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