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