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

#ifndef VECTOR_HPP_
#define VECTOR_HPP_

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
35
36
37

#define NO_ID false
#define ID true

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

incardon's avatar
incardon committed
40
// Perform a ghost get or a ghost put
incardon's avatar
incardon committed
41
42
43
#define GET	1
#define PUT 2

incardon's avatar
incardon committed
44
45
46
47
// Write the particles with ghost
#define NO_GHOST 0
#define WITH_GHOST 2

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

73
template<unsigned int dim, typename St, typename prop, typename Decomposition = CartDecomposition<dim,St>, typename Memory = HeapMemory>
Pietro Incardona's avatar
Pietro Incardona committed
74
class vector_dist : public vector_dist_comm<dim,St,prop,Decomposition,Memory>
incardon's avatar
incardon committed
75
76
77
{
private:

78
	//! Ghost marker, all the particle with id > g_m are ghost all with g_m < are real particle
79
80
	size_t g_m = 0;

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

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

89
	//! Virtual cluster
incardon's avatar
incardon committed
90
91
	Vcluster & v_cl;

92

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

incardon's avatar
incardon committed
110
		// resize the position vector
111
		v_pos.resize(p_np);
incardon's avatar
incardon committed
112
113

		// resize the properties vector
114
115
116
		v_prp.resize(p_np);

		g_m = p_np;
Pietro Incardona's avatar
Pietro Incardona committed
117
	}
incardon's avatar
incardon committed
118

incardon's avatar
incardon committed
119
120
121
122
123
124
125
126
127
128
129
130
131
	/*! \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
132
133
134
135
136
137
138
139
public:

	//! space type
	typedef St stype;

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

Pietro Incardona's avatar
Pietro Incardona committed
140
141
142
143
	/*! \brief Operator= for distributed vector
	 *
	 * \param v vector to copy
	 *
Pietro Incardona's avatar
Pietro Incardona committed
144
145
	 * \return itself
	 *
Pietro Incardona's avatar
Pietro Incardona committed
146
147
148
	 */
	vector_dist<dim,St,prop,Decomposition,Memory> & operator=(const vector_dist<dim,St,prop,Decomposition,Memory> & v)
	{
Pietro Incardona's avatar
Pietro Incardona committed
149
		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
150
151
152
153
154
155
156
157
158
159
160
161

		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
162
163
	 * \return itself
	 *
Pietro Incardona's avatar
Pietro Incardona committed
164
165
166
	 */
	vector_dist<dim,St,prop,Decomposition,Memory> & operator=(vector_dist<dim,St,prop,Decomposition,Memory> && v)
	{
Pietro Incardona's avatar
Pietro Incardona committed
167
		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
168
169
170
171
172
173
174
175

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

		this->operator=(v);
	}

Pietro Incardona's avatar
Pietro Incardona committed
191
	/*! \brief Copy constructor
Pietro Incardona's avatar
Pietro Incardona committed
192
	 *
Pietro Incardona's avatar
Pietro Incardona committed
193
	 * \param v vector to copy
Pietro Incardona's avatar
Pietro Incardona committed
194
195
	 *
	 */
incardon's avatar
incardon committed
196
	vector_dist(vector_dist<dim,St,prop,Decomposition,Memory> && v) noexcept
Pietro Incardona's avatar
Pietro Incardona committed
197
198
199
200
201
202
203
204
	: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
205

Pietro Incardona's avatar
Pietro Incardona committed
206
	/*! \brief Constructor with predefined decomposition
Pietro Incardona's avatar
Pietro Incardona committed
207
	 *
Pietro Incardona's avatar
Pietro Incardona committed
208
209
	 * \param dec is the decomposition
	 * \param np number of particles
Pietro Incardona's avatar
Pietro Incardona committed
210
211
212
	 *
	 */
	vector_dist(const Decomposition & dec, size_t np) :
Pietro Incardona's avatar
Pietro Incardona committed
213
	vector_dist_comm<dim,St,prop,Decomposition,Memory>(dec), v_cl(create_vcluster())
Pietro Incardona's avatar
Pietro Incardona committed
214
215
216
217
218
219
220
221
222
223
224
225
226
	{
#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
227
	 * \param bc boundary conditions
Pietro Incardona's avatar
Pietro Incardona committed
228
	 * \param g Ghost margins
229
230
231
232
	 * \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
233
234
	 *
	 */
235
	vector_dist(size_t np, Box<dim, St> box, const size_t (&bc)[dim], const Ghost<dim, St> & g, size_t opt = 0)
Pietro Incardona's avatar
Pietro Incardona committed
236
	:v_cl(create_vcluster())
Pietro Incardona's avatar
Pietro Incardona committed
237
238
239
240
241
	{
#ifdef SE_CLASS2
		check_new(this,8,VECTOR_DIST_EVENT,4);
#endif

incardon's avatar
incardon committed
242
243
244
		if (opt >> 32 != 0)
			this->setDecompositionGranularity(opt >> 32);

incardon's avatar
incardon committed
245
246
		check_parameters(box);

Pietro Incardona's avatar
Pietro Incardona committed
247
		init_structures(np);
248
		this->init_decomposition(box,bc,g,opt);
Pietro Incardona's avatar
Pietro Incardona committed
249
250
	}

251
252
253
254
255
256
257
	~vector_dist()
	{
#ifdef SE_CLASS2
		check_delete(this);
#endif
	}

Pietro Incardona's avatar
Pietro Incardona committed
258
	/*! \brief return the local size of the vector
259
	 *
Pietro Incardona's avatar
Pietro Incardona committed
260
	 * \return local size
261
262
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
263
	size_t size_local()
264
	{
Pietro Incardona's avatar
Pietro Incardona committed
265
		return g_m;
266
267
	}

incardon's avatar
incardon committed
268
269
270
271
272
	/*! \brief return the local size of the vector
	 *
	 * \return local size
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
273
	size_t size_local_with_ghost()
incardon's avatar
incardon committed
274
	{
Pietro Incardona's avatar
Pietro Incardona committed
275
		return v_pos.size();
incardon's avatar
incardon committed
276
277
	}

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

292
293
294
295
296
297
298
299
300
301
302
303
304
305
	/*! \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());
	}

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

Pietro Incardona's avatar
Pietro Incardona committed
321
322
323
324
325
326
327
328
329
330
	/*! \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
331
	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
332
333
334
335
	{
		return v_prp.template get<id>(vec_key.getKey());
	}

336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
	/*! \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
357
	inline auto getLastPos() -> decltype(v_pos.template get<0>(0))
358
	{
Pietro Incardona's avatar
Pietro Incardona committed
359
		return v_pos.template get<0>(g_m - 1);
360
361
362
363
364
365
366
367
368
369
370
	}

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

incardon's avatar
incardon committed
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
	/*! \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
389
390
		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
391

incardon's avatar
incardon committed
392
		// Processor bounding box
incardon's avatar
incardon committed
393
394
		Box<dim, St> pbox = getDecomposition().getProcessorBounds();

incardon's avatar
incardon committed
395
396
397
		// Ghost padding extension
		Ghost<dim,size_t> g_ext(0);
		cell_list.Initialize(cd_sm,pbox,pad);
incardon's avatar
incardon committed
398

incardon's avatar
incardon committed
399
		updateCellListSym(cell_list);
incardon's avatar
incardon committed
400
401

		return cell_list;
incardon's avatar
incardon committed
402
403
	}

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

		return getCellList(r_cut, g);
420
421
	}

422
423
424
425
426
427
428
429
430
431
432
433
	/*! \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%
Pietro Incardona's avatar
Pietro Incardona committed
434
		Ghost<dim,St> g = getDecomposition().getGhost();
Pietro Incardona's avatar
Pietro Incardona committed
435
		g.magnify(1.013);
436
437
438
439

		return getCellList_hilb(r_cut, g);
	}

Pietro Incardona's avatar
Pietro Incardona committed
440
441
442
443
444
445
446
447
448
	/*! \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
449
		populate_cell_list(v_pos,cell_list,g_m,CL_NON_SYMMETRIC);
Pietro Incardona's avatar
Pietro Incardona committed
450

incardon's avatar
incardon committed
451
452
453
454
455
456
457
458
459
460
461
462
		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
463
		populate_cell_list(v_pos,cell_list,g_m,CL_SYMMETRIC);
incardon's avatar
incardon committed
464

Pietro Incardona's avatar
Pietro Incardona committed
465
		cell_list.set_gm(g_m);
Pietro Incardona's avatar
Pietro Incardona committed
466
467
	}

468
469
470
471
472
473
474
475
476
477
478
479
	/*! \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
480
481
	 * \return the CellList
	 *
482
	 */
483
	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
484
485
486
	{
		CellL cell_list;

487
488
		// Division array
		size_t div[dim];
489

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

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

496
497
498
499
500
501
		cell_list.Initialize(pbox, div);

		updateCellList(cell_list);

		return cell_list;
	}
502

503
504
505
506
507
508
509
510
511
512
513
514
	/*! \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
515
516
	 * \return The Cell-list
	 *
517
518
519
520
521
522
	 */
	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
523
524
		size_t div[dim];

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

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

531
		cell_list.Initialize(pbox, div, g_m);
532

Pietro Incardona's avatar
Pietro Incardona committed
533
		updateCellList(cell_list);
Pietro Incardona's avatar
Pietro Incardona committed
534
535
536

		return cell_list;
	}
incardon's avatar
incardon committed
537

incardon's avatar
incardon committed
538
539
540
541
	/*! \brief for each particle get the symmetric verlet list
	 *
	 * \param r_cut cut-off radius
	 *
incardon's avatar
incardon committed
542
543
	 * \return the verlet list
	 *
incardon's avatar
incardon committed
544
545
546
547
	 */
	VerletList<dim,St,FAST,shift<dim,St> > getVerletSym(St r_cut)
	{
		VerletList<dim,St,FAST,shift<dim,St>> ver;
548

incardon's avatar
incardon committed
549
550
551
552
553
554
555
		// Processor bounding box
		Box<dim, St> pbox = getDecomposition().getProcessorBounds();

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

		return ver;
	}
Pietro Incardona's avatar
Pietro Incardona committed
556

incardon's avatar
incardon committed
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
588
589
590
591
592
593
	/*! \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;
	}

594
595
596
597
	/*! \brief for each particle get the verlet list
	 *
	 * \param r_cut cut-off radius
	 *
incardon's avatar
incardon committed
598
599
	 * \return a VerletList object
	 *
600
	 */
incardon's avatar
incardon committed
601
	VerletList<dim,St,FAST,shift<dim,St> > getVerlet(St r_cut)
602
603
604
605
	{
		VerletList<dim,St,FAST,shift<dim,St>> ver;

		// get the processor bounding box
Pietro Incardona's avatar
Pietro Incardona committed
606
607
608
609
610
		Box<dim, St> bt = getDecomposition().getProcessorBounds();

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

Pietro Incardona's avatar
Pietro Incardona committed
612
613
		// enlarge the box where the Verlet is defined
		bt.enlarge(g);
614

incardon's avatar
incardon committed
615
		ver.Initialize(bt,getDecomposition().getProcessorBounds(),r_cut,v_pos,g_m,VL_NON_SYMMETRIC);
616
617
618
619

		return ver;
	}

Pietro Incardona's avatar
Pietro Incardona committed
620
621
622
623
624
625
626
627
628
629
	/*! \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
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
		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);
Pietro Incardona's avatar
Pietro Incardona committed
651
652
	}

653
654
655
656
657
	/*! \brief for each particle get the verlet list
	 *
	 * \param verlet output verlet list for each particle
	 * \param r_cut cut-off radius
	 *
658
659
	 * \deprecated
	 *
660
	 */
661
	void getVerletDeprecated(openfpm::vector<openfpm::vector<size_t>> & verlet, St r_cut)
662
663
664
665
666
667
668
669
	{
		// 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
670
		St r_cut2 = r_cut * r_cut;
671
672

		// iterate the particles
673
674
675
676
677
678
679
		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
680
			Point<dim, St> p = this->getPos(key);
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697

			// 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
698
				Point<dim, St> q = this->getPos(nnp);
699
700
701
702
703
704
705
706
707
708
709

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

				// Next particle
				++NN;
			}

			// next particle
			++it_p;
		}
710
711
	}

Yaroslav's avatar
Yaroslav committed
712
713
714
715
716
717
718
719
720
721
722
	/*! \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
723
		reorder(m,getDecomposition().getGhost());
Yaroslav's avatar
Yaroslav committed
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
	}


	/*! \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
751
		Box<dim,St> pbox = getDecomposition().getProcessorBounds();
Yaroslav's avatar
Yaroslav committed
752
753
754
755
756
757
758
759
760
761
762
		// 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;
		}

763
		cell_list.Initialize(pbox,div);
Yaroslav's avatar
Yaroslav committed
764
765
766
767
768
769
770
771
772

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

		auto it = getIterator();

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

773
			cell_list.add(this->getPos(key),key.getKey());
Yaroslav's avatar
Yaroslav committed
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821

			++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);
	}

822
823
	/*! \brief It return the number of particles contained by the previous processors
	 *
Pietro Incardona's avatar
Pietro Incardona committed
824
	 * \warning It only work with the initial decomposition
825
826
827
828
829
830
831
832
833
	 *
	 * 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
834
835
	 * \return number of particles contained by the previous processors
	 *
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
	 */
	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
857

858
	/*! \brief Get an iterator that traverse domain and ghost particles
incardon's avatar
incardon committed
859
860
861
862
	 *
	 * \return an iterator
	 *
	 */
863
	vector_dist_iterator getIterator()
incardon's avatar
incardon committed
864
	{
865
		return vector_dist_iterator(0, v_pos.size());
incardon's avatar
incardon committed
866
867
	}

868
869
870
871
	/*! /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
872
873
	 * \param sz size of the grid
	 *
874
875
876
	 * \return a Grid iterator
	 *
	 */
877
	inline grid_dist_id_iterator_dec<Decomposition> getGridIterator(const size_t (&sz)[dim])
878
879
880
	{
		grid_key_dx<dim> start;
		grid_key_dx<dim> stop;
881
		for (size_t i = 0; i < dim; i++)
882
		{
883
			start.set_d(i, 0);
884
			stop.set_d(i, sz[i] - 1);
885
886
		}

Pietro Incardona's avatar
Pietro Incardona committed
887
		grid_dist_id_iterator_dec<Decomposition> it_dec(getDecomposition(), sz, start, stop);
888
889
890
		return it_dec;
	}

891
892
893
894
895
	/*! \brief Get the iterator across the position of the ghost particles
	 *
	 * \return an iterator
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
896
	vector_dist_iterator getGhostIterator() const
897
	{
898
		return vector_dist_iterator(g_m, v_pos.size());
899
900
	}

901
	/*! \brief Get an iterator that traverse the particles in the domain
902
903
904
905
	 *
	 * \return an iterator
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
906
	vector_dist_iterator getDomainIterator() const
907
	{
908
		return vector_dist_iterator(0, g_m);
909
910
	}

Pietro Incardona's avatar
Pietro Incardona committed
911
912
913
914
915
916
917
918
919
920
	/*! \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
921
922
923
924
925
	/*! \brief Get the decomposition
	 *
	 * \return
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
926
927
928
929
930
	inline Decomposition & getDecomposition()
	{
		return vector_dist_comm<dim,St,prop,Decomposition,Memory>::getDecomposition();
	}

incardon's avatar
incardon committed
931
932
933
934
935
936
937
938
939
940
	/*! \brief Get the decomposition
	 *
	 * \return
	 *
	 */
	inline const Decomposition & getDecomposition() const
	{
		return vector_dist_comm<dim,St,prop,Decomposition,Memory>::getDecomposition();
	}

Pietro Incardona's avatar
Pietro Incardona committed
941
942
943
944
945
946
947
948
949
950
951
952
953
	/*! \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
954
	{
Pietro Incardona's avatar
Pietro Incardona committed
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
		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
984
	}
incardon's avatar
incardon committed
985

Pietro Incardona's avatar
Pietro Incardona committed
986
987
988
989
990
	/*! \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
	 *
incardon's avatar
incardon committed
991
992
993
	 * \param opt_ options. It is an optional parameter.
	 *             If set to NO_CHANGE_ELEMENTS the communication has lower latencies. This option has some usage limitations, so please refere to the examples
	 *             for further explanations
incardon's avatar
incardon committed
994
	 *
Pietro Incardona's avatar
Pietro Incardona committed
995
996
	 *
	 */
incardon's avatar
incardon committed
997
	template<template<typename,typename> class op, int ... prp> inline void ghost_put(size_t opt_ = NONE)
Pietro Incardona's avatar
Pietro Incardona committed
998
	{
incardon's avatar
incardon committed
999
		this->template ghost_put_<op,prp...>(v_pos,v_prp,g_m,opt_);
Pietro Incardona's avatar
Pietro Incardona committed
1000
1001
	}

Pietro Incardona's avatar
Pietro Incardona committed
1002
1003
1004
1005
1006
1007
1008
1009
	/*! \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
	 *
	 */
1010
1011
1012
1013
	void remove(openfpm::vector<size_t> & keys, size_t start = 0)
	{
		v_pos.remove(keys, start);
		v_prp.remove(keys, start);
1014
1015

		g_m -= keys.size();
1016
1017
	}

Pietro Incardona's avatar
Pietro Incardona committed
1018
1019
	/*! \brief Remove one element from the distributed vector
	 *
Pietro Incardona's avatar
Pietro Incardona committed
1020
	 * \param key remove one element from the vector
Pietro Incardona's avatar
Pietro Incardona committed
1021
1022
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
1023
1024
1025
1026
1027
1028
1029
1030
	void remove(size_t key)
	{
		v_pos.remove(key);
		v_prp.remove(key);

		g_m--;
	}

incardon's avatar
incardon committed
1031
	/*! \brief Add the computation cost on the decomposition coming from the particles
Pietro Incardona's avatar
Pietro Incardona committed
1032
1033
	 *
	 */
incardon's avatar
incardon committed
1034
	template <typename Model=ModelLin>inline void addComputationCosts(Model md=Model())
1035
	{
incardon's avatar
incardon committed
1036
		CellDecomposer_sm<dim, St, shift<dim,St>> cdsm;
1037

Pietro Incardona's avatar
Pietro Incardona committed
1038
1039
		Decomposition & dec = getDecomposition();

1040
1041
		cdsm.setDimensions(dec.getDomain(), dec.getGrid().getSize(), 0);

Pietro Incardona's avatar
Pietro Incardona committed
1042
		for (size_t i = 0; i < getDecomposition().getNSubSubDomains(); i++)
Pietro Incardona's avatar
Pietro Incardona committed
1043
			dec.setSubSubDomainComputationCost(i, 1);
1044

tonynsyde's avatar
tonynsyde committed
1045
1046
		auto it = getDomainIterator();

1047
1048
		while (it.isNext())
		{
Pietro Incardona's avatar
Pietro Incardona committed
1049
			size_t v = cdsm.getCell(this->getPos(it.get()));
1050

incardon's avatar
incardon committed
1051
			md.addComputation(dec,*this,v,it.get().getKey());
1052
1053
1054
1055

			++it;
		}

incardon's avatar
incardon committed
1056
1057
1058
1059
		// Go throught all the sub-sub-domains and apply the model

		for (size_t i = 0 ; i < dec.getDistribution().getNSubSubDomains(); i++)
			md.applyModel(dec,i);
1060
1061
	}

incardon's avatar
incardon committed
1062
1063
	/*! \brief Output particle position and properties
	 *
1064
	 * \param out output
1065
	 * \param opt VTK_WRITER or CSV_WRITER
incardon's avatar
incardon committed
1066
	 *
1067
	 * \return true if the file has been written without error
incardon's avatar
incardon committed
1068
1069
	 *
	 */
1070
	inline bool write(std::string out, int opt = NO_GHOST | VTK_WRITER )
incardon's avatar
incardon committed
1071
1072
	{

Pietro Incardona's avatar
Pietro Incardona committed
1073
1074
1075
1076
		if ((opt & 0xFFFF0000) == CSV_WRITER)
		{
			// CSVWriter test
			CSVWriter<openfpm::vector<Point<dim,St>>, openfpm::vector<prop> > csv_writer;
incardon's avatar
incardon committed
1077

Pietro Incardona's avatar
Pietro Incardona committed
1078
			std::string output = std::to_string(out + "_" + std::to_string(v_cl.getProcessUnitID()) + std::to_string(".csv"));
Pietro Incardona's avatar
Pietro Incardona committed
1079
1080
1081
1082

			// Write the CSV
			return csv_writer.write(output,v_pos,v_prp);
		}
1083
		else if ((opt & 0xFFFF0000) == VTK_WRITER)
Pietro Incardona's avatar
Pietro Incardona committed
1084
		{
1085
			// VTKWriter for a set of points
Pietro Incardona's avatar
Pietro Incardona committed
1086
			VTKWriter<boost::mpl::pair<openfpm::vector<Point<dim,St>>, openfpm::vector<prop>>, VECTOR_POINTS> vtk_writer;
1087
			vtk_writer.add(v_pos,v_prp,g_m);
incardon's avatar
incardon committed
1088

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

1091
1092
			// Write the VTK file
			return vtk_writer.write(output);
Pietro Incardona's avatar
Pietro Incardona committed
1093
		}
Pietro Incardona's avatar
Pietro Incardona committed
1094
1095

		return false;
1096
1097
	}

Pietro Incardona's avatar
Pietro Incardona committed
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
	/*! \brief Delete the particles on the ghost
	 *
	 *
	 */
	void deleteGhost()
	{
		v_pos.resize(g_m);
		v_prp.resize(g_m);
	}

1108
1109
	/*! \brief Resize the vector (locally)
	 *
Pietro Incardona's avatar
Pietro Incardona committed
1110
	 * \warning It automatically delete the ghosts
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
	 *
	 * \param rs
	 *
	 */
	void resize(size_t rs)
	{
		deleteGhost();

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

		g_m = rs;
	}

1125
1126
1127
	/*! \brief Output particle position and properties
	 *
	 * \param out output
Pietro Incardona's avatar
Pietro Incardona committed
1128
	 * \param iteration (we can append the number at the end of the file_name)
1129
1130
1131
1132
1133
1134
1135
	 * \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
1136
1137
1138
1139
		if ((opt & 0xFFFF0000) == CSV_WRITER)
		{
			// CSVWriter test
			CSVWriter<openfpm::vector<Point<dim, St>>, openfpm::vector<prop> > csv_writer;
1140

Pietro Incardona's avatar
Pietro Incardona committed
1141
			std::string output = std::to_string(out + "_" + std::to_string(v_cl.getProcessUnitID()) + "_" + std::to_string(iteration) + std::to_string(".csv"));
1142

Pietro Incardona's avatar
Pietro Incardona committed
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
			// 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
1157
	}
1158

Pietro Incardona's avatar
Pietro Incardona committed
1159
1160
1161
	/*! \brief Get the Celllist parameters
	 *
	 * \param r_cut spacing of the cell-list
Pietro Incardona's avatar
Pietro Incardona committed
1162
	 * \param div division required for the cell-list
Pietro Incardona's avatar
Pietro Incardona committed
1163
	 * \param box where the Cell list must be defined (In general Processor domain + Ghost)
Pietro Incardona's avatar
Pietro Incardona committed
1164
	 * \param enlarge Optionally a request to make the space a littler bit larger than Processor domain + Ghost
Pietro Incardona's avatar
Pietro Incardona committed
1165
1166
1167
1168
1169
	 *        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
1170
1171
1172
1173
1174
1175
1176
1177
		// 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
1178
1179
1180

		// output the fixed domain
		box = pbox;
Pietro Incardona's avatar
Pietro Incardona committed
1181
1182
	}

Pietro Incardona's avatar
Pietro Incardona committed
1183
	/*! \brief It return the id of structure in the allocation list
1184
1185
1186
	 *
	 * \see print_alloc and SE_CLASS2
	 *
Pietro Incardona's avatar
Pietro Incardona committed
1187
1188
	 * \return the id
	 *
1189
1190
1191
1192
1193
	 */
	long int who()
	{
#ifdef SE_CLASS2
		return check_whoami(this,8);
Pietro Incardona's avatar
Pietro Incardona committed
1194
1195
#else
		return -1;
1196
1197
#endif
	}
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211

	/*! \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
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221

	/*! \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
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243

	/*! \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;
	}
incardon's avatar
incardon committed
1244
1245
1246
1247
1248
1249

	/*! \brief Get a special particle iterator able to iterate across particles using
	 *         symmetric crossing scheme
	 *
	 * \param NN Cell-List neighborhood
	 *
incardon's avatar
incardon committed
1250
1251
	 * \return Particle iterator
	 *
incardon's avatar
incardon committed
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
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
	 */
	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
1300
1301
};

Pietro Incardona's avatar
Pietro Incardona committed
1302

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