trunk/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 <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 "rad_option.hh"
  23
  24 namespace voro {
  25
  26 /** \brief Pure virtual class from which wall objects are derived.
  27  *
  28  * This is a pure virtual class for a generic wall object. A wall object
  29  * can be specified by deriving a new class from this and specifying the
  30  * functions.*/
  31 class wall {
  32         public:
  33                 virtual ~wall() {}
  34                 /** A pure virtual function for testing whether a point is
  35                  * inside the wall object. */
  36                 virtual bool point_inside(double x,double y,double z) = 0;
  37                 /** A pure virtual function for cutting a cell without
  38                  * neighbor-tracking with a wall. */
  39                 virtual bool cut_cell(voronoicell &c,double x,double y,double z) = 0;
  40                 /** A pure virtual function for cutting a cell with
  41                  * neighbor-tracking enabled with a wall. */
  42                 virtual bool cut_cell(voronoicell_neighbor &c,double x,double y,double z) = 0;
  43 };
  44
  45 /** \brief A class for storing a list of pointers to walls.
  46  *
  47  * This class stores a list of pointers to wall classes. It contains several
  48  * simple routines that make use of the wall classes (such as telling whether a
  49  * given position is inside all of the walls or not). It can be used by itself,
  50  * but also forms part of container_base, for associating walls with this
  51  * class. */
  52 class wall_list {
  53         public:
  54                 /** An array holding pointers to wall objects. */
  55                 wall **walls;
  56                 /** A pointer to the next free position to add a wall pointer.
  57                  */
  58                 wall **wep;
  59                 wall_list();
  60                 ~wall_list();
  61                 /** Adds a wall to the list.
  62                  * \param[in] w the wall to add. */
  63                 inline void add_wall(wall *w) {
  64                         if(wep==wel) increase_wall_memory();
  65                         *(wep++)=w;
  66                 }
  67                 /** Adds a wall to the list.
  68                  * \param[in] w a reference to the wall to add. */
  69                 inline void add_wall(wall &w) {add_wall(&w);}
  70                 void add_wall(wall_list &wl);
  71                 /** Determines whether a given position is inside all of the
  72                  * walls on the list.
  73                  * \param[in] (x,y,z) the position to test.
  74                  * \return True if it is inside, false if it is outside. */
  75                 inline bool point_inside_walls(double x,double y,double z) {
  76                         for(wall **wp=walls;wp<wep;wp++) if(!((*wp)->point_inside(x,y,z))) return false;
  77                         return true;
  78                 }
  79                 /** Cuts a Voronoi cell by all of the walls currently on
  80                  * the list.
  81                  * \param[in] c a reference to the Voronoi cell class.
  82                  * \param[in] (x,y,z) the position of the cell.
  83                  * \return True if the cell still exists, false if the cell is
  84                  * deleted. */
  85                 template<class c_class>
  86                 bool apply_walls(c_class &c,double x,double y,double z) {
  87                         for(wall **wp=walls;wp<wep;wp++) if(!((*wp)->cut_cell(c,x,y,z))) return false;
  88                         return true;
  89                 }
  90                 void deallocate();
  91         protected:
  92                 void increase_wall_memory();
  93                 /** A pointer to the limit of the walls array, used to
  94                  * determine when array is full. */
  95                 wall **wel;
  96                 /** The current amount of memory allocated for walls. */
  97                 int current_wall_size;
  98 };
  99
 100 /** \brief Class for representing a particle system in a three-dimensional
 101  * rectangular box.
 102  *
 103  * This class represents a system of particles in a three-dimensional
 104  * rectangular box. Any combination of non-periodic and periodic coordinates
 105  * can be used in the three coordinate directions. The class is not intended
 106  * for direct use, but instead forms the base of the container and
 107  * container_poly classes that add specialized routines for computing the
 108  * regular and radical Voronoi tessellations respectively. It contains routines
 109  * that are commonly between these two classes, such as those for drawing the
 110  * domain, and placing particles within the internal data structure.
 111  *
 112  * The class is derived from the wall_list class, which encapsulates routines
 113  * for associating walls with the container, and the voro_base class, which
 114  * encapsulates routines about the underlying computational grid. */
 115 class container_base : public voro_base, public wall_list {
 116         public:
 117                 /** The minimum x coordinate of the container. */
 118                 const double ax;
 119                 /** The maximum x coordinate of the container. */
 120                 const double bx;
 121                 /** The minimum y coordinate of the container. */
 122                 const double ay;
 123                 /** The maximum y coordinate of the container. */
 124                 const double by;
 125                 /** The minimum z coordinate of the container. */
 126                 const double az;
 127                 /** The maximum z coordinate of the container. */
 128                 const double bz;
 129                 /** The maximum length squared that could be encountered in the
 130                  * Voronoi cell calculation. */
 131                 const double max_len_sq;
 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                 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(*this);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(*this);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(*this);
 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(*this);
 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                 /** Computes the Voronoi cell for a ghost particle at a given
 480                  * location.
 481                  * \param[out] c a Voronoi cell class in which to store the
 482                  *               computed cell.
 483                  * \param[in] (x,y,z) the location of the ghost particle.
 484                  * \return True if the cell was computed. If the cell cannot be
 485                  * computed, if it is removed entirely by a wall or boundary
 486                  * condition, then the routine returns false. */
 487                 template<class v_cell>
 488                 inline bool compute_ghost_cell(v_cell &c,double x,double y,double z) {
 489                         int ijk;
 490                         if(put_locate_block(ijk,x,y,z)) {
 491                                 double *pp=p[ijk]+3*co[ijk]++;
 492                                 *(pp++)=x;*(pp++)=y;*pp=z;
 493                                 bool q=compute_cell(c,ijk,co[ijk]-1);
 494                                 co[ijk]--;
 495                                 return q;
 496                         }
 497                         return false;
 498                 }
 499         private:
 500                 voro_compute<container> vc;
 501                 friend class voro_compute<container>;
 502 };
 503
 504 /** \brief Extension of the container_base class for computing radical Voronoi
 505  * tessellations.
 506  *
 507  * This class is an extension of container_base class that has routines
 508  * specifically for computing the radical Voronoi tessellation that depends on
 509  * the particle radii. */
 510 class container_poly : public container_base, public radius_poly {
 511         public:
 512                 container_poly(double ax_,double bx_,double ay_,double by_,double az_,double bz_,
 513                                 int nx_,int ny_,int nz_,bool xperiodic_,bool yperiodic_,bool zperiodic_,int init_mem);
 514                 void clear();
 515                 void put(int n,double x,double y,double z,double r);
 516                 void put(particle_order &vo,int n,double x,double y,double z,double r);
 517                 void import(FILE *fp=stdin);
 518                 void import(particle_order &vo,FILE *fp=stdin);
 519                 /** Imports a list of particles from an open file stream into
 520                  * the container_poly class. Entries of five numbers (Particle
 521                  * ID, x position, y position, z position, radius) are searched
 522                  * for. If the file cannot be successfully read, then the
 523                  * routine causes a fatal error.
 524                  * \param[in] filename the name of the file to open and read
 525                  *                     from. */
 526                 inline void import(const char* filename) {
 527                         FILE *fp=safe_fopen(filename,"r");
 528                         import(fp);
 529                         fclose(fp);
 530                 }
 531                 /** Imports a list of particles from an open file stream into
 532                  * the container_poly class. Entries of five numbers (Particle
 533                  * ID, x position, y position, z position, radius) are searched
 534                  * for. In addition, the order in which particles are read is
 535                  * saved into an ordering class. If the file cannot be
 536                  * successfully read, then the routine causes a fatal error.
 537                  * \param[in,out] vo the ordering class to use.
 538                  * \param[in] filename the name of the file to open and read
 539                  *                     from. */
 540                 inline void import(particle_order &vo,const char* filename) {
 541                         FILE *fp=safe_fopen(filename,"r");
 542                         import(vo,fp);
 543                         fclose(fp);
 544                 }
 545                 void compute_all_cells();
 546                 double sum_cell_volumes();
 547                 /** Dumps particle IDs, positions and radii to a file.
 548                  * \param[in] vl the loop class to use.
 549                  * \param[in] fp a file handle to write to. */
 550                 template<class c_loop>
 551                 void draw_particles(c_loop &vl,FILE *fp) {
 552                         double *pp;
 553                         if(vl.start()) do {
 554                                 pp=p[vl.ijk]+4*vl.q;
 555                                 fprintf(fp,"%d %g %g %g %g\n",id[vl.ijk][vl.q],*pp,pp[1],pp[2],pp[3]);
 556                         } while(vl.inc());
 557                 }
 558                 /** Dumps all of the particle IDs, positions and radii to a
 559                  * file.
 560                  * \param[in] fp a file handle to write to. */
 561                 inline void draw_particles(FILE *fp=stdout) {
 562                         c_loop_all vl(*this);
 563                         draw_particles(vl,fp);
 564                 }
 565                 /** Dumps all of the particle IDs, positions and radii to a
 566                  * file.
 567                  * \param[in] filename the name of the file to write to. */
 568                 inline void draw_particles(const char *filename) {
 569                         FILE *fp=safe_fopen(filename,"w");
 570                         draw_particles(fp);
 571                         fclose(fp);
 572                 }
 573                 /** Dumps particle positions in POV-Ray format.
 574                  * \param[in] vl the loop class to use.
 575                  * \param[in] fp a file handle to write to. */
 576                 template<class c_loop>
 577                 void draw_particles_pov(c_loop &vl,FILE *fp) {
 578                         double *pp;
 579                         if(vl.start()) do {
 580                                 pp=p[vl.ijk]+4*vl.q;
 581                                 fprintf(fp,"// id %d\nsphere{<%g,%g,%g>,%g}\n",
 582                                                 id[vl.ijk][vl.q],*pp,pp[1],pp[2],pp[3]);
 583                         } while(vl.inc());
 584                 }
 585                 /** Dumps all the particle positions in POV-Ray format.
 586                  * \param[in] fp a file handle to write to. */
 587                 inline void draw_particles_pov(FILE *fp=stdout) {
 588                         c_loop_all vl(*this);
 589                         draw_particles_pov(vl,fp);
 590                 }
 591                 /** Dumps all the particle positions in POV-Ray format.
 592                  * \param[in] filename the name of the file to write to. */
 593                 inline void draw_particles_pov(const char *filename) {
 594                         FILE *fp=safe_fopen(filename,"w");
 595                         draw_particles_pov(fp);
 596                         fclose(fp);
 597                 }
 598                 /** Computes Voronoi cells and saves the output in gnuplot
 599                  * format.
 600                  * \param[in] vl the loop class to use.
 601                  * \param[in] fp a file handle to write to. */
 602                 template<class c_loop>
 603                 void draw_cells_gnuplot(c_loop &vl,FILE *fp) {
 604                         voronoicell c;double *pp;
 605                         if(vl.start()) do if(compute_cell(c,vl)) {
 606                                 pp=p[vl.ijk]+ps*vl.q;
 607                                 c.draw_gnuplot(*pp,pp[1],pp[2],fp);
 608                         } while(vl.inc());
 609                 }
 610                 /** Compute all Voronoi cells and saves the output in gnuplot
 611                  * format.
 612                  * \param[in] fp a file handle to write to. */
 613                 inline void draw_cells_gnuplot(FILE *fp=stdout) {
 614                         c_loop_all vl(*this);
 615                         draw_cells_gnuplot(vl,fp);
 616                 }
 617                 /** Compute all Voronoi cells and saves the output in gnuplot
 618                  * format.
 619                  * \param[in] filename the name of the file to write to. */
 620                 inline void draw_cells_gnuplot(const char *filename) {
 621                         FILE *fp=safe_fopen(filename,"w");
 622                         draw_cells_gnuplot(fp);
 623                         fclose(fp);
 624                 }
 625                 /** Computes Voronoi cells and saves the output in POV-Ray
 626                  * format.
 627                  * \param[in] vl the loop class to use.
 628                  * \param[in] fp a file handle to write to. */
 629                 template<class c_loop>
 630                 void draw_cells_pov(c_loop &vl,FILE *fp) {
 631                         voronoicell c(*this);double *pp;
 632                         if(vl.start()) do if(compute_cell(c,vl)) {
 633                                 fprintf(fp,"// cell %d\n",id[vl.ijk][vl.q]);
 634                                 pp=p[vl.ijk]+ps*vl.q;
 635                                 c.draw_pov(*pp,pp[1],pp[2],fp);
 636                         } while(vl.inc());
 637                 }
 638                 /** Computes all Voronoi cells and saves the output in POV-Ray
 639                  * format.
 640                  * \param[in] fp a file handle to write to. */
 641                 inline void draw_cells_pov(FILE *fp=stdout) {
 642                         c_loop_all vl(*this);
 643                         draw_cells_pov(vl,fp);
 644                 }
 645                 /** Computes all Voronoi cells and saves the output in POV-Ray
 646                  * format.
 647                  * \param[in] filename the name of the file to write to. */
 648                 inline void draw_cells_pov(const char *filename) {
 649                         FILE *fp=safe_fopen(filename,"w");
 650                         draw_cells_pov(fp);
 651                         fclose(fp);
 652                 }
 653                 /** Computes the Voronoi cells and saves customized information
 654                  * about them.
 655                  * \param[in] vl the loop class to use.
 656                  * \param[in] format the custom output string to use.
 657                  * \param[in] fp a file handle to write to. */
 658                 template<class c_loop>
 659                 void print_custom(c_loop &vl,const char *format,FILE *fp) {
 660                         int ijk,q;double *pp;
 661                         if(contains_neighbor(format)) {
 662                                 voronoicell_neighbor c(*this);
 663                                 if(vl.start()) do if(compute_cell(c,vl)) {
 664                                         ijk=vl.ijk;q=vl.q;pp=p[ijk]+ps*q;
 665                                         c.output_custom(format,id[ijk][q],*pp,pp[1],pp[2],pp[3],fp);
 666                                 } while(vl.inc());
 667                         } else {
 668                                 voronoicell c(*this);
 669                                 if(vl.start()) do if(compute_cell(c,vl)) {
 670                                         ijk=vl.ijk;q=vl.q;pp=p[ijk]+ps*q;
 671                                         c.output_custom(format,id[ijk][q],*pp,pp[1],pp[2],pp[3],fp);
 672                                 } while(vl.inc());
 673                         }
 674                 }
 675                 /** Computes the Voronoi cell for a particle currently being
 676                  * referenced by a loop class.
 677                  * \param[out] c a Voronoi cell class in which to store the
 678                  *               computed cell.
 679                  * \param[in] vl the loop class to use.
 680                  * \return True if the cell was computed. If the cell cannot be
 681                  * computed, if it is removed entirely by a wall or boundary
 682                  * condition, then the routine returns false. */
 683                 template<class v_cell,class c_loop>
 684                 inline bool compute_cell(v_cell &c,c_loop &vl) {
 685                         return vc.compute_cell(c,vl.ijk,vl.q,vl.i,vl.j,vl.k);
 686                 }
 687                 /** Computes the Voronoi cell for given particle.
 688                  * \param[out] c a Voronoi cell class in which to store the
 689                  *               computed cell.
 690                  * \param[in] ijk the block that the particle is within.
 691                  * \param[in] q the index of the particle within the block.
 692                  * \return True if the cell was computed. If the cell cannot be
 693                  * computed, if it is removed entirely by a wall or boundary
 694                  * condition, then the routine returns false. */
 695                 template<class v_cell>
 696                 inline bool compute_cell(v_cell &c,int ijk,int q) {
 697                         int k=ijk/nxy,ijkt=ijk-nxy*k,j=ijkt/nx,i=ijkt-j*nx;
 698                         return vc.compute_cell(c,ijk,q,i,j,k);
 699                 }
 700                 /** Computes the Voronoi cell for a ghost particle at a given
 701                  * location.
 702                  * \param[out] c a Voronoi cell class in which to store the
 703                  *               computed cell.
 704                  * \param[in] (x,y,z) the location of the ghost particle.
 705                  * \param[in] r the radius of the ghost particle.
 706                  * \return True if the cell was computed. If the cell cannot be
 707                  * computed, if it is removed entirely by a wall or boundary
 708                  * condition, then the routine returns false. */
 709                 template<class v_cell>
 710                 inline bool compute_ghost_cell(v_cell &c,double x,double y,double z,double r) {
 711                         int ijk;
 712                         if(put_locate_block(ijk,x,y,z)) {
 713                                 double *pp=p[ijk]+4*co[ijk]++,tm=max_radius;
 714                                 *(pp++)=x;*(pp++)=y;*(pp++)=z;*pp=r;
 715                                 if(r>max_radius) max_radius=r;
 716                                 bool q=compute_cell(c,ijk,co[ijk]-1);
 717                                 co[ijk]--;max_radius=tm;
 718                                 return q;
 719                         }
 720                         return false;
 721                 }
 722                 void print_custom(const char *format,FILE *fp=stdout);
 723                 void print_custom(const char *format,const char *filename);
 724                 bool find_voronoi_cell(double x,double y,double z,double &rx,double &ry,double &rz,int &pid);
 725         private:
 726                 voro_compute<container_poly> vc;
 727                 friend class voro_compute<container_poly>;
 728 };
 729
 730 }
 731
 732 #endif