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. */ particle_order(int init_size=init_ordering_size)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. */ ~particle_order()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. */ add(int ijk,int q)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> c_loop_base(c_class & con)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. */ pos(double & x,double & y,double & z)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. */ pos(int & pid,double & x,double & y,double & z,double & r)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. */ x()155 inline double x() {return p[ijk][ps*q];} 156 /** Returns the y position of the particle currently being 157 * considered by the loop. */ y()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. */ z()161 inline double z() {return p[ijk][ps*q+2];} 162 /** Returns the ID of the particle currently being considered 163 * by the loop. */ pid()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> c_loop_all(c_class & con)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. */ start()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. */ inc()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. */ next_block()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> c_loop_subset(c_class & con)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. */ inc()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; step_mod(int a,int b)258 inline int step_mod(int a,int b) {return a>=0?a%b:b-1-(b-1-a)%b;} step_div(int a,int b)259 inline int step_div(int a,int b) {return a>=0?a/b:-1+(a+1)/b;} step_int(double a)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> c_loop_order(c_class & con,particle_order & vo_)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. */ start()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. */ inc()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. */ decode()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> c_loop_all_periodic(c_class & con)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. */ start()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. */ inc()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. */ next_block()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> c_loop_order_periodic(c_class & con,particle_order & vo_)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. */ start()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. */ inc()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. */ decode()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 463