1#ifndef HYPERCUBE_HPP
2#define HYPERCUBE_HPP
3
4#include "Grid/grid_sm.hpp"
5#include "Grid/comb.hpp"
6#include "util/mathutil.hpp"
7
8template<unsigned int dim, unsigned int subdim> class SubHyperCube;
9
10/*! \brief output stream overload for printing
11 *
12 * \param str ostream
13 * \param c combination to print
14 *
15 */
16template<unsigned int dim> std::ostream& operator<<(std::ostream& str, const comb<dim> & c)
17{
18 // print the combination of ostream
19 for (size_t i = 0 ; i < dim-1 ; i++)
20 {
21 /* coverity[dead_error_line] */
22 str << c.c[i] << ";";
23 }
24
25 str << c.c[dim-1];
26
27 return str;
28}
29
30/*! \brief This class calculate elements of the hyper-cube
31 *
32 * This class give you a set of utility functions for the hyper-cube like getting
33 * number of faces, number of edge, number of vertex, or in general number of elements
34 * of dimension d, position of each element
35 *
36 * * 0d Hyper-cube vertex
37 * * 1d Hypercube segment
38 * * 2d Hypercube square
39 * * 3d Hypercube Cube
40 * ...
41 *
42 * \tparam dim dimensionality of the Hyper-cube
43 *
44 * ### Get vertex and edge on a line
45 * \snippet HyperCube_unit_test.hpp Get vertex and edge on a line
46 * ### Get vertex edge and surfaces of a square
47 * \snippet HyperCube_unit_test.hpp Get vertex edge and surfaces of a square
48 * ### Get vertex edge surfaces and volumes of a cube
49 * \snippet HyperCube_unit_test.hpp Get vertex edge surfaces and volumes of a cube
50 *
51 * hyper-cube define only the features of an N-dimensional hyper-cube, does not define
52 * where is is located and its size, use Box for that purpose
53 *
54 */
55
56template<unsigned int dim>
57class HyperCube
58{
59public:
60
61 /*! \brief Get the number of Elements of dimension d
62 *
63 * \param d dimensionality of the element
64 * \return the number of elements of dimension d
65 *
66 */
67 static size_t getNumberOfElements_R(size_t d)
68 {
69 // The formula to calculate the number of element of size d is given by
70 //
71 // 2^(dim-d) * C(dim,d) with C(dim,d) = dim!/(d!(dim-d)!)
72
73 return pow(2,dim-d) * openfpm::math::C(dim,d);
74 }
75
76 /*! \brief Get the sum of the number of elements from d to d_t (included)
77 *
78 * \param d_t
79 * \return the sum of the number of elements from d to d_t
80 *
81 */
82
83 static size_t getNumberOfElementsTo_R(size_t d_t)
84 {
85 size_t N_ele = 0;
86
87 for (size_t i = 0 ; i <= d_t ; i++)
88 N_ele += getNumberOfElements_R(i);
89
90 return N_ele;
91 }
92
93 /*! brief Calculate the position (combinations) of all the elements of size d
94 *
95 * \param d dimensionality of the object returned in the combinations
96 * \return all the combinations
97 *
98 */
99
100 static std::vector<comb<dim>> getCombinations_R(size_t d)
101 {
102#ifdef SE_CLASS1
103 if (d > dim)
104 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " " << d << " must be smaller than " << "\n";
105#endif
106
107 // Create an Iterator_g_const
108 // And a vector that store all the combination
109
110 std::vector<comb<dim>> v;
111 Iterator_g_const it(dim-d,dim);
112
113 // for each non-zero elements
114 // basically the key will store the position of the
115 // non zero elements, while BinPermutations will
116 // fill the array of all the permutations
117 //
118
119 while (it.isNext())
120 {
121 grid_key_dx_r & key = it.get();
122
123 // Calculate the permutation
124
125 BinPermutations(key,v);
126 ++it;
127 }
128
129 // case when d == dim
130
131 if (d == dim)
132 {
133 comb<dim> c;
134 c.zero();
135
136 v.push_back(c);
137 }
138
139 // return the combinations
140
141 return v;
142 }
143
144 /*! brief Calculate the position (combinations) of all the elements of size d
145 *
146 * \param d dimensionality of the object returned in the combinations
147 * \return all the combinations
148 *
149 */
150
151 static std::vector<comb<dim>> getCombinations_R_bc(size_t d, const size_t (& bc)[dim])
152 {
153#ifdef SE_CLASS1
154 if (d > dim)
155 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " " << d << " must be smaller than " << "\n";
156#endif
157
158 size_t dg[dim];
159
160 size_t k = 0;
161 // get the indexes of the free degree of freedom
162 for (size_t i = 0 ; i < dim ; i++)
163 {
164 if (bc[i] == PERIODIC)
165 {
166 dg[k] = i;
167 k++;
168 }
169 }
170
171 // Create an Iterator_g_const
172 // And a vector that store all the combination
173
174 std::vector<comb<dim>> v;
175 Iterator_g_const it(k-(d-(dim - k)),k);
176
177 // for each non-zero elements
178 // basically the key will store the position of the
179 // non zero elements, while BinPermutations will
180 // fill the array of all the permutations
181 //
182
183 while (it.isNext())
184 {
185 grid_key_dx_r key = it.get();
186
187 // Now we adjust the non zero position
188
189 for (size_t i = 0 ; i < key.getDim() ; i++)
190 {key.set_d(i,dg[key.get(i)]);}
191
192 // Calculate the permutation
193
194 BinPermutations(key,v);
195 ++it;
196 }
197
198 // case when d == dim
199
200 if (d == dim)
201 {
202 comb<dim> c;
203 c.zero();
204
205 v.push_back(c);
206 }
207
208 // return the combinations
209
210 return v;
211 }
212
213 /*! \brief Binary permutations
214 *
215 * Fill v with all the possible binary permutations
216 * it produce 2^(pos.getDim()) Permutations
217 *
218 * Example
219 *
220 * if getDim() is 2
221 *
222 * it produce 4 configuration
223 *
224 * (1,1) (1,-1) (-1,1) (-1,-1)
225 *
226 * and fill the number in the position indicated by Iterator_g_const
227 *
228 * \param pos slots inside comb to fill with all permutations
229 * \param v vector to fill with the permutations
230 *
231 */
232 static void BinPermutations(grid_key_dx_r & pos, std::vector<comb<dim>> & v)
233 {
234 size_t p_stop = pow(2,pos.getDim());
235 for (size_t p = 0 ; p < p_stop ; p++)
236 {
237 // Create a new permutation
238 struct comb<dim> c;
239
240 // zero the combination
241 c.zero();
242
243 // the bit of p give the permutations 0 mean bottom 1 mean up
244 // Fill the combination based on the bit of p
245
246 for (size_t k = 0 ; k < pos.getDim() ; k++)
247 {
248 // bit of p
249 bool bit = (p >> k) & 0x01;
250
251 // Fill the combination
252 c.c[pos.get(k)] = (bit == 0)? 1 : -1;
253 }
254
255 // save the combination
256
257 v.push_back(c);
258 }
259 }
260
261 /*! \brief Binary permutations
262 *
263 * Fill v with all the possible binary permutations
264 * it produce 2^(pos.getDim()) Permutations
265 *
266 * Example
267 *
268 * if getDim() is 2
269 *
270 * it produce 4 configuration
271 *
272 * (1,1) (1,-0) (0,1) (0,0)
273 *
274 * from another prospective given
275 *
276 * \verbatim
277 *
278 * +----#----+
279 * | |
280 * | |
281 * # * #
282 * | |
283 * | |
284 * +----#----+
285 *
286 * \endverbatim
287 *
288 * combination in the center (*) the down-left vertex (+). down and left edge (#)
289 *
290 * \param v vector to fill with the permutations
291 *
292 */
293 static void BinPermutationsSt(std::vector<comb<dim>> & v)
294 {
295 size_t p_stop = pow(2,dim);
296 for (size_t p = 0 ; p < p_stop ; p++)
297 {
298 // Create a new permutation
299 struct comb<dim> c;
300
301 // zero the combination
302 c.zero();
303
304 // the bit of p give the permutations 0 mean bottom 1 mean up
305 // Fill the combination based on the bit of p
306
307 for (size_t k = 0 ; k < dim ; k++)
308 {
309 // bit of p
310 bool bit = (p >> k) & 0x01;
311
312 // Fill the combination
313 c.c[k] = (bit == 0)? 0 : -1;
314 }
315
316 // save the combination
317
318 v.push_back(c);
319 }
320 }
321
322 /*! \brief Linearize the permutation given by BinPermutationSt
323 *
324 * Suppose BinPermutation return the following combination
325 *
326 * (-1,-1) (-1,0) (0,-1) (0,0)
327 *
328 * giving (0,-1) it return 2
329 *
330 * \param c combination to linearize
331 *
332 * \return the linearized permutation
333 *
334 * \see BinPermitationSt
335 *
336 */
337 static size_t LinPerm(comb<dim> & c)
338 {
339 size_t id = 0;
340
341 for (size_t i = 0 ; i < dim ; i++)
342 {
343 if (c.c[i] == -1)
344 {id = id | (1 << i);}
345 }
346
347 return id;
348 }
349
350 /*
351 static SubHyperCube<dim,dim-1> getSubHyperCube(int d)
352 {
353 SubHyperCube<dim,dim-1> sub;
354
355 return sub;
356 }*/
357
358 /** \brief Linearize the combination
359 *
360 * It map the combination into a linear id, in particular given the vector of combinations
361 * with get getCombinations_R, given the combination it give where is located in the vector
362 *
363 * \param c given combination
364 * \return the linearized combination
365 *
366 */
367 static size_t LinId(comb<dim> & c)
368 {
369 // calculate the numbers of non-zero
370 size_t d = 0;
371
372 size_t pos_n_zero[dim];
373
374 for (size_t i = 0 ; i < dim ; i++)
375 {
376 if (c.c[i] != 0)
377 {d++;}
378 }
379
380 // Get the position of the non-zero
381 size_t pn_zero = 0;
382 for (size_t i = 0 ; i < dim ; i++)
383 {
384 if (c.c[i] != 0)
385 {
386 pos_n_zero[d-pn_zero-1] = i;
387 pn_zero++;
388 }
389 }
390
391 size_t lin_id = 0;
392
393 // Cumulative value
394 long int val = 0;
395 long int cum_val = 0;
396 for(long int i = d - 1; i >= 0 ; i--)
397 {
398 // check the out-of-bound, outside is assumed to be -1 so (- pos_n_zero[i+1] - 1) = 0
399 if (i+1 < (long int)d)
400 {
401 /* coverity[dead_error_line] */
402 val = pos_n_zero[i] - pos_n_zero[i+1] - 1;
403 }
404 else
405 {
406 /* coverty[uninit_use] */
407 val = pos_n_zero[i];
408 }
409
410 for (long int j = 0 ; j < (long int)val; j++)
411 {
412 // C is not safe, check the limit
413 /* coverity[dead_error_line] */
414 if (((long int)dim)-cum_val-j-1 >= 0 && i > 0 && ((long int)dim)-cum_val-j >= i)
415 lin_id += openfpm::math::C(dim-cum_val-j-1,i);
416 else
417 lin_id += 1;
418 }
419
420 cum_val += (val + 1);
421 }
422
423 // multiply for the permutation
424 lin_id *= pow(2,d);
425
426 // calculate the permutation position
427 size_t id = 0;
428
429 for (size_t i = 0 ; i < d ; i++)
430 {
431 if (c.c[pos_n_zero[i]] == -1)
432 {id = id | (1 << i);}
433 }
434
435 // return the correct id
436
437 return lin_id + id;
438 }
439
440 /** \brief isPositive return if the combination d is a positive or a negative
441 *
442 * For an hyper-cube of dimension dim we have 2*dim faces combinations half on positive direction
443 * half on negative direction, the function check
444 * if the d combination is negative or positive
445 *
446 * \param d
447 *
448 * \return true if the combination is in positive direction
449 *
450 */
451 static bool isPositive(size_t d)
452 {
453 return (d % 2) == 0;
454 }
455
456 /*! \brief return the combination of the positive face on direction d
457 *
458 * \param d direction
459 *
460 * \return id of the combination
461 *
462 */
463 static int positiveFace(int d)
464 {
465 return d * 2;
466 }
467
468 /*! \brief Return the combination of the negative face on direction d
469 *
470 * \param d direction
471 *
472 * \return id of the combination
473 *
474 */
475 static int negativeFace(int d)
476 {
477 return d * 2 + 1;
478 }
479};
480
481/*! \brief This represent a sub-hyper-cube of an hyper-cube like a face or an edge of a cube
482 *
483 * It give a set of utility function to work with sub-hyper-cubes like the hyper-cube
484 *
485 * \tparam dimensionality of the hyper-cube
486 * \tparam dimensionality of the sub-hyper-cube
487 *
488 * ### First we get the surfaces of an hyper-cube
489 * \snippet HyperCube_unit_test.hpp Getting the surfaces of the cube
490 * ### Getting the vertices of one particular surface of the cube
491 * \snippet HyperCube_unit_test.hpp Getting the vertices of one surface of the cube
492 * ### Getting the edges of one particular surface of the cube
493 * \snippet HyperCube_unit_test.hpp Getting the edges of the surfaces of the cube
494 *
495 */
496
497template<unsigned int dim, unsigned int subdim>
498class SubHyperCube : public HyperCube<subdim>
499{
500public:
501
502 /*! brief Calculate the position (combinations) of all the elements of size d in the sub-hyper-cube
503 *
504 * \param c identify the position of the sub-hyper-cube in the hypercube
505 * \param d dimensionality of the objects
506 * \return all the combinations
507 *
508 */
509 static std::vector<comb<dim>> getCombinations_R(comb<dim> c, int d)
510 {
511#ifdef DEBUG
512 if (c.n_zero() < d)
513 {
514 std::cerr << "Error SubHyperCube: " << __FILE__ << " " << __LINE__ << " the hyper-cube selected must have dimensionality bigger than the dimensionality of the requested combinations, or the number of zero in c must me bigger than d" << "\n";
515 }
516#endif
517
518 // if sub-dim == 0 return c
519
520 if (subdim == 0)
521 {
522 std::vector<comb<dim>> vc;
523 vc.push_back(c);
524
525 return vc;
526 }
527
528 // Create an Iterator_g_const
529 // And a vector that store all the combination
530
531 std::vector<comb<subdim>> v = HyperCube<subdim>::getCombinations_R(d);
532
533 // Create new combinations
534 std::vector<comb<dim>> vc(v.size());
535
536 // for each combination
537 for (size_t i = 0 ; i < v.size() ; i++)
538 {
539 // sub j counter
540 int sub_j = 0;
541 // for each zero (direction spanned by the sub-hyper-cube)
542 for (size_t j = 0 ; j < dim ; j++)
543 {
544 if (c.c[j] == 0)
545 {
546 // take the combination from the sub-hyper-cube
547 vc[i].c[j] = v[i].c[sub_j];
548 sub_j++;
549 }
550 else
551 {
552 // take the combination from the hyper-cube position
553 vc[i].c[j] = c.c[j];
554 }
555 }
556 }
557
558 // return the combinations
559 return vc;
560 }
561};
562
563#endif
564