branches/2d_boundary/src/container_2d.hh

   1 /** \file container_2d.hh
   2  * \brief Header file for the container_2d class. */
   3
   4 #ifndef VOROPP_CONTAINER_2D_HH
   5 #define VOROPP_CONTAINER_2D_HH
   6
   7 #include <cstdio>
   8 #include <cstdlib>
   9 #include <cmath>
  10 using namespace std;
  11
  12 #include "config.hh"
  13 #include "cell_2d.hh"
  14
  15 class voropp_loop_2d;
  16
  17 /** \brief A class representing the whole 2D simulation region.
  18  *
  19  * The container class represents the whole simulation region. The
  20  * container constructor sets up the geometry and periodicity, and divides
  21  * the geometry into rectangular grid of blocks, each of which handles the
  22  * particles in a particular area. Routines exist for putting in particles,
  23  * importing particles from standard input, and carrying out Voronoi
  24  * calculations. */
  25 class container_2d {
  26         public:
  27                 /** degbugging aid. */
  28                 bool debug;
  29                 /** The minimum x coordinate of the container. */
  30                 const double ax;
  31                 /** The maximum x coordinate of the container. */
  32                 const double bx;
  33                 /** The minimum y coordinate of the container. */
  34                 const double ay;
  35                 /** The maximum y coordinate of the container. */
  36                 const double by;
  37                 /** The box length in the x direction, set to (bx-ax)/nx. */
  38                 const double boxx;
  39                 /** The box length in the y direction, set to (by-ay)/ny. */
  40                 const double boxy;
  41                 /** The inverse box length in the x direction. */
  42                 const double xsp;
  43                 /** The inverse box length in the y direction. */
  44                 const double ysp;
  45                 /** The number of boxes in the x direction. */
  46                 const int nx;
  47                 /** The number of boxes in the y direction. */
  48                 const int ny;
  49                 /** A constant, set to the value of nx multiplied by ny, which
  50                  * is used in the routines which step through boxes in
  51                  * sequence. */
  52                 const int nxy;
  53                 /** A boolean value that determines if the x coordinate in
  54                  * periodic or not. */
  55                 const bool xperiodic;
  56                 /** A boolean value that determines if the y coordinate in
  57                  * periodic or not. */
  58                 const bool yperiodic;
  59                 /** A boolean value that specifies whether the domain will be convex or non-convex.
  60                 * affected methods include cell_2d:plane, container_2d:import, container_2d:initialize_vo
  61                 ronoicell
  62                 * and container_2d:compute_cellsphere. */
  63                 const bool convex;
  64                 /** This array holds the number of particles within each
  65                  * computational box of the container. */
  66                 int *co;
  67
  68                 /** This array holds the maximum amount of particle memory for
  69                  * each computational box of the container. If the number of
  70                  * particles in a particular box ever approaches this limit,
  71                  * more is allocated using the add_particle_memory() function.
  72                  */
  73                 int *mem;
  74                 /** This array holds the numerical IDs of each particle in each
  75                  * computational box. */
  76                 int **id;
  77                 /** For each computational box, this holds the boundary walls that pass through that box.
  78                 */
  79                 int **wid;
  80                 /** A two dimensional array holding particle positions. For the
  81                  * derived container_poly class, this also holds particle
  82                  * radii. */
  83                 double **p;
  84                 /** A two dimensional array for holding the number of labels per
  85                  * particle, plus other status information. */
  86                 unsigned int **nlab;
  87                 /** A two dimensional array of pointers to label information.
  88                  */
  89                 int ***plab;
  90                 /** An array created from *tmp and is referenced my *soip. soi
  91                  * stands for "spheres of influence" */
  92                 int *soi;
  93                 /** The current number of boundary points. */
  94                 int noofbnds;
  95                 /** Specifies the memory allocation for the boundary points. */
  96                 int bnds_size;
  97                 /** This array holds the vertices of the non-convex domain. If
  98                  * the domain is convex this variable is not initialized. The
  99                  * format is (bnds[0],bnds[1])=(x1,y1) for some (x1,y1)
 100                  * boundary point and (bnds[2],bnds[3])=(x2,y2) is the next
 101                  * vertex in the boundary from (x1,y1) if we are going
 102                  * clockwise. */
 103                 double *bnds;
 104                 /** Information about how the boundary vertices are connected.
 105                  * A single integer is stored for each vertex, giving the index
 106                  * of the connecting vertex in the counter-clockwise sense. */
 107                 int *edb;
 108                 /** Contains information about which points are boundary points,
 109                 i.e. lie at a nonconvexity. **/
 110                 int **bndpts;
 111                 /** Contains information about which boundary points are problem points **/
 112                 bool *probpts;
 113                 bool *extpts;
 114                 container_2d(double xa,double xb,double ya,double yb,int xn,int yn,bool xper,bool yper,bool convex_,int memi);
 115                 ~container_2d();
 116                 void setup();
 117                 int crossproductz(double x1, double y1, double x2, double y2);
 118                 void semi_circle_labelling(double x1, double y1, double x2, double y2, int wid);
 119                 void add_edb_memory(int i);
 120                 void initial_cut(voronoicell_2d &c, double x, double y, int bid);
 121                 void initial_cut_nonconvex(voronoicell_2d &c, double x, double y, int bid);
 122                 void import(FILE *fp=stdin);
 123                 void container_projection(double x, double y, double &ccy, double &ccy, double &nx, double &ny);
 124                 /** Imports a list of particles from a file.
 125                  * \param[in] filename the file to read from. */
 126                 inline void import(const char *filename) {
 127                         FILE *fp(voropp_safe_fopen(filename,"r"));
 128                         import(fp);
 129                         fclose(fp);
 130                 }
 131                 void Perp_vector(double x1, double y1, double &x1, double &y1, int dcp);
 132                 void tag_walls(double x1, double y1, double x2, double y2, int wid);
 133                 inline void tag(int wallid, int box) {
 134                         int length=wid[box][0], *a;
 135                         if(length%6==0){
 136                                 a=new int[length+6];
 137                                 for(int j=0; j<length; j++){
 138                                         a[j]=wid[box][j];
 139                                 }
 140                                 delete [] wid[box]; wid[box]=a;
 141                         }
 142                         wid[box][length]=wallid;
 143                         wid[box][0]++;
 144                 }
 145                 void draw_particles(FILE *fp=stdout);
 146                 /** Dumps all the particle positions and IDs to a file.
 147                  * \param[in] filename the file to write to. */
 148                 inline void draw_particles(const char *filename) {
 149                         FILE *fp(voropp_safe_fopen(filename,"w"));
 150                         draw_particles(fp);
 151                         fclose(fp);
 152                 }
 153                 void draw_particles_pov(FILE *fp=stdout);
 154                 /** Dumps all the particles positions in POV-Ray format.
 155                  * \param[in] filename the file to write to. */
 156                 inline void draw_particles_pov(const char *filename) {
 157                         FILE *fp(voropp_safe_fopen(filename,"w"));
 158                         draw_particles_pov(fp);
 159                         fclose(fp);
 160                 }
 161                 void draw_cells_gnuplot(FILE *fp=stdout);
 162                 /** Computes the Voronoi cells for all particles and saves the
 163                  * output in gnuplot format.
 164                  * \param[in] filename the file to write to. */
 165                 inline void draw_cells_gnuplot(const char *filename) {
 166                         FILE *fp(voropp_safe_fopen(filename,"w"));
 167                         draw_cells_gnuplot(fp);
 168                         fclose(fp);
 169                 }
 170                 void draw_cells_pov(FILE *fp=stdout);
 171                 /** Computes the Voronoi cells for all particles and saves the
 172                  * output in POV-Ray format.
 173                  * \param[in] filename the file to write to. */
 174                 inline void draw_cells_pov(const char *filename) {
 175                         FILE *fp(voropp_safe_fopen(filename,"w"));
 176                         draw_cells_pov(fp);
 177                         fclose(fp);
 178                 }
 179                 void draw_boundary(FILE *fp=stdout);
 180                 inline void draw_boundary(const char *filename) {
 181                         FILE *fp(voropp_safe_fopen(filename,"w"));
 182                         draw_boundary(fp);
 183                         fclose(fp);
 184                 }
 185                 void print_custom(const char *format,FILE *fp=stdout);
 186                 /** Computes the Voronoi cells for all particles in the
 187                  * container, and for each cell, outputs a line containing
 188                  * custom information about the cell structure. The output
 189                  * format is specified using an input string with control
 190                  * sequences similar to the standard C printf() routine.
 191                  * \param[in] format the format of the output lines, using
 192                  *                   control sequences to denote the different
 193                  *                   cell statistics.
 194                  * \param[in] filename the file to write to. */
 195                 inline void print_custom(const char *format,const char *filename) {
 196                         FILE *fp(voropp_safe_fopen(filename,"w"));
 197                         print_custom(format,fp);
 198                         fclose(fp);
 199                 }
 200                 double sum_cell_areas();
 201                 void compute_all_cells();
 202                 void debug_output();
 203                 /** An overloaded version of the compute_cell_sphere routine,
 204                  * that sets up the x and y variables.
 205                  *\param[in,out] c a reference to a voronoicell object.
 206                  * \param[in] (i,j) the coordinates of the block that the test
 207                  *                  particle is in.
 208                  * \param[in] ij the index of the block that the test particle
 209                  *               is in, set to i+nx*j.
 210                  * \param[in] s the index of the particle within the test
 211                  *              block.
 212                  * \return False if the Voronoi cell was completely removed
 213                  * during the computation and has zero volume, true otherwise.
 214                  */
 215                 inline bool compute_cell_sphere(voronoicell_2d &c,int i,int j,int ij,int s) {
 216                         double x=p[s][2*ij],y=p[s][2*ij+1];
 217                         return compute_cell_sphere(c,i,j,ij,s,x,y);
 218                 }
 219                 bool compute_cell_sphere(voronoicell_2d &c,int i,int j,int ij,int s,double x,double y);
 220                 bool initialize_voronoicell(voronoicell_2d &c,double x,double y);
 221                 bool initialize_voronoicell_nonconvex(voronoicell_2d &c,double x, double y, int bid);
 222                 bool initialize_voronoicell_boundary(voronoicell_2d &c, double x, double y, int bid);
 223                 bool OKCuttingParticle(double gx, double gy, int gbox, int gindex, double cx, double cy,
 224                                         int cbox, int cindex, bool boundary, int bid);
 225                 void put(int n,double x,double y, int bndloc);
 226                 void clear();
 227         private:
 228                 /** A temporary array used in label creation. */
 229                 int *tmp;
 230                 /** The next free position in the temporary label array. */
 231                 int *tmpp;
 232                 /** The end of the temporary label array. */
 233                 int *tmpe;
 234                 inline bool put_locate_block(int &ij,double &x,double &y);
 235                 inline bool put_remap(int &ij,double &x,double &y);
 236                 inline double max(double a, double b) {return b<a?a:b;}
 237                 inline double min(double a, double b) {return (b<a)?b:a;}
 238                 inline double dist(double x1,double y1,double x2,double y2) {
 239                         return pow((((y2-y1)*(y2-y1))+((x2-x1)*(x2-x1))),.5);
 240                 }
 241                 void create_label_table();
 242                 void add_temporary_label_memory();
 243                 void add_particle_memory(int i);
 244                 void add_boundary_memory();
 245                 /** Custom int function, that gives consistent stepping for
 246                  * negative numbers. With normal int, we have
 247                  * (-1.5,-0.5,0.5,1.5) -> (-1,0,0,1). With this routine, we
 248                  * have (-1.5,-0.5,0.5,1.5) -> (-2,-1,0,1). */
 249                 inline int step_int(double a) {return a<0?int(a)-1:int(a);}
 250                 /** Custom modulo function, that gives consistent stepping for
 251                  * negative numbers. */
 252                 inline int step_mod(int a,int b) {return a>=0?a%b:b-1-(b-1-a)%b;}
 253                 /** Custom integer division function, that gives consistent
 254                  * stepping for negative numbers. */
 255                 inline int step_div(int a,int b) {return a>=0?a/b:-1+(a+1)/b;}
 256                 friend class voropp_loop_2d;
 257 };
 258
 259
 260 /** \brief A class to handle loops on regions of the container handling
 261  * non-periodic and periodic boundary conditions.
 262  *
 263  * Many of the container routines require scanning over a rectangular sub-grid
 264  * of blocks, and the routines for handling this are stored in the
 265  * voropp_loop_2d class. A voropp_loop_2d class can first be initialized to
 266  * either calculate the subgrid which is within a distance r of a vector
 267  * (vx,vy,vz), or a subgrid corresponding to a rectangular box. The routine
 268  * inc() can then be successively called to step through all the blocks within
 269  * this subgrid.
 270  */
 271 class voropp_loop_2d {
 272         public:
 273                 voropp_loop_2d(container_2d &con);
 274                 int init(double vx,double vy,double r,double &px,double &py);
 275                 int init(double xmin,double xmax,double ymin,double ymax,double &px,double &py);
 276                 int inc(double &px,double &py);
 277                 /** The current block index in the x direction, referencing a
 278                  * real cell in the range 0 to nx-1. */
 279                 int ip;
 280                 /** The current block index in the y direction, referencing a
 281                  * real cell in the range 0 to ny-1. */
 282                 int jp;
 283         private:
 284                 const double boxx,boxy,xsp,ysp,ax,ay;
 285                 const int nx,ny,nxy;
 286                 const bool xperiodic,yperiodic;
 287                 double apx,apy;
 288                 int i,j,ai,bi,aj,bj,s;
 289                 int aip,ajp,inc1;
 290                 /** Custom modulo function, that gives consistent stepping for
 291                  * negative numbers. */
 292                 inline int step_mod(int a,int b) {return a>=0?a%b:b-1-(b-1-a)%b;}
 293                 /** Custom integer division function, that gives consistent
 294                  * stepping for negative numbers. */
 295                 inline int step_div(int a,int b) {return a>=0?a/b:-1+(a+1)/b;}
 296                 /** Custom int function, that gives consistent stepping for
 297                  * negative numbers. With normal int, we have
 298                  * (-1.5,-0.5,0.5,1.5) -> (-1,0,0,1). With this routine, we
 299                  * have (-1.5,-0.5,0.5,1.5) -> (-2,-1,0,1). */
 300                 inline int step_int(double a) {return a<0?int(a)-1:int(a);}
 301 };
 302
 303 #endif