vector_dist.hpp 38.5 KB
Newer Older
Pietro Incardona's avatar
Pietro Incardona 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"
Pietro Incardona's avatar
Pietro Incardona committed
12
#include "VCluster.hpp"
Pietro Incardona's avatar
Pietro Incardona committed
13
#include "Space/Shape/Point.hpp"
Pietro Incardona's avatar
Pietro Incardona committed
14
15
16
#include "Vector/vector_dist_iterator.hpp"
#include "Space/Shape/Box.hpp"
#include "Vector/vector_dist_key.hpp"
Pietro Incardona's avatar
Pietro Incardona committed
17
18
#include "memory/PreAllocHeapMemory.hpp"
#include "memory/PtrMemory.hpp"
Pietro Incardona's avatar
Pietro Incardona committed
19
#include "NN/CellList/CellList.hpp"
20
#include "util/common.hpp"
Pietro Incardona's avatar
Pietro Incardona committed
21
#include "util/object_util.hpp"
Pietro Incardona's avatar
Pietro Incardona 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"
26
#include "Grid/grid_dist_id_iterator_dec.hpp"
Pietro Incardona's avatar
Pietro Incardona committed
27
#include "Vector/vector_dist_ofb.hpp"
28
29

#define V_SUB_UNIT_FACTOR 64
Pietro Incardona's avatar
Pietro Incardona committed
30
31
32
33

#define NO_ID false
#define ID true

Pietro Incardona's avatar
Pietro Incardona committed
34
// Perform a ghost get or a ghost put
Pietro Incardona's avatar
Pietro Incardona committed
35
36
37
#define GET	1
#define PUT 2

Pietro Incardona's avatar
Pietro Incardona committed
38
#define NO_POSITION 1
39
#define WITH_POSITION 2
Pietro Incardona's avatar
Pietro Incardona committed
40

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

Pietro Incardona's avatar
Pietro Incardona 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 ...
Pietro Incardona's avatar
Pietro Incardona committed
67
68
69
 *
 */

70
template<unsigned int dim, typename St, typename prop, typename Decomposition, typename Memory = HeapMemory>
Pietro Incardona's avatar
Pietro Incardona committed
71
72
73
74
class vector_dist
{
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;

Pietro Incardona's avatar
Pietro Incardona committed
78
79
80
	//! Space Decomposition
	Decomposition dec;

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;
Pietro Incardona's avatar
Pietro Incardona 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;
Pietro Incardona's avatar
Pietro Incardona committed
88

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

92
	// definition of the send vector for position
93
	typedef openfpm::vector<Point<dim, St>, ExtPreAlloc<Memory>, openfpm::grow_policy_identity> send_pos_vector;
94
95
96
97
98
99
100
101
102

	//////////////////////////////
	// COMMUNICATION variables
	//////////////////////////////

	//! It map the processor id with the communication request into map procedure
	openfpm::vector<size_t> p_map_req;

	//! For each near processor, outgoing particle id and shift vector
103
104
105
106
	openfpm::vector<openfpm::vector<size_t>> opart;

	//! For each near processor, particle shift vector
	openfpm::vector<openfpm::vector<size_t>> oshift;
107
108
109
110
111

	//! For each adjacent processor the size of the ghost sending buffer
	openfpm::vector<size_t> ghost_prc_sz;

	//! Sending buffer for the ghost particles properties
112
	BHeapMemory g_prp_mem;
113
114

	//! Sending buffer for the ghost particles position
115
	BHeapMemory g_pos_mem;
116
117
118
119

	//! For each adjacent processor it store the size of the receiving message in byte
	openfpm::vector<size_t> recv_sz;

120
	//! For each adjacent processor it store the received message for ghost get
121
	openfpm::vector<BHeapMemory> recv_mem_gg;
122

123
	//! For each processor it store the received message for global map
124
	openfpm::vector<BHeapMemory> recv_mem_gm;
125

126
127
128
129
130
131
132
	/*! \brief It store for each processor the position and properties vector of the particles
	 *
	 *
	 */
	struct pos_prop
	{
		//! position vector
133
		openfpm::vector<Point<dim, St>, PreAllocHeapMemory<2>, openfpm::grow_policy_identity> pos;
134
		//! properties vector
135
		openfpm::vector<prop, PreAllocHeapMemory<2>, openfpm::grow_policy_identity> prp;
136
	};
137

138
139
140
	/*! \brief Label particles for mappings
	 *
	 * \param lbl_p Particle labeled
141
	 * \param prc_sz For each processor the number of particles to send
142
143
144
	 * \param opart id of the particles to send
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
145
	template<typename obp> void labelParticleProcessor(openfpm::vector<openfpm::vector<size_t>> & lbl_p, openfpm::vector<size_t> & prc_sz, openfpm::vector<size_t> & opart)
146
	{
147
148
		// reset lbl_p
		lbl_p.resize(v_cl.getProcessingUnits());
149
		for (size_t i = 0; i < lbl_p.size(); i++)
150
151
			lbl_p.get(i).clear();

152
153
154
155
156
157
158
159
160
161
162
163
164
		// resize the label buffer
		prc_sz.resize(v_cl.getProcessingUnits());

		auto it = v_pos.getIterator();

		// Label all the particles with the processor id where they should go
		while (it.isNext())
		{
			auto key = it.get();

			// Apply the boundary conditions
			dec.applyPointBC(v_pos.get(key));

Pietro Incardona's avatar
Pietro Incardona committed
165
166
167
			size_t p_id = 0;

			// Check if the particle is inside the domain
168
			if (dec.getDomain().isInside(v_pos.get(key)) == true)
Pietro Incardona's avatar
Pietro Incardona committed
169
				p_id = dec.processorIDBC(v_pos.get(key));
170
			else
171
				p_id = obp::out(key, v_cl.getProcessUnitID());
172

173
			// Particle to move
174
175
			if (p_id != v_cl.getProcessUnitID())
			{
176
				if ((long int) p_id != -1)
177
				{
178
179
					prc_sz.get(p_id)++;lbl_p
					.get(p_id).add(key);
180
181
182
183
184
185
					opart.add(key);
				}
				else
				{
					opart.add(key);
				}
186
187
188
189
190
191
192
193
			}

			// Add processors and add size

			++it;
		}
	}

194
195
196
197
	/*! \brief Label the particles
	 *
	 * It count the number of particle to send to each processors and save its ids
	 *
198
199
200
201
202
203
	 * \param prc_sz For each processor the number of particles to send
	 * \param opart id if of the particles to send
	 * \param shift_id shift correction id
	 *
	 * \see nn_prcs::getShiftvectors()
	 *
204
	 */
205
	void labelParticlesGhost()
206
207
208
209
210
211
212
213
214
	{
		// Buffer that contain the number of elements to send for each processor
		ghost_prc_sz.clear();
		ghost_prc_sz.resize(dec.getNNProcessors());

		// Buffer that contain for each processor the id of the particle to send
		opart.clear();
		opart.resize(dec.getNNProcessors());

215
216
217
218
		// Buffer that contain for each processor the id of the shift vector
		oshift.clear();
		oshift.resize(dec.getNNProcessors());

219
		// Iterate over all particles
220
		auto it = v_pos.getIteratorTo(g_m);
221
222
223
224
225
226
		while (it.isNext())
		{
			auto key = it.get();

			// Given a particle, it return which processor require it (first id) and shift id, second id
			// For an explanation about shifts vectors please consult getShiftVector in ie_ghost
227
			const openfpm::vector<std::pair<size_t, size_t>> & vp_id = dec.template ghost_processorID_pair<typename Decomposition::lc_processor_id, typename Decomposition::shift_id>(v_pos.get(key), UNIQUE);
228

229
			for (size_t i = 0; i < vp_id.size(); i++)
230
231
232
233
234
			{
				// processor id
				size_t p_id = vp_id.get(i).first;

				// add particle to communicate
Pietro Incardona's avatar
Pietro Incardona committed
235
236
				ghost_prc_sz.get(p_id)++;
				opart.get(p_id).add(key);
237
				oshift.get(p_id).add(vp_id.get(i).second);
238
239
240
241
242
243
244
245
246
247
			}

			++it;
		}
	}

	/*! \brief Add local particles based on the boundary conditions
	 *
	 * In order to understand what this function use the following
	 *
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
	 \verbatim

	 [1,1]
	 +---------+------------------------+---------+
	 | (1,-1)  |                        | (1,1)   |
	 |   |     |    (1,0) --> 7         |   |     |
	 |   v     |                        |   v     |
	 |   6     |                        |   8     |
	 +--------------------------------------------+
	 |         |                        |         |
	 |         |                        |         |
	 |         |                        |         |
	 | (-1,0)  |                        | (1,0)   |
	 |    |    |                        |   |     |
	 |    v    |      (0,0) --> 4       |   v     |
	 |    3    |                        |   5     |
	 |         |                        |         |
	 B	|         |                        |     A   |
	 *	|         |                        |    *    |
	 |         |                        |         |
	 |         |                        |         |
	 |         |                        |         |
	 +--------------------------------------------+
	 | (-1,-1) |                        | (-1,1)  |
	 |    |    |   (-1,0) --> 1         |    |    |
	 |    v    |                        |    v    |
	 |    0    |                        |    2    |
	 +---------+------------------------+---------+


	 \endverbatim
279
280
281
282
283
284
285
286
287
288
289

	 *
	 *  The box is the domain, while all boxes at the border (so not (0,0) ) are the
	 *  ghost part at the border of the domain. If a particle A is in the position in figure
	 *  a particle B must be created. This function duplicate the particle A, if A and B are
	 *  local
	 *
	 */
	void add_loc_particles_bc()
	{
		// get the shift vectors
290
		const openfpm::vector<Point<dim, St>> & shifts = dec.getShiftVectors();
291

292
293
294
		// this map is used to check if a combination is already present
		std::unordered_map<size_t, size_t> map_cmb;

295
		// Add local particles coming from periodic boundary, the only boxes that count are the one
296
		// touching the border, filter them
297

298
299
		// The boxes touching the border of the domain are divided in groups (first vector)
		// each group contain internal ghost coming from sub-domain of the same section
300
		openfpm::vector_std<openfpm::vector_std<Box<dim, St>>>box_f;
301
302
		openfpm::vector_std<comb<dim>> box_cmb;

303
		for (size_t i = 0; i < dec.getNLocalSub(); i++)
304
305
306
		{
			size_t Nl = dec.getLocalNIGhost(i);

307
			for (size_t j = 0; j < Nl; j++)
308
			{
309
310
				// If the ghost does not come from the intersection with an out of
				// border sub-domain the combination is all zero and n_zero return dim
311
				if (dec.getLocalIGhostPos(i, j).n_zero() == dim)
312
313
					continue;

314
				// Check if we already have boxes with such combination
315
				auto it = map_cmb.find(dec.getLocalIGhostPos(i, j).lin());
316
317
318
319
				if (it == map_cmb.end())
				{
					// we do not have it
					box_f.add();
320
321
322
					box_f.last().add(dec.getLocalIGhostBox(i, j));
					box_cmb.add(dec.getLocalIGhostPos(i, j));
					map_cmb[dec.getLocalIGhostPos(i, j).lin()] = box_f.size() - 1;
323
324
325
326
				}
				else
				{
					// we have it
327
					box_f.get(it->second).add(dec.getLocalIGhostBox(i, j));
328
329
				}

330
331
332
333
334
335
336
337
			}
		}

		if (box_f.size() == 0)
			return;
		else
		{
			// Label the internal (assigned) particles
338
			auto it = v_pos.getIteratorTo(g_m);
339
340
341
342
343
344

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

				// If particles are inside these boxes
345
				for (size_t i = 0; i < box_f.size(); i++)
346
				{
347
					for (size_t j = 0; j < box_f.get(i).size(); j++)
348
					{
349
350
						if (box_f.get(i).get(j).isInside(v_pos.get(key)) == true)
						{
351
							Point<dim, St> p = v_pos.get(key);
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
							// shift
							p -= shifts.get(box_cmb.get(i).lin());

							// add this particle shifting its position
							v_pos.add(p);
							v_prp.add();
							v_prp.last() = v_prp.get(key);

							// boxes in one group can be overlapping
							// we do not have to search for the other
							// boxes otherwise we will have duplicate particles
							//
							// A small note overlap of boxes across groups is fine
							// (and needed) because each group has different shift
							// producing non overlapping particles
							//
							break;
						}
370
371
372
373
374
375
376
377
378
379
380
381
382
383
					}
				}

				++it;
			}
		}
	}

	/*! \brief This function fill the send buffer for the particle position after the particles has been label with labelParticles
	 *
	 * \param g_pos_send Send buffer to fill
	 * \param prAlloc_pos Memory object for the send buffer
	 *
	 */
384
	void fill_send_ghost_pos_buf(openfpm::vector<send_pos_vector> & g_pos_send, ExtPreAlloc<Memory> * prAlloc_pos)
385
386
	{
		// get the shift vectors
387
		const openfpm::vector<Point<dim, St>> & shifts = dec.getShiftVectors();
388
389
390

		// create a number of send buffers equal to the near processors
		g_pos_send.resize(ghost_prc_sz.size());
391
		for (size_t i = 0; i < g_pos_send.size(); i++)
392
393
394
395
396
397
398
399
400
		{
			// set the preallocated memory to ensure contiguity
			g_pos_send.get(i).setMemory(*prAlloc_pos);

			// resize the sending vector (No allocation is produced)
			g_pos_send.get(i).resize(ghost_prc_sz.get(i));
		}

		// Fill the send buffer
401
		for (size_t i = 0; i < opart.size(); i++)
402
		{
403
			for (size_t j = 0; j < opart.get(i).size(); j++)
404
			{
405
				Point<dim, St> s = v_pos.get(opart.get(i).get(j));
406
				s -= shifts.get(oshift.get(i).get(j));
407
				g_pos_send.get(i).set(j, s);
408
409
410
411
412
413
414
415
416
417
418
419
420
421
			}
		}
	}

	/*! \brief This function fill the send buffer for properties after the particles has been label with labelParticles
	 *
	 * \tparam send_vector type used to send data
	 * \tparam prp_object object containing only the properties to send
	 * \tparam prp set of properties to send
	 *
	 * \param g_send_prp Send buffer to fill
	 * \param prAlloc_prp Memory object for the send buffer
	 *
	 */
422
	template<typename send_vector, typename prp_object, int ... prp> void fill_send_ghost_prp_buf(openfpm::vector<send_vector> & g_send_prp, ExtPreAlloc<Memory> * prAlloc_prp)
423
424
425
	{
		// create a number of send buffers equal to the near processors
		g_send_prp.resize(ghost_prc_sz.size());
426
		for (size_t i = 0; i < g_send_prp.size(); i++)
427
428
429
430
431
432
433
434
435
		{
			// set the preallocated memory to ensure contiguity
			g_send_prp.get(i).setMemory(*prAlloc_prp);

			// resize the sending vector (No allocation is produced)
			g_send_prp.get(i).resize(ghost_prc_sz.get(i));
		}

		// Fill the send buffer
436
		for (size_t i = 0; i < opart.size(); i++)
437
		{
438
			for (size_t j = 0; j < opart.get(i).size(); j++)
439
440
			{
				// source object type
441
				typedef encapc<1, prop, typename openfpm::vector<prop>::memory_conf> encap_src;
442
				// destination object type
443
				typedef encapc<1, prp_object, typename openfpm::vector<prp_object>::memory_conf> encap_dst;
444
445

				// Copy only the selected properties
446
				object_si_d<encap_src, encap_dst, OBJ_ENCAP, prp...>(v_prp.get(opart.get(i).get(j)), g_send_prp.get(i).get(j));
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
			}
		}
	}

	/*! \brief allocate and fill the send buffer for the map function
	 *
	 * \param prc_r List of processor rank involved in the send
	 * \param prc_r_sz For each processor in the list the size of the message to send
	 * \param pb send buffer
	 *
	 */
	void fill_send_map_buf(openfpm::vector<size_t> & prc_r, openfpm::vector<size_t> & prc_sz_r, openfpm::vector<pos_prop> & pb)
	{
		pb.resize(prc_r.size());

462
		for (size_t i = 0; i < prc_r.size(); i++)
463
464
		{
			// Create the size required to store the particles position and properties to communicate
465
466
			size_t s1 = openfpm::vector<Point<dim, St>, HeapMemory, openfpm::grow_policy_identity>::calculateMem(prc_sz_r.get(i), 0);
			size_t s2 = openfpm::vector<prop, HeapMemory, openfpm::grow_policy_identity>::calculateMem(prc_sz_r.get(i), 0);
467
468

			// Preallocate the memory
469
			size_t sz[2] = { s1, s2 };
470
471
472
473
474
475
476
477
478
479
480
481
482
			PreAllocHeapMemory<2> * mem = new PreAllocHeapMemory<2>(sz);

			// Set the memory allocator
			pb.get(i).pos.setMemory(*mem);
			pb.get(i).prp.setMemory(*mem);

			// set the size and allocate, using mem warant that pos and prp is contiguous
			pb.get(i).pos.resize(prc_sz_r.get(i));
			pb.get(i).prp.resize(prc_sz_r.get(i));
		}

		// Run through all the particles and fill the sending buffer

483
		for (size_t i = 0; i < opart.size(); i++)
484
485
486
487
488
489
490
491
492
		{
			auto it = opart.get(i).getIterator();
			size_t lbl = p_map_req.get(i);

			while (it.isNext())
			{
				size_t key = it.get();
				size_t id = opart.get(i).get(key);

493
494
				pb.get(lbl).pos.set(key, v_pos.get(id));
				pb.get(lbl).prp.set(key, v_prp.get(id));
495
496

				++it;
497
498
499
500
501
502
503
504
505
506
507
			}
		}
	}

	/*! \brief This function process the receiced data for the properties and populate the ghost
	 *
	 * \tparam send_vector type used to send data
	 * \tparam prp_object object containing only the properties to send
	 * \tparam prp set of properties to send
	 *
	 */
508
	template<typename send_vector, typename prp_object, int ... prp> void process_received_ghost_prp()
509
510
511
512
513
	{
		// Mark the ghost part
		g_m = v_prp.size();

		// Process the received data (recv_mem_gg != 0 if you have data)
514
		for (size_t i = 0; i < dec.getNNProcessors() && recv_mem_gg.size() != 0; i++)
515
516
517
518
519
		{
			// calculate the number of received elements
			size_t n_ele = recv_sz.get(i) / sizeof(prp_object);

			// add the received particles to the vector
520
			PtrMemory * ptr1 = new PtrMemory(recv_mem_gg.get(i).getPointer(), recv_sz.get(i));
521
522

			// create vector representation to a piece of memory already allocated
523
			openfpm::vector<prp_object, PtrMemory, openfpm::grow_policy_identity> v2;
524
525
526
527
528
529
530

			v2.setMemory(*ptr1);

			// resize with the number of elements
			v2.resize(n_ele);

			// Add the ghost particle
531
			v_prp.template add_prp<prp_object, PtrMemory, openfpm::grow_policy_identity, prp...>(v2);
532
533
534
535
536
537
		}
	}

	/*! \brief This function process the received data for the properties and populate the ghost
	 *
	 */
538
	void process_received_ghost_pos()
539
540
	{
		// Process the received data (recv_mem_gg != 0 if you have data)
541
		for (size_t i = 0; i < dec.getNNProcessors() && recv_mem_gg.size() != 0; i++)
542
543
		{
			// calculate the number of received elements
544
			size_t n_ele = recv_sz.get(i) / sizeof(Point<dim, St> );
545
546

			// add the received particles to the vector
547
			PtrMemory * ptr1 = new PtrMemory(recv_mem_gg.get(i).getPointer(), recv_sz.get(i));
548
549
550

			// create vector representation to a piece of memory already allocated

551
			openfpm::vector<Point<dim, St>, PtrMemory, openfpm::grow_policy_identity> v2;
552
553
554
555
556
557
558

			v2.setMemory(*ptr1);

			// resize with the number of elements
			v2.resize(n_ele);

			// Add the ghost particle
559
			v_pos.template add<PtrMemory, openfpm::grow_policy_identity>(v2);
560
561
562
		}
	}

563
564
565
566
567
568
569
570
571
	/*! \brief Process the received particles
	 *
	 * \param list of the out-going particles
	 *
	 */
	void process_received_map(openfpm::vector<size_t> & out_part)
	{
		size_t o_p_id = 0;

572
		for (size_t i = 0; i < recv_mem_gm.size(); i++)
573
574
575
		{
			// Get the number of elements

576
			size_t n_ele = recv_mem_gm.get(i).size() / (sizeof(Point<dim, St> ) + sizeof(prop));
577
578

			// Pointer of the received positions for each near processor
579
			void * ptr_pos = (unsigned char *) recv_mem_gm.get(i).getPointer();
580
			// Pointer of the received properties for each near processor
581
			void * ptr_prp = (unsigned char *) recv_mem_gm.get(i).getPointer() + n_ele * sizeof(Point<dim, St> );
582

583
584
			PtrMemory * ptr1 = new PtrMemory(ptr_pos, n_ele * sizeof(Point<dim, St> ));
			PtrMemory * ptr2 = new PtrMemory(ptr_prp, n_ele * sizeof(prop));
585
586
587

			// create vector representation to a piece of memory already allocated

588
589
			openfpm::vector<Point<dim, St>, PtrMemory, openfpm::grow_policy_identity> vpos;
			openfpm::vector<prop, PtrMemory, openfpm::grow_policy_identity> vprp;
590
591
592
593
594
595
596
597
598
599

			vpos.setMemory(*ptr1);
			vprp.setMemory(*ptr2);

			vpos.resize(n_ele);
			vprp.resize(n_ele);

			// Add the received particles to v_pos and v_prp

			size_t j = 0;
600
			for (; j < vpos.size() && o_p_id < out_part.size(); j++, o_p_id++)
601
			{
602
603
				v_pos.set(out_part.get(o_p_id), vpos.get(j));
				v_prp.set(out_part.get(o_p_id), vprp.get(j));
604
605
			}

606
			for (; j < vpos.size(); j++)
607
608
			{
				v_pos.add();
609
				v_pos.set(v_pos.size() - 1, vpos.get(j));
610
				v_prp.add();
611
				v_prp.set(v_prp.size() - 1, vprp.get(j));
612
613
614
615
616
			}
		}

		// remove the (out-going particles) in the vector

617
618
		v_pos.remove(out_part, o_p_id);
		v_prp.remove(out_part, o_p_id);
619
620
	}

621
622
623
624
625
626
627
628
629
630
	/*! \brief Calculate send buffers total size and allocation
	 *
	 * \tparam prp_object object containing only the properties to send
	 *
	 * \param size_byte_prp total size for the property buffer
	 * \param size_byte_pos total size for the position buffer
	 * \param pap_prp allocation sequence for the property buffer
	 * \param pap_pos allocation sequence for the position buffer
	 *
	 */
631
	template<typename prp_object> void calc_send_ghost_buf(size_t & size_byte_prp, size_t & size_byte_pos, std::vector<size_t> & pap_prp, std::vector<size_t> & pap_pos)
632
633
	{
		// Calculate the total size required for the sending buffer
634
		for (size_t i = 0; i < ghost_prc_sz.size(); i++)
635
		{
636
			size_t alloc_ele = openfpm::vector<prp_object, HeapMemory, openfpm::grow_policy_identity>::calculateMem(ghost_prc_sz.get(i), 0);
637
638
639
			pap_prp.push_back(alloc_ele);
			size_byte_prp += alloc_ele;

640
			alloc_ele = openfpm::vector<Point<dim, St>, HeapMemory, openfpm::grow_policy_identity>::calculateMem(ghost_prc_sz.get(i), 0);
641
642
643
644
645
			pap_pos.push_back(alloc_ele);
			size_byte_pos += alloc_ele;
		}
	}

Pietro Incardona's avatar
Pietro Incardona committed
646
647
648
649
public:

	/*! \brief Constructor
	 *
650
651
	 * \param np number of elements
	 * \param box domain where the vector of elements live
652
	 * \param boundary conditions
653
	 * \param g Ghost margins
Pietro Incardona's avatar
Pietro Incardona committed
654
655
	 *
	 */
656
657
	vector_dist(size_t np, Box<dim, St> box, const size_t (&bc)[dim], const Ghost<dim, St> & g) :
			dec(*global_v_cluster), v_cl(*global_v_cluster)
Pietro Incardona's avatar
Pietro Incardona committed
658
	{
659
660
661
662
#ifdef SE_CLASS2
		check_new(this,8,VECTOR_DIST_EVENT,4);
#endif

663
		// convert to a local number of elements
664
665
666
667
668
669
670
671
		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++;
672

Pietro Incardona's avatar
Pietro Incardona committed
673
		// resize the position vector
674
		v_pos.resize(p_np);
Pietro Incardona's avatar
Pietro Incardona committed
675
676

		// resize the properties vector
677
678
679
		v_prp.resize(p_np);

		g_m = p_np;
Pietro Incardona's avatar
Pietro Incardona committed
680
681
682
683
684

		// Create a valid decomposition of the space
		// Get the number of processor and calculate the number of sub-domain
		// for decomposition
		size_t n_proc = v_cl.getProcessingUnits();
685
		size_t n_sub = n_proc * getDefaultNsubsub();
Pietro Incardona's avatar
Pietro Incardona committed
686
687
688

		// Calculate the maximum number (before merging) of sub-domain on
		// each dimension
689
		size_t div[dim];
690
691
692
693
		for (size_t i = 0; i < dim; i++)
		{
			div[i] = openfpm::math::round_big_2(pow(n_sub, 1.0 / dim));
		}
Pietro Incardona's avatar
Pietro Incardona committed
694
695

		// Create the sub-domains
696
		dec.setParameters(div, box, bc, g);
Pietro Incardona's avatar
Pietro Incardona committed
697
		dec.decompose();
Pietro Incardona's avatar
Pietro Incardona committed
698
699
	}

700
701
702
703
704
705
706
	~vector_dist()
	{
#ifdef SE_CLASS2
		check_delete(this);
#endif
	}

707
708
709
710
711
712
713
	/*! \brief Get the number of minimum sub-domain
	 *
	 * \return minimum number
	 *
	 */
	static size_t getDefaultNsubsub()
	{
714
		return V_SUB_UNIT_FACTOR;
715
716
	}

Pietro Incardona's avatar
Pietro Incardona committed
717
718
719
720
721
722
723
	/*! \brief return the local size of the vector
	 *
	 * \return local size
	 *
	 */
	size_t size_local()
	{
Pietro Incardona's avatar
Pietro Incardona committed
724
		return g_m;
Pietro Incardona's avatar
Pietro Incardona committed
725
726
	}

727
	/*! \brief Get the position of an element
Pietro Incardona's avatar
Pietro Incardona committed
728
	 *
729
730
731
732
733
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \param vec_key element
	 *
	 * \return the position of the element in space
Pietro Incardona's avatar
Pietro Incardona committed
734
735
	 *
	 */
736
	template<unsigned int id> inline auto getPos(vect_dist_key_dx vec_key) -> decltype(v_pos.template get<id>(vec_key.getKey()))
Pietro Incardona's avatar
Pietro Incardona committed
737
	{
738
		return v_pos.template get<id>(vec_key.getKey());
Pietro Incardona's avatar
Pietro Incardona committed
739
740
	}

741
	/*! \brief Get the property of an element
Pietro Incardona's avatar
Pietro Incardona committed
742
	 *
743
744
745
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \tparam id property id
Pietro Incardona's avatar
Pietro Incardona committed
746
747
	 * \param vec_key vector element
	 *
748
749
	 * \return return the selected property of the vector element
	 *
Pietro Incardona's avatar
Pietro Incardona committed
750
	 */
751
	template<unsigned int id> inline auto getProp(vect_dist_key_dx vec_key) -> decltype(v_prp.template get<id>(vec_key.getKey()))
Pietro Incardona's avatar
Pietro Incardona committed
752
	{
753
		return v_prp.template get<id>(vec_key.getKey());
Pietro Incardona's avatar
Pietro Incardona committed
754
755
	}

756
	/*! \brief It move all the particles that does not belong to the local processor to the respective processor
Pietro Incardona's avatar
Pietro Incardona committed
757
758
	 *
	 * \tparam out of bound policy it specify what to do when the particles are detected out of bound
759
760
761
	 *
	 * 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
762
	 * contain non local particles
Pietro Incardona's avatar
Pietro Incardona committed
763
764
	 *
	 */
765
	template<typename obp = KillParticle> void map()
Pietro Incardona's avatar
Pietro Incardona committed
766
	{
Pietro Incardona's avatar
Pietro Incardona committed
767
		// outgoing particles-id
768
		openfpm::vector<size_t> out_part;
Pietro Incardona's avatar
Pietro Incardona committed
769
770
771
772
773

		// Processor communication size
		openfpm::vector<size_t> prc_sz(v_cl.getProcessingUnits());

		// It contain the list of the processors this processor should to communicate with
Pietro Incardona's avatar
Pietro Incardona committed
774
775
		openfpm::vector<size_t> p_list;

776
777
778
779
780
		// map completely reset the ghost part
		v_pos.resize(g_m);
		v_prp.resize(g_m);

		// Contain the processor id of each particle (basically where they have to go)
781
		labelParticleProcessor<obp>(opart, prc_sz, out_part);
Pietro Incardona's avatar
Pietro Incardona committed
782

783
784
		// Calculate the sending buffer size for each processor, put this information in
		// a contiguous buffer
785
		p_map_req.resize(v_cl.getProcessingUnits());
Pietro Incardona's avatar
Pietro Incardona committed
786
787
788
		openfpm::vector<size_t> prc_sz_r;
		openfpm::vector<size_t> prc_r;

789
		for (size_t i = 0; i < v_cl.getProcessingUnits(); i++)
Pietro Incardona's avatar
Pietro Incardona committed
790
		{
791
			if (prc_sz.get(i) != 0)
Pietro Incardona's avatar
Pietro Incardona committed
792
			{
793
				p_map_req.get(i) = prc_r.size();
Pietro Incardona's avatar
Pietro Incardona committed
794
795
796
797
798
				prc_r.add(i);
				prc_sz_r.add(prc_sz.get(i));
			}
		}

799
		// Allocate the send buffers
Pietro Incardona's avatar
Pietro Incardona committed
800

801
		openfpm::vector<pos_prop> pb;
Pietro Incardona's avatar
Pietro Incardona committed
802

803
		// fill the send buffers
804
		fill_send_map_buf(prc_r, prc_sz_r, pb);
Pietro Incardona's avatar
Pietro Incardona committed
805
806
807

		// Create the set of pointers
		openfpm::vector<void *> ptr(prc_r.size());
808
		for (size_t i = 0; i < prc_r.size(); i++)
Pietro Incardona's avatar
Pietro Incardona committed
809
810
811
812
813
		{
			ptr.get(i) = pb.get(i).pos.getPointer();
		}

		// convert the particle number to buffer size
814
		for (size_t i = 0; i < prc_sz_r.size(); i++)
Pietro Incardona's avatar
Pietro Incardona committed
815
		{
816
			prc_sz_r.get(i) = prc_sz_r.get(i) * (sizeof(prop) + sizeof(Point<dim, St> ));
Pietro Incardona's avatar
Pietro Incardona committed
817
818
819
820
		}

		// Send and receive the particles

821
		recv_mem_gm.clear();
822
		v_cl.sendrecvMultipleMessagesNBX(prc_sz_r.size(), (size_t *) prc_sz_r.getPointer(), (size_t *) prc_r.getPointer(), (void **) ptr.getPointer(), vector_dist::message_alloc_map, this, NEED_ALL_SIZE);
Pietro Incardona's avatar
Pietro Incardona committed
823

824
		// Process the incoming particles
Pietro Incardona's avatar
Pietro Incardona committed
825

826
		process_received_map(out_part);
Pietro Incardona's avatar
Pietro Incardona committed
827

828
		// mark the ghost part
Pietro Incardona's avatar
Pietro Incardona committed
829

830
831
		g_m = v_pos.size();
	}
Pietro Incardona's avatar
Pietro Incardona committed
832

833
	/*! \brief It synchronize the properties and position of the ghost particles
Pietro Incardona's avatar
Pietro Incardona committed
834
	 *
835
836
	 * \tparam prp list of properties to get synchronize
	 * \param opt options WITH_POSITION, it send also the positional information of the particles
Pietro Incardona's avatar
Pietro Incardona committed
837
838
	 *
	 */
839
	template<int ... prp> void ghost_get(size_t opt = WITH_POSITION)
Pietro Incardona's avatar
Pietro Incardona committed
840
	{
Pietro Incardona's avatar
Pietro Incardona committed
841
842
843
844
		// Unload receive buffer
		for (size_t i = 0 ; i < recv_mem_gg.size() ; i++)
			recv_sz.get(i) = 0;

845
		// Sending property object
846
		typedef object<typename object_creator<typename prop::type, prp...>::type> prp_object;
847
848

		// send vector for each processor
849
		typedef openfpm::vector<prp_object, ExtPreAlloc<Memory>, openfpm::grow_policy_identity> send_vector;
850

851
852
853
		// reset the ghost part
		v_pos.resize(g_m);
		v_prp.resize(g_m);
Pietro Incardona's avatar
Pietro Incardona committed
854

855
		// Label all the particles
856
		labelParticlesGhost();
Pietro Incardona's avatar
Pietro Incardona committed
857

858
		// Calculate memory and allocation for the send buffers
Pietro Incardona's avatar
Pietro Incardona committed
859

860
		// Total size
Pietro Incardona's avatar
Pietro Incardona committed
861
862
		size_t size_byte_prp = 0;
		size_t size_byte_pos = 0;
Pietro Incardona's avatar
Pietro Incardona committed
863

864
		// allocation patterns for property and position send buffer
Pietro Incardona's avatar
Pietro Incardona committed
865
866
		std::vector<size_t> pap_prp;
		std::vector<size_t> pap_pos;
Pietro Incardona's avatar
Pietro Incardona committed
867

868
		calc_send_ghost_buf<prp_object>(size_byte_prp, size_byte_pos, pap_prp, pap_pos);
Pietro Incardona's avatar
Pietro Incardona committed
869

870
		// Create memory for the send buffer
Pietro Incardona's avatar
Pietro Incardona committed
871

Pietro Incardona's avatar
Pietro Incardona committed
872
		g_prp_mem.resize(size_byte_prp);
873
874
		if (opt != NO_POSITION)
			g_pos_mem.resize(size_byte_pos);
Pietro Incardona's avatar
Pietro Incardona committed
875

876
		// Create and fill send buffer for particle properties
Pietro Incardona's avatar
Pietro Incardona committed
877

878
		ExtPreAlloc<Memory> * prAlloc_prp = new ExtPreAlloc<Memory>(pap_prp, g_prp_mem);
Pietro Incardona's avatar
Pietro Incardona committed
879
		openfpm::vector<send_vector> g_send_prp;
880
		fill_send_ghost_prp_buf<send_vector, prp_object, prp...>(g_send_prp, prAlloc_prp);
Pietro Incardona's avatar
Pietro Incardona committed
881

882
		// Create and fill the send buffer for the particle position
Pietro Incardona's avatar
Pietro Incardona committed
883

884
		ExtPreAlloc<Memory> * prAlloc_pos;
Pietro Incardona's avatar
Pietro Incardona committed
885
		openfpm::vector<send_pos_vector> g_pos_send;
Pietro Incardona's avatar
Pietro Incardona committed
886
887
		if (opt != NO_POSITION)
		{
888
889
			prAlloc_pos = new ExtPreAlloc<Memory>(pap_pos, g_pos_mem);
			fill_send_ghost_pos_buf(g_pos_send, prAlloc_pos);
Pietro Incardona's avatar
Pietro Incardona committed
890
891
		}

892
		// Create processor list
Pietro Incardona's avatar
Pietro Incardona committed
893
		openfpm::vector<size_t> prc;
894
		for (size_t i = 0; i < opart.size(); i++)
Pietro Incardona's avatar
Pietro Incardona committed
895
			prc.add(dec.IDtoProc(i));
Pietro Incardona's avatar
Pietro Incardona committed
896

897
		// Send/receive the particle properties information
898
899
		v_cl.sendrecvMultipleMessagesNBX(prc, g_send_prp, msg_alloc_ghost_get, this);
		process_received_ghost_prp<send_vector, prp_object, prp...>();
900
901
902

		if (opt != NO_POSITION)
		{
903
			// Send/receive the particle properties information
904
			v_cl.sendrecvMultipleMessagesNBX(prc, g_pos_send, msg_alloc_ghost_get, this);
905
			process_received_ghost_pos();
906
907
		}

908
909
		add_loc_particles_bc();
	}
910

911
	/*! \brief Call-back to allocate buffer to receive incoming elements (particles)
912
913
914
915
916
917
918
919
920
921
	 *
	 * \param msg_i message size required to receive from i
	 * \param total_msg message size to receive from all the processors
	 * \param total_p the total number of processor want to communicate with you
	 * \param i processor id
	 * \param ri request id (it is an id that goes from 0 to total_p, and is unique
	 *           every time message_alloc is called)
	 * \param ptr void pointer parameter for additional data to pass to the call-back
	 *
	 */
922
	static void * msg_alloc_ghost_get(size_t msg_i, size_t total_msg, size_t total_p, size_t i, size_t ri, void * ptr)
923
	{
924
		vector_dist<dim, St, prop, Decomposition, Memory> * v = static_cast<vector_dist<dim, St, prop, Decomposition, Memory> *>(ptr);
925
926
927
928
929
930

		v->recv_sz.resize(v->dec.getNNProcessors());
		v->recv_mem_gg.resize(v->dec.getNNProcessors());

		// Get the local processor id
		size_t lc_id = v->dec.ProctoID(i);
Pietro Incardona's avatar
Pietro Incardona committed
931

932
933
934
		// resize the receive buffer
		v->recv_mem_gg.get(lc_id).resize(msg_i);
		v->recv_sz.get(lc_id) = msg_i;
Pietro Incardona's avatar
Pietro Incardona committed
935

936
		return v->recv_mem_gg.get(lc_id).getPointer();
Pietro Incardona's avatar
Pietro Incardona committed
937
938
	}

939
	/*! \brief Call-back to allocate buffer to receive incoming elements (particles)
Pietro Incardona's avatar
Pietro Incardona committed
940
	 *
941
942
943
	 * \param msg_i size required to receive the message from i
	 * \param total_msg total size to receive from all the processors
	 * \param total_p the total number of processor that want to communicate with you
Pietro Incardona's avatar
Pietro Incardona committed
944
	 * \param i processor id
Pietro Incardona's avatar
Pietro Incardona committed
945
946
	 * \param ri request id (it is an id that goes from 0 to total_p, and is unique
	 *           every time message_alloc is called)
Pietro Incardona's avatar
Pietro Incardona committed
947
948
	 * \param ptr a pointer to the vector_dist structure
	 *
949
	 * \return the pointer where to store the message for the processor i
Pietro Incardona's avatar
Pietro Incardona committed
950
951
	 *
	 */
952
	static void * message_alloc_map(size_t msg_i, size_t total_msg, size_t total_p, size_t i, size_t ri, void * ptr)
Pietro Incardona's avatar
Pietro Incardona committed
953
954
	{
		// cast the pointer
955
		vector_dist<dim, St, prop, Decomposition, Memory> * vd = static_cast<vector_dist<dim, St, prop, Decomposition, Memory> *>(ptr);
Pietro Incardona's avatar
Pietro Incardona committed
956

957
958
		vd->recv_mem_gm.resize(vd->v_cl.getProcessingUnits());
		vd->recv_mem_gm.get(i).resize(msg_i);
Pietro Incardona's avatar
Pietro Incardona committed
959

960
		return vd->recv_mem_gm.get(i).getPointer();
Pietro Incardona's avatar
Pietro Incardona committed
961
962
	}

963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
	/*! \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
	 *
	 */
	template<unsigned int id> inline auto getLastPos() -> decltype(v_pos.template get<id>(0))
	{
986
		return v_pos.template get<id>(g_m - 1);
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
	}

	/*! \brief Get the property of the last element
	 *
	 * see the vector_dist iterator usage to get an element key
	 *
	 * \tparam id property id
	 * \param vec_key vector element
	 *
	 * \return return the selected property of the vector element
	 *
	 */
	template<unsigned int id> inline auto getLastProp() -> decltype(v_prp.template get<id>(0))
	{
1001
		return v_prp.template get<id>(g_m - 1);
1002
1003
	}

Pietro Incardona's avatar
Pietro Incardona committed
1004
1005
1006
1007
	/*! \brief Construct a cell list starting from the stored particles
	 *
	 * \tparam CellL CellList type to construct
	 *
1008
1009
1010
1011
	 * \param r_cut interation radius, or size of each cell
	 *
	 * \return the Cell list
	 *
Pietro Incardona's avatar
Pietro Incardona committed
1012
	 */
1013
	template<typename CellL = CellList<dim, St, FAST, shift<dim, St> > > CellL getCellList(St r_cut)
Pietro Incardona's avatar