vector_dist.hpp 64.4 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 "config.h"
incardon's avatar
incardon committed
12
#include "HDF5_wr/HDF5_wr.hpp"
incardon's avatar
incardon committed
13
#include "VCluster/VCluster.hpp"
incardon's avatar
incardon committed
14
#include "Space/Shape/Point.hpp"
incardon's avatar
incardon committed
15
#include "Vector/Iterators/vector_dist_iterator.hpp"
incardon's avatar
incardon committed
16 17
#include "Space/Shape/Box.hpp"
#include "Vector/vector_dist_key.hpp"
incardon's avatar
incardon committed
18
#include "memory/PtrMemory.hpp"
incardon's avatar
incardon committed
19
#include "NN/CellList/CellList.hpp"
incardon's avatar
incardon committed
20
#include "NN/CellList/CellListFast_gen.hpp"
21
#include "util/common.hpp"
incardon's avatar
incardon committed
22
#include "util/object_util.hpp"
incardon's avatar
incardon committed
23
#include "memory/ExtPreAlloc.hpp"
24
#include "CSVWriter/CSVWriter.hpp"
25
#include "VTKWriter/VTKWriter.hpp"
26
#include "Decomposition/common.hpp"
incardon's avatar
incardon committed
27
#include "Grid/Iterators/grid_dist_id_iterator_dec.hpp"
28
#include "Grid/grid_key_dx_iterator_hilbert.hpp"
29
#include "Vector/vector_dist_ofb.hpp"
30
#include "Decomposition/CartDecomposition.hpp"
31
#include "data_type/aggregate.hpp"
32
#include "NN/VerletList/VerletList.hpp"
33
#include "vector_dist_comm.hpp"
incardon's avatar
incardon committed
34
#include "DLB/LB_Model.hpp"
incardon's avatar
incardon committed
35
#include "Vector/vector_map_iterator.hpp"
incardon's avatar
incardon committed
36 37
#include "NN/CellList/ParticleIt_Cells.hpp"
#include "NN/CellList/ProcKeys.hpp"
incardon's avatar
incardon committed
38 39
#include "Vector/vector_dist_kernel.hpp"
#include "NN/CellList/cuda/CellList_gpu.hpp"
incardon's avatar
incardon committed
40

incardon's avatar
incardon committed
41 42
#define DEC_GRAN(gr) ((size_t)gr << 32)

incardon's avatar
incardon committed
43 44 45 46 47 48 49 50 51 52 53 54 55
#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
56 57 58
#define NO_ID false
#define ID true

incardon's avatar
incardon committed
59
// Perform a ghost get or a ghost put
incardon's avatar
incardon committed
60 61 62
#define GET	1
#define PUT 2

incardon's avatar
incardon committed
63 64 65 66
// Write the particles with ghost
#define NO_GHOST 0
#define WITH_GHOST 2

incardon's avatar
incardon committed
67 68 69 70
#define GCL_NON_SYMMETRIC 0
#define GCL_SYMMETRIC 1
#define GCL_HILBERT 2

incardon's avatar
incardon committed
71

incardon's avatar
incardon committed
72
//! General function t get a cell-list
incardon's avatar
incardon committed
73
template<unsigned int dim, typename St, typename CellL, typename Vector, unsigned int impl>
incardon's avatar
incardon committed
74 75 76 77 78 79 80 81 82 83 84 85 86
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
87
		return vd.template getCellList<CellL>(r_cut);
incardon's avatar
incardon committed
88 89 90 91
	}
};

//! General function t get a cell-list
incardon's avatar
incardon committed
92 93
template<unsigned int dim, typename St, typename CellL, typename Vector>
struct gcl<dim,St,CellL,Vector,GCL_HILBERT>
incardon's avatar
incardon committed
94 95 96 97 98 99 100 101 102 103
{
	/*! \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
104
	static inline CellL get(Vector & vd, const St & r_cut, const Ghost<dim,St> & g)
incardon's avatar
incardon committed
105 106 107 108 109
	{
		return vd.getCellList_hilb(r_cut,g);
	}
};

incardon's avatar
incardon committed
110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128
//! 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
129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149
/////////////////// GCL anisotropic ///////////////////////////////////////////

//! General function t get a cell-list
template<unsigned int dim, typename St, typename CellL, typename Vector, unsigned int impl>
struct gcl_An
{
	/*! \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 size_t (& div)[dim], const size_t (& pad)[dim], const Ghost<dim,St> & g)
	{
		return vd.template getCellListSym<CellL>(div,pad);
	}
};

150
#define CELL_MEMFAST(dim,St) CellList_gen<dim, St, Process_keys_lin, Mem_fast<>, shift<dim, St> >
151 152
#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
153

154
#define CELL_MEMFAST_HILB(dim,St) CellList_gen<dim, St, Process_keys_hilb, Mem_fast<>, shift<dim, St> >
155 156
#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
157

158 159 160 161
#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> >

incardon's avatar
incardon committed
162
#define VERLET_MEMFAST_INT(dim,St) VerletList<dim,St,Mem_fast<HeapMemory,unsigned int>,shift<dim,St> >
163 164 165 166 167 168 169 170 171
#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
172

incardon's avatar
incardon committed
173
/*! \brief Distributed vector
174
 *
175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194
 * 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 ...
195
 * \tparam Memory layout
incardon's avatar
incardon committed
196 197 198
 *
 */

199 200 201 202 203 204 205
template<unsigned int dim,
         typename St,
         typename prop,
         typename Decomposition = CartDecomposition<dim,St>,
         typename Memory = HeapMemory,
         template<typename> class layout_base = memory_traits_lin>
class vector_dist : public vector_dist_comm<dim,St,prop,Decomposition,Memory,layout_base>
incardon's avatar
incardon committed
206
{
incardon's avatar
incardon committed
207 208 209
public:

	//! Self type
incardon's avatar
incardon committed
210
	typedef vector_dist<dim,St,prop,Decomposition,Memory,layout_base> self;
incardon's avatar
incardon committed
211 212 213 214

	//! property object
	typedef prop value_type;

incardon's avatar
incardon committed
215 216
private:

217
	//! Ghost marker, all the particle with id > g_m are ghost all with g_m < are real particle
218 219
	size_t g_m = 0;

220 221
	//! Particle position vector, (It has 2 elements) the first has real particles assigned to a processor
	//! the second element contain unassigned particles
222
	openfpm::vector<Point<dim, St>,Memory,typename layout_base<Point<dim,St>>::type,layout_base> v_pos;
incardon's avatar
incardon committed
223

224 225
	//! Particle properties vector, (It has 2 elements) the first has real particles assigned to a processor
	//! the second element contain unassigned particles
226
	openfpm::vector<prop,Memory,typename layout_base<prop>::type,layout_base> v_prp;
incardon's avatar
incardon committed
227

incardon's avatar
incardon committed
228 229 230 231 232 233 234 235
#ifdef CUDA_GPU

	openfpm::vector<prop,Memory,typename layout_base<prop>::type,layout_base> v_prp_out;

	openfpm::vector<Point<dim, St>,Memory,typename layout_base<Point<dim,St>>::type,layout_base> v_pos_out;

#endif

236
	//! Virtual cluster
incardon's avatar
incardon committed
237 238
	Vcluster & v_cl;

incardon's avatar
incardon committed
239 240 241
	//! option used to create this vector
	size_t opt = 0;

242 243 244
	//! Name of the properties
	openfpm::vector<std::string> prp_names;

incardon's avatar
incardon committed
245 246 247 248 249
#ifdef SE_CLASS3

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

#endif
250

Pietro Incardona's avatar
Pietro Incardona committed
251
	/*! \brief Initialize the structures
incardon's avatar
incardon committed
252
	 *
Pietro Incardona's avatar
Pietro Incardona committed
253
	 * \param np number of particles
incardon's avatar
incardon committed
254 255
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
256
	void init_structures(size_t np)
incardon's avatar
incardon committed
257
	{
258
		// convert to a local number of elements
259 260 261 262 263 264 265 266
		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++;
267

incardon's avatar
incardon committed
268
		// resize the position vector
269
		v_pos.resize(p_np);
incardon's avatar
incardon committed
270 271

		// resize the properties vector
272 273 274
		v_prp.resize(p_np);

		g_m = p_np;
Pietro Incardona's avatar
Pietro Incardona committed
275
	}
incardon's avatar
incardon committed
276

incardon's avatar
incardon committed
277
	/*! \brief Check if the parameters describe a valid vector. In case it does not report an error
278 279 280 281 282 283 284 285
	 *
	 * \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
286
		{
287
			std::cerr << __FILE__ << ":" << __LINE__ << "  Error the domain is not valid " << box.toString() << std::endl;
incardon's avatar
incardon committed
288 289
			ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
		}
290 291 292

	}

incardon's avatar
incardon committed
293 294 295 296 297
	/*! \brief It check that the r_cut is not bugger than the ghost
	 *
	 * \param r_cut cut-off radius
	 *
	 */
incardon's avatar
incardon committed
298
	void check_ghost_compatible_rcut(St r_cut)
incardon's avatar
incardon committed
299 300 301
	{
		for (size_t i = 0 ; i < dim ; i++)
		{
incardon's avatar
incardon committed
302 303 304 305 306
			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
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
	/*! \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;
		}
	}
355

Pietro Incardona's avatar
Pietro Incardona committed
356 357 358 359 360 361 362 363
public:

	//! space type
	typedef St stype;

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

364 365 366 367
	/*! \brief Operator= for distributed vector
	 *
	 * \param v vector to copy
	 *
Pietro Incardona's avatar
Pietro Incardona committed
368 369
	 * \return itself
	 *
370
	 */
incardon's avatar
incardon committed
371 372
	vector_dist<dim,St,prop,Decomposition,Memory,layout_base> &
	operator=(const vector_dist<dim,St,prop,Decomposition,Memory,layout_base> & v)
373
	{
incardon's avatar
incardon committed
374
		static_cast<vector_dist_comm<dim,St,prop,Decomposition,Memory,layout_base> *>(this)->operator=(static_cast<vector_dist_comm<dim,St,prop,Decomposition,Memory,layout_base>>(v));
375 376 377 378 379

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

incardon's avatar
incardon committed
380 381 382 383 384 385
#ifdef SE_CLASS3
		se3 = v.se3;
#endif

		opt = v.opt;

386 387 388 389 390 391 392
		return *this;
	}

	/*! \brief Operator= for distributed vector
	 *
	 * \param v vector to copy
	 *
Pietro Incardona's avatar
Pietro Incardona committed
393 394
	 * \return itself
	 *
395
	 */
incardon's avatar
incardon committed
396 397
	vector_dist<dim,St,prop,Decomposition,Memory,layout_base> &
	operator=(vector_dist<dim,St,prop,Decomposition,Memory,layout_base> && v)
398
	{
incardon's avatar
incardon committed
399
		static_cast<vector_dist_comm<dim,St,prop,Decomposition,Memory,layout_base> *>(this)->operator=(static_cast<vector_dist_comm<dim,St,prop,Decomposition,Memory,layout_base> >(v));
400 401 402 403 404

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

incardon's avatar
incardon committed
405 406 407 408 409 410
#ifdef SE_CLASS3
		se3 = v.se3;
#endif

		opt = v.opt;

411 412 413
		return *this;
	}

incardon's avatar
incardon committed
414

Pietro Incardona's avatar
Pietro Incardona committed
415
	/*! \brief Copy Constructor
416
	 *
Pietro Incardona's avatar
Pietro Incardona committed
417
	 * \param v vector to copy
418 419
	 *
	 */
incardon's avatar
incardon committed
420 421
	vector_dist(const vector_dist<dim,St,prop,Decomposition,Memory,layout_base> & v)
	:vector_dist_comm<dim,St,prop,Decomposition,Memory,layout_base>(v.getDecomposition()),v_cl(v.v_cl) SE_CLASS3_VDIST_CONSTRUCTOR
422 423 424 425 426 427 428 429
	{
#ifdef SE_CLASS2
		check_new(this,8,VECTOR_DIST_EVENT,4);
#endif

		this->operator=(v);
	}

Pietro Incardona's avatar
Pietro Incardona committed
430
	/*! \brief Copy constructor
431
	 *
Pietro Incardona's avatar
Pietro Incardona committed
432
	 * \param v vector to copy
433 434
	 *
	 */
incardon's avatar
incardon committed
435
	vector_dist(vector_dist<dim,St,prop,Decomposition,Memory,layout_base> && v) noexcept
incardon's avatar
incardon committed
436
	:v_cl(v.v_cl) SE_CLASS3_VDIST_CONSTRUCTOR
437 438 439 440 441 442
	{
#ifdef SE_CLASS2
		check_new(this,8,VECTOR_DIST_EVENT,4);
#endif

		this->operator=(v);
incardon's avatar
incardon committed
443 444 445 446

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

Pietro Incardona's avatar
Pietro Incardona committed
449
	/*! \brief Constructor with predefined decomposition
Pietro Incardona's avatar
Pietro Incardona committed
450
	 *
Pietro Incardona's avatar
Pietro Incardona committed
451 452
	 * \param dec is the decomposition
	 * \param np number of particles
Pietro Incardona's avatar
Pietro Incardona committed
453 454 455
	 *
	 */
	vector_dist(const Decomposition & dec, size_t np) :
incardon's avatar
incardon committed
456
	vector_dist_comm<dim,St,prop,Decomposition,Memory,layout_base>(dec), v_cl(create_vcluster()) SE_CLASS3_VDIST_CONSTRUCTOR
Pietro Incardona's avatar
Pietro Incardona committed
457 458 459 460 461 462
	{
#ifdef SE_CLASS2
		check_new(this,8,VECTOR_DIST_EVENT,4);
#endif

		init_structures(np);
incardon's avatar
incardon committed
463 464 465 466

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


incardon's avatar
incardon committed
470
	/*! \brief Constructor of a distributed vector
Pietro Incardona's avatar
Pietro Incardona committed
471 472 473
	 *
	 * \param np number of elements
	 * \param box domain where the vector of elements live
Pietro Incardona's avatar
Pietro Incardona committed
474
	 * \param bc boundary conditions
Pietro Incardona's avatar
Pietro Incardona committed
475
	 * \param g Ghost margins
476
	 * \param opt [Optional] additional options. BIND_DEC_TO_GHOST Bind the decomposition to be multiple of the
477 478
	 *          ghost size. This is required if we want to use symmetric to eliminate
	 *          ghost communications.
479
	 * \param gdist [Optional] override the default distribution grid
Pietro Incardona's avatar
Pietro Incardona committed
480 481
	 *
	 */
482
	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
483
	:v_cl(create_vcluster()),opt(opt) SE_CLASS3_VDIST_CONSTRUCTOR
Pietro Incardona's avatar
Pietro Incardona committed
484 485 486 487 488
	{
#ifdef SE_CLASS2
		check_new(this,8,VECTOR_DIST_EVENT,4);
#endif

incardon's avatar
incardon committed
489 490 491
		if (opt >> 32 != 0)
			this->setDecompositionGranularity(opt >> 32);

492 493
		check_parameters(box);

Pietro Incardona's avatar
Pietro Incardona committed
494
		init_structures(np);
495
		this->init_decomposition(box,bc,g,opt,gdist);
incardon's avatar
incardon committed
496 497 498 499

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

502 503 504 505 506 507 508
	~vector_dist()
	{
#ifdef SE_CLASS2
		check_delete(this);
#endif
	}

509 510 511 512 513 514 515 516 517
	/*! \brief remove all the elements
	 *
	 *
	 */
	void clear()
	{
		resize(0);
	}

518
	/*! \brief return the local size of the vector
519
	 *
520
	 * \return local size
521 522
	 *
	 */
incardon's avatar
incardon committed
523
	size_t size_local() const
524
	{
525
		return g_m;
526 527
	}

incardon's avatar
incardon committed
528 529 530 531 532
	/*! \brief return the local size of the vector
	 *
	 * \return local size
	 *
	 */
incardon's avatar
incardon committed
533
	size_t size_local_with_ghost() const
incardon's avatar
incardon committed
534
	{
535
		return v_pos.size();
incardon's avatar
incardon committed
536 537
	}

incardon's avatar
incardon committed
538 539
#ifndef ONLY_READWRITE_GETTER

540
	/*! \brief Get the position of an element
incardon's avatar
incardon committed
541
	 *
542 543 544 545 546
	 * 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
547 548
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
549
	inline auto getPos(vect_dist_key_dx vec_key) -> decltype(v_pos.template get<0>(vec_key.getKey()))
incardon's avatar
incardon committed
550
	{
551 552 553 554
#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
555
		return v_pos.template get<0>(vec_key.getKey());
incardon's avatar
incardon committed
556 557
	}

558 559 560 561 562 563 564 565 566 567 568
	/*! \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()))
	{
569 570 571
#ifdef SE_CLASS3
		check_for_pos_nan_inf<prop::max_prop_real,prop::max_prop>(*this,vec_key.getKey());
#endif
572 573 574
		return v_pos.template get<0>(vec_key.getKey());
	}

575 576 577 578 579 580 581 582 583
	/*! \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
	 *
	 */
584 585
	inline auto getPos(size_t vec_key) -> decltype(v_pos.template get<0>(vec_key))
	{
586 587 588
#ifdef SE_CLASS3
		check_for_pos_nan_inf<prop::max_prop_real,prop::max_prop>(*this,vec_key);
#endif
589 590 591 592 593 594 595 596 597 598 599 600
		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
601
	inline auto getPos(size_t vec_key) const -> decltype(v_pos.template get<0>(vec_key))
602
	{
603 604 605
#ifdef SE_CLASS3
		check_for_pos_nan_inf<prop::max_prop_real,prop::max_prop>(*this,vec_key);
#endif
606 607 608
		return v_pos.template get<0>(vec_key);
	}

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

627 628 629 630 631 632 633 634 635 636
	/*! \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
637
	template<unsigned int id> inline auto getProp(vect_dist_key_dx vec_key) const -> decltype(v_prp.template get<id>(vec_key.getKey()))
638
	{
639 640 641
#ifdef SE_CLASS3
		check_for_prop_nan_inf<id,prop::max_prop+SE3_STATUS>(*this,vec_key.getKey());
#endif
642 643 644
		return v_prp.template get<id>(vec_key.getKey());
	}

645 646 647 648 649 650 651 652 653 654 655 656
	/*! \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))
	{
657 658 659
#ifdef SE_CLASS3
		check_for_prop_nan_inf<id,prop::max_prop+SE3_STATUS>(*this,vec_key);
#endif
660 661 662
		return v_prp.template get<id>(vec_key);
	}

663 664 665 666 667 668 669 670 671 672
	/*! \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
673
	template<unsigned int id> inline auto getProp(size_t vec_key) const -> decltype(v_prp.template get<id>(vec_key))
674
	{
675 676 677
#ifdef SE_CLASS3
		check_for_prop_nan_inf<id,prop::max_prop+SE3_STATUS>(*this,vec_key);
#endif
678 679 680
		return v_prp.template get<id>(vec_key);
	}

incardon's avatar
incardon committed
681 682
#endif

incardon's avatar
incardon committed
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 757 758 759 760 761 762 763 764 765
///////////////////// 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
766
	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
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
	{
		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
796
	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
797 798 799 800
	{
		return v_prp.template get<id>(vec_key);
	}

incardon's avatar
incardon committed
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 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
///////////////////// 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
868
	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
869 870 871 872 873 874 875 876 877 878
	{
#ifdef SE_CLASS3
		se3.template read<id>(*this,vec_key.getKey());
#endif

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

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

879 880 881 882 883 884 885 886 887 888 889 890 891 892
	/*! \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
893 894 895 896 897

#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
898 899
	}

900 901
#ifndef ONLY_READWRITE_GETTER

902 903 904 905 906
	/*! \brief Get the position of the last element
	 *
	 * \return the position of the element in space
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
907
	inline auto getLastPos() -> decltype(v_pos.template get<0>(0))
908
	{
Pietro Incardona's avatar
Pietro Incardona committed
909
		return v_pos.template get<0>(g_m - 1);
910 911 912 913 914 915 916 917 918 919 920
	}

	/*! \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))
	{
921
		return v_prp.template get<id>(g_m - 1);
922 923
	}

924 925
#endif

incardon's avatar
incardon committed
926 927 928 929 930 931 932 933 934 935
////////////////////////////// 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
936
		se3.template read<prop::max_prop_real>(*this,g_m-1);
incardon's avatar
incardon committed
937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966
#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
967
		se3.template write<prop::max_prop_real>(*this,g_m-1);
incardon's avatar
incardon committed
968 969 970 971 972 973 974 975 976 977 978 979 980 981 982
#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
983
		se3.template write<id>(*this,g_m-1);
incardon's avatar
incardon committed
984 985 986 987 988 989 990
#endif

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

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

991 992 993 994 995 996 997 998 999
	/*! \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
	 *
	 */
1000
	template<typename CellL = CellList<dim, St, Mem_fast<>, shift<dim, St> > > CellL getCellListSym(St r_cut)
1001
	{
incardon's avatar
incardon committed
1002
#ifdef SE_CLASS1
incardon's avatar
incardon committed
1003
		if (!(opt & BIND_DEC_TO_GHOST))
incardon's avatar
incardon committed
1004
		{
incardon's avatar
incardon committed
1005
			if (getDecomposition().getGhost().getLow(dim-1) == 0.0)
1006 1007 1008 1009
			{
				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
1010 1011 1012
		}
#endif

1013 1014 1015 1016
		// Cell list
		CellL cell_list;

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

incardon's avatar
incardon committed
1020
		// Processor bounding box
1021 1022
		Box<dim, St> pbox = getDecomposition().getProcessorBounds();

incardon's avatar
incardon committed
1023 1024 1025
		// Ghost padding extension
		Ghost<dim,size_t> g_ext(0);
		cell_list.Initialize(cd_sm,pbox,pad);
1026
		cell_list.set_ndec(getDecomposition().get_ndec());
1027

1028
		updateCellListSym(cell_list);
incardon's avatar
incardon committed
1029 1030

		return cell_list;
1031 1032
	}

incardon's avatar
incardon committed
1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078
	/*! \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
	 *
	 */
	template<typename CellL = CellList<dim, St, Mem_fast<>, shift<dim, St> > >
	CellL getCellListSym(const size_t (& div)[dim],
						 const size_t (& pad)[dim])
	{
#ifdef SE_CLASS1
		if (!(opt & BIND_DEC_TO_GHOST))
		{
			if (getDecomposition().getGhost().getLow(dim-1) == 0.0)
			{
				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);
			}
		}
#endif

		size_t pad_max = pad[0];
		for (size_t i = 1 ; i < dim ; i++)
		{if (pad[i] > pad_max)	{pad_max = pad[i];}}

		// Cell list
		CellL cell_list;

		CellDecomposer_sm<dim,St,shift<dim,St>> cd_sm;
		cd_sm.setDimensions(getDecomposition().getDomain(),div,pad_max);

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

		// Ghost padding extension
		Ghost<dim,size_t> g_ext(0);
		cell_list.Initialize(cd_sm,pbox,pad_max);
		cell_list.set_ndec(getDecomposition().get_ndec());

		updateCellListSym(cell_list);

		return cell_list;
	}
1079

1080 1081 1082 1083
	/*! \brief Construct a cell list starting from the stored particles
	 *
	 * \tparam CellL CellList type to construct
	 *
1084
	 * \param r_cut interation radius, or size of each cell
incardon's avatar
incardon committed
1085
	 * \param no_se3 avoid SE_CLASS3 checking
1086 1087 1088
	 *
	 * \return the Cell list
	 *
1089
	 */
1090
	template<typename CellL = CellList_gen<dim, St, Process_keys_lin, Mem_fast<>, shift<dim, St> > >
incardon's avatar
incardon committed
1091
	CellL getCellList(St r_cut, bool no_se3 = false)
1092
	{
incardon's avatar
incardon committed
1093
#ifdef SE_CLASS3
1094 1095
		if (no_se3 == false)
		{se3.getNN();}
incardon's avatar
incardon committed
1096 1097 1098 1099 1100
#endif
#ifdef SE_CLASS1
		check_ghost_compatible_rcut(r_cut);
#endif

Pietro Incardona's avatar
Pietro Incardona committed
1101
		// Get ghost and anlarge by 1%
1102
		Ghost<dim,St> g = getDecomposition().getGhost();
1103
		g.magnify(1.013);
Pietro Incardona's avatar
Pietro Incardona committed
1104

1105
		return getCellList<CellL>(r_cut, g,no_se3);
1106 1107
	}

incardon's avatar
incardon committed
1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177
#ifdef CUDA_GPU

	/*! \brief Construct a cell list starting from the stored particles, this version of cell-list can be offloaded with
	 *         the function toGPU
	 *
	 * \param r_cut interation radius, or size of each cell
	 *
	 * \return the Cell list
	 *
	 */
	CellList_gpu<dim,St,CudaMemory,shift_only<dim, St>> getCellListGPU(St r_cut, bool no_se3 = false)
	{
#ifdef SE_CLASS3
		if (no_se3 == false)
		{se3.getNN();}
#endif
#ifdef SE_CLASS1
		check_ghost_compatible_rcut(r_cut);
#endif

		// Get ghost and anlarge by 1%
		Ghost<dim,St> g = getDecomposition().getGhost();
		g.magnify(1.013);

		return getCellListGPU(r_cut, g,no_se3);
	}

	/*! \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
	 *
	 * \return the CellList
	 *
	 */
	CellList_gpu<dim,St,CudaMemory,shift_only<dim, St>> getCellListGPU(St r_cut, const Ghost<dim, St> & enlarge, bool no_se3 = false)
	{
#ifdef SE_CLASS3
		if (no_se3 == false)
		{se3.getNN();}
#endif

		// Division array
		size_t div[dim];

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

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

		CellList_gpu<dim,St,CudaMemory,shift_only<dim, St>> cell_list(pbox,div);

		v_prp_out.resize(v_pos.size());
		v_pos_out.resize(v_pos.size());

		cell_list.template construct<decltype(v_pos),decltype(v_prp)>(v_pos,v_pos_out,v_prp,v_prp_out);

		return cell_list;
	}

#endif

1178 1179 1180 1181 1182 1183 1184 1185 1186
	/*! \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
	 *
	 */
1187
	template<typename CellL = CellList_gen<dim, St, Process_keys_hilb, Mem_fast<>, shift<dim, St> > >
incardon's avatar
incardon committed
1188
	CellL getCellList_hilb(St r_cut)
1189
	{
incardon's avatar
incardon committed
1190 1191 1192 1193 1194 1195 1196
#ifdef SE_CLASS3
		se3.getNN();
#endif
#ifdef SE_CLASS1
		check_ghost_compatible_rcut(r_cut);
#endif

1197
		// Get ghost and anlarge by 1%
1198
		Ghost<dim,St> g = getDecomposition().getGhost();
Pietro Incardona's avatar
Pietro Incardona committed
1199
		g.magnify(1.013);
1200 1201 1202 1203

		return getCellList_hilb(r_cut, g);
	}

Pietro Incardona's avatar
Pietro Incardona committed
1204 1205 1206 1207 1208
	/*! \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
1209
	 * \param no_se3 avoid se class 3 checking
Pietro Incardona's avatar
Pietro Incardona committed
1210 1211
	 *
	 */
1212
	template<typename CellL> void updateCellList(CellL & cell_list, bool no_se3 = false)
Pietro Incardona's avatar
Pietro Incardona committed
1213
	{
incardon's avatar
incardon committed
1214
#ifdef SE_CLASS3
1215 1216
		if (no_se3 == false)
		{se3.getNN();}
incardon's avatar
incardon committed
1217 1218
#endif

incardon's avatar
incardon committed
1219 1220 1221 1222
		// 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
1223
		{r_cut = std::max(r_cut,cell_list.getCellBox().getHigh(i));}
incardon's avatar
incardon committed
1224

incardon's avatar
incardon committed
1225 1226 1227
		// 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
1228 1229 1230 1231 1232 1233 1234 1235 1236

		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
1237
			CellL cli_tmp = gcl<dim,St,CellL,self,GCL_NON_SYMMETRIC>::get(*this,r_cut,getDecomposition().getGhost());
incardon's avatar
incardon committed
1238 1239 1240

			cell_list.swap(cli_tmp);
		}
incardon's avatar
incardon committed
1241 1242 1243 1244 1245 1246 1247 1248 1249
	}

	/*! \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
1250 1251
	template<typename CellL = CellList<dim, St, Mem_fast<>, shift<dim, St> > >
	void updateCellListSym(CellL & cell_list)
incardon's avatar
incardon committed
1252
	{
incardon's avatar
incardon committed
1253 1254 1255
#ifdef SE_CLASS3
		se3.getNN();
#endif
incardon's avatar
incardon committed
1256

incardon's avatar
incardon committed
1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268
		// 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
1269 1270 1271 1272
			CellL cli_tmp = gcl_An<dim,St,CellL,self,GCL_SYMMETRIC>::get(*this,
																		 cell_list.getDivWP(),
																		 cell_list.getPadding(),
																		 getDecomposition().getGhost());
incardon's avatar
incardon committed
1273 1274 1275

			cell_list.swap(cli_tmp);
		}
Pietro Incardona's avatar
Pietro Incardona committed
1276 1277
	}

1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288
	/*! \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
1289
	 * \param no_se3 avoid se_class3 cheking default false
1290
	 *
Pietro Incardona's avatar
Pietro Incardona committed
1291 1292
	 * \return the CellList
	 *
1293
	 */
1294
	template<typename CellL = CellList_gen<dim, St, Process_keys_lin, Mem_fast<>, shift<dim, St> > >
incardon's avatar
incardon committed
1295
	CellL getCellList(St r_cut, const Ghost<dim, St> & enlarge, bool no_se3 = false)
1296
	{
incardon's avatar
incardon committed
1297
#ifdef SE_CLASS3
1298 1299
		if (no_se3 == false)
		{se3.getNN();}
incardon's avatar
incardon committed
1300 1301
#endif

1302 1303
		CellL cell_list;

1304 1305
		// Division array
		size_t div[dim];
1306

1307
		// get the processor bounding box
1308
		Box<dim, St> pbox = getDecomposition().getProcessorBounds();
1309

1310
		// Processor bounding box
1311
		cl_param_calculate(pbox, div, r_cut, enlarge);
1312

incardon's avatar
incardon committed
1313 1314
		cell_list.Initialize(pbox, div);
		cell_list.set_gm(g_m);
incardon's avatar
incardon committed
1315
		cell_list.set_ndec(getDecomposition().get_ndec());
1316

1317
		updateCellList(cell_list,no_se3);
1318 1319 1320

		return cell_list;
	}
1321

1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333
	/*! \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
1334 1335
	 * \return The Cell-list
	 *
1336
	 */
1337
	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)
1338
	{
incardon's avatar
incardon committed
1339 1340 1341 1342
#ifdef SE_CLASS3
		se3.getNN();
#endif

1343 1344 1345
		CellL cell_list;

		// Division array
1346 1347
		size_t div[dim];

1348
		// get the processor bounding box
1349
		Box<dim, St> pbox = getDecomposition().getProcessorBounds();
1350

1351
		// Processor bounding box
1352
		cl_param_calculate(pbox,div, r_cut, enlarge);
1353

incardon's avatar
incardon committed
1354 1355
		cell_list.Initialize(pbox, div);
		cell_list.set_gm(g_m);
incardon's avatar
incardon committed
1356
		cell_list.set_ndec(getDecomposition().get_ndec());
1357

Pietro Incardona's avatar
Pietro Incardona committed
1358
		updateCellList(cell_list);
1359 1360 1361

		return cell_list;
	}
incardon's avatar
incardon committed
1362

incardon's avatar
incardon committed
1363 1364 1365 1366
	/*! \brief for each particle get the symmetric verlet list
	 *
	 * \param r_cut cut-off radius
	 *
incardon's avatar
incardon committed
1367 1368
	 * \return the verlet list
	 *
incardon's avatar
incardon committed
1369
	 */
1370 1371
	template <typename VerletL = VerletList<dim,St,Mem_fast<>,shift<dim,St> >>
	VerletL getVerletSym(St r_cut)
incardon's avatar
incardon committed
1372
	{
incardon's avatar
incardon committed
1373 1374 1375 1376
#ifdef SE_CLASS3
		se3.getNN();
#endif

1377
		VerletL ver;
1378

incardon's avatar
incardon committed
1379 1380 1381 1382 1383
		// 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
1384 1385
		ver.set_ndec(getDecomposition().get_ndec());

incardon's avatar
incardon committed
1386 1387
		return ver;
	}
1388

incardon's avatar
incardon committed
1389 1390 1391 1392 1393 1394 1395
	/*! \brief for each particle get the symmetric verlet list
	 *
	 * \param r_cut cut-off radius
	 *
	 * \return the verlet list
	 *
	 */
1396 1397
	template <typename VerletL = VerletList<dim,St,Mem_fast<>,shift<dim,St> >>
	VerletL getVerletCrs(St r_cut)
incardon's avatar
incardon committed
1398
	{
1399
#ifdef SE_CLASS1
incardon's avatar
incardon committed
1400
		if (!(opt & BIND_DEC_TO_GHOST))
1401 1402 1403 1404 1405 1406
		{
			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
1407 1408 1409 1410
#ifdef SE_CLASS3
		se3.getNN();
#endif

1411
		VerletL ver;
incardon's avatar
incardon committed
1412 1413 1414 1415 1416 1417 1418 1419 1420 1421 1422

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

		// Initialize the verlet list
		ver.InitializeCrs(getDecomposition().getDomain(),pbox,getDecomposition().getGhost(),r_cut,v_pos,g_m);

		// Get the internal cell list
		auto & NN = ver.getInternalCellList();

		// Shift
1423
		grid_key_dx<dim> shift;
incardon's avatar
incardon committed
1424 1425 1426

		// Add padding
		for (size_t i = 0 ; i < dim ; i++)
1427
			shift.set_d(i,NN.getPadding(i));
incardon's avatar
incardon committed
1428 1429 1430

		grid_sm<dim,void> gs = NN.getInternalGrid();

1431 1432
		getDecomposition().setNNParameters(shift,gs);

incardon's avatar
incardon committed
1433
		ver.createVerletCrs(r_cut,g_m,v_pos,
1434 1435
				            getDecomposition().getCRSDomainCells(),
							getDecomposition().getCRSAnomDomainCells());
incardon's avatar
incardon committed
1436

incardon's avatar
incardon committed
1437 1438
		ver.set_ndec(getDecomposition().get_ndec());

incardon's avatar
incardon committed
1439 1440 1441
		return ver;
	}

1442 1443 1444 1445
	/*! \brief for each particle get the verlet list
	 *
	 * \param r_cut cut-off radius
	 *
incardon's avatar
incardon committed
1446 1447
	 * \return a VerletList object
	 *
1448
	 */
1449 1450
	template <typename VerletL = VerletList<dim,St,Mem_fast<>,shift<dim,St> >>
	VerletL getVerlet(St r_cut)
1451
	{
incardon's avatar
incardon committed
1452 1453 1454 1455
#ifdef SE_CLASS3
		se3.getNN();
#endif

1456
		VerletL ver;
1457 1458

		// get the processor bounding box
1459 1460 1461 1462 1463
		Box<dim, St> bt = getDecomposition().getProcessorBounds();

		// Get the ghost
		Ghost<dim,St> g = getDecomposition().getGhost();
		g.magnify(1.013);
1464

1465 1466
		// enlarge the box where the Verlet is defined
		bt.enlarge(g);
1467

incardon's avatar
incardon committed
1468
		ver.Initialize(bt,getDecomposition().getProcessorBounds(),r_cut,v_pos,g_m,VL_NON_SYMMETRIC);
1469

incardon's avatar
incardon committed
1470 1471
		ver.set_ndec(getDecomposition().get_ndec());

1472 1473 1474
		return ver;
	}

1475 1476 1477 1478 1479
	/*! \brief for each particle get the verlet list
	 *
	 * \param r_cut cut-off radius
	 * \param ver Verlet to update
	 * \param r_cut cutoff radius
incardon's avatar
incardon committed
1480
	 * \param opt option like VL_SYMMETRIC and VL_NON_SYMMETRIC or VL_CRS_SYMMETRIC
1481 1482
	 *
	 */
1483
	template<typename Mem_type> void updateVerlet(VerletList<dim,St,Mem_type,shift<dim,St> > & ver, St r_cut, size_t opt = VL_NON_SYMMETRIC)
1484
	{
incardon's avatar
incardon committed
1485 1486 1487 1488
#ifdef SE_CLASS3
		se3.getNN();
#endif

incardon's avatar
incardon committed
1489 1490 1491 1492 1493 1494
		if (opt == VL_SYMMETRIC)
		{
			auto & NN = ver.getInternalCellList();

			// Here we have to check that the Box defined by the Cell-list is the same as the domain box of this
			// processor. if it is not like that we have to completely reconstruct from stratch
incardon's avatar
incardon committed
1495
			bool to_reconstruct = NN.get_ndec() != getDecomposition().get_ndec();
incardon's avatar
incardon committed
1496 1497 1498 1499 1500

			if (to_reconstruct == false)
				ver.update(getDecomposition().getDomain(),r_cut,v_pos,g_m, opt);
			else
			{
1501
				VerletList<dim,St,Mem_type,shift<dim,St> > ver_tmp;
incardon's avatar
incardon committed
1502

1503
				ver_tmp = getVerlet<VerletList<dim,St,Mem_type,shift<dim,St> >>(r_cut);
incardon's avatar
incardon committed
1504 1505 1506 1507
				ver.swap(ver);
			}
		}
		else if (opt == VL_CRS_SYMMETRIC)
incardon's avatar
incardon committed
1508
		{
1509 1510 1511 1512 1513 1514 1515 1516
#ifdef SE_CLASS1
			if ((opt & BIND_DEC_TO_GHOST))
			{
				std::cerr << __FILE__ << ":" << __LINE__ << " Error the vector has been constructed without BIND_DEC_TO_GHOST, updateVerlet with the option VL_CRS_SYMMETRIC 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
1517 1518
			auto & NN = ver.getInternalCellList();

incardon's avatar
incardon committed
1519 1520
			// Here we have to check that the Box defined by the Cell-list is the same as the domain box of this
			// processor. if it is not like that we have to completely reconstruct from stratch
incardon's avatar
incardon committed
1521
			bool to_reconstruct = NN.get_ndec() != getDecomposition().get_ndec();
incardon's avatar
incardon committed
1522 1523 1524 1525

			if (to_reconstruct == false)
			{
				// Shift
1526
				grid_key_dx<dim> shift;
incardon's avatar
incardon committed
1527 1528 1529

				// Add padding
				for (size_t i = 0 ; i < dim ; i++)
1530
					shift.set_d(i,NN.getPadding(i));
incardon's avatar
incardon committed
1531 1532 1533

				grid_sm<dim,void> gs = NN.getInternalGrid();

1534 1535
				getDecomposition().setNNParameters(shift,gs);

incardon's avatar
incardon committed
1536
				ver.updateCrs(getDecomposition().getDomain(),r_cut,v_pos,g_m,
1537 1538
						      getDecomposition().getCRSDomainCells(),
							  getDecomposition().getCRSAnomDomainCells());
incardon's avatar
incardon committed
1539 1540 1541
			}
			else
			{
1542
				VerletList<dim,St,Mem_type,shift<dim,St> > ver_tmp;
incardon's avatar
incardon committed
1543

1544
				ver_tmp = getVerletCrs<VerletList<dim,St,Mem_type,shift<dim,St> >>(r_cut);
incardon's avatar
incardon committed
1545 1546 1547 1548 1549 1550 1551 1552 1553
				ver.swap(ver_tmp);
			}
		}
		else
		{
			auto & NN = ver.getInternalCellList();

			// Here we have to check that the Box defined by the Cell-list is the same as the domain box of this
			// processor. if it is not like that we have to completely reconstruct from stratch
incardon's avatar
incardon committed
1554
			bool to_reconstruct = NN.get_ndec() != getDecomposition().get_ndec();
incardon's avatar
incardon committed
1555 1556 1557 1558 1559

			if (to_reconstruct == false)
				ver.update(getDecomposition().getDomain(),r_cut,v_pos,g_m, opt);
			else
			{
1560
				VerletList<dim,St,Mem_type,shift<dim,St> > ver_tmp;
incardon's avatar
incardon committed
1561

1562
				ver_tmp = getVerlet<VerletList<dim,St,Mem_type,shift<dim,St> >>(r_cut);
incardon's avatar
incardon committed
1563 1564 1565
				ver.swap(ver_tmp);
			}
		}
1566 1567
	}

incardon's avatar
incardon committed
1568

1569 1570 1571 1572 1573 1574 1575
	/*! \brief Construct a cell list starting from the stored particles and reorder a vector according to the Hilberts curve
	 *
	 * \tparam CellL CellList type to construct
	 *
	 * \param m an order of a hilbert curve
	 *
	 */
1576 1577
	template<typename CellL=CellList_gen<dim,St,Process_keys_lin,Mem_bal<>,shift<dim,St> > >
	void reorder (int32_t m, reorder_opt opt = reorder_opt::HILBERT)
1578
	{
1579
		reorder<CellL>(m,getDecomposition().getGhost(),opt);
1580 1581 1582 1583 1584 1585 1586 1587 1588 1589 1590 1591 1592 1593 1594
	}


	/*! \brief Construct a cell list starting from the stored particles and reorder a vector according to the Hilberts curve
	 *
	 *
	 *It differs from the reorder(m) 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)
	 *
	 * \param m order of a curve
	 * \param enlarge In case of padding particles the cell list must be enlarged, like a ghost this parameter say how much must be enlarged
	 *
	 */
1595 1596
	template<typename CellL=CellList_gen<dim,St,Process_keys_lin,Mem_bal<>,shift<dim,St> > >
	void reorder(int32_t m, const Ghost<dim,St> & enlarge, reorder_opt opt = reorder_opt::HILBERT)