branches/exact/src/c_loops.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.hh
   8  * \brief Header file for the loop classes. */
   9
  10 #ifndef VOROPP_C_LOOPS_HH
  11 #define VOROPP_C_LOOPS_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 class, determining what type of
  24  * geometrical region to loop over. */
  25 enum c_loop_subset_mode {
  26         sphere,
  27         box,
  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 {
  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                 /** The number of blocks in the z direction. */
  91                 const int nz;
  92                 /** A constant, set to the value of nx multiplied by ny, which
  93                  * is used in the routines that step through blocks in
  94                  * sequence. */
  95                 const int nxy;
  96                 /** A constant, set to the value of nx*ny*nz, which is used in
  97                  * the routines that step through blocks in sequence. */
  98                 const int nxyz;
  99                 /** The number of floating point numbers per particle in the
 100                  * associated container data structure. */
 101                 const int ps;
 102                 /** A pointer to the particle position information in the
 103                  * associated container data structure. */
 104                 double **p;
 105                 /** A pointer to the particle ID information in the associated
 106                  * container data structure. */
 107                 int **id;
 108                 /** A pointer to the particle counts in the associated
 109                  * container data structure. */
 110                 int *co;
 111                 /** The current x-index of the block under consideration by the
 112                  * loop. */
 113                 int i;
 114                 /** The current y-index of the block under consideration by the
 115                  * loop. */
 116                 int j;
 117                 /** The current z-index of the block under consideration by the
 118                  * loop. */
 119                 int k;
 120                 /** The current index of the block under consideration by the
 121                  * loop. */
 122                 int ijk;
 123                 /** The index of the particle under consideration within the current
 124                  * block. */
 125                 int q;
 126                 /** The constructor copies several necessary constants from the
 127                  * base container class.
 128                  * \param[in] con the container class to use. */
 129                 template<class c_class>
 130                 c_loop_base(c_class &con) : nx(con.nx), ny(con.ny), nz(con.nz),
 131                                             nxy(con.nxy), nxyz(con.nxyz), ps(con.ps),
 132                                             p(con.p), id(con.id), co(con.co) {}
 133                 /** Returns the position vector of the particle currently being
 134                  * considered by the loop.
 135                  * \param[out] (x,y,z) the position vector of the particle. */
 136                 inline void pos(double &x,double &y,double &z) {
 137                         double *pp=p[ijk]+ps*q;
 138                         x=*(pp++);y=*(pp++);z=*pp;
 139                 }
 140                 /** Returns the ID, position vector, and radius of the particle
 141                  * currently being considered by the loop.
 142                  * \param[out] pid the particle ID.
 143                  * \param[out] (x,y,z) the position vector of the particle.
 144                  * \param[out] r the radius of the particle. If no radius
 145                  *               information is available the default radius
 146                  *               value is returned. */
 147                 inline void pos(int &pid,double &x,double &y,double &z,double &r) {
 148                         pid=id[ijk][q];
 149                         double *pp=p[ijk]+ps*q;
 150                         x=*(pp++);y=*(pp++);z=*pp;
 151                         r=ps==3?default_radius:*(++pp);
 152                 }
 153                 /** Returns the x position of the particle currently being
 154                  * considered by the loop. */
 155                 inline double x() {return p[ijk][ps*q];}
 156                 /** Returns the y position of the particle currently being
 157                  * considered by the loop. */
 158                 inline double y() {return p[ijk][ps*q+1];}
 159                 /** Returns the z position of the particle currently being
 160                  * considered by the loop. */
 161                 inline double z() {return p[ijk][ps*q+2];}
 162                 /** Returns the ID of the particle currently being considered
 163                  * by the loop. */
 164                 inline int pid() {return id[ijk][q];}
 165 };
 166
 167 /** \brief Class for looping over all of the particles in a container.
 168  *
 169  * This is one of the simplest loop classes, that scans the computational
 170  * blocks in order, and scans all the particles within each block in order. */
 171 class c_loop_all : public c_loop_base {
 172         public:
 173                 /** The constructor copies several necessary constants from the
 174                  * base container class.
 175                  * \param[in] con the container class to use. */
 176                 template<class c_class>
 177                 c_loop_all(c_class &con) : c_loop_base(con) {}
 178                 /** Sets the class to consider the first particle.
 179                  * \return True if there is any particle to consider, false
 180                  * otherwise. */
 181                 inline bool start() {
 182                         i=j=k=ijk=q=0;
 183                         while(co[ijk]==0) if(!next_block()) return false;
 184                         return true;
 185                 }
 186                 /** Finds the next particle to test.
 187                  * \return True if there is another particle, false if no more
 188                  * particles are available. */
 189                 inline bool inc() {
 190                         q++;
 191                         if(q>=co[ijk]) {
 192                                 q=0;
 193                                 do {
 194                                         if(!next_block()) return false;
 195                                 } while(co[ijk]==0);
 196                         }
 197                         return true;
 198                 }
 199         private:
 200                 /** Updates the internal variables to find the next
 201                  * computational block with any particles.
 202                  * \return True if another block is found, false if there are
 203                  * no more blocks. */
 204                 inline bool next_block() {
 205                         ijk++;
 206                         i++;
 207                         if(i==nx) {
 208                                 i=0;j++;
 209                                 if(j==ny) {
 210                                         j=0;k++;
 211                                         if(ijk==nxyz) return false;
 212                                 }
 213                         }
 214                         return true;
 215                 }
 216 };
 217
 218 /** \brief Class for looping over a subset of particles in a container.
 219  *
 220  * This class can loop over a subset of particles in a certain geometrical
 221  * region within the container. The class can be set up to loop over a
 222  * rectangular box or sphere. It can also rectangular group of internal
 223  * computational blocks. */
 224 class c_loop_subset : public c_loop_base {
 225         public:
 226                 /** The current mode of operation, determining whether tests
 227                  * should be applied to particles to ensure they are within a
 228                  * certain geometrical object. */
 229                 c_loop_subset_mode mode;
 230                 /** The constructor copies several necessary constants from the
 231                  * base container class.
 232                  * \param[in] con the container class to use. */
 233                 template<class c_class>
 234                 c_loop_subset(c_class &con) : c_loop_base(con), ax(con.ax), ay(con.ay), az(con.az),
 235                         sx(con.bx-ax), sy(con.by-ay), sz(con.bz-az), xsp(con.xsp), ysp(con.ysp), zsp(con.zsp),
 236                         xperiodic(con.xperiodic), yperiodic(con.yperiodic), zperiodic(con.zperiodic) {}
 237                 void setup_sphere(double vx,double vy,double vz,double r,bool bounds_test=true);
 238                 void setup_box(double xmin,double xmax,double ymin,double ymax,double zmin,double zmax,bool bounds_test=true);
 239                 void setup_intbox(int ai_,int bi_,int aj_,int bj_,int ak_,int bk_);
 240                 bool start();
 241                 /** Finds the next particle to test.
 242                  * \return True if there is another particle, false if no more
 243                  * particles are available. */
 244                 inline bool inc() {
 245                         do {
 246                                 q++;
 247                                 while(q>=co[ijk]) {q=0;if(!next_block()) return false;}
 248                         } while(mode!=no_check&&out_of_bounds());
 249                         return true;
 250                 }
 251         private:
 252                 const double ax,ay,az,sx,sy,sz,xsp,ysp,zsp;
 253                 const bool xperiodic,yperiodic,zperiodic;
 254                 double px,py,pz,apx,apy,apz;
 255                 double v0,v1,v2,v3,v4,v5;
 256                 int ai,bi,aj,bj,ak,bk,s;
 257                 int ci,cj,ck,di,dj,dk,inc1,inc2;
 258                 inline int step_mod(int a,int b) {return a>=0?a%b:b-1-(b-1-a)%b;}
 259                 inline int step_div(int a,int b) {return a>=0?a/b:-1+(a+1)/b;}
 260                 inline int step_int(double a) {return a<0?int(a)-1:int(a);}
 261                 void setup_common();
 262                 bool next_block();
 263                 bool out_of_bounds();
 264 };
 265
 266 /** \brief Class for looping over all of the particles specified in a
 267  * pre-assembled particle_order class.
 268  *
 269  * The particle_order class can be used to create a specific order of particles
 270  * within the container. This class can then loop over these particles in this
 271  * order. The class is particularly useful in cases where the ordering of the
 272  * output must match the ordering of particles as they were inserted into the
 273  * container. */
 274 class c_loop_order : public c_loop_base {
 275         public:
 276                 /** A reference to the ordering class to use. */
 277                 particle_order &vo;
 278                 /** A pointer to the current position in the ordering class. */
 279                 int *cp;
 280                 /** A pointer to the end position in the ordering class. */
 281                 int *op;
 282                 /** The constructor copies several necessary constants from the
 283                  * base class, and sets up a reference to the ordering class to
 284                  * use.
 285                  * \param[in] con the container class to use.
 286                  * \param[in] vo_ the ordering class to use. */
 287                 template<class c_class>
 288                 c_loop_order(c_class &con,particle_order &vo_)
 289                 : c_loop_base(con), vo(vo_), nx(con.nx), nxy(con.nxy) {}
 290                 /** Sets the class to consider the first particle.
 291                  * \return True if there is any particle to consider, false
 292                  * otherwise. */
 293                 inline bool start() {
 294                         cp=vo.o;op=vo.op;
 295                         if(cp!=op) {
 296                                 ijk=*(cp++);decode();
 297                                 q=*(cp++);
 298                                 return true;
 299                         } else return false;
 300                 }
 301                 /** Finds the next particle to test.
 302                  * \return True if there is another particle, false if no more
 303                  * particles are available. */
 304                 inline bool inc() {
 305                         if(cp==op) return false;
 306                         ijk=*(cp++);decode();
 307                         q=*(cp++);
 308                         return true;
 309                 }
 310         private:
 311                 /** The number of computational blocks in the x direction. */
 312                 const int nx;
 313                 /** The number of computational blocks in a z-slice. */
 314                 const int nxy;
 315                 /** Takes the current block index and computes indices in the
 316                  * x, y, and z directions. */
 317                 inline void decode() {
 318                         k=ijk/nxy;
 319                         int ijkt=ijk-nxy*k;
 320                         j=ijkt/nx;
 321                         i=ijkt-j*nx;
 322                 }
 323 };
 324
 325 /** \brief A class for looping over all particles in a container_periodic or
 326  * container_periodic_poly class.
 327  *
 328  * Since the container_periodic and container_periodic_poly classes have a
 329  * fundamentally different memory organization, the regular loop classes cannot
 330  * be used with them. */
 331 class c_loop_all_periodic : public c_loop_base {
 332         public:
 333                 /** The constructor copies several necessary constants from the
 334                  * base periodic container class.
 335                  * \param[in] con the periodic container class to use. */
 336                 template<class c_class>
 337                 c_loop_all_periodic(c_class &con) : c_loop_base(con), ey(con.ey), ez(con.ez), wy(con.wy), wz(con.wz),
 338                         ijk0(nx*(ey+con.oy*ez)), inc2(2*nx*con.ey+1) {}
 339                 /** Sets the class to consider the first particle.
 340                  * \return True if there is any particle to consider, false
 341                  * otherwise. */
 342                 inline bool start() {
 343                         i=0;
 344                         j=ey;
 345                         k=ez;
 346                         ijk=ijk0;
 347                         q=0;
 348                         while(co[ijk]==0) if(!next_block()) return false;
 349                         return true;
 350                 }
 351                 /** Finds the next particle to test.
 352                  * \return True if there is another particle, false if no more
 353                  * particles are available. */
 354                 inline bool inc() {
 355                         q++;
 356                         if(q>=co[ijk]) {
 357                                 q=0;
 358                                 do {
 359                                         if(!next_block()) return false;
 360                                 } while(co[ijk]==0);
 361                         }
 362                         return true;
 363                 }
 364         private:
 365                 /** The lower y index (inclusive) of the primary domain within
 366                  * the block structure. */
 367                 int ey;
 368                 /** The lower y index (inclusive) of the primary domain within
 369                  * the block structure. */
 370                 int ez;
 371                 /** The upper y index (exclusive) of the primary domain within
 372                  * the block structure. */
 373                 int wy;
 374                 /** The upper z index (exclusive) of the primary domain within
 375                  * the block structure. */
 376                 int wz;
 377                 /** The index of the (0,0,0) block within the block structure.
 378                  */
 379                 int ijk0;
 380                 /** A value to increase ijk by when the z index is increased.
 381                  */
 382                 int inc2;
 383                 /** Updates the internal variables to find the next
 384                  * computational block with any particles.
 385                  * \return True if another block is found, false if there are
 386                  * no more blocks. */
 387                 inline bool next_block() {
 388                         i++;
 389                         if(i==nx) {
 390                                 i=0;j++;
 391                                 if(j==wy) {
 392                                         j=ey;k++;
 393                                         if(k==wz) return false;
 394                                         ijk+=inc2;
 395                                 } else ijk++;
 396                         } else ijk++;
 397                         return true;
 398                 }
 399 };
 400
 401 /** \brief Class for looping over all of the particles specified in a
 402  * pre-assembled particle_order class, for use with container_periodic classes.
 403  *
 404  * The particle_order class can be used to create a specific order of particles
 405  * within the container. This class can then loop over these particles in this
 406  * order. The class is particularly useful in cases where the ordering of the
 407  * output must match the ordering of particles as they were inserted into the
 408  * container. */
 409 class c_loop_order_periodic : public c_loop_base {
 410         public:
 411                 /** A reference to the ordering class to use. */
 412                 particle_order &vo;
 413                 /** A pointer to the current position in the ordering class. */
 414                 int *cp;
 415                 /** A pointer to the end position in the ordering class. */
 416                 int *op;
 417                 /** The constructor copies several necessary constants from the
 418                  * base class, and sets up a reference to the ordering class to
 419                  * use.
 420                  * \param[in] con the container class to use.
 421                  * \param[in] vo_ the ordering class to use. */
 422                 template<class c_class>
 423                 c_loop_order_periodic(c_class &con,particle_order &vo_)
 424                 : c_loop_base(con), vo(vo_), nx(con.nx), oxy(con.nx*con.oy) {}
 425                 /** Sets the class to consider the first particle.
 426                  * \return True if there is any particle to consider, false
 427                  * otherwise. */
 428                 inline bool start() {
 429                         cp=vo.o;op=vo.op;
 430                         if(cp!=op) {
 431                                 ijk=*(cp++);decode();
 432                                 q=*(cp++);
 433                                 return true;
 434                         } else return false;
 435                 }
 436                 /** Finds the next particle to test.
 437                  * \return True if there is another particle, false if no more
 438                  * particles are available. */
 439                 inline bool inc() {
 440                         if(cp==op) return false;
 441                         ijk=*(cp++);decode();
 442                         q=*(cp++);
 443                         return true;
 444                 }
 445         private:
 446                 /** The number of computational blocks in the x direction. */
 447                 const int nx;
 448                 /** The number of computational blocks in a z-slice. */
 449                 const int oxy;
 450                 /** Takes the current block index and computes indices in the
 451                  * x, y, and z directions. */
 452                 inline void decode() {
 453                         k=ijk/oxy;
 454                         int ijkt=ijk-oxy*k;
 455                         j=ijkt/nx;
 456                         i=ijkt-j*nx;
 457                 }
 458 };
 459
 460 }
 461
 462 #endif