1/*
2 * CellDecomposer.hpp
3 *
4 * Created on: Apr 1, 2015
5 * Author: Pietro Incardona
6 */
7
8#ifndef CELLDECOMPOSER_HPP_
9#define CELLDECOMPOSER_HPP_
10
11#include "Space/SpaceBox.hpp"
12#include "Space/Matrix.hpp"
13#include "util/copy_compare/meta_compare.hpp"
14#include "Grid/grid_sm.hpp"
15
16#define CELL_DECOMPOSER 8001lu
17
18
19
20
21/*! This Class apply a shift transformation before converting to Cell-ID
22 *
23 *
24 */
25template<unsigned int dim, typename T>
26class shift
27{
28 //! Shift point
29 Point<dim,T> sh;
30
31 //! Matrix transformation
32 Matrix<dim,T> mat;
33
34public:
35
36 /*! \brief Constructor
37 *
38 * \param t Matrix transformation
39 * \param s shift
40 *
41 */
42 __device__ __host__ shift(const Matrix<dim,T> & t, const Point<dim,T> & s)
43 :sh(s)
44 {
45 mat.identity();
46 }
47
48 /*! \brief Shift the point transformation
49 *
50 * \param s source point
51 * \param i coordinate
52 *
53 * \return the transformed coordinate
54 *
55 */
56 __device__ __host__ inline T transform(const T(&s)[dim], const size_t i) const
57 {
58 return s[i] - sh.get(i);
59 }
60
61 /*! \brief Shift the point transformation
62 *
63 * \param s source point
64 * \param i coordinate
65 *
66 * \return the transformed coordinate
67 *
68 */
69 __device__ __host__ inline T transform(const Point<dim,T> & s, const size_t i) const
70 {
71 return s.get(i) - sh.get(i);
72 }
73
74 /*! \brief Shift the point transformation
75 *
76 * \param s source point
77 * \param i coordinate
78 *
79 * \return the transformed coordinate
80 *
81 */
82 template<typename Mem> inline __device__ __host__ T transform(const encapc<1,Point<dim,T>,Mem> & s, const size_t i) const
83 {
84 return s.template get<0>()[i] - sh.get(i);
85 }
86
87 /*! \brief Set the transformation Matrix and shift
88 *
89 * \param mat Matrix transformation
90 * \param orig origin point
91 *
92 */
93 __device__ __host__ inline void setTransform(const Matrix<dim,T> & mat, const Point<dim,T> & orig)
94 {
95 for (size_t i = 0 ; i < dim ; i++)
96 sh.get(i) = orig.get(i);
97 }
98
99 /*! \brief Get the shift vector
100 *
101 * \return the shift vector
102 *
103 */
104 inline const Point<dim,T> & getOrig() const
105 {
106 return sh;
107 }
108
109 /*! \brief Get the transformation Matrix
110 *
111 * \return the transformation matrix
112 *
113 */
114 inline const Matrix<dim,T> & getMat() const
115 {
116 return mat;
117 }
118
119 /*! \brief It return true if the shift match
120 *
121 * \param s shift to compare with
122 *
123 * \return true if it match
124 *
125 */
126 inline bool operator==(const shift<dim,T> & s)
127 {
128 return sh == s.sh;
129 }
130
131 /*! \brief It return true if the shift is different
132 *
133 * \param s shift to compare with
134 *
135 * \return true if the shift is different
136 *
137 */
138 inline bool operator!=(const shift<dim,T> & s)
139 {
140 return !this->operator==(s);
141 }
142};
143
144/*! This Class apply a shift transformation before converting to Cell-ID
145 *
146 *
147 */
148template<unsigned int dim, typename T>
149class shift_only
150{
151 //! Shift point
152 Point<dim,T> sh;
153
154public:
155
156 /*! \brief Default constructor
157 *
158 */
159 shift_only()
160 {
161 sh.zero();
162 }
163
164 /*! \brief Constructor
165 *
166 * \param t Matrix transformation
167 * \param s shift
168 *
169 */
170 shift_only(const Matrix<dim,T> & t, const Point<dim,T> & s)
171 :sh(s)
172 {
173 }
174
175 /*! \brief Shift the point transformation
176 *
177 * \param s source point
178 * \param i coordinate
179 *
180 * \return the transformed coordinate
181 *
182 */
183 __device__ __host__ inline T transform(const T * s, const int i) const
184 {
185 return s[i] - sh.get(i);
186 }
187
188 /*! \brief Shift the point transformation
189 *
190 * \param s source point
191 * \param i coordinate
192 *
193 * \return the transformed coordinate
194 *
195 */
196 __device__ __host__ inline T transform(const T(&s)[dim], const int i) const
197 {
198 return s[i] - sh.get(i);
199 }
200
201 /*! \brief Shift the point transformation
202 *
203 * \param s source point
204 * \param i coordinate
205 *
206 * \return the transformed coordinate
207 *
208 */
209 __device__ __host__ inline T transform(const Point<dim,T> & s, const int i) const
210 {
211 return s.get(i) - sh.get(i);
212 }
213
214 /*! \brief Shift the point transformation
215 *
216 * \param s source point
217 * \param i coordinate
218 *
219 * \return the transformed coordinate
220 *
221 */
222 template<typename Mem> __device__ __host__ inline T transform(const encapc<1,Point<dim,T>,Mem> & s, const int i) const
223 {
224 return s.template get<0>()[i] - sh.get(i);
225 }
226
227 /*! \brief Set the transformation Matrix and shift
228 *
229 * \param mat Matrix transformation
230 * \param orig origin point
231 *
232 */
233 inline void setTransform(const Matrix<dim,T> & mat, const Point<dim,T> & orig)
234 {
235 for (size_t i = 0 ; i < dim ; i++)
236 sh.get(i) = orig.get(i);
237 }
238
239 /*! \brief Get the shift vector
240 *
241 * \return the shift vector
242 *
243 */
244 __device__ __host__ inline const Point<dim,T> & getOrig() const
245 {
246 return sh;
247 }
248
249
250 /*! \brief It return true if the shift match
251 *
252 * \param s shift to compare with
253 *
254 * \return true if it match
255 *
256 */
257 inline bool operator==(const shift<dim,T> & s)
258 {
259 return sh == s.sh;
260 }
261
262 /*! \brief It return true if the shift is different
263 *
264 * \param s shift to compare with
265 *
266 * \return true if the shift is different
267 *
268 */
269 inline bool operator!=(const shift<dim,T> & s)
270 {
271 return !this->operator==(s);
272 }
273};
274
275//! No transformation
276template<unsigned int dim, typename T>
277class no_transform
278{
279 //! shift transform
280 Point<dim,T> orig;
281
282 //! Matrix transform
283 Matrix<dim,T> mat;
284
285public:
286
287 /*! \brief Constructor
288 *
289 * \param t Matrix transformation
290 * \param s shift
291 *
292 */
293 __device__ __host__ no_transform(const Matrix<dim,T> & t, const Point<dim,T> & s)
294 {
295 orig.zero();
296 mat.identity();
297 }
298
299 /*! \brief Shift the point transformation
300 *
301 * \param s source point
302 * \param i coordinate
303 *
304 * \return the transformed coordinate
305 *
306 */
307 __device__ __host__ inline T transform(const T(&s)[dim], const size_t i) const
308 {
309 return s[i];
310 }
311
312 /*! \brief No transformation
313 *
314 * \param s source point
315 * \param i coordinate
316 *
317 * \return the source point coordinate
318 *
319 */
320 __device__ __host__ inline T transform(const Point<dim,T> & s, const size_t i) const
321 {
322 return s.get(i);
323 }
324
325 /*! \brief No transformation
326 *
327 * \param s source point
328 * \param i coordinate
329 *
330 * \return the point coordinate
331 *
332 */
333 template<typename Mem> __device__ __host__ inline T transform(const encapc<1,Point<dim,T>,Mem> & s, const size_t i) const
334 {
335 return s.template get<Point<dim,T>::x>()[i];
336 }
337
338 /*! \brief Set the transformation Matrix and shift
339 *
340 * \param mat Matrix transformation
341 * \param orig origin point
342 *
343 */
344 inline void setTransform(const Matrix<dim,T> & mat, const Point<dim,T> & orig)
345 {
346
347 }
348
349 /*! \brief It return always true true
350 *
351 * There is nothing to compare
352 *
353 * \param nt unused
354 *
355 * \return true
356 *
357 */
358 inline bool operator==(const no_transform<dim,T> & nt)
359 {
360 return true;
361 }
362
363 /*! \brief It return always false
364 *
365 * There is nothing to compare they cannot be differents
366 *
367 * \param nt unused
368 *
369 * \return false
370 *
371 */
372 inline bool operator!=(const no_transform<dim,T> & nt)
373 {
374 return false;
375 }
376
377 /*! \brief Get the origin
378 *
379 * \return the shift vector
380 *
381 */
382 inline const Point<dim,T> & getOrig() const
383 {
384 return orig;
385 }
386
387 /*! \brief Get the transformation Matrix
388 *
389 * \return get the return the matrix
390 *
391 */
392 inline const Matrix<dim,T> & getMat() const
393 {
394 return mat;
395 }
396};
397
398//! No transformation
399template<unsigned int dim, typename T>
400class no_transform_only
401{
402
403public:
404
405 /*! \brief Constructor
406 *
407 * \param t Matrix transformation
408 * \param s shift
409 *
410 */
411 no_transform_only(const Matrix<dim,T> & t, const Point<dim,T> & s)
412 {
413 }
414
415 /*! \brief Shift the point transformation
416 *
417 * \param s source point
418 * \param i coordinate
419 *
420 * \return the transformed coordinate
421 *
422 */
423 __device__ __host__ inline T transform(const T * s, const int i) const
424 {
425 return s[i];
426 }
427
428 /*! \brief Shift the point transformation
429 *
430 * \param s source point
431 * \param i coordinate
432 *
433 * \return the transformed coordinate
434 *
435 */
436 __device__ __host__ inline T transform(const T(&s)[dim], const int i) const
437 {
438 return s[i];
439 }
440
441 /*! \brief No transformation
442 *
443 * \param s source point
444 * \param i coordinate
445 *
446 * \return the source point coordinate
447 *
448 */
449 __device__ __host__ inline T transform(const Point<dim,T> & s, const int i) const
450 {
451 return s.get(i);
452 }
453
454 /*! \brief No transformation
455 *
456 * \param s source point
457 * \param i coordinate
458 *
459 * \return the point coordinate
460 *
461 */
462 template<typename Mem> __device__ __host__ inline T transform(const encapc<1,Point<dim,T>,Mem> & s, const int i) const
463 {
464 return s.template get<Point<dim,T>::x>()[i];
465 }
466
467 /*! \brief Set the transformation Matrix and shift
468 *
469 * \param mat Matrix transformation
470 * \param orig origin point
471 *
472 */
473 inline void setTransform(const Matrix<dim,T> & mat, const Point<dim,T> & orig)
474 {
475
476 }
477
478 /*! \brief It return always true true
479 *
480 * There is nothing to compare
481 *
482 * \param nt unused
483 *
484 * \return true
485 *
486 */
487 inline bool operator==(const no_transform<dim,T> & nt)
488 {
489 return true;
490 }
491
492 /*! \brief It return always false
493 *
494 * There is nothing to compare they cannot be differents
495 *
496 * \param nt unused
497 *
498 * \return false
499 *
500 */
501 inline bool operator!=(const no_transform<dim,T> & nt)
502 {
503 return false;
504 }
505};
506
507/*! \brief Decompose a space into cells
508 *
509 * It is a convenient class for cell decomposition of an N dimensional space into cells
510 * and index linearization with getCell. Transform parameter is used to shift the space
511 * from the origin where the cell list is defined. It basically apply a transformation to the
512 * point every time we call getCell of getCellGrid. The shift is given by the starting point (p1)
513 * of the box that define where the Cell decomposition live
514 *
515 * \tparam dim dimension of the space divided by the cell
516 * \tparam T information the cell contain
517 * \tparam transform type of transformation (no_transformation shift ... ), carefully if p1 it different from {0, ... 0}
518 * shift class must be used
519 *
520 * Example of a 2D space decompised into Cells, 6x6 structure with padding 1 without shift, cell indicated with p are padding cell
521 * the origin of the cell or point (0,0) is marked with cell number 9
522 *
523 * \verbatim
524 * +-----------------------+
525 * |p |p |p |p |p |p |p |p |
526 * +-----------------------+
527 * |p | | | | | | |p |
528 * +-----------------------+
529 * |p | | | | | | |p |
530 * +-----------------------+
531 * |p | | | | | | |p |
532 * +-----------------------+
533 * |p |9 | | | | | |p |
534 * +-----------------------+
535 * |p |p |p |p |p |p |p |p |
536 * +-----------------------+
537 * \endverbatim
538 *
539 * ### Cell decomposer without shift
540 * \snippet CellDecomposer_unit_tests.hpp Cell decomposer use without shift
541 * ### Cell decomposer with padding
542 * \snippet CellDecomposer_unit_tests.hpp Test Cell decomposer with padding
543 * ### Cell decomposer with shift
544 * \snippet CellDecompose_unit_tests.hpp Test Cell decomposer with shift
545 *
546 */
547template<unsigned int dim,typename T, typename transform = no_transform<dim,T>>
548class CellDecomposer_sm
549{
550 // Point in the middle of a cell
551 //
552 // \verbatim
553 //
554 // C (0.1,0.1)
555 // +-----+
556 // | |
557 // | .P |
558 // | |
559 // +-----+
560 //
561 // \endverbatim
562 //
563 // C is the cell and P is the point inside the middle of the cell
564 // for example if the cell is (0.1,0.1) P is (0.05,0.05)
565 //
566 //
567 Point<dim,T> p_middle;
568
569 // Point transformation before get the Cell object (useful for example to shift the cell list)
570 transform t;
571
572 /*! \brief Convert the coordinates into id
573 *
574 * \param x coordinate
575 * \param s dimension
576 *
577 */
578 inline size_t ConvertToID(const T (&x)[dim] ,size_t s) const
579 {
580 size_t id = openfpm::math::size_t_floor(t.transform(x,s) / box_unit.getHigh(s)) + off[s];
581 id = (id >= gr_cell.size(s))?(gr_cell.size(s)-1-cell_shift.get(s)):id-cell_shift.get(s);
582 return id;
583 }
584
585 /*! \brief Convert the coordinates into id
586 *
587 * \param x point
588 * \param s dimension
589 *
590 */
591 inline size_t ConvertToID(const Point<dim,T> & x ,size_t s, size_t sc = 0) const
592 {
593 size_t id = openfpm::math::size_t_floor(t.transform(x,s) / box_unit.getHigh(s)) + off[s];
594 id = (id >= gr_cell.size(s))?(gr_cell.size(s)-1-cell_shift.get(s)):id-cell_shift.get(s);
595 return id;
596 }
597
598 /*! \brief Convert the coordinates into id with negative machine precision expansion
599 *
600 * \param x point
601 * \param s dimension
602 *
603 */
604 inline size_t ConvertToID_me(const Point<dim,T> & x ,size_t s, size_t sc = 0) const
605 {
606 T cc = t.transform(x,s) / box_unit.getHigh(s) - 0.015625;
607 size_t id = openfpm::math::size_t_floor(cc) + off[s];
608 id = (id >= gr_cell.size(s))?(gr_cell.size(s)-1-cell_shift.get(s)):id-cell_shift.get(s);
609 return id;
610 }
611
612 /*! \brief Convert the coordinates into id with positive machine precision expansion
613 *
614 * \param x point
615 * \param s dimension
616 *
617 */
618 inline size_t ConvertToID_pe(const Point<dim,T> & x ,size_t s, size_t sc = 0) const
619 {
620 T cc = t.transform(x,s) / box_unit.getHigh(s) + 0.015625;
621 size_t id = openfpm::math::size_t_floor(cc) + off[s];
622 id = (id >= gr_cell.size(s))?(gr_cell.size(s)-1-cell_shift.get(s)):id-cell_shift.get(s);
623 return id;
624 }
625
626 /*! \brief Convert the coordinates into id without apply shift
627 *
628 * \param x coordinate
629 * \param s dimension
630 *
631 */
632 inline size_t ConvertToID_ns(const T (&x)[dim] ,size_t s) const
633 {
634 size_t id = openfpm::math::size_t_floor(t.transform(x,s) / box_unit.getHigh(s)) + off[s];
635 id = (id >= gr_cell.size(s))?(gr_cell.size(s)-1):id;
636 return id;
637 }
638
639 /*! \brief Convert the coordinates into id without apply shift
640 *
641 * \param x point
642 * \param s dimension
643 *
644 */
645 inline size_t ConvertToID_ns(const Point<dim,T> & x ,size_t s, size_t sc = 0) const
646 {
647 size_t id = openfpm::math::size_t_floor(t.transform(x,s) / box_unit.getHigh(s)) + off[s];
648 id = (id >= gr_cell.size(s))?(gr_cell.size(s)-1):id;
649 return id;
650 }
651
652 /*! \brief Convert the coordinates into id
653 *
654 * \param x point
655 * \param s dimension
656 *
657 */
658 template <typename Mem> inline size_t ConvertToID_(const encapc<1,Point<dim,T>,Mem> & x ,size_t s, size_t sc = 0) const
659 {
660 size_t id = (size_t)(t.transform(x,s) / box_unit.getHigh(s)) + off[s];
661 id = (id >= gr_cell.size(s))?(gr_cell.size(s)-1-cell_shift.get(s)):id-cell_shift.get(s);
662 return id;
663 }
664
665 /*! \Check if a particle is outside the domain of the cell-list
666 *
667 * \param pos position of the particle
668 * \param s coordinate to check
669 *
670 */
671 template<typename Ele> inline void check_and_print_error(const Ele & pos ,size_t s) const
672 {
673#ifdef SE_CLASS1
674 if (tot_n_cell == 0)
675 {
676 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " using an uninitialized CellDecomposer";
677 ACTION_ON_ERROR(CELL_DECOMPOSER);
678 }
679
680 if (pos[s] < box.getLow(s) - off[s]*box_unit.getP2()[s] || pos[s] > box.getHigh(s) + off[s]*box_unit.getP2()[s])
681 {
682 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point " << Point<dim,T>(pos).toString() << " is not inside the cell space";
683 ACTION_ON_ERROR(CELL_DECOMPOSER);
684 }
685#endif
686 }
687
688
689 template<typename Ele> inline size_t getCellDom_impl(const Ele & pos) const
690 {
691 check_and_print_error(pos,0);
692
693 size_t cell_id = ConvertToID_ns(pos,0);
694
695 cell_id = (cell_id == gr_cell.size(0) - off[0])?gr_cell.size(0) - off[0] - 1:cell_id;
696 cell_id = (cell_id == off[0]-1)?off[0]:cell_id;
697
698 cell_id -= cell_shift.get(0);
699
700 for (size_t s = 1 ; s < dim ; s++)
701 {
702 /* coverity[dead_error_begin] */
703 check_and_print_error(pos,s);
704
705 size_t cell_idt = ConvertToID_ns(pos,s);
706
707 cell_idt = (cell_idt == gr_cell.size(s) - off[s])?gr_cell.size(s) - off[s] - 1:cell_idt;
708 cell_idt = (cell_idt == off[s]-1)?off[s]:cell_idt;
709
710 cell_idt -= cell_shift.get(s);
711
712 cell_id += gr_cell2.size_s(s-1) * cell_idt;
713 }
714
715 return cell_id;
716 }
717
718 template<typename Ele> inline size_t getCellPad_impl(const Ele & pos) const
719 {
720 check_and_print_error(pos,0);
721
722 size_t cell_id = ConvertToID_ns(pos,0);
723 cell_id = (cell_id == off[0])?off[0]-1:cell_id;
724 cell_id = (cell_id == gr_cell.size(0) - off[0] - 1)?gr_cell.size(0) - off[0]:cell_id;
725
726 cell_id -= cell_shift.get(0);
727
728 for (size_t s = 1 ; s < dim ; s++)
729 {
730 check_and_print_error(pos,s);
731
732 size_t cell_idt = ConvertToID_ns(pos,s);
733 cell_idt = (cell_idt == off[s])?off[s]-1:cell_idt;
734 cell_idt = (cell_idt == gr_cell.size(s) - off[s] - 1)?gr_cell.size(s) - off[s]:cell_idt;
735
736 cell_idt -= cell_shift.get(s);
737
738 cell_id += gr_cell2.size_s(s-1) * cell_idt;
739 }
740
741 return cell_id;
742 }
743
744
745protected:
746
747 //! Total number of cell
748 size_t tot_n_cell;
749
750 //! Domain of the cell list
751 SpaceBox<dim,T> box;
752
753 //! Unit box of the Cell list
754 SpaceBox<dim,T> box_unit;
755
756 //! Grid structure of the Cell list
757 grid_sm<dim,void> gr_cell;
758
759 //! Grid structure of the internal Cell list
760 grid_sm<dim,void> gr_cell2;
761
762 //! Box in continuum for the gr_cell2
763 Box<dim,T> box_gr_cell2;
764
765 //! cell padding on each dimension
766 size_t off[dim];
767
768 //! cell_shift
769 Point<dim,long int> cell_shift;
770
771 // temporal buffer
772 mutable size_t div_wp[dim];
773
774 /*! \brief Initialize all the structures
775 *
776 */
777 void Initialize(const size_t pad, const size_t (& div)[dim])
778 {
779#ifdef SE_CLASS1
780
781 for (size_t i = 0 ; i < dim ; i++)
782 {
783 if (div[i] == 0)
784 {
785 std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " the number of cells on each dimension must be different from zero\n";
786 ACTION_ON_ERROR(CELL_DECOMPOSER)
787 }
788 }
789
790#endif
791
792 // set div_wp to zero
793 for (size_t i = 0 ; i < dim ; i++)
794 {div_wp[i] = 0;}
795
796 // created a padded div
797 size_t div_p[dim];
798
799 for (size_t i = 0 ; i < dim ; i++)
800 div_p[i] = div[i] + 2*pad;
801
802 gr_cell.setDimensions(div_p);
803 gr_cell2.setDimensions(div_p);
804
805 tot_n_cell = 1;
806
807 // Total number of cells and calculate the unit cell size
808
809 for (size_t i = 0 ; i < dim ; i++)
810 {
811 tot_n_cell *= gr_cell.size(i);
812
813 // Cell are padded by
814 box_unit.setHigh(i,(box.getHigh(i) - box.getLow(i)) / (gr_cell.size(i)- 2*pad) );
815 }
816
817 for (size_t i = 0; i < dim ; i++)
818 off[i] = pad;
819
820 // Initialize p_middle
821
822 p_middle = box_unit.getP2();
823 p_middle = p_middle / 2;
824 }
825
826public:
827
828 /*! \brief Return the transformation
829 *
830 * \return the transform
831 *
832 */
833 const transform & getTransform()
834 {
835 return t;
836 }
837
838 /*! \brief Return the underlying grid information of the cell list
839 *
840 * \return the grid infos
841 *
842 */
843 const grid_sm<dim,void> & getGrid() const
844 {
845#ifdef DEBUG
846 if (tot_n_cell == 0)
847 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " using an uninitialized CellDecomposer";
848#endif
849
850 return gr_cell;
851 }
852
853 /*! \brief Get the cell-ids
854 *
855 * Convert the point coordinates into the cell ids
856 *
857 * \param pos Point position
858 *
859 * \return the cell-ids ad a grid_key_dx<dim>
860 *
861 */
862 inline grid_key_dx<dim> getCellGrid(const T (& pos)[dim]) const
863 {
864#ifdef SE_CLASS1
865 if (tot_n_cell == 0)
866 {
867 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " using an uninitialized CellDecomposer";
868 ACTION_ON_ERROR(CELL_DECOMPOSER);
869 }
870#endif
871
872 grid_key_dx<dim> key;
873 key.set_d(0,ConvertToID(pos,0));
874
875 for (size_t s = 1 ; s < dim ; s++)
876 {
877#ifdef SE_CLASS1
878 if ((size_t)(t.transform(pos,s) / box_unit.getHigh(s)) + off[s] < 0)
879 {
880 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point is not inside the cell space\n";
881 ACTION_ON_ERROR(CELL_DECOMPOSER);
882 }
883#endif
884 key.set_d(s,ConvertToID(pos,s));
885
886 }
887
888 return key;
889 }
890
891 /*! \brief Get the cell-id
892 *
893 * Convert the point coordinates into the cell id
894 *
895 * \param cell-id in grid units
896 *
897 * \return the cell-id
898 *
899 */
900 inline size_t getCellLin(grid_key_dx<dim> && k) const
901 {
902#ifdef SE_CLASS1
903 if (tot_n_cell == 0)
904 {
905 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " using an uninitialized CellDecomposer";
906 ACTION_ON_ERROR(CELL_DECOMPOSER);
907 }
908
909 if (gr_cell2.size(0) < k.get(0) + off[0] - cell_shift.get(0))
910 {
911 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;
912 ACTION_ON_ERROR(CELL_DECOMPOSER);
913 }
914#endif
915
916 size_t cell_id = k.get(0) + off[0] - cell_shift.get(0);
917
918 for (size_t s = 1 ; s < dim ; s++)
919 {
920#ifdef SE_CLASS1
921 if (gr_cell2.size(s) < k.get(s) + off[s] - cell_shift.get(s))
922 {
923 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;
924 ACTION_ON_ERROR(CELL_DECOMPOSER);
925 }
926#endif
927 cell_id += gr_cell2.size_s(s-1) * (k.get(s) + off[s] -cell_shift.get(s));
928 }
929
930 return cell_id;
931 }
932
933 /*! \brief Get the cell-id
934 *
935 * Convert the point coordinates into the cell id
936 *
937 * \param cell-id in grid units
938 *
939 * \return the cell-id
940 *
941 */
942 inline size_t getCellLin(const grid_key_dx<dim> & k) const
943 {
944#ifdef SE_CLASS1
945 if (tot_n_cell == 0)
946 {
947 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " using an uninitialized CellDecomposer";
948 ACTION_ON_ERROR(CELL_DECOMPOSER);
949 }
950
951 if (gr_cell2.size(0) < k.get(0) + off[0] - cell_shift.get(0))
952 {
953 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;
954 ACTION_ON_ERROR(CELL_DECOMPOSER);
955 }
956#endif
957
958 size_t cell_id = k.get(0) + off[0] -cell_shift.get(0);
959
960 for (size_t s = 1 ; s < dim ; s++)
961 {
962#ifdef SE_CLASS1
963 if (gr_cell2.size(s) < k.get(s) + off[s] - cell_shift.get(s))
964 {
965 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;
966 ACTION_ON_ERROR(CELL_DECOMPOSER);
967 }
968#endif
969 cell_id += gr_cell2.size_s(s-1) * (k.get(s) + off[s] -cell_shift.get(s));
970 }
971
972 return cell_id;
973 }
974
975 /*! \brief Get the cell-ids with negative machine precision expansion
976 *
977 * Convert the point coordinates into the cell ids (Careful it include padding)
978 *
979 * \param pos Point position
980 *
981 * \return the cell-ids ad a grid_key_dx<dim>
982 *
983 */
984 inline grid_key_dx<dim> getCellGrid_me(const Point<dim,T> & pos) const
985 {
986#ifdef SE_CLASS1
987 if (tot_n_cell == 0)
988 {
989 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " using an uninitialized CellDecomposer" << std::endl;
990 ACTION_ON_ERROR(CELL_DECOMPOSER);
991 }
992#endif
993
994 grid_key_dx<dim> key;
995 key.set_d(0,ConvertToID_me(pos,0));
996
997 for (size_t s = 1 ; s < dim ; s++)
998 {
999#ifdef SE_CLASS1
1000 if ((size_t)(t.transform(pos,s) / box_unit.getHigh(s)) + off[s] < 0)
1001 {
1002 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point is not inside the cell space" << std::endl;
1003 ACTION_ON_ERROR(CELL_DECOMPOSER);
1004 }
1005#endif
1006 /* coverity[dead_error_line] */
1007 key.set_d(s,ConvertToID_me(pos,s));
1008 }
1009
1010 return key;
1011 }
1012
1013 /*! \brief Get the cell-ids with positive machine precision expansion
1014 *
1015 * Convert the point coordinates into the cell ids (Careful it include padding)
1016 *
1017 * \param pos Point position
1018 *
1019 * \return the cell-ids ad a grid_key_dx<dim>
1020 *
1021 */
1022 inline grid_key_dx<dim> getCellGrid_pe(const Point<dim,T> & pos) const
1023 {
1024#ifdef SE_CLASS1
1025 if (tot_n_cell == 0)
1026 {
1027 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " using an uninitialized CellDecomposer" << std::endl;
1028 ACTION_ON_ERROR(CELL_DECOMPOSER);
1029 }
1030#endif
1031
1032 grid_key_dx<dim> key;
1033 key.set_d(0,ConvertToID_pe(pos,0));
1034
1035 for (size_t s = 1 ; s < dim ; s++)
1036 {
1037#ifdef SE_CLASS1
1038 if ((size_t)(t.transform(pos,s) / box_unit.getHigh(s)) + off[s] < 0)
1039 {
1040 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point is not inside the cell space" << std::endl;
1041 ACTION_ON_ERROR(CELL_DECOMPOSER);
1042 }
1043#endif
1044 /* coverity[dead_error_line] */
1045 key.set_d(s,ConvertToID_pe(pos,s));
1046 }
1047
1048 return key;
1049 }
1050
1051 /*! \brief Get the cell-ids
1052 *
1053 * Convert the point coordinates into the cell ids (Careful it include padding)
1054 *
1055 * \param pos Point position
1056 *
1057 * \return the cell-ids ad a grid_key_dx<dim>
1058 *
1059 */
1060 inline grid_key_dx<dim> getCellGrid(const Point<dim,T> & pos) const
1061 {
1062#ifdef SE_CLASS1
1063 if (tot_n_cell == 0)
1064 {
1065 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " using an uninitialized CellDecomposer" << std::endl;
1066 ACTION_ON_ERROR(CELL_DECOMPOSER);
1067 }
1068#endif
1069
1070 grid_key_dx<dim> key;
1071 key.set_d(0,ConvertToID(pos,0));
1072
1073 for (size_t s = 1 ; s < dim ; s++)
1074 {
1075#ifdef SE_CLASS1
1076 if ((size_t)(t.transform(pos,s) / box_unit.getHigh(s)) + off[s] < 0)
1077 {
1078 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point is not inside the cell space" << std::endl;
1079 ACTION_ON_ERROR(CELL_DECOMPOSER);
1080 }
1081#endif
1082 /* coverity[dead_error_line] */
1083 key.set_d(s,ConvertToID(pos,s));
1084 }
1085
1086 return key;
1087 }
1088
1089 /*! \brief Get the cell-id enforcing that is NOT a cell from the padding
1090 *
1091 * Convert the point coordinates into the cell id
1092 *
1093 * \note this function is in general used to bypass round-off error
1094 *
1095 * \param pos Point position
1096 *
1097 * \return the cell-id
1098 *
1099 */
1100 inline size_t getCellDom(const Point<dim,T> & pos) const
1101 {
1102 return getCellDom_impl<Point<dim,T>>(pos);
1103 }
1104
1105
1106 /*! \brief Get the cell-id enforcing that is NOT a cell from the padding
1107 *
1108 * Convert the point coordinates into the cell id
1109 *
1110 * \note this function is in general used to bypass round-off error
1111 *
1112 * \param pos Point position
1113 *
1114 * \return the cell-id
1115 *
1116 */
1117 inline size_t getCellDom(const T (& pos)[dim]) const
1118 {
1119 return getCellDom_impl<T[dim]>(pos);
1120 }
1121
1122 /*! \brief Get the cell-id enforcing that is from a padding cell
1123 *
1124 * Convert the point coordinates into the cell id
1125 *
1126 * \note this function is in general used to bypass round-off error
1127 *
1128 * \param pos Point position
1129 *
1130 * \return the cell-id
1131 *
1132 */
1133 inline size_t getCellPad(const Point<dim,T> & pos) const
1134 {
1135 return getCellPad_impl<Point<dim,T>>(pos);
1136 }
1137
1138 /*! \brief Get the cell-id enforcing that is from a padding cell
1139 *
1140 * Convert the point coordinates into the cell id
1141 *
1142 * \note this function is in general used to bypass round-off error
1143 *
1144 * \param pos Point position
1145 *
1146 * \return the cell-id
1147 *
1148 */
1149 inline size_t getCellPad(const T (& pos)[dim]) const
1150 {
1151 return getCellPad_impl<T[dim]>(pos);
1152 }
1153
1154 /*! \brief Get the cell-id
1155 *
1156 * Convert the point coordinates into the cell id
1157 *
1158 * \param pos Point position
1159 *
1160 * \return the cell-id
1161 *
1162 */
1163 inline size_t getCell(const T (& pos)[dim]) const
1164 {
1165#ifdef SE_CLASS1
1166 if (tot_n_cell == 0)
1167 {
1168 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " using an uninitialized CellDecomposer";
1169 ACTION_ON_ERROR(CELL_DECOMPOSER);
1170 }
1171
1172 if (pos[0] < box.getLow(0) - off[0]*box_unit.getP2()[0] || pos[0] > box.getHigh(0) + off[0]*box_unit.getP2()[0])
1173 {
1174 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point " << toPointString(pos) << " is not inside the cell space";
1175 ACTION_ON_ERROR(CELL_DECOMPOSER);
1176 }
1177#endif
1178
1179 size_t cell_id = ConvertToID(pos,0);
1180
1181 for (size_t s = 1 ; s < dim ; s++)
1182 {
1183#ifdef SE_CLASS1
1184 if (pos[s] < box.getLow(s) - off[s]*box_unit.getP2()[s] || pos[s] > box.getHigh(s) + off[s]*box_unit.getP2()[s])
1185 {
1186 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point " << toPointString(pos) << " is not inside the cell space";
1187 ACTION_ON_ERROR(CELL_DECOMPOSER);
1188 }
1189#endif
1190 cell_id += gr_cell2.size_s(s-1) * ConvertToID(pos,s);
1191 }
1192
1193 return cell_id;
1194 }
1195
1196 /*! \brief Get the cell-id
1197 *
1198 * Convert the point coordinates into the cell id
1199 *
1200 * \param pos Point position
1201 *
1202 * \return the cell-id
1203 *
1204 */
1205 inline size_t getCell(const Point<dim,T> & pos) const
1206 {
1207#ifdef SE_CLASS1
1208 if (tot_n_cell == 0)
1209 {
1210 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " using an uninitialized CellDecomposer";
1211 ACTION_ON_ERROR(CELL_DECOMPOSER);
1212 }
1213
1214 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])
1215 {
1216 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point " << pos.toPointString() << " is not inside the cell space";
1217 ACTION_ON_ERROR(CELL_DECOMPOSER);
1218 }
1219#endif
1220
1221 size_t cell_id = ConvertToID(pos,0);
1222
1223 for (size_t s = 1 ; s < dim ; s++)
1224 {
1225#ifdef SE_CLASS1
1226 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])
1227 {
1228 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point " << pos.toPointString() << " is not inside the cell space";
1229 ACTION_ON_ERROR(CELL_DECOMPOSER);
1230 }
1231#endif
1232 /* coverity[dead_error_line] */
1233 cell_id += gr_cell2.size_s(s-1) * ConvertToID(pos,s);
1234 }
1235
1236 return cell_id;
1237 }
1238
1239 /*! \brief Get the cell-id
1240 *
1241 * Convert the point coordinates into the cell id
1242 *
1243 * \param pos Point position
1244 *
1245 * \return the cell-id
1246 *
1247 */
1248 template<typename Mem> inline size_t getCell(const encapc<1,Point<dim,T>,Mem> & pos) const
1249 {
1250
1251#ifdef SE_CLASS1
1252 if (tot_n_cell == 0)
1253 {
1254 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " using an uninitialized CellDecomposer";
1255 ACTION_ON_ERROR(CELL_DECOMPOSER);
1256 }
1257
1258 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])
1259 {
1260 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point " << toPointString(pos) << " is not inside the cell space " << box.toString() << std::endl;
1261 ACTION_ON_ERROR(CELL_DECOMPOSER);
1262 }
1263#endif
1264
1265 size_t cell_id = ConvertToID_(pos,0);
1266
1267 for (size_t s = 1 ; s < dim ; s++)
1268 {
1269#ifdef SE_CLASS1
1270
1271 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])
1272 {
1273 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point " << toPointString(pos) << " is not inside the cell space";
1274 ACTION_ON_ERROR(CELL_DECOMPOSER);
1275 }
1276#endif
1277 cell_id += gr_cell2.size_s(s-1) * ConvertToID_(pos,s);
1278 }
1279
1280 return cell_id;
1281 }
1282
1283 /*! \brief Return the smallest box containing the grid points
1284 *
1285 * Suppose a grid 5x5 defined on a Box<2,float> box({0.0,0.0},{1.0,1.0})
1286 * and feeding to the function a Box<2,float>({0.4,0.4},{0.8,0.8}), it will return
1287 * a Box<2,size_t> (2,2) and (3,3). A visualization it is shown in the
1288 * picture below. (the grid points are centered on each cell)
1289 *
1290 * \verbatim
1291 *
1292 +----------+
1293 |. . . . . |
1294 | |
1295 |. . . . . |
1296 | +---+ |
1297 |. .|. .|. |
1298 | | | |
1299 |. .|. .|. |
1300 | +---+ |
1301 |. . . . . |
1302 +----------+
1303
1304 \endverbatim
1305 *
1306 *
1307 * \div grid size on each dimension
1308 * \box Small box in the picture
1309 *
1310 * The big box is defined by "this" box
1311 *
1312 * \return the box containing the grid points
1313 *
1314 */
1315 inline Box<dim,size_t> getGridPoints(const Box<dim,T> & s_box) const
1316 {
1317 // Box with inside grid
1318 Box<dim,size_t> bx;
1319
1320 // Point p2
1321 Point<dim,T> p2 = s_box.getP2();
1322 p2 = p2 - p_middle;
1323
1324 // Point p1
1325 Point<dim,T> p1 = s_box.getP1();
1326 p1 = p1 + p_middle;
1327
1328 bx.setP2(getCellGrid(p2));
1329 bx.setP1(getCellGrid(p1));
1330
1331 return bx;
1332 }
1333
1334 /*! \brief Set the domain to decompose
1335 *
1336 * \param box Domain to decompose
1337 * \param div array with the number of cells on each dimensions
1338 * \param pad padding cell
1339 *
1340 */
1341 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)
1342 {
1343 Matrix<dim,T> mat;
1344 mat.identity();
1345 t.setTransform(mat,box.getP1());
1346 this->box = box;
1347
1348 Initialize(pad,div);
1349
1350 size_t cells[dim];
1351
1352 for (size_t i = 0 ; i < dim ; i++)
1353 cells[i] = div2[i] + 2*pad;
1354
1355 gr_cell2.setDimensions(cells);
1356
1357 for (size_t i = 0 ; i < dim ; i++)
1358 this->cell_shift.get(i) = cell_shift.get(i) - off[i];
1359 }
1360
1361 /*! \brief Set the domain to decompose
1362 *
1363 * \param box Domain to decompose
1364 * \param div array with the number of cells on each dimensions
1365 * \param pad padding cell
1366 *
1367 */
1368 inline void setDimensions(const Box<dim,T> & box, const size_t (&div)[dim], const size_t pad)
1369 {
1370 Matrix<dim,T> mat;
1371 mat.identity();
1372 t.setTransform(mat,box.getP1());
1373 this->box = box;
1374 Initialize(pad,div);
1375 this->cell_shift = 0;
1376 }
1377
1378 /*! \brief Set the cell decomposition parameters + the nested
1379 *
1380 * \param box Domain to decompose
1381 * \param div array with the number of cells on each dimensions
1382 * \param div2 array with the number of cells on each dimension for the nested decomposer
1383 * \param mat transformation matrix the cell space is transformed by p' = A * p
1384 * \param orig origin of the cell decomposition
1385 * \param pad padding cell
1386 *
1387 */
1388 inline void setDimensions(const Box<dim,T> & box, const size_t (&div)[dim], const size_t (&div2)[dim], const Matrix<dim,T> & mat, const size_t pad, Point<dim,long int> cell_shift)
1389 {
1390 t.setTransform(mat,box.getP1());
1391 this->box = box;
1392
1393 Initialize(pad,div);
1394
1395 // The nested cell is big div2 + 2*off
1396
1397 size_t div_with_off[dim];
1398
1399 for(size_t i = 0 ; i < dim ; i++)
1400 div_with_off[i] = div2[i] + 2*off[i];
1401
1402 gr_cell2.setDimensions(div_with_off);
1403
1404 for (size_t i = 0 ; i < dim ; i++)
1405 this->cell_shift.get(i) = cell_shift.get(i) - off[i];
1406 }
1407
1408 /*! \brief Set the cell decomposition parameters
1409 *
1410 * \param box Domain to decompose
1411 * \param div array with the number of cells on each dimensions
1412 * \param mat transformation matrix the cell space is transformed by p' = A * p
1413 * \param orig origin of the cell decomposition
1414 * \param pad padding cell
1415 *
1416 */
1417 inline void setDimensions(const Box<dim,T> & box, const size_t (&div)[dim], const Matrix<dim,T> & mat, const size_t pad)
1418 {
1419 t.setTransform(mat,box.getP1());
1420 this->box = box;
1421 Initialize(pad,div);
1422 this->cell_shift = 0;
1423 }
1424
1425
1426 /*! \brief Set a consistent cell decomposer
1427 *
1428 * When we want to create a new extended CellDecomposer consistent with the old one
1429 * this function can be used to guarantee such costrain.
1430 * With consistent we mean that for each cell of the old CellDecomposer exist
1431 * a Cell in the new CellDecomposer that perfectly overlap
1432 *
1433 * \param cd OLD CellDecomposer
1434 * \param cell_extension extension of the CellDecomposer in term of cells to add in each directions (like Ghost)
1435 *
1436 */
1437 inline void setDimensions(const CellDecomposer_sm<dim,T,transform> & cd, const Box<dim,size_t> & cell_extension)
1438 {
1439 this->cell_shift = 0;
1440
1441 // Get the space transformation
1442
1443 t.setTransform(cd.getMat(),cd.getOrig());
1444
1445 // The domain is equivalent to the old one
1446 this->box = cd.box;
1447
1448 // The padding must be calculated
1449
1450 size_t pad = 0;
1451
1452 for (size_t i = 0 ; i < dim ; i++)
1453 {
1454 if (pad < cell_extension.getLow(i))
1455 pad = cell_extension.getLow(i);
1456 else if (pad > cell_extension.getHigh(i))
1457 pad = cell_extension.getHigh(i);
1458 }
1459
1460 // We have to give the old division
1461
1462 size_t sz_div[dim];
1463
1464 for (size_t i = 0 ; i < dim ; i++)
1465 sz_div[i] = cd.gr_cell.size(i) - 2*cd.off[i];
1466
1467 Initialize(pad,sz_div);
1468 }
1469
1470 /*! \brief Constructor
1471 *
1472 * \param box Space where is defined the cell list
1473 * \param div Reference array to the number of divisions on each dimensions
1474 * \param mat Transformation matrix, the point is transformed as p' = mat * p
1475 * \param pad cell padding
1476 *
1477 * Example for div = {6,6} and pad = 1
1478 *
1479 * \verbatim
1480 * +-----------------------+
1481 * |p |p |p |p |p |p |p |p |
1482 * +-----------------------+
1483 * |p | | | | | | |p |
1484 * +-----------------------+
1485 * |p | | | | | | |p |
1486 * +-----------------------+
1487 * |p | | | | | | |p |
1488 * +-----------------------+
1489 * |p |9 | | | | | |p |
1490 * +-----------------------+
1491 * |p |p |p |p |p |p |p |p |
1492 * +-----------------------+
1493 * \endverbatim
1494 *
1495 * Cell with p are padding cell cell that are around but external the box, the cell number 9 that
1496 * is at the origin of the box is identified with 9
1497 *
1498 */
1499 CellDecomposer_sm(const SpaceBox<dim,T> & box, const size_t (&div)[dim], Matrix<dim,T> & mat, const size_t pad)
1500 :t(Matrix<dim,T>::identity(),box.getP1()),box(box),gr_cell()
1501 {
1502 Initialize(pad);
1503 }
1504
1505 /*! \brief Constructor
1506 *
1507 * \param box Space where is defined the cell list (it is assumed p1 = {0, .... 0})
1508 * \param div Reference array to the number of divisions on each dimensions
1509 * \param pad cell padding
1510 *
1511 * Example for div = {7,7} and pad = 1
1512 *
1513 * \verbatim
1514 * +-----------------------+
1515 * |p |p |p |p |p |p |p |p |
1516 * +-----------------------+
1517 * |p | | | | | | |p |
1518 * +-----------------------+
1519 * |p | | | | | | |p |
1520 * +-----------------------+
1521 * |p | | | | | | |p |
1522 * +-----------------------+
1523 * |p |9 | | | | | |p |
1524 * +-----------------------+
1525 * |p |p |p |p |p |p |p |p |
1526 * +-----------------------+
1527 *
1528 * \endverbatim
1529 *
1530 * Cell with p are padding cell cell that are around but external the box, the cell number 9 that
1531 * is at the origin of the box is identified with 9
1532 *
1533 */
1534 CellDecomposer_sm(const SpaceBox<dim,T> & box, const size_t (&div)[dim], const size_t pad)
1535 :t(Matrix<dim,T>::identity(),box.getP1()),box(box),gr_cell(div)
1536 {
1537 Initialize(pad,div);
1538 }
1539
1540 /*! \brief Constructor for a consistent CellDecomposer construction
1541 *
1542 * \param cd Old CellDecomposer
1543 * \param ext Extension box
1544 *
1545 *
1546 */
1547 CellDecomposer_sm(const CellDecomposer_sm<dim,T,transform> & cd, Box<dim,size_t> & ext)
1548 :t(Matrix<dim,T>::identity(),cd.getOrig())
1549 {
1550 setDimensions(cd,ext);
1551 }
1552
1553
1554 //! default constructor
1555 CellDecomposer_sm()
1556 :t(Matrix<dim,T>::identity(),Point<dim,T>::zero_p()),tot_n_cell(0)
1557 {
1558
1559 }
1560
1561 /*! \brief Return the box that represent the cell
1562 *
1563 * Can be considered the spacing between vertices of the cells
1564 *
1565 * \return the box
1566 *
1567 */
1568 inline const Box<dim,T> & getCellBox() const
1569 {
1570 return box_unit;
1571 }
1572
1573 /*! \brief Get the transformation Matrix of the cell decomposer
1574 *
1575 * \return the transformation Matrix
1576 *
1577 */
1578 inline const Matrix<dim,T> & getMat() const
1579 {
1580 return t.getMat();
1581 }
1582
1583 /*! \brief Get the origin of the cell decomposer
1584 *
1585 * \return the origin
1586 *
1587 */
1588 inline const Point<dim,T> & getOrig() const
1589 {
1590 return t.getOrig();
1591 }
1592
1593 /*! \brief Convert a Box in the domain space into cell units (Negative contour included, positive contour excluded)
1594 *
1595 * Given the following
1596 *
1597 * \verbatim
1598 *
1599 +-----+-----+-----+-----+-----+-----+ (1.0. 1.0) Domain box
1600 | | | | | | |
1601 | | | | | | |
1602 | +-----------------+ | | |
1603 +-----+-----+-----+-----+-----+-----+
1604 | | | | | | | | |
1605Box "b" <-----------------+ | | | | | | Grid (7, 6)
1606(0.1 , 0.42) | | | | | | | | | Cells (6,5)
1607(0.64, 0.85) +-----+-----+-----+-----+-----+-----+
1608 | | | | | | | | |
1609 | | | | | | | | |
1610 | +-----------------+ | | |
1611 +-----+-----+-----+-----+-----+-----+
1612 | | | | | | |
1613 | | | | | | |
1614 +-----+-----+-----+-----+-----+-----+
1615 | | | | | | |
1616 | | | | | | |
1617 | | | | | | |
1618 +-----+-----+-----+-----+-----+-----+
1619 (0.0, 0.0)
1620
1621
1622 + = grid points
1623
1624 \verbatim
1625
1626 It return a Box with P1 = (0,2), P2 = (3,4)
1627
1628 *
1629 * \param b Box in domain space
1630 * \param bc boundary conditions
1631 *
1632 * \return Box in grid units, if P2 < P1 the box does not include any grid points
1633 *
1634 */
1635 inline Box<dim,long int> convertDomainSpaceIntoCellUnits(const Box<dim,T> & b_d, const size_t (& bc)[dim]) const
1636 {
1637 Box<dim,long int> g_box;
1638 Box<dim,T> b = b_d;
1639 b -= getOrig();
1640
1641 // Convert b into grid units
1642 b /= getCellBox().getP2();
1643
1644 // Considering that we are interested in a box decomposition of the space
1645 // where each box does not intersect any other boxes in the decomposition we include the negative
1646 // countour and exclude the positive one. So ceilP1 do the job for P1 while ceilP2 - 1
1647 // do the job for P2
1648
1649 b.floorP1();
1650 b.ceilP2();
1651
1652 g_box = b;
1653
1654 // Translate the box by the offset
1655
1656 for (size_t i = 0 ; i < dim ; i++)
1657 {
1658 g_box.setLow(i,g_box.getLow(i) + off[i]);
1659 g_box.setHigh(i,g_box.getHigh(i) + off[i]);
1660 }
1661
1662 // on the other hand with non periodic boundary condition, the positive border of the
1663 // sub-domain at the edge of the domain must be included
1664
1665 Point<dim,size_t> p_move;
1666
1667 for (size_t i = 0 ; i < dim ; i++)
1668 {
1669 // we are at the positive border (We are assuming that there are not rounding error in the decomposition)
1670 if (b_d.getHigh(i) == box.getHigh(i))
1671 {
1672 if (bc[i] == NON_PERIODIC)
1673 {
1674 // Here the positive boundary is included
1675 g_box.setHigh(i,gr_cell.size(i) - off[i]);
1676 }
1677 else
1678 {
1679 // Carefull in periodic gr_cell is one bigger than the non-periodic
1680 // and the positive boundary is excluded
1681 g_box.setHigh(i,gr_cell.size(i)-1 - off[i]);
1682 }
1683 }
1684
1685
1686 if (b_d.getLow(i) == box.getHigh(i))
1687 {
1688 if (bc[i] == NON_PERIODIC)
1689 {
1690 // The instruction is the same but the meaning is different
1691 // for this reason there is anyway a branch
1692 // Here the border is not included
1693 g_box.setLow(i,gr_cell.size(i) - off[i]);
1694 }
1695 else
1696 {
1697 // Carefull in periodic gr_cell is one bigger than the non-periodic
1698 // Here the border is included
1699 g_box.setLow(i,gr_cell.size(i) - off[i]);
1700 }
1701 }
1702
1703 /////////// Low boundary
1704
1705 // we are at the positive border (We are assuming that there are not rounding error in the decomposition)
1706 /* coverity[copy_paste_error] */
1707 if (b_d.getHigh(i) == box.getLow(i))
1708 g_box.setHigh(i,off[i]);
1709
1710
1711 if (b_d.getLow(i) == box.getLow(i))
1712 g_box.setLow(i,off[i]);
1713 }
1714
1715 return g_box;
1716 }
1717
1718
1719 /*! \brief Convert a Box in the domain space into grid units (Negative contour included, positive contour excluded)
1720 *
1721 * Given the following
1722 *
1723 * \verbatim
1724 *
1725 +-----+-----+-----+-----+-----+-----+ (1.0. 1.0) Domain box
1726 | | | | | | |
1727 | | | | | | |
1728 | +-----------------+ | | |
1729 +-----+-----+-----+-----+-----+-----+
1730 | | | | | | | | |
1731Box "b" <-----------------+ | | | | | | Grid (7, 6)
1732(0.1 , 0.42) | | | | | | | | |
1733(0.64, 0.85) +-----+-----+-----+-----+-----+-----+
1734 | | | | | | | | |
1735 | | | | | | | | |
1736 | +-----------------+ | | |
1737 +-----+-----+-----+-----+-----+-----+
1738 | | | | | | |
1739 | | | | | | |
1740 +-----+-----+-----+-----+-----+-----+
1741 | | | | | | |
1742 | | | | | | |
1743 | | | | | | |
1744 +-----+-----+-----+-----+-----+-----+
1745 (0.0, 0.0)
1746
1747
1748 + = grid points
1749
1750 \verbatim
1751
1752 It return a Box with P1 = (1,3), P2 = (3,4)
1753
1754 *
1755 * \param b Box in domain space
1756 * \param bc boundary conditions
1757 *
1758 * \return Box in grid units, if P2 < P1 the box does not include any grid points
1759 *
1760 */
1761 inline Box<dim,long int> convertDomainSpaceIntoGridUnits(const Box<dim,T> & b_d, const size_t (& bc)[dim]) const
1762 {
1763 Box<dim,long int> g_box;
1764 Box<dim,T> b = b_d;
1765 b -= getOrig();
1766
1767 // Convert b into grid units
1768 b /= getCellBox().getP2();
1769
1770 // Considering that we are interested in a box decomposition of the space
1771 // where each box does not intersect any other boxes in the decomposition we include the negative
1772 // countour and exclude the positive one. So ceilP1 do the job for P1 while ceilP2 - 1
1773 // do the job for P2
1774
1775 b.ceilP1();
1776
1777 // (we do -1 later)
1778 b.ceilP2();
1779 for (size_t i = 0 ; i < dim ; i++) {b.setHigh(i,b.getHigh(i)-1);}
1780
1781 g_box = b;
1782
1783 // on the other hand with non periodic boundary condition, the positive border of the
1784 // sub-domain at the edge of the domain must be included
1785
1786 Point<dim,size_t> p_move;
1787
1788 for (size_t i = 0 ; i < dim ; i++)
1789 {
1790 // we are at the positive border (We are assuming that there are not rounding error in the decomposition)
1791 if (b_d.getHigh(i) == box.getHigh(i))
1792 {
1793 if (bc[i] == NON_PERIODIC)
1794 {
1795 // Here the positive boundary is included
1796 g_box.setHigh(i,gr_cell.size(i) - off[i]);
1797 }
1798 else
1799 {
1800 // Carefull in periodic gr_cell is one bigger than the non-periodic
1801 // and the positive boundary is excluded
1802 g_box.setHigh(i,gr_cell.size(i)-1 - off[i]);
1803 }
1804 }
1805
1806
1807 if (b_d.getLow(i) == box.getHigh(i))
1808 {
1809 if (bc[i] == NON_PERIODIC)
1810 {
1811 // The instruction is the same but the meaning is different
1812 // for this reason there is anyway a branch
1813 // Here the border is not included
1814 g_box.setLow(i,gr_cell.size(i) - off[i]);
1815 }
1816 else
1817 {
1818 // Carefull in periodic gr_cell is one bigger than the non-periodic
1819 // Here the border is included
1820 g_box.setLow(i,gr_cell.size(i) - off[i]);
1821 }
1822 }
1823
1824 /////////// Low boundary
1825
1826 // we are at the positive border (We are assuming that there are not rounding error in the decomposition)
1827 /* coverity[copy_paste_error] */
1828 if (b_d.getHigh(i) == box.getLow(i))
1829 g_box.setHigh(i,off[i]);
1830
1831
1832 if (b_d.getLow(i) == box.getLow(i))
1833 g_box.setLow(i,off[i]);
1834 }
1835
1836 return g_box;
1837 }
1838
1839
1840 /*! \brief Convert a Box from grid units into domain space (Negative contour included, positive contour excluded)
1841 *
1842 * Given the following
1843 *
1844 * \verbatim
1845
1846 +-----+-----+-----+-----+-----+-----+ (1.0. 1.0) Domain box
1847 | | | | | | |
1848 | | | | | | |
1849 | | | | | | |
1850 +-----+-----+-----+-----+----{2}----+
1851 | | | | | | |
1852 | | | | | | | |
1853 | | | | | | |
1854 +-----+-----+-----+-----+-----+-----+
1855 | | | | | | |
1856 | | | | | | |
1857 | | | | | |
1858 +-----+----{1}----+-----+-----+-----+
1859 | | | | | | |
1860 | | | | | | |
1861 +-----+-----+-----+-----+-----+-----+
1862 | | | | | | |
1863 | | | | | | |
1864 | | | | | | |
1865 +-----+-----+-----+-----+-----+-----+
1866 (0.0, 0.0)
1867
1868
1869 + = grid points
1870
1871 \verbatim
1872
1873 Giving a box P1 = (2,2), P2 = (5,4)
1874 it return a box (0.333333,0.33333333) (0.8333333,0.8)
1875
1876 *
1877 * \param b Box in grid units
1878 *
1879 * \return the box in domain space
1880 *
1881 */
1882 inline Box<dim,T> convertCellUnitsIntoDomainSpace(const Box<dim,long int> & b_d) const
1883 {
1884 Box<dim,T> be;
1885
1886 for (size_t i = 0 ; i < dim ; i++)
1887 {
1888 if ((long int)gr_cell.size(i) - (long int)off[i] == b_d.getLow(i))
1889 be.setLow(i,box.getHigh(i));
1890 else if ((long int)off[i] == b_d.getLow(i))
1891 be.setLow(i,box.getLow(i));
1892 else
1893 be.setLow(i,(b_d.getLow(i) - off[i]) * box_unit.getP2()[i] + box.getLow(i));
1894
1895 if ((long int)gr_cell.size(i) - (long int)off[i] == b_d.getHigh(i))
1896 be.setHigh(i,box.getHigh(i));
1897 else if ((long int)off[i] == b_d.getHigh(i))
1898 be.setHigh(i,box.getLow(i));
1899 else
1900 be.setHigh(i,(b_d.getHigh(i) - off[i]) * box_unit.getP2()[i] + box.getLow(i));
1901 }
1902
1903 return be;
1904 }
1905
1906
1907 /*! \brief Convert a Box from grid units into domain space (Negative contour included, positive contour excluded)
1908 *
1909 * Given the following
1910 *
1911 * \verbatim
1912
1913 +-----+-----+-----+-----+-----+-----+ (1.0. 1.0) Domain box
1914 | | | | | | |
1915 | | *-----------------------* |
1916 | | | | | | | | |
1917 +-----+--|--+-----+-----+----{2}-|--+
1918 | | | | | | | | |
1919 | | | | | | | | |
1920 | | | | | | | | |
1921 +-----+--|--+-----+-----+-----+--|--+
1922 | | | | | | | | |
1923 | | | | | | | | |
1924 | | | | | | | | |
1925 +-----+--|-{1}----+-----+-----+--|--+
1926 | | | | | | | | |
1927 | | *-----------------------* |
1928 | | | | | | |
1929 +-----+-----+-----+-----+-----+-----+
1930 | | | | | | |
1931 | | | | | | |
1932 | | | | | | |
1933 +-----+-----+-----+-----+-----+-----+
1934 (0.0, 0.0)
1935
1936
1937 + = grid points
1938
1939 \verbatim
1940
1941 Giving a box P1 = (2,2), P2 = (5,4)
1942 it return a box (0.25,0.25) (0.916666,0.9)
1943
1944 *
1945 * \param b Box in grid units
1946 *
1947 * The box is never allowed to be bigger than the domain and is always cropped
1948 * by the domain size
1949 *
1950 * \return the box in domain space
1951 *
1952 */
1953 inline Box<dim,T> convertCellUnitsIntoDomainSpaceMiddle(const Box<dim,long int> & b_d) const
1954 {
1955 Box<dim,T> be;
1956
1957 for (size_t i = 0 ; i < dim ; i++)
1958 {
1959 if ((long int)gr_cell.size(i) - (long int)off[i] == b_d.getLow(i))
1960 {be.setLow(i,box.getHigh(i));}
1961 else if ((long int)off[i] == b_d.getLow(i))
1962 {be.setLow(i,box.getLow(i));}
1963 else
1964 {be.setLow(i,(b_d.getLow(i) - off[i]) * box_unit.getP2()[i] + box.getLow(i) - box_unit.getP2()[i] / 2.0);}
1965
1966 if ((long int)gr_cell.size(i) - (long int)off[i] == b_d.getHigh(i))
1967 {be.setHigh(i,box.getHigh(i));}
1968 else if ((long int)off[i] == b_d.getHigh(i))
1969 {be.setHigh(i,box.getLow(i));}
1970 else
1971 {be.setHigh(i,(b_d.getHigh(i) - off[i]) * box_unit.getP2()[i] + box.getLow(i) + box_unit.getP2()[i] / 2.0);}
1972 }
1973
1974 return be;
1975 }
1976
1977 /*! \brief Return the number of divisions of the Cell Decomposer (including padding)
1978 *
1979 * \return the number of divisions
1980 *
1981 */
1982 const size_t (& getDiv() const)[dim]
1983 {
1984 return gr_cell.getSize();
1985 }
1986
1987
1988 /*! \brief Return the number of divisions of the Cell Decomposer (without padding)
1989 *
1990 * \return the number of divisions
1991 *
1992 */
1993 const size_t (& getDivWP() const)[dim]
1994 {
1995 for (size_t i = 0 ; i < dim ; i++)
1996 {div_wp[i] = gr_cell.getSize()[i] - 2*getPadding(i);}
1997
1998 return div_wp;
1999 }
2000
2001 /*! \brief Return the domain where the CellDecomposer is defined
2002 *
2003 * \return The domain of the remote machine
2004 *
2005 */
2006 const Box<dim,T> & getDomain() const
2007 {
2008 return box;
2009 }
2010
2011 /*! \brief it swap the content of two Cell Decomposer
2012 *
2013 * \param cd CellDecomposer to swap with
2014 *
2015 */
2016 inline void swap(CellDecomposer_sm<dim,T,transform> & cd)
2017 {
2018 // swap all the members
2019 p_middle.swap(cd.p_middle);
2020
2021 // Point transformation before get the Cell object (useful for example to shift the cell list)
2022 transform t_t = t;
2023 t = cd.t;
2024 cd.t = t_t;
2025
2026 // Total number of cell
2027 size_t tot_n_cell_t = tot_n_cell;
2028 tot_n_cell = cd.tot_n_cell;
2029 cd.tot_n_cell = tot_n_cell_t;
2030
2031 box.swap(cd.box);
2032 box_unit.swap(cd.box_unit);
2033 gr_cell.swap(cd.gr_cell);
2034 gr_cell2.swap(cd.gr_cell2);
2035
2036 for (size_t i = 0 ; i < dim ; i++)
2037 {
2038 size_t off_t = off[i];
2039 off[i] = cd.off[i];
2040 cd.off[i] = off_t;
2041
2042 size_t cs_t = cell_shift.get(i);
2043 cell_shift.get(i) = cd.cell_shift.get(i);
2044 cd.cell_shift.get(i) = cs_t;
2045 }
2046 }
2047
2048 /*! \brief Check that the CellDecomposer is the same
2049 *
2050 * \param cd Cell decomposer
2051 *
2052 * \return true if the two CellDecomposer are the same
2053 *
2054 */
2055 inline bool operator==(const CellDecomposer_sm<dim,T,transform> & cd)
2056 {
2057 if (meta_compare<Point<dim,T>>::meta_compare_f(p_middle,cd.p_middle) == false)
2058 return false;
2059
2060 if (t != cd.t)
2061 return false;
2062
2063 if (tot_n_cell != cd.tot_n_cell)
2064 return false;
2065
2066 if (box != cd.box)
2067 return false;
2068
2069 if (box_unit != cd.box_unit)
2070 return false;
2071
2072 if (gr_cell != cd.gr_cell)
2073 return false;
2074
2075 if (gr_cell2 != cd.gr_cell2)
2076 return false;
2077
2078 for (size_t i = 0 ; i < dim ; i++)
2079 {
2080 if (off[i] != cd.off[i])
2081 return false;
2082
2083 if (cell_shift.get(i) != cd.cell_shift.get(i))
2084 return false;
2085 }
2086
2087 return true;
2088 }
2089
2090 /*! \brief Check that the CellDecomposer is the same
2091 *
2092 * \param cd Cell decomposer
2093 *
2094 * \return true if the two CellDecomposer is different
2095 *
2096 */
2097 inline bool operator!=(const CellDecomposer_sm<dim,T,transform> & cd)
2098 {
2099 return ! this->operator==(cd);
2100 }
2101
2102 /*! \brief Return the number of padding cells of the Cell decomposer
2103 *
2104 * \param i dimension
2105 *
2106 * \return the number of padding cells
2107 *
2108 */
2109 size_t getPadding(size_t i) const
2110 {
2111 return off[i];
2112 }
2113
2114 /*! \brief Return the number of padding cells of the Cell decomposer as an array
2115 *
2116 *
2117 * \return the number of padding cells
2118 *
2119 */
2120 size_t (& getPadding())[dim]
2121 {
2122 return off;
2123 }
2124
2125 /*! \brief Return the index of the first cell in the domain
2126 *
2127 * This function make sense if the CellDecomposer has been
2128 * created with setDimensions with div and div2
2129 *
2130 * \return the first domain cell
2131 *
2132 */
2133 grid_key_dx<dim> getStartDomainCell() const
2134 {
2135 grid_key_dx<dim> key;
2136
2137 for (size_t i = 0 ; i < dim ; i++)
2138 {
2139 key.set_d(i, cell_shift.get(i));
2140 }
2141
2142 return key;
2143 }
2144
2145 /*! \brief Return the index of the last cell in the domain
2146 *
2147 * This function make sense if the CellDecomposer has been
2148 * created with setDimensions with div and div2
2149 *
2150 * \return the last domain cell
2151 *
2152 */
2153 grid_key_dx<dim> getStopDomainCell() const
2154 {
2155 grid_key_dx<dim> key;
2156
2157 for (size_t i = 0 ; i < dim ; i++)
2158 {
2159 key.set_d(i,cell_shift.get(i) + gr_cell2.size(i) - 2*getPadding(i) - 1);
2160 }
2161
2162 return key;
2163 }
2164
2165 /*! \brief Return the cell shift
2166 *
2167 * \return the cell shifting
2168 *
2169 */
2170 grid_key_dx<dim> getShift() const
2171 {
2172 grid_key_dx<dim> k;
2173
2174 for (size_t i = 0 ; i < dim ; i++)
2175 k.set_d(i,cell_shift.get(i));
2176
2177 return k;
2178 }
2179
2180 /*! \brief Return the grid structure of the internal cell list
2181 *
2182 * \return the grid information
2183 *
2184 */
2185 const grid_sm<dim,void> & getInternalGrid() const
2186 {
2187 return gr_cell2;
2188 }
2189};
2190
2191
2192#endif /* CELLDECOMPOSER_HPP_ */
2193