1/*
2 * CellList.hpp
3 *
4 * Created on: Mar 21, 2015
5 * Author: Pietro Incardona
6 */
7
8#ifndef CELLLIST_HPP_
9#define CELLLIST_HPP_
10
11#include "CellList_def.hpp"
12#include "Vector/map_vector.hpp"
13#include "CellDecomposer.hpp"
14#include "Space/SpaceBox.hpp"
15#include "util/mathutil.hpp"
16#include "CellNNIterator.hpp"
17#include "Space/Shape/HyperCube.hpp"
18#include "CellListNNIteratorRadius.hpp"
19#include <unordered_map>
20
21#include "CellListIterator.hpp"
22#include "ParticleIt_Cells.hpp"
23#include "ParticleItCRS_Cells.hpp"
24#include "util/common.hpp"
25
26#include "NN/Mem_type/MemFast.hpp"
27#include "NN/Mem_type/MemBalanced.hpp"
28#include "NN/Mem_type/MemMemoryWise.hpp"
29#include "NN/CellList/NNc_array.hpp"
30#include "cuda/CellList_cpu_ker.cuh"
31
32//! Wrapper of the unordered map
33template<typename key,typename val>
34class wrap_unordered_map: public std::unordered_map<key,val>
35{
36};
37
38#ifdef HAVE_LIBQUADMATH
39
40#include <boost/multiprecision/float128.hpp>
41
42
43//! Wrapper of the unordered map
44template<typename val>
45class wrap_unordered_map<boost::multiprecision::float128,val>
46{
47};
48
49#endif
50
51//! Point at witch the cell do a reallocation (it should the the maximum for all the implementations)
52#define CELL_REALLOC 16ul
53
54#define STARTING_NSLOT 16
55
56/*! \brief Calculate the the Neighborhood for symmetric interactions CSR scheme
57 *
58 * \param cNN calculated cross neighborhood
59 * \param div Number of divisions in each direction
60 *
61 */
62template<unsigned int dim> void NNcalc_csr(openfpm::vector<std::pair<grid_key_dx<dim>,grid_key_dx<dim>>> & cNN)
63{
64 // Calculate the NNc_full array, it is a structure to get the neighborhood array
65
66 // compile-time array {0,0,0,....} {2,2,2,...} {1,1,1,...}
67
68 typedef typename generate_array<size_t,dim, Fill_zero>::result NNzero;
69 typedef typename generate_array<size_t,dim, Fill_two>::result NNtwo;
70
71 // Generate the sub-grid iterator
72
73 size_t div[dim];
74
75 // Calculate the divisions
76
77 for (size_t i = 0 ; i < dim ; i++)
78 div[i] = 4;
79
80 grid_sm<dim,void> gs(div);
81 grid_key_dx_iterator_sub<dim> gr_sub3(gs,NNzero::data,NNtwo::data);
82
83 grid_key_dx<dim> src_; // Source cell
84 for (size_t i = 0; i < dim; i++)
85 src_.set_d(i,1);
86
87 size_t middle = gs.LinId(src_);
88
89 // Calculate the symmetric crs array
90 while (gr_sub3.isNext())
91 {
92 auto dst = gr_sub3.get();
93 grid_key_dx<dim> src = src_;
94
95 if ((long int)middle > gs.LinId(dst))
96 {
97 ++gr_sub3;
98 continue;
99 }
100
101 // Here we adjust src and dst to be in the positive quadrant
102
103 for (size_t i = 0 ; i < dim; i++)
104 {
105 if (dst.get(i) == 0)
106 {
107 src.set_d(i,src.get(i) + 1);
108 dst.set_d(i,dst.get(i) + 1);
109 }
110 }
111
112 src -= src_;
113 dst -= src_;
114
115 cNN.add(std::pair<grid_key_dx<dim>,grid_key_dx<dim>>(src,dst));
116
117 ++gr_sub3;
118
119 }
120};
121
122/*! \brief Calculate the the Neighborhood for symmetric interactions
123 *
124 * \param cNN calculated cross neighborhood
125 * \param div Number of divisions in each direction
126 *
127 */
128template<unsigned int dim> void NNcalc_sym(openfpm::vector<grid_key_dx<dim>> & cNN)
129{
130 // Calculate the NNc_full array, it is a structure to get the neighborhood array
131
132 // compile-time array {0,0,0,....} {2,2,2,...} {1,1,1,...}
133
134 typedef typename generate_array<size_t,dim, Fill_zero>::result NNzero;
135 typedef typename generate_array<size_t,dim, Fill_two>::result NNtwo;
136
137 // Generate the sub-grid iterator
138
139 size_t div[dim];
140
141 // Calculate the divisions
142
143 for (size_t i = 0 ; i < dim ; i++)
144 div[i] = 4;
145
146 grid_sm<dim,void> gs(div);
147 grid_key_dx_iterator_sub<dim> gr_sub3(gs,NNzero::data,NNtwo::data);
148
149 grid_key_dx<dim> src_; // Source cell
150 for (size_t i = 0; i < dim; i++)
151 src_.set_d(i,1);
152
153 size_t middle = gs.LinId(src_);
154
155 // Calculate the symmetric array
156 while (gr_sub3.isNext())
157 {
158 auto dst = gr_sub3.get();
159
160 if ((long int)middle > gs.LinId(dst))
161 {
162 ++gr_sub3;
163 continue;
164 }
165
166 cNN.add(dst - src_);
167
168 ++gr_sub3;
169
170 }
171};
172
173/*! \brief Calculate the Neighborhood cells
174 *
175 * \param cNN calculated cross neighborhood
176 * \param div Number of divisions in each direction
177 *
178 */
179template<unsigned int dim> void NNcalc_full(openfpm::vector<grid_key_dx<dim>> & cNN)
180{
181 // Calculate the NNc_full array, it is a structure to get the neighborhood array
182
183 // compile-time array {0,0,0,....} {2,2,2,...} {1,1,1,...}
184
185 typedef typename generate_array<size_t,dim, Fill_zero>::result NNzero;
186 typedef typename generate_array<size_t,dim, Fill_two>::result NNtwo;
187
188 // Generate the sub-grid iterator
189
190 size_t div[dim];
191
192 // Calculate the divisions
193
194 for (size_t i = 0 ; i < dim ; i++)
195 {div[i] = 4;}
196
197 grid_sm<dim,void> gs(div);
198 grid_key_dx_iterator_sub<dim> gr_sub3(gs,NNzero::data,NNtwo::data);
199
200 grid_key_dx<dim> src_; // Source cell
201 for (size_t i = 0; i < dim; i++)
202 src_.set_d(i,1);
203
204 // Calculate the symmetric crs array
205 while (gr_sub3.isNext())
206 {
207 auto dst = gr_sub3.get();
208
209 cNN.add(dst - src_);
210
211 ++gr_sub3;
212 }
213};
214
215
216/*! Calculate the neighborhood cells based on the radius
217 *
218 * \note To the calculated neighborhood cell you have to add the id of the central cell
219 *
220 \verbatim
221 +-----------------------+
222 |p |p |p |p |p |p |p |p |
223 +-----------------------+
224 |p | | | | | | |p |
225 +-----------------------+
226 |p | | |7 |8 |9 | |p |
227 +-----------------------+
228 |p | | |-1|0 |1 | |p |
229 +-----------------------+
230 |p |9 | |-9|-8|-7| |p |
231 +-----------------------+
232 |p |p |p |p |p |p |p |p |
233 +-----------------------+
234 \endverbatim
235 *
236 * The number indicate the cell id calculated
237 *
238 * -9,-8,-7,-1,0,1,7,8,9
239 *
240 * The cell 0 has id = 22 in the big cell matrix, so to calculate the
241 * neighborhood cells you have to sum the id of the center cell
242 *
243 * 13,14,15,21,22,23,29,30,31
244 *
245 * \param r_cut Cutoff-radius
246 * \param NNcell vector containing the neighborhood cells ids
247 *
248 */
249template<unsigned int dim, typename T>
250void NNcalc_rad(T r_cut, openfpm::vector<long int> & NNcell, const Box<dim,T> & cell_box, const grid_sm<dim,void> & gs)
251{
252 size_t n_cell[dim];
253 size_t n_cell_mid[dim];
254
255 Point<dim,T> spacing = cell_box.getP2();
256
257 for (size_t i = 0 ; i < dim ; i++)
258 {
259 n_cell[i] = 2*(std::ceil(r_cut / spacing.get(i)))+1;
260 n_cell_mid[i] = n_cell[i] / 2;
261 }
262
263 grid_sm<dim,void> gsc(n_cell);
264 grid_key_dx_iterator<dim> gkdi(gsc);
265
266 Box<dim,T> cell_zero;
267
268 for (unsigned int i = 0 ; i < dim ; i++)
269 {
270 cell_zero.setLow(i,n_cell_mid[i]*spacing.get(i));
271 cell_zero.setHigh(i,(n_cell_mid[i]+1)*spacing.get(i));
272 }
273
274 NNcell.clear();
275 while (gkdi.isNext())
276 {
277 auto key = gkdi.get();
278
279 Box<dim,T> cell;
280
281 for (unsigned int i = 0 ; i < dim ; i++)
282 {
283 cell.setLow(i,key.get(i)*spacing.get(i));
284 cell.setHigh(i,(key.get(i)+1)*spacing.get(i));
285 }
286
287 // here we check if the cell is in the radius.
288 T min_distance = cell.min_distance(cell_zero);
289 if (min_distance > r_cut)
290 {++gkdi;continue;}
291
292 for (size_t i = 0 ; i < dim ; i++)
293 {key.set_d(i,key.get(i) - n_cell_mid[i]);}
294
295 NNcell.add(gs.LinId(key));
296
297 ++gkdi;
298 }
299}
300
301/* NOTE all the implementations
302 *
303 * has complexity O(1) in getting the cell id and the elements in a cell
304 * but with different constants
305 *
306 */
307
308/*! \brief Class for FAST cell list implementation
309 *
310 * This class implement the FAST cell list, fast but memory
311 * expensive. The memory allocation is (M * N_cell_max)*sizeof(ele) + M*8
312 *
313 * * M = number of cells
314 * * N_cell_max = maximum number of elements in a cell
315 * * ele = element the structure is storing
316 *
317 * \note Because N_cell_max >= N/M then M * N_cell_max >= O(N)
318 *
319 * \warning Not not use for high asymmetric distribution
320 *
321 * Example of a 2D Cell list 6x6 structure with padding 1 without shift, cell indicated with p are padding cell
322 * the origin of the cell or point (0,0) is marked with cell number 9
323 *
324 * \verbatim
325 * +-----------------------+
326 * |p |p |p |p |p |p |p |p |
327 * +-----------------------+
328 * |p | | | | | | |p |
329 * +-----------------------+
330 * |p | | | | | | |p |
331 * +-----------------------+
332 * |p | | | | | | |p |
333 * +-----------------------+
334 * |p |9 | | | | | |p |
335 * +-----------------------+
336 * |p |p |p |p |p |p |p |p |
337 * +-----------------------+
338 * \endverbatim
339 *
340 *
341 * \tparam dim Dimensionality of the space
342 * \tparam T type of the space float, double ...
343 * \tparam base Base structure that store the information
344 *
345 * ### Declaration of a cell list
346 * \snippet CellList_test.hpp Declare a cell list
347 * ### Usage of cell list [CellS == CellList<3,double,FAST>]
348 * \snippet CellList_test.hpp Usage of cell list
349 * ### Remove one particle from each cell
350 * \snippet CellList_test.hpp remove one particle from each cell
351 * ### Usage of the neighborhood iterator
352 * \snippet CellList_test.hpp Usage of the neighborhood iterator
353 *
354 */
355template<unsigned int dim, typename T, typename Mem_type, typename transform = no_transform<dim,T>, typename vector_pos_type = openfpm::vector<Point<dim,T>>>
356class CellList : public CellDecomposer_sm<dim,T,transform>, public Mem_type
357{
358protected:
359 //! The array contain the neighborhood of the cell-id in case of asymmetric interaction
360 //
361 // * * *
362 // * x *
363 // * * *
364
365
366 NNc_array<dim,(unsigned int)openfpm::math::pow(3,dim)> NNc_full;
367
368 //! The array contain the neighborhood of the cell-id in case of symmetric interaction
369 //
370 // * * *
371 // x *
372 //
373// long int NNc_sym[openfpm::math::pow(3,dim)/2+1];
374 NNc_array<dim,(unsigned int)openfpm::math::pow(3,dim)/2+1> NNc_sym;
375
376private:
377
378 //! Caching of r_cutoff radius
379 wrap_unordered_map<T,openfpm::vector<long int>> rcache;
380
381 //! True if has been initialized from CellDecomposer
382 bool from_cd;
383
384 //! Additional information in general (used to understand if the cell-list)
385 //! has been constructed from an old decomposition
386 size_t n_dec;
387
388 //! Cells for the neighborhood radius
389 openfpm::vector<long int> nnc_rad;
390
391
392
393 //! Initialize the structures of the data structure
394 void InitializeStructures(const size_t (& div)[dim], size_t tot_n_cell, size_t slot=STARTING_NSLOT)
395 {
396 Mem_type::init_to_zero(slot,tot_n_cell);
397
398 NNc_full.set_size(div);
399 NNc_full.init_full();
400
401 NNc_sym.set_size(div);
402 NNc_sym.init_sym();
403 }
404
405 void setCellDecomposer(CellDecomposer_sm<dim,T,transform> & cd, const CellDecomposer_sm<dim,T,transform> & cd_sm, const Box<dim,T> & dom_box, size_t pad) const
406 {
407 size_t bc[dim];
408
409 for (size_t i = 0 ; i < dim ; i++)
410 {bc[i] = NON_PERIODIC;}
411
412 Box<dim,long int> bx = cd_sm.convertDomainSpaceIntoCellUnits(dom_box,bc);
413
414 size_t div[dim];
415 size_t div_big[dim];
416
417 for (size_t i = 0 ; i < dim ; i++)
418 {
419 div[i] = bx.getHigh(i) - bx.getLow(i);
420 div_big[i] = cd_sm.getDiv()[i] - 2*cd_sm.getPadding(i);
421 }
422
423 cd.setDimensions(cd_sm.getDomain(),div_big,div, pad, bx.getP1());
424 }
425
426public:
427
428 //! Type of internal memory structure
429 typedef Mem_type Mem_type_type;
430
431 typedef CellNNIteratorSym<dim,CellList<dim,T,Mem_type,transform,vector_pos_type>,vector_pos_type,RUNTIME,NO_CHECK> SymNNIterator;
432
433 //! Object type that the structure store
434 typedef typename Mem_type::local_index_type value_type;
435
436 //! Type of the coordinate space (double float)
437 typedef T stype;
438
439 //!
440 typedef vector_pos_type internal_vector_pos_type;
441
442 /*! \brief Return the underlying grid information of the cell list
443 *
444 * \return the grid infos
445 *
446 */
447 const grid_sm<dim,void> & getGrid()
448 {
449 return CellDecomposer_sm<dim,T,transform>::getGrid();
450 }
451
452 /*! Initialize the cell list from a well-define Cell-decomposer
453 *
454 * In some cases is needed to have a Cell-list with Cells consistent
455 * with a well predefined CellDecomposer. In this case we use this function.
456 * Using this initialization the Cell-list maintain the Cells defined by this
457 * Cell-decomposer consistently
458 *
459 * \param cd_sm Cell-Decomposer
460 * \param dom_box domain box (carefully this is going to be adjusted)
461 * \param pad cell-list padding
462 * \param slot slots for each cell
463 *
464 */
465 void Initialize(CellDecomposer_sm<dim,T,transform> & cd_sm, const Box<dim,T> & dom_box,const size_t pad = 1, size_t slot=STARTING_NSLOT)
466 {
467 size_t bc[dim];
468 for (size_t i = 0 ; i < dim ; i++) {bc[i] = NON_PERIODIC;}
469
470 Box<dim,long int> bx = cd_sm.convertDomainSpaceIntoCellUnits(dom_box,bc);
471
472 setCellDecomposer(*this,cd_sm,dom_box,pad);
473
474 size_t div_w_pad[dim];
475 size_t tot_cell = 1;
476
477 for (size_t i = 0 ; i < dim ; i++)
478 {
479 div_w_pad[i] = bx.getHigh(i) - bx.getLow(i) + 2*pad;
480 tot_cell *= div_w_pad[i];
481 }
482
483 // here we set the cell-shift for the CellDecomposer
484
485 InitializeStructures(div_w_pad,tot_cell);
486
487 // Initialized from CellDecomposer
488 from_cd = true;
489 }
490
491 /*! Initialize the cell list
492 *
493 * \param box Domain where this cell list is living
494 * \param div grid size on each dimension
495 * \param pad padding cell
496 * \param slot maximum number of slot
497 *
498 */
499 void Initialize(const Box<dim,T> & box, const size_t (&div)[dim], const size_t pad = 1, size_t slot=STARTING_NSLOT)
500 {
501 SpaceBox<dim,T> sbox(box);
502
503 // Initialize point transformation
504
505 Initialize(sbox,div,pad,slot);
506 }
507
508 /*! Initialize the cell list constructor
509 *
510 * \param box Domain where this cell list is living
511 * \param div grid size on each dimension
512 * \param pad padding cell
513 * \param slot maximum number of slot
514 *
515 */
516 void Initialize(const SpaceBox<dim,T> & box, const size_t (&div)[dim], const size_t pad = 1, size_t slot=STARTING_NSLOT)
517 {
518 Matrix<dim,T> mat;
519
520 CellDecomposer_sm<dim,T,transform>::setDimensions(box,div, mat, pad);
521 Mem_type::set_slot(slot);
522
523 // create the array that store the number of particle on each cell and se it to 0
524 InitializeStructures(this->gr_cell.getSize(),this->gr_cell.size());
525
526 from_cd = false;
527 }
528
529 //! Default Constructor
530 CellList()
531
532 :Mem_type(STARTING_NSLOT)
533 {};
534
535 //! Copy constructor
536 CellList(const CellList<dim,T,Mem_type,transform,vector_pos_type> & cell)
537 :Mem_type(STARTING_NSLOT)
538 {
539 this->operator=(cell);
540 }
541
542 //! Copy constructor
543 CellList(CellList<dim,T,Mem_type,transform,vector_pos_type> && cell)
544 :Mem_type(STARTING_NSLOT)
545 {
546 this->operator=(cell);
547 }
548
549
550 /*! \brief Cell list constructor
551 *
552 * \param box Domain where this cell list is living
553 * \param div grid size on each dimension
554 * \param mat Matrix transformation
555 * \param pad Cell padding
556 * \param slot maximum number of slot
557 *
558 */
559 CellList(Box<dim,T> & box, const size_t (&div)[dim], Matrix<dim,T> mat, const size_t pad = 1, size_t slot=STARTING_NSLOT)
560 :Mem_type(slot),CellDecomposer_sm<dim,T,transform>(box,div,mat,box.getP1(),pad)
561 {
562 SpaceBox<dim,T> sbox(box);
563 Initialize(sbox,div,pad,slot);
564 }
565
566 /*! \brief Cell list constructor
567 *
568 * \param box Domain where this cell list is living
569 * \param div grid size on each dimension
570 * \param pad Cell padding
571 * \param slot maximum number of slot
572 *
573 */
574 CellList(Box<dim,T> & box, const size_t (&div)[dim], const size_t pad = 1, size_t slot=STARTING_NSLOT)
575 :Mem_type(slot),n_dec(0)
576 {
577 SpaceBox<dim,T> sbox(box);
578 Initialize(sbox,div,pad,slot);
579 }
580
581 /*! \brief Cell list constructor
582 *
583 * \param box Domain where this cell list is living
584 * \param div grid size on each dimension
585 * \param pad Cell padding
586 * \param slot maximum number of slot
587 *
588 */
589 CellList(SpaceBox<dim,T> & box, const size_t (&div)[dim], const size_t pad = 1, size_t slot=STARTING_NSLOT)
590 :Mem_type(slot)
591 {
592 Initialize(box,div,pad,slot);
593 }
594
595 /*! \brief Cell list constructor from a cell decomposer
596 *
597 * \see Initialize
598 *
599 * \param cd_sm Cell-Decomposer
600 * \param box domain box (carefully this is going to be adjusted)
601 * \param pad Cell list padding
602 * \param slot number of slot for each cell
603 *
604 */
605 CellList(CellDecomposer_sm<dim,T,transform> & cd_sm, const Box<dim,T> & box, const size_t pad = 1, size_t slot=STARTING_NSLOT)
606 :Mem_type(slot),n_dec(0)
607 {
608 Initialize(cd_sm,box,pad,slot);
609 }
610
611 /*! \brief Destructor
612 *
613 *
614 */
615 ~CellList()
616 {}
617
618 /*! \brief Constructor from a temporal object
619 *
620 * \param cell Cell list structure
621 *
622 * \return itself
623 *
624 */
625 CellList<dim,T,Mem_type,transform> & operator=(CellList<dim,T,Mem_type,transform> && cell)
626 {
627 std::copy(&cell.NNc_full[0],&cell.NNc_full[openfpm::math::pow(3,dim)],&NNc_full[0]);
628 std::copy(&cell.NNc_sym[0],&cell.NNc_sym[openfpm::math::pow(3,dim)/2+1],&NNc_sym[0]);
629
630 Mem_type::swap(static_cast<Mem_type &&>(cell));
631
632 static_cast<CellDecomposer_sm<dim,T,transform> &>(*this).swap(cell);
633
634 n_dec = cell.n_dec;
635 from_cd = cell.from_cd;
636
637 return *this;
638 }
639
640 /*! \brief Constructor from a temporal object
641 *
642 * \param cell Cell list structure
643 *
644 * \return itself
645 *
646 */
647 CellList<dim,T,Mem_type,transform,vector_pos_type> & operator=(const CellList<dim,T,Mem_type,transform,vector_pos_type> & cell)
648 {
649 NNc_full = cell.NNc_full;
650 NNc_sym = cell.NNc_sym;
651
652 Mem_type::operator=(static_cast<const Mem_type &>(cell));
653
654 static_cast<CellDecomposer_sm<dim,T,transform> &>(*this) = static_cast<const CellDecomposer_sm<dim,T,transform> &>(cell);
655
656 n_dec = cell.n_dec;
657 from_cd = cell.from_cd;
658
659 return *this;
660 }
661
662 /*! \brief Constructor from a temporal object
663 *
664 * \param cell Cell list structure
665 *
666 * \return itself
667 *
668 */
669 template<typename Mem_type2>
670 CellList<dim,T,Mem_type,transform,vector_pos_type> & operator=(const CellList<dim,T,Mem_type2,transform,vector_pos_type> & cell)
671 {
672 NNc_full = cell.private_get_NNc_full();
673 NNc_sym = cell.private_get_NNc_sym();
674
675 Mem_type::copy_general(static_cast<const Mem_type2 &>(cell));
676
677 static_cast<CellDecomposer_sm<dim,T,transform> &>(*this) = static_cast<const CellDecomposer_sm<dim,T,transform> &>(cell);
678
679 n_dec = cell.get_ndec();
680 from_cd = cell.private_get_from_cd();
681
682 return *this;
683 }
684
685 /*! \brief Set the radius for the getNNIteratorRadius
686 *
687 * \param radius
688 *
689 */
690 void setRadius(T radius)
691 {
692 NNcalc_rad(radius,nnc_rad,this->getCellBox(),this->getGrid());
693 }
694
695 /*! \brief Get an iterator over particles following the cell structure
696 *
697 * \param dom_cells cells in the domain
698 *
699 * \return a particle iterator
700 *
701 */
702 ParticleIt_Cells<dim,CellList<dim,T,Mem_fast<>,transform>> getDomainIterator(openfpm::vector<size_t> & dom_cells)
703 {
704 ParticleIt_Cells<dim,CellList<dim,T,Mem_fast<>,transform>> it(*this,dom_cells);
705
706 return it;
707 }
708
709 /*! \brief Add to the cell
710 *
711 * \param cell_id Cell id where to add
712 * \param ele element to add
713 *
714 */
715 inline void addCell(size_t cell_id, typename Mem_type::local_index_type ele)
716 {
717 Mem_type::addCell(cell_id,ele);
718 }
719
720 /*! \brief Add an element in the cell list
721 *
722 * \param pos array that contain the coordinate
723 * \param ele element to store
724 *
725 */
726 inline void add(const T (& pos)[dim], typename Mem_type::local_index_type ele)
727 {
728 // calculate the Cell id
729 size_t cell_id = this->getCell(pos);
730
731 Mem_type::add(cell_id,ele);
732 }
733
734 /*! \brief Add an element in the cell list
735 *
736 * \param pos array that contain the coordinate
737 * \param ele element to store
738 *
739 */
740 inline void add(const Point<dim,T> & pos, typename Mem_type::local_index_type ele)
741 {
742 // calculate the Cell id
743 size_t cell_id = this->getCell(pos);
744
745 Mem_type::add(cell_id,ele);
746 }
747
748
749 /*! \brief Add an element in the cell list forcing to be in the domain cells
750 *
751 * \warning careful is intended to be used ONLY to avoid round-off problems
752 *
753 * \param pos array that contain the coordinate
754 * \param ele element to store
755 *
756 */
757 inline void addDom(const T (& pos)[dim], typename Mem_type::local_index_type ele)
758 {
759 // calculate the Cell id
760
761 size_t cell_id = this->getCellDom(pos);
762
763 // add the element to the cell
764
765 addCell(cell_id,ele);
766 }
767
768 /*! \brief Add an element in the cell list forcing to be in the domain cells
769 *
770 * \warning careful is intended to be used ONLY to avoid round-off problems
771 *
772 * \param pos array that contain the coordinate
773 * \param ele element to store
774 *
775 */
776 inline void addDom(const Point<dim,T> & pos, typename Mem_type::local_index_type ele)
777 {
778 // calculate the Cell id
779
780 size_t cell_id = this->getCellDom(pos);
781
782 // add the element to the cell
783
784 addCell(cell_id,ele);
785 }
786
787 /*! \brief Add an element in the cell list forcing to be in the padding cells
788 *
789 * \warning careful is intended to be used ONLY to avoid round-off problems
790 *
791 * \param pos array that contain the coordinate
792 * \param ele element to store
793 *
794 */
795 inline void addPad(const T (& pos)[dim], typename Mem_type::local_index_type ele)
796 {
797 // calculate the Cell id
798
799 size_t cell_id = this->getCellPad(pos);
800
801 // add the element to the cell
802
803 addCell(cell_id,ele);
804 }
805
806 /*! \brief Add an element in the cell list forcing to be in the padding cells
807 *
808 * \warning careful is intended to be used ONLY to avoid round-off problems
809 *
810 * \param pos array that contain the coordinate
811 * \param ele element to store
812 *
813 */
814 inline void addPad(const Point<dim,T> & pos, typename Mem_type::local_index_type ele)
815 {
816 // calculate the Cell id
817
818 size_t cell_id = this->getCell(pos);
819
820 // add the element to the cell
821
822 addCell(cell_id,ele);
823 }
824
825 /*! \brief remove an element from the cell
826 *
827 * \param cell cell id
828 * \param ele element id
829 *
830 */
831 inline void remove(size_t cell, size_t ele)
832 {
833 Mem_type::remove(cell,ele);
834 }
835
836 /*! \brief Get the number of cells this cell-list contain
837 *
838 * \return number of cells
839 */
840 inline size_t getNCells() const
841 {
842 return Mem_type::size();
843 }
844
845 /*! \brief Return the number of elements in the cell
846 *
847 * \param cell_id id of the cell
848 *
849 * \return number of elements in the cell
850 *
851 */
852 inline size_t getNelements(const size_t cell_id) const
853 {
854 return Mem_type::getNelements(cell_id);
855 }
856
857 /*! \brief Get an element in the cell
858 *
859 * \tparam i property to get
860 *
861 * \param cell cell id
862 * \param ele element id
863 *
864 * \return The element value
865 *
866 */
867 inline auto get(size_t cell, size_t ele) -> decltype(this->Mem_type::get(cell,ele))
868 {
869 return Mem_type::get(cell,ele);
870 }
871
872 /*! \brief Get an element in the cell
873 *
874 * \tparam i property to get
875 *
876 * \param cell cell id
877 * \param ele element id
878 *
879 * \return The element value
880 *
881 */
882 inline auto get(size_t cell, size_t ele) const -> decltype(this->Mem_type::get(cell,ele))
883 {
884 return Mem_type::get(cell,ele);
885 }
886
887 /*! \brief Swap the memory
888 *
889 * \param cl Cell list with witch you swap the memory
890 *
891 */
892 inline void swap(CellList<dim,T,Mem_type,transform,vector_pos_type> & cl)
893 {
894 NNc_full.swap(cl.NNc_full);
895 NNc_sym.swap(cl.NNc_sym);
896
897 Mem_type::swap(static_cast<Mem_type &>(cl));
898
899 static_cast<CellDecomposer_sm<dim,T,transform> &>(*this).swap(static_cast<CellDecomposer_sm<dim,T,transform> &>(cl));
900
901 n_dec = cl.n_dec;
902 from_cd = cl.from_cd;
903 }
904
905 /*! \brief Get the Cell iterator
906 *
907 * \param cell
908 *
909 * \return the iterator to the elements inside cell
910 *
911 */
912 CellIterator<CellList<dim,T,Mem_type,transform>> getCellIterator(size_t cell)
913 {
914 return CellIterator<CellList<dim,T,Mem_type,transform>>(cell,*this);
915 }
916
917 /*! \brief Get the Neighborhood iterator
918 *
919 * It iterate across all the element of the selected cell and the near cells
920 *
921 * \verbatim
922
923 * * *
924 * x *
925 * * *
926
927 \endverbatim
928 *
929 * * x is the selected cell
930 * * * are the near cell
931 *
932 * \param cell cell id
933 *
934 * \return An iterator across the neighhood particles
935 *
936 */
937 template<unsigned int impl=NO_CHECK> inline CellNNIteratorRadius<dim,CellList<dim,T,Mem_type,transform>,impl> getNNIteratorRadius(size_t cell)
938 {
939 CellNNIteratorRadius<dim,CellList<dim,T,Mem_type,transform>,impl> cln(cell,nnc_rad,*this);
940 return cln;
941 }
942
943 /*! \brief Get the Neighborhood iterator
944 *
945 * It iterate across all the element of the selected cell and the near cells
946 *
947 * \verbatim
948
949 * * *
950 * x *
951 * * *
952
953 \endverbatim
954 *
955 * * x is the selected cell
956 * * * are the near cell
957 *
958 * \param cell cell id
959 *
960 * \return An iterator across the neighhood particles
961 *
962 */
963 template<unsigned int impl=NO_CHECK>
964 __attribute__((always_inline)) inline CellNNIterator<dim,CellList<dim,T,Mem_type,transform,vector_pos_type>,(int)FULL,impl> getNNIterator(size_t cell)
965 {
966 CellNNIterator<dim,CellList<dim,T,Mem_type,transform,vector_pos_type>,(int)FULL,impl> cln(cell,NNc_full,*this);
967 return cln;
968
969 }
970
971 /*! \brief Get the symmetric Neighborhood iterator
972 *
973 * It iterate across all the element of the selected cell and the near cells up to some selected radius
974 *
975 * \param cell cell id
976 * \param r_cut radius
977 *
978 * \return An iterator across the neighborhood particles
979 *
980 */
981 template<unsigned int impl=NO_CHECK>
982 __attribute__((always_inline)) inline CellNNIteratorRadius<dim,CellList<dim,T,Mem_type,transform,vector_pos_type>,impl> getNNIteratorRadius(size_t cell, T r_cut)
983 {
984 openfpm::vector<long int> & NNc = rcache[r_cut];
985
986 if (NNc.size() == 0)
987 {NNcalc_rad(r_cut,NNc,this->getCellBox(),this->getGrid());}
988
989 CellNNIteratorRadius<dim,CellList<dim,T,Mem_type,transform,vector_pos_type>,impl> cln(cell,NNc,*this);
990
991 return cln;
992 }
993
994
995
996 /*! \brief Get the symmetric Neighborhood iterator
997 *
998 * It iterate across all the element of the selected cell and the near cells
999 *
1000 * \verbatim
1001
1002 * * *
1003 x *
1004
1005 \endverbatim
1006 *
1007 * * x is the selected cell
1008 * * * are the near cell
1009 *
1010 * \param cell cell id
1011 * \param p particle id
1012 *
1013 * \return An aiterator across the neighborhood particles
1014 *
1015 */
1016 template<unsigned int impl>
1017 __attribute__((always_inline)) inline CellNNIteratorSym<dim,CellList<dim,T,Mem_type,transform,vector_pos_type>,vector_pos_type,(unsigned int)SYM,impl>
1018 getNNIteratorSym(size_t cell, size_t p, const vector_pos_type & v)
1019 {
1020#ifdef SE_CLASS1
1021 if (from_cd == false)
1022 {std::cerr << __FILE__ << ":" << __LINE__ << " Warning when you try to get a symmetric neighborhood iterator, you must construct the Cell-list in a symmetric way" << std::endl;}
1023#endif
1024
1025 CellNNIteratorSym<dim,CellList<dim,T,Mem_type,transform,vector_pos_type>,vector_pos_type,SYM,impl> cln(cell,p,NNc_sym,*this,v);
1026 return cln;
1027 }
1028
1029 /*! \brief Get the symmetric Neighborhood iterator
1030 *
1031 * It iterate across all the element of the selected cell and the near cells
1032 *
1033 * \verbatim
1034
1035 * * *
1036 x *
1037
1038 \endverbatim
1039 *
1040 * * x is the selected cell
1041 * * * are the near cell
1042 *
1043 * \param cell cell id
1044 * \param p particle id
1045 * \param v_p1 first phase for particle p
1046 * \param v_p2 second phase for particle q
1047 *
1048 * \return An aiterator across the neighborhood particles
1049 *
1050 */
1051 template<unsigned int impl, typename vector_pos_type2>
1052 __attribute__((always_inline)) inline CellNNIteratorSymMP<dim,CellList<dim,T,Mem_type,transform,vector_pos_type>,vector_pos_type2,(unsigned int)SYM,impl>
1053 getNNIteratorSymMP(size_t cell, size_t p, const vector_pos_type2 & v_p1, const vector_pos_type2 & v_p2)
1054 {
1055#ifdef SE_CLASS1
1056 if (from_cd == false)
1057 std::cerr << __FILE__ << ":" << __LINE__ << " Warning when you try to get a symmetric neighborhood iterator, you must construct the Cell-list in a symmetric way" << std::endl;
1058#endif
1059
1060 CellNNIteratorSymMP<dim,CellList<dim,T,Mem_type,transform,vector_pos_type>,vector_pos_type2,SYM,impl> cln(cell,p,NNc_sym,*this,v_p1,v_p2);
1061 return cln;
1062 }
1063
1064 /*! \brief Get the symmetric neighborhood
1065 *
1066 * \return the symmetric neighborhood
1067 *
1068 */
1069 const NNc_array<dim,(unsigned int)openfpm::math::pow(3,dim)/2+1> & getNNc_sym() const
1070 {
1071 return NNc_sym;
1072 }
1073
1074 /*! \brief Return the number of padding cells of the Cell decomposer
1075 *
1076 * \param i dimension
1077 *
1078 * \return the number of padding cells
1079 *
1080 */
1081 size_t getPadding(size_t i) const
1082 {
1083 return CellDecomposer_sm<dim,T,transform>::getPadding(i);
1084 }
1085
1086 /*! \brief Return the number of padding cells of the Cell decomposer as an array
1087 *
1088 *
1089 * \return the number of padding cells
1090 *
1091 */
1092 size_t (& getPadding())[dim]
1093 {
1094 return CellDecomposer_sm<dim,T,transform>::getPadding();
1095 }
1096
1097 /*! \brief Clear the cell list
1098 *
1099 */
1100 void clear()
1101 {
1102 Mem_type::clear();
1103 }
1104
1105 /*! \brief Litterary destroy the memory of the cell list, including the retained one
1106 *
1107 */
1108 void destroy()
1109 {
1110 Mem_type::destroy();
1111 }
1112
1113 /*! \brief Return the starting point of the cell p
1114 *
1115 * \param cell_id cell id
1116 *
1117 * \return the index
1118 *
1119 */
1120 __attribute__((always_inline)) inline const typename Mem_type::local_index_type & getStartId(typename Mem_type::local_index_type cell_id) const
1121 {
1122 return Mem_type::getStartId(cell_id);
1123 }
1124
1125 /*! \brief Return the end point of the cell p
1126 *
1127 * \param cell_id cell id
1128 *
1129 * \return the stop index
1130 *
1131 */
1132 __attribute__((always_inline)) inline const typename Mem_type::local_index_type & getStopId(typename Mem_type::local_index_type cell_id) const
1133 {
1134 return Mem_type::getStopId(cell_id);
1135 }
1136
1137 /*! \brief Return the neighborhood id
1138 *
1139 * \param part_id particle id
1140 *
1141 * \return the neighborhood id
1142 *
1143 */
1144 __attribute__((always_inline)) inline const typename Mem_type::local_index_type & get_lin(const typename Mem_type::local_index_type * part_id) const
1145 {
1146 return Mem_type::get_lin(part_id);
1147 }
1148
1149//////////////////////////////// POINTLESS BUT REQUIRED TO RESPECT THE INTERFACE //////////////////
1150
1151 //! Ghost marker
1152 size_t g_m = 0;
1153
1154 /*! \brief return the ghost marker
1155 *
1156 * \return ghost marker
1157 *
1158 */
1159 inline size_t get_gm()
1160 {
1161 return g_m;
1162 }
1163
1164 /*! \brief Set the ghost marker
1165 *
1166 * \param g_m marker
1167 *
1168 */
1169 inline void set_gm(size_t g_m)
1170 {
1171 this->g_m = g_m;
1172 }
1173
1174#ifdef CUDA_GPU
1175
1176 CellList_cpu_ker<dim,T,typename Mem_type::toKernel_type,transform> toKernel()
1177 {
1178 typedef typename Mem_type::local_index_type ids_type;
1179 typedef typename Mem_type::local_index_type cnt_type;
1180
1181 openfpm::array<T,dim,cnt_type> spacing_c;
1182 openfpm::array<ids_type,dim,cnt_type> div_c;
1183 openfpm::array<ids_type,dim,cnt_type> off;
1184
1185 for (size_t i = 0 ; i < dim ; i++)
1186 {
1187 spacing_c[i] = CellDecomposer_sm<dim,T,transform>::getCellBox().getHigh(i);
1188 div_c[i] = CellDecomposer_sm<dim,T,transform>::getDiv()[i];
1189 off[i] = CellDecomposer_sm<dim,T,transform>::getPadding(i);
1190 }
1191
1192 CellList_cpu_ker<dim,T,typename Mem_type::toKernel_type,transform> cl(Mem_type::toKernel(),
1193 spacing_c,
1194 div_c,
1195 off,
1196 CellDecomposer_sm<dim,T,transform>::getTransform());
1197
1198 return cl;
1199 }
1200
1201
1202 void hostToDevice()
1203 {
1204 Mem_type::hostToDevice();
1205 }
1206
1207#endif
1208
1209 const NNc_array<dim,(unsigned int)openfpm::math::pow(3,dim)> & private_get_NNc_full () const
1210 {
1211 return NNc_full;
1212 }
1213
1214 const NNc_array<dim,(unsigned int)openfpm::math::pow(3,dim)/2+1> & private_get_NNc_sym () const
1215 {
1216 return NNc_sym;
1217 }
1218
1219 bool private_get_from_cd() const
1220 {
1221 return from_cd;
1222 }
1223
1224/////////////////////////////////////
1225
1226 void re_setBoxNN()
1227 {}
1228
1229/////////////////////////////////////
1230
1231 /*! \brief Set the n_dec number
1232 *
1233 * \param n_dec
1234 *
1235 */
1236 void set_ndec(size_t n_dec)
1237 {
1238 this->n_dec = n_dec;
1239 }
1240
1241 /*! \brief Set the n_dec number
1242 *
1243 * \return n_dec
1244 *
1245 */
1246 size_t get_ndec() const
1247 {
1248 return n_dec;
1249 }
1250
1251/////////////////////////////////////
1252};
1253
1254/*! \brief Calculate parameters for the cell list
1255 *
1256 * \param div Division array
1257 * \param r_cut interation radius or size of each cell
1258 * \param enlarge In case of padding particles the cell list must be enlarged, like a ghost. This parameter says how much must be enlarged
1259 *
1260 * \return the processor bounding box
1261 */
1262template<unsigned int dim, typename St> static inline void cl_param_calculate(Box<dim, St> & pbox, size_t (&div)[dim], St r_cut, const Ghost<dim, St> & enlarge)
1263{
1264 // calculate the parameters of the cell list
1265
1266 // extend by the ghost
1267 pbox.enlarge(enlarge);
1268
1269 // Calculate the division array and the cell box
1270 for (size_t i = 0; i < dim; i++)
1271 {
1272 div[i] = static_cast<size_t>((pbox.getP2().get(i) - pbox.getP1().get(i)) / r_cut);
1273 div[i]++;
1274 pbox.setHigh(i,pbox.getLow(i) + div[i]*r_cut);
1275 }
1276}
1277
1278/*! \brief Calculate parameters for the symmetric cell list
1279 *
1280 * \param[in] dom Simulation domain
1281 * \param[output] cd_sm This cell-decomposer is set according to the needed division
1282 * \param[in] g Ghost dimensions
1283 * \param[output] pad required padding for the cell-list
1284 *
1285 * \return the processor bounding box
1286 */
1287template<unsigned int dim, typename St> static
1288inline void cl_param_calculateSym(const Box<dim,St> & dom,
1289 CellDecomposer_sm<dim,St,shift<dim,St>> & cd_sm,
1290 Ghost<dim,St> g,
1291 St r_cut,
1292 size_t & pad)
1293{
1294 size_t div[dim];
1295
1296 for (size_t i = 0 ; i < dim ; i++)
1297 div[i] = (dom.getHigh(i) - dom.getLow(i)) / r_cut;
1298
1299 g.magnify(1.013);
1300
1301 // Calculate the maximum padding
1302 for (size_t i = 0 ; i < dim ; i++)
1303 {
1304 size_t tmp = std::ceil(fabs(g.getLow(i)) / r_cut);
1305 pad = (pad > tmp)?pad:tmp;
1306
1307 tmp = std::ceil(fabs(g.getHigh(i)) / r_cut);
1308 pad = (pad > tmp)?pad:tmp;
1309 }
1310
1311 cd_sm.setDimensions(dom,div,pad);
1312}
1313
1314
1315/*! \brief Calculate parameters for the symmetric cell list
1316 *
1317 * \param[in] dom Simulation domain
1318 * \param[out] cd_sm This cell-decomposer is set according to the needed division
1319 * \param[in] g Ghost part extension
1320 * \param[out] pad required padding for the cell-list
1321 *
1322 * \return the processor bounding box
1323 */
1324template<unsigned int dim, typename St>
1325static inline void cl_param_calculateSym(const Box<dim,St> & dom,
1326 CellDecomposer_sm<dim,St,shift<dim,St>> & cd_sm,
1327 Ghost<dim,St> g,
1328 size_t & pad)
1329{
1330 size_t div[dim];
1331
1332 for (size_t i = 0 ; i < dim ; i++)
1333 div[i] = (dom.getHigh(i) - dom.getLow(i)) / g.getHigh(i);
1334
1335 g.magnify(1.013);
1336
1337 pad = 1;
1338
1339 cd_sm.setDimensions(dom,div,pad);
1340}
1341
1342
1343#endif /* CELLLIST_HPP_ */
1344