branches/2d/src/v_connect/ctr_boundary_2d.hh

   1 // Voro++, a cell-based Voronoi library
   2 //
   3 // Authors  : Chris H. Rycroft (LBL / UC Berkeley)
   4 //            Cody Robert Dance (UC Berkeley)
   5 // Email    : chr@alum.mit.edu
   6 // Date     : August 30th 2011
   7
   8 /** \file ctr_boundary_2d.hh
   9  * \brief Header file for the container_boundary_2d and related classes. */
  10
  11 #ifndef VOROPP_CTR_BOUNDARY_2D_HH
  12 #define VOROPP_CTR_BOUNDARY_2D_HH
  13
  14 #include <cstdio>
  15 #include <cstdlib>
  16 #include <cmath>
  17 #include <vector>
  18 using namespace std;
  19
  20 #include "config.hh"
  21 #include "common.hh"
  22 #include "v_base_2d.hh"
  23 #include "cell_2d.hh"
  24 #include "c_loops_2d.hh"
  25 #include "rad_option.hh"
  26 #include "v_compute_2d.hh"
  27
  28 namespace voro {
  29
  30 /** \brief Class for representing a particle system in a three-dimensional
  31  * rectangular box.
  32  *
  33  * This class represents a system of particles in a three-dimensional
  34  * rectangular box. Any combination of non-periodic and periodic coordinates
  35  * can be used in the three coordinate directions. The class is not intended
  36  * for direct use, but instead forms the base of the container and
  37  * container_poly classes that add specialized routines for computing the
  38  * regular and radical Voronoi tessellations respectively. It contains routines
  39  * that are commonly between these two classes, such as those for drawing the
  40  * domain, and placing particles within the internal data structure.
  41  *
  42  * The class is derived from the wall_list class, which encapsulates routines
  43  * for associating walls with the container, and the voro_base class, which
  44  * encapsulates routines about the underlying computational grid. */
  45 class container_boundary_2d : public voro_base_2d, public radius_mono {
  46         public:
  47
  48                 /** The minimum x coordinate of the container. */
  49                 const double ax;
  50                 /** The maximum x coordinate of the container. */
  51                 const double bx;
  52                 /** The minimum y coordinate of the container. */
  53                 const double ay;
  54                 /** The maximum y coordinate of the container. */
  55                 const double by;
  56                 /** A boolean value that determines if the x coordinate in
  57                  * periodic or not. */
  58                 const bool xperiodic;
  59                 /** A boolean value that determines if the y coordinate in
  60                  * periodic or not. */
  61                 const bool yperiodic;
  62                 /** This array holds the numerical IDs of each particle in each
  63                  * computational box. */
  64                 int **id;
  65                 /** A two dimensional array holding particle positions. For the
  66                  * derived container_poly class, this also holds particle
  67                  * radii. */
  68                 double **p;
  69                 /** This array holds the number of particles within each
  70                  * computational box of the container. */
  71                 int *co;
  72                 /** This array holds the maximum amount of particle memory for
  73                  * each computational box of the container. If the number of
  74                  * particles in a particular box ever approaches this limit,
  75                  * more is allocated using the add_particle_memory() function.
  76                  */
  77                 int *mem;
  78                 int **wid;
  79                 int **nlab;
  80                 int ***plab;
  81                 int **bndpts;
  82                 int boundary_track;
  83                 int edbc;
  84                 int edbm;
  85                 int *edb;
  86                 double *bnds;
  87
  88                 /** The amount of memory in the array structure for each
  89                  * particle. This is set to 2 when the basic class is
  90                  * initialized, so that the array holds (x,y) positions. If the
  91                  * 2D container class is initialized as part of the derived class
  92                  * container_poly_2d, then this is set to 3, to also hold the
  93                  * particle radii. */
  94                 const int ps;
  95                 container_boundary_2d(double ax_,double bx_,double ay_,double by_,
  96                              int nx_,int ny_,bool xperiodic_,bool yperiodic_,int init_mem);
  97                 ~container_boundary_2d();
  98                 void region_count();
  99                 /** Initializes the Voronoi cell prior to a compute_cell
 100                  * operation for a specific particle being carried out by a
 101                  * voro_compute class. The cell is initialized to fill the
 102                  * entire container. For non-periodic coordinates, this is set
 103                  * by the position of the walls. For periodic coordinates, the
 104                  * space is equally divided in either direction from the
 105                  * particle's initial position. Plane cuts made by any walls
 106                  * that have been added are then applied to the cell.
 107                  o* \param[in,out] c a reference to a voronoicell_nonconvex_2d object.
 108                  * \param[in] ij the block that the particle is within.
 109                  * \param[in] q the index of the particle within its block.
 110                  * \param[in] (ci,cj) the coordinates of the block in the
 111                  *                    container coordinate system.
 112                  * \param[out] (i,j) the coordinates of the test block relative
 113                  *                   to the voro_compute coordinate system.
 114                  * \param[out] (x,y) the position of the particle.
 115                  * \param[out] disp a block displacement used internally by the
 116                  *                  compute_cell routine.
 117                  * \return False if the plane cuts applied by walls completely
 118                  * removed the cell, true otherwise. */
 119                 template<class v_cell_2d>
 120                 inline bool initialize_voronoicell(v_cell_2d &c,int ij,int q,int ci,int cj,
 121                                 int &i,int &j,double &x,double &y,int &disp) {
 122                         double x1,x2,y1,y2,*pp=p[ij]+ps*q;
 123                                 x=*(pp++);y=*(pp++);
 124                         if(xperiodic) {x1=-(x2=0.5*(bx-ax));i=nx;} else {x1=ax-x;x2=bx-x;i=ci;}
 125                         if(yperiodic) {y1=-(y2=0.5*(by-ay));j=ny;} else {y1=ay-y;y2=by-y;j=cj;}
 126                         if(bndpts[ij][q]==-1) c.init(x1,x2,y1,y2);
 127                         else {
 128                                 int &bid=bndpts[ij][q];
 129                                 double cx=bnds[2*bid],cy=bnds[2*bid+1];
 130                                 int nwid=edb[2*bid],lwid=edb[2*bid+1];
 131                                 double lx=bnds[2*lwid],ly=bnds[2*lwid+1];
 132                                 double nx=bnds[2*nwid],ny=bnds[2*nwid+1];
 133                                 c.init_nonconvex(x1,x2,y1,y2,nx-cx,ny-cy,lx-cx,ly-cy);
 134                         }
 135                         disp=ij-i-nx*j;
 136                         return true;
 137                 }
 138                 template<class v_cell_2d>
 139                 inline void connect_to_cell(v_cell_2d &c,int ij, int q){
 140                         if(full_connect){
 141                                 c.set_id(id[ij][q]);
 142                                 c.full_connect_on();
 143                         }
 144                 }
 145
 146                 bool point_inside(double x,double y);
 147                 template<class v_cell_2d>
 148                 bool boundary_cuts(v_cell_2d &c,int ij,double x,double y);
 149                 /** Initializes parameters for a find_voronoi_cell call within
 150                  * the voro_compute template.
 151                  * \param[in] (ci,cj) the coordinates of the test block in
 152                  *                    the container coordinate system.
 153                  * \param[in] ij the index of the test block
 154                  * \param[out] (i,j) the coordinates of the test block relative
 155                  *                   to the voro_compute coordinate system.
 156                  * \param[out] disp a block displacement used internally by the
 157                  *                  find_voronoi_cell routine. */
 158                 inline void initialize_search(int ci,int cj,int ij,int &i,int &j,int &disp) {
 159                         i=xperiodic?nx:ci;
 160                         j=yperiodic?ny:cj;
 161                         disp=ij-i-nx*j;
 162                 }
 163                 /** Returns the position of a particle currently being computed
 164                  * relative to the computational block that it is within. It is
 165                  * used to select the optimal worklist entry to use.
 166                  * \param[in] (x,y) the position of the particle.
 167                  * \param[in] (ci,cj) the block that the particle is within.
 168                  * \param[out] (fx,fy) the position relative to the block.
 169                  */
 170                 inline void frac_pos(double x,double y,double ci,double cj,
 171                                 double &fx,double &fy) {
 172                         fx=x-ax-boxx*ci;
 173                         fy=y-ay-boxy*cj;
 174                 }
 175                 /** Calculates the index of block in the container structure
 176                  * corresponding to given coordinates.
 177                  * \param[in] (ci,cj) the coordinates of the original block in
 178                  *                    the current computation, relative to the
 179                  *                    container coordinate system.
 180                  * \param[in] (ei,ej) the displacement of the current block
 181                  *                    from the original block.
 182                  * \param[in,out] (qx,qy) the periodic displacement that must
 183                  *                        be added to the particles within the
 184                  *                        computed block.
 185                  * \param[in] disp a block displacement used internally by the
 186                  *                 find_voronoi_cell andcompute_cell routines.
 187                  * \return The block index. */
 188                 inline int region_index(int ci,int cj,int ei,int ej,double &qx,double &qy,int &disp) {
 189                         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;}
 190                         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;}
 191                         return disp+ei+nx*ej;
 192                 }
 193                 void draw_domain_gnuplot(FILE *fp=stdout);
 194                 /** Draws an outline of the domain in Gnuplot format.
 195                  * \param[in] filename the filename to write to. */
 196                 inline void draw_domain_gnuplot(const char* filename) {
 197                         FILE *fp=safe_fopen(filename,"w");
 198                         draw_domain_gnuplot(fp);
 199                         fclose(fp);
 200                 }
 201                 void draw_domain_pov(FILE *fp=stdout);
 202                 /** Draws an outline of the domain in Gnuplot format.
 203                  * \param[in] filename the filename to write to. */
 204                 inline void draw_domain_pov(const char* filename) {
 205                         FILE *fp=safe_fopen(filename,"w");
 206                         draw_domain_pov(fp);
 207                         fclose(fp);
 208                 }
 209                 void draw_boundary_gnuplot(FILE *fp=stdout);
 210                 inline void draw_boundary_gnuplot(const char* filename) {
 211                         FILE *fp=safe_fopen(filename,"w");
 212                         draw_boundary_gnuplot(fp);
 213                         fclose(fp);
 214                 }
 215                 /** Sums up the total number of stored particles.
 216                  * \return The number of particles. */
 217                 inline int total_particles() {
 218                         int tp=*co;
 219                         for(int *cop=co+1;cop<co+nxy;cop++) tp+=*cop;
 220                         return tp;
 221                 }
 222                 inline void start_boundary() {boundary_track=edbc;}
 223                 void end_boundary();
 224                 void register_boundary(double x,double y);
 225                 void clear();
 226                 void put(int n,double x,double y);
 227                 void put(particle_order &vo,int n,double x,double y);
 228
 229
 230
 231
 232
 233
 234
 235                 void import(FILE *fp=stdin);
 236                 /** Imports a list of particles from an open file stream into
 237                  * the container. Entries of three numbers (Particle ID, x
 238                  * position, y position) are searched for. If the file cannot
 239                  * be successfully read, then the routine causes a fatal error.
 240                  * \param[in] filename the name of the file to open and read
 241                  *                     from. */
 242                 inline void import(const char* filename) {
 243                         FILE *fp=safe_fopen(filename,"r");
 244                         import(fp);
 245                         fclose(fp);
 246                 }
 247                 void compute_all_cells();
 248                 double sum_cell_areas();
 249                 /** Dumps particle IDs and positions to a file.
 250                  * \param[in] vl the loop class to use.
 251                  * \param[in] fp a file handle to write to. */
 252                 template<class c_loop_2d>
 253                 void draw_particles(c_loop_2d &vl,FILE *fp) {
 254                         double *pp;
 255                         if(vl.start()) do {
 256                                 pp=p[vl.ij]+2*vl.q;
 257                                 fprintf(fp,"%d %g %g\n",id[vl.ij][vl.q],*pp,pp[1]);
 258                         } while(vl.inc());
 259                 }
 260                 /** Dumps all of the particle IDs and positions to a file.
 261                  * \param[in] fp a file handle to write to. */
 262                 inline void draw_particles(FILE *fp=stdout) {
 263                         c_loop_all_2d vl(*this);
 264                         draw_particles(vl,fp);
 265                 }
 266                 /** Dumps all of the particle IDs and positions to a file.
 267                  * \param[in] filename the name of the file to write to. */
 268                 inline void draw_particles(const char *filename) {
 269                         FILE *fp=safe_fopen(filename,"w");
 270                         draw_particles(fp);
 271                         fclose(fp);
 272                 }
 273                 /** Dumps particle positions in POV-Ray format.
 274                  * \param[in] vl the loop class to use.
 275                  * \param[in] fp a file handle to write to. */
 276                 template<class c_loop_2d>
 277                 void draw_particles_pov(c_loop_2d &vl,FILE *fp) {
 278                         double *pp;
 279                         if(vl.start()) do {
 280                                 pp=p[vl.ij]+2*vl.q;
 281                                 fprintf(fp,"// id %d\nsphere{<%g,%g,0>,s}\n",
 282                                                 id[vl.ij][vl.q],*pp,pp[1]);
 283                         } while(vl.inc());
 284                 }
 285                 /** Dumps all particle positions in POV-Ray format.
 286                  * \param[in] fp a file handle to write to. */
 287                 inline void draw_particles_pov(FILE *fp=stdout) {
 288                         c_loop_all_2d vl(*this);
 289                         draw_particles_pov(vl,fp);
 290                 }
 291                 /** Dumps all particle positions in POV-Ray format.
 292                  * \param[in] filename the name of the file to write to. */
 293                 inline void draw_particles_pov(const char *filename) {
 294                         FILE *fp=safe_fopen(filename,"w");
 295                         draw_particles_pov(fp);
 296                         fclose(fp);
 297                 }
 298                 /** Computes Voronoi cells and saves the output in gnuplot
 299                  * format.
 300                  * \param[in] vl the loop class to use.
 301                  * \param[in] fp a file handle to write to. */
 302                 template<class c_loop_2d>
 303                 void draw_cells_gnuplot(c_loop_2d &vl,FILE *fp) {
 304                         voronoicell_nonconvex_2d c;double *pp;
 305                         if(vl.start()) do if(compute_cell(c,vl)) {
 306                                 pp=p[vl.ij]+ps*vl.q;
 307                                 fprintf(fp,"# [%d]\n",id[vl.ij][vl.q]);
 308                                 c.draw_gnuplot(*pp,pp[1],fp);
 309                                 fputs("\n",fp);
 310                         } while(vl.inc());
 311                 }
 312                 /** Computes all Voronoi cells and saves the output in gnuplot
 313                  * format.
 314                  * \param[in] fp a file handle to write to. */
 315                 inline void draw_cells_gnuplot(FILE *fp=stdout) {
 316                         c_loop_all_2d vl(*this);
 317                         draw_cells_gnuplot(vl,fp);
 318                 }
 319                 /** Computes all Voronoi cells and saves the output in gnuplot
 320                  * format.
 321                  * \param[in] filename the name of the file to write to. */
 322                 inline void draw_cells_gnuplot(const char *filename) {
 323                         FILE *fp=safe_fopen(filename,"w");
 324                         draw_cells_gnuplot(fp);
 325                         fclose(fp);
 326                 }
 327                 /** Computes Voronoi cells and saves the output in POV-Ray
 328                  * format.
 329                  * \param[in] vl the loop class to use.
 330                  * \param[in] fp a file handle to write to. */
 331                 template<class c_loop_2d>
 332                 void draw_cells_pov(c_loop_2d &vl,FILE *fp) {
 333                         voronoicell_nonconvex_2d c;double *pp;
 334                         if(vl.start()) do if(compute_cell(c,vl)) {
 335                                 fprintf(fp,"// cell %d\n",id[vl.ij][vl.q]);
 336                                 pp=p[vl.ij]+ps*vl.q;
 337                                 c.draw_pov(*pp,pp[1],fp);
 338                         } while(vl.inc());
 339                 }
 340                 /** Computes all Voronoi cells and saves the output in POV-Ray
 341                  * format.
 342                  * \param[in] fp a file handle to write to. */
 343                 inline void draw_cells_pov(FILE *fp=stdout) {
 344                         c_loop_all_2d vl(*this);
 345                         draw_cells_pov(vl,fp);
 346                 }
 347                 /** Computes all Voronoi cells and saves the output in POV-Ray
 348                  * format.
 349                  * \param[in] filename the name of the file to write to. */
 350                 inline void draw_cells_pov(const char *filename) {
 351                         FILE *fp=safe_fopen(filename,"w");
 352                         draw_cells_pov(fp);
 353                         fclose(fp);
 354                 }
 355                 /** Computes the Voronoi cells and saves customized information
 356                  * about them.
 357                  * \param[in] vl the loop class to use.
 358                  * \param[in] format the custom output string to use.
 359                  * \param[in] fp a file handle to write to. */
 360                 template<class c_loop_2d>
 361                 void print_custom(c_loop_2d &vl,const char *format,FILE *fp) {
 362                         int ij,q;double *pp;
 363                 //      bool glob=false, loc=false;
 364                 //      if(contains_neighbor_global(format)){
 365                 //              init_globne();
 366                 //              glob=true;
 367                 //      }
 368                         if(contains_neighbor(format)){
 369                 //              loc=true;
 370                 //      }
 371                 //      if(glob || loc) {
 372                                 voronoicell_nonconvex_neighbor_2d c;
 373                                 if(vl.start()) do if(compute_cell(c,vl)) {
 374
 375                                         ij=vl.ij;q=vl.q;pp=p[ij]+ps*q;
 376                         //              if(glob) add_globne_info(id[ij][q], c.ne, c.p);
 377                                         c.output_custom(format,id[ij][q],*pp,pp[1],default_radius_2d,fp);
 378                                 } while(vl.inc());
 379                         //      if(glob) print_globne(fp);
 380                         } else {
 381                                 voronoicell_nonconvex_2d c;
 382                                 if(vl.start()) do if(compute_cell(c,vl)) {
 383                                         ij=vl.ij;q=vl.q;pp=p[ij]+ps*q;
 384                                         c.output_custom(format,id[ij][q],*pp,pp[1],default_radius_2d,fp);
 385                                 } while(vl.inc());
 386                         }
 387                 }
 388                 void print_custom(const char *format,FILE *fp=stdout);
 389                 void print_custom(const char *format,const char *filename);
 390                 //bool find_voronoi_cell(double x,double y,double &rx,double &ry,int &pid);
 391                 /** Computes the Voronoi cell for a particle currently being
 392                  * referenced by a loop class.
 393                  * \param[out] c a Voronoi cell class in which to store the
 394                  *               computed cell.
 395                  * \param[in] vl the loop class to use.
 396                  * \return True if the cell was computed. If the cell cannot be
 397                  * computed, if it is removed entirely by a wall or boundary
 398                  * condition, then the routine returns false. */
 399                 template<class v_cell_2d,class c_loop_2d>
 400                 inline bool compute_cell(v_cell_2d &c,c_loop_2d &vl) {
 401                         return vc.compute_cell(c,vl.ij,vl.q,vl.i,vl.j);
 402                 }
 403                 /** Computes the Voronoi cell for given particle.
 404                  * \param[out] c a Voronoi cell class in which to store the
 405                  *               computed cell.
 406                  * \param[in] ij the block that the particle is within.
 407                  * \param[in] q the index of the particle within the block.
 408                  * \return True if the cell was computed. If the cell cannot be
 409                  * computed, if it is removed entirely by a wall or boundary
 410                  * condition, then the routine returns false. */
 411                 template<class v_cell_2d>
 412                 inline bool compute_cell(v_cell_2d &c,int ij,int q) {
 413                         int j=ij/nx,i=ij-j*nx;
 414                         return vc.compute_cell(c,ij,q,i,j);
 415                 }
 416                 void setup();
 417                 bool skip(int ij,int l,double x,double y);
 418         private:
 419                 inline void draw_block(int ij) {
 420                         int i=ij%nx,j=ij/nx;
 421                         double lx=ax+i*boxx,ly=ay+j*boxy,ux=lx+boxx,uy=ly+boxy;
 422                         printf("%g %g\n%g %g\n%g %g\n%g %g\n%g %g\n\n\n",lx,ly,ux,ly,ux,uy,lx,uy,lx,ly);
 423                 }
 424                 int *soi;
 425                 int *tmp;
 426                 int *tmpp;
 427                 int *tmpe;
 428                 void create_label_table();
 429                 void add_particle_memory(int i);
 430                 inline bool put_locate_block(int &ij,double &x,double &y);
 431                 inline bool put_remap(int &ij,double &x,double &y);
 432                 inline bool remap(int &ai,int &aj,int &ci,int &cj,double &x,double &y,int &ij);
 433                 void add_temporary_label_memory();
 434                 void add_boundary_memory();
 435                 void tag_line(int &ij,int ije,int wid);
 436                 inline void tag(int ij,int wid);
 437                 void tag_walls(double x1,double y1,double x2,double y2,int wid);
 438                 inline bool cross_product(double x1,double y1,double x2,double y2) {
 439                         return x1*y2>x2*y1;
 440                 }
 441                 void semi_circle_labeling(double x1,double y1,double x2,double y2,int bid);
 442                 voro_compute_2d<container_boundary_2d> vc;
 443                 friend class voro_compute_2d<container_boundary_2d>;
 444 };
 445
 446 }
 447
 448 #endif