8 #ifndef CELLDECOMPOSER_HPP_
9 #define CELLDECOMPOSER_HPP_
11 #include "Space/SpaceBox.hpp"
12 #include "Space/Matrix.hpp"
13 #include "util/copy_compare/meta_compare.hpp"
14 #include "Grid/grid_sm.hpp"
16 #define CELL_DECOMPOSER 8001lu
25 template<
unsigned int dim,
typename T>
56 inline T
transform(
const T(&s)[dim],
const size_t i)
const
58 return s[i] -
sh.get(i);
71 return s.
get(i) -
sh.get(i);
84 return s.template get<0>()[i] -
sh.get(i);
95 for (
size_t i = 0 ; i < dim ; i++)
96 sh.get(i) = orig.
get(i);
147 template<
unsigned int dim,
typename T>
178 inline T
transform(
const T(&s)[dim],
const size_t i)
const
206 return s.template get<Point<dim,T>::x>()[i];
310 template<
unsigned int dim,
typename T,
typename transform = no_transform<dim,T>>
311 class CellDecomposer_sm
341 inline size_t ConvertToID(
const T (&x)[dim] ,
size_t s)
const
343 size_t id = openfpm::math::size_t_floor(t.transform(x,s) / box_unit.getHigh(s)) + off[s];
344 id = (
id >= gr_cell.size(s))?(gr_cell.size(s)-1-cell_shift.get(s)):
id-cell_shift.get(s);
354 inline size_t ConvertToID(
const Point<dim,T> & x ,
size_t s,
size_t sc = 0)
const
356 size_t id = openfpm::math::size_t_floor(t.transform(x,s) / box_unit.getHigh(s)) + off[s];
357 id = (
id >= gr_cell.size(s))?(gr_cell.size(s)-1-cell_shift.get(s)):
id-cell_shift.get(s);
367 inline size_t ConvertToID_ns(
const T (&x)[dim] ,
size_t s)
const
369 size_t id = openfpm::math::size_t_floor(t.transform(x,s) / box_unit.getHigh(s)) + off[s];
370 id = (
id >= gr_cell.size(s))?(gr_cell.size(s)-1):
id;
380 inline size_t ConvertToID_ns(
const Point<dim,T> & x ,
size_t s,
size_t sc = 0)
const
382 size_t id = openfpm::math::size_t_floor(t.transform(x,s) / box_unit.getHigh(s)) + off[s];
383 id = (
id >= gr_cell.size(s))?(gr_cell.size(s)-1):
id;
393 template <
typename Mem>
inline size_t ConvertToID_(
const encapc<1,
Point<dim,T>,Mem> & x ,
size_t s,
size_t sc = 0)
const
395 size_t id = (size_t)(t.transform(x,s) / box_unit.getHigh(s)) + off[s];
396 id = (
id >= gr_cell.size(s))?(gr_cell.size(s)-1-cell_shift.get(s)):
id-cell_shift.get(s);
406 template<
typename Ele>
inline void check_and_print_error(
const Ele & pos ,
size_t s)
const
411 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" using an uninitialized CellDecomposer";
412 ACTION_ON_ERROR(CELL_DECOMPOSER);
415 if (pos[0] < box.getLow(s) - off[s]*box_unit.getP2()[s] || pos[s] > box.getHigh(s) + off[s]*box_unit.getP2()[s])
417 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" point " <<
Point<dim,T>(pos).toString() <<
" is not inside the cell space";
418 ACTION_ON_ERROR(CELL_DECOMPOSER);
424 template<
typename Ele>
inline size_t getCellDom_impl(
const Ele & pos)
const
426 check_and_print_error(pos,0);
428 size_t cell_id = ConvertToID_ns(pos,0);
430 cell_id = (cell_id == gr_cell.size(0) - off[0])?gr_cell.size(0) - off[0] - 1:cell_id;
431 cell_id = (cell_id == off[0]-1)?off[0]:cell_id;
433 cell_id -= cell_shift.get(0);
435 for (
size_t s = 1 ; s < dim ; s++)
438 check_and_print_error(pos,s);
440 size_t cell_idt = ConvertToID_ns(pos,s);
442 cell_idt = (cell_idt == gr_cell.size(s) - off[s])?gr_cell.size(s) - off[s] - 1:cell_idt;
443 cell_idt = (cell_idt == off[s]-1)?off[s]:cell_idt;
445 cell_idt -= cell_shift.get(s);
447 cell_id += gr_cell2.size_s(s-1) * cell_idt;
453 template<
typename Ele>
inline size_t getCellPad_impl(
const Ele & pos)
const
455 check_and_print_error(pos,0);
457 size_t cell_id = ConvertToID_ns(pos,0);
458 cell_id = (cell_id == off[0])?off[0]-1:cell_id;
459 cell_id = (cell_id == gr_cell.size(0) - off[0] - 1)?gr_cell.size(0) - off[0]:cell_id;
461 cell_id -= cell_shift.get(0);
463 for (
size_t s = 1 ; s < dim ; s++)
465 check_and_print_error(pos,s);
467 size_t cell_idt = ConvertToID_ns(pos,s);
468 cell_idt = (cell_idt == off[s])?off[s]-1:cell_idt;
469 cell_idt = (cell_idt == gr_cell.size(s) - off[s] - 1)?gr_cell.size(s) - off[s]:cell_idt;
471 cell_idt -= cell_shift.get(s);
473 cell_id += gr_cell2.size_s(s-1) * cell_idt;
507 mutable size_t div_wp[dim];
512 void Initialize(
const size_t pad,
const size_t (& div)[dim])
516 for (
size_t i = 0 ; i < dim ; i++)
520 std::cerr <<
"Error " << __FILE__ <<
":" << __LINE__ <<
" the number of cells on each dimension must be different from zero\n";
521 ACTION_ON_ERROR(CELL_DECOMPOSER)
530 for (
size_t i = 0 ; i < dim ; i++)
531 div_p[i] = div[i] + 2*pad;
533 gr_cell.setDimensions(div_p);
534 gr_cell2.setDimensions(div_p);
540 for (
size_t i = 0 ; i < dim ; i++)
542 tot_n_cell *= gr_cell.size(i);
545 box_unit.setHigh(i,(box.getHigh(i) - box.getLow(i)) / (gr_cell.size(i)- 2*pad) );
548 for (
size_t i = 0; i < dim ; i++)
553 p_middle = box_unit.getP2();
554 p_middle = p_middle / 2;
568 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" using an uninitialized CellDecomposer";
588 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" using an uninitialized CellDecomposer";
589 ACTION_ON_ERROR(CELL_DECOMPOSER);
594 key.
set_d(0,ConvertToID(pos,0));
596 for (
size_t s = 1 ; s < dim ; s++)
599 if ((
size_t)(t.transform(pos,s) / box_unit.getHigh(s)) + off[s] < 0)
601 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" point is not inside the cell space\n";
602 ACTION_ON_ERROR(CELL_DECOMPOSER);
605 key.
set_d(s,ConvertToID(pos,s));
626 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" using an uninitialized CellDecomposer";
627 ACTION_ON_ERROR(CELL_DECOMPOSER);
630 if (gr_cell2.size(0) < k.get(0) + off[0] - cell_shift.get(0))
632 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" cell coordinate 0 = " << k.get(0) + off[0] - cell_shift.get(0) <<
" bigger than cell space " << gr_cell2.size(0) << std::endl;
633 ACTION_ON_ERROR(CELL_DECOMPOSER);
637 size_t cell_id = k.get(0) + off[0] - cell_shift.get(0);
639 for (
size_t s = 1 ; s < dim ; s++)
642 if (gr_cell2.size(s) < k.get(s) + off[s] - cell_shift.get(s))
644 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" cell coordinate 0 = " << k.get(0) + off[0] - cell_shift.get(0) <<
" bigger than cell space " << gr_cell2.size(s) << std::endl;
645 ACTION_ON_ERROR(CELL_DECOMPOSER);
648 cell_id += gr_cell2.size_s(s-1) * (k.get(s) + off[s] -cell_shift.get(s));
668 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" using an uninitialized CellDecomposer";
669 ACTION_ON_ERROR(CELL_DECOMPOSER);
672 if (gr_cell2.size(0) < k.
get(0) + off[0] - cell_shift.get(0))
674 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" cell coordinate 0 = " << k.
get(0) + off[0] - cell_shift.get(0) <<
" bigger than cell space " << gr_cell2.size(0) << std::endl;
675 ACTION_ON_ERROR(CELL_DECOMPOSER);
679 size_t cell_id = k.
get(0) + off[0] -cell_shift.get(0);
681 for (
size_t s = 1 ; s < dim ; s++)
684 if (gr_cell2.size(s) < k.
get(s) + off[s] - cell_shift.get(s))
686 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" cell coordinate 0 = " << k.
get(0) + off[0] - cell_shift.get(0) <<
" bigger than cell space " << gr_cell2.size(s) << std::endl;
687 ACTION_ON_ERROR(CELL_DECOMPOSER);
690 cell_id += gr_cell2.size_s(s-1) * (k.
get(s) + off[s] -cell_shift.get(s));
710 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" using an uninitialized CellDecomposer" << std::endl;
711 ACTION_ON_ERROR(CELL_DECOMPOSER);
716 key.
set_d(0,ConvertToID(pos,0));
718 for (
size_t s = 1 ; s < dim ; s++)
721 if ((
size_t)(t.transform(pos,s) / box_unit.getHigh(s)) + off[s] < 0)
723 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" point is not inside the cell space" << std::endl;
724 ACTION_ON_ERROR(CELL_DECOMPOSER);
728 key.
set_d(s,ConvertToID(pos,s));
745 inline size_t getCellDom(
const Point<dim,T> & pos)
const
747 return getCellDom_impl<Point<dim,T>>(pos);
762 inline size_t getCellDom(
const T (& pos)[dim])
const
764 return getCellDom_impl<T[dim]>(pos);
778 inline size_t getCellPad(
const Point<dim,T> & pos)
const
780 return getCellPad_impl<Point<dim,T>>(pos);
794 inline size_t getCellPad(
const T (& pos)[dim])
const
796 return getCellPad_impl<T[dim]>(pos);
808 inline size_t getCell(
const T (& pos)[dim])
const
813 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" using an uninitialized CellDecomposer";
814 ACTION_ON_ERROR(CELL_DECOMPOSER);
817 if (pos[0] < box.getLow(0) - off[0]*box_unit.getP2()[0] || pos[0] > box.getHigh(0) + off[0]*box_unit.getP2()[0])
819 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" point " << toPointString(pos) <<
" is not inside the cell space";
820 ACTION_ON_ERROR(CELL_DECOMPOSER);
824 size_t cell_id = ConvertToID(pos,0);
826 for (
size_t s = 1 ; s < dim ; s++)
829 if (pos[s] < box.getLow(s) - off[s]*box_unit.getP2()[s] || pos[s] > box.getHigh(s) + off[s]*box_unit.getP2()[s])
831 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" point " << toPointString(pos) <<
" is not inside the cell space";
832 ACTION_ON_ERROR(CELL_DECOMPOSER);
835 cell_id += gr_cell2.size_s(s-1) * ConvertToID(pos,s);
855 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" using an uninitialized CellDecomposer";
856 ACTION_ON_ERROR(CELL_DECOMPOSER);
859 if (pos.
get(0) < box.getLow(0) - off[0]*box_unit.getP2()[0] || pos.
get(0) > box.getHigh(0) + off[0]*box_unit.getP2()[0])
861 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" point " << pos.
toPointString() <<
" is not inside the cell space";
862 ACTION_ON_ERROR(CELL_DECOMPOSER);
866 size_t cell_id = ConvertToID(pos,0);
868 for (
size_t s = 1 ; s < dim ; s++)
871 if (pos.
get(s) < box.getLow(s) - off[s]*box_unit.getP2()[s] || pos.
get(s) > box.getHigh(s) + off[s]*box_unit.getP2()[s])
873 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" point " << pos.
toPointString() <<
" is not inside the cell space";
874 ACTION_ON_ERROR(CELL_DECOMPOSER);
878 cell_id += gr_cell2.size_s(s-1) * ConvertToID(pos,s);
893 template<
typename Mem>
inline size_t getCell(
const encapc<1,
Point<dim,T>,Mem> & pos)
const
899 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" using an uninitialized CellDecomposer";
900 ACTION_ON_ERROR(CELL_DECOMPOSER);
903 if (pos.template get<0>()[0] < box.getLow(0) - off[0]*box_unit.getP2()[0] || pos.template get<0>()[0] > box.getHigh(0) + off[0]*box_unit.getP2()[0])
905 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" point " << toPointString(pos) <<
" is not inside the cell space " << box.toString() << std::endl;
906 ACTION_ON_ERROR(CELL_DECOMPOSER);
910 size_t cell_id = ConvertToID_(pos,0);
912 for (
size_t s = 1 ; s < dim ; s++)
916 if (pos.template get<0>()[s] < box.getLow(s) - off[s]*box_unit.getP2()[s] || pos.template get<0>()[s] > box.getHigh(s) + off[s]*box_unit.getP2()[s])
918 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" point " << toPointString(pos) <<
" is not inside the cell space";
919 ACTION_ON_ERROR(CELL_DECOMPOSER);
922 cell_id += gr_cell2.size_s(s-1) * ConvertToID_(pos,s);
973 bx.
setP2(getCellGrid(p2));
974 bx.
setP1(getCellGrid(p1));
986 inline void setDimensions(
const Box<dim,T> & box,
const size_t (&div)[dim],
const size_t (&div2)[dim],
const size_t pad,
Point<dim,long int> cell_shift)
990 t.setTransform(mat,box.
getP1());
1006 for (
size_t i = 0 ; i < dim ; i++)
1007 cells[i] = div2[i] + 2*pad;
1009 gr_cell2.setDimensions(cells);
1011 for (
size_t i = 0 ; i < dim ; i++)
1012 this->cell_shift.
get(i) = cell_shift.
get(i) - off[i];
1022 inline void setDimensions(
const Box<dim,T> & box,
const size_t (&div)[dim],
const size_t pad)
1026 t.setTransform(mat,box.
getP1());
1028 Initialize(pad,div);
1029 this->cell_shift = 0;
1044 t.setTransform(mat,box.
getP1());
1047 Initialize(pad,div);
1051 size_t div_with_off[dim];
1053 for(
size_t i = 0 ; i < dim ; i++)
1054 div_with_off[i] = div2[i] + 2*off[i];
1056 gr_cell2.setDimensions(div_with_off);
1058 for (
size_t i = 0 ; i < dim ; i++)
1059 this->cell_shift.
get(i) = cell_shift.
get(i) - off[i];
1071 inline void setDimensions(
const Box<dim,T> & box,
const size_t (&div)[dim],
const Matrix<dim,T> & mat,
const size_t pad)
1073 t.setTransform(mat,box.
getP1());
1075 Initialize(pad,div);
1076 this->cell_shift = 0;
1091 inline void setDimensions(
const CellDecomposer_sm<dim,T,transform> & cd,
const Box<dim,size_t> & cell_extension)
1093 this->cell_shift = 0;
1097 t.setTransform(cd.getMat(),cd.getOrig());
1106 for (
size_t i = 0 ; i < dim ; i++)
1108 if (pad < cell_extension.
getLow(i))
1109 pad = cell_extension.
getLow(i);
1110 else if (pad > cell_extension.
getHigh(i))
1111 pad = cell_extension.
getHigh(i);
1118 for (
size_t i = 0 ; i < dim ; i++)
1119 sz_div[i] = cd.gr_cell.size(i) - 2*cd.off[i];
1121 Initialize(pad,sz_div);
1154 :t(
Matrix<dim,T>::identity(),box.getP1()),box(box),gr_cell()
1188 CellDecomposer_sm(
const SpaceBox<dim,T> & box,
const size_t (&div)[dim],
const size_t pad)
1189 :t(
Matrix<dim,T>::identity(),box.getP1()),box(box),gr_cell(div)
1191 Initialize(pad,div);
1201 CellDecomposer_sm(
const CellDecomposer_sm<dim,T,transform> & cd,
Box<dim,size_t> & ext)
1202 :t(
Matrix<dim,T>::identity(),cd.getOrig())
1204 setDimensions(cd,ext);
1210 :t(
Matrix<dim,T>::identity(),
Point<dim,T>::zero_p()),tot_n_cell(0)
1296 b /= getCellBox().
getP2();
1310 for (
size_t i = 0 ; i < dim ; i++)
1321 for (
size_t i = 0 ; i < dim ; i++)
1326 if (bc[i] == NON_PERIODIC)
1329 g_box.
setHigh(i,gr_cell.size(i) - off[i]);
1335 g_box.
setHigh(i,gr_cell.size(i)-1 - off[i]);
1342 if (bc[i] == NON_PERIODIC)
1347 g_box.
setLow(i,gr_cell.size(i) - off[i]);
1353 g_box.
setLow(i,gr_cell.size(i) - off[i]);
1422 b /= getCellBox().
getP2();
1433 for (
size_t i = 0 ; i < dim ; i++) {b.
setHigh(i,b.
getHigh(i)-1);}
1442 for (
size_t i = 0 ; i < dim ; i++)
1447 if (bc[i] == NON_PERIODIC)
1450 g_box.
setHigh(i,gr_cell.size(i) - off[i]);
1456 g_box.
setHigh(i,gr_cell.size(i)-1 - off[i]);
1463 if (bc[i] == NON_PERIODIC)
1468 g_box.
setLow(i,gr_cell.size(i) - off[i]);
1474 g_box.
setLow(i,gr_cell.size(i) - off[i]);
1500 b /= getCellBox().
getP2();
1507 for (
size_t i = 0 ; i < dim ; i++)
1520 for (
size_t i = 0 ; i < dim ; i++)
1575 for (
size_t i = 0 ; i < dim ; i++)
1577 if ((
long int)gr_cell.size(i) - (
long int)off[i] == b_d.
getLow(i))
1579 else if ((
long int)off[i] == b_d.
getLow(i))
1584 if ((
long int)gr_cell.size(i) - (
long int)off[i] == b_d.
getHigh(i))
1586 else if ((
long int)off[i] == b_d.
getHigh(i))
1600 const size_t (& getDiv()
const)[dim]
1602 return gr_cell.getSize();
1611 const size_t (& getDivWP()
const)[dim]
1613 for (
size_t i = 0 ; i < dim ; i++)
1614 {div_wp[i] = gr_cell.getSize()[i] - 2*getPadding(i);}
1634 inline void swap(CellDecomposer_sm<dim,T,transform> & cd)
1637 p_middle.
swap(cd.p_middle);
1645 size_t tot_n_cell_t = tot_n_cell;
1646 tot_n_cell = cd.tot_n_cell;
1647 cd.tot_n_cell = tot_n_cell_t;
1650 box_unit.swap(cd.box_unit);
1651 gr_cell.swap(cd.gr_cell);
1652 gr_cell2.swap(cd.gr_cell2);
1654 for (
size_t i = 0 ; i < dim ; i++)
1656 size_t off_t = off[i];
1660 size_t cs_t = cell_shift.
get(i);
1661 cell_shift.
get(i) = cd.cell_shift.get(i);
1662 cd.cell_shift.get(i) = cs_t;
1673 inline bool operator==(
const CellDecomposer_sm<dim,T,transform> & cd)
1681 if (tot_n_cell != cd.tot_n_cell)
1687 if (box_unit != cd.box_unit)
1690 if (gr_cell != cd.gr_cell)
1693 if (gr_cell2 != cd.gr_cell2)
1696 for (
size_t i = 0 ; i < dim ; i++)
1698 if (off[i] != cd.off[i])
1701 if (cell_shift.
get(i) != cd.cell_shift.get(i))
1715 inline bool operator!=(
const CellDecomposer_sm<dim,T,transform> & cd)
1717 return ! this->operator==(cd);
1727 size_t getPadding(
size_t i)
const
1738 size_t (& getPadding())[dim]
1755 for (
size_t i = 0 ; i < dim ; i++)
1775 for (
size_t i = 0 ; i < dim ; i++)
1777 key.
set_d(i,cell_shift.
get(i) + gr_cell2.size(i) - 2*getPadding(i) - 1);
1792 for (
size_t i = 0 ; i < dim ; i++)
bool operator!=(const shift< dim, T > &s)
It return true if the shift is different.
This class represent an N-dimensional box.
void setP1(const grid_key_dx< dim > &p1)
Set the point P1 of the box.
T getLow(int i) const
get the i-coordinate of the low bound interval of the box
const Point< dim, T > & getOrig() const
Get the shift vector.
grid_key_dx is the key to access any element in the grid
T getHigh(int i) const
get the high interval of the box
void setHigh(int i, T val)
set the high interval of the box
void floorP2()
Apply the ceil operation to the point P2.
shift(const Matrix< dim, T > &t, const Point< dim, T > &s)
Constructor.
This class implement the point shape in an N-dimensional space.
Point< dim, T > getP1() const
Get the point p1.
void swap(Box< dim, T > &b)
exchange the data of two boxes
void setTransform(const Matrix< dim, T > &mat, const Point< dim, T > &orig)
Set the transformation Matrix and shift.
T transform(const T(&s)[dim], const size_t i) const
Shift the point transformation.
mem_id get(size_t i) const
Get the i index.
void setP2(const grid_key_dx< dim > &p2)
Set the point P2 of the box.
void floorP1()
Apply the ceil operation to the point P1.
Matrix< dim, T > mat
Matrix transformation.
void ceilP1()
Apply the ceil operation to the point P1.
bool operator==(const shift< dim, T > &s)
It return true if the shift match.
This class implement an NxN (dense) matrix.
const T & get(size_t i) const
Get coordinate.
void setLow(int i, T val)
set the low interval of the box
This class represent an N-dimensional box.
This class is a trick to indicate the compiler a specific specialization pattern. ...
T transform(const encapc< 1, Point< dim, T >, Mem > &s, const size_t i) const
Shift the point transformation.
static Matrix< dim, T > identity()
Identity matrix.
Point< dim, T > sh
Shift point.
T transform(const Point< dim, T > &s, const size_t i) const
Shift the point transformation.
void ceilP2()
Apply the ceil operation to the point P2.
void set_d(size_t i, mem_id id)
Set the i index.
std::string toPointString() const
Convert the point into a string.
const Matrix< dim, T > & getMat() const
Get the transformation Matrix.
Point< dim, T > getP2() const
Get the point p2.