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