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

#ifndef VECTOR_HPP_
#define VECTOR_HPP_

incardon's avatar
incardon committed
11
#include "HDF5_wr/HDF5_wr.hpp"
incardon's avatar
incardon committed
12
#include "VCluster/VCluster.hpp"
incardon's avatar
incardon committed
13
#include "Space/Shape/Point.hpp"
incardon's avatar
incardon committed
14
#include "Vector/Iterators/vector_dist_iterator.hpp"
incardon's avatar
incardon committed
15
16
#include "Space/Shape/Box.hpp"
#include "Vector/vector_dist_key.hpp"
incardon's avatar
incardon committed
17
#include "memory/PtrMemory.hpp"
incardon's avatar
incardon committed
18
#include "NN/CellList/CellList.hpp"
incardon's avatar
incardon committed
19
#include "NN/CellList/CellListFast_gen.hpp"
20
#include "util/common.hpp"
incardon's avatar
incardon committed
21
#include "util/object_util.hpp"
incardon's avatar
incardon committed
22
#include "memory/ExtPreAlloc.hpp"
23
#include "CSVWriter/CSVWriter.hpp"
Pietro Incardona's avatar
Pietro Incardona committed
24
#include "VTKWriter/VTKWriter.hpp"
25
#include "Decomposition/common.hpp"
incardon's avatar
incardon committed
26
#include "Grid/Iterators/grid_dist_id_iterator_dec.hpp"
Yaroslav's avatar
Yaroslav committed
27
#include "Grid/grid_key_dx_iterator_hilbert.hpp"
Pietro Incardona's avatar
Pietro Incardona committed
28
#include "Vector/vector_dist_ofb.hpp"
29
#include "Decomposition/CartDecomposition.hpp"
30
#include "data_type/aggregate.hpp"
31
#include "NN/VerletList/VerletList.hpp"
Pietro Incardona's avatar
Pietro Incardona committed
32
#include "vector_dist_comm.hpp"
incardon's avatar
incardon committed
33
#include "DLB/LB_Model.hpp"
incardon's avatar
incardon committed
34
#include "Vector/vector_map_iterator.hpp"
incardon's avatar
incardon committed
35
36
#include "NN/CellList/ParticleIt_Cells.hpp"
#include "NN/CellList/ProcKeys.hpp"
incardon's avatar
incardon committed
37

incardon's avatar
incardon committed
38
39
#define DEC_GRAN(gr) ((size_t)gr << 32)

incardon's avatar
incardon committed
40
41
42
43
44
45
46
47
48
49
50
51
52
#define VECTOR_DIST_ERROR_OBJECT std::runtime_error("Runtime vector distributed error");

#ifdef SE_CLASS3
#include "se_class3_vector.hpp"
#endif

#ifdef SE_CLASS3
	#define SE_CLASS3_VDIST_CONSTRUCTOR ,se3(getDecomposition(),*this)
#else
	#define SE_CLASS3_VDIST_CONSTRUCTOR
#endif


incardon's avatar
incardon committed
53
54
55
#define NO_ID false
#define ID true

incardon's avatar
incardon committed
56
// Perform a ghost get or a ghost put
incardon's avatar
incardon committed
57
58
59
#define GET	1
#define PUT 2

incardon's avatar
incardon committed
60
61
62
63
// Write the particles with ghost
#define NO_GHOST 0
#define WITH_GHOST 2

incardon's avatar
incardon committed
64
65
66
67
#define GCL_NON_SYMMETRIC 0
#define GCL_SYMMETRIC 1
#define GCL_HILBERT 2

incardon's avatar
incardon committed
68

incardon's avatar
incardon committed
69
//! General function t get a cell-list
incardon's avatar
incardon committed
70
template<unsigned int dim, typename St, typename CellL, typename Vector, unsigned int impl>
incardon's avatar
incardon committed
71
72
73
74
75
76
77
78
79
80
81
82
83
struct gcl
{
	/*! \brief Get the Cell list based on the type
	 *
	 * \param vd Distributed vector
	 * \param r_cut Cut-off radius
	 * \param g Ghost
	 *
	 * \return the constructed cell-list
	 *
	 */
	static inline CellL get(Vector & vd, const St & r_cut, const Ghost<dim,St> & g)
	{
incardon's avatar
incardon committed
84
		return vd.template getCellList<CellL>(r_cut);
incardon's avatar
incardon committed
85
86
87
88
	}
};

//! General function t get a cell-list
incardon's avatar
incardon committed
89
90
template<unsigned int dim, typename St, typename CellL, typename Vector>
struct gcl<dim,St,CellL,Vector,GCL_HILBERT>
incardon's avatar
incardon committed
91
92
93
94
95
96
97
98
99
100
{
	/*! \brief Get the Cell list based on the type
	 *
	 * \param vd Distributed vector
	 * \param r_cut Cut-off radius
	 * \param g Ghost
	 *
	 * \return the constructed cell-list
	 *
	 */
incardon's avatar
incardon committed
101
	static inline CellL get(Vector & vd, const St & r_cut, const Ghost<dim,St> & g)
incardon's avatar
incardon committed
102
103
104
105
106
	{
		return vd.getCellList_hilb(r_cut,g);
	}
};

incardon's avatar
incardon committed
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
//! General function t get a cell-list
template<unsigned int dim, typename St, typename CellL, typename Vector>
struct gcl<dim,St,CellL,Vector,GCL_SYMMETRIC>
{
	/*! \brief Get the Cell list based on the type
	 *
	 * \param vd Distributed vector
	 * \param r_cut Cut-off radius
	 * \param g Ghost
	 *
	 * \return the constructed cell-list
	 *
	 */
	static inline CellL get(Vector & vd, const St & r_cut, const Ghost<dim,St> & g)
	{
		return vd.getCellListSym(r_cut);
	}
};

incardon's avatar
incardon committed
126
#define CELL_MEMFAST(dim,St) CellList_gen<dim, St, Process_keys_lin, Mem_fast<>, shift<dim, St> >
incardon's avatar
incardon committed
127
128
#define CELL_MEMBAL(dim,St) CellList_gen<dim, St, Process_keys_lin, Mem_bal<>, shift<dim, St> >
#define CELL_MEMMW(dim,St) CellList_gen<dim, St, Process_keys_lin, Mem_mw<>, shift<dim, St> >
incardon's avatar
incardon committed
129

incardon's avatar
incardon committed
130
#define CELL_MEMFAST_HILB(dim,St) CellList_gen<dim, St, Process_keys_hilb, Mem_fast<>, shift<dim, St> >
incardon's avatar
incardon committed
131
132
#define CELL_MEMBAL_HILB(dim,St) CellList_gen<dim, St, Process_keys_hilb, Mem_bal<>, shift<dim, St> >
#define CELL_MEMMW_HILB(dim,St) CellList_gen<dim, St, Process_keys_hilb, Mem_mw<>, shift<dim, St> >
incardon's avatar
incardon committed
133

incardon's avatar
incardon committed
134
135
136
137
138
139
140
141
142
143
144
145
146
147
#define VERLET_MEMFAST(dim,St) VerletList<dim,St,Mem_fast<>,shift<dim,St> >
#define VERLET_MEMBAL(dim,St)  VerletList<dim,St,Mem_bal<>,shift<dim,St> >
#define VERLET_MEMMW(dim,St)   VerletList<dim,St,Mem_mw<>,shift<dim,St> >

#define VERLET_MEMFAST_INT(dim,St) VerletList<dim,St,Mem_fast<unsigned int>,shift<dim,St> >
#define VERLET_MEMBAL_INT(dim,St)  VerletList<dim,St,Mem_bal<unsigned int>,shift<dim,St> >
#define VERLET_MEMMW_INT(dim,St)   VerletList<dim,St,Mem_mw<unsigned int>,shift<dim,St> >

enum reorder_opt
{
	NO_REORDER = 0,
	HILBERT = 1,
	LINEAR = 2
};
incardon's avatar
incardon committed
148

incardon's avatar
incardon committed
149
/*! \brief Distributed vector
150
 *
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
 * This class reppresent a distributed vector, the distribution of the structure
 * is based on the positional information of the elements the vector store
 *
 * ## Create a vector of random elements on each processor 2D
 * \snippet vector_dist_unit_test.hpp Create a vector of random elements on each processor 2D
 *
 * ## Create a vector of random elements on each processor 3D
 * \snippet vector_dist_unit_test.hpp Create a vector of random elements on each processor 3D
 *
 * ## Create a vector of elements distributed on a grid like way
 * \snippet vector_dist_unit_test.hpp Create a vector of elements distributed on a grid like way
 *
 * ## Redistribute the particles and sync the ghost properties
 * \snippet vector_dist_unit_test.hpp Redistribute the particles and sync the ghost properties
 *
 * \tparam dim Dimensionality of the space where the elements lives
 * \tparam St type of space float, double ...
 * \tparam prop properties the vector element store in OpenFPM data structure format
 * \tparam Decomposition Decomposition strategy to use CartDecomposition ...
 * \tparam Memory Memory pool where store the information HeapMemory ...
incardon's avatar
incardon committed
171
172
173
 *
 */

incardon's avatar
incardon committed
174
175
176
177
178
179
180
template<unsigned int dim,
         typename St,
		 typename prop,
		 typename layout = typename memory_traits_lin<prop>::type,
		 template <typename> class layout_base = memory_traits_lin,
		 typename Decomposition = CartDecomposition<dim,St>,
		 typename Memory = HeapMemory>
incardon's avatar
incardon committed
181
class vector_dist : public vector_dist_comm<dim,St,prop,layout,layout_base,Decomposition,Memory>
incardon's avatar
incardon committed
182
{
incardon's avatar
incardon committed
183
184
185
public:

	//! Self type
incardon's avatar
incardon committed
186
	typedef vector_dist<dim,St,prop,layout,layout_base,Decomposition,Memory> self;
incardon's avatar
incardon committed
187
188
189
190

	//! property object
	typedef prop value_type;

incardon's avatar
incardon committed
191
192
private:

193
	//! Ghost marker, all the particle with id > g_m are ghost all with g_m < are real particle
194
195
	size_t g_m = 0;

196
197
	//! Particle position vector, (It has 2 elements) the first has real particles assigned to a processor
	//! the second element contain unassigned particles
incardon's avatar
incardon committed
198
	openfpm::vector<Point<dim, St>,Memory> v_pos;
incardon's avatar
incardon committed
199

200
201
	//! Particle properties vector, (It has 2 elements) the first has real particles assigned to a processor
	//! the second element contain unassigned particles
incardon's avatar
incardon committed
202
	openfpm::vector<prop,Memory,layout,layout_base> v_prp;
incardon's avatar
incardon committed
203

204
	//! Virtual cluster
incardon's avatar
incardon committed
205
206
	Vcluster & v_cl;

incardon's avatar
incardon committed
207
208
209
	//! option used to create this vector
	size_t opt = 0;

210
211
212
	//! Name of the properties
	openfpm::vector<std::string> prp_names;

incardon's avatar
incardon committed
213
214
215
216
217
#ifdef SE_CLASS3

	se_class3_vector<prop::max_prop,dim,St,Decomposition,self> se3;

#endif
218

Pietro Incardona's avatar
Pietro Incardona committed
219
	/*! \brief Initialize the structures
incardon's avatar
incardon committed
220
	 *
Pietro Incardona's avatar
Pietro Incardona committed
221
	 * \param np number of particles
incardon's avatar
incardon committed
222
223
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
224
	void init_structures(size_t np)
incardon's avatar
incardon committed
225
	{
226
		// convert to a local number of elements
227
228
229
230
231
232
233
234
		size_t p_np = np / v_cl.getProcessingUnits();

		// Get non divisible part
		size_t r = np % v_cl.getProcessingUnits();

		// Distribute the remain particles
		if (v_cl.getProcessUnitID() < r)
			p_np++;
235

incardon's avatar
incardon committed
236
		// resize the position vector
237
		v_pos.resize(p_np);
incardon's avatar
incardon committed
238
239

		// resize the properties vector
240
241
242
		v_prp.resize(p_np);

		g_m = p_np;
Pietro Incardona's avatar
Pietro Incardona committed
243
	}
incardon's avatar
incardon committed
244

incardon's avatar
incardon committed
245
	/*! \brief Check if the parameters describe a valid vector. In case it does not report an error
incardon's avatar
incardon committed
246
247
248
249
250
251
252
253
	 *
	 * \param box Box to check
	 *
	 */
	void check_parameters(Box<dim,St> & box)
	{
		// if the box is not valid return an error
		if (box.isValid() == false)
incardon's avatar
incardon committed
254
		{
incardon's avatar
incardon committed
255
			std::cerr << __FILE__ << ":" << __LINE__ << "  Error the domain is not valid " << box.toString() << std::endl;
incardon's avatar
incardon committed
256
257
			ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
		}
incardon's avatar
incardon committed
258
259
260

	}

incardon's avatar
incardon committed
261
262
263
264
265
	/*! \brief It check that the r_cut is not bugger than the ghost
	 *
	 * \param r_cut cut-off radius
	 *
	 */
incardon's avatar
incardon committed
266
	void check_ghost_compatible_rcut(St r_cut)
incardon's avatar
incardon committed
267
268
269
	{
		for (size_t i = 0 ; i < dim ; i++)
		{
incardon's avatar
incardon committed
270
271
272
273
274
			if (fabs(getDecomposition().getGhost().getLow(i)) < r_cut)
			{
				std::cerr << __FILE__ << ":" << __LINE__ << " Error the cut off radius " << r_cut << " is bigger that the ghost layer on the dimension " << i << " lower=" << getDecomposition().getGhost().getLow(i) << std::endl;
				ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
			}
incardon's avatar
incardon committed
275
276
277
		}
	}

incardon's avatar
incardon committed
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
	/*! \brief Reorder based on hilbert space filling curve
	 *
	 * \param v_pos_dest reordered vector of position
	 * \param v_prp_dest reordered vector of properties
	 * \param m order of the space filling curve
	 * \param cell_list cell-list
	 *
	 */
	template<typename CellL, typename sfc_it>
	void reorder_sfc(openfpm::vector<Point<dim,St>> & v_pos_dest,
						 openfpm::vector<prop> & v_prp_dest,
						 sfc_it & h_it,
						 CellL & cell_list)
	{
		v_pos_dest.resize(v_pos.size());
		v_prp_dest.resize(v_prp.size());

		//Index for v_pos_dest
		size_t count = 0;

		grid_key_dx<dim> ksum;

		for (size_t i = 0; i < dim ; i++)
		{ksum.set_d(i,cell_list.getPadding(i));}

		while (h_it.isNext())
		{
		  auto key = h_it.get();
		  key += ksum;

		  size_t lin = cell_list.getGrid().LinId(key);

		  // for each particle in the Cell "lin"
		  for (size_t i = 0; i < cell_list.getNelements(lin); i++)
		  {
			  //reorder
			  auto v = cell_list.get(lin,i);
			  v_pos_dest.get(count) = v_pos.get(v);
			  v_prp_dest.get(count) = v_prp.get(v);

			  count++;
		  }
		  ++h_it;
		}
	}
incardon's avatar
incardon committed
323

Pietro Incardona's avatar
Pietro Incardona committed
324
325
326
327
328
329
330
331
public:

	//! space type
	typedef St stype;

	//! dimensions of space
	static const unsigned int dims = dim;

Pietro Incardona's avatar
Pietro Incardona committed
332
333
334
335
	/*! \brief Operator= for distributed vector
	 *
	 * \param v vector to copy
	 *
Pietro Incardona's avatar
Pietro Incardona committed
336
337
	 * \return itself
	 *
Pietro Incardona's avatar
Pietro Incardona committed
338
	 */
incardon's avatar
incardon committed
339
	vector_dist<dim,St,prop,layout,layout_base,Decomposition,Memory> & operator=(const vector_dist<dim,St,prop,layout,layout_base,Decomposition,Memory> & v)
Pietro Incardona's avatar
Pietro Incardona committed
340
	{
incardon's avatar
incardon committed
341
		static_cast<vector_dist_comm<dim,St,prop,layout,layout_base,Decomposition,Memory> *>(this)->operator=(static_cast<vector_dist_comm<dim,St,prop,layout,layout_base,Decomposition,Memory>>(v));
Pietro Incardona's avatar
Pietro Incardona committed
342
343
344
345
346

		g_m = v.g_m;
		v_pos = v.v_pos;
		v_prp = v.v_prp;

incardon's avatar
incardon committed
347
348
349
350
351
352
#ifdef SE_CLASS3
		se3 = v.se3;
#endif

		opt = v.opt;

Pietro Incardona's avatar
Pietro Incardona committed
353
354
355
356
357
358
359
		return *this;
	}

	/*! \brief Operator= for distributed vector
	 *
	 * \param v vector to copy
	 *
Pietro Incardona's avatar
Pietro Incardona committed
360
361
	 * \return itself
	 *
Pietro Incardona's avatar
Pietro Incardona committed
362
	 */
incardon's avatar
incardon committed
363
	vector_dist<dim,St,prop,layout,layout_base,Decomposition,Memory> & operator=(vector_dist<dim,St,prop,layout,layout_base,Decomposition,Memory> && v)
Pietro Incardona's avatar
Pietro Incardona committed
364
	{
incardon's avatar
incardon committed
365
		static_cast<vector_dist_comm<dim,St,prop,layout,layout_base,Decomposition,Memory> *>(this)->operator=(static_cast<vector_dist_comm<dim,St,prop,layout,layout_base,Decomposition,Memory> >(v));
Pietro Incardona's avatar
Pietro Incardona committed
366
367
368
369
370

		g_m = v.g_m;
		v_pos.swap(v.v_pos);
		v_prp.swap(v.v_prp);

incardon's avatar
incardon committed
371
372
373
374
375
376
#ifdef SE_CLASS3
		se3 = v.se3;
#endif

		opt = v.opt;

Pietro Incardona's avatar
Pietro Incardona committed
377
378
379
		return *this;
	}

incardon's avatar
incardon committed
380

Pietro Incardona's avatar
Pietro Incardona committed
381
	/*! \brief Copy Constructor
Pietro Incardona's avatar
Pietro Incardona committed
382
	 *
Pietro Incardona's avatar
Pietro Incardona committed
383
	 * \param v vector to copy
Pietro Incardona's avatar
Pietro Incardona committed
384
385
	 *
	 */
incardon's avatar
incardon committed
386
	vector_dist(const vector_dist<dim,St,prop,layout,layout_base,Decomposition,Memory> & v)
incardon's avatar
incardon committed
387
	:vector_dist_comm<dim,St,prop,layout,layout_base,Decomposition,Memory>(v.getDecomposition()),v_cl(v.v_cl) SE_CLASS3_VDIST_CONSTRUCTOR
Pietro Incardona's avatar
Pietro Incardona committed
388
389
390
391
392
393
394
395
	{
#ifdef SE_CLASS2
		check_new(this,8,VECTOR_DIST_EVENT,4);
#endif

		this->operator=(v);
	}

Pietro Incardona's avatar
Pietro Incardona committed
396
	/*! \brief Copy constructor
Pietro Incardona's avatar
Pietro Incardona committed
397
	 *
Pietro Incardona's avatar
Pietro Incardona committed
398
	 * \param v vector to copy
Pietro Incardona's avatar
Pietro Incardona committed
399
400
	 *
	 */
incardon's avatar
incardon committed
401
	vector_dist(vector_dist<dim,St,prop,layout,layout_base,Decomposition,Memory> && v) noexcept
incardon's avatar
incardon committed
402
	:v_cl(v.v_cl) SE_CLASS3_VDIST_CONSTRUCTOR
Pietro Incardona's avatar
Pietro Incardona committed
403
404
405
406
407
408
	{
#ifdef SE_CLASS2
		check_new(this,8,VECTOR_DIST_EVENT,4);
#endif

		this->operator=(v);
incardon's avatar
incardon committed
409
410
411
412

#ifdef SE_CLASS3
		se3.Initialize();
#endif
Pietro Incardona's avatar
Pietro Incardona committed
413
	}
Pietro Incardona's avatar
Pietro Incardona committed
414

Pietro Incardona's avatar
Pietro Incardona committed
415
	/*! \brief Constructor with predefined decomposition
Pietro Incardona's avatar
Pietro Incardona committed
416
	 *
Pietro Incardona's avatar
Pietro Incardona committed
417
418
	 * \param dec is the decomposition
	 * \param np number of particles
Pietro Incardona's avatar
Pietro Incardona committed
419
420
421
	 *
	 */
	vector_dist(const Decomposition & dec, size_t np) :
incardon's avatar
incardon committed
422
	vector_dist_comm<dim,St,prop,layout,layout_base,Decomposition,Memory>(dec), v_cl(create_vcluster()) SE_CLASS3_VDIST_CONSTRUCTOR
Pietro Incardona's avatar
Pietro Incardona committed
423
424
425
426
427
428
	{
#ifdef SE_CLASS2
		check_new(this,8,VECTOR_DIST_EVENT,4);
#endif

		init_structures(np);
incardon's avatar
incardon committed
429
430
431
432

#ifdef SE_CLASS3
		se3.Initialize();
#endif
Pietro Incardona's avatar
Pietro Incardona committed
433
434
435
	}


incardon's avatar
incardon committed
436
	/*! \brief Constructor of a distributed vector
Pietro Incardona's avatar
Pietro Incardona committed
437
438
439
	 *
	 * \param np number of elements
	 * \param box domain where the vector of elements live
Pietro Incardona's avatar
Pietro Incardona committed
440
	 * \param bc boundary conditions
Pietro Incardona's avatar
Pietro Incardona committed
441
	 * \param g Ghost margins
incardon's avatar
incardon committed
442
	 * \param opt [Optional] additional options. BIND_DEC_TO_GHOST Bind the decomposition to be multiple of the
443
444
	 *          ghost size. This is required if we want to use symmetric to eliminate
	 *          ghost communications.
incardon's avatar
incardon committed
445
	 * \param gdist [Optional] override the default distribution grid
Pietro Incardona's avatar
Pietro Incardona committed
446
447
	 *
	 */
448
	vector_dist(size_t np, Box<dim, St> box, const size_t (&bc)[dim], const Ghost<dim, St> & g, size_t opt = 0, const grid_sm<dim,void> & gdist = grid_sm<dim,void>())
incardon's avatar
incardon committed
449
	:v_cl(create_vcluster()),opt(opt) SE_CLASS3_VDIST_CONSTRUCTOR
Pietro Incardona's avatar
Pietro Incardona committed
450
451
452
453
454
	{
#ifdef SE_CLASS2
		check_new(this,8,VECTOR_DIST_EVENT,4);
#endif

incardon's avatar
incardon committed
455
456
457
		if (opt >> 32 != 0)
			this->setDecompositionGranularity(opt >> 32);

incardon's avatar
incardon committed
458
459
		check_parameters(box);

Pietro Incardona's avatar
Pietro Incardona committed
460
		init_structures(np);
461
		this->init_decomposition(box,bc,g,opt,gdist);
incardon's avatar
incardon committed
462
463
464
465

#ifdef SE_CLASS3
		se3.Initialize();
#endif
Pietro Incardona's avatar
Pietro Incardona committed
466
467
	}

468
469
470
471
472
473
474
	~vector_dist()
	{
#ifdef SE_CLASS2
		check_delete(this);
#endif
	}

475
476
477
478
479
480
481
482
483
	/*! \brief remove all the elements
	 *
	 *
	 */
	void clear()
	{
		resize(0);
	}

Pietro Incardona's avatar
Pietro Incardona committed
484
	/*! \brief return the local size of the vector
485
	 *
Pietro Incardona's avatar
Pietro Incardona committed
486
	 * \return local size
487
488
	 *
	 */
incardon's avatar
incardon committed
489
	size_t size_local() const
490
	{
Pietro Incardona's avatar
Pietro Incardona committed
491
		return g_m;
492
493
	}

incardon's avatar
incardon committed
494
495
496
497
498
	/*! \brief return the local size of the vector
	 *
	 * \return local size
	 *
	 */
incardon's avatar
incardon committed
499
	size_t size_local_with_ghost() const
incardon's avatar
incardon committed
500
	{
Pietro Incardona's avatar
Pietro Incardona committed
501
		return v_pos.size();
incardon's avatar
incardon committed
502
503
	}

incardon's avatar
incardon committed
504
505
#ifndef ONLY_READWRITE_GETTER

506
	/*! \brief Get the position of an element
incardon's avatar
incardon committed
507
	 *
508
509
510
511
512
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \param vec_key element
	 *
	 * \return the position of the element in space
incardon's avatar
incardon committed
513
514
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
515
	inline auto getPos(vect_dist_key_dx vec_key) -> decltype(v_pos.template get<0>(vec_key.getKey()))
incardon's avatar
incardon committed
516
	{
incardon's avatar
incardon committed
517
518
519
520
#ifdef SE_CLASS3
		check_for_pos_nan_inf<prop::max_prop_real,prop::max_prop>(*this,vec_key.getKey());
#endif

Pietro Incardona's avatar
Pietro Incardona committed
521
		return v_pos.template get<0>(vec_key.getKey());
incardon's avatar
incardon committed
522
523
	}

524
525
526
527
528
529
530
531
532
533
534
	/*! \brief Get the position of an element
	 *
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \param vec_key element
	 *
	 * \return the position of the element in space
	 *
	 */
	inline auto getPos(vect_dist_key_dx vec_key) const -> decltype(v_pos.template get<0>(vec_key.getKey()))
	{
incardon's avatar
incardon committed
535
536
537
#ifdef SE_CLASS3
		check_for_pos_nan_inf<prop::max_prop_real,prop::max_prop>(*this,vec_key.getKey());
#endif
538
539
540
		return v_pos.template get<0>(vec_key.getKey());
	}

541
542
543
544
545
546
547
548
549
	/*! \brief Get the position of an element
	 *
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \param vec_key element
	 *
	 * \return the position of the element in space
	 *
	 */
incardon's avatar
incardon committed
550
551
	inline auto getPos(size_t vec_key) -> decltype(v_pos.template get<0>(vec_key))
	{
incardon's avatar
incardon committed
552
553
554
#ifdef SE_CLASS3
		check_for_pos_nan_inf<prop::max_prop_real,prop::max_prop>(*this,vec_key);
#endif
incardon's avatar
incardon committed
555
556
557
558
559
560
561
562
563
564
565
566
		return v_pos.template get<0>(vec_key);
	}

	/*! \brief Get the position of an element
	 *
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \param vec_key element
	 *
	 * \return the position of the element in space
	 *
	 */
incardon's avatar
incardon committed
567
	inline auto getPos(size_t vec_key) const -> decltype(v_pos.template get<0>(vec_key))
568
	{
incardon's avatar
incardon committed
569
570
571
#ifdef SE_CLASS3
		check_for_pos_nan_inf<prop::max_prop_real,prop::max_prop>(*this,vec_key);
#endif
572
573
574
		return v_pos.template get<0>(vec_key);
	}

575
	/*! \brief Get the property of an element
incardon's avatar
incardon committed
576
	 *
577
578
579
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \tparam id property id
incardon's avatar
incardon committed
580
581
	 * \param vec_key vector element
	 *
582
583
	 * \return return the selected property of the vector element
	 *
incardon's avatar
incardon committed
584
	 */
585
	template<unsigned int id> inline auto getProp(vect_dist_key_dx vec_key) -> decltype(v_prp.template get<id>(vec_key.getKey()))
incardon's avatar
incardon committed
586
	{
incardon's avatar
incardon committed
587
588
589
#ifdef SE_CLASS3
		check_for_prop_nan_inf<id,prop::max_prop+SE3_STATUS>(*this,vec_key.getKey());
#endif
590
		return v_prp.template get<id>(vec_key.getKey());
incardon's avatar
incardon committed
591
592
	}

Pietro Incardona's avatar
Pietro Incardona committed
593
594
595
596
597
598
599
600
601
602
	/*! \brief Get the property of an element
	 *
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \tparam id property id
	 * \param vec_key vector element
	 *
	 * \return return the selected property of the vector element
	 *
	 */
incardon's avatar
incardon committed
603
	template<unsigned int id> inline auto getProp(vect_dist_key_dx vec_key) const -> decltype(v_prp.template get<id>(vec_key.getKey()))
Pietro Incardona's avatar
Pietro Incardona committed
604
	{
incardon's avatar
incardon committed
605
606
607
#ifdef SE_CLASS3
		check_for_prop_nan_inf<id,prop::max_prop+SE3_STATUS>(*this,vec_key.getKey());
#endif
Pietro Incardona's avatar
Pietro Incardona committed
608
609
610
		return v_prp.template get<id>(vec_key.getKey());
	}

incardon's avatar
incardon committed
611
612
613
614
615
616
617
618
619
620
621
622
	/*! \brief Get the property of an element
	 *
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \tparam id property id
	 * \param vec_key vector element
	 *
	 * \return return the selected property of the vector element
	 *
	 */
	template<unsigned int id> inline auto getProp(size_t vec_key) -> decltype(v_prp.template get<id>(vec_key))
	{
incardon's avatar
incardon committed
623
624
625
#ifdef SE_CLASS3
		check_for_prop_nan_inf<id,prop::max_prop+SE3_STATUS>(*this,vec_key);
#endif
incardon's avatar
incardon committed
626
627
628
		return v_prp.template get<id>(vec_key);
	}

629
630
631
632
633
634
635
636
637
638
	/*! \brief Get the property of an element
	 *
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \tparam id property id
	 * \param vec_key vector element
	 *
	 * \return return the selected property of the vector element
	 *
	 */
incardon's avatar
incardon committed
639
	template<unsigned int id> inline auto getProp(size_t vec_key) const -> decltype(v_prp.template get<id>(vec_key))
640
	{
incardon's avatar
incardon committed
641
642
643
#ifdef SE_CLASS3
		check_for_prop_nan_inf<id,prop::max_prop+SE3_STATUS>(*this,vec_key);
#endif
644
645
646
		return v_prp.template get<id>(vec_key);
	}

incardon's avatar
incardon committed
647
648
#endif

incardon's avatar
incardon committed
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
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
///////////////////// Read and write with no check

	/*! \brief Get the position of an element
	 *
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \param vec_key element
	 *
	 * \return the position of the element in space
	 *
	 */
	inline auto getPosNC(vect_dist_key_dx vec_key) -> decltype(v_pos.template get<0>(vec_key.getKey()))
	{
		return v_pos.template get<0>(vec_key.getKey());
	}

	/*! \brief Get the position of an element
	 *
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \param vec_key element
	 *
	 * \return the position of the element in space
	 *
	 */
	inline auto getPosNC(vect_dist_key_dx vec_key) const -> decltype(v_pos.template get<0>(vec_key.getKey()))
	{
		return v_pos.template get<0>(vec_key.getKey());
	}

	/*! \brief Get the position of an element
	 *
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \param vec_key element
	 *
	 * \return the position of the element in space
	 *
	 */
	inline auto getPosNC(size_t vec_key) -> decltype(v_pos.template get<0>(vec_key))
	{
		return v_pos.template get<0>(vec_key);
	}

	/*! \brief Get the position of an element
	 *
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \param vec_key element
	 *
	 * \return the position of the element in space
	 *
	 */
	inline auto getPosNC(size_t vec_key) const -> decltype(v_pos.template get<0>(vec_key))
	{
		return v_pos.template get<0>(vec_key);
	}

	/*! \brief Get the property of an element
	 *
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \tparam id property id
	 * \param vec_key vector element
	 *
	 * \return return the selected property of the vector element
	 *
	 */
	template<unsigned int id> inline auto getPropNC(vect_dist_key_dx vec_key) -> decltype(v_prp.template get<id>(vec_key.getKey()))
	{
		return v_prp.template get<id>(vec_key.getKey());
	}

	/*! \brief Get the property of an element
	 *
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \tparam id property id
	 * \param vec_key vector element
	 *
	 * \return return the selected property of the vector element
	 *
	 */
incardon's avatar
incardon committed
732
	template<unsigned int id> inline auto getPropNC(vect_dist_key_dx vec_key) const -> decltype(v_prp.template get<id>(vec_key.getKey()))
incardon's avatar
incardon committed
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
	{
		return v_prp.template get<id>(vec_key.getKey());
	}

	/*! \brief Get the property of an element
	 *
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \tparam id property id
	 * \param vec_key vector element
	 *
	 * \return return the selected property of the vector element
	 *
	 */
	template<unsigned int id> inline auto getPropNC(size_t vec_key) -> decltype(v_prp.template get<id>(vec_key))
	{
		return v_prp.template get<id>(vec_key);
	}

	/*! \brief Get the property of an element
	 *
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \tparam id property id
	 * \param vec_key vector element
	 *
	 * \return return the selected property of the vector element
	 *
	 */
incardon's avatar
incardon committed
762
	template<unsigned int id> inline auto getPropNC(size_t vec_key) const -> decltype(v_prp.template get<id>(vec_key))
incardon's avatar
incardon committed
763
764
765
766
	{
		return v_prp.template get<id>(vec_key);
	}

incardon's avatar
incardon committed
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
///////////////////// Read and Write function

	/*! \brief Get the position of an element
	 *
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \param vec_key element
	 *
	 * \return the position of the element in space
	 *
	 */
	inline auto getPosWrite(vect_dist_key_dx vec_key) -> decltype(v_pos.template get<0>(vec_key.getKey()))
	{
#ifdef SE_CLASS3
		se3.template write<prop::max_prop_real>(*this,vec_key.getKey());
#endif

		return v_pos.template get<0>(vec_key.getKey());
	}

	/*! \brief Get the position of an element
	 *
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \param vec_key element
	 *
	 * \return the position of the element in space
	 *
	 */
	inline auto getPosRead(vect_dist_key_dx vec_key) const -> decltype(v_pos.template get<0>(vec_key.getKey()))
	{
#ifdef SE_CLASS3
		se3.template read<prop::max_prop_real>(*this,vec_key.getKey());
#endif

		return v_pos.template get<0>(vec_key.getKey());
	}

	/*! \brief Get the property of an element
	 *
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \tparam id property id
	 * \param vec_key vector element
	 *
	 * \return return the selected property of the vector element
	 *
	 */
	template<unsigned int id> inline auto getPropWrite(vect_dist_key_dx vec_key) -> decltype(v_prp.template get<id>(vec_key.getKey()))
	{
#ifdef SE_CLASS3
		se3.template write<id>(*this,vec_key.getKey());
#endif

		return v_prp.template get<id>(vec_key.getKey());
	}

	/*! \brief Get the property of an element
	 *
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \tparam id property id
	 * \param vec_key vector element
	 *
	 * \return return the selected property of the vector element
	 *
	 */
incardon's avatar
incardon committed
834
	template<unsigned int id> inline auto getPropRead(vect_dist_key_dx vec_key) const -> decltype(v_prp.template get<id>(vec_key.getKey()))
incardon's avatar
incardon committed
835
836
837
838
839
840
841
842
843
844
	{
#ifdef SE_CLASS3
		se3.template read<id>(*this,vec_key.getKey());
#endif

		return v_prp.template get<id>(vec_key.getKey());
	}

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

845
846
847
848
849
850
851
852
853
854
855
856
857
858
	/*! \brief Add local particle
	 *
	 * It add a local particle, with "local" we mean in this processor
	 * the particle can be also created out of the processor domain, in this
	 * case a call to map is required. Added particles are always created at the
	 * end and can be accessed with getLastPos and getLastProp
	 *
	 */
	void add()
	{
		v_prp.insert(g_m);
		v_pos.insert(g_m);

		g_m++;
incardon's avatar
incardon committed
859
860
861
862
863

#ifdef SE_CLASS3
		for (size_t i = 0 ; i < prop::max_prop_real+1 ; i++)
			v_prp.template get<prop::max_prop_real>(g_m-1)[i] = UNINITIALIZED;
#endif
864
865
	}

incardon's avatar
incardon committed
866
867
#ifndef ONLY_READWRITE_GETTER

868
869
870
871
872
	/*! \brief Get the position of the last element
	 *
	 * \return the position of the element in space
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
873
	inline auto getLastPos() -> decltype(v_pos.template get<0>(0))
874
	{
Pietro Incardona's avatar
Pietro Incardona committed
875
		return v_pos.template get<0>(g_m - 1);
876
877
878
879
880
881
882
883
884
885
886
	}

	/*! \brief Get the property of the last element
	 *
	 * \tparam id property id
	 *
	 * \return return the selected property of the vector element
	 *
	 */
	template<unsigned int id> inline auto getLastProp() -> decltype(v_prp.template get<id>(0))
	{
887
		return v_prp.template get<id>(g_m - 1);
888
889
	}

incardon's avatar
incardon committed
890
891
#endif

incardon's avatar
incardon committed
892
893
894
895
896
897
898
899
900
901
////////////////////////////// READ AND WRITE VERSION //////////

	/*! \brief Get the position of the last element
	 *
	 * \return the position of the element in space
	 *
	 */
	inline auto getLastPosRead() -> decltype(v_pos.template get<0>(0))
	{
#ifdef SE_CLASS3
incardon's avatar
incardon committed
902
		se3.template read<prop::max_prop_real>(*this,g_m-1);
incardon's avatar
incardon committed
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
#endif

		return v_pos.template get<0>(g_m - 1);
	}

	/*! \brief Get the property of the last element
	 *
	 * \tparam id property id
	 *
	 * \return return the selected property of the vector element
	 *
	 */
	template<unsigned int id> inline auto getLastPropRead() -> decltype(v_prp.template get<id>(0))
	{
#ifdef SE_CLASS3
		se3.read<id>(*this,g_m-1);
#endif

		return v_prp.template get<id>(g_m - 1);
	}


	/*! \brief Get the position of the last element
	 *
	 * \return the position of the element in space
	 *
	 */
	inline auto getLastPosWrite() -> decltype(v_pos.template get<0>(0))
	{
#ifdef SE_CLASS3
incardon's avatar
incardon committed
933
		se3.template write<prop::max_prop_real>(*this,g_m-1);
incardon's avatar
incardon committed
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
#endif

		return v_pos.template get<0>(g_m - 1);
	}

	/*! \brief Get the property of the last element
	 *
	 * \tparam id property id
	 *
	 * \return return the selected property of the vector element
	 *
	 */
	template<unsigned int id> inline auto getLastPropWrite() -> decltype(v_prp.template get<id>(0))
	{
#ifdef SE_CLASS3
incardon's avatar
incardon committed
949
		se3.template write<id>(*this,g_m-1);
incardon's avatar
incardon committed
950
951
952
953
954
955
956
#endif

		return v_prp.template get<id>(g_m - 1);
	}

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

incardon's avatar
incardon committed
957
958
959
960
961
962
963
964
965
	/*! \brief Construct a cell list symmetric based on a cut of radius
	 *
	 * \tparam CellL CellList type to construct
	 *
	 * \param r_cut interation radius, or size of each cell
	 *
	 * \return the Cell list
	 *
	 */
incardon's avatar
incardon committed
966
	template<typename CellL = CellList<dim, St, Mem_fast<>, shift<dim, St> > > CellL getCellListSym(St r_cut)
incardon's avatar
incardon committed
967
	{
incardon's avatar
incardon committed
968
#ifdef SE_CLASS1
incardon's avatar
incardon committed
969
		if (!(opt & BIND_DEC_TO_GHOST))
incardon's avatar
incardon committed
970
		{
incardon's avatar
incardon committed
971
			if (getDecomposition().getGhost().getLow(dim-1) == 0.0)
972
973
974
975
			{
				std::cerr << __FILE__ << ":" << __LINE__ << " Error the vector has been constructed without BIND_DEC_TO_GHOST, If you construct a vector without BIND_DEC_TO_GHOST the ghost must be full without reductions " << std::endl;
				ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
			}
incardon's avatar
incardon committed
976
977
978
		}
#endif

incardon's avatar
incardon committed
979
980
981
982
		// Cell list
		CellL cell_list;

		size_t pad = 0;
incardon's avatar
incardon committed
983
984
		CellDecomposer_sm<dim,St,shift<dim,St>> cd_sm;
		cl_param_calculateSym(getDecomposition().getDomain(),cd_sm,getDecomposition().getGhost(),r_cut,pad);
incardon's avatar
incardon committed
985

incardon's avatar
incardon committed
986
		// Processor bounding box
incardon's avatar
incardon committed
987
988
		Box<dim, St> pbox = getDecomposition().getProcessorBounds();

incardon's avatar
incardon committed
989
990
991
		// Ghost padding extension
		Ghost<dim,size_t> g_ext(0);
		cell_list.Initialize(cd_sm,pbox,pad);
992
		cell_list.set_ndec(getDecomposition().get_ndec());
incardon's avatar
incardon committed
993

incardon's avatar
incardon committed
994
		updateCellListSym(cell_list);
incardon's avatar
incardon committed
995
996

		return cell_list;
incardon's avatar
incardon committed
997
998
	}

999

Pietro Incardona's avatar
Pietro Incardona committed
1000
1001
1002
1003
	/*! \brief Construct a cell list starting from the stored particles
	 *
	 * \tparam CellL CellList type to construct
	 *
1004
	 * \param r_cut interation radius, or size of each cell
incardon's avatar
incardon committed
1005
	 * \param no_se3 avoid SE_CLASS3 checking
1006
1007
1008
	 *
	 * \return the Cell list
	 *
Pietro Incardona's avatar
Pietro Incardona committed
1009
	 */
incardon's avatar
incardon committed
1010
	template<typename CellL = CellList_gen<dim, St, Process_keys_lin, Mem_fast<>, shift<dim, St> > >
incardon's avatar
incardon committed
1011
	CellL getCellList(St r_cut, bool no_se3 = false)
1012
	{
incardon's avatar
incardon committed
1013
#ifdef SE_CLASS3
incardon's avatar
incardon committed
1014
1015
		if (no_se3 == false)
		{se3.getNN();}
incardon's avatar
incardon committed
1016
1017
1018
1019
1020
#endif
#ifdef SE_CLASS1
		check_ghost_compatible_rcut(r_cut);
#endif

Pietro Incardona's avatar
Pietro Incardona committed
1021
		// Get ghost and anlarge by 1%
Pietro Incardona's avatar
Pietro Incardona committed
1022
		Ghost<dim,St> g = getDecomposition().getGhost();
1023
		g.magnify(1.013);
Pietro Incardona's avatar
Pietro Incardona committed
1024

incardon's avatar
incardon committed
1025
		return getCellList<CellL>(r_cut, g,no_se3);
1026
1027
	}

1028
1029
1030
1031
1032
1033
1034
1035
1036
	/*! \brief Construct an hilbert cell list starting from the stored particles
	 *
	 * \tparam CellL CellList type to construct
	 *
	 * \param r_cut interation radius, or size of each cell
	 *
	 * \return the Cell list
	 *
	 */
incardon's avatar
incardon committed
1037
	template<typename CellL = CellList_gen<dim, St, Process_keys_hilb, Mem_fast<>, shift<dim, St> > >
incardon's avatar
incardon committed
1038
	CellL getCellList_hilb(St r_cut)
1039
	{
incardon's avatar
incardon committed
1040
1041
1042
1043
1044
1045
1046
#ifdef SE_CLASS3
		se3.getNN();
#endif
#ifdef SE_CLASS1
		check_ghost_compatible_rcut(r_cut);
#endif

1047
		// Get ghost and anlarge by 1%
Pietro Incardona's avatar
Pietro Incardona committed
1048
		Ghost<dim,St> g = getDecomposition().getGhost();
Pietro Incardona's avatar
Pietro Incardona committed
1049
		g.magnify(1.013);
1050
1051
1052
1053

		return getCellList_hilb(r_cut, g);
	}

Pietro Incardona's avatar
Pietro Incardona committed
1054
1055
1056
1057
1058
	/*! \brief Update a cell list using the stored particles
	 *
	 * \tparam CellL CellList type to construct
	 *
	 * \param cell_list Cell list to update
incardon's avatar
incardon committed
1059
	 * \param no_se3 avoid se class 3 checking
Pietro Incardona's avatar
Pietro Incardona committed
1060
1061
	 *
	 */
incardon's avatar
incardon committed
1062
	template<typename CellL> void updateCellList(CellL & cell_list, bool no_se3 = false)
Pietro Incardona's avatar
Pietro Incardona committed
1063
	{
incardon's avatar
incardon committed
1064
#ifdef SE_CLASS3
incardon's avatar
incardon committed
1065
1066
		if (no_se3 == false)
		{se3.getNN();}
incardon's avatar
incardon committed
1067
1068
#endif

incardon's avatar
incardon committed
1069
1070
1071
1072
1073
1074
		// This function assume equal spacing in all directions
		// but in the worst case we take the maximum
		St r_cut = 0;
		for (size_t i = 0 ; i < dim ; i++)
			r_cut = std::max(r_cut,cell_list.getCellBox().getHigh(i));

incardon's avatar
incardon committed
1075
1076
1077
		// Here we have to check that the Cell-list has been constructed
		// from the same decomposition
		bool to_reconstruct = cell_list.get_ndec() != getDecomposition().get_ndec();
incardon's avatar
incardon committed
1078
1079
1080
1081
1082
1083
1084
1085
1086

		if (to_reconstruct == false)
		{
			populate_cell_list(v_pos,cell_list,g_m,CL_NON_SYMMETRIC);

			cell_list.set_gm(g_m);
		}
		else
		{
incardon's avatar
incardon committed
1087
			CellL cli_tmp = gcl<dim,St,CellL,self,GCL_NON_SYMMETRIC>::get(*this,r_cut,getDecomposition().getGhost());
incardon's avatar
incardon committed
1088
1089
1090

			cell_list.swap(cli_tmp);
		}
incardon's avatar
incardon committed
1091
1092
1093
1094
1095
1096
1097
1098
1099
	}

	/*! \brief Update a cell list using the stored particles
	 *
	 * \tparam CellL CellList type to construct
	 *
	 * \param cell_list Cell list to update
	 *
	 */
incardon's avatar
incardon committed
1100
	template<typename CellL = CellList<dim, St, Mem_fast<>, shift<dim, St> > > void updateCellListSym(CellL & cell_list)
incardon's avatar
incardon committed
1101
	{
incardon's avatar
incardon committed
1102
1103
1104
#ifdef SE_CLASS3
		se3.getNN();
#endif
incardon's avatar
incardon committed
1105

incardon's avatar
incardon committed
1106
1107
1108
1109
		// This function assume equal spacing in all directions
		// but in the worst case we take the maximum
		St r_cut = 0;
		for (size_t i = 0 ; i < dim ; i++)
incardon's avatar
incardon committed
1110
		{r_cut = std::max(r_cut,cell_list.getCellBox().getHigh(i));}
incardon's avatar
incardon committed
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123

		// Here we have to check that the Cell-list has been constructed
		// from the same decomposition
		bool to_reconstruct = cell_list.get_ndec() != getDecomposition().get_ndec();

		if (to_reconstruct == false)
		{
			populate_cell_list(v_pos,cell_list,g_m,CL_SYMMETRIC);

			cell_list.set_gm(g_m);
		}
		else
		{
incardon's avatar
incardon committed
1124
			CellL cli_tmp = gcl<dim,St,CellL,self,GCL_SYMMETRIC>::get(*this,r_cut,getDecomposition().getGhost());
incardon's avatar
incardon committed
1125
1126
1127

			cell_list.swap(cli_tmp);
		}
Pietro Incardona's avatar
Pietro Incardona committed
1128
1129
	}

1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
	/*! \brief Construct a cell list starting from the stored particles
	 *
	 * It differ from the get getCellList for an additional parameter, in case the
	 * domain + ghost is not big enough to contain additional padding particles, a Cell list
	 * with bigger space can be created
	 * (padding particles in general are particles added by the user out of the domains)
	 *
	 * \tparam CellL CellList type to construct
	 *
	 * \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 say how much must be enlarged
incardon's avatar
incardon committed
1141
	 * \param no_se3 avoid se_class3 cheking default false
1142
	 *
Pietro Incardona's avatar
Pietro Incardona committed
1143
1144
	 * \return the CellList
	 *
1145
	 */
incardon's avatar
incardon committed
1146
	template<typename CellL = CellList_gen<dim, St, Process_keys_lin, Mem_fast<>, shift<dim, St> > >
incardon's avatar
incardon committed
1147
	CellL getCellList(St r_cut, const Ghost<dim, St> & enlarge, bool no_se3 = false)
Pietro Incardona's avatar
Pietro Incardona committed
1148
	{
incardon's avatar
incardon committed
1149
#ifdef SE_CLASS3
incardon's avatar
incardon committed
1150
1151
		if (no_se3 == false)
		{se3.getNN();}
incardon's avatar
incardon committed
1152
1153
#endif

Pietro Incardona's avatar
Pietro Incardona committed
1154
1155
		CellL cell_list;

1156
1157
		// Division array
		size_t div[dim];
1158

1159
		// get the processor bounding box
Pietro Incardona's avatar
Pietro Incardona committed
1160
		Box<dim, St> pbox = getDecomposition().getProcessorBounds();
1161

1162
		// Processor bounding box
1163
		cl_param_calculate(pbox, div, r_cut, enlarge);
1164

incardon's avatar
incardon committed
1165
1166
		cell_list.Initialize(pbox, div);
		cell_list.set_gm(g_m);
incardon's avatar
incardon committed
1167
		cell_list.set_ndec(getDecomposition().get_ndec());
1168

incardon's avatar
incardon committed
1169
		updateCellList(cell_list,no_se3);
1170
1171
1172

		return cell_list;
	}
1173

1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
	/*! \brief Construct an hilbert cell list starting from the stored particles
	 *
	 * It differ from the get getCellList for an additional parameter, in case the
	 * domain + ghost is not big enough to contain additional padding particles, a Cell list
	 * with bigger space can be created
	 * (padding particles in general are particles added by the user out of the domains)
	 *
	 * \tparam CellL CellList type to construct
	 *
	 * \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 say how much must be enlarged
	 *
Pietro Incardona's avatar
Pietro Incardona committed
1186
1187
	 * \return The Cell-list
	 *
1188
	 */
incardon's avatar
incardon committed
1189
	template<typename CellL = CellList_gen<dim, St, Process_keys_hilb, Mem_fast<>, shift<dim, St> > > CellL getCellList_hilb(St r_cut, const Ghost<dim, St> & enlarge)
1190
	{
incardon's avatar
incardon committed
1191
1192
1193
1194
#ifdef SE_CLASS3
		se3.getNN();
#endif

1195
1196
1197
		CellL cell_list;

		// Division array
1198
1199
		size_t div[dim];

1200
		// get the processor bounding box
Pietro Incardona's avatar
Pietro Incardona committed
1201
		Box<dim, St> pbox = getDecomposition().getProcessorBounds();
1202

1203
		// Processor bounding box
1204
		cl_param_calculate(pbox,div, r_cut, enlarge);
1205

incardon's avatar
incardon committed
1206
1207
		cell_list.Initialize(pbox, div);
		cell_list.set_gm(g_m);
incardon's avatar
incardon committed
1208
		cell_list.set_ndec(getDecomposition().get_ndec());
1209

Pietro Incardona's avatar
Pietro Incardona committed
1210
		updateCellList(cell_list);
Pietro Incardona's avatar
Pietro Incardona committed
1211
1212
1213

		return cell_list;
	}
incardon's avatar
incardon committed
1214

incardon's avatar
incardon committed
1215
1216
1217
1218
	/*! \brief for each particle get the symmetric verlet list
	 *
	 * \param r_cut cut-off radius
	 *
incardon's avatar
incardon committed
1219
1220
	 * \return the verlet list
	 *
incardon's avatar
incardon committed
1221
	 */
incardon's avatar
incardon committed
1222
1223
	template <typename VerletL = VerletList<dim,St,Mem_fast<>,shift<dim,St> >>
	VerletL getVerletSym(St r_cut)
incardon's avatar
incardon committed
1224
	{
incardon's avatar
incardon committed
1225
1226
1227
1228
#ifdef SE_CLASS3
		se3.getNN();
#endif

incardon's avatar
incardon committed
1229
		VerletL ver;
1230

incardon's avatar
incardon committed
1231
1232
1233
1234
1235
		// Processor bounding box
		Box<dim, St> pbox = getDecomposition().getProcessorBounds();

		ver.InitializeSym(getDecomposition().getDomain(),pbox,getDecomposition().getGhost(),r_cut,v_pos,g_m);

incardon's avatar
incardon committed
1236
1237
		ver.set_ndec(getDecomposition().get_ndec());

incardon's avatar
incardon committed
1238
1239
		return ver;
	}
Pietro Incardona's avatar
Pietro Incardona committed
1240

incardon's avatar
incardon committed
1241
1242
1243
1244
1245
1246
1247
	/*! \brief for each particle get the symmetric verlet list
	 *
	 * \param r_cut cut-off radius
	 *
	 * \return the verlet list
	 *
	 */
incardon's avatar
incardon committed
1248
1249
	template <typename VerletL = VerletList<dim,St,Mem_fast<>,shift<dim,St> >>
	VerletL getVerletCrs(St r_cut)
incardon's avatar
incardon committed
1250
	{
1251
#ifdef SE_CLASS1
incardon's avatar
incardon committed
1252
		if (!(opt & BIND_DEC_TO_GHOST))
1253
1254
1255
1256
1257
1258
		{
			std::cerr << __FILE__ << ":" << __LINE__ << " Error the vector has been constructed without BIND_DEC_TO_GHOST, getVerletCrs require the vector to be constructed with BIND_DEC_TO_GHOST option " << std::endl;
			ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
		}
#endif

incardon's avatar
incardon committed
1259
1260
1261
1262
#ifdef SE_CLASS3
		se3.getNN();
#endif

incardon's avatar
incardon committed
1263
		VerletL ver;
incardon's avatar
incardon committed
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274

		// Processor bounding box
		Box<dim, St> pbox = getDecomposition().getProcessorBounds();

		// Initialize the verlet list
		ver.InitializeCrs(getDecomposition().getDomain(),pbox,getDecomposition().getGhost(),<