grid_dist_id.hpp 30.3 KB
Newer Older
incardon's avatar
incardon committed
1
2
3
4
#ifndef COM_UNIT_HPP
#define COM_UNIT_HPP

#include <vector>
5
#include <unordered_map>
incardon's avatar
incardon committed
6
7
8
#include "Grid/map_grid.hpp"
#include "VCluster.hpp"
#include "Space/SpaceBox.hpp"
9
#include "util/mathutil.hpp"
incardon's avatar
incardon committed
10
11
#include "grid_dist_id_iterator.hpp"
#include "grid_dist_key.hpp"
12
#include "NN/CellList/CellDecomposer.hpp"
13
14
15
#include "util/object_util.hpp"
#include "memory/ExtPreAlloc.hpp"
#include "VTKWriter.hpp"
16
17
#include "Packer.hpp"
#include "Unpacker.hpp"
incardon's avatar
incardon committed
18
19
20
21
22
23

#define SUB_UNIT_FACTOR 64


/*! \brief This is a distributed grid
 *
24
25
26
 * Implementation of a distributed grid with decomposition on the ids.
 * A distributed grid is a grid distributed across processors.
 * The decomposition is performed on the ids of the elements
incardon's avatar
incardon committed
27
28
29
 *
 *
 * \param dim Dimensionality of the grid
30
31
 * \param St Type of space where the grid is living
 * \param T object the grid is storing
incardon's avatar
incardon committed
32
33
34
35
36
 * \param Decomposition Class that decompose the grid for example CartDecomposition
 * \param Mem Is the allocator
 * \param device type of base structure is going to store the data
 *
 */
37
template<unsigned int dim, typename St, typename T, typename Decomposition,typename Memory=HeapMemory , typename device_grid=grid_cpu<dim,T> >
incardon's avatar
incardon committed
38
39
class grid_dist_id
{
40
41
42
	// Domain
	Box<dim,St> domain;

incardon's avatar
incardon committed
43
	// Ghost expansion
44
	Ghost<dim,St> ghost;
incardon's avatar
incardon committed
45
46
47
48
49
50
51

	//! Local grids
	Vcluster_object_array<device_grid> loc_grid;

	//! Space Decomposition
	Decomposition dec;

incardon's avatar
incardon committed
52
53
54
	//! Extension of each grid: Domain and ghost + domain
	openfpm::vector<GBoxes<device_grid::dims>> gdb_ext;

incardon's avatar
incardon committed
55
56
57
	//! Size of the grid on each dimension
	size_t g_sz[dim];

58
59
	//! Structure that divide the space into cells
	CellDecomposer_sm<dim,St> cd_sm;
incardon's avatar
incardon committed
60

61
	//! Communicator class
incardon's avatar
incardon committed
62
63
	Vcluster & v_cl;

64
65
66
	//! It map a global ghost id (g_id) to the external ghost box information
	std::unordered_map<size_t,size_t> g_id_to_external_ghost_box;

incardon's avatar
incardon committed
67
68
	/*! \brief Get the grid size
	 *
incardon's avatar
incardon committed
69
	 * Given a domain, the resolution of the grid on it and another spaceBox contained in the domain
incardon's avatar
incardon committed
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
	 * it give the size on all directions of the local grid
	 *
	 * \param sp SpaceBox enclosing the local grid
	 * \param domain Space box enclosing the physical domain or part of it
	 * \param v_size grid size on this physical domain
	 *
	 * \return An std::vector representing the local grid on each dimension
	 *
	 */
	std::vector<size_t> getGridSize(SpaceBox<dim,typename Decomposition::domain_type> & sp, Box<dim,typename Decomposition::domain_type> & domain, size_t (& v_size)[dim])
	{
		std::vector<size_t> tmp;
		for (size_t d = 0 ; d < dim ; d++)
		{
			//! Get the grid size compared to the domain space and its resolution
			typename Decomposition::domain_type dim_sz = (sp.getHigh(d) - sp.getLow(d)) / ((domain.getHigh(d) - domain.getLow(d)) / v_size[d]) + 0.5;

			// push the size of the local grid
			tmp.push_back(dim_sz);
		}
		return tmp;
	}

	/*! \brief Get the grid size
	 *
	 * Get the grid size, given a spaceBox
	 * it give the size on all directions of the local grid
	 *
	 * \param sp SpaceBox enclosing the local grid
	 * \param sz array to fill with the local grid size on each dimension
	 *
	 */
	void getGridSize(SpaceBox<dim,size_t> & sp, size_t (& v_size)[dim])
	{
		for (size_t d = 0 ; d < dim ; d++)
		{
			// push the size of the local grid
incardon's avatar
incardon committed
107
			v_size[d] = sp.getHigh(d) - sp.getLow(d);
incardon's avatar
incardon committed
108
109
110
		}
	}

111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
	// Receiving size
	openfpm::vector<size_t> recv_sz;

	// Receiving buffer for particles ghost get
	openfpm::vector<HeapMemory> recv_mem_gg;

	/*! \brief Call-back to allocate buffer to receive incoming objects (external ghost boxes)
	 *
	 * \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
	 *
	 */
	static void * msg_alloc_external_box(size_t msg_i ,size_t total_msg, size_t total_p, size_t i, size_t ri, void * ptr)
	{
		grid_dist_id<dim,St,T,Decomposition,Memory,device_grid> * g = static_cast<grid_dist_id<dim,St,T,Decomposition,Memory,device_grid> *>(ptr);

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

		// Get the local processor id
		size_t lc_id = g->dec.ProctoID(i);

		// resize the receive buffer
		g->recv_mem_gg.get(lc_id).resize(msg_i);
		g->recv_sz.get(lc_id) = msg_i;

		return g->recv_mem_gg.get(lc_id).getPointer();
	}

145
	/*! \brief Create per-processor internal ghost box list in grid units
146
147
	 *
	 * It also create g_id_to_external_ghost_box
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
	 *
	 */
	void create_ig_box()
	{
		// Get the grid info
		auto g = cd_sm.getGrid();

		if (init_i_g_box == true)	return;

		// Get the number of near processors
		for (size_t i = 0 ; i < dec.getNNProcessors() ; i++)
		{
			ig_box.add();
			auto&& pib = ig_box.last();

			pib.prc = dec.IDtoProc(i);
			for (size_t j = 0 ; j < dec.getProcessorNIGhost(i) ; j++)
			{
				// Get the internal ghost boxes and transform into grid units
incardon's avatar
incardon committed
167
168
169
170
171
172
				::Box<dim,St> ib_dom = dec.getProcessorIGhostBox(i,j);
				::Box<dim,size_t> ib = cd_sm.convertDomainSpaceIntoGridUnits(ib_dom);

				// Check if ib is valid if not it mean that the internal ghost does not contain information so skip it
				if (ib.isValid() == false)
					continue;
173
174
175
176

				// save the box and the sub-domain id (it is calculated as the linearization of P1)
				::Box<dim,size_t> cvt = ib;

177
178
				i_box_id bid_t;
				bid_t.box = cvt;
179
				bid_t.g_id = g.LinId(bid_t.box.middle().asArray());
180
				bid_t.sub = dec.getProcessorIGhostSub(i,j);
181
182
183
184
185
186
187
				pib.bid.add(bid_t);
			}
		}

		init_i_g_box = true;
	}

188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
	/*! \brief Create per-processor internal ghost box list in grid units
	 *
	 */
	void create_eg_box()
	{
		// Get the grid info
		auto g = cd_sm.getGrid();

		if (init_e_g_box == true)	return;

		// Get the number of near processors
		for (size_t i = 0 ; i < dec.getNNProcessors() ; i++)
		{
			eg_box.add();
			auto&& pib = eg_box.last();

			pib.prc = dec.IDtoProc(i);
			for (size_t j = 0 ; j < dec.getProcessorNEGhost(i) ; j++)
			{
incardon's avatar
incardon committed
207
208
209
210
211
212
213
				// Get the external ghost boxes and transform into grid units
				::Box<dim,St> ib_dom = dec.getProcessorEGhostBox(i,j);
				::Box<dim,size_t> ib = cd_sm.convertDomainSpaceIntoGridUnits(ib_dom);

				// Check if ib is valid if not it mean that the internal ghost does not contain information so skip it
				if (ib.isValid() == false)
					continue;
214
215
216

				// save the box and the unique external ghost box id (linearization of P1)
				// It is (locally) unique because it is ensured that external ghost boxes does not overlap
incardon's avatar
incardon committed
217
				// Carefull it is not unique from the internal ghost box
218
219
220
221
222
223

				// sub domain id at which belong the external ghost box
				size_t sub_id = dec.getProcessorEGhostSub(i,j);

				e_box_id bid_t;
				bid_t.sub = sub_id;
224
225
				bid_t.g_e_box = ib;
				bid_t.l_e_box = ib;
226
				// Translate in local coordinate
227
				Box<dim,long int> tb = ib;
incardon's avatar
incardon committed
228
229
				tb -= gdb_ext.get(sub_id).origin;
				bid_t.l_e_box = tb;
230
231

				pib.bid.add(bid_t);
incardon's avatar
incardon committed
232
233

				// Add the map between the global ghost box id and id of the external box in the vector
234
235
				size_t g_id = g.LinId(ib.middle().asArray());
				g_id_to_external_ghost_box[g_id] = pib.bid.size()-1;
236
237
238
239
240
241
			}
		}

		init_e_g_box = true;
	}

incardon's avatar
incardon committed
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
	bool init_local_i_g_box = false;

	/*! \brief Create local internal ghost box in grid units
	 *
	 */
	void create_local_ig_box()
	{
		// Get the grid info
		auto g = cd_sm.getGrid();

		if (init_local_i_g_box == true)	return;

		// Get the number of near processors
		for (size_t i = 0 ; i < dec.getNLocalHyperCube() ; i++)
		{
			loc_ig_box.add();
			auto&& pib = loc_ig_box.last();

			for (size_t j = 0 ; j < dec.getLocalNIGhost(i) ; j++)
			{
				// Get the internal ghost boxes and transform into grid units
				::Box<dim,St> ib_dom = dec.getLocalIGhostBox(i,j);
				::Box<dim,size_t> ib = cd_sm.convertDomainSpaceIntoGridUnits(ib_dom);

				// Check if ib is valid if not it mean that the internal ghost does not contain information so skip it
				if (ib.isValid() == false)
					continue;

incardon's avatar
incardon committed
270
271
272
273
				pib.bid.add();
				pib.bid.last().box = ib;
				pib.bid.last().sub = dec.getLocalIGhostSub(i,j);
				pib.bid.last().k = dec.getLocalIGhostE(i,j);
incardon's avatar
incardon committed
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
			}
		}

		init_local_i_g_box = true;
	}

	bool init_local_e_g_box = false;

	/*! \brief Create per-processor internal ghost box list in grid units
	 *
	 */
	void create_local_eg_box()
	{
		// Get the grid info
		auto g = cd_sm.getGrid();

		if (init_local_e_g_box == true)	return;

		// Get the number of near processors
		for (size_t i = 0 ; i < dec.getNLocalHyperCube() ; i++)
		{
			loc_eg_box.add();
			auto&& pib = loc_eg_box.last();

			for (size_t j = 0 ; j < dec.getLocalNEGhost(i) ; j++)
			{
				// Get the internal ghost boxes and transform into grid units
				::Box<dim,St> ib_dom = dec.getLocalEGhostBox(i,j);
				::Box<dim,size_t> ib = cd_sm.convertDomainSpaceIntoGridUnits(ib_dom);

incardon's avatar
incardon committed
304
305
306
				// Warning even if the ib is not a valid in grid unit we are forced to keep it
				// otherwise the value returned from dec.getLocalEGhostSub(i,j) will point to an
				// invalid or wrong box
incardon's avatar
incardon committed
307

incardon's avatar
incardon committed
308
309
310
				pib.bid.add();
				pib.bid.last().box = ib;
				pib.bid.last().sub = dec.getLocalEGhostSub(i,j);
incardon's avatar
incardon committed
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
			}
		}

		init_local_e_g_box = true;
	}

	/*! \brief Sync the local ghost part
	 *
	 * \tparam prp... properties to sync
	 *
	 */
	template<int... prp> void ghost_get_local()
	{
		//! For all the sub-domains
		for (size_t i = 0 ; i < loc_ig_box.size() ; i++)
		{
			//! For all the internal ghost boxes of each sub-domain
			for (size_t j = 0 ; j < loc_ig_box.get(i).bid.size() ; j++)
			{
incardon's avatar
incardon committed
330
331
332
				Box<dim,size_t> bx_src = loc_ig_box.get(i).bid.get(j).box;
				// convert into local
				bx_src -= gdb_ext.get(i).origin;
incardon's avatar
incardon committed
333

incardon's avatar
incardon committed
334
335
				// sub domain connected with external box
				size_t sub_id_dst = loc_ig_box.get(i).bid.get(j).sub;
incardon's avatar
incardon committed
336
337

				// local external ghost box connected
incardon's avatar
incardon committed
338
				size_t k = loc_ig_box.get(i).bid.get(j).k;
incardon's avatar
incardon committed
339

incardon's avatar
incardon committed
340
341
342
343
				Box<dim,size_t> bx_dst = loc_eg_box.get(sub_id_dst).bid.get(k).box;

				// convert into local
				bx_dst -= gdb_ext.get(sub_id_dst).origin;
incardon's avatar
incardon committed
344
345
346

				// create 2 sub grid iterator
				grid_key_dx_iterator_sub<dim> sub_src(loc_grid.get(i).getGrid(),bx_src.getKP1(),bx_src.getKP2());
incardon's avatar
incardon committed
347
				grid_key_dx_iterator_sub<dim> sub_dst(loc_grid.get(sub_id_dst).getGrid(),bx_dst.getKP1(),bx_dst.getKP2());
incardon's avatar
incardon committed
348
349
350

#ifdef DEBUG

incardon's avatar
incardon committed
351
352
353
				if (loc_eg_box.get(sub_id_dst).bid.get(k).sub != i)
					std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " source and destination are not correctly linked" << "\n";

incardon's avatar
incardon committed
354
355
356
357
358
359
				if (sub_src.getVolume() != sub_dst.getVolume())
					std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " source and destination does not match in size" << "\n";

#endif

				const auto & gs = loc_grid.get(i);
incardon's avatar
incardon committed
360
				auto & gd = loc_grid.get(sub_id_dst);
incardon's avatar
incardon committed
361
362
363
364

				while (sub_src.isNext())
				{
					// Option 1
incardon's avatar
incardon committed
365
					gd.set(sub_dst.get(),gs,sub_src.get());
incardon's avatar
incardon committed
366
367

					// Option 2
incardon's avatar
incardon committed
368
//					gd.get_o(sub_dst.get()) = gs.get_o(sub_src.get());
incardon's avatar
incardon committed
369
370
371
372
373
374
375
376

					++sub_src;
					++sub_dst;
				}
			}
		}
	}

incardon's avatar
incardon committed
377
378
379
380
381
382
383
384
385
386
387
388
389
	/*! \brief Check the the grid has valid size
	 *
	 * Distributed grids with size < 2 on each dimension are not supported
	 *
	 */
	inline void check_size(const size_t (& g_sz)[dim])
	{
		for (size_t i = 0 ; i < dim ; i++)
		{
			if (g_sz[i] < 2)
				std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " distrobuted grids with size smaller than 2 are not supported\n";
		}
	}
incardon's avatar
incardon committed
390

incardon's avatar
incardon committed
391
392
393
public:

	//! constructor
394
	grid_dist_id(Vcluster v_cl, Decomposition & dec, const size_t (& g_sz)[dim], const Box<dim,St> & domain, const Ghost<dim,T> & ghost)
incardon's avatar
incardon committed
395
	:domain(domain),ghost(ghost),loc_grid(NULL),v_cl(v_cl),dec(dec)
incardon's avatar
incardon committed
396
	{
incardon's avatar
incardon committed
397
398
		check_size(g_sz);

incardon's avatar
incardon committed
399
400
401
402
403
404
405
		// For a 5x5 grid you have 4x4 Cell
		size_t c_g[dim];
		for (size_t i = 0 ; i < dim ; i++)	{c_g[i] = g_sz[i]-1;}

		// Initialize the cell decomposer
		cd_sm.setDimensions(domain,c_g,0);

incardon's avatar
incardon committed
406
407
408
409
410
411
412
413
414
415
416
417
		// fill the global size of the grid
		for (int i = 0 ; i < dim ; i++)	{this->g_sz[i] = g_sz[i];}

		// Get the number of processor and calculate the number of sub-domain
		// for decomposition
		size_t n_proc = v_cl.getProcessingUnits();
		size_t n_sub = n_proc * SUB_UNIT_FACTOR;

		// Calculate the maximum number (before merging) of sub-domain on
		// each dimension
		size_t div[dim];
		for (int i = 0 ; i < dim ; i++)
418
		{div[i] = openfpm::math::round_big_2(pow(n_sub,1.0/dim));}
incardon's avatar
incardon committed
419
420
421

		// Create the sub-domains
		dec.setParameters(div);
incardon's avatar
incardon committed
422
423
424

		// Create local grid
		Create();
425
426
427

		// Calculate ghost boxes
		dec.calculateGhostBoxes(ghost);
incardon's avatar
incardon committed
428
429
	}

430
431
432
433
434
435
	/*! \brief Constrcuctor
	 *
	 * \param g_sz array with the grid size on each dimension
	 * \param domain
	 *
	 */
436
	grid_dist_id(const size_t (& g_sz)[dim],const Box<dim,St> & domain, const Ghost<dim,St> & g)
incardon's avatar
incardon committed
437
	:domain(domain),ghost(g),dec(Decomposition(*global_v_cluster)),v_cl(*global_v_cluster)
incardon's avatar
incardon committed
438
	{
incardon's avatar
incardon committed
439
440
441
		// check that the grid has valid size
		check_size(g_sz);

incardon's avatar
incardon committed
442
443
		// For a 5x5 grid you have 4x4 Cell
		size_t c_g[dim];
incardon's avatar
incardon committed
444
		for (size_t i = 0 ; i < dim ; i++)	{c_g[i] = (g_sz[i]-1 > 0)?(g_sz[i]-1):1;}
incardon's avatar
incardon committed
445
446
447
448

		// Initialize the cell decomposer
		cd_sm.setDimensions(domain,c_g,0);

incardon's avatar
incardon committed
449
		// fill the global size of the grid
450
		for (size_t i = 0 ; i < dim ; i++)	{this->g_sz[i] = g_sz[i];}
incardon's avatar
incardon committed
451

incardon's avatar
incardon committed
452
453
454
455
456
457
458
459
		// Get the number of processor and calculate the number of sub-domain
		// for decomposition
		size_t n_proc = v_cl.getProcessingUnits();
		size_t n_sub = n_proc * SUB_UNIT_FACTOR;

		// Calculate the maximum number (before merging) of sub-domain on
		// each dimension
		size_t div[dim];
460
		for (size_t i = 0 ; i < dim ; i++)
461
		{div[i] = openfpm::math::round_big_2(pow(n_sub,1.0/dim));}
incardon's avatar
incardon committed
462
463

		// Create the sub-domains
incardon's avatar
incardon committed
464
		dec.setParameters(div,domain,ghost);
incardon's avatar
incardon committed
465
466

		// Create local grid
incardon's avatar
incardon committed
467
		Create();
468
469

		// Calculate ghost boxes
incardon's avatar
incardon committed
470
		dec.calculateGhostBoxes();
incardon's avatar
incardon committed
471
472
473
474
475
476
477
478
479
480
481
482
483
	}

	/*! \brief Get the object that store the decomposition information
	 *
	 * \return the decomposition object
	 *
	 */

	Decomposition & getDecomposition()
	{
		return dec;
	}

incardon's avatar
incardon committed
484
485
486
487
488
489
490
491
492
493
494
	/*! \brief Return the cell decomposer
	 *
	 * \return the cell decomposer
	 *
	 */
	const CellDecomposer_sm<dim,St> & getCellDecomposer()
	{
		return cd_sm;
	}

	/*! \brief Create the grids on memory
incardon's avatar
incardon committed
495
496
497
498
499
	 *
	 */

	void Create()
	{
500
501
		Box<dim,St> g_rnd_box;
		for (size_t i = 0 ; i < dim ; i++)	{g_rnd_box.setHigh(i,0.5); g_rnd_box.setLow(i,-0.5);}
502

incardon's avatar
incardon committed
503
504
505
506
507
508
509
510
511
512
513
514
		// Get the number of local grid needed
		size_t n_grid = dec.getNLocalHyperCube();

		// create local grids for each hyper-cube
		loc_grid = v_cl.allocate<device_grid>(n_grid);

		// Size of the grid on each dimension
		size_t l_res[dim];

		// Allocate the grids
		for (size_t i = 0 ; i < n_grid ; i++)
		{
incardon's avatar
incardon committed
515
516
			gdb_ext.add();

incardon's avatar
incardon committed
517
			// Get the local hyper-cube
518
			SpaceBox<dim,St> sp = dec.getLocalHyperCube(i);
incardon's avatar
incardon committed
519
			SpaceBox<dim,St> sp_g = dec.getSubDomainWithGhost(i);
incardon's avatar
incardon committed
520

incardon's avatar
incardon committed
521
			// Convert from SpaceBox<dim,St> to SpaceBox<dim,long int>
incardon's avatar
incardon committed
522
			SpaceBox<dim,long int> sp_t = cd_sm.convertDomainSpaceIntoGridUnits(sp);
incardon's avatar
incardon committed
523
			SpaceBox<dim,long int> sp_tg = cd_sm.convertDomainSpaceIntoGridUnits(sp_g);
524

525
			//! Save the origin of the local grid
incardon's avatar
incardon committed
526
			gdb_ext.last().origin = sp_tg.getP1();
incardon's avatar
incardon committed
527

528
529
			// save information about the local grid: domain box seen inside the domain + ghost box (see GDBoxes for a visual meaning)
			// and where the GDBox start, or the origin of the local grid (+ghost) in global coordinate
incardon's avatar
incardon committed
530
			gdb_ext.last().Dbox = sp_t;
incardon's avatar
incardon committed
531
			gdb_ext.last().Dbox -= sp_tg.getP1();
incardon's avatar
incardon committed
532

incardon's avatar
incardon committed
533
534
			// center to zero
			sp_tg -= sp_tg.getP1();
incardon's avatar
incardon committed
535

536
			// Get the size of the local grid
incardon's avatar
incardon committed
537
			for (size_t i = 0 ; i < dim ; i++) {l_res[i] = (sp_tg.getHigh(i) >= 0)?(sp_tg.getHigh(i)+1):0;}
538
539

			// Set the dimensions of the local grid
incardon's avatar
incardon committed
540
541
542
543
			loc_grid.get(i).template resize<Memory>(l_res);
		}
	}

incardon's avatar
incardon committed
544
	/*! \brief Check that the global grid key is inside the grid domain
incardon's avatar
incardon committed
545
	 *
incardon's avatar
incardon committed
546
	 * \return true if is inside
incardon's avatar
incardon committed
547
	 *
incardon's avatar
incardon committed
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
	 */
	bool isInside(const grid_key_dx<dim> & gk) const
	{
		for (size_t i = 0 ; i < dim ; i++)
		{
			if (gk.get(i) < 0 || gk.get(i) >= (long int)g_sz[i])
				return false;
		}

		return true;
	}

	/*! \brief It return an iterator that span the full grid domain (each processor span its local domain)
	 *
	 * \return the iterator
	 *
	 */
	grid_dist_iterator<dim,device_grid,FREE> getDomainIterator()
	{
		grid_dist_iterator<dim,device_grid,FREE> it(loc_grid,gdb_ext);

		return it;
	}

	/*! \brief It return an iterator that span the grid domain + ghost part
incardon's avatar
incardon committed
573
574
575
	 *
	 *
	 */
incardon's avatar
incardon committed
576
	grid_dist_iterator<dim,device_grid,FIXED> getDomainGhostIterator()
incardon's avatar
incardon committed
577
	{
incardon's avatar
incardon committed
578
		grid_dist_iterator<dim,device_grid,FIXED> it(loc_grid,gdb_ext);
incardon's avatar
incardon committed
579
580
581
582
583
584
585
586

		return it;
	}

	//! Destructor
	~grid_dist_id()
	{
	}
incardon's avatar
incardon committed
587
588
589
590
591
592
593
594
595
596
597

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

	Vcluster & getVC()
	{
		return v_cl;
	}
incardon's avatar
incardon committed
598
599
600
601
602
603
604
605
606
607
608
609

	/*! \brief Get the reference of the selected element
	 *
	 *
	 * \param p property to get (is an integer)
	 * \param v1 grid_key that identify the element in the grid
	 *
	 */
	template <unsigned int p>inline auto get(grid_dist_key_dx<dim> & v1) -> typename std::add_lvalue_reference<decltype(loc_grid.get(v1.getSub()).template get<p>(v1.getKey()))>::type
	{
		return loc_grid.get(v1.getSub()).template get<p>(v1.getKey());
	}
610

611
	/*! \brief it store a box, its unique id and the sub-domain from where it come from
612
613
	 *
	 */
614
	struct i_box_id
615
616
617
618
619
	{
		//! Box
		::Box<dim,size_t> box;

		//! id
620
621
622
623
624
625
		size_t g_id;

		//! sub
		size_t sub;
	};

incardon's avatar
incardon committed
626
627
628
629
630
631
632
633
634
635
636
637
638
	/*! \brief it store an internal ghost box, the linked external ghost box and the sub-domain from where
	 *  it come from as internal ghost box
	 *
	 */
	struct i_lbox_id
	{
		//! Box
		::Box<dim,size_t> box;

		//! sub-domain id
		size_t sub;

		//! external box
incardon's avatar
incardon committed
639
		size_t k;
incardon's avatar
incardon committed
640
641
	};

642
643
644
645
646
647
648
649
650
651
652
653
654
655
	/*! \brief It store the information about the external ghost box
	 *
	 *
	 */
	struct e_box_id
	{
		//! Box defining the external ghost box in global coordinates
		::Box<dim,size_t> g_e_box;

		//! Box defining the external ghost box in local coordinates
		::Box<dim,size_t> l_e_box;

		//! sub_id in which sub-domain this box live
		size_t sub;
656
657
	};

incardon's avatar
incardon committed
658
659
660
661
662
663
664
665
666
667
668
669
670
	/*! \brief It store the information about the external ghost box
	 *
	 *
	 */
	struct e_lbox_id
	{
		//! Box defining the external ghost box in local coordinates
		::Box<dim,size_t> box;

		//! sub_id in which sub-domain this box live
		size_t sub;
	};

671
	/*! \brief Per-processor Internal ghost box
672
673
	 *
	 */
674
	struct ip_box_grid
675
	{
676
677
		// ghost in grid units
		openfpm::vector<i_box_id> bid;
678
679
680
681
682

		//! processor id
		size_t prc;
	};

incardon's avatar
incardon committed
683
684
685
686
687
688
689
690
691
	/*! \brief local Internal ghost box
	 *
	 */
	struct i_lbox_grid
	{
		// ghost in grid units
		openfpm::vector<i_lbox_id> bid;
	};

692
693
694
695
696
697
698
699
700
701
702
703
	/*! \brief Per-processor external ghost box
	 *
	 */
	struct ep_box_grid
	{
		// ghost in grid units
		openfpm::vector<e_box_id> bid;

		//! processor id
		size_t prc;
	};

incardon's avatar
incardon committed
704
705
706
707
708
709
710
711
712
	/*! \brief Per-processor external ghost box
	 *
	 */
	struct e_lbox_grid
	{
		// ghost in grid units
		openfpm::vector<e_lbox_id> bid;
	};

713
714
715
716
717
718
719
720
721
722
	//! Memory for the ghost sending buffer
	Memory g_send_prp_mem;

	//! Memory for the ghost sending buffer
	Memory g_recv_prp_mem;

	//! Flag that indicate if the external ghost box has been initialized
	bool init_e_g_box = false;

	//! Flag that indicate if the internal ghost box has been initialized
723
724
725
	bool init_i_g_box = false;

	//! Internal ghost boxes in grid units
726
727
728
729
	openfpm::vector<ip_box_grid> ig_box;

	//! External ghost boxes in grid units
	openfpm::vector<ep_box_grid> eg_box;
730

incardon's avatar
incardon committed
731
732
733
734
735
736
	//! Local internal ghost boxes in grid units
	openfpm::vector<i_lbox_grid> loc_ig_box;

	//! Local external ghost boxes in grid units
	openfpm::vector<e_lbox_grid> loc_eg_box;

737
738
739
740
741
742
743
744
745
746
747
	/*! \brief It synchronize getting the ghost part of the grid
	 *
	 * \tparam prp Properties to get (sequence of properties ids)
	 * \opt options (unused)
	 *
	 */
	template<int... prp> void ghost_get()
	{
		// Sending property object
		typedef object<typename object_creator<typename T::type,prp...>::type> prp_object;

748
		// Convert the ghost  internal boxes into grid unit boxes
749
750
		create_ig_box();

751
752
		// Convert the ghost external boxes into grid unit boxes
		create_eg_box();
753

incardon's avatar
incardon committed
754
755
756
757
758
759
		// Convert the local ghost internal boxes into grid unit boxes
		create_local_ig_box();

		// Convert the local external ghost boxes into grid unit boxes
		create_local_eg_box();

760
		// total number of sending vector
761
		std::vector<size_t> pap_prp;
762

763
764
		// Create a packing request vector
		for ( size_t i = 0 ; i < ig_box.size() ; i++ )
765
766
		{
			// for each ghost box
767
			for (size_t j = 0 ; j < ig_box.get(i).bid.size() ; j++)
768
			{
769
770
771
772
773
774
775
776
777
778
779
780
				// And linked sub-domain
				size_t sub_id = ig_box.get(i).bid.get(j).sub;
				// Internal ghost box
				Box<dim,size_t> g_ig_box = ig_box.get(i).bid.get(j).box;
				g_ig_box -= gdb_ext.get(sub_id).origin.template convertPoint<size_t>();

				// Pack a size_t for the internal ghost id
				Packer<size_t,HeapMemory>::packRequest(pap_prp);
				// Create a sub grid iterator spanning the internal ghost layer
				grid_key_dx_iterator_sub<dim> sub_it(loc_grid.get(sub_id).getGrid(),g_ig_box.getKP1(),g_ig_box.getKP2());
				// and pack the internal ghost grid
				Packer<device_grid,HeapMemory>::template packRequest<prp...>(loc_grid.get(sub_id),sub_it,pap_prp);
781
782
783
784
			}
		}

		// resize the property buffer memory
785
		g_send_prp_mem.resize(ExtPreAlloc<Memory>::calculateMem(pap_prp));
786
787

		// Create an object of preallocated memory for properties
788
789
		ExtPreAlloc<Memory> & prAlloc_prp = *(new ExtPreAlloc<Memory>(pap_prp,g_send_prp_mem));
		prAlloc_prp.incRef();
790

791
792
		// Pack information
		Pack_stat sts;
793

794
795
		// Pack the information for each processor and send it
		for ( size_t i = 0 ; i < ig_box.size() ; i++ )
796
		{
797
798
			sts.mark();

799
			// for each ghost box
800
			for (size_t j = 0 ; j < ig_box.get(i).bid.size() ; j++)
801
			{
802
803
804
805
806
807
808
809
810
811
812
813
814
815
				// And linked sub-domain
				size_t sub_id = ig_box.get(i).bid.get(j).sub;
				// Internal ghost box
				Box<dim,size_t> g_ig_box = ig_box.get(i).bid.get(j).box;
				g_ig_box -= gdb_ext.get(sub_id).origin.template convertPoint<size_t>();
				// Ghost box global id
				size_t g_id = ig_box.get(i).bid.get(j).g_id;

				// Pack a size_t for the internal ghost id
				Packer<size_t,HeapMemory>::pack(prAlloc_prp,g_id,sts);
				// Create a sub grid iterator spanning the internal ghost layer
				grid_key_dx_iterator_sub<dim> sub_it(loc_grid.get(sub_id).getGrid(),g_ig_box.getKP1(),g_ig_box.getKP2());
				// and pack the internal ghost grid
				Packer<device_grid,HeapMemory>::template pack<prp...>(prAlloc_prp,loc_grid.get(sub_id),sub_it,sts);
816
			}
817
818
			// send the request
			v_cl.send(ig_box.get(i).prc,0,sts.getMarkPointer(prAlloc_prp),sts.getMarkSize(prAlloc_prp));
819
820
		}

821
822
		// Calculate the total information to receive from each processors
		std::vector<size_t> prp_recv;
823

824
825
826
827
		//! Receive the information from each processors
		for ( size_t i = 0 ; i < eg_box.size() ; i++ )
		{
			prp_recv.push_back(0);
828

829
830
831
832
833
834
			// for each external ghost box
			for (size_t j = 0 ; j < eg_box.get(i).bid.size() ; j++)
			{
				// External ghost box
				Box<dim,size_t> g_eg_box = eg_box.get(i).bid.get(j).g_e_box;
				prp_recv[prp_recv.size()-1] += g_eg_box.getVolumeKey() * sizeof(prp_object) + sizeof(size_t);
835
836
837
			}
		}

838
839
		//! Resize the receiving buffer
		g_recv_prp_mem.resize(ExtPreAlloc<Memory>::calculateMem(prp_recv));
840

841
		// Create an object of preallocated memory for properties
842
		ExtPreAlloc<Memory> & prRecv_prp = *(new ExtPreAlloc<Memory>(prp_recv,g_recv_prp_mem));
843
		prRecv_prp.incRef();
844

incardon's avatar
incardon committed
845
		// queue the receives
846
847
		for ( size_t i = 0 ; i < eg_box.size() ; i++ )
		{
848
849
			prRecv_prp.allocate(prp_recv[i]);
			v_cl.recv(eg_box.get(i).prc,0,prRecv_prp.getPointer(),prp_recv[i]);
850
		}
851

incardon's avatar
incardon committed
852
853
854
855
856
		// Before wait for the communication to complete we sync the local ghost
		// in order to overlap with communication

		ghost_get_local<prp...>();

857
858
859
		// wait to receive communication
		v_cl.execute();

incardon's avatar
incardon committed
860
		Unpack_stat ps;
861
862
863

		// Unpack the object
		for ( size_t i = 0 ; i < eg_box.size() ; i++ )
864
		{
865
866
			// for each external ghost box
			for (size_t j = 0 ; j < eg_box.get(i).bid.size() ; j++)
867
			{
868
				// Unpack the ghost box global-id
869

870
871
872
873
874
875
876
877
878
879
880
				size_t g_id;
				Unpacker<size_t,HeapMemory>::unpack(prRecv_prp,g_id,ps);

				size_t l_id = 0;
				// convert the global id into local id
				auto key = g_id_to_external_ghost_box.find(g_id);
				if (key != g_id_to_external_ghost_box.end()) // FOUND
					l_id = key->second;
				else
				{
					// NOT FOUND
881

882
883
884
					// It must be always found, if not it mean that the processor has no-idea of
					// what is stored and conseguently do not know how to unpack, print a critical error
					// and return
885

886
887
888
889
890
891
892
893
					std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Critical, cannot unpack object, because received data cannot be interpreted\n";

					return;
				}

				// Get the external ghost box associated with the packed information
				Box<dim,size_t> box = eg_box.get(i).bid.get(l_id).l_e_box;
				size_t sub_id = eg_box.get(i).bid.get(l_id).sub;
894

895
896
				// sub-grid where to unpack
				grid_key_dx_iterator_sub<dim> sub2(loc_grid.get(sub_id).getGrid(),box.getKP1(),box.getKP2());
897

898
899
900
901
				// Unpack
				Unpacker<device_grid,HeapMemory>::template unpack<prp...>(prRecv_prp,sub2,loc_grid.get(sub_id),ps);
			}
		}
902
903
	}

incardon's avatar
incardon committed
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
	/*! \brief Convert a g_dist_key_dx into a global key
	 *
	 * \see grid_dist_key_dx
	 * \see grid_dist_key_dx_iterator
	 *
	 * \return the global position in the grid
	 *
	 */
	inline grid_key_dx<dim> getGKey(const grid_dist_key_dx<dim> & k)
	{
		// Get the sub-domain id
		size_t sub_id = k.getSub();

		grid_key_dx<dim> k_glob = k.getKey();

		// shift
		k_glob = k_glob + gdb_ext.get(sub_id).origin;

		return k_glob;
	}

925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
	/*! \brief Write the grid_dist_id information as VTK file
	 *
	 * The function generate several files
	 *
	 * 1)
	 *
	 * where X is the processor number
	 *
	 * \param output directory where to write the files
	 *
	 */
	bool write(std::string output) const
	{
		VTKWriter<openfpm::vector<::Box<dim,size_t>>,VECTOR_BOX> vtk_box1;

		openfpm::vector< openfpm::vector< ::Box<dim,size_t> > > boxes;

		//! Carefully we have to ensure that boxes does not reallocate inside the for loop
		boxes.reserve(ig_box.size());

		//! Write internal ghost in grid units (Color encoded)
		for (size_t p = 0 ; p < ig_box.size() ; p++)
		{
			boxes.add();

			// Create a vector of boxes
			for (size_t j = 0 ; j < ig_box.get(p).bid.size() ; j++)
			{
				boxes.last().add(ig_box.get(p).bid.get(j).box);
			}

			vtk_box1.add(boxes.last());
		}
		vtk_box1.write(output + std::string("internal_ghost_") + std::to_string(v_cl.getProcessUnitID()) + std::string(".vtk"));

		return true;
	}
incardon's avatar
incardon committed
962
963
964
965
966
967
968
};

/*! \brief This is a distributed grid
 *
 * Implementation of a distributed grid with id decomposition. A distributed grid is a grid distributed
 * across processors. The decomposition is performed on the id of the elements
 *
969
 * 1D specialization
incardon's avatar
incardon committed
970
971
972
973
974
975
976
977
978
979
980
981
 *
 * \param dim Dimensionality of the grid
 * \param T type of grid
 * \param Decomposition Class that decompose the grid for example CartDecomposition
 * \param Mem Is the allocator
 * \param device type of base structure is going to store the data
 *
 */

template<typename T, typename Decomposition,typename Memory , typename device_grid >
class grid_dist_id<1,T,Decomposition,Memory,device_grid>
{
982
983
	// Ghost
	Ghost<1,T> ghost;
incardon's avatar
incardon committed
984
985
986
987
988
989
990
991
992
993

	//! Local grids
	Vcluster_object_array<device_grid> loc_grid;

	//! Size of the grid on each dimension
	size_t g_sz[1];

	//! Communicator class
	Vcluster & v_cl;

incardon's avatar
incardon committed
994
995
996
	//! Extension of each grid: Domain and ghost + domain
	openfpm::vector<GBoxes<device_grid::dims>> gdb_ext;

incardon's avatar
incardon committed
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
	/*! \brief Get the grid size
	 *
	 * Get the grid size, given a domain, the resolution on it and another spaceBox
	 * it give the size on all directions of the local grid
	 *
	 * \param sp SpaceBox enclosing the local grid
	 * \param domain Space box enclosing the physical domain or part of it
	 * \param v_size grid size on this physical domain
	 *
	 * \return An std::vector representing the local grid on each dimension
	 *
	 */
	std::vector<size_t> getGridSize(SpaceBox<1,typename Decomposition::domain_type> & sp, Box<1,typename Decomposition::domain_type> & domain, size_t (& v_size)[1])
	{
		std::vector<size_t> tmp;
		for (size_t d = 0 ; d < 1 ; d++)
		{
			// push the size of the local grid
			tmp.push_back(g_sz[0]);
		}
		return tmp;
	}

	/*! \brief Get the grid size
	 *
	 * Get the grid size, given a spaceBox
	 * it give the size on all directions of the local grid
	 *
	 * \param sp SpaceBox enclosing the local grid
	 * \param sz array to fill with the local grid size on each dimension
	 *
	 */
	void getGridSize(SpaceBox<1,size_t> & sp, size_t (& v_size)[1])
	{
		for (size_t d = 0 ; d < 1 ; d++)
		{
			// push the size of the local grid
incardon's avatar
incardon committed
1034
			v_size[d] = sp.getHigh(d) - sp.getLow(d);
incardon's avatar
incardon committed
1035
1036
1037
1038
1039
1040
		}
	}

public:

	//! constructor
1041
	grid_dist_id(Vcluster v_cl, Decomposition & dec, size_t (& g_sz)[1], Box<1,T> & ghost)
incardon's avatar
incardon committed
1042
1043
	:ghost(ghost),loc_grid(NULL),v_cl(v_cl)
	{
incardon's avatar
incardon committed
1044
1045
1046
		// All this code is completely untested to assume broken
		std::cout << "Error: " << __FILE__ << ":" << __LINE__ << " this structure is untested to assume broken ";

incardon's avatar
incardon committed
1047
1048
1049
1050
1051
1052
1053
1054
1055
		// fill the global size of the grid
		for (int i = 0 ; i < 1 ; i++)	{this->g_sz[i] = g_sz[i];}

		// Create local memory
		Create();
	}

	//! constructor
	grid_dist_id(size_t (& g_sz)[1])
1056
	:v_cl(*global_v_cluster),ghost(0)
incardon's avatar
incardon committed
1057
	{
incardon's avatar
incardon committed
1058
1059
1060
		// All this code is completely untested to assume broken
		std::cout << "Error: " << __FILE__ << ":" << __LINE__ << " this structure is untested to assume broken ";

incardon's avatar
incardon committed
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
		// fill the global size of the grid
		for (int i = 0 ; i < 1 ; i++)	{this->g_sz[i] = g_sz[i];}

		// Create local memory
		Create();
	}

	/*! \brief Create the grid on memory
	 *
	 */

	void Create()
	{
		size_t n_grid = 1;

		// create local grids for each hyper-cube
		loc_grid = v_cl.allocate<device_grid>(n_grid);

		// Size of the grid on each dimension
		size_t l_res[1];

		// Calculate the local grid size
		l_res[0] = g_sz[0] / v_cl.getProcessingUnits();

		// Distribute the remaining
		size_t residual = g_sz[0] % v_cl.getProcessingUnits();
		if (v_cl.getProcessUnitID() < residual)
			l_res[0]++;

		// Set the dimensions of the local grid
		loc_grid.get(0).template resize<Memory>(l_res);
	}

	/*! \brief It return an iterator of the bulk part of the grid with a specified margin
	 *
	 * For margin we mean that every point is at least m points far from the border
	 *
	 * \param m margin
	 *
	 * \return An iterator to a grid with specified margins
	 *
	 */
incardon's avatar
incardon committed
1103
	grid_dist_iterator<1,device_grid,FREE> getDomainIterator()
incardon's avatar
incardon committed
1104
	{
incardon's avatar
incardon committed
1105
		grid_dist_iterator<1,device_grid,FREE> it(loc_grid,gdb_ext);
incardon's avatar
incardon committed
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136

		return it;
	}

	//! Destructor
	~grid_dist_id()
	{
	}

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

	Vcluster & getVC()
	{
		return v_cl;
	}

	/*! \brief Get the reference of the selected element
	 *
	 *
	 * \param p property to get (is an integer)
	 * \param v1 grid_key that identify the element in the grid
	 *
	 */
	template <unsigned int p>inline auto get(grid_dist_key_dx<1> & v1) -> typename std::add_lvalue_reference<decltype(loc_grid.get(v1.getSub()).template get<p>(v1.getKey()))>::type
	{
		return loc_grid.get(0).template get<p>(v1.getKey());
	}
incardon's avatar
incardon committed
1137
1138
1139
};

#endif