CellList.hpp 27 KB
Newer Older
incardon's avatar
incardon committed
1
2
3
4
5
6
7
8
9
10
/*
 * CellList.hpp
 *
 *  Created on: Mar 21, 2015
 *      Author: Pietro Incardona
 */

#ifndef CELLLIST_HPP_
#define CELLLIST_HPP_

11
#include "Vector/map_vector.hpp"
12
#include "CellDecomposer.hpp"
incardon's avatar
incardon committed
13
14
#include "Space/SpaceBox.hpp"
#include "util/mathutil.hpp"
15
#include "CellNNIterator.hpp"
incardon's avatar
incardon committed
16
#include "Space/Shape/HyperCube.hpp"
17
18
#include "CellListNNIteratorRadius.hpp"
#include <unordered_map>
19

incardon's avatar
incardon committed
20
21
22
23
24
#include "CellListIterator.hpp"
#include "ParticleIt_Cells.hpp"
#include "ParticleItCRS_Cells.hpp"
#include "util/common.hpp"

25
26
27
#include "NN/Mem_type/MemFast.hpp"
#include "NN/Mem_type/MemBalanced.hpp"
#include "NN/Mem_type/MemMemoryWise.hpp"
incardon's avatar
incardon committed
28
#include "NN/CellList/NNc_array.hpp"
incardon's avatar
incardon committed
29

30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
//! Wrapper of the unordered map
template<typename key,typename val>
class wrap_unordered_map: public std::unordered_map<key,val>
{
};

#ifdef HAVE_LIBQUADMATH

#include <boost/multiprecision/float128.hpp>


//! Wrapper of the unordered map
template<typename val>
class wrap_unordered_map<boost::multiprecision::float128,val>
{
};

#endif
48

49
//! Point at witch the cell do a reallocation (it should the the maximum for all the implementations)
50
#define CELL_REALLOC 16ul
51

52
53
#define STARTING_NSLOT 16

incardon's avatar
incardon committed
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
/*! \brief Calculate the the Neighborhood for symmetric interactions CSR scheme
 *
 * \param cNN calculated cross neighborhood
 * \param div Number of divisions in each direction
 *
 */
template<unsigned int dim> void NNcalc_csr(openfpm::vector<std::pair<grid_key_dx<dim>,grid_key_dx<dim>>> & cNN)
{
	// Calculate the NNc_full array, it is a structure to get the neighborhood array

	// compile-time array {0,0,0,....}  {2,2,2,...} {1,1,1,...}

	typedef typename generate_array<size_t,dim, Fill_zero>::result NNzero;
	typedef typename generate_array<size_t,dim, Fill_two>::result NNtwo;

	// Generate the sub-grid iterator

	size_t div[dim];

	// Calculate the divisions

	for (size_t i = 0 ; i < dim ; i++)
		div[i] = 4;

	grid_sm<dim,void> gs(div);
	grid_key_dx_iterator_sub<dim> gr_sub3(gs,NNzero::data,NNtwo::data);

	grid_key_dx<dim> src_; // Source cell
	for (size_t i = 0; i < dim; i++)
		src_.set_d(i,1);

	size_t middle = gs.LinId(src_);

	// Calculate the symmetric crs array
	while (gr_sub3.isNext())
	{
		auto dst = gr_sub3.get();
		grid_key_dx<dim> src = src_;

		if ((long int)middle > gs.LinId(dst))
		{
			++gr_sub3;
			continue;
		}

		// Here we adjust src and dst to be in the positive quadrant

		for (size_t i = 0 ; i < dim; i++)
		{
			if (dst.get(i) == 0)
			{
				src.set_d(i,src.get(i) + 1);
				dst.set_d(i,dst.get(i) + 1);
			}
		}

		src -= src_;
		dst -= src_;

		cNN.add(std::pair<grid_key_dx<dim>,grid_key_dx<dim>>(src,dst));

		++gr_sub3;

	}
};

/*! \brief Calculate the the Neighborhood for symmetric interactions
 *
 * \param cNN calculated cross neighborhood
 * \param div Number of divisions in each direction
 *
 */
template<unsigned int dim> void NNcalc_sym(openfpm::vector<grid_key_dx<dim>> & cNN)
{
	// Calculate the NNc_full array, it is a structure to get the neighborhood array

	// compile-time array {0,0,0,....}  {2,2,2,...} {1,1,1,...}

	typedef typename generate_array<size_t,dim, Fill_zero>::result NNzero;
	typedef typename generate_array<size_t,dim, Fill_two>::result NNtwo;

	// Generate the sub-grid iterator

	size_t div[dim];

	// Calculate the divisions

	for (size_t i = 0 ; i < dim ; i++)
		div[i] = 4;

	grid_sm<dim,void> gs(div);
	grid_key_dx_iterator_sub<dim> gr_sub3(gs,NNzero::data,NNtwo::data);

	grid_key_dx<dim> src_; // Source cell
	for (size_t i = 0; i < dim; i++)
		src_.set_d(i,1);

	size_t middle = gs.LinId(src_);

	// Calculate the symmetric array
	while (gr_sub3.isNext())
	{
		auto dst = gr_sub3.get();

		if ((long int)middle > gs.LinId(dst))
		{
			++gr_sub3;
			continue;
		}

		cNN.add(dst - src_);

		++gr_sub3;

	}
};

/*! \brief Calculate the Neighborhood cells
 *
 * \param cNN calculated cross neighborhood
 * \param div Number of divisions in each direction
 *
 */
template<unsigned int dim> void NNcalc_full(openfpm::vector<grid_key_dx<dim>> & cNN)
{
	// Calculate the NNc_full array, it is a structure to get the neighborhood array

	// compile-time array {0,0,0,....}  {2,2,2,...} {1,1,1,...}

	typedef typename generate_array<size_t,dim, Fill_zero>::result NNzero;
	typedef typename generate_array<size_t,dim, Fill_two>::result NNtwo;

	// Generate the sub-grid iterator

	size_t div[dim];

	// Calculate the divisions

	for (size_t i = 0 ; i < dim ; i++)
incardon's avatar
incardon committed
193
	{div[i] = 4;}
incardon's avatar
incardon committed
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213

	grid_sm<dim,void> gs(div);
	grid_key_dx_iterator_sub<dim> gr_sub3(gs,NNzero::data,NNtwo::data);

	grid_key_dx<dim> src_; // Source cell
	for (size_t i = 0; i < dim; i++)
		src_.set_d(i,1);

	// Calculate the symmetric crs array
	while (gr_sub3.isNext())
	{
		auto dst = gr_sub3.get();

		cNN.add(dst - src_);

		++gr_sub3;
	}
};


incardon's avatar
incardon committed
214
215
216
217
218
219
220
/* NOTE all the implementations
 *
 * has complexity O(1) in getting the cell id and the elements in a cell
 * but with different constants
 *
 */

221
222
223
224
/*! \brief Class for FAST cell list implementation
 *
 * This class implement the FAST cell list, fast but memory
 * expensive. The memory allocation is (M * N_cell_max)*sizeof(ele) + M*8
incardon's avatar
incardon committed
225
 *
226
227
228
 * * M = number of cells
 * * N_cell_max = maximum number of elements in a cell
 * * ele = element the structure is storing
incardon's avatar
incardon committed
229
 *
230
 * \note Because N_cell_max >= N/M then M * N_cell_max >= O(N)
incardon's avatar
incardon committed
231
 *
232
 * \warning Not not use for high asymmetric distribution
incardon's avatar
incardon committed
233
 *
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
 * Example of a 2D Cell list 6x6 structure with padding 1 without shift, cell indicated with p are padding cell
 * the origin of the cell or point (0,0) is marked with cell number 9
 *
 * \verbatim
 * +-----------------------+
 * |p |p |p |p |p |p |p |p |
 * +-----------------------+
 * |p |  |  |  |  |  |  |p |
 * +-----------------------+
 * |p |  |  |  |  |  |  |p |
 * +-----------------------+
 * |p |  |  |  |  |  |  |p |
 * +-----------------------+
 * |p |9 |  |  |  |  |  |p |
 * +-----------------------+
 * |p |p |p |p |p |p |p |p |
 * +-----------------------+
 * \endverbatim
 *
 *
 * \tparam dim Dimensionality of the space
 * \tparam T type of the space float, double ...
 * \tparam base Base structure that store the information
 *
 * ### Declaration of a cell list
 * \snippet CellList_test.hpp Declare a cell list
 * ### Usage of cell list [CellS == CellList<3,double,FAST>]
 * \snippet CellList_test.hpp Usage of cell list
 * ### Remove one particle from each cell
 * \snippet CellList_test.hpp remove one particle from each cell
 * ### Usage of the neighborhood iterator
 * \snippet CellList_test.hpp Usage of the neighborhood iterator
incardon's avatar
incardon committed
266
267
 *
 */
268
template<unsigned int dim, typename T,  typename Mem_type, typename transform = no_transform<dim,T>, typename base=openfpm::vector<size_t>>
incardon's avatar
incardon committed
269
class CellList : public CellDecomposer_sm<dim,T,transform>, public Mem_type
incardon's avatar
incardon committed
270
{
271
272
273
274
275
276
277
protected:
	//! The array contain the neighborhood of the cell-id in case of asymmetric interaction
	//
	//    * * *
	//    * x *
	//    * * *

incardon's avatar
incardon committed
278
279
280

	NNc_array<dim,openfpm::math::pow(3,dim)> NNc_full;
//	long int NNc_full[openfpm::math::pow(3,dim)];
281
282
283
284
285
286

	//! The array contain the neighborhood of the cell-id in case of symmetric interaction
	//
	//   * * *
	//     x *
	//
incardon's avatar
incardon committed
287
288
//	long int NNc_sym[openfpm::math::pow(3,dim)/2+1];
	NNc_array<dim,openfpm::math::pow(3,dim)/2+1> NNc_sym;
289
290
291
292

private:

	//! Caching of r_cutoff radius
293
	wrap_unordered_map<T,openfpm::vector<long int>> rcache;
294

incardon's avatar
incardon committed
295
296
297
298
299
300
	//! True if has been initialized from CellDecomposer
	bool from_cd;

	//! Additional information in general (used to understand if the cell-list)
	//! has been constructed from an old decomposition
	size_t n_dec;
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369

	/*! Calculate the neighborhood cells based on the radius
	 *
	 * \note To the calculated neighborhood cell you have to add the id of the central cell
	 *
		\verbatim
       +-----------------------+
       |p |p |p |p |p |p |p |p |
       +-----------------------+
       |p |  |  |  |  |  |  |p |
       +-----------------------+
       |p |  |  |7 |8 |9 |  |p |
       +-----------------------+
       |p |  |  |-1|0 |1 |  |p |
       +-----------------------+
       |p |9 |  |-9|-8|-7|  |p |
       +-----------------------+
       |p |p |p |p |p |p |p |p |
       +-----------------------+
		\endverbatim
	 *
	 * The number indicate the cell id calculated
	 *
	 * -9,-8,-7,-1,0,1,7,8,9
	 *
	 * The cell 0 has id = 22 in the big cell matrix, so to calculate the
	 * neighborhood cells you have to sum the id of the center cell
	 *
	 * 13,14,15,21,22,23,29,30,31
	 *
	 * \param r_cut Cutoff-radius
	 * \param NNcell vector containing the neighborhood cells ids
	 *
	 */
	void NNcalc(T r_cut, openfpm::vector<long int> & NNcell)
	{
		size_t n_cell[dim];
		size_t n_cell_mid[dim];

		Point<dim,T> spacing = this->getCellBox().getP2();
		const grid_sm<dim,void> & gs = this->getGrid();

		for (size_t i = 0 ; i < dim ; i++)
		{
			n_cell[i] = 2*(std::ceil(r_cut / spacing.get(i)))+1;
			n_cell_mid[i] = n_cell[i] / 2;
		}

		grid_sm<dim,void> gsc(n_cell);
		grid_key_dx_iterator<dim> gkdi(gsc);

		while (gkdi.isNext())
		{
			auto key = gkdi.get();

			for (size_t i = 0 ; i < dim ; i++)
				key.set_d(i,key.get(i) - n_cell_mid[i]);

			NNcell.add(gs.LinId(key));

			++gkdi;
		}
	}

	//! Initialize the structures of the data structure
	void InitializeStructures(const size_t (& div)[dim], size_t tot_n_cell, size_t slot=STARTING_NSLOT)
	{
		Mem_type::init_to_zero(slot,tot_n_cell);

incardon's avatar
incardon committed
370
371
		NNc_full.set_size(div);
		NNc_full.init_full();
372

incardon's avatar
incardon committed
373
374
		NNc_sym.set_size(div);
		NNc_sym.init_sym();
incardon's avatar
incardon committed
375
	}
376

incardon's avatar
incardon committed
377
378
379
	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
	{
		size_t bc[dim];
380

incardon's avatar
incardon committed
381
382
		for (size_t i = 0 ; i < dim ; i++)
			bc[i] = NON_PERIODIC;
383

incardon's avatar
incardon committed
384
		Box<dim,long int> bx = cd_sm.convertDomainSpaceIntoCellUnits(dom_box,bc);
385

incardon's avatar
incardon committed
386
387
		size_t div[dim];
		size_t div_big[dim];
388

incardon's avatar
incardon committed
389
390
391
392
		for (size_t i = 0 ; i < dim ; i++)
		{
			div[i] = bx.getHigh(i) - bx.getLow(i);
			div_big[i] = cd_sm.getDiv()[i] - 2*cd_sm.getPadding(i);
393
		}
incardon's avatar
incardon committed
394
395

		cd.setDimensions(cd_sm.getDomain(),div_big,div, pad, bx.getP1());
396
397
398
399
	}

public:

incardon's avatar
incardon committed
400
401
402
	//! Type of internal memory structure
	typedef Mem_type Mem_type_type;

incardon's avatar
incardon committed
403
404
	typedef CellNNIteratorSym<dim,CellList<dim,T,Mem_type,transform,base>,RUNTIME,NO_CHECK> SymNNIterator;

405
406
407
	//! Object type that the structure store
	typedef typename base::value_type value_type;

incardon's avatar
incardon committed
408
	//! Type of the coordinate space (double float)
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
	typedef T stype;

	/*! \brief Return the underlying grid information of the cell list
	 *
	 * \return the grid infos
	 *
	 */
	const grid_sm<dim,void> & getGrid()
	{
		return CellDecomposer_sm<dim,T,transform>::getGrid();
	}

	/*! Initialize the cell list from a well-define Cell-decomposer
	 *
	 * In some cases is needed to have a Cell-list with Cells consistent
	 * with a well predefined CellDecomposer. In this case we use this function.
	 * Using this initialization the Cell-list maintain the Cells defined by this
	 * Cell-decomposer consistently
	 *
	 * \param cd_sm Cell-Decomposer
	 * \param dom_box domain box (carefully this is going to be adjusted)
	 * \param pad cell-list padding
	 * \param slot slots for each cell
	 *
	 */
	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)
	{
		size_t bc[dim];
		for (size_t i = 0 ; i < dim ; i++)	{bc[i] = NON_PERIODIC;}

		Box<dim,long int> bx = cd_sm.convertDomainSpaceIntoCellUnits(dom_box,bc);

incardon's avatar
incardon committed
441
442
		setCellDecomposer(*this,cd_sm,dom_box,pad);

443
444
445
446
447
448
449
450
451
		size_t div_w_pad[dim];
		size_t tot_cell = 1;

		for (size_t i = 0 ; i < dim ; i++)
		{
			div_w_pad[i] = bx.getHigh(i) - bx.getLow(i) + 2*pad;
			tot_cell *= div_w_pad[i];
		}

incardon's avatar
incardon committed
452
		// here we set the cell-shift for the CellDecomposer
453
454

		InitializeStructures(div_w_pad,tot_cell);
incardon's avatar
incardon committed
455
456
457

		// Initialized from CellDecomposer
		from_cd = true;
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
	}

	/*! Initialize the cell list
	 *
	 * \param box Domain where this cell list is living
	 * \param div grid size on each dimension
	 * \param pad padding cell
	 * \param slot maximum number of slot
	 *
	 */
	void Initialize(const Box<dim,T> & box, const size_t (&div)[dim], const size_t pad = 1, size_t slot=STARTING_NSLOT)
	{
		SpaceBox<dim,T> sbox(box);

		// Initialize point transformation

		Initialize(sbox,div,pad,slot);
	}

	/*! Initialize the cell list constructor
	 *
	 * \param box Domain where this cell list is living
	 * \param div grid size on each dimension
	 * \param pad padding cell
	 * \param slot maximum number of slot
	 *
	 */
	void Initialize(const SpaceBox<dim,T> & box, const size_t (&div)[dim], const size_t pad = 1, size_t slot=STARTING_NSLOT)
	{
		Matrix<dim,T> mat;

		CellDecomposer_sm<dim,T,transform>::setDimensions(box,div, mat, pad);
incardon's avatar
incardon committed
490
		Mem_type::set_slot(slot);
491
492
493
494

		// create the array that store the number of particle on each cell and se it to 0
		InitializeStructures(this->gr_cell.getSize(),this->gr_cell.size());

incardon's avatar
incardon committed
495
		from_cd = false;
496
497
498
499
	}

	//! Default Constructor
	CellList()
Yaroslav's avatar
Yaroslav committed
500

501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
	:Mem_type(STARTING_NSLOT)
	{};

	//! Copy constructor
	CellList(const CellList<dim,T,Mem_type,transform,base> & cell)
	:Mem_type(STARTING_NSLOT)
	{
		this->operator=(cell);
	}

	//! Copy constructor
	CellList(CellList<dim,T,Mem_type,transform,base> && cell)
	:Mem_type(STARTING_NSLOT)
	{
		this->operator=(cell);
	}


	/*! \brief Cell list constructor
	 *
	 * \param box Domain where this cell list is living
	 * \param div grid size on each dimension
	 * \param mat Matrix transformation
	 * \param pad Cell padding
	 * \param slot maximum number of slot
	 *
	 */
	CellList(Box<dim,T> & box, const size_t (&div)[dim], Matrix<dim,T> mat, const size_t pad = 1, size_t slot=STARTING_NSLOT)
	:Mem_type(slot),CellDecomposer_sm<dim,T,transform>(box,div,mat,box.getP1(),pad)
	{
		SpaceBox<dim,T> sbox(box);
		Initialize(sbox,div,pad,slot);
	}

	/*! \brief Cell list constructor
	 *
	 * \param box Domain where this cell list is living
	 * \param div grid size on each dimension
	 * \param pad Cell padding
	 * \param slot maximum number of slot
	 *
	 */
	CellList(Box<dim,T> & box, const size_t (&div)[dim], const size_t pad = 1, size_t slot=STARTING_NSLOT)
incardon's avatar
incardon committed
544
	:Mem_type(slot),n_dec(0)
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
	{
		SpaceBox<dim,T> sbox(box);
		Initialize(sbox,div,pad,slot);
	}

	/*! \brief Cell list constructor
	 *
	 * \param box Domain where this cell list is living
	 * \param div grid size on each dimension
	 * \param pad Cell padding
	 * \param slot maximum number of slot
	 *
	 */
	CellList(SpaceBox<dim,T> & box, const size_t (&div)[dim], const size_t pad = 1, size_t slot=STARTING_NSLOT)
	:Mem_type(slot)
	{
		Initialize(box,div,pad,slot);
	}

	/*! \brief Cell list constructor from a cell decomposer
	 *
	 * \see Initialize
	 *
	 * \param cd_sm Cell-Decomposer
	 * \param box domain box (carefully this is going to be adjusted)
	 * \param pad Cell list padding
	 * \param slot number of slot for each cell
	 *
	 */
	CellList(CellDecomposer_sm<dim,T,transform> & cd_sm, const Box<dim,T> & box, const size_t pad = 1, size_t slot=STARTING_NSLOT)
incardon's avatar
incardon committed
575
	:Mem_type(slot),n_dec(0)
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
	{
		Initialize(cd_sm,box,pad,slot);
	}

	/*! \brief Destructor
	 *
	 *
	 */
	~CellList()
	{}

	/*! \brief Constructor from a temporal object
	 *
	 * \param cell Cell list structure
	 *
	 * \return itself
	 *
	 */
	CellList<dim,T,Mem_type,transform,base> & operator=(CellList<dim,T,Mem_type,transform,base> && cell)
	{
		std::copy(&cell.NNc_full[0],&cell.NNc_full[openfpm::math::pow(3,dim)],&NNc_full[0]);
		std::copy(&cell.NNc_sym[0],&cell.NNc_sym[openfpm::math::pow(3,dim)/2+1],&NNc_sym[0]);

599
		Mem_type::swap(static_cast<Mem_type &&>(cell));
600
601
602

		static_cast<CellDecomposer_sm<dim,T,transform> &>(*this).swap(cell);

incardon's avatar
incardon committed
603
604
		n_dec = cell.n_dec;

605
606
607
608
609
610
611
612
613
614
615
616
		return *this;
	}

	/*! \brief Constructor from a temporal object
	 *
	 * \param cell Cell list structure
	 *
	 * \return itself
	 *
	 */
	CellList<dim,T,Mem_type,transform,base> & operator=(const CellList<dim,T,Mem_type,transform,base> & cell)
	{
incardon's avatar
incardon committed
617
618
		NNc_full = cell.NNc_full;
		NNc_sym = cell.NNc_sym;
619

620
		Mem_type::operator=(static_cast<const Mem_type &>(cell));
621
622
623

		static_cast<CellDecomposer_sm<dim,T,transform> &>(*this) = static_cast<const CellDecomposer_sm<dim,T,transform> &>(cell);

incardon's avatar
incardon committed
624
625
		n_dec = cell.n_dec;

626
627
628
		return *this;
	}

incardon's avatar
incardon committed
629
	/*! \brief Get an iterator over particles following the cell structure
630
631
	 *
	 * \param dom_cells cells in the domain
incardon's avatar
incardon committed
632
633
634
635
	 *
	 * \return a particle iterator
	 *
	 */
incardon's avatar
incardon committed
636
	ParticleIt_Cells<dim,CellList<dim,T,Mem_fast<>,transform,base>> getDomainIterator(openfpm::vector<size_t> & dom_cells)
incardon's avatar
incardon committed
637
	{
incardon's avatar
incardon committed
638
		ParticleIt_Cells<dim,CellList<dim,T,Mem_fast<>,transform,base>> it(*this,dom_cells);
incardon's avatar
incardon committed
639
640
641
642

		return it;
	}

643
644
645
646
647
648
649
650
	/*! \brief Add to the cell
	 *
	 * \param cell_id Cell id where to add
	 * \param ele element to add
	 *
	 */
	inline void addCell(size_t cell_id, typename base::value_type ele)
	{
651
		Mem_type::addCell(cell_id,ele);
652
653
654
655
656
657
658
659
660
661
	}

	/*! \brief Add an element in the cell list
	 *
	 * \param pos array that contain the coordinate
	 * \param ele element to store
	 *
	 */
	inline void add(const T (& pos)[dim], typename base::value_type ele)
	{
incardon's avatar
incardon committed
662
663
664
665
		// calculate the Cell id
		size_t cell_id = this->getCell(pos);

		Mem_type::add(cell_id,ele);
666
667
668
669
670
671
672
673
674
675
	}

	/*! \brief Add an element in the cell list
	 *
	 * \param pos array that contain the coordinate
	 * \param ele element to store
	 *
	 */
	inline void add(const Point<dim,T> & pos, typename base::value_type ele)
	{
incardon's avatar
incardon committed
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
		// calculate the Cell id
		size_t cell_id = this->getCell(pos);

		Mem_type::add(cell_id,ele);
	}


	/*! \brief Add an element in the cell list forcing to be in the domain cells
	 *
	 * \warning careful is intended to be used ONLY to avoid round-off problems
	 *
	 * \param pos array that contain the coordinate
	 * \param ele element to store
	 *
	 */
	inline void addDom(const T (& pos)[dim], typename base::value_type ele)
	{
		// calculate the Cell id

		size_t cell_id = this->getCellDom(pos);

		// add the element to the cell

		addCell(cell_id,ele);
	}

	/*! \brief Add an element in the cell list forcing to be in the domain cells
	 *
	 * \warning careful is intended to be used ONLY to avoid round-off problems
	 *
	 * \param pos array that contain the coordinate
	 * \param ele element to store
	 *
	 */
	inline void addDom(const Point<dim,T> & pos, typename base::value_type ele)
	{
		// calculate the Cell id

		size_t cell_id = this->getCellDom(pos);

		// add the element to the cell

		addCell(cell_id,ele);
	}

	/*! \brief Add an element in the cell list forcing to be in the padding cells
	 *
	 * \warning careful is intended to be used ONLY to avoid round-off problems
	 *
	 * \param pos array that contain the coordinate
	 * \param ele element to store
	 *
	 */
	inline void addPad(const T (& pos)[dim], typename base::value_type ele)
	{
		// calculate the Cell id

		size_t cell_id = this->getCellPad(pos);

		// add the element to the cell

		addCell(cell_id,ele);
	}

	/*! \brief Add an element in the cell list forcing to be in the padding cells
	 *
	 * \warning careful is intended to be used ONLY to avoid round-off problems
	 *
	 * \param pos array that contain the coordinate
	 * \param ele element to store
	 *
	 */
	inline void addPad(const Point<dim,T> & pos, typename base::value_type ele)
	{
		// calculate the Cell id

		size_t cell_id = this->getCell(pos);

		// add the element to the cell

		addCell(cell_id,ele);
757
758
759
760
761
762
763
764
765
766
	}

	/*! \brief remove an element from the cell
	 *
	 * \param cell cell id
	 * \param ele element id
	 *
	 */
	inline void remove(size_t cell, size_t ele)
	{
767
		Mem_type::remove(cell,ele);
768
769
770
771
772
773
774
775
776
777
778
	}

	/*! \brief Return the number of elements in the cell
	 *
	 * \param cell_id id of the cell
	 *
	 * \return number of elements in the cell
	 *
	 */
	inline size_t getNelements(const size_t cell_id) const
	{
779
		return Mem_type::getNelements(cell_id);
780
781
782
783
784
785
786
787
788
789
790
791
	}

	/*! \brief Get an element in the cell
	 *
	 * \tparam i property to get
	 *
	 * \param cell cell id
	 * \param ele element id
	 *
	 * \return The element value
	 *
	 */
792
	inline auto get(size_t cell, size_t ele) -> decltype(this->Mem_type::get(cell,ele))
793
	{
794
		return Mem_type::get(cell,ele);
795
796
	}

incardon's avatar
incardon committed
797
798
799
800
801
802
803
804
805
806
	/*! \brief Get an element in the cell
	 *
	 * \tparam i property to get
	 *
	 * \param cell cell id
	 * \param ele element id
	 *
	 * \return The element value
	 *
	 */
incardon's avatar
incardon committed
807
	inline auto get(size_t cell, size_t ele) const -> decltype(this->Mem_type::get(cell,ele))
incardon's avatar
incardon committed
808
809
810
	{
		return Mem_type::get(cell,ele);
	}
811
812
813
814
815
816
817
818

	/*! \brief Swap the memory
	 *
	 * \param cl Cell list with witch you swap the memory
	 *
	 */
	inline void swap(CellList<dim,T,Mem_type,transform,base> & cl)
	{
incardon's avatar
incardon committed
819
820
		NNc_full.swap(cl.NNc_full);
		NNc_sym.swap(cl.NNc_sym);
821

822
		Mem_type::swap(static_cast<Mem_type &>(cl));
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909

		static_cast<CellDecomposer_sm<dim,T,transform> &>(*this) = static_cast<const CellDecomposer_sm<dim,T,transform> &>(cl);
	}

	/*! \brief Get the Cell iterator
	 *
	 * \param cell
	 *
	 * \return the iterator to the elements inside cell
	 *
	 */
	CellIterator<CellList<dim,T,Mem_type,transform,base>> getCellIterator(size_t cell)
	{
		return CellIterator<CellList<dim,T,Mem_type,transform,base>>(cell,*this);
	}

	/*! \brief Get the Neighborhood iterator
	 *
	 * It iterate across all the element of the selected cell and the near cells
	 *
	 *  \verbatim

	     * * *
	     * x *
	     * * *

	   \endverbatim
	 *
	 * * x is the selected cell
	 * * * are the near cell
	 *
	 * \param cell cell id
	 *
	 * \return An iterator across the neighhood particles
	 *
	 */
	template<unsigned int impl=NO_CHECK> inline CellNNIterator<dim,CellList<dim,T,Mem_type,transform,base>,FULL,impl> getNNIterator(size_t cell)
	{
		CellNNIterator<dim,CellList<dim,T,Mem_type,transform,base>,FULL,impl> cln(cell,NNc_full,*this);
		return cln;

	}

	/*! \brief Get the Neighborhood iterator
	 *
	 * It iterate across all the element of the selected cell and the near cells up to some selected radius
	 *
	 * \param cell cell id
	 * \param r_cut radius
	 *
	 * \return An iterator across the neighborhood particles
	 *
	 */
	template<unsigned int impl=NO_CHECK> inline CellNNIteratorRadius<dim,CellList<dim,T,Mem_type,transform,base>,impl> getNNIteratorRadius(size_t cell, T r_cut)
	{
		openfpm::vector<long int> & NNc = rcache[r_cut];

		if (NNc.size() == 0)
			NNcalc(r_cut,NNc);

		CellNNIteratorRadius<dim,CellList<dim,T,Mem_type,transform,base>,impl> cln(cell,NNc,*this);

		return cln;
	}

	/*! \brief Get the Neighborhood iterator
	 *
	 * It iterate across all the element of the selected cell and the near cells
	 *
	 *  \verbatim

	   * * *
	     x *

	   \endverbatim
	 *
	 * * x is the selected cell
	 * * * are the near cell
	 *
	 * \param cell cell id
	 * \param p particle id
	 *
	 * \return An aiterator across the neighborhood particles
	 *
	 */
	template<unsigned int impl> inline CellNNIteratorSym<dim,CellList<dim,T,Mem_type,transform,base>,SYM,impl> getNNIteratorSym(size_t cell, size_t p, const openfpm::vector<Point<dim,T>> & v)
	{
incardon's avatar
incardon committed
910
911
912
913
914
#ifdef SE_CLASS1
		if (from_cd == false)
			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;
#endif

915
916
917
918
		CellNNIteratorSym<dim,CellList<dim,T,Mem_type,transform,base>,SYM,impl> cln(cell,p,NNc_sym,*this,v);
		return cln;
	}

incardon's avatar
incardon committed
919
	/*! \brief Get the symmetric neighborhood
920
	 *
incardon's avatar
incardon committed
921
	 * \return the symmetric neighborhood
922
923
	 *
	 */
incardon's avatar
incardon committed
924
	const NNc_array<dim,(unsigned int)openfpm::math::pow(3,dim)/2+1> & getNNc_sym() const
925
	{
incardon's avatar
incardon committed
926
		return NNc_sym;
927
928
929
930
931
932
933
934
935
	}

	/*! \brief Return the number of padding cells of the Cell decomposer
	 *
	 * \param i dimension
	 *
	 * \return the number of padding cells
	 *
	 */
incardon's avatar
incardon committed
936
	size_t getPadding(size_t i) const
937
938
939
940
	{
		return CellDecomposer_sm<dim,T,transform>::getPadding(i);
	}

Yaroslav's avatar
Yaroslav committed
941

942
943
944
945
946
	/*! \brief Clear the cell list
	 *
	 */
	void clear()
	{
947
		Mem_type::clear();
948
949
	}

950
	/*! \brief Return the starting point of the cell p
951
	 *
952
	 * \param cell_id cell id
953
954
955
956
	 *
	 * \return the index
	 *
	 */
incardon's avatar
incardon committed
957
	inline const typename Mem_type::loc_index & getStartId(typename Mem_type::loc_index cell_id) const
958
	{
959
		return Mem_type::getStartId(cell_id);
960
961
	}

962
	/*! \brief Return the end point of the cell p
963
	 *
964
	 * \param cell_id cell id
965
966
967
968
	 *
	 * \return the stop index
	 *
	 */
incardon's avatar
incardon committed
969
	inline const typename Mem_type::loc_index & getStopId(typename Mem_type::loc_index cell_id) const
970
	{
971
		return Mem_type::getStopId(cell_id);
972
973
974
975
976
977
978
979
980
	}

	/*! \brief Return the neighborhood id
	 *
	 * \param part_id particle id
	 *
	 * \return the neighborhood id
	 *
	 */
incardon's avatar
incardon committed
981
	inline const typename Mem_type::loc_index & get_lin(const typename Mem_type::loc_index * part_id) const
982
	{
983
		return Mem_type::get_lin(part_id);
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
	}

//////////////////////////////// POINTLESS BUT REQUIRED TO RESPECT THE INTERFACE //////////////////

	//! Ghost marker
	size_t g_m = 0;

	/*! \brief return the ghost marker
	 *
	 * \return ghost marker
	 *
	 */
	inline size_t get_gm()
	{
		return g_m;
	}

	/*! \brief Set the ghost marker
	 *
	 * \param g_m marker
	 *
	 */
	inline void set_gm(size_t g_m)
	{
		this->g_m = g_m;
	}

/////////////////////////////////////
incardon's avatar
incardon committed
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035

	/////////////////////////////////////

		/*! \brief Set the n_dec number
		 *
		 * \param n_dec
		 *
		 */
		void set_ndec(size_t n_dec)
		{
			this->n_dec = n_dec;
		}

		/*! \brief Set the n_dec number
		 *
		 * \return n_dec
		 *
		 */
		size_t get_ndec() const
		{
			return n_dec;
		}

	/////////////////////////////////////
incardon's avatar
incardon committed
1036
1037
};

Pietro Incardona's avatar
Pietro Incardona committed
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
/*! \brief Calculate parameters for the cell list
 *
 * \param div Division array
 * \param r_cut interation radius or size of each cell
 * \param enlarge In case of padding particles the cell list must be enlarged, like a ghost. This parameter says how much must be enlarged
 *
 * \return the processor bounding box
 */
template<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)
{
	// calculate the parameters of the cell list

	// extend by the ghost
	pbox.enlarge(enlarge);

	// Calculate the division array and the cell box
	for (size_t i = 0; i < dim; i++)
	{
		div[i] = static_cast<size_t>((pbox.getP2().get(i) - pbox.getP1().get(i)) / r_cut);
		div[i]++;
		pbox.setHigh(i,pbox.getLow(i) + div[i]*r_cut);
	}
}

1062
1063
1064
1065
/*! \brief Calculate parameters for the symmetric cell list
 *
 * \param[in] dom Simulation domain
 * \param[output] cd_sm This cell-decomposer is set according to the needed division
incardon's avatar
incardon committed
1066
 * \param[in] g Ghost dimensions
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
 * \param[output] pad required padding for the cell-list
 *
 * \return the processor bounding box
 */
template<unsigned int dim, typename St> static inline void cl_param_calculateSym(const Box<dim,St> & dom, CellDecomposer_sm<dim,St,shift<dim,St>> & cd_sm, Ghost<dim,St> g, St r_cut, size_t & pad)
{
	size_t div[dim];

	for (size_t i = 0 ; i < dim ; i++)
		div[i] = (dom.getHigh(i) - dom.getLow(i)) / r_cut;

	g.magnify(1.013);

	// Calculate the maximum padding
	for (size_t i = 0 ; i < dim ; i++)
	{
		size_t tmp = std::ceil(fabs(g.getLow(i)) / r_cut);
		pad = (pad > tmp)?pad:tmp;
incardon's avatar
incardon committed
1085
1086
1087

		tmp = std::ceil(fabs(g.getHigh(i)) / r_cut);
		pad = (pad > tmp)?pad:tmp;
1088
1089
1090
1091
1092
	}

	cd_sm.setDimensions(dom,div,pad);
}

incardon's avatar
incardon committed
1093

incardon's avatar
incardon committed
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
/*! \brief Calculate parameters for the symmetric cell list
 *
 * \param[in] dom Simulation domain
 * \param[out] cd_sm This cell-decomposer is set according to the needed division
 * \param[in] g Ghost part extension
 * \param[out] pad required padding for the cell-list
 *
 * \return the processor bounding box
 */
template<unsigned int dim, typename St> static inline void cl_param_calculateSym(const Box<dim,St> & dom, CellDecomposer_sm<dim,St,shift<dim,St>> & cd_sm, Ghost<dim,St> g, size_t & pad)
{
	size_t div[dim];

	for (size_t i = 0 ; i < dim ; i++)
		div[i] = (dom.getHigh(i) - dom.getLow(i)) / g.getHigh(i);

	g.magnify(1.013);

	pad = 1;

	cd_sm.setDimensions(dom,div,pad);
}

incardon's avatar
incardon committed
1117
1118

#endif /* CELLLIST_HPP_ */