branches/2d/src/c_loops_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 c_loops_2d.hh
   8  * \brief Header file for the 2D loop classes. */
   9
  10 #ifndef VOROPP_C_LOOPS_2D_HH
  11 #define VOROPP_C_LOOPS_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
  21 namespace voro {
  22
  23 /** A type associated with a c_loop_subset_2d class, determining what type of
  24  * geometrical region to loop over. */
  25 enum c_loop_subset_mode_2d {
  26         circle,
  27         rectangle,
  28         no_check
  29 };
  30
  31 /** \brief A class for storing ordering information when particles are added to
  32  * a container.
  33  *
  34  * When particles are added to a container class, they are sorted into an
  35  * internal computational grid of blocks. The particle_order class provides a
  36  * mechanism for remembering which block particles were sorted into. The import
  37  * and put routines in the container class have variants that also take a
  38  * particle_order class. Each time they are called, they will store the block
  39  * that the particle was sorted into, plus the position of the particle within
  40  * the block. The particle_order class can used by the c_loop_order class to
  41  * specifically loop over the particles that have their information stored
  42  * within it. */
  43 class particle_order {
  44         public:
  45                 /** A pointer to the array holding the ordering. */
  46                 int *o;
  47                 /** A pointer to the next position in the ordering array in
  48                  * which to store an entry. */
  49                 int *op;
  50                 /** The current memory allocation for the class, set to the
  51                  * number of entries which can be stored. */
  52                 int size;
  53                 /** The particle_order constructor allocates memory to store the
  54                  * ordering information.
  55                  * \param[in] init_size the initial amount of memory to
  56                  *                      allocate. */
  57                 particle_order(int init_size=init_ordering_size)
  58                         : o(new int[init_size<<1]),op(o),size(init_size) {}
  59                 /** The particle_order destructor frees the dynamically allocated
  60                  * memory used to store the ordering information. */
  61                 ~particle_order() {
  62                         delete [] o;
  63                 }
  64                 /** Adds a record to the order, corresponding to the memory
  65                  * address of where a particle was placed into the container.
  66                  * \param[in] ijk the block into which the particle was placed.
  67                  * \param[in] q the position within the block where the
  68                  *              particle was placed. */
  69                 inline void add(int ijk,int q) {
  70                         if(op==o+size) add_ordering_memory();
  71                         *(op++)=ijk;*(op++)=q;
  72                 }
  73         private:
  74                 void add_ordering_memory();
  75 };
  76
  77 /** \brief Base class for looping over particles in a container.
  78  *
  79  * This class forms the base of all classes that can loop over a subset of
  80  * particles in a contaner in some order. When initialized, it stores constants
  81  * about the corresponding container geometry. It also contains a number of
  82  * routines for interrogating which particle currently being considered by the
  83  * loop, which are common between all of the derived classes. */
  84 class c_loop_base_2d {
  85         public:
  86                 /** The number of blocks in the x direction. */
  87                 const int nx;
  88                 /** The number of blocks in the y direction. */
  89                 const int ny;
  90                 /** A constant, set to the value of nx multiplied by ny, which
  91                  * is used in the routines that step through blocks in
  92                  * sequence. */
  93                 const int nxy;
  94                 /** The number of floating point numbers per particle in the
  95                  * associated container data structure. */
  96                 const int ps;
  97                 /** A pointer to the particle position information in the
  98                  * associated container data structure. */
  99                 double **p;
 100                 /** A pointer to the particle ID information in the associated
 101                  * container data structure. */
 102                 int **id;
 103                 /** A pointer to the particle counts in the associated
 104                  * container data structure. */
 105                 int *co;
 106                 /** The current x-index of the block under consideration by the
 107                  * loop. */
 108                 int i;
 109                 /** The current y-index of the block under consideration by the
 110                  * loop. */
 111                 int j;
 112                 /** The current index of the block under consideration by the
 113                  * loop. */
 114                 int ij;
 115                 /** The index of the particle under consideration within the current
 116                  * block. */
 117                 int q;
 118                 /** The constructor copies several necessary constants from the
 119                  * base container class.
 120                  * \param[in] con the container class to use. */
 121                 template<class c_class_2d>
 122                 c_loop_base_2d(c_class_2d &con) : nx(con.nx), ny(con.ny), nxy(con.nxy),
 123                                                   ps(con.ps), p(con.p), id(con.id),
 124                                                   co(con.co) {}
 125                 /** Returns the position vector of the particle currently being
 126                  * considered by the loop.
 127                  * \param[out] (x,y) the position vector of the particle. */
 128                 inline void pos(double &x,double &y) {
 129                         double *pp=p[ij]+ps*q;
 130                         x=*(pp++);y=*pp;
 131                 }
 132                 /** Returns the ID, position vector, and radius of the particle
 133                  * currently being considered by the loop.
 134                  * \param[out] pid the particle ID.
 135                  * \param[out] (x,y) the position vector of the particle.
 136                  * \param[out] r the radius of the particle. If no radius
 137                  *               information is available the default radius
 138                  *               value is returned. */
 139                 inline void pos(int &pid,double &x,double &y,double &r) {
 140                         pid=id[ij][q];
 141                         double *pp=p[ij]+ps*q;
 142                         x=*(pp++);y=*pp;
 143                         r=ps==2?default_radius_2d:*(++pp);
 144                 }
 145                 /** Returns the x position of the particle currently being
 146                  * considered by the loop. */
 147                 inline double x() {return p[ij][ps*q];}
 148                 /** Returns the y position of the particle currently being
 149                  * considered by the loop. */
 150                 inline double y() {return p[ij][ps*q+1];}
 151                 /** Returns the ID of the particle currently being considered
 152                  * by the loop. */
 153                 inline int pid() {return id[ij][q];}
 154 };
 155
 156 /** \brief Class for looping over all of the particles in a container.
 157  *
 158  * This is one of the simplest loop classes, that scans the computational
 159  * blocks in order, and scans all the particles within each block in order. */
 160 class c_loop_all_2d : public c_loop_base_2d {
 161         public:
 162                 /** The constructor copies several necessary constants from the
 163                  * base container class.
 164                  * \param[in] con the container class to use. */
 165                 template<class c_class_2d>
 166                 c_loop_all_2d(c_class_2d &con) : c_loop_base_2d(con) {}
 167                 /** Sets the class to consider the first particle.
 168                  * \return True if there is any particle to consider, false
 169                  * otherwise. */
 170                 inline bool start() {
 171                         i=j=ij=q=0;
 172                         while(co[ij]==0) if(!next_block()) return false;
 173                         return true;
 174                 }
 175                 /** Finds the next particle to test.
 176                  * \return True if there is another particle, false if no more
 177                  * particles are available. */
 178                 inline bool inc() {
 179                         q++;
 180                         if(q>=co[ij]) {
 181                                 q=0;
 182                                 do {
 183                                         if(!next_block()) return false;
 184                                 } while(co[ij]==0);
 185                         }
 186                         return true;
 187                 }
 188         private:
 189                 /** Updates the internal variables to find the next
 190                  * computational block with any particles.
 191                  * \return True if another block is found, false if there are
 192                  * no more blocks. */
 193                 inline bool next_block() {
 194                         ij++;
 195                         i++;
 196                         if(i==nx) {
 197                                 i=0;j++;
 198                                 if(j==ny) return false;
 199                         }
 200                         return true;
 201                 }
 202 };
 203
 204 /** \brief Class for looping over a subset of particles in a container.
 205  *
 206  * This class can loop over a subset of particles in a certain geometrical
 207  * region within the container. The class can be set up to loop over a
 208  * rectangle or circle. It can also rectangular group of internal computational
 209  * blocks. */
 210 class c_loop_subset_2d : public c_loop_base_2d {
 211         public:
 212                 /** The current mode of operation, determining whether tests
 213                  * should be applied to particles to ensure they are within a
 214                  * certain geometrical object. */
 215                 c_loop_subset_mode_2d mode;
 216                 /** The constructor copies several necessary constants from the
 217                  * base container class.
 218                  * \param[in] con the container class to use. */
 219                 template<class c_class_2d>
 220                 c_loop_subset_2d(c_class_2d &con) : c_loop_base_2d(con), ax(con.ax), ay(con.ay),
 221                         sx(con.bx-ax), sy(con.by-ay), xsp(con.xsp), ysp(con.ysp),
 222                         xperiodic(con.xperiodic), yperiodic(con.yperiodic) {}
 223                 void setup_circle(double vx,double vy,double r,bool bounds_test=true);
 224                 void setup_box(double xmin,double xmax,double ymin,double ymax,bool bounds_test=true);
 225                 void setup_intbox(int ai_,int bi_,int aj_,int bj_);
 226                 bool start();
 227                 /** Finds the next particle to test.
 228                  * \return True if there is another particle, false if no more
 229                  * particles are available. */
 230                 inline bool inc() {
 231                         do {
 232                                 q++;
 233                                 while(q>=co[ij]) {q=0;if(!next_block()) return false;}
 234                         } while(mode!=no_check&&out_of_bounds());
 235                         return true;
 236                 }
 237         private:
 238                 const double ax,ay,sx,sy,xsp,ysp;
 239                 const bool xperiodic,yperiodic;
 240                 double px,py,apx,apy;
 241                 double v0,v1,v2,v3;
 242                 int ai,bi,aj,bj;
 243                 int ci,cj,di,dj,inc1;
 244                 inline int step_mod(int a,int b) {return a>=0?a%b:b-1-(b-1-a)%b;}
 245                 inline int step_div(int a,int b) {return a>=0?a/b:-1+(a+1)/b;}
 246                 inline int step_int(double a) {return a<0?int(a)-1:int(a);}
 247                 void setup_common();
 248                 bool next_block();
 249                 bool out_of_bounds();
 250 };
 251
 252 /** \brief Class for looping over all of the particles specified in a
 253  * pre-assembled particle_order class.
 254  *
 255  * The particle_order class can be used to create a specific order of particles
 256  * within the container. This class can then loop over these particles in this
 257  * order. The class is particularly useful in cases where the ordering of the
 258  * output must match the ordering of particles as they were inserted into the
 259  * container. */
 260 class c_loop_order_2d : public c_loop_base_2d {
 261         public:
 262                 /** A reference to the ordering class to use. */
 263                 particle_order &vo;
 264                 /** A pointer to the current position in the ordering class. */
 265                 int *cp;
 266                 /** A pointer to the end position in the ordering class. */
 267                 int *op;
 268                 /** The constructor copies several necessary constants from the
 269                  * base class, and sets up a reference to the ordering class to
 270                  * use.
 271                  * \param[in] con the container class to use.
 272                  * \param[in] vo_ the ordering class to use. */
 273                 template<class c_class_2d>
 274                 c_loop_order_2d(c_class_2d &con,particle_order &vo_)
 275                 : c_loop_base_2d(con), vo(vo_), nx(con.nx) {}
 276                 /** Sets the class to consider the first particle.
 277                  * \return True if there is any particle to consider, false
 278                  * otherwise. */
 279                 inline bool start() {
 280                         cp=vo.o;op=vo.op;
 281                         if(cp!=op) {
 282                                 ij=*(cp++);decode();
 283                                 q=*(cp++);
 284                                 return true;
 285                         } else return false;
 286                 }
 287                 /** Finds the next particle to test.
 288                  * \return True if there is another particle, false if no more
 289                  * particles are available. */
 290                 inline bool inc() {
 291                         if(cp==op) return false;
 292                         ij=*(cp++);decode();
 293                         q=*(cp++);
 294                         return true;
 295                 }
 296         private:
 297                 /** The number of computational blocks in the x direction. */
 298                 const int nx;
 299                 /** Takes the current block index and computes indices in the
 300                  * x, y directions. */
 301                 inline void decode() {
 302                         j=ij/nx;
 303                         i=ij-j*nx;
 304                 }
 305 };
 306
 307 }
 308
 309 #endif