vector_dist.hpp 43.6 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

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

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

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

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

incardon's avatar
incardon committed
61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98
//! 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
99
/*! \brief Distributed vector
100
 *
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
 * 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
121 122 123
 *
 */

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

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

	//! property object
	typedef prop value_type;

incardon's avatar
incardon committed
135 136
private:

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

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

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

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

incardon's avatar
incardon committed
151 152 153 154 155 156 157 158
	//! 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
159

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

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

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

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

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

	}

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

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

	//! space type
	typedef St stype;

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

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

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

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

		opt = v.opt;

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

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

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

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

		opt = v.opt;

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

incardon's avatar
incardon committed
275

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

		this->operator=(v);
	}

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

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

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

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

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

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


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

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

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

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

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

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

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

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

incardon's avatar
incardon committed
389 390
#ifndef ONLY_READWRITE_GETTER

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

405 406 407 408 409 410 411 412 413 414 415 416 417 418
	/*! \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());
	}

419
	/*! \brief Get the property of an element
incardon's avatar
incardon committed
420
	 *
421 422 423
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \tparam id property id
incardon's avatar
incardon committed
424 425
	 * \param vec_key vector element
	 *
426 427
	 * \return return the selected property of the vector element
	 *
incardon's avatar
incardon committed
428
	 */
429
	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
430
	{
431
		return v_prp.template get<id>(vec_key.getKey());
incardon's avatar
incardon committed
432 433
	}

Pietro Incardona's avatar
Pietro Incardona committed
434 435 436 437 438 439 440 441 442 443
	/*! \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
444
	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
445 446 447 448
	{
		return v_prp.template get<id>(vec_key.getKey());
	}

incardon's avatar
incardon committed
449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528
#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());
	}

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

529 530 531 532 533 534 535 536 537 538 539 540 541 542
	/*! \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
543 544 545 546 547

#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
548 549 550 551 552 553 554
	}

	/*! \brief Get the position of the last element
	 *
	 * \return the position of the element in space
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
555
	inline auto getLastPos() -> decltype(v_pos.template get<0>(0))
556
	{
Pietro Incardona's avatar
Pietro Incardona committed
557
		return v_pos.template get<0>(g_m - 1);
558 559 560 561 562 563 564 565 566 567 568
	}

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

incardon's avatar
incardon committed
572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636
////////////////////////////// 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
637 638 639 640 641 642 643 644 645 646 647
	/*! \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
648 649 650 651 652 653 654 655
#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
656 657 658 659
		// Cell list
		CellL cell_list;

		size_t pad = 0;
incardon's avatar
incardon committed
660 661
		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
662

incardon's avatar
incardon committed
663
		// Processor bounding box
incardon's avatar
incardon committed
664 665
		Box<dim, St> pbox = getDecomposition().getProcessorBounds();

incardon's avatar
incardon committed
666 667 668
		// Ghost padding extension
		Ghost<dim,size_t> g_ext(0);
		cell_list.Initialize(cd_sm,pbox,pad);
669
		cell_list.set_ndec(getDecomposition().get_ndec());
incardon's avatar
incardon committed
670

incardon's avatar
incardon committed
671
		updateCellListSym(cell_list);
incardon's avatar
incardon committed
672 673

		return cell_list;
incardon's avatar
incardon committed
674 675
	}

Pietro Incardona's avatar
Pietro Incardona committed
676 677 678 679
	/*! \brief Construct a cell list starting from the stored particles
	 *
	 * \tparam CellL CellList type to construct
	 *
680 681 682 683
	 * \param r_cut interation radius, or size of each cell
	 *
	 * \return the Cell list
	 *
Pietro Incardona's avatar
Pietro Incardona committed
684
	 */
685
	template<typename CellL = CellList<dim, St, FAST, shift<dim, St> > > CellL getCellList(St r_cut)
686
	{
incardon's avatar
incardon committed
687 688 689 690 691 692 693
#ifdef SE_CLASS3
		se3.getNN();
#endif
#ifdef SE_CLASS1
		check_ghost_compatible_rcut(r_cut);
#endif

Pietro Incardona's avatar
Pietro Incardona committed
694
		// Get ghost and anlarge by 1%
Pietro Incardona's avatar
Pietro Incardona committed
695
		Ghost<dim,St> g = getDecomposition().getGhost();
696
		g.magnify(1.013);
Pietro Incardona's avatar
Pietro Incardona committed
697 698

		return getCellList(r_cut, g);
699 700
	}

701 702 703 704 705 706 707 708 709 710 711
	/*! \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
712 713 714 715 716 717 718
#ifdef SE_CLASS3
		se3.getNN();
#endif
#ifdef SE_CLASS1
		check_ghost_compatible_rcut(r_cut);
#endif

719
		// Get ghost and anlarge by 1%
Pietro Incardona's avatar
Pietro Incardona committed
720
		Ghost<dim,St> g = getDecomposition().getGhost();
Pietro Incardona's avatar
Pietro Incardona committed
721
		g.magnify(1.013);
722 723 724 725

		return getCellList_hilb(r_cut, g);
	}

Pietro Incardona's avatar
Pietro Incardona committed
726 727 728 729 730 731 732 733 734
	/*! \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
735 736 737 738
#ifdef SE_CLASS3
		se3.getNN();
#endif

incardon's avatar
incardon committed
739 740 741 742 743 744
		// 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
745 746 747
		// 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
748 749 750 751 752 753 754 755 756 757 758 759 760

		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
761 762 763 764 765 766 767 768 769 770 771
	}

	/*! \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
772 773 774
#ifdef SE_CLASS3
		se3.getNN();
#endif
incardon's avatar
incardon committed
775

incardon's avatar
incardon committed
776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797
		// 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
798 799
	}

800 801 802 803 804 805 806 807 808 809 810 811
	/*! \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
812 813
	 * \return the CellList
	 *
814
	 */
815
	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
816
	{
incardon's avatar
incardon committed
817 818 819 820
#ifdef SE_CLASS3
		se3.getNN();
#endif

Pietro Incardona's avatar
Pietro Incardona committed
821 822
		CellL cell_list;

823 824
		// Division array
		size_t div[dim];
825

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

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

832
		cell_list.Initialize(pbox, div);
incardon's avatar
incardon committed
833
		cell_list.set_ndec(getDecomposition().get_ndec());
834 835 836 837 838

		updateCellList(cell_list);

		return cell_list;
	}
839

840 841 842 843 844 845 846 847 848 849 850 851
	/*! \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
852 853
	 * \return The Cell-list
	 *
854 855 856
	 */
	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
857 858 859 860
#ifdef SE_CLASS3
		se3.getNN();
#endif

861 862 863
		CellL cell_list;

		// Division array
864 865
		size_t div[dim];

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

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

872
		cell_list.Initialize(pbox, div, g_m);
incardon's avatar
incardon committed
873
		cell_list.set_ndec(getDecomposition().get_ndec());
874

Pietro Incardona's avatar
Pietro Incardona committed
875
		updateCellList(cell_list);
Pietro Incardona's avatar
Pietro Incardona committed
876 877 878

		return cell_list;
	}
incardon's avatar
incardon committed
879

incardon's avatar
incardon committed
880 881 882 883
	/*! \brief for each particle get the symmetric verlet list
	 *
	 * \param r_cut cut-off radius
	 *
incardon's avatar
incardon committed
884 885
	 * \return the verlet list
	 *
incardon's avatar
incardon committed
886 887 888
	 */
	VerletList<dim,St,FAST,shift<dim,St> > getVerletSym(St r_cut)
	{
incardon's avatar
incardon committed
889 890 891 892
#ifdef SE_CLASS3
		se3.getNN();
#endif

incardon's avatar
incardon committed
893
		VerletList<dim,St,FAST,shift<dim,St>> ver;
894

incardon's avatar
incardon committed
895 896 897 898 899
		// 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
900 901
		ver.set_ndec(getDecomposition().get_ndec());

incardon's avatar
incardon committed
902 903
		return ver;
	}
Pietro Incardona's avatar
Pietro Incardona committed
904

incardon's avatar
incardon committed
905 906 907 908 909 910 911 912 913
	/*! \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
914 915 916 917
#ifdef SE_CLASS3
		se3.getNN();
#endif

incardon's avatar
incardon committed
918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940
		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
		grid_key_dx<dim> cell_shift = NN.getShift();

		// Shift
		grid_key_dx<dim> shift = NN.getShift();

		// Add padding
		for (size_t i = 0 ; i < dim ; i++)
			shift.set_d(i,shift.get(i) + NN.getPadding(i));

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

incardon's avatar
incardon committed
941 942 943
		ver.createVerletCrs(r_cut,g_m,v_pos,
				            getDecomposition().getCRSDomainCells(shift,cell_shift,gs),
							getDecomposition().getCRSAnomDomainCells(shift,cell_shift,gs));
incardon's avatar
incardon committed
944

incardon's avatar
incardon committed
945 946
		ver.set_ndec(getDecomposition().get_ndec());

incardon's avatar
incardon committed
947 948 949
		return ver;
	}

950 951 952 953
	/*! \brief for each particle get the verlet list
	 *
	 * \param r_cut cut-off radius
	 *
incardon's avatar
incardon committed
954 955
	 * \return a VerletList object
	 *
956
	 */
incardon's avatar
incardon committed
957
	VerletList<dim,St,FAST,shift<dim,St> > getVerlet(St r_cut)
958
	{
incardon's avatar
incardon committed
959 960 961 962
#ifdef SE_CLASS3
		se3.getNN();
#endif

963 964 965
		VerletList<dim,St,FAST,shift<dim,St>> ver;

		// get the processor bounding box
Pietro Incardona's avatar
Pietro Incardona committed
966 967 968 969 970
		Box<dim, St> bt = getDecomposition().getProcessorBounds();

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

Pietro Incardona's avatar
Pietro Incardona committed
972 973
		// enlarge the box where the Verlet is defined
		bt.enlarge(g);
974

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

incardon's avatar
incardon committed
977 978
		ver.set_ndec(getDecomposition().get_ndec());

979 980 981
		return ver;
	}

Pietro Incardona's avatar
Pietro Incardona committed
982 983 984 985 986
	/*! \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
987
	 * \param opt option like VL_SYMMETRIC and VL_NON_SYMMETRIC or VL_CRS_SYMMETRIC
Pietro Incardona's avatar
Pietro Incardona committed
988 989 990 991
	 *
	 */
	void updateVerlet(VerletList<dim,St,FAST,shift<dim,St> > & ver, St r_cut, size_t opt = VL_NON_SYMMETRIC)
	{
incardon's avatar
incardon committed
992 993 994 995
#ifdef SE_CLASS3
		se3.getNN();
#endif

incardon's avatar
incardon committed
996 997 998 999 1000 1001
		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
1002
			bool to_reconstruct = NN.get_ndec() != getDecomposition().get_ndec();
incardon's avatar
incardon committed
1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014

			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
1015 1016 1017
		{
			auto & NN = ver.getInternalCellList();

incardon's avatar
incardon committed
1018 1019
			// 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
1020
			bool to_reconstruct = NN.get_ndec() != getDecomposition().get_ndec();
incardon's avatar
incardon committed
1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035

			if (to_reconstruct == false)
			{
				// Shift
				grid_key_dx<dim> cell_shift = NN.getShift();

				// Shift
				grid_key_dx<dim> shift = NN.getShift();

				// Add padding
				for (size_t i = 0 ; i < dim ; i++)
					shift.set_d(i,shift.get(i) + NN.getPadding(i));

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

incardon's avatar
incardon committed
1036 1037 1038
				ver.updateCrs(getDecomposition().getDomain(),r_cut,v_pos,g_m,
						      getDecomposition().getCRSDomainCells(shift,cell_shift,gs),
							  getDecomposition().getCRSAnomDomainCells(shift,cell_shift,gs));
incardon's avatar
incardon committed
1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053
			}
			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
1054
			bool to_reconstruct = NN.get_ndec() != getDecomposition().get_ndec();
incardon's avatar
incardon committed
1055 1056 1057 1058 1059 1060

			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
1061

incardon's avatar
incardon committed
1062 1063 1064 1065
				ver_tmp = getVerlet(r_cut);
				ver.swap(ver_tmp);
			}
		}
Pietro Incardona's avatar
Pietro Incardona committed
1066 1067
	}

incardon's avatar
incardon committed
1068

Yaroslav's avatar
Yaroslav committed
1069 1070 1071 1072 1073 1074 1075 1076 1077
	/*! \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
	 *
	 */
	template<typename CellL=CellList<dim,St,FAST,shift<dim,St> > > void reorder (int32_t m)
	{
Pietro Incardona's avatar
Pietro Incardona committed
1078
		reorder(m,getDecomposition().getGhost());
Yaroslav's avatar
Yaroslav committed
1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105
	}


	/*! \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
	 *
	 */
	template<typename CellL=CellList<dim,St,FAST,shift<dim,St> > > void reorder(int32_t m, const Ghost<dim,St> & enlarge)
	{
		// reset the ghost part
		v_pos.resize(g_m);
		v_prp.resize(g_m);


		CellL cell_list;

		// calculate the parameters of the cell list

		// get the processor bounding box
Pietro Incardona's avatar
Pietro Incardona committed
1106
		Box<dim,St> pbox = getDecomposition().getProcessorBounds();
Yaroslav's avatar
Yaroslav committed
1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117
		// extend by the ghost
		pbox.enlarge(enlarge);

		size_t div[dim];

		// Calculate the division array and the cell box
		for (size_t i = 0 ; i < dim ; i++)
		{
			div[i] = 1 << m;
		}

1118
		cell_list.Initialize(pbox,div);
Yaroslav's avatar
Yaroslav committed
1119 1120 1121 1122 1123 1124 1125 1126 1127

		// for each particle add the particle to the cell list

		auto it = getIterator();

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

incardon's avatar
incardon committed
1128 1129 1130
			Point<dim,St> xp = this->getPos(key);

			cell_list.add(xp,key.getKey());
Yaroslav's avatar
Yaroslav committed
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 1178

			++it;