vector_dist.hpp 31.7 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.hpp"
incardon's avatar
incardon committed
13
#include "Space/Shape/Point.hpp"
incardon's avatar
incardon committed
14 15 16
#include "Vector/vector_dist_iterator.hpp"
#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"
24
#include "VTKWriter/VTKWriter.hpp"
25
#include "Decomposition/common.hpp"
26
#include "Grid/grid_dist_id_iterator_dec.hpp"
27
#include "Grid/grid_key_dx_iterator_hilbert.hpp"
28
#include "Vector/vector_dist_ofb.hpp"
29
#include "Decomposition/CartDecomposition.hpp"
30
#include "data_type/aggregate.hpp"
31
#include "NN/VerletList/VerletList.hpp"
32
#include "vector_dist_comm.hpp"
incardon's avatar
incardon committed
33 34 35 36

#define NO_ID false
#define ID true

incardon's avatar
incardon committed
37
// Perform a ghost get or a ghost put
incardon's avatar
incardon committed
38 39 40
#define GET	1
#define PUT 2

incardon's avatar
incardon committed
41 42 43 44
// Write the particles with ghost
#define NO_GHOST 0
#define WITH_GHOST 2

incardon's avatar
incardon committed
45
/*! \brief Distributed vector
46
 *
47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
 * 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
67 68 69
 *
 */

70
template<unsigned int dim, typename St, typename prop, typename Decomposition = CartDecomposition<dim,St>, typename Memory = HeapMemory>
71
class vector_dist : public vector_dist_comm<dim,St,prop,Decomposition,Memory>
incardon's avatar
incardon committed
72 73 74
{
private:

75
	//! Ghost marker, all the particle with id > g_m are ghost all with g_m < are real particle
76 77
	size_t g_m = 0;

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

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

86
	//! Virtual cluster
incardon's avatar
incardon committed
87 88
	Vcluster & v_cl;

89

Pietro Incardona's avatar
Pietro Incardona committed
90
	/*! \brief Initialize the structures
incardon's avatar
incardon committed
91
	 *
Pietro Incardona's avatar
Pietro Incardona committed
92
	 * \param np number of particles
incardon's avatar
incardon committed
93 94
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
95
	void init_structures(size_t np)
incardon's avatar
incardon committed
96
	{
97
		// convert to a local number of elements
98 99 100 101 102 103 104 105
		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++;
106

incardon's avatar
incardon committed
107
		// resize the position vector
108
		v_pos.resize(p_np);
incardon's avatar
incardon committed
109 110

		// resize the properties vector
111 112 113
		v_prp.resize(p_np);

		g_m = p_np;
Pietro Incardona's avatar
Pietro Incardona committed
114
	}
incardon's avatar
incardon committed
115

116 117 118 119 120 121 122 123 124 125 126 127 128
	/*! \brief Checl if the parameters describe a valid vector. In case it does not report an error
	 *
	 * \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)
			std::cerr << __FILE__ << ":" << __LINE__ << "  Error the domain is not valid " << box.toString() << std::endl;

	}

Pietro Incardona's avatar
Pietro Incardona committed
129 130 131 132 133 134 135 136
public:

	//! space type
	typedef St stype;

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

137 138 139 140
	/*! \brief Operator= for distributed vector
	 *
	 * \param v vector to copy
	 *
Pietro Incardona's avatar
Pietro Incardona committed
141 142
	 * \return itself
	 *
143 144 145
	 */
	vector_dist<dim,St,prop,Decomposition,Memory> & operator=(const vector_dist<dim,St,prop,Decomposition,Memory> & v)
	{
146
		static_cast<vector_dist_comm<dim,St,prop,Decomposition,Memory> *>(this)->operator=(static_cast<vector_dist_comm<dim,St,prop,Decomposition,Memory>>(v));
147 148 149 150 151 152 153 154 155 156 157 158

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

		return *this;
	}

	/*! \brief Operator= for distributed vector
	 *
	 * \param v vector to copy
	 *
Pietro Incardona's avatar
Pietro Incardona committed
159 160
	 * \return itself
	 *
161 162 163
	 */
	vector_dist<dim,St,prop,Decomposition,Memory> & operator=(vector_dist<dim,St,prop,Decomposition,Memory> && v)
	{
164
		static_cast<vector_dist_comm<dim,St,prop,Decomposition,Memory> *>(this)->operator=(static_cast<vector_dist_comm<dim,St,prop,Decomposition,Memory> >(v));
165 166 167 168 169 170 171 172

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

		return *this;
	}

Pietro Incardona's avatar
Pietro Incardona committed
173
	/*! \brief Copy Constructor
174
	 *
Pietro Incardona's avatar
Pietro Incardona committed
175
	 * \param v vector to copy
176 177 178
	 *
	 */
	vector_dist(const vector_dist<dim,St,prop,Decomposition,Memory> & v)
incardon's avatar
incardon committed
179
	:vector_dist_comm<dim,St,prop,Decomposition,Memory>(v.getDecomposition()),v_cl(v.v_cl)
180 181 182 183 184 185 186 187
	{
#ifdef SE_CLASS2
		check_new(this,8,VECTOR_DIST_EVENT,4);
#endif

		this->operator=(v);
	}

Pietro Incardona's avatar
Pietro Incardona committed
188
	/*! \brief Copy constructor
189
	 *
Pietro Incardona's avatar
Pietro Incardona committed
190
	 * \param v vector to copy
191 192
	 *
	 */
incardon's avatar
incardon committed
193
	vector_dist(vector_dist<dim,St,prop,Decomposition,Memory> && v) noexcept
194 195 196 197 198 199 200 201
	:v_cl(v.v_cl)
	{
#ifdef SE_CLASS2
		check_new(this,8,VECTOR_DIST_EVENT,4);
#endif

		this->operator=(v);
	}
Pietro Incardona's avatar
Pietro Incardona committed
202

Pietro Incardona's avatar
Pietro Incardona committed
203
	/*! \brief Constructor with predefined decomposition
Pietro Incardona's avatar
Pietro Incardona committed
204
	 *
Pietro Incardona's avatar
Pietro Incardona committed
205 206
	 * \param dec is the decomposition
	 * \param np number of particles
Pietro Incardona's avatar
Pietro Incardona committed
207 208 209
	 *
	 */
	vector_dist(const Decomposition & dec, size_t np) :
210
	vector_dist_comm<dim,St,prop,Decomposition,Memory>(dec), v_cl(create_vcluster())
Pietro Incardona's avatar
Pietro Incardona committed
211 212 213 214 215 216 217 218 219 220 221 222 223
	{
#ifdef SE_CLASS2
		check_new(this,8,VECTOR_DIST_EVENT,4);
#endif

		init_structures(np);
	}


	/*! \brief Constructor
	 *
	 * \param np number of elements
	 * \param box domain where the vector of elements live
Pietro Incardona's avatar
Pietro Incardona committed
224
	 * \param bc boundary conditions
Pietro Incardona's avatar
Pietro Incardona committed
225
	 * \param g Ghost margins
226 227 228 229
	 * \param opt additional options.
	 *        * BIND_DEC_TO_GHOST Bind the decomposition to be multiple of the
	 *          ghost size. This is required if we want to use symmetric to eliminate
	 *          ghost communications.
Pietro Incardona's avatar
Pietro Incardona committed
230 231
	 *
	 */
232
	vector_dist(size_t np, Box<dim, St> box, const size_t (&bc)[dim], const Ghost<dim, St> & g, size_t opt = 0)
233
	:v_cl(create_vcluster())
Pietro Incardona's avatar
Pietro Incardona committed
234 235 236 237 238
	{
#ifdef SE_CLASS2
		check_new(this,8,VECTOR_DIST_EVENT,4);
#endif

239 240
		check_parameters(box);

Pietro Incardona's avatar
Pietro Incardona committed
241
		init_structures(np);
242
		this->init_decomposition(box,bc,g,opt);
Pietro Incardona's avatar
Pietro Incardona committed
243 244
	}

245 246 247 248 249 250 251
	~vector_dist()
	{
#ifdef SE_CLASS2
		check_delete(this);
#endif
	}

252
	/*! \brief return the local size of the vector
253
	 *
254
	 * \return local size
255 256
	 *
	 */
257
	size_t size_local()
258
	{
259
		return g_m;
260 261
	}

incardon's avatar
incardon committed
262 263 264 265 266
	/*! \brief return the local size of the vector
	 *
	 * \return local size
	 *
	 */
267
	size_t size_local_with_ghost()
incardon's avatar
incardon committed
268
	{
269
		return v_pos.size();
incardon's avatar
incardon committed
270 271
	}

272
	/*! \brief Get the position of an element
incardon's avatar
incardon committed
273
	 *
274 275 276 277 278
	 * 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
279 280
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
281
	inline auto getPos(vect_dist_key_dx vec_key) -> decltype(v_pos.template get<0>(vec_key.getKey()))
incardon's avatar
incardon committed
282
	{
Pietro Incardona's avatar
Pietro Incardona committed
283
		return v_pos.template get<0>(vec_key.getKey());
incardon's avatar
incardon committed
284 285
	}

286 287 288 289 290 291 292 293 294 295 296 297 298 299
	/*! \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());
	}

300
	/*! \brief Get the property of an element
incardon's avatar
incardon committed
301
	 *
302 303 304
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \tparam id property id
incardon's avatar
incardon committed
305 306
	 * \param vec_key vector element
	 *
307 308
	 * \return return the selected property of the vector element
	 *
incardon's avatar
incardon committed
309
	 */
310
	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
311
	{
312
		return v_prp.template get<id>(vec_key.getKey());
incardon's avatar
incardon committed
313 314
	}

315 316 317 318 319 320 321 322 323 324
	/*! \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
	 *
	 */
325
	template<unsigned int id> inline auto getProp(vect_dist_key_dx vec_key) const -> const decltype(v_prp.template get<id>(vec_key.getKey()))
326 327 328 329
	{
		return v_prp.template get<id>(vec_key.getKey());
	}

330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350
	/*! \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++;
	}

	/*! \brief Get the position of the last element
	 *
	 * \return the position of the element in space
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
351
	inline auto getLastPos() -> decltype(v_pos.template get<0>(0))
352
	{
Pietro Incardona's avatar
Pietro Incardona committed
353
		return v_pos.template get<0>(g_m - 1);
354 355 356 357 358 359 360 361 362 363 364
	}

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

368 369 370 371 372 373 374 375 376 377 378 379 380 381 382
	/*! \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)
	{
		// Cell list
		CellL cell_list;

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

incardon's avatar
incardon committed
386
		// Processor bounding box
387 388
		Box<dim, St> pbox = getDecomposition().getProcessorBounds();

incardon's avatar
incardon committed
389 390 391
		// Ghost padding extension
		Ghost<dim,size_t> g_ext(0);
		cell_list.Initialize(cd_sm,pbox,pad);
392 393

		updateCellList(cell_list);
incardon's avatar
incardon committed
394 395

		return cell_list;
396 397
	}

398 399 400 401
	/*! \brief Construct a cell list starting from the stored particles
	 *
	 * \tparam CellL CellList type to construct
	 *
402 403 404 405
	 * \param r_cut interation radius, or size of each cell
	 *
	 * \return the Cell list
	 *
406
	 */
407
	template<typename CellL = CellList<dim, St, FAST, shift<dim, St> > > CellL getCellList(St r_cut)
408
	{
Pietro Incardona's avatar
Pietro Incardona committed
409
		// Get ghost and anlarge by 1%
410
		Ghost<dim,St> g = getDecomposition().getGhost();
411
		g.magnify(1.013);
Pietro Incardona's avatar
Pietro Incardona committed
412 413

		return getCellList(r_cut, g);
414 415
	}

416 417 418 419 420 421 422 423 424 425 426 427
	/*! \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)
	{
		// Get ghost and anlarge by 1%
428
		Ghost<dim,St> g = getDecomposition().getGhost();
Pietro Incardona's avatar
Pietro Incardona committed
429
		g.magnify(1.013);
430 431 432 433

		return getCellList_hilb(r_cut, g);
	}

Pietro Incardona's avatar
Pietro Incardona committed
434 435 436 437 438 439 440 441 442
	/*! \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
443
		populate_cell_list(v_pos,cell_list,g_m,CL_NON_SYMMETRIC);
Pietro Incardona's avatar
Pietro Incardona committed
444

incardon's avatar
incardon committed
445 446 447 448 449 450 451 452 453 454 455 456
		cell_list.set_gm(g_m);
	}

	/*! \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
457
		populate_cell_list(v_pos,cell_list,g_m,CL_SYMMETRIC);
incardon's avatar
incardon committed
458

Pietro Incardona's avatar
Pietro Incardona committed
459
		cell_list.set_gm(g_m);
Pietro Incardona's avatar
Pietro Incardona committed
460 461
	}

462 463 464 465 466 467 468 469 470 471 472 473
	/*! \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
474 475
	 * \return the CellList
	 *
476
	 */
477
	template<typename CellL = CellList<dim, St, FAST, shift<dim, St> > > CellL getCellList(St r_cut, const Ghost<dim, St> & enlarge)
478 479 480
	{
		CellL cell_list;

481 482
		// Division array
		size_t div[dim];
483

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

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

490 491 492 493 494 495
		cell_list.Initialize(pbox, div);

		updateCellList(cell_list);

		return cell_list;
	}
496

497 498 499 500 501 502 503 504 505 506 507 508
	/*! \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
509 510
	 * \return The Cell-list
	 *
511 512 513 514 515 516
	 */
	template<typename CellL = CellList_hilb<dim, St, FAST, shift<dim, St> > > CellL getCellList_hilb(St r_cut, const Ghost<dim, St> & enlarge)
	{
		CellL cell_list;

		// Division array
517 518
		size_t div[dim];

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

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

525
		cell_list.Initialize(pbox, div, g_m);
526

Pietro Incardona's avatar
Pietro Incardona committed
527
		updateCellList(cell_list);
528 529 530

		return cell_list;
	}
incardon's avatar
incardon committed
531

incardon's avatar
incardon committed
532 533 534 535
	/*! \brief for each particle get the symmetric verlet list
	 *
	 * \param r_cut cut-off radius
	 *
incardon's avatar
incardon committed
536 537
	 * \return the verlet list
	 *
incardon's avatar
incardon committed
538 539 540 541
	 */
	VerletList<dim,St,FAST,shift<dim,St> > getVerletSym(St r_cut)
	{
		VerletList<dim,St,FAST,shift<dim,St>> ver;
542

incardon's avatar
incardon committed
543 544 545 546 547 548 549
		// Processor bounding box
		Box<dim, St> pbox = getDecomposition().getProcessorBounds();

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

		return ver;
	}
550

incardon's avatar
incardon committed
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
	/*! \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)
	{
		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();

		ver.createVerletCrs(r_cut,g_m,v_pos,getDecomposition().getDomainCells(shift,cell_shift,gs),getDecomposition().getAnomDomainCells(shift,cell_shift,gs));

		return ver;
	}

588 589 590 591
	/*! \brief for each particle get the verlet list
	 *
	 * \param r_cut cut-off radius
	 *
incardon's avatar
incardon committed
592 593
	 * \return a VerletList object
	 *
594
	 */
incardon's avatar
incardon committed
595
	VerletList<dim,St,FAST,shift<dim,St> > getVerlet(St r_cut)
596 597 598 599
	{
		VerletList<dim,St,FAST,shift<dim,St>> ver;

		// get the processor bounding box
600 601 602 603 604
		Box<dim, St> bt = getDecomposition().getProcessorBounds();

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

606 607
		// enlarge the box where the Verlet is defined
		bt.enlarge(g);
608

incardon's avatar
incardon committed
609
		ver.Initialize(bt,getDecomposition().getProcessorBounds(),r_cut,v_pos,g_m,VL_NON_SYMMETRIC);
610 611 612 613

		return ver;
	}

614 615 616 617 618 619 620 621 622 623
	/*! \brief for each particle get the verlet list
	 *
	 * \param r_cut cut-off radius
	 * \param ver Verlet to update
	 * \param r_cut cutoff radius
	 * \param opt option like VL_SYMMETRIC and VL_NON_SYMMETRIC
	 *
	 */
	void updateVerlet(VerletList<dim,St,FAST,shift<dim,St> > & ver, St r_cut, size_t opt = VL_NON_SYMMETRIC)
	{
incardon's avatar
incardon committed
624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644
		if (opt == VL_CRS_SYMMETRIC)
		{
			// 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();

			ver.updateCrs(getDecomposition().getDomain(),r_cut,v_pos,g_m,getDecomposition().getDomainCells(shift,cell_shift,gs),getDecomposition().getAnomDomainCells(shift,cell_shift,gs));
		}
		else
			ver.update(getDecomposition().getDomain(),r_cut,v_pos,g_m, opt);
645 646
	}

647 648 649 650 651
	/*! \brief for each particle get the verlet list
	 *
	 * \param verlet output verlet list for each particle
	 * \param r_cut cut-off radius
	 *
652 653
	 * \deprecated
	 *
654
	 */
655
	void getVerletDeprecated(openfpm::vector<openfpm::vector<size_t>> & verlet, St r_cut)
656 657 658 659 660 661 662 663
	{
		// resize verlet to store the number of particles
		verlet.resize(size_local());

		// get the cell-list
		auto cl = getCellList(r_cut);

		// square of the cutting radius
664
		St r_cut2 = r_cut * r_cut;
665 666

		// iterate the particles
667 668 669 670 671 672 673
		auto it_p = this->getDomainIterator();
		while (it_p.isNext())
		{
			// key
			vect_dist_key_dx key = it_p.get();

			// Get the position of the particles
Pietro Incardona's avatar
Pietro Incardona committed
674
			Point<dim, St> p = this->getPos(key);
675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691

			// Clear the neighborhood of the particle
			verlet.get(key.getKey()).clear();

			// Get the neighborhood of the particle
			auto NN = cl.template getNNIterator<NO_CHECK>(cl.getCell(p));
			while (NN.isNext())
			{
				auto nnp = NN.get();

				// p != q
				if (nnp == key.getKey())
				{
					++NN;
					continue;
				}

Pietro Incardona's avatar
Pietro Incardona committed
692
				Point<dim, St> q = this->getPos(nnp);
693 694 695 696 697 698 699 700 701 702 703

				if (p.distance2(q) < r_cut2)
					verlet.get(key.getKey()).add(nnp);

				// Next particle
				++NN;
			}

			// next particle
			++it_p;
		}
704 705
	}

706 707 708 709 710 711 712 713 714 715 716
	/*! \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)
	{
717
		reorder(m,getDecomposition().getGhost());
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
	}


	/*! \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
745
		Box<dim,St> pbox = getDecomposition().getProcessorBounds();
746 747 748 749 750 751 752 753 754 755 756
		// 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;
		}

757
		cell_list.Initialize(pbox,div);
758 759 760 761 762 763 764 765 766

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

		auto it = getIterator();

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

767
			cell_list.add(this->getPos(key),key.getKey());
768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815

			++it;
		}

		// Use cell_list to reorder v_pos

		//destination vector
		openfpm::vector<Point<dim,St>> v_pos_dest;
		openfpm::vector<prop> v_prp_dest;

		v_pos_dest.resize(v_pos.size());
		v_prp_dest.resize(v_prp.size());

		//hilberts curve iterator
		grid_key_dx_iterator_hilbert<dim> h_it(m);

		//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;
		}

		v_pos.swap(v_pos_dest);
		v_prp.swap(v_prp_dest);
	}

816 817
	/*! \brief It return the number of particles contained by the previous processors
	 *
Pietro Incardona's avatar
Pietro Incardona committed
818
	 * \warning It only work with the initial decomposition
819 820 821 822 823 824 825 826 827
	 *
	 * Given 1000 particles and 3 processors, you will get
	 *
	 * * Processor 0: 0
	 * * Processor 1: 334
	 * * Processor 2: 667
	 *
	 * \param np initial number of particles
	 *
Pietro Incardona's avatar
Pietro Incardona committed
828 829
	 * \return number of particles contained by the previous processors
	 *
830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850
	 */
	size_t init_size_accum(size_t np)
	{
		size_t accum = 0;

		// convert to a local number of elements
		size_t p_np = np / v_cl.getProcessingUnits();

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

		accum = p_np * v_cl.getProcessUnitID();

		// Distribute the remain particles
		if (v_cl.getProcessUnitID() <= r)
			accum += v_cl.getProcessUnitID();
		else
			accum += r;

		return accum;
	}
incardon's avatar
incardon committed
851

852
	/*! \brief Get an iterator that traverse domain and ghost particles
incardon's avatar
incardon committed
853 854 855 856
	 *
	 * \return an iterator
	 *
	 */
857
	vector_dist_iterator getIterator()
incardon's avatar
incardon committed
858
	{
859
		return vector_dist_iterator(0, v_pos.size());
incardon's avatar
incardon committed
860 861
	}

862 863 864 865
	/*! /brief Get a grid Iterator
	 *
	 * Usefull function to place particles on a grid or grid-like (grid + noise)
	 *
Pietro Incardona's avatar
Pietro Incardona committed
866 867
	 * \param sz size of the grid
	 *
868 869 870
	 * \return a Grid iterator
	 *
	 */
871
	inline grid_dist_id_iterator_dec<Decomposition> getGridIterator(const size_t (&sz)[dim])
872 873 874
	{
		grid_key_dx<dim> start;
		grid_key_dx<dim> stop;
875
		for (size_t i = 0; i < dim; i++)
876
		{
877
			start.set_d(i, 0);
878
			stop.set_d(i, sz[i] - 1);
879 880
		}

881
		grid_dist_id_iterator_dec<Decomposition> it_dec(getDecomposition(), sz, start, stop);
882 883 884
		return it_dec;
	}

885 886 887 888 889
	/*! \brief Get the iterator across the position of the ghost particles
	 *
	 * \return an iterator
	 *
	 */
890
	vector_dist_iterator getGhostIterator() const
891
	{
892
		return vector_dist_iterator(g_m, v_pos.size());
893 894
	}

895
	/*! \brief Get an iterator that traverse the particles in the domain
896 897 898 899
	 *
	 * \return an iterator
	 *
	 */
900
	vector_dist_iterator getDomainIterator() const
901
	{
902
		return vector_dist_iterator(0, g_m);
903 904
	}

905 906 907 908 909 910 911 912 913 914
	/*! \brief Get an iterator that traverse the particles in the domain
	 *
	 * \return an iterator
	 *
	 */
	vector_dist_iterator getDomainAndGhostIterator() const
	{
		return vector_dist_iterator(0, v_pos.size());
	}

incardon's avatar
incardon committed
915 916 917 918 919
	/*! \brief Get the decomposition
	 *
	 * \return
	 *
	 */
920 921 922 923 924
	inline Decomposition & getDecomposition()
	{
		return vector_dist_comm<dim,St,prop,Decomposition,Memory>::getDecomposition();
	}

incardon's avatar
incardon committed
925 926 927 928 929 930 931 932 933 934
	/*! \brief Get the decomposition
	 *
	 * \return
	 *
	 */
	inline const Decomposition & getDecomposition() const
	{
		return vector_dist_comm<dim,St,prop,Decomposition,Memory>::getDecomposition();
	}

935 936 937 938 939 940 941 942 943 944 945 946 947
	/*! \brief It move all the particles that does not belong to the local processor to the respective processor
	 *
	 * \tparam out of bound policy it specify what to do when the particles are detected out of bound
	 *
	 * In general this function is called after moving the particles to move the
	 * elements out the local processor. Or just after initialization if each processor
	 * contain non local particles
	 *
	 * \tparam prp properties to communicate
	 *
	 *
	 */
	template<unsigned int ... prp> void map_list()
incardon's avatar
incardon committed
948
	{
949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977
		this->template map_list_<prp...>(v_pos,v_prp,g_m);
	}


	/*! \brief It move all the particles that does not belong to the local processor to the respective processor
	 *
	 * \tparam out of bound policy it specify what to do when the particles are detected out of bound
	 *
	 * In general this function is called after moving the particles to move the
	 * elements out the local processor. Or just after initialization if each processor
	 * contain non local particles
	 *
	 *
	 */
	template<typename obp = KillParticle> void map()
	{
		this->template map_<obp>(v_pos,v_prp,g_m);
	}

	/*! \brief It synchronize the properties and position of the ghost particles
	 *
	 * \tparam prp list of properties to get synchronize
	 *
	 * \param opt options WITH_POSITION, it send also the positional information of the particles
	 *
	 */
	template<int ... prp> inline void ghost_get(size_t opt = WITH_POSITION)
	{
		this->template ghost_get_<prp...>(v_pos,v_prp,g_m,opt);
incardon's avatar
incardon committed
978
	}
incardon's avatar
incardon committed
979

Pietro Incardona's avatar
Pietro Incardona committed
980 981 982 983 984 985 986 987 988 989 990 991
	/*! \brief It synchronize the properties and position of the ghost particles
	 *
	 * \tparam op which kind of operation to apply
	 * \tparam prp list of properties to get synchronize
	 *
	 *
	 */
	template<template<typename,typename> class op, int ... prp> inline void ghost_put()
	{
		this->template ghost_put_<op,prp...>(v_pos,v_prp,g_m);
	}

Pietro Incardona's avatar
Pietro Incardona committed
992 993 994 995 996 997 998 999
	/*! \brief Remove a set of elements from the distributed vector
	 *
	 * \warning keys must be sorted
	 *
	 * \param keys vector of elements to eliminate
	 * \param start from where to eliminate
	 *
	 */
1000 1001 1002 1003
	void remove(openfpm::vector<size_t> & keys, size_t start = 0)
	{
		v_pos.remove(keys, start);
		v_prp.remove(keys, start);
1004 1005

		g_m -= keys.size();
1006 1007
	}

Pietro Incardona's avatar
Pietro Incardona committed
1008 1009
	/*! \brief Remove one element from the distributed vector
	 *
Pietro Incardona's avatar
Pietro Incardona committed
1010
	 * \param key remove one element from the vector
Pietro Incardona's avatar
Pietro Incardona committed
1011 1012
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
1013 1014 1015 1016 1017 1018 1019 1020
	void remove(size_t key)
	{
		v_pos.remove(key);
		v_prp.remove(key);

		g_m--;
	}

Pietro Incardona's avatar
Pietro Incardona committed
1021 1022 1023
	/*! \brief Add the computation cost on the decomposition comming from the particles
	 *
	 */
1024 1025 1026 1027
	inline void addComputationCosts()
	{
		CellDecomposer_sm<dim, St> cdsm;

1028 1029
		Decomposition & dec = getDecomposition();

1030 1031
		cdsm.setDimensions(dec.getDomain(), dec.getGrid().getSize(), 0);

1032
		for (size_t i = 0; i < getDecomposition().getNSubSubDomains(); i++)
1033
		{
Pietro Incardona's avatar
Pietro Incardona committed
1034
			dec.setSubSubDomainComputationCost(i, 1);
1035 1036
		}

tonynsyde's avatar
tonynsyde committed
1037 1038
		auto it = getDomainIterator();

1039 1040
		while (it.isNext())
		{
Pietro Incardona's avatar
Pietro Incardona committed
1041
			size_t v = cdsm.getCell(this->getPos(it.get()));
1042 1043 1044 1045 1046 1047 1048 1049

			dec.addComputationCost(v, 1);

			++it;
		}

	}

incardon's avatar
incardon committed
1050 1051
	/*! \brief Output particle position and properties
	 *
1052
	 * \param out output
1053
	 * \param opt VTK_WRITER or CSV_WRITER
incardon's avatar
incardon committed
1054
	 *
1055
	 * \return true if the file has been written without error
incardon's avatar
incardon committed
1056 1057
	 *
	 */
1058
	inline bool write(std::string out, int opt = NO_GHOST | VTK_WRITER )
incardon's avatar
incardon committed
1059 1060
	{

1061 1062 1063 1064
		if ((opt & 0xFFFF0000) == CSV_WRITER)
		{
			// CSVWriter test
			CSVWriter<openfpm::vector<Point<dim,St>>, openfpm::vector<prop> > csv_writer;
incardon's avatar
incardon committed
1065

Pietro Incardona's avatar
Pietro Incardona committed
1066
			std::string output = std::to_string(out + "_" + std::to_string(v_cl.getProcessUnitID()) + std::to_string(".csv"));
1067 1068 1069 1070

			// Write the CSV
			return csv_writer.write(output,v_pos,v_prp);
		}
1071
		else if ((opt & 0xFFFF0000) == VTK_WRITER)
1072
		{
1073
			// VTKWriter for a set of points
1074
			VTKWriter<boost::mpl::pair<openfpm::vector<Point<dim,St>>, openfpm::vector<prop>>, VECTOR_POINTS> vtk_writer;
1075
			vtk_writer.add(v_pos,v_prp,g_m);
incardon's avatar
incardon committed
1076

Pietro Incardona's avatar
Pietro Incardona committed
1077
			std::string output = std::to_string(out + "_" + std::to_string(v_cl.getProcessUnitID()) + std::to_string(".vtk"));
incardon's avatar
incardon committed
1078

1079 1080
			// Write the VTK file
			return vtk_writer.write(output);
1081
		}
1082 1083

		return false;
1084 1085
	}

Pietro Incardona's avatar
Pietro Incardona committed
1086 1087 1088 1089 1090 1091 1092 1093 1094 1095
	/*! \brief Delete the particles on the ghost
	 *
	 *
	 */
	void deleteGhost()
	{
		v_pos.resize(g_m);
		v_prp.resize(g_m);
	}

1096 1097
	/*! \brief Resize the vector (locally)
	 *
Pietro Incardona's avatar
Pietro Incardona committed
1098
	 * \warning It automatically delete the ghosts
1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112
	 *
	 * \param rs
	 *
	 */
	void resize(size_t rs)
	{
		deleteGhost();

		v_pos.resize(rs);
		v_prp.resize(rs);

		g_m = rs;
	}

1113 1114 1115
	/*! \brief Output particle position and properties
	 *
	 * \param out output
Pietro Incardona's avatar
Pietro Incardona committed
1116
	 * \param iteration (we can append the number at the end of the file_name)
1117 1118 1119 1120 1121 1122 1123
	 * \param opt NO_GHOST or WITH_GHOST
	 *
	 * \return if the file has been written correctly
	 *
	 */
	inline bool write(std::string out, size_t iteration, int opt = NO_GHOST)
	{
Pietro Incardona's avatar
Pietro Incardona committed
1124 1125 1126 1127
		if ((opt & 0xFFFF0000) == CSV_WRITER)
		{
			// CSVWriter test
			CSVWriter<openfpm::vector<Point<dim, St>>, openfpm::vector<prop> > csv_writer;
1128

1129
			std::string output = std::to_string(out + "_" + std::to_string(v_cl.getProcessUnitID()) + "_" + std::to_string(iteration) + std::to_string(".csv"));
1130

Pietro Incardona's avatar
Pietro Incardona committed
1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144
			// Write the CSV
			return csv_writer.write(output, v_pos, v_prp);
		}
		else
		{
			// VTKWriter for a set of points
			VTKWriter<boost::mpl::pair<openfpm::vector<Point<dim,St>>, openfpm::vector<prop>>, VECTOR_POINTS> vtk_writer;
			vtk_writer.add(v_pos,v_prp,g_m);

			std::string output = std::to_string(out + "_" + std::to_string(v_cl.getProcessUnitID()) + "_" + std::to_string(iteration) + std::to_string(".vtk"));

			// Write the VTK file
			return vtk_writer.write(output);
		}
incardon's avatar
incardon committed
1145
	}
1146

1147 1148 1149
	/*! \brief Get the Celllist parameters
	 *
	 * \param r_cut spacing of the cell-list
Pietro Incardona's avatar
Pietro Incardona committed
1150
	 * \param div division required for the cell-list
1151
	 * \param box where the Cell list must be defined (In general Processor domain + Ghost)
Pietro Incardona's avatar
Pietro Incardona committed
1152
	 * \param enlarge Optionally a request to make the space a littler bit larger than Processor domain + Ghost
1153 1154 1155 1156 1157
	 *        keeping the cell list consistent with the requests
	 *
	 */
	void getCellListParams(St r_cut, size_t (&div)[dim],Box<dim, St> & box, Ghost<dim,St> enlarge = Ghost<dim,St>(0.0))
	{
incardon's avatar
incardon committed
1158 1159 1160 1161 1162 1163 1164 1165
		// get the processor bounding box
		Box<dim, St> pbox = getDecomposition().getProcessorBounds();

		// enlarge the processor bounding box by the ghost
		Ghost<dim,St> g = getDecomposition().getGhost();
		pbox.enlarge(g);

		cl_param_calculate(pbox, div,r_cut,enlarge);
incardon's avatar
incardon committed
1166 1167 1168

		// output the fixed domain
		box = pbox;
1169 1170
	}

Pietro Incardona's avatar
Pietro Incardona committed
1171
	/*! \brief It return the id of structure in the allocation list
1172 1173 1174
	 *
	 * \see print_alloc and SE_CLASS2
	 *
Pietro Incardona's avatar
Pietro Incardona committed
1175 1176
	 * \return the id
	 *
1177 1178 1179 1180 1181
	 */
	long int who()
	{
#ifdef SE_CLASS2
		return check_whoami(this,8);
1182 1183
#else
		return -1;
1184 1185
#endif
	}
1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199

	/*! \brief Get the Virtual Cluster machine
	 *
	 * \return the Virtual cluster machine
	 *
	 */

	Vcluster & getVC()
	{
#ifdef SE_CLASS2
		check_valid(this,8);
#endif
		return v_cl;
	}
incardon's avatar
incardon committed
1200 1201 1202 1203 1204 1205 1206 1207 1208 1209

	/*! \brief return the position vector of all the particles
	 *
	 * \return the particle position vector
	 *
	 */
	const openfpm::vector<Point<dim,St>> & getPosVector() const
	{
		return v_pos;
	}
incardon's avatar
incardon committed
1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231

	/*! \brief It return the sum of the particles in the previous processors
	 *
	 * \return the particles number
	 *
	 */
	size_t accum()
	{
		openfpm::vector<size_t> accu;

		size_t sz = size_local();

		v_cl.allGather(sz,accu);
		v_cl.execute();

		sz = 0;

		for (size_t i = 0 ; i < v_cl.getProcessUnitID() ; i++)
			sz += accu.get(i);

		return sz;
	}
1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285

	/*! \brief Get a special particle iterator able to iterate across particles using
	 *         symmetric crossing scheme
	 *
	 * \param NN Cell-List neighborhood
	 *
	 */
	template<typename cli> ParticleItCRS_Cells<dim,cli> getParticleIteratorCRS(cli & NN)
	{
		// 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();

		// First we check that
		return ParticleItCRS_Cells<dim,cli>(NN,getDecomposition().getDomainCells(shift,cell_shift,gs),getDecomposition().getAnomDomainCells(shift,cell_shift,gs),NN.getNNc_sym());
	}

	/*! \brief Return from which cell we have to start in case of CRS interation
	 *         scheme
	 *
	 * \param NN cell-list
	 *
	 * \return The starting cell point
	 *
	 */
	template<typename Celllist> grid_key_dx<dim> getCRSStart(Celllist & NN)
	{
		return NN.getStartDomainCell();
	}

	/*! \brief Return from which cell we have to stop in case of CRS interation
	 *         scheme
	 *
	 * \param NN cell-list
	 *
	 * \return The stop cell point
	 *
	 */
	template<typename Celllist> grid_key_dx<dim> getCRSStop(Celllist & NN)
	{
		grid_key_dx<dim> key = NN.getStopDomainCell();

		for (size_t i = 0 ; i < dim ; i++)
			key.set_d(i,key.get(i) + 1);
		return key;
	}
incardon's avatar
incardon committed
1286 1287
};

1288

incardon's avatar
incardon committed
1289
#endif /* VECTOR_HPP_ */