vector_dist.hpp 46 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_

11
#include "HDF5_XdmfWriter/HDF5_XdmfWriter.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"
19
#include "NN/CellList/CellListFast_hilb.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

incardon's avatar
incardon committed
36 37 38 39 40 41 42 43 44 45 46 47 48
#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
49 50 51
#define NO_ID false
#define ID true

incardon's avatar
incardon committed
52 53
#define DEC_GRAN(gr) ((size_t)gr << 32)

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

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

incardon's avatar
incardon committed
62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99
//! General function t get a cell-list
template<unsigned int dim, typename St, typename CellL, typename Vector>
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)
	{
		return vd.getCellList(r_cut);
	}
};

//! General function t get a cell-list
template<unsigned int dim, typename St, typename Vector>
struct gcl<dim,St,CellList_hilb<dim, St, FAST, shift<dim, St> >,Vector>
{
	/*! \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 CellList_hilb<dim, St, FAST, shift<dim, St> > get(Vector & vd, const St & r_cut, const Ghost<dim,St> & g)
	{
		return vd.getCellList_hilb(r_cut,g);
	}
};

incardon's avatar
incardon committed
100
/*! \brief Distributed vector
101
 *
102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121
 * 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
122 123 124
 *
 */

125
template<unsigned int dim, typename St, typename prop, typename Decomposition = CartDecomposition<dim,St>, typename Memory = HeapMemory>
Pietro Incardona's avatar
Pietro Incardona committed
126
class vector_dist : public vector_dist_comm<dim,St,prop,Decomposition,Memory>
incardon's avatar
incardon committed
127
{
incardon's avatar
incardon committed
128 129 130 131 132 133 134 135
public:

	//! Self type
	typedef vector_dist<dim,St,prop,Decomposition,Memory> self;

	//! property object
	typedef prop value_type;

incardon's avatar
incardon committed
136 137
private:

138
	//! Ghost marker, all the particle with id > g_m are ghost all with g_m < are real particle
139 140
	size_t g_m = 0;

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

145 146
	//! Particle properties vector, (It has 2 elements) the first has real particles assigned to a processor
	//! the second element contain unassigned particles
147
	openfpm::vector<prop> v_prp;
incardon's avatar
incardon committed
148

149
	//! Virtual cluster
incardon's avatar
incardon committed
150 151
	Vcluster & v_cl;

incardon's avatar
incardon committed
152 153 154 155 156 157 158 159
	//! option used to create this vector
	size_t opt = 0;

#ifdef SE_CLASS3

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

#endif
160

Pietro Incardona's avatar
Pietro Incardona committed
161
	/*! \brief Initialize the structures
incardon's avatar
incardon committed
162
	 *
Pietro Incardona's avatar
Pietro Incardona committed
163
	 * \param np number of particles
incardon's avatar
incardon committed
164 165
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
166
	void init_structures(size_t np)
incardon's avatar
incardon committed
167
	{
168
		// convert to a local number of elements
169 170 171 172 173 174 175 176
		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++;
177

incardon's avatar
incardon committed
178
		// resize the position vector
179
		v_pos.resize(p_np);
incardon's avatar
incardon committed
180 181

		// resize the properties vector
182 183 184
		v_prp.resize(p_np);

		g_m = p_np;
Pietro Incardona's avatar
Pietro Incardona committed
185
	}
incardon's avatar
incardon committed
186

incardon's avatar
incardon committed
187
	/*! \brief Check if the parameters describe a valid vector. In case it does not report an error
incardon's avatar
incardon committed
188 189 190 191 192 193 194 195
	 *
	 * \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
196
		{
incardon's avatar
incardon committed
197
			std::cerr << __FILE__ << ":" << __LINE__ << "  Error the domain is not valid " << box.toString() << std::endl;
incardon's avatar
incardon committed
198 199
			ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
		}
incardon's avatar
incardon committed
200 201 202

	}

incardon's avatar
incardon committed
203 204 205 206 207
	/*! \brief It check that the r_cut is not bugger than the ghost
	 *
	 * \param r_cut cut-off radius
	 *
	 */
incardon's avatar
incardon committed
208
	void check_ghost_compatible_rcut(St r_cut)
incardon's avatar
incardon committed
209 210 211
	{
		for (size_t i = 0 ; i < dim ; i++)
		{
incardon's avatar
incardon committed
212 213 214 215 216
			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
217 218 219
		}
	}

Pietro Incardona's avatar
Pietro Incardona committed
220 221 222 223 224 225 226 227
public:

	//! space type
	typedef St stype;

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

Pietro Incardona's avatar
Pietro Incardona committed
228 229 230 231
	/*! \brief Operator= for distributed vector
	 *
	 * \param v vector to copy
	 *
Pietro Incardona's avatar
Pietro Incardona committed
232 233
	 * \return itself
	 *
Pietro Incardona's avatar
Pietro Incardona committed
234 235 236
	 */
	vector_dist<dim,St,prop,Decomposition,Memory> & operator=(const vector_dist<dim,St,prop,Decomposition,Memory> & v)
	{
Pietro Incardona's avatar
Pietro Incardona committed
237
		static_cast<vector_dist_comm<dim,St,prop,Decomposition,Memory> *>(this)->operator=(static_cast<vector_dist_comm<dim,St,prop,Decomposition,Memory>>(v));
Pietro Incardona's avatar
Pietro Incardona committed
238 239 240 241 242

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

incardon's avatar
incardon committed
243 244 245 246 247 248
#ifdef SE_CLASS3
		se3 = v.se3;
#endif

		opt = v.opt;

Pietro Incardona's avatar
Pietro Incardona committed
249 250 251 252 253 254 255
		return *this;
	}

	/*! \brief Operator= for distributed vector
	 *
	 * \param v vector to copy
	 *
Pietro Incardona's avatar
Pietro Incardona committed
256 257
	 * \return itself
	 *
Pietro Incardona's avatar
Pietro Incardona committed
258 259 260
	 */
	vector_dist<dim,St,prop,Decomposition,Memory> & operator=(vector_dist<dim,St,prop,Decomposition,Memory> && v)
	{
Pietro Incardona's avatar
Pietro Incardona committed
261
		static_cast<vector_dist_comm<dim,St,prop,Decomposition,Memory> *>(this)->operator=(static_cast<vector_dist_comm<dim,St,prop,Decomposition,Memory> >(v));
Pietro Incardona's avatar
Pietro Incardona committed
262 263 264 265 266

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

incardon's avatar
incardon committed
267 268 269 270 271 272
#ifdef SE_CLASS3
		se3 = v.se3;
#endif

		opt = v.opt;

Pietro Incardona's avatar
Pietro Incardona committed
273 274 275
		return *this;
	}

incardon's avatar
incardon committed
276

Pietro Incardona's avatar
Pietro Incardona committed
277
	/*! \brief Copy Constructor
Pietro Incardona's avatar
Pietro Incardona committed
278
	 *
Pietro Incardona's avatar
Pietro Incardona committed
279
	 * \param v vector to copy
Pietro Incardona's avatar
Pietro Incardona committed
280 281 282
	 *
	 */
	vector_dist(const vector_dist<dim,St,prop,Decomposition,Memory> & v)
incardon's avatar
incardon committed
283
	:vector_dist_comm<dim,St,prop,Decomposition,Memory>(v.getDecomposition()),v_cl(v.v_cl) SE_CLASS3_VDIST_CONSTRUCTOR
Pietro Incardona's avatar
Pietro Incardona committed
284 285 286 287 288 289 290 291
	{
#ifdef SE_CLASS2
		check_new(this,8,VECTOR_DIST_EVENT,4);
#endif

		this->operator=(v);
	}

Pietro Incardona's avatar
Pietro Incardona committed
292
	/*! \brief Copy constructor
Pietro Incardona's avatar
Pietro Incardona committed
293
	 *
Pietro Incardona's avatar
Pietro Incardona committed
294
	 * \param v vector to copy
Pietro Incardona's avatar
Pietro Incardona committed
295 296
	 *
	 */
incardon's avatar
incardon committed
297
	vector_dist(vector_dist<dim,St,prop,Decomposition,Memory> && v) noexcept
incardon's avatar
incardon committed
298
	:v_cl(v.v_cl) SE_CLASS3_VDIST_CONSTRUCTOR
Pietro Incardona's avatar
Pietro Incardona committed
299 300 301 302 303 304
	{
#ifdef SE_CLASS2
		check_new(this,8,VECTOR_DIST_EVENT,4);
#endif

		this->operator=(v);
incardon's avatar
incardon committed
305 306 307 308

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

Pietro Incardona's avatar
Pietro Incardona committed
311
	/*! \brief Constructor with predefined decomposition
Pietro Incardona's avatar
Pietro Incardona committed
312
	 *
Pietro Incardona's avatar
Pietro Incardona committed
313 314
	 * \param dec is the decomposition
	 * \param np number of particles
Pietro Incardona's avatar
Pietro Incardona committed
315 316 317
	 *
	 */
	vector_dist(const Decomposition & dec, size_t np) :
incardon's avatar
incardon committed
318
	vector_dist_comm<dim,St,prop,Decomposition,Memory>(dec), v_cl(create_vcluster()) SE_CLASS3_VDIST_CONSTRUCTOR
Pietro Incardona's avatar
Pietro Incardona committed
319 320 321 322 323 324
	{
#ifdef SE_CLASS2
		check_new(this,8,VECTOR_DIST_EVENT,4);
#endif

		init_structures(np);
incardon's avatar
incardon committed
325 326 327 328

#ifdef SE_CLASS3
		se3.Initialize();
#endif
Pietro Incardona's avatar
Pietro Incardona committed
329 330 331
	}


incardon's avatar
incardon committed
332
	/*! \brief Constructor of a distributed vector
Pietro Incardona's avatar
Pietro Incardona committed
333 334 335
	 *
	 * \param np number of elements
	 * \param box domain where the vector of elements live
Pietro Incardona's avatar
Pietro Incardona committed
336
	 * \param bc boundary conditions
Pietro Incardona's avatar
Pietro Incardona committed
337
	 * \param g Ghost margins
incardon's avatar
incardon committed
338
	 * \param opt additional options. BIND_DEC_TO_GHOST Bind the decomposition to be multiple of the
339 340
	 *          ghost size. This is required if we want to use symmetric to eliminate
	 *          ghost communications.
Pietro Incardona's avatar
Pietro Incardona committed
341 342
	 *
	 */
343
	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
344
	:v_cl(create_vcluster()),opt(opt) SE_CLASS3_VDIST_CONSTRUCTOR
Pietro Incardona's avatar
Pietro Incardona committed
345 346 347 348 349
	{
#ifdef SE_CLASS2
		check_new(this,8,VECTOR_DIST_EVENT,4);
#endif

incardon's avatar
incardon committed
350 351 352
		if (opt >> 32 != 0)
			this->setDecompositionGranularity(opt >> 32);

incardon's avatar
incardon committed
353 354
		check_parameters(box);

Pietro Incardona's avatar
Pietro Incardona committed
355
		init_structures(np);
356
		this->init_decomposition(box,bc,g,opt,gdist);
incardon's avatar
incardon committed
357 358 359 360

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

363 364 365 366 367 368 369
	~vector_dist()
	{
#ifdef SE_CLASS2
		check_delete(this);
#endif
	}

Pietro Incardona's avatar
Pietro Incardona committed
370
	/*! \brief return the local size of the vector
371
	 *
Pietro Incardona's avatar
Pietro Incardona committed
372
	 * \return local size
373 374
	 *
	 */
incardon's avatar
incardon committed
375
	size_t size_local() const
376
	{
Pietro Incardona's avatar
Pietro Incardona committed
377
		return g_m;
378 379
	}

incardon's avatar
incardon committed
380 381 382 383 384
	/*! \brief return the local size of the vector
	 *
	 * \return local size
	 *
	 */
incardon's avatar
incardon committed
385
	size_t size_local_with_ghost() const
incardon's avatar
incardon committed
386
	{
Pietro Incardona's avatar
Pietro Incardona committed
387
		return v_pos.size();
incardon's avatar
incardon committed
388 389
	}

incardon's avatar
incardon committed
390 391
#ifndef ONLY_READWRITE_GETTER

392
	/*! \brief Get the position of an element
incardon's avatar
incardon committed
393
	 *
394 395 396 397 398
	 * 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
399 400
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
401
	inline auto getPos(vect_dist_key_dx vec_key) -> decltype(v_pos.template get<0>(vec_key.getKey()))
incardon's avatar
incardon committed
402
	{
Pietro Incardona's avatar
Pietro Incardona committed
403
		return v_pos.template get<0>(vec_key.getKey());
incardon's avatar
incardon committed
404 405
	}

406 407 408 409 410 411 412 413 414 415 416 417 418 419
	/*! \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()))
	{
		return v_pos.template get<0>(vec_key.getKey());
	}

420 421 422 423 424 425 426 427 428
	/*! \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
429 430 431 432 433 434 435 436 437 438 439 440 441 442
	inline auto getPos(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
	 *
	 */
incardon's avatar
incardon committed
443
	inline auto getPos(size_t vec_key) const -> decltype(v_pos.template get<0>(vec_key))
444 445 446 447
	{
		return v_pos.template get<0>(vec_key);
	}

448
	/*! \brief Get the property of an element
incardon's avatar
incardon committed
449
	 *
450 451 452
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \tparam id property id
incardon's avatar
incardon committed
453 454
	 * \param vec_key vector element
	 *
455 456
	 * \return return the selected property of the vector element
	 *
incardon's avatar
incardon committed
457
	 */
458
	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
459
	{
460
		return v_prp.template get<id>(vec_key.getKey());
incardon's avatar
incardon committed
461 462
	}

Pietro Incardona's avatar
Pietro Incardona committed
463 464 465 466 467 468 469 470 471 472
	/*! \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
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
473
	template<unsigned int id> inline auto getProp(vect_dist_key_dx vec_key) const -> const decltype(v_prp.template get<id>(vec_key.getKey()))
Pietro Incardona's avatar
Pietro Incardona committed
474 475 476 477
	{
		return v_prp.template get<id>(vec_key.getKey());
	}

incardon's avatar
incardon committed
478 479 480 481 482 483 484 485 486 487 488 489 490 491 492
	/*! \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))
	{
		return v_prp.template get<id>(vec_key);
	}

493 494 495 496 497 498 499 500 501 502 503 504 505 506 507
	/*! \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) const -> const decltype(v_prp.template get<id>(vec_key))
	{
		return v_prp.template get<id>(vec_key);
	}

incardon's avatar
incardon committed
508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587
#endif

///////////////////// 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
	 *
	 */
	template<unsigned int id> inline auto getPropRead(vect_dist_key_dx vec_key) const -> const decltype(v_prp.template get<id>(vec_key.getKey()))
	{
#ifdef SE_CLASS3
		se3.template read<id>(*this,vec_key.getKey());
#endif

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

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

588 589 590 591 592 593 594 595 596 597 598 599 600 601
	/*! \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
602 603 604 605 606

#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
607 608 609 610 611 612 613
	}

	/*! \brief Get the position of the last element
	 *
	 * \return the position of the element in space
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
614
	inline auto getLastPos() -> decltype(v_pos.template get<0>(0))
615
	{
Pietro Incardona's avatar
Pietro Incardona committed
616
		return v_pos.template get<0>(g_m - 1);
617 618 619 620 621 622 623 624 625 626 627
	}

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

incardon's avatar
incardon committed
631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 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
////////////////////////////// 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
		se3.read<prop::max_prop_real>(*this,g_m-1);
#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
		se3.write<prop::max_prop_real>(*this,g_m-1);
#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
		se3.write<id>(*this,g_m-1);
#endif

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

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

incardon's avatar
incardon committed
696 697 698 699 700 701 702 703 704 705 706
	/*! \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, FAST, shift<dim, St> > > CellL getCellListSym(St r_cut)
	{
incardon's avatar
incardon committed
707 708 709 710 711 712 713 714
#ifdef SE_CLASS1
		if ((opt & BIND_DEC_TO_GHOST))
		{
			std::cerr << __FILE__ << ":" << __LINE__ << " error to get symmetric cell-list you must construct the vector with the option BIND_DEC_TO_GHOST " << std::endl;
			ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
		}
#endif

incardon's avatar
incardon committed
715 716 717 718
		// Cell list
		CellL cell_list;

		size_t pad = 0;
incardon's avatar
incardon committed
719 720
		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
721

incardon's avatar
incardon committed
722
		// Processor bounding box
incardon's avatar
incardon committed
723 724
		Box<dim, St> pbox = getDecomposition().getProcessorBounds();

incardon's avatar
incardon committed
725 726 727
		// Ghost padding extension
		Ghost<dim,size_t> g_ext(0);
		cell_list.Initialize(cd_sm,pbox,pad);
728
		cell_list.set_ndec(getDecomposition().get_ndec());
incardon's avatar
incardon committed
729

incardon's avatar
incardon committed
730
		updateCellListSym(cell_list);
incardon's avatar
incardon committed
731 732

		return cell_list;
incardon's avatar
incardon committed
733 734
	}

735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750

	/*! \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, FAST, shift<dim, St> > > CellL getCellListSymNoBind(St r_cut)
	{
		return getCellList(r_cut);
	}


Pietro Incardona's avatar
Pietro Incardona committed
751 752 753 754
	/*! \brief Construct a cell list starting from the stored particles
	 *
	 * \tparam CellL CellList type to construct
	 *
755 756 757 758
	 * \param r_cut interation radius, or size of each cell
	 *
	 * \return the Cell list
	 *
Pietro Incardona's avatar
Pietro Incardona committed
759
	 */
760
	template<typename CellL = CellList<dim, St, FAST, shift<dim, St> > > CellL getCellList(St r_cut)
761
	{
incardon's avatar
incardon committed
762 763 764 765 766 767 768
#ifdef SE_CLASS3
		se3.getNN();
#endif
#ifdef SE_CLASS1
		check_ghost_compatible_rcut(r_cut);
#endif

Pietro Incardona's avatar
Pietro Incardona committed
769
		// Get ghost and anlarge by 1%
Pietro Incardona's avatar
Pietro Incardona committed
770
		Ghost<dim,St> g = getDecomposition().getGhost();
771
		g.magnify(1.013);
Pietro Incardona's avatar
Pietro Incardona committed
772 773

		return getCellList(r_cut, g);
774 775
	}

776 777 778 779 780 781 782 783 784 785 786
	/*! \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
	 *
	 */
	template<typename CellL = CellList_hilb<dim, St, FAST, shift<dim, St> > > CellL getCellList_hilb(St r_cut)
	{
incardon's avatar
incardon committed
787 788 789 790 791 792 793
#ifdef SE_CLASS3
		se3.getNN();
#endif
#ifdef SE_CLASS1
		check_ghost_compatible_rcut(r_cut);
#endif

794
		// Get ghost and anlarge by 1%
Pietro Incardona's avatar
Pietro Incardona committed
795
		Ghost<dim,St> g = getDecomposition().getGhost();
Pietro Incardona's avatar
Pietro Incardona committed
796
		g.magnify(1.013);
797 798 799 800

		return getCellList_hilb(r_cut, g);
	}

Pietro Incardona's avatar
Pietro Incardona committed
801 802 803 804 805 806 807 808 809
	/*! \brief Update a cell list using the stored particles
	 *
	 * \tparam CellL CellList type to construct
	 *
	 * \param cell_list Cell list to update
	 *
	 */
	template<typename CellL = CellList<dim, St, FAST, shift<dim, St> > > void updateCellList(CellL & cell_list)
	{
incardon's avatar
incardon committed
810 811 812 813
#ifdef SE_CLASS3
		se3.getNN();
#endif

incardon's avatar
incardon committed
814 815 816 817 818 819
		// 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
820 821 822
		// 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
823 824 825 826 827 828 829 830 831 832 833 834 835

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

			cell_list.set_gm(g_m);
		}
		else
		{
			CellL cli_tmp = gcl<dim,St,CellL,self>::get(*this,r_cut,getDecomposition().getGhost());

			cell_list.swap(cli_tmp);
		}
incardon's avatar
incardon committed
836 837 838 839 840 841 842 843 844 845 846
	}

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

incardon's avatar
incardon committed
851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872
		// 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));

		// 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
		{
			CellL cli_tmp = gcl<dim,St,CellL,self>::get(*this,r_cut,getDecomposition().getGhost());

			cell_list.swap(cli_tmp);
		}
Pietro Incardona's avatar
Pietro Incardona committed
873 874
	}

875 876 877 878 879 880 881 882 883 884 885 886
	/*! \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
	 *
Pietro Incardona's avatar
Pietro Incardona committed
887 888
	 * \return the CellList
	 *
889
	 */
890
	template<typename CellL = CellList<dim, St, FAST, shift<dim, St> > > CellL getCellList(St r_cut, const Ghost<dim, St> & enlarge)
Pietro Incardona's avatar
Pietro Incardona committed
891
	{
incardon's avatar
incardon committed
892 893 894 895
#ifdef SE_CLASS3
		se3.getNN();
#endif

Pietro Incardona's avatar
Pietro Incardona committed
896 897
		CellL cell_list;

898 899
		// Division array
		size_t div[dim];
900

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

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

907
		cell_list.Initialize(pbox, div);
incardon's avatar
incardon committed
908
		cell_list.set_ndec(getDecomposition().get_ndec());
909 910 911 912 913

		updateCellList(cell_list);

		return cell_list;
	}
914

915 916 917 918 919 920 921 922 923 924 925 926
	/*! \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
927 928
	 * \return The Cell-list
	 *
929 930 931
	 */
	template<typename CellL = CellList_hilb<dim, St, FAST, shift<dim, St> > > CellL getCellList_hilb(St r_cut, const Ghost<dim, St> & enlarge)
	{
incardon's avatar
incardon committed
932 933 934 935
#ifdef SE_CLASS3
		se3.getNN();
#endif

936 937 938
		CellL cell_list;

		// Division array
939 940
		size_t div[dim];

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

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

947
		cell_list.Initialize(pbox, div, g_m);
incardon's avatar
incardon committed
948
		cell_list.set_ndec(getDecomposition().get_ndec());
949

Pietro Incardona's avatar
Pietro Incardona committed
950
		updateCellList(cell_list);
Pietro Incardona's avatar
Pietro Incardona committed
951 952 953

		return cell_list;
	}
incardon's avatar
incardon committed
954

incardon's avatar
incardon committed
955 956 957 958
	/*! \brief for each particle get the symmetric verlet list
	 *
	 * \param r_cut cut-off radius
	 *
incardon's avatar
incardon committed
959 960
	 * \return the verlet list
	 *
incardon's avatar
incardon committed
961 962 963
	 */
	VerletList<dim,St,FAST,shift<dim,St> > getVerletSym(St r_cut)
	{
incardon's avatar
incardon committed
964 965 966 967
#ifdef SE_CLASS3
		se3.getNN();
#endif

incardon's avatar
incardon committed
968
		VerletList<dim,St,FAST,shift<dim,St>> ver;
969

incardon's avatar
incardon committed
970 971 972 973 974
		// 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
975 976
		ver.set_ndec(getDecomposition().get_ndec());

incardon's avatar
incardon committed
977 978
		return ver;
	}
Pietro Incardona's avatar
Pietro Incardona committed
979

incardon's avatar
incardon committed
980 981 982 983 984 985 986 987 988
	/*! \brief for each particle get the symmetric verlet list
	 *
	 * \param r_cut cut-off radius
	 *
	 * \return the verlet list
	 *
	 */
	VerletList<dim,St,FAST,shift<dim,St> > getVerletCrs(St r_cut)
	{
incardon's avatar
incardon committed
989 990 991 992
#ifdef SE_CLASS3
		se3.getNN();
#endif

incardon's avatar
incardon committed
993 994 995 996 997 998 999 1000 1001 1002 1003 1004
		VerletList<dim,St,FAST,shift<dim,St>> ver;

		// 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
1005
		grid_key_dx<dim> shift;
incardon's avatar
incardon committed
1006 1007 1008

		// Add padding
		for (size_t i = 0 ; i < dim ; i++)
1009
			shift.set_d(i,NN.getPadding(i));
incardon's avatar
incardon committed
1010 1011 1012

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

1013 1014
		getDecomposition().setNNParameters(shift,gs);

incardon's avatar
incardon committed
1015
		ver.createVerletCrs(r_cut,g_m,v_pos,
1016 1017
				            getDecomposition().getCRSDomainCells(),
							getDecomposition().getCRSAnomDomainCells());
incardon's avatar
incardon committed
1018

incardon's avatar
incardon committed
1019 1020
		ver.set_ndec(getDecomposition().get_ndec());

incardon's avatar
incardon committed
1021 1022 1023
		return ver;
	}

1024 1025 1026 1027
	/*! \brief for each particle get the verlet list
	 *
	 * \param r_cut cut-off radius
	 *
incardon's avatar
incardon committed
1028 1029
	 * \return a VerletList object
	 *
1030
	 */
incardon's avatar
incardon committed
1031
	VerletList<dim,St,FAST,shift<dim,St> > getVerlet(St r_cut)
1032
	{
incardon's avatar
incardon committed
1033 1034 1035 1036
#ifdef SE_CLASS3
		se3.getNN();
#endif

1037 1038 1039
		VerletList<dim,St,FAST,shift<dim,St>> ver;

		// get the processor bounding box
Pietro Incardona's avatar
Pietro Incardona committed
1040 1041 1042 1043 1044
		Box<dim, St> bt = getDecomposition().getProcessorBounds();

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

Pietro Incardona's avatar
Pietro Incardona committed
1046 1047
		// enlarge the box where the Verlet is defined
		bt.enlarge(g);
1048

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

incardon's avatar
incardon committed
1051 1052
		ver.set_ndec(getDecomposition().get_ndec());

1053 1054 1055
		return ver;
	}

Pietro Incardona's avatar
Pietro Incardona committed
1056 1057 1058 1059 1060
	/*! \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
1061
	 * \param opt option like VL_SYMMETRIC and VL_NON_SYMMETRIC or VL_CRS_SYMMETRIC
Pietro Incardona's avatar
Pietro Incardona committed
1062 1063 1064 1065
	 *
	 */
	void updateVerlet(VerletList<dim,St,FAST,shift<dim,St> > & ver, St r_cut, size_t opt = VL_NON_SYMMETRIC)
	{
incardon's avatar
incardon committed
1066 1067 1068 1069
#ifdef SE_CLASS3
		se3.getNN();
#endif

incardon's avatar
incardon committed
1070 1071 1072 1073 1074 1075
		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
1076
			bool to_reconstruct = NN.get_ndec() != getDecomposition().get_ndec();
incardon's avatar
incardon committed
1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088

			if (to_reconstruct == false)
				ver.update(getDecomposition().getDomain(),r_cut,v_pos,g_m, opt);
			else
			{
				VerletList<dim,St,FAST,shift<dim,St> > ver_tmp;

				ver_tmp = getVerlet(r_cut);
				ver.swap(ver);
			}
		}
		else if (opt == VL_CRS_SYMMETRIC)
incardon's avatar
incardon committed
1089 1090 1091
		{
			auto & NN = ver.getInternalCellList();

incardon's avatar
incardon committed
1092 1093
			// 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
1094
			bool to_reconstruct = NN.get_ndec() != getDecomposition().get_ndec();
incardon's avatar
incardon committed
1095 1096 1097 1098

			if (to_reconstruct == false)
			{
				// Shift
1099
				grid_key_dx<dim> shift;
incardon's avatar
incardon committed
1100 1101 1102

				// Add padding
				for (size_t i = 0 ; i < dim ; i++)
1103
					shift.set_d(i,NN.getPadding(i));
incardon's avatar
incardon committed
1104 1105 1106

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

1107 1108
				getDecomposition().setNNParameters(shift,gs);

incardon's avatar
incardon committed
1109
				ver.updateCrs(getDecomposition().getDomain(),r_cut,v_pos,g_m,
1110 1111
						      getDecomposition().getCRSDomainCells(),
							  getDecomposition().getCRSAnomDomainCells());
incardon's avatar
incardon committed
1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126
			}
			else
			{
				VerletList<dim,St,FAST,shift<dim,St> > ver_tmp;

				ver_tmp = getVerletCrs(r_cut);
				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
1127
			bool to_reconstruct = NN.get_ndec() != getDecomposition().get_ndec();
incardon's avatar
incardon committed
1128 1129 1130 1131 1132 1133

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

incardon's avatar
incardon committed