grid_dist_id.hpp source code [openfpm/src/Grid/grid_dist_id.hpp]

1	#ifndef COM_UNIT_HPP
2	#define COM_UNIT_HPP
3
4	#include <vector>
5	#include <unordered_map>
6	#include "Grid/map_grid.hpp"
7	#include "VCluster/VCluster.hpp"
8	#include "Space/SpaceBox.hpp"
9	#include "util/mathutil.hpp"
10	#include "VTKWriter/VTKWriter.hpp"
11	#ifdef __NVCC__
12	#include "SparseGridGpu/SparseGridGpu.hpp"
13	#endif
14	#include "Iterators/grid_dist_id_iterator_dec.hpp"
15	#include "Iterators/grid_dist_id_iterator.hpp"
16	#include "Iterators/grid_dist_id_iterator_sub.hpp"
17	#include "grid_dist_key.hpp"
18	#include "NN/CellList/CellDecomposer.hpp"
19	#include "util/object_util.hpp"
20	#include "memory/ExtPreAlloc.hpp"
21	#include "Packer_Unpacker/Packer.hpp"
22	#include "Packer_Unpacker/Unpacker.hpp"
23	#include "Decomposition/CartDecomposition.hpp"
24	#include "data_type/aggregate.hpp"
25	#include "hdf5.h"
26	#include "grid_dist_id_comm.hpp"
27	#include "HDF5_wr/HDF5_wr.hpp"
28	#include "SparseGrid/SparseGrid.hpp"
29	#include "lib/pdata.hpp"
30	#ifdef __NVCC__
31	#include "cuda/grid_dist_id_kernels.cuh"
32	#include "Grid/cuda/grid_dist_id_iterator_gpu.cuh"
33	#endif
34
35	/! \brief It contain the offset necessary to move to coarser and finer level grids*
36	*
37	*/
38	template<unsigned int dim>
39	struct offset_mv
40	{
41	//! offset to move up on an upper grid (coarse)
42	Point<dim,long int> up;
43
44	//! offset to move on the lower grid (finer)
45	Point<dim,long int> dw;
46	};
47
48	//! Internal ghost box sent to construct external ghost box into the other processors
49	template<unsigned int dim>
50	struct Box_fix
51	{
52	//! Box in global unit
53	Box<dim,size_t> bx;
54	//! In which sector live the box
55	comb<dim> cmb;
56	//! Global id of the internal ghost box
57	size_t g_id;
58	//! from which sub-domain this internal ghost box is generated (or with which sub-domain is overlapping)
59	size_t r_sub;
60	};
61
62	#define NO_GDB_EXT_SWITCH 0x1000
63
64	/! \brief This is a distributed grid*
65	*
66	* Implementation of a distributed grid the decomposition is geometrical, grid
67	* is splitted across several processor
68	*
69	* \param dim Dimensionality of the grid
70	* \param St Type of space where the grid is living
71	* \param T object the grid is storing
72	* \param Decomposition Class that decompose the grid for example CartDecomposition
73	* \param Mem Is the allocator
74	* \param device type of base structure is going to store the data
75	*
76	* ### Create a distributed grid and access it
77	* \snippet grid_dist_id_unit_test.cpp Create and access a distributed grid
78	* ### Synchronize the ghosts and check the information
79	* \snippet grid_dist_id_unit_test.cpp Synchronize the ghost and check the information
80	* ### Create and access a distributed grid for complex structures
81	* \snippet grid_dist_id_unit_test.cpp Create and access a distributed grid complex
82	* ### Synchronize a distributed grid for complex structures
83	* \snippet grid_dist_id_unit_test.cpp Synchronized distributed grid complex
84	* ### Usage of a grid dist iterator sub
85	* \snippet grid_dist_id_iterators_unit_tests.hpp Usage of a sub_grid iterator
86	* ### Construct two grid with the same decomposition
87	* \snippet grid_dist_id_unit_test.cpp Construct two grid with the same decomposition
88	*
89	*/
90	template<unsigned int dim,
91	typename St,
92	typename T,
93	typename Decomposition = CartDecomposition<dim,St>,
94	typename Memory=HeapMemory ,
95	typename device_grid=grid_cpu<dim,T> >
96	class grid_dist_id : public grid_dist_id_comm<dim,St,T,Decomposition,Memory,device_grid>
97	{
98	typedef grid_dist_id<dim,St,T,Decomposition,Memory,device_grid> self;
99
100	//! Domain
101	Box<dim,St> domain;
102
103	//! Ghost expansion
104	Ghost<dim,St> ghost;
105
106	//! Ghost expansion
107	Ghost<dim,long int> ghost_int;
108
109	//! Local grids
110	mutable openfpm::vector<device_grid> loc_grid;
111
112	//! Old local grids
113	mutable openfpm::vector<device_grid> loc_grid_old;
114
115	//! Space Decomposition
116	Decomposition dec;
117
118	//! gdb_ext markers
119	//! In the case where the grid is defined everywhere
120	//! gdb_ext_marker is useless and so is empty
121	//! in the case we have a grid defined on a smaller set
122	//! of boxes gbd_ext_markers indicate the division across
123	//! subdomains. For example Sub-domain 0 produce 2 grid
124	//! Sub-domain 1 produce 3 grid Sub-domain 2 produce 2 grid
125	//! Sub-domain 3 produce 1 grid
126	//! gdb_ext_markers contain 0,2,5,7,8
127	openfpm::vector<size_t> gdb_ext_markers;
128
129	//! Extension of each grid: Domain and ghost + domain
130	openfpm::vector<GBoxes<device_grid::dims>> gdb_ext;
131
132	//! Global gdb_ext
133	mutable openfpm::vector<GBoxes<device_grid::dims>> gdb_ext_global;
134
135	//! Extension of each old grid (old): Domain and ghost + domain
136	openfpm::vector<GBoxes<device_grid::dims>> gdb_ext_old;
137
138	//! Size of the grid on each dimension
139	size_t g_sz[dim];
140
141	//! Structure that divide the space into cells
142	CellDecomposer_sm<dim,St,shift<dim,St>> cd_sm;
143
144	//! number of insert each GPU thread does
145	size_t gpu_n_insert_thread;
146
147	//! Communicator class
148	Vcluster<> & v_cl;
149
150	//! properties names
151	openfpm::vector<std::string> prp_names;
152
153	//! It map a global ghost id (g_id) to the external ghost box information
154	//! It is unique across all the near processor
155	std::unordered_map<size_t,size_t> g_id_to_external_ghost_box;
156
157	/! Link a received external ghost box to the linked eg_box.*
158	* When the grid is defined everywhere for each received external ghost box
159	* exists one eg_box linked to it that contain the information
160	* on how to transfer the information to the associated sub-domain grid.
161	* Unfortunately when we specify where the grid is defined, a received external
162	* ghost box can be linked to multiple sub-domain grids (one sub-domain can have multiple
163	* sub grids).
164	* So in standard situation (grid defined everywhere)
165	* a received external ghost box is linked to a single
166	* eg_box entry and eb_gid_list play mainly no role.
167	* (play no role but must be filled ghost_get expect this structure to be
168	* filled consistently, it will be clear later how to do it in this case).
169	* When the grid is not defined everywhere a received ghost box can be linked
170	* to multiple external ghost boxes. (Like in figure)
171	*
172	* \verbatim
173	+--------------------------------------------+------
174	\| Sub-domain \| Another sub-domain
175	\| +------+ +---------+ \|
176	\| \| \| \| \| \|
177	\| \| \| \| \| \|
178	\| \| \| \| \| \|
179	\| \| 3 \| \| 4 \| \|
180	\| empty \| \| empty \| \| \|
181	\| \| \| \| \| \|
182	\| \| \| \| \| \|
183	\| \| \| \| \| \| 1
184	\| \| \| \| \| \|
185	+-+-----+----+------+-----------+---------+----+-----+----- Processor bound
186	\|##########*****#############\|*\|
187	\| \| 0
188	\| \|
189	\| \|
190	\| 9 \|
191	\| \|
192	\| \|
193	\| \|
194	+--------------------------------------------+
195
196	* \endverbatim
197	*
198	* As we can see here the grid number 9 on processo 0 has an internal ghost box
199	* The internal ghost-box is sent to processor 1 and is a received external
200	* ghost box. This external ghost box is partially shared in two separated grids.
201	* It is important to note that 3 and 4 are grid defined externally and are
202	* not defined by the sub-domain border. It is important also to note that the sub-domain
203	* granularity in processor 1 define the granularity of the internal ghost box in
204	* processor 0 and consequently every external ghost box in processor 1 is linked
205	* uniquely with one internal ghost box in processor 0. On the other hand if we have
206	* a secondary granularity define by external boxes like 3 and 4 this is not anymore
207	* true and one internal ghost box in 0 can be linked with multiple grids.
208	* The granularity of the space division is different from the granularity of where
209	* the grid is defined. Space decomposition exist independently from the data-structure
210	* and can be shared across multiple data-structure this mean that cannot be redefined
211	* based on where is the grid definitions.
212	* The internal ghost box could be redefined in order to respect the granularity.
213	* We do not do this for 3 main reason.
214	*
215	* 1) The definition box must be communicated across processors.
216	* 2) An interprocessor global-id link must be established with lower sub-domain
217	* granularty
218	* 3) Despite the points * are not linked, but must be anyway sent
219	* to processor 1, this mean that make not too much sense to increase
220	* the granularity in advance on processor 0, but it is better receive
221	* the information an than solve the lower granularity locally
222	* on processor 1
223	*
224	*/
225	openfpm::vector<e_box_multi<dim>> eb_gid_list;
226
227	//! It map a global ghost id (g_id) to the internal ghost box information
228	//! (is unique for processor), it is not unique across all the near processor
229	openfpm::vector<std::unordered_map<size_t,size_t>> g_id_to_internal_ghost_box;
230
231	//! Receiving size
232	openfpm::vector<size_t> recv_sz;
233
234	//! Receiving buffer for particles ghost get
235	openfpm::vector<HeapMemory> recv_mem_gg;
236
237	//! Grid informations object
238	grid_sm<dim,T> ginfo;
239
240	//! Grid informations object without type
241	grid_sm<dim,void> ginfo_v;
242
243	//! Set of boxes that define where the grid is defined
244	openfpm::vector<Box<dim,long int>> bx_def;
245
246	//! Indicate if we have to use bx_def to define the grid
247	bool use_bx_def = false;
248
249	//! Indicate if the local internal ghost box has been initialized
250	bool init_local_i_g_box = false;
251
252	//! Indicate if the local external ghost box has been initialized
253	bool init_local_e_g_box = false;
254
255	//! Flag that indicate if the external ghost box has been initialized
256	bool init_e_g_box = false;
257
258	//! Flag that indicate if the internal ghost box has been initialized
259	bool init_i_g_box = false;
260
261	//! Flag that indicate if the internal and external ghost box has been fixed
262	bool init_fix_ie_g_box = false;
263
264	//! Internal ghost boxes in grid units
265	openfpm::vector<ip_box_grid<dim>> ig_box;
266
267	//! External ghost boxes in grid units
268	openfpm::vector<ep_box_grid<dim>> eg_box;
269
270	//! Local internal ghost boxes in grid units
271	openfpm::vector<i_lbox_grid<dim>> loc_ig_box;
272
273	//! Local external ghost boxes in grid units
274	openfpm::vector<e_lbox_grid<dim>> loc_eg_box;
275
276	//! Number of sub-sub-domain for each processor
277	size_t v_sub_unit_factor = `64`;
278
279	/! \brief Call-back to allocate buffer to receive incoming objects (external ghost boxes)*
280	*
281	* \param msg_i message size required to receive from i
282	* \param total_msg message size to receive from all the processors
283	* \param total_p the total number of processor want to communicate with you
284	* \param i processor id
285	* \param ri request id (it is an id that goes from 0 to total_p, and is unique
286	* every time message_alloc is called)
287	* \param ptr void pointer parameter for additional data to pass to the call-back
288	*
289	*/
290	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)
291	{
292	grid_dist_id<dim,St,T,Decomposition,Memory,device_grid> * g = static_cast<grid_dist_id<dim,St,T,Decomposition,Memory,device_grid> *>(ptr);
293
294	g->recv_sz.resize(g->dec.getNNProcessors());
295	g->recv_mem_gg.resize(g->dec.getNNProcessors());
296
297	// Get the local processor id
298	size_t lc_id = g->dec.ProctoID(i);
299
300	// resize the receive buffer
301	g->recv_mem_gg.get(lc_id).resize(msg_i);
302	g->recv_sz.get(lc_id) = msg_i;
303
304	return g->recv_mem_gg.get(lc_id).getPointer();
305	}
306
307	/! \brief this function is for optimization of the ghost size*
308	*
309	* Because the decomposition work in continuum and discrete ghost is
310	* converted in continuum, in some case continuum ghost because of
311	* rounding-off error can produce ghost bigger than the discrete selected
312	* one. This function adjust for this round-off error
313	*
314	* \param sub_id sub-domain id
315	* \param sub_domain_other the other sub-domain
316	* \param ib internal ghost box to adjust
317	*
318	*/
319	void set_for_adjustment(size_t sub_id,
320	const Box<dim,St> & sub_domain_other,
321	const comb<dim> & cmb,
322	Box<dim,long int> & ib,
323	Ghost<dim,long int> & g)
324	{
325	if (g.isInvalidGhost() == true \|\| use_bx_def == true)
326	{return;}
327
328	Box<dim,long int> sub_domain;
329
330	if (use_bx_def == false)
331	{
332	sub_domain = gdb_ext.get(sub_id).Dbox;
333	sub_domain += gdb_ext.get(sub_id).origin;
334	}
335
336	// Convert from SpaceBox<dim,St> to SpaceBox<dim,long int>
337	Box<dim,long int> sub_domain_other_exp = cd_sm.convertDomainSpaceIntoGridUnits(sub_domain_other,dec.periodicity());
338
339	// translate sub_domain_other based on cmb
340	for (size_t i = `0` ; i < dim ; i++)
341	{
342	if (cmb.c[i] == `1`)
343	{
344	sub_domain_other_exp.setLow(i,sub_domain_other_exp.getLow(i) - ginfo.size(i));
345	sub_domain_other_exp.setHigh(i,sub_domain_other_exp.getHigh(i) - ginfo.size(i));
346	}
347	else if (cmb.c[i] == -`1`)
348	{
349	sub_domain_other_exp.setLow(i,sub_domain_other_exp.getLow(i) + ginfo.size(i));
350	sub_domain_other_exp.setHigh(i,sub_domain_other_exp.getHigh(i) + ginfo.size(i));
351	}
352	}
353
354	sub_domain_other_exp.enlarge(g);
355	if (sub_domain_other_exp.Intersect(sub_domain,ib) == false)
356	{
357	for (size_t i = `0` ; i < dim ; i++)
358	{ib.setHigh(i,ib.getLow(i) - `1`);}
359	}
360	}
361
362
363
364	/! \brief Create per-processor internal ghost boxes list in grid units and g_id_to_external_ghost_box*
365	*
366	*/
367	void create_ig_box()
368	{
369	// temporal vector used for computation
370	openfpm::vector_std<result_box<dim>> ibv;
371
372	if (init_i_g_box == true) return;
373
374	// Get the grid info
375	auto g = cd_sm.getGrid();
376
377	g_id_to_internal_ghost_box.resize(dec.getNNProcessors());
378
379	// Get the number of near processors
380	for (size_t i = `0` ; i < dec.getNNProcessors() ; i++)
381	{
382	ig_box.add();
383	auto&& pib = ig_box.last();
384
385	pib.prc = dec.IDtoProc(i);
386	for (size_t j = `0` ; j < dec.getProcessorNIGhost(i) ; j++)
387	{
388	// Get the internal ghost boxes and transform into grid units
389	::Box<dim,St> ib_dom = dec.getProcessorIGhostBox(i,j);
390	::Box<dim,long int> ib = cd_sm.convertDomainSpaceIntoGridUnits(ib_dom,dec.periodicity());
391
392	size_t sub_id = dec.getProcessorIGhostSub(i,j);
393	size_t r_sub = dec.getProcessorIGhostSSub(i,j);
394
395	auto & n_box = dec.getNearSubdomains(dec.IDtoProc(i));
396
397	set_for_adjustment(sub_id,
398	n_box.get(r_sub),dec.getProcessorIGhostPos(i,j),
399	ib,ghost_int);
400
401	// Here we intersect the internal ghost box with the definition boxes
402	// this operation make sense when the grid is not defined in the full
403	// domain and we have to intersect the internal ghost box with all the box
404	// that define where the grid is defined
405	bx_intersect<dim>(bx_def,use_bx_def,ib,ibv);
406
407	for (size_t k = `0` ; k < ibv.size() ; k++)
408	{
409	// Check if ib is valid if not it mean that the internal ghost does not contain information so skip it
410	if (ibv.get(k).bx.isValid() == false)
411	{continue;}
412
413	// save the box and the sub-domain id (it is calculated as the linearization of P1)
414	::Box<dim,size_t> cvt = ibv.get(k).bx;
415
416	i_box_id<dim> bid_t;
417	bid_t.box = cvt;
418	bid_t.g_id = dec.getProcessorIGhostId(i,j) \| (k) << `52`;
419
420	bid_t.sub = convert_to_gdb_ext(dec.getProcessorIGhostSub(i,j),
421	ibv.get(k).id,
422	gdb_ext,
423	gdb_ext_markers);
424
425	bid_t.cmb = dec.getProcessorIGhostPos(i,j);
426	bid_t.r_sub = dec.getProcessorIGhostSSub(i,j);
427	pib.bid.add(bid_t);
428
429	g_id_to_internal_ghost_box.get(i)[bid_t.g_id \| (k) << `52` ] = pib.bid.size()-`1`;
430	}
431	}
432	}
433
434	init_i_g_box = true;
435	}
436
437
438	/! \brief Create per-processor internal ghost box list in grid units*
439	*
440	*/
441	void create_eg_box()
442	{
443	// Get the grid info
444	auto g = cd_sm.getGrid();
445
446	if (init_e_g_box == true) return;
447
448	// Here we collect all the calculated internal ghost box in the sector different from 0 that this processor has
449
450	openfpm::vector<size_t> prc;
451	openfpm::vector<size_t> prc_recv;
452	openfpm::vector<size_t> sz_recv;
453	openfpm::vector<openfpm::vector<Box_fix<dim>>> box_int_send(dec.getNNProcessors());
454	openfpm::vector<openfpm::vector<Box_fix<dim>>> box_int_recv;
455
456	for(size_t i = `0` ; i < dec.getNNProcessors() ; i++)
457	{
458	for (size_t j = `0` ; j < ig_box.get(i).bid.size() ; j++)
459	{
460	box_int_send.get(i).add();
461	box_int_send.get(i).last().bx = ig_box.get(i).bid.get(j).box;
462	box_int_send.get(i).last().g_id = ig_box.get(i).bid.get(j).g_id;
463	box_int_send.get(i).last().r_sub = ig_box.get(i).bid.get(j).r_sub;
464	box_int_send.get(i).last().cmb = ig_box.get(i).bid.get(j).cmb;
465	}
466	prc.add(dec.IDtoProc(i));
467	}
468
469	v_cl.SSendRecv(box_int_send,box_int_recv,prc,prc_recv,sz_recv);
470
471	eg_box.resize(dec.getNNProcessors());
472
473	for (size_t i = `0` ; i < eg_box.size() ; i++)
474	{eg_box.get(i).prc = dec.IDtoProc(i);}
475
476	for (size_t i = `0` ; i < box_int_recv.size() ; i++)
477	{
478	size_t pib_cnt = `0`;
479	size_t p_id = dec.ProctoID(prc_recv.get(i));
480	auto&& pib = eg_box.get(p_id);
481	pib.prc = prc_recv.get(i);
482
483	eg_box.get(p_id).recv_pnt = `0`;
484	eg_box.get(p_id).n_r_box = box_int_recv.get(i).size();
485
486	// For each received internal ghost box
487	for (size_t j = `0` ; j < box_int_recv.get(i).size() ; j++)
488	{
489	size_t vol_recv = box_int_recv.get(i).get(j).bx.getVolumeKey();
490
491	eg_box.get(p_id).recv_pnt += vol_recv;
492	size_t send_list_id = box_int_recv.get(i).get(j).r_sub;
493
494	if (use_bx_def == true)
495	{
496	// First we understand if the internal ghost box sent intersect
497	// some local extended sub-domain.
498
499	// eb_gid_list, for an explanation check the declaration
500	eb_gid_list.add();
501	eb_gid_list.last().e_id = p_id;
502
503	// Now we have to check if a received external ghost box intersect one
504	// or more sub-grids
505	for (size_t k = `0` ; k < gdb_ext.size() ; k++)
506	{
507	Box<dim,long int> bx = gdb_ext.get(k).GDbox;
508	bx += gdb_ext.get(k).origin;
509
510	Box<dim,long int> output;
511	Box<dim,long int> flp_i = flip_box(box_int_recv.get(i).get(j).bx,box_int_recv.get(i).get(j).cmb,ginfo);
512
513	// it intersect one sub-grid
514	if (bx.Intersect(flp_i,output))
515	{
516	// link
517
518	size_t g_id = box_int_recv.get(i).get(j).g_id;
519	add_eg_box<dim>(k,box_int_recv.get(i).get(j).cmb,output,
520	g_id,
521	gdb_ext.get(k).origin,
522	box_int_recv.get(i).get(j).bx.getP1(),
523	pib.bid);
524
525	eb_gid_list.last().eb_list.add(pib.bid.size() - `1`);
526
527	g_id_to_external_ghost_box [g_id] = eb_gid_list.size() - `1`;
528	}
529	}
530
531	// now we check if exist a full match across the full intersected
532	// ghost parts
533
534	bool no_match = true;
535	for (size_t k = `0` ; k < eb_gid_list.last().eb_list.size() ; k++)
536	{
537	size_t eb_id = eb_gid_list.last().eb_list.get(k);
538
539	if (pib.bid.get(eb_id).g_e_box == box_int_recv.get(i).get(j).bx)
540	{
541	// full match found
542
543	eb_gid_list.last().full_match = eb_id;
544	no_match = false;
545
546	break;
547	}
548	}
549
550	// This is the case where a full match has not been found. In this case we
551	// generate an additional gdb_ext and local grid with only external ghost
552
553	if (no_match == true)
554	{
555	// Create a grid with the same size of the external ghost
556
557	size_t sz[dim];
558	for (size_t s = `0` ; s < dim ; s++)
559	{sz[s] = box_int_recv.get(i).get(j).bx.getHigh(s) - box_int_recv.get(i).get(j).bx.getLow(s) + `1`;}
560
561	// Add an unlinked gdb_ext
562	// An unlinked gdb_ext is an empty domain with only a external ghost
563	// part
564	Box<dim,long int> output = flip_box(box_int_recv.get(i).get(j).bx,box_int_recv.get(i).get(j).cmb,ginfo);
565
566	GBoxes<dim> tmp;
567	tmp.GDbox = box_int_recv.get(i).get(j).bx;
568	tmp.GDbox -= tmp.GDbox.getP1();
569	tmp.origin = output.getP1();
570	for (size_t s = `0` ; s < dim ; s++)
571	{
572	// we set an invalid box, there is no-domain
573	tmp.Dbox.setLow(s,`0`);
574	tmp.Dbox.setHigh(s,-`1`);
575	}
576	tmp.k = (size_t)-`1`;
577	gdb_ext.add(tmp);
578
579	// create the local grid
580
581	loc_grid.add();
582	loc_grid.last().resize(sz);
583
584	// Add an external ghost box
585
586	size_t g_id = box_int_recv.get(i).get(j).g_id;
587	add_eg_box<dim>(gdb_ext.size()-`1`,box_int_recv.get(i).get(j).cmb,output,
588	g_id,
589	gdb_ext.get(gdb_ext.size()-`1`).origin,
590	box_int_recv.get(i).get(j).bx.getP1(),
591	pib.bid);
592
593	// now we map the received ghost box to the information of the
594	// external ghost box created
595	eb_gid_list.last().full_match = pib.bid.size() - `1`;
596	eb_gid_list.last().eb_list.add(pib.bid.size() - `1`);
597	g_id_to_external_ghost_box [g_id] = eb_gid_list.size() - `1`;
598	}
599
600	// now we create lr_e_box
601
602	size_t fm = eb_gid_list.last().full_match;
603	size_t sub_id = pib.bid.get(fm).sub;
604
605	for ( ; pib_cnt < pib.bid.size() ; pib_cnt++)
606	{
607	pib.bid.get(pib_cnt).lr_e_box -= gdb_ext.get(sub_id).origin;
608	}
609
610	}
611	else
612	{
613	// Get the list of the sent sub-domains
614	// and recover the id of the sub-domain from
615	// the sent list
616	const openfpm::vector<size_t> & s_sub = dec.getSentSubdomains(p_id);
617
618	size_t sub_id = s_sub.get(send_list_id);
619
620	e_box_id<dim> bid_t;
621	bid_t.sub = sub_id;
622	bid_t.cmb = box_int_recv.get(i).get(j).cmb;
623	bid_t.cmb.sign_flip();
624	::Box<dim,long int> ib = flip_box(box_int_recv.get(i).get(j).bx,box_int_recv.get(i).get(j).cmb,ginfo);
625	bid_t.g_e_box = ib;
626	bid_t.g_id = box_int_recv.get(i).get(j).g_id;
627	// Translate in local coordinate
628	Box<dim,long int> tb = ib;
629	tb -= gdb_ext.get(sub_id).origin;
630	bid_t.l_e_box = tb;
631
632	pib.bid.add(bid_t);
633	eb_gid_list.add();
634	eb_gid_list.last().eb_list.add(pib.bid.size()-`1`);
635	eb_gid_list.last().e_id = p_id;
636	eb_gid_list.last().full_match = pib.bid.size()-`1`;
637
638	g_id_to_external_ghost_box[bid_t.g_id] = eb_gid_list.size()-`1`;
639	}
640	}
641	}
642
643	init_e_g_box = true;
644	}
645
646	/! \brief Create local internal ghost box in grid units*
647	*
648	*/
649	void create_local_ig_box()
650	{
651	openfpm::vector_std<result_box<dim>> ibv;
652
653	// Get the grid info
654	auto g = cd_sm.getGrid();
655
656	if (init_local_i_g_box == true) return;
657
658	// Get the number of sub-domains
659	for (size_t i = `0` ; i < dec.getNSubDomain() ; i++)
660	{
661	loc_ig_box.add();
662	auto&& pib = loc_ig_box.last();
663
664	for (size_t j = `0` ; j < dec.getLocalNIGhost(i) ; j++)
665	{
666	if (use_bx_def == true)
667	{
668	// Get the internal ghost boxes and transform into grid units
669	::Box<dim,St> ib_dom = dec.getLocalIGhostBox(i,j);
670	::Box<dim,long int> ib = cd_sm.convertDomainSpaceIntoGridUnits(ib_dom,dec.periodicity());
671
672	// Here we intersect the internal ghost box with the definition boxes
673	// this operation make sense when the grid is not defined in the full
674	// domain and we have to intersect the internal ghost box with all the box
675	// that define where the grid is defined
676	bx_intersect<dim>(bx_def,use_bx_def,ib,ibv);
677
678	for (size_t k = `0` ; k < ibv.size() ; k++)
679	{
680	// Check if ib is valid if not it mean that the internal ghost does not contain information so skip it
681	if (ibv.get(k).bx.isValid() == false)
682	{continue;}
683
684	pib.bid.add();
685	pib.bid.last().box = ibv.get(k).bx;
686
687	pib.bid.last().sub_gdb_ext = convert_to_gdb_ext(i,
688	ibv.get(k).id,
689	gdb_ext,
690	gdb_ext_markers);
691
692	pib.bid.last().sub = dec.getLocalIGhostSub(i,j);
693
694	// It will be filled later
695	pib.bid.last().cmb = dec.getLocalIGhostPos(i,j);
696	}
697	}
698	else
699	{
700	// Get the internal ghost boxes and transform into grid units
701	::Box<dim,St> ib_dom = dec.getLocalIGhostBox(i,j);
702	::Box<dim,long int> ib = cd_sm.convertDomainSpaceIntoGridUnits(ib_dom,dec.periodicity());
703
704	// Check if ib is valid if not it mean that the internal ghost does not contain information so skip it
705	if (ib.isValid() == false)
706	continue;
707
708	size_t sub_id = i;
709	size_t r_sub = dec.getLocalIGhostSub(i,j);
710
711	set_for_adjustment(sub_id,dec.getSubDomain(r_sub),
712	dec.getLocalIGhostPos(i,j),ib,ghost_int);
713
714	// Check if ib is valid if not it mean that the internal ghost does not contain information so skip it
715	if (ib.isValid() == false)
716	continue;
717
718	pib.bid.add();
719	pib.bid.last().box = ib;
720	pib.bid.last().sub = dec.getLocalIGhostSub(i,j);
721	pib.bid.last().k.add(dec.getLocalIGhostE(i,j));
722	pib.bid.last().cmb = dec.getLocalIGhostPos(i,j);
723	pib.bid.last().sub_gdb_ext = i;
724	}
725	}
726
727	if (use_bx_def == true)
728	{
729	// unfortunately boxes that define where the grid is located can generate
730	// additional internal ghost boxes
731
732	for (size_t j = gdb_ext_markers.get(i) ; j < gdb_ext_markers.get(i+`1`) ; j++)
733	{
734	// intersect within box in the save sub-domain
735
736	for (size_t k = gdb_ext_markers.get(i) ; k < gdb_ext_markers.get(i+`1`) ; k++)
737	{
738	if (j == k) {continue;}
739
740	// extend k and calculate the internal ghost box
741	Box<dim,long int> bx_e = gdb_ext.get(k).GDbox;
742	bx_e += gdb_ext.get(k).origin;
743	Box<dim,long int> bx = gdb_ext.get(j).Dbox;
744	bx += gdb_ext.get(j).origin;
745
746	Box<dim,long int> output;
747	if (bx.Intersect(bx_e, output) == true)
748	{
749	pib.bid.add();
750
751	pib.bid.last().box = output;
752
753	pib.bid.last().sub_gdb_ext = j;
754	pib.bid.last().sub = i;
755
756	// these ghost always in the quadrant zero
757	pib.bid.last().cmb.zero();
758
759	}
760	}
761	}
762	}
763
764	}
765
766
767	init_local_i_g_box = true;
768	}
769
770	/! \brief Create per-processor external ghost boxes list in grid units*
771	*
772	*/
773	void create_local_eg_box()
774	{
775	// Get the grid info
776	auto g = cd_sm.getGrid();
777
778	if (init_local_e_g_box == true) return;
779
780	loc_eg_box.resize(dec.getNSubDomain());
781
782	// Get the number of sub-domain
783	for (size_t i = `0` ; i < dec.getNSubDomain() ; i++)
784	{
785	for (size_t j = `0` ; j < loc_ig_box.get(i).bid.size() ; j++)
786	{
787	long int volume_linked = `0`;
788
789	size_t le_sub = loc_ig_box.get(i).bid.get(j).sub;
790	auto & pib = loc_eg_box.get(le_sub);
791
792	if (use_bx_def == true)
793	{
794
795	// We check if an external local ghost box intersect one
796	// or more sub-grids
797	for (size_t k = `0` ; k < gdb_ext.size() ; k++)
798	{
799	Box<dim,long int> bx = gdb_ext.get(k).Dbox;
800	bx += gdb_ext.get(k).origin;
801
802	Box<dim,long int> gbx = gdb_ext.get(k).GDbox;
803	gbx += gdb_ext.get(k).origin;
804
805	Box<dim,long int> output;
806	Box<dim,long int> flp_i = flip_box(loc_ig_box.get(i).bid.get(j).box,loc_ig_box.get(i).bid.get(j).cmb,ginfo);
807
808	bool intersect_domain = bx.Intersect(flp_i,output);
809	bool intersect_gdomain = gbx.Intersect(flp_i,output);
810
811	// it intersect one sub-grid
812	if (intersect_domain == false && intersect_gdomain == true)
813	{
814	// fill the link variable
815	loc_ig_box.get(i).bid.get(j).k.add(pib.bid.size());
816	size_t s = loc_ig_box.get(i).bid.get(j).k.last();
817
818	comb<dim> cmb = loc_ig_box.get(i).bid.get(j).cmb;
819	cmb.sign_flip();
820
821	add_loc_eg_box(le_sub,
822	loc_ig_box.get(i).bid.get(j).sub_gdb_ext,
823	j,
824	k,
825	pib.bid,
826	output,
827	cmb);
828
829
830	volume_linked += pib.bid.last().ebox.getVolumeKey();
831	}
832	}
833	}
834	else
835	{
836	size_t k = loc_ig_box.get(i).bid.get(j).sub;
837	auto & pib = loc_eg_box.get(k);
838
839	size_t s = loc_ig_box.get(i).bid.get(j).k.get(`0`);
840	pib.bid.resize(dec.getLocalNEGhost(k));
841
842	pib.bid.get(s).ebox = flip_box(loc_ig_box.get(i).bid.get(j).box,loc_ig_box.get(i).bid.get(j).cmb,ginfo);
843	pib.bid.get(s).sub = dec.getLocalEGhostSub(k,s);
844	pib.bid.get(s).cmb = loc_ig_box.get(i).bid.get(j).cmb;
845	pib.bid.get(s).cmb.sign_flip();
846	pib.bid.get(s).k = j;
847	pib.bid.get(s).initialized = true;
848	pib.bid.get(s).sub_gdb_ext = k;
849	}
850	}
851	}
852
853	init_local_e_g_box = true;
854	}
855
856
857	/! \brief Check the grid has a valid size*
858	*
859	* \param g_sz size of the grid
860	*
861	*/
862	inline void check_size(const size_t (& g_sz)[dim])
863	{
864	for (size_t i = `0` ; i < dim ; i++)
865	{
866	if (g_sz[i] < `2`)
867	{std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " distributed grids with size smaller than 2 are not supported\n";}
868	}
869	}
870
871	/! \brief Check the domain is valid*
872	*
873	* \param dom domain is valid
874	*
875	*/
876	inline void check_domain(const Box<dim,St> & dom)
877	{
878	for (size_t i = `0` ; i < dim ; i++)
879	{
880	if (dom.getLow(i) >= dom.getHigh(i))
881	{std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " error the simulation domain is invalid\n";}
882	}
883	}
884
885	/! \brief Create the grids on memory*
886	*
887	* \param bx_def Where the grid is defined
888	* \param use_bx_def use the array that define where the grid is defined
889	*
890	*/
891	void Create(openfpm::vector<Box<dim,long int>> & bx_def,
892	const Ghost<dim,long int> & g,
893	bool use_bx_def)
894	{
895	// create gdb
896	create_gdb_ext<dim,Decomposition>(gdb_ext,gdb_ext_markers,dec,cd_sm,bx_def,g,use_bx_def);
897
898	size_t n_grid = gdb_ext.size();
899
900	// create local grids for each hyper-cube
901	loc_grid.resize(n_grid);
902
903	// Size of the grid on each dimension
904	size_t l_res[dim];
905
906	// Allocate the grids
907	for (size_t i = `0` ; i < n_grid ; i++)
908	{
909	SpaceBox<dim,long int> sp_tg = gdb_ext.get(i).GDbox;
910
911	// Get the size of the local grid
912	// The boxes indicate the extension of the index the size
913	// is this extension +1
914	// for example a 1D box (interval) from 0 to 3 in one dimension have
915	// the points 0,1,2,3 = so a total of 4 points
916	for (size_t j = `0` ; j < dim ; j++)
917	{l_res[j] = (sp_tg.getHigh(j) >= `0`)?(sp_tg.getHigh(j)+`1`):`0`;}
918
919	// Set the dimensions of the local grid
920	loc_grid.get(i).resize(l_res);
921	}
922	}
923
924
925	/! \brief Initialize the Cell decomposer of the grid enforcing perfect overlap of the cells*
926	*
927	* \param cd_old the CellDecomposer we are trying to mach
928	* \param ext extension of the domain
929	*
930	*/
931	inline void InitializeCellDecomposer(const CellDecomposer_sm<dim,St,shift<dim,St>> & cd_old, const Box<dim,size_t> & ext)
932	{
933	// Initialize the cell decomposer
934	cd_sm.setDimensions(cd_old,ext);
935	}
936
937	/! \brief Initialize the Cell decomposer of the grid*
938	*
939	* \param g_sz Size of the grid
940	* \param bc boundary conditions
941	*
942	*/
943	inline void InitializeCellDecomposer(const size_t (& g_sz)[dim], const size_t (& bc)[dim])
944	{
945	// check that the grid has valid size
946	check_size(g_sz);
947
948	// get the size of the cell decomposer
949	size_t c_g[dim];
950	getCellDecomposerPar<dim>(c_g,g_sz,bc);
951
952	// Initialize the cell decomposer
953	cd_sm.setDimensions(domain,c_g,`0`);
954	}
955
956	/! \brief Initialize the grid*
957	*
958	* \param g_sz Global size of the grid
959	* \param bc boundary conditions
960	*
961	*/
962	inline void InitializeDecomposition(const size_t (& g_sz)[dim], const size_t (& bc)[dim], const grid_sm<dim,void> & g_dist = grid_sm<dim,void>())
963	{
964	// fill the global size of the grid
965	for (size_t i = `0` ; i < dim ; i++) {this->g_sz[i] = g_sz[i];}
966
967	// Get the number of processor and calculate the number of sub-domain
968	// for decomposition
969	size_t n_proc = v_cl.getProcessingUnits();
970	size_t n_sub = n_proc * v_sub_unit_factor;
971
972	// Calculate the maximum number (before merging) of sub-domain on
973	// each dimension
974	size_t div[dim];
975	for (size_t i = `0` ; i < dim ; i++)
976	{div[i] = openfpm::math::round_big_2(pow(n_sub,`1.0`/dim));}
977
978	if (g_dist.size(`0`) != `0`)
979	{
980	for (size_t i = `0` ; i < dim ; i++)
981	{div[i] = g_dist.size(i);}
982	}
983
984	// Create the sub-domains
985	dec.setParameters(div,domain,bc,ghost);
986	dec.decompose(dec_options::DEC_SKIP_ICELL);
987	}
988
989	/! \brief Initialize the grid*
990	*
991	* \param g_sz Global size of the grid
992	*
993	*/
994	inline void InitializeStructures(const size_t (& g_sz)[dim])
995	{
996	// an empty
997	openfpm::vector<Box<dim,long int>> empty;
998
999	// Ghost zero
1000	Ghost<dim,long int> zero;
1001
1002	InitializeStructures(g_sz,empty,zero,false);
1003	}
1004
1005	/! \brief Initialize the grid*
1006	*
1007	* \param g_sz Global size of the grid
1008	* \param g ghost extension of the grid in integer unit
1009	* \param bx set of boxes that define where is defined the grid
1010	*
1011	*/
1012	inline void InitializeStructures(const size_t (& g_sz)[dim],
1013	openfpm::vector<Box<dim,long int>> & bx,
1014	const Ghost<dim,long int> & g,
1015	bool use_bx_def)
1016	{
1017	// fill the global size of the grid
1018	for (size_t i = `0` ; i < dim ; i++) {this->g_sz[i] = g_sz[i];}
1019
1020	// Create local grid
1021	Create(bx,g,use_bx_def);
1022	}
1023
1024	// Ghost as integer
1025	Ghost<dim,long int> gint = Ghost<dim,long int>(`0`);
1026
1027	protected:
1028
1029	/! \brief Given a local sub-domain i with a local grid Domain + ghost return the part of the local grid that is domain*
1030	*
1031	* \param i sub-domain
1032	*
1033	* \return the Box defining the domain in the local grid
1034	*
1035	*/
1036	Box<dim,size_t> getDomain(size_t i)
1037	{
1038	return gdb_ext.get(i).Dbox;
1039	}
1040
1041	/! \brief Convert a ghost from grid point units into continus space*
1042	*
1043	* \param gd Ghost in continuous space
1044	* \param cd_sm CellDecomposer of the grid
1045	*
1046	* \return the ghost in continuous unit
1047	*
1048	*/
1049	static inline Ghost<dim,float> convert_ghost(const Ghost<dim,long int> & gd, const CellDecomposer_sm<dim,St,shift<dim,St>> & cd_sm)
1050	{
1051	Ghost<dim,float> gc;
1052
1053	// get the grid spacing
1054	Box<dim,St> sp = cd_sm.getCellBox();
1055
1056	// enlarge 0.001 of the spacing
1057	sp.magnify_fix_P1(`1.1`);
1058
1059	// set the ghost
1060	for (size_t i = `0` ; i < dim ; i++)
1061	{
1062	gc.setLow(i,gd.getLow(i)*(sp.getHigh(i)));
1063	gc.setHigh(i,gd.getHigh(i)*(sp.getHigh(i)));
1064	}
1065
1066	return gc;
1067	}
1068
1069	/! \brief Set the minimum number of sub-domain per processor*
1070	*
1071	* \param n_sub
1072	*
1073	*/
1074	void setDecompositionGranularity(size_t n_sub)
1075	{
1076	this->v_sub_unit_factor = n_sub;
1077	}
1078
1079	void reset_ghost_structures()
1080	{
1081	g_id_to_internal_ghost_box.clear();
1082	ig_box.clear();
1083	init_i_g_box = false;
1084
1085	eg_box.clear();
1086	eb_gid_list.clear();
1087	init_e_g_box = false;
1088
1089	init_local_i_g_box = false;
1090	loc_ig_box.clear();
1091
1092	init_local_e_g_box = false;
1093	loc_eg_box.clear();
1094	}
1095
1096	public:
1097
1098	//! Which kind of grid the structure store
1099	typedef device_grid d_grid;
1100
1101	//! Decomposition used
1102	typedef Decomposition decomposition;
1103
1104	//! value_type
1105	typedef T value_type;
1106
1107	//! Type of space
1108	typedef St stype;
1109
1110	//! Type of Memory
1111	typedef Memory memory_type;
1112
1113	//! Type of device grid
1114	typedef device_grid device_grid_type;
1115
1116	//! Number of dimensions
1117	static const unsigned int dims = dim;
1118
1119	/! \brief Get the domain where the grid is defined*
1120	*
1121	* \return the domain of the grid
1122	*
1123	*/
1124	inline const Box<dim,St> getDomain() const
1125	{
1126	return domain;
1127	}
1128
1129	/! \brief Get the point where it start the origin of the grid of the sub-domain i*
1130	*
1131	* \param i sub-domain
1132	*
1133	* \return the point
1134	*
1135	*/
1136	Point<dim,St> getOffset(size_t i)
1137	{
1138	return pmul(Point<dim,St>(gdb_ext.get(i).origin), cd_sm.getCellBox().getP2()) + getDomain().getP1();
1139	}
1140
1141	/! \brief Get the spacing of the grid in direction i*
1142	*
1143	* \param i dimension
1144	*
1145	* \return the spacing
1146	*
1147	*/
1148	inline St spacing(size_t i) const
1149	{
1150	return cd_sm.getCellBox().getHigh(i);
1151	}
1152
1153	/! \brief Return the total number of points in the grid*
1154	*
1155	* \return number of points
1156	*
1157	*/
1158	size_t size() const
1159	{
1160	return ginfo_v.size();
1161	}
1162
1163	/! \brief set the background value*
1164	*
1165	* You can use this function make sense in case of sparse in case of dense
1166	* it does nothing
1167	*
1168	*/
1169	void setBackgroundValue(T & bv)
1170	{
1171	setBackground_impl<T,decltype(loc_grid)> func(bv,loc_grid);
1172	boost::mpl::for_each_ref< boost::mpl::range_c<int,`0`,T::max_prop>>(func);
1173	}
1174
1175	/! \brief set the background value*
1176	*
1177	* You can use this function make sense in case of sparse in case of dense
1178	* it does nothing
1179	*
1180	*/
1181	template<unsigned int p>
1182	void setBackgroundValue(const typename boost::mpl::at<typename T::type,boost::mpl::int_<p>>::type & bv)
1183	{
1184	for (size_t i = `0` ; i < loc_grid.size() ; i++)
1185	{loc_grid.get(i).template setBackgroundValue<p>(bv);}
1186	}
1187
1188	/! \brief Return the local total number of points inserted in the grid*
1189	*
1190	* in case of dense grid it return the number of local points, in case of
1191	* sparse it return the number of inserted points
1192	*
1193	* \return number of points
1194	*
1195	*/
1196	size_t size_local_inserted() const
1197	{
1198	size_t lins = `0`;
1199
1200	for (size_t i = `0` ; i < loc_grid.size() ; i++)
1201	{lins += loc_grid.get(i).size_inserted();}
1202
1203	return lins;
1204	}
1205
1206	/! \brief Return the total number of points in the grid*
1207	*
1208	* \param i direction
1209	*
1210	* \return number of points on direction i
1211	*
1212	*/
1213	size_t size(size_t i) const
1214	{
1215	return ginfo_v.size(i);
1216	}
1217
1218	/! \brief Copy constructor*
1219	*
1220	* \param g grid to copy
1221	*
1222	*/
1223	grid_dist_id(const grid_dist_id<dim,St,T,Decomposition,Memory,device_grid> & g)
1224	:grid_dist_id_comm<dim,St,T,Decomposition,Memory,device_grid>(g),
1225	domain(g.domain),
1226	ghost(g.ghost),
1227	ghost_int(g.ghost_int),
1228	loc_grid(g.loc_grid),
1229	loc_grid_old(g.loc_grid_old),
1230	dec(g.dec),
1231	gdb_ext_markers(g.gdb_ext_markers),
1232	gdb_ext(g.gdb_ext),
1233	gdb_ext_global(g.gdb_ext_global),
1234	gdb_ext_old(g.gdb_ext_old),
1235	cd_sm(g.cd_sm),
1236	v_cl(g.v_cl),
1237	prp_names(g.prp_names),
1238	g_id_to_external_ghost_box(g.g_id_to_external_ghost_box),
1239	g_id_to_internal_ghost_box(g.g_id_to_internal_ghost_box),
1240	ginfo(g.ginfo),
1241	ginfo_v(g.ginfo_v),
1242	init_local_i_g_box(g.init_local_i_g_box),
1243	init_local_e_g_box(g.init_local_e_g_box)
1244	{
1245	#ifdef SE_CLASS2
1246	check_new(this,`8`,GRID_DIST_EVENT,`4`);
1247	#endif
1248
1249	for (size_t i = `0` ; i < dim ; i++)
1250	{g_sz[i] = g.g_sz[i];}
1251	}
1252
1253	/! \brief This constructor is special, it construct an expanded grid that perfectly overlap with the previous*
1254	*
1255	* The key-word here is "perfectly overlap". Using the default constructor you could create
1256	* something similar, but because of rounding-off error it can happen that it is not perfectly overlapping
1257	*
1258	* \param g previous grid
1259	* \param gh Ghost part in grid units
1260	* \param ext extension of the grid (must be positive on every direction)
1261	*
1262	*/
1263	template<typename H>
1264	grid_dist_id(const grid_dist_id<dim,St,H,typename Decomposition::base_type,Memory,grid_cpu<dim,H>> & g,
1265	const Ghost<dim,long int> & gh,
1266	Box<dim,size_t> ext)
1267	:ghost_int(gh),dec(create_vcluster()),v_cl(create_vcluster())
1268	{
1269	#ifdef SE_CLASS2
1270	check_new(this,`8`,GRID_DIST_EVENT,`4`);
1271	#endif
1272
1273	size_t ext_dim[dim];
1274	for (size_t i = `0` ; i < dim ; i++) {ext_dim[i] = g.getGridInfoVoid().size(i) + ext.getKP1().get(i) + ext.getKP2().get(i);}
1275
1276	// Set the grid info of the extended grid
1277	ginfo.setDimensions(ext_dim);
1278	ginfo_v.setDimensions(ext_dim);
1279
1280	InitializeCellDecomposer(g.getCellDecomposer(),ext);
1281
1282	ghost = convert_ghost(gh,cd_sm);
1283
1284	// Extend the grid by the extension part and calculate the domain
1285
1286	for (size_t i = `0` ; i < dim ; i++)
1287	{
1288	g_sz[i] = g.size(i) + ext.getLow(i) + ext.getHigh(i);
1289
1290	if (g.getDecomposition().periodicity(i) == NON_PERIODIC)
1291	{
1292	this->domain.setLow(i,g.getDomain().getLow(i) - ext.getLow(i) * g.spacing(i) - g.spacing(i) / `2.0`);
1293	this->domain.setHigh(i,g.getDomain().getHigh(i) + ext.getHigh(i) * g.spacing(i) + g.spacing(i) / `2.0`);
1294	}
1295	else
1296	{
1297	this->domain.setLow(i,g.getDomain().getLow(i) - ext.getLow(i) * g.spacing(i));
1298	this->domain.setHigh(i,g.getDomain().getHigh(i) + ext.getHigh(i) * g.spacing(i));
1299	}
1300	}
1301
1302	dec.setParameters(g.getDecomposition(),ghost,this->domain);
1303
1304	// an empty
1305	openfpm::vector<Box<dim,long int>> empty;
1306
1307	InitializeStructures(g.getGridInfoVoid().getSize(),empty,gh,false);
1308	}
1309
1310	/! It constructs a grid of a specified size, defined on a specified Box space, forcing to follow a specified decomposition and with a specified ghost size*
1311	*
1312	* \param dec Decomposition
1313	* \param g_sz grid size on each dimension
1314	* \param ghost Ghost part
1315	*
1316	*/
1317	grid_dist_id(const Decomposition & dec,
1318	const size_t (& g_sz)[dim],
1319	const Ghost<dim,St> & ghost)
1320	:domain(dec.getDomain()),ghost(ghost),ghost_int(INVALID_GHOST),dec(dec),v_cl(create_vcluster()),
1321	ginfo(g_sz),ginfo_v(g_sz)
1322	{
1323	#ifdef SE_CLASS2
1324	check_new(this,`8`,GRID_DIST_EVENT,`4`);
1325	#endif
1326
1327	InitializeCellDecomposer(g_sz,dec.periodicity());
1328
1329	this->dec = dec.duplicate(ghost);
1330
1331	InitializeStructures(g_sz);
1332	}
1333
1334	/! It constructs a grid of a specified size, defined on a specified Box space, forcing to follow a specified decomposition and with a specified ghost size*
1335	*
1336	* \param dec Decomposition
1337	* \param g_sz grid size on each dimension
1338	* \param ghost Ghost part
1339	*
1340	*/
1341	grid_dist_id(Decomposition && dec, const size_t (& g_sz)[dim],
1342	const Ghost<dim,St> & ghost)
1343	:domain(dec.getDomain()),ghost(ghost),dec(dec),ginfo(g_sz),
1344	ginfo_v(g_sz),v_cl(create_vcluster()),ghost_int(INVALID_GHOST)
1345	{
1346	#ifdef SE_CLASS2
1347	check_new(this,`8`,GRID_DIST_EVENT,`4`);
1348	#endif
1349
1350	InitializeCellDecomposer(g_sz,dec.periodicity());
1351
1352	this->dec = dec.duplicate(ghost);
1353
1354	InitializeStructures(g_sz);
1355	}
1356
1357	/! It constructs a grid of a specified size, defined on a specified Box space, forcing to follow a specified decomposition, and having a specified ghost size*
1358	*
1359	* \param dec Decomposition
1360	* \param g_sz grid size on each dimension
1361	* \param g Ghost part (given in grid units)
1362	*
1363	* \warning In very rare case the ghost part can be one point bigger than the one specified
1364	*
1365	*/
1366	grid_dist_id(const Decomposition & dec, const size_t (& g_sz)[dim],
1367	const Ghost<dim,long int> & g)
1368	:domain(dec.getDomain()),ghost_int(g),dec(create_vcluster()),v_cl(create_vcluster()),
1369	ginfo(g_sz),ginfo_v(g_sz)
1370	{
1371	#ifdef SE_CLASS2
1372	check_new(this,`8`,GRID_DIST_EVENT,`4`);
1373	#endif
1374
1375	InitializeCellDecomposer(g_sz,dec.periodicity());
1376
1377	ghost = convert_ghost(g,cd_sm);
1378	this->dec = dec.duplicate(ghost);
1379
1380	// an empty
1381	openfpm::vector<Box<dim,long int>> empty;
1382
1383	// Initialize structures
1384	InitializeStructures(g_sz,empty,g,false);
1385	}
1386
1387	/! It construct a grid of a specified size, defined on a specified Box space, forcing to follow a specified decomposition, and having a specified ghost size*
1388	*
1389	* \param dec Decomposition
1390	* \param g_sz grid size on each dimension
1391	* \param g Ghost part (given in grid units)
1392	*
1393	* \warning In very rare case the ghost part can be one point bigger than the one specified
1394	*
1395	*/
1396	grid_dist_id(Decomposition && dec, const size_t (& g_sz)[dim],
1397	const Ghost<dim,long int> & g)
1398	:domain(dec.getDomain()),dec(dec),v_cl(create_vcluster()),ginfo(g_sz),
1399	ginfo_v(g_sz),ghost_int(g)
1400	{
1401	#ifdef SE_CLASS2
1402	check_new(this,`8`,GRID_DIST_EVENT,`4`);
1403	#endif
1404	InitializeCellDecomposer(g_sz,dec.periodicity());
1405
1406	ghost = convert_ghost(g,cd_sm);
1407	this->dec = dec.duplicate(ghost);
1408
1409	// an empty
1410	openfpm::vector<Box<dim,long int>> empty;
1411
1412	// Initialize structures
1413	InitializeStructures(g_sz,empty,g,false);
1414	}
1415
1416	/! It construct a grid of a specified size, defined on a specified Box space, and having a specified ghost size*
1417	*
1418	* \param g_sz grid size on each dimension
1419	* \param domain Box that contain the grid
1420	* \param g Ghost part (given in grid units)
1421	*
1422	* \warning In very rare case the ghost part can be one point bigger than the one specified
1423	*
1424	*/
1425	grid_dist_id(const size_t (& g_sz)[dim],const Box<dim,St> & domain, const Ghost<dim,St> & g, size_t opt = `0`)
1426	:grid_dist_id(g_sz,domain,g,create_non_periodic<dim>(),opt)
1427	{
1428	}
1429
1430	/! It construct a grid of a specified size, defined on a specified Box space, having a specified ghost size and periodicity*
1431	*
1432	* \param g_sz grid size on each dimension
1433	* \param domain Box that contain the grid
1434	* \param g Ghost part of the domain (given in grid units)
1435	*
1436	* \warning In very rare case the ghost part can be one point bigger than the one specified
1437	*
1438	*/
1439	grid_dist_id(const size_t (& g_sz)[dim],const Box<dim,St> & domain, const Ghost<dim,long int> & g, size_t opt = `0`)
1440	:grid_dist_id(g_sz,domain,g,create_non_periodic<dim>(),opt)
1441	{
1442	}
1443
1444	/! It construct a grid of a specified size, defined on a specified Box space, having a specified ghost size, and specified periodicity*
1445	*
1446	* \param g_sz grid size on each dimension
1447	* \param domain Box that contain the grid
1448	* \param g Ghost part (given in grid units)
1449	* \param p Boundary conditions
1450	*
1451	* \warning In very rare case the ghost part can be one point bigger than the one specified
1452	*
1453	*/
1454	grid_dist_id(const size_t (& g_sz)[dim],const Box<dim,St> & domain, const Ghost<dim,St> & g, const periodicity<dim> & p, size_t opt = `0`)
1455	:domain(domain),ghost(g),ghost_int(INVALID_GHOST),dec(create_vcluster()),v_cl(create_vcluster()),ginfo(g_sz),ginfo_v(g_sz)
1456	{
1457	#ifdef SE_CLASS2
1458	check_new(this,`8`,GRID_DIST_EVENT,`4`);
1459	#endif
1460
1461	if (opt >> `32` != `0`)
1462	{this->setDecompositionGranularity(opt >> `32`);}
1463
1464	check_domain(domain);
1465	InitializeCellDecomposer(g_sz,p.bc);
1466	InitializeDecomposition(g_sz, p.bc);
1467	InitializeStructures(g_sz);
1468	}
1469
1470
1471	/! It construct a grid of a specified size, defined on a specified Box space, having a specified ghost size and periodicity*
1472	*
1473	* \param g_sz grid size on each dimension
1474	* \param domain Box that contain the grid
1475	* \param g Ghost part of the domain (given in grid units)
1476	* \param p periodicity
1477	*
1478	* \warning In very rare case the ghost part can be one point bigger than the one specified
1479	*
1480	*/
1481	grid_dist_id(const size_t (& g_sz)[dim],const Box<dim,St> & domain, const Ghost<dim,long int> & g, const periodicity<dim> & p, size_t opt = `0`, const grid_sm<dim,void> & g_dec = grid_sm<dim,void>())
1482	:domain(domain),ghost_int(g),dec(create_vcluster()),v_cl(create_vcluster()),ginfo(g_sz),ginfo_v(g_sz)
1483	{
1484	#ifdef SE_CLASS2
1485	check_new(this,`8`,GRID_DIST_EVENT,`4`);
1486	#endif
1487
1488	if (opt >> `32` != `0`)
1489	{this->setDecompositionGranularity(opt >> `32`);}
1490
1491	check_domain(domain);
1492	InitializeCellDecomposer(g_sz,p.bc);
1493
1494	ghost = convert_ghost(g,cd_sm);
1495
1496	InitializeDecomposition(g_sz,p.bc,g_dec);
1497
1498	// an empty
1499	openfpm::vector<Box<dim,long int>> empty;
1500
1501	// Initialize structures
1502	InitializeStructures(g_sz,empty,g,false);
1503	}
1504
1505
1506	/! \brief It construct a grid on the full domain restricted*
1507	* to the set of boxes specified
1508	*
1509	* In particular the grid is defined in the space equal to the
1510	* domain intersected the boxes defined by bx
1511	*
1512	* \param g_sz grid size on each dimension
1513	* \param domain where the grid is constructed
1514	* \param g ghost size
1515	* \param p periodicity of the grid
1516	* \param bx set of boxes where the grid is defined
1517	*
1518	*
1519	*/
1520	grid_dist_id(const size_t (& g_sz)[dim],
1521	const Box<dim,St> & domain,
1522	const Ghost<dim,long int> & g,
1523	const periodicity<dim> & p,
1524	openfpm::vector<Box<dim,long int>> & bx_def)
1525	:domain(domain),dec(create_vcluster()),v_cl(create_vcluster()),ginfo(g_sz),ginfo_v(g_sz),gint(g)
1526	{
1527	#ifdef SE_CLASS2
1528	check_new(this,`8`,GRID_DIST_EVENT,`4`);
1529	#endif
1530
1531	check_domain(domain);
1532	InitializeCellDecomposer(g_sz,p.bc);
1533
1534	ghost = convert_ghost(g,cd_sm);
1535
1536	InitializeDecomposition(g_sz, p.bc);
1537	InitializeStructures(g_sz,bx_def,g,true);
1538	this->bx_def = bx_def;
1539	this->use_bx_def = true;
1540	}
1541
1542	/! \brief Get an object containing the grid informations*
1543	*
1544	* \return an information object about this grid
1545	*
1546	*/
1547	const grid_sm<dim,T> & getGridInfo() const
1548	{
1549	#ifdef SE_CLASS2
1550	check_valid(this,`8`);
1551	#endif
1552	return ginfo;
1553	}
1554
1555	/! \brief Get an object containing the grid informations without type*
1556	*
1557	* \return an information object about this grid
1558	*
1559	*/
1560	const grid_sm<dim,void> & getGridInfoVoid() const
1561	{
1562	#ifdef SE_CLASS2
1563	check_valid(this,`8`);
1564	#endif
1565	return ginfo_v;
1566	}
1567
1568	/! \brief Get the object that store the information about the decomposition*
1569	*
1570	* \return the decomposition object
1571	*
1572	*/
1573	Decomposition & getDecomposition()
1574	{
1575	#ifdef SE_CLASS2
1576	check_valid(this,`8`);
1577	#endif
1578	return dec;
1579	}
1580
1581	/! \brief Get the object that store the information about the decomposition*
1582	*
1583	* \return the decomposition object
1584	*
1585	*/
1586	const Decomposition & getDecomposition() const
1587	{
1588	#ifdef SE_CLASS2
1589	check_valid(this,`8`);
1590	#endif
1591	return dec;
1592	}
1593
1594	/! \brief Return the cell decomposer*
1595	*
1596	* \return the cell decomposer
1597	*
1598	*/
1599	const CellDecomposer_sm<dim,St,shift<dim,St>> & getCellDecomposer() const
1600	{
1601	#ifdef SE_CLASS2
1602	check_valid(this,`8`);
1603	#endif
1604	return cd_sm;
1605	}
1606
1607	/! \brief Check that the global grid key is inside the grid domain*
1608	*
1609	* \param gk point to check
1610	*
1611	* \return true if is inside
1612	*
1613	*/
1614	bool isInside(const grid_key_dx<dim> & gk) const
1615	{
1616	#ifdef SE_CLASS2
1617	check_valid(this,`8`);
1618	#endif
1619	for (size_t i = `0` ; i < dim ; i++)
1620	{
1621	if (gk.get(i) < `0` \|\| gk.get(i) >= (long int)g_sz[i])
1622	{return false;}
1623	}
1624
1625	return true;
1626	}
1627
1628	/! \brief Get the total number of grid points for the calling processor*
1629	*
1630	* \return The number of grid points
1631	*
1632	*/
1633	size_t getLocalDomainSize() const
1634	{
1635	#ifdef SE_CLASS2
1636	check_valid(this,`8`);
1637	#endif
1638	size_t total = `0`;
1639
1640	for (size_t i = `0` ; i < gdb_ext.size() ; i++)
1641	{
1642	total += gdb_ext.get(i).Dbox.getVolumeKey();
1643	}
1644
1645	return total;
1646	}
1647
1648	/! \brief Get the total number of grid points with ghost for the calling processor*
1649	*
1650	* \return The number of grid points
1651	*
1652	*/
1653	size_t getLocalDomainWithGhostSize() const
1654	{
1655	#ifdef SE_CLASS2
1656	check_valid(this,`8`);
1657	#endif
1658	size_t total = `0`;
1659
1660	for (size_t i = `0` ; i < gdb_ext.size() ; i++)
1661	{
1662	total += gdb_ext.get(i).GDbox.getVolumeKey();
1663	}
1664
1665	return total;
1666	}
1667
1668
1669	/! \brief It return the informations about the local grids*
1670	*
1671	* \return The information about the local grids
1672	*
1673	*/
1674	const openfpm::vector<GBoxes<device_grid::dims>> & getLocalGridsInfo()
1675	{
1676	#ifdef SE_CLASS2
1677	check_valid(this,`8`);
1678	#endif
1679	return gdb_ext;
1680	}
1681
1682	/! \brief It gathers the information about local grids for all of the processors*
1683	*
1684	* \param gdb_ext_global where to store the grid infos
1685	*
1686	*/
1687	void getGlobalGridsInfo(openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext_global) const
1688	{
1689	#ifdef SE_CLASS2
1690	check_valid(this,`8`);
1691	#endif
1692	gdb_ext_global.clear();
1693
1694	v_cl.SGather(gdb_ext,gdb_ext_global,`0`);
1695	v_cl.execute();
1696
1697	size_t size_r;
1698	size_t size = gdb_ext_global.size();
1699
1700	if (v_cl.getProcessUnitID() == `0`)
1701	{
1702	for (size_t i = `1`; i < v_cl.getProcessingUnits(); i++)
1703	v_cl.send(i,`0`,&size,sizeof(size_t));
1704
1705	size_r = size;
1706	}
1707	else
1708	v_cl.recv(`0`,`0`,&size_r,sizeof(size_t));
1709
1710	v_cl.execute();
1711
1712	gdb_ext_global.resize(size_r);
1713
1714
1715	if (v_cl.getProcessUnitID() == `0`)
1716	{
1717	for (size_t i = `1`; i < v_cl.getProcessingUnits(); i++)
1718	v_cl.send(i,`0`,gdb_ext_global);
1719	}
1720	else
1721	v_cl.recv(`0`,`0`,gdb_ext_global);
1722
1723	v_cl.execute();
1724	}
1725
1726
1727	/! \brief It return an iterator that span the full grid domain (each processor span its local domain)*
1728	*
1729	* \return the iterator
1730	*
1731	*/
1732	grid_dist_iterator<dim,device_grid,
1733	decltype(device_grid::type_of_subiterator()),
1734	FREE>
1735	getOldDomainIterator() const
1736	{
1737	#ifdef SE_CLASS2
1738	check_valid(this,`8`);
1739	#endif
1740
1741	grid_key_dx<dim> stop(ginfo_v.getSize());
1742	grid_key_dx<dim> one;
1743	one.one();
1744	stop = stop - one;
1745
1746	grid_dist_iterator<dim,device_grid,
1747	decltype(device_grid::type_of_subiterator()),
1748	FREE> it(loc_grid_old,gdb_ext_old,stop);
1749
1750	return it;
1751	}
1752
1753	/! /brief Get a grid Iterator*
1754	*
1755	* In case of dense grid getGridIterator is equivalent to getDomainIterator
1756	* in case of sparse grid getDomainIterator go across all the
1757	* inserted point get grid iterator run across all grid points independently
1758	* that the point has been insert or not
1759	*
1760	* \return a Grid iterator
1761	*
1762	*/
1763	inline grid_dist_id_iterator_dec<Decomposition> getGridIterator(const grid_key_dx<dim> & start, const grid_key_dx<dim> & stop)
1764	{
1765	grid_dist_id_iterator_dec<Decomposition> it_dec(getDecomposition(), g_sz, start, stop);
1766	return it_dec;
1767	}
1768
1769	#ifdef __NVCC__
1770
1771	/! \brief Insert point in the grid*
1772	*
1773	* \param f1 lambda function to insert point
1774	* \param f2 lambda function to set points
1775	*/
1776	template<typename lambda_t1, typename lambda_t2>
1777	void addPoints(lambda_t1 f1, lambda_t2 f2)
1778	{
1779	auto it = getGridIteratorGPU();
1780	it.setGPUInsertBuffer(`1`);
1781
1782	it.template launch<`1`>(launch_insert_sparse(),f1,f2);
1783	}
1784
1785	/! \brief Insert point in the grid between start and stop*
1786	*
1787	* \param start point
1788	* \param stop point
1789	* \param f1 lambda function to insert point
1790	* \param f2 lambda function to set points
1791	*/
1792	template<typename lambda_t1, typename lambda_t2>
1793	void addPoints(grid_key_dx<dim> k1, grid_key_dx<dim> k2, lambda_t1 f1, lambda_t2 f2)
1794	{
1795	auto it = getGridIteratorGPU(k1,k2);
1796	it.setGPUInsertBuffer(`1`);
1797
1798	it.template launch<`1`>(launch_insert_sparse(),f1,f2);
1799	}
1800
1801	/! /brief Get a grid Iterator in GPU*
1802	*
1803	* In case of dense grid getGridIterator is equivalent to getDomainIteratorGPU
1804	* in case of sparse distributed grid getDomainIterator go across all the
1805	* inserted point getGridIteratorGPU run across all grid points independently
1806	* that the point has been insert or not
1807	*
1808	* \param start point
1809	* \param stop point
1810	*
1811	* \return a Grid iterator
1812	*
1813	*/
1814	inline grid_dist_id_iterator_gpu<Decomposition,openfpm::vector<device_grid>>
1815	getGridIteratorGPU(const grid_key_dx<dim> & start, const grid_key_dx<dim> & stop)
1816	{
1817	grid_dist_id_iterator_gpu<Decomposition,openfpm::vector<device_grid>> it_dec(loc_grid,getDecomposition(), g_sz, start, stop);
1818	return it_dec;
1819	}
1820
1821	/! /brief Get a grid Iterator in GPU*
1822	*
1823	* In case of dense grid getGridIterator is equivalent to getDomainIteratorGPU
1824	* in case of sparse distributed grid getDomainIterator go across all the
1825	* inserted point getGridIteratorGPU run across all grid points independently
1826	* that the point has been insert or not
1827	*
1828	* \return a Grid iterator
1829	*
1830	*/
1831	inline grid_dist_id_iterator_gpu<Decomposition,openfpm::vector<device_grid>>
1832	getGridIteratorGPU()
1833	{
1834	grid_dist_id_iterator_gpu<Decomposition,openfpm::vector<device_grid>> it_dec(loc_grid,getDecomposition(), g_sz);
1835	return it_dec;
1836	}
1837
1838	#endif
1839
1840	/! /brief Get a grid Iterator running also on ghost area*
1841	*
1842	* In case of dense grid getGridIterator is equivalent to getDomainIterator
1843	* in case if sparse distributed grid getDomainIterator go across all the
1844	* inserted point get grid iterator run across all grid points independently
1845	* that the point has been insert or not
1846	*
1847	* \return a Grid iterator
1848	*
1849	*/
1850	inline grid_dist_id_iterator_dec<Decomposition,true> getGridGhostIterator(const grid_key_dx<dim> & start, const grid_key_dx<dim> & stop)
1851	{
1852	grid_dist_id_iterator_dec<Decomposition,true> it_dec(getDecomposition(), g_sz, start, stop);
1853	return it_dec;
1854	}
1855
1856	/! /brief Get a grid Iterator*
1857	*
1858	* In case of dense grid getGridIterator is equivalent to getDomainIterator
1859	* in case if sparse distributed grid getDomainIterator go across all the
1860	* inserted point get grid iterator run across all grid points independently
1861	* that the point has been insert or not
1862	*
1863	* \return a Grid iterator
1864	*
1865	*/
1866	inline grid_dist_id_iterator_dec<Decomposition> getGridIterator()
1867	{
1868	grid_key_dx<dim> start;
1869	grid_key_dx<dim> stop;
1870	for (size_t i = `0`; i < dim; i++)
1871	{
1872	start.set_d(i, `0`);
1873	stop.set_d(i, g_sz[i] - `1`);
1874	}
1875
1876	grid_dist_id_iterator_dec<Decomposition> it_dec(getDecomposition(), g_sz, start, stop);
1877	return it_dec;
1878	}
1879
1880	/! \brief It return an iterator that span the full grid domain (each processor span its local domain)*
1881	*
1882	* \return the iterator
1883	*
1884	*/
1885	grid_dist_iterator<dim,device_grid,
1886	decltype(device_grid::type_of_subiterator()),FREE>
1887	getDomainIterator() const
1888	{
1889	#ifdef SE_CLASS2
1890	check_valid(this,`8`);
1891	#endif
1892
1893	grid_key_dx<dim> stop(ginfo_v.getSize());
1894	grid_key_dx<dim> one;
1895	one.one();
1896	stop = stop - one;
1897
1898	grid_dist_iterator<dim,device_grid,
1899	decltype(device_grid::type_of_subiterator()),
1900	FREE> it(loc_grid,gdb_ext,stop);
1901
1902	return it;
1903	}
1904
1905	/! \brief It return an iterator that span the full grid domain (each processor span its local domain)*
1906	*
1907	* \param stencil_pnt stencil points
1908	*
1909	* \return the iterator
1910	*
1911	*/
1912	template<unsigned int Np>
1913	grid_dist_iterator<dim,device_grid,
1914	decltype(device_grid::template type_of_subiterator<stencil_offset_compute<dim,Np>>()),
1915	FREE,
1916	stencil_offset_compute<dim,Np> >
1917	getDomainIteratorStencil(const grid_key_dx<dim> (& stencil_pnt)[Np]) const
1918	{
1919	#ifdef SE_CLASS2
1920	check_valid(this,`8`);
1921	#endif
1922
1923	grid_key_dx<dim> stop(ginfo_v.getSize());
1924	grid_key_dx<dim> one;
1925	one.one();
1926	stop = stop - one;
1927
1928	grid_dist_iterator<dim,device_grid,
1929	decltype(device_grid::template type_of_subiterator<stencil_offset_compute<dim,Np>>()),
1930	FREE,
1931	stencil_offset_compute<dim,Np>> it(loc_grid,gdb_ext,stop,stencil_pnt);
1932
1933	return it;
1934	}
1935
1936	/! \brief It return an iterator that span the grid domain + ghost part*
1937	*
1938	* \return the iterator
1939	*
1940	*/
1941	grid_dist_iterator<dim,device_grid,
1942	decltype(device_grid::type_of_iterator()),
1943	FIXED>
1944	getDomainGhostIterator() const
1945	{
1946	#ifdef SE_CLASS2
1947	check_valid(this,`8`);
1948	#endif
1949	grid_key_dx<dim> stop;
1950	for (size_t i = `0` ; i < dim ; i++)
1951	{stop.set_d(i,`0`);}
1952
1953	grid_dist_iterator<dim,device_grid,
1954	decltype(device_grid::type_of_iterator()),
1955	FIXED> it(loc_grid,gdb_ext,stop);
1956
1957	return it;
1958	}
1959
1960	/! \brief It return an iterator that span the grid domain only in the specified*
1961	* part
1962	*
1963	* The key spanned are the one inside the box spanned by the start point and the end
1964	* point included
1965	*
1966	* \param start point
1967	* \param stop point
1968	*
1969	* \return the sub-domain iterator
1970	*
1971	*/
1972	grid_dist_iterator_sub<dim,device_grid>
1973	getSubDomainIterator(const grid_key_dx<dim> & start,
1974	const grid_key_dx<dim> & stop) const
1975	{
1976	#ifdef SE_CLASS2
1977	check_valid(this,`8`);
1978	#endif
1979	grid_dist_iterator_sub<dim,device_grid> it(start,stop,loc_grid,gdb_ext);
1980
1981	return it;
1982	}
1983
1984	/! \brief It return an iterator that span the grid domain only in the specified*
1985	* part
1986	*
1987	* The key spanned are the one inside the box spanned by the start point and the end
1988	* point included
1989	*
1990	* \param start point
1991	* \param stop point
1992	*
1993	* \return an iterator on the sub-part of the grid
1994	*
1995	*/
1996	grid_dist_iterator_sub<dim,device_grid> getSubDomainIterator(const long int (& start)[dim], const long int (& stop)[dim]) const
1997	{
1998	grid_dist_iterator_sub<dim,device_grid> it(grid_key_dx<dim>(start),grid_key_dx<dim>(stop),loc_grid,gdb_ext);
1999
2000	return it;
2001	}
2002
2003	//! Destructor
2004	~grid_dist_id()
2005	{
2006	#ifdef SE_CLASS2
2007	check_delete(this);
2008	#endif
2009	dec.decRef();
2010	}
2011
2012	/! \brief Get the Virtual Cluster machine*
2013	*
2014	* \return the Virtual cluster machine
2015	*
2016	*/
2017	Vcluster<> & getVC()
2018	{
2019	#ifdef SE_CLASS2
2020	check_valid(this,`8`);
2021	#endif
2022	return v_cl;
2023	}
2024
2025	/! \brief Eliminate many internal temporary buffer you can use this between flushes if you get some out of memory*
2026	*
2027	*
2028	*/
2029	void removeUnusedBuffers()
2030	{
2031	for (int i = `0` ; i < loc_grid.size() ; i++)
2032	{
2033	loc_grid.get(i).removeUnusedBuffers();
2034	}
2035	}
2036
2037	/! \brief Indicate that this grid is not staggered*
2038	*
2039	* \return false
2040	*
2041	*/
2042	bool is_staggered()
2043	{
2044	return false;
2045	}
2046
2047	/! \brief remove an element in the grid*
2048	*
2049	* In case of dense grid this function print a warning, in case of sparse
2050	* grid this function remove a grid point.
2051	*
2052	* \param v1 grid_key that identify the element in the grid
2053	*
2054	* \return a reference to the inserted element
2055	*
2056	*/
2057	template <typename bg_key> inline void remove(const grid_dist_key_dx<dim,bg_key> & v1)
2058	{
2059	#ifdef SE_CLASS2
2060	check_valid(this,`8`);
2061	#endif
2062	return loc_grid.get(v1.getSub()).remove(v1.getKey());
2063	}
2064
2065	/! \brief remove an element in the grid*
2066	*
2067	* In case of dense grid this function print a warning, in case of sparse
2068	* grid this function remove a grid point.
2069	*
2070	* \param v1 grid_key that identify the element in the grid
2071	*
2072	* \return a reference to the inserted element
2073	*
2074	*/
2075	template <typename bg_key> inline void remove_no_flush(const grid_dist_key_dx<dim,bg_key> & v1)
2076	{
2077	#ifdef SE_CLASS2
2078	check_valid(this,`8`);
2079	#endif
2080	return loc_grid.get(v1.getSub()).remove_no_flush(v1.getKey());
2081	}
2082
2083
2084	template<typename ... v_reduce>
2085	void flush(flush_type opt = flush_type::FLUSH_ON_HOST)
2086	{
2087	for (size_t i = `0` ; i < loc_grid.size() ; i++)
2088	{
2089	loc_grid.get(i).template flush<v_reduce ...>(v_cl.getmgpuContext(),opt);
2090	}
2091	}
2092
2093	/! \brief remove an element in the grid*
2094	*
2095	* In case of dense grid this function print a warning, in case of sparse
2096	* grid this function remove a grid point.
2097	*
2098	* \param v1 grid_key that identify the element in the grid
2099	*
2100	* \return a reference to the inserted element
2101	*
2102	*/
2103	inline void flush_remove()
2104	{
2105	#ifdef SE_CLASS2
2106	check_valid(this,`8`);
2107	#endif
2108	for (size_t i = `0` ; i < loc_grid.size() ; i++)
2109	{loc_grid.get(i).flush_remove();}
2110	}
2111
2112	/! \brief insert an element in the grid*
2113	*
2114	* In case of dense grid this function is equivalent to get, in case of sparse
2115	* grid this function insert a grid point. When the point already exist it return
2116	* a reference to the already existing point. In case of massive insert Sparse grids
2117	* it give a reference to the inserted element in the insert buffer
2118	*
2119	* \tparam p property to get (is an integer)
2120	* \param v1 grid_key that identify the element in the grid
2121	*
2122	* \return a reference to the inserted element
2123	*
2124	*/
2125	template <unsigned int p,typename bg_key>inline auto insert(const grid_dist_key_dx<dim,bg_key> & v1)
2126	-> decltype(loc_grid.get(v1.getSub()).template insert<p>(v1.getKey()))
2127	{
2128	#ifdef SE_CLASS2
2129	check_valid(this,`8`);
2130	#endif
2131
2132	return loc_grid.get(v1.getSub()).template insert<p>(v1.getKey());
2133	}
2134
2135
2136	/! \brief insert an element in the grid*
2137	*
2138	* In case of dense grid this function is equivalent to get, in case of sparse
2139	* grid this function insert a grid point. When the point already exist it return
2140	* a reference to the already existing point. In case of massive insert Sparse grids
2141	* The point is inserted immediately and a reference to the inserted element is returned
2142	*
2143	* \warning This function is not fast an unlucky insert can potentially cost O(N) where N is the number
2144	* of points (worst case)
2145	*
2146	* \tparam p property to get (is an integer)
2147	* \param v1 grid_key that identify the element in the grid
2148	*
2149	* \return a reference to the inserted element
2150	*
2151	*/
2152	template <unsigned int p,typename bg_key>inline auto insertFlush(const grid_dist_key_dx<dim,bg_key> & v1)
2153	-> decltype(loc_grid.get(v1.getSub()).template insertFlush<p>(v1.getKey()))
2154	{
2155	#ifdef SE_CLASS2
2156	check_valid(this,`8`);
2157	#endif
2158
2159	return loc_grid.get(v1.getSub()).template insertFlush<p>(v1.getKey());
2160	}
2161
2162	/! \brief Get the reference of the selected element*
2163	*
2164	* \tparam p property to get (is an integer)
2165	* \param v1 grid_key that identify the element in the grid
2166	*
2167	* \return the selected element
2168	*
2169	*/
2170	template <unsigned int p, typename bg_key>
2171	inline auto get(const grid_dist_key_dx<dim,bg_key> & v1) const
2172	-> typename std::add_lvalue_reference<decltype(loc_grid.get(v1.getSub()).template get<p>(v1.getKey()))>::type
2173	{
2174	#ifdef SE_CLASS2
2175	check_valid(this,`8`);
2176	#endif
2177	return loc_grid.get(v1.getSub()).template get<p>(v1.getKey());
2178	}
2179
2180
2181	/! \brief Get the reference of the selected element*
2182	*
2183	* \tparam p property to get (is an integer)
2184	* \param v1 grid_key that identify the element in the grid
2185	*
2186	* \return the selected element
2187	*
2188	*/
2189	template <unsigned int p, typename bg_key>
2190	inline auto get(const grid_dist_key_dx<dim,bg_key> & v1)
2191	-> decltype(loc_grid.get(v1.getSub()).template get<p>(v1.getKey()))
2192	{
2193	#ifdef SE_CLASS2
2194	check_valid(this,`8`);
2195	#endif
2196	return loc_grid.get(v1.getSub()).template get<p>(v1.getKey());
2197	}
2198
2199	/! \brief Get the reference of the selected element*
2200	*
2201	* \tparam p property to get (is an integer)
2202	* \param v1 grid_key that identify the element in the grid
2203	*
2204	* \return the selected element
2205	*
2206	*/
2207	template <unsigned int p = `0`>
2208	inline auto get(const grid_dist_g_dx<device_grid> & v1) const
2209	-> decltype(v1.getSub()->template get<p>(v1.getKey()))
2210	{
2211	#ifdef SE_CLASS2
2212	check_valid(this,`8`);
2213	#endif
2214	return v1.getSub()->template get<p>(v1.getKey());
2215	}
2216
2217	/! \brief Get the reference of the selected element*
2218	*
2219	* \tparam p property to get (is an integer)
2220	* \param v1 grid_key that identify the element in the grid
2221	*
2222	* \return the selected element
2223	*
2224	*/
2225	template <unsigned int p = `0`>
2226	inline auto get(const grid_dist_g_dx<device_grid> & v1) -> decltype(v1.getSub()->template get<p>(v1.getKey()))
2227	{
2228	#ifdef SE_CLASS2
2229	check_valid(this,`8`);
2230	#endif
2231	return v1.getSub()->template get<p>(v1.getKey());
2232	}
2233
2234	/! \brief Get the reference of the selected element*
2235	*
2236	* \tparam p property to get (is an integer)
2237	* \param v1 grid_key that identify the element in the grid
2238	*
2239	* \return the selected element
2240	*
2241	*/
2242	template <unsigned int p = `0`>
2243	inline auto get(const grid_dist_lin_dx & v1) const -> decltype(loc_grid.get(v1.getSub()).template get<p>(v1.getKey()))
2244	{
2245	#ifdef SE_CLASS2
2246	check_valid(this,`8`);
2247	#endif
2248	return loc_grid.get(v1.getSub()).template get<p>(v1.getKey());
2249	}
2250
2251	/! \brief Get the reference of the selected element*
2252	*
2253	* \tparam p property to get (is an integer)
2254	* \param v1 grid_key that identify the element in the grid
2255	*
2256	* \return the selected element
2257	*
2258	*/
2259	template <unsigned int p = `0`>
2260	inline auto get(const grid_dist_lin_dx & v1) -> decltype(loc_grid.get(v1.getSub()).template get<p>(v1.getKey()))
2261	{
2262	#ifdef SE_CLASS2
2263	check_valid(this,`8`);
2264	#endif
2265	return loc_grid.get(v1.getSub()).template get<p>(v1.getKey());
2266	}
2267
2268	/! \brief Get the reference of the selected element*
2269	*
2270	* \tparam p property to get (is an integer)
2271	* \param v1 grid_key that identify the element in the grid
2272	*
2273	* \return the selected element
2274	*
2275	*/
2276	template <typename bg_key>
2277	inline Point<dim,St> getPos(const grid_dist_key_dx<dim,bg_key> & v1)
2278	{
2279	#ifdef SE_CLASS2
2280	check_valid(this,`8`);
2281	#endif
2282	Point<dim,St> p;
2283
2284	for (int i = `0` ; i < dim ; i++)
2285	{
2286	p.get(i) = (gdb_ext.get(v1.getSub()).origin.get(i) + v1.getKeyRef().get(i)) * this->spacing(i);
2287	}
2288
2289	return p;
2290	}
2291
2292	/! \brief Check if the point exist*
2293	*
2294	* \param v1 grid_key that identify the element in the grid
2295	*
2296	* \return the true if the point exist
2297	*
2298	*/
2299	template<typename bg_key>
2300	inline bool existPoint(const grid_dist_key_dx<dim,bg_key> & v1) const
2301	{
2302	return loc_grid.get(v1.getSub()).existPoint(v1.getKey());
2303	}
2304
2305	/! \brief Get the reference of the selected element*
2306	*
2307	* \tparam p property to get (is an integer)
2308	* \param v1 grid_key that identify the element in the grid
2309	*
2310	* \return the selected element
2311	*
2312	*/
2313	template <unsigned int p = `0`>
2314	inline auto getProp(const grid_dist_key_dx<dim> & v1) const -> decltype(this->template get<p>(v1))
2315	{
2316	return this->template get<p>(v1);
2317	}
2318
2319	/! \brief Get the reference of the selected element*
2320	*
2321	* \tparam p property to get (is an integer)
2322	* \param v1 grid_key that identify the element in the grid
2323	*
2324	* \return the selected element
2325	*
2326	*/
2327	template <unsigned int p = `0`>
2328	inline auto getProp(const grid_dist_key_dx<dim> & v1) -> decltype(this->template get<p>(v1))
2329	{
2330	return this->template get<p>(v1);
2331	}
2332
2333	/! \brief It synchronize the ghost parts*
2334	*
2335	* \tparam prp... Properties to synchronize
2336	*
2337	*/
2338	template<int... prp> void ghost_get(size_t opt = `0`)
2339	{
2340	#ifdef SE_CLASS2
2341	check_valid(this,`8`);
2342	#endif
2343
2344	// Convert the ghost internal boxes into grid unit boxes
2345	create_ig_box();
2346
2347	// Convert the ghost external boxes into grid unit boxes
2348	create_eg_box();
2349
2350	// Convert the local ghost internal boxes into grid unit boxes
2351	create_local_ig_box();
2352
2353	// Convert the local external ghost boxes into grid unit boxes
2354	create_local_eg_box();
2355
2356	grid_dist_id_comm<dim,St,T,Decomposition,Memory,device_grid>::template ghost_get_<prp...>(ig_box,
2357	eg_box,
2358	loc_ig_box,
2359	loc_eg_box,
2360	gdb_ext,
2361	eb_gid_list,
2362	use_bx_def,
2363	loc_grid,
2364	ginfo_v,
2365	g_id_to_external_ghost_box,
2366	opt);
2367	}
2368
2369	/! \brief It synchronize the ghost parts*
2370	*
2371	* \tparam prp... Properties to synchronize
2372	*
2373	*/
2374	template<template<typename,typename> class op,int... prp> void ghost_put()
2375	{
2376	#ifdef SE_CLASS2
2377	check_valid(this,`8`);
2378	#endif
2379
2380	// Convert the ghost internal boxes into grid unit boxes
2381	create_ig_box();
2382
2383	// Convert the ghost external boxes into grid unit boxes
2384	create_eg_box();
2385
2386	// Convert the local ghost internal boxes into grid unit boxes
2387	create_local_ig_box();
2388
2389	// Convert the local external ghost boxes into grid unit boxes
2390	create_local_eg_box();
2391
2392	grid_dist_id_comm<dim,St,T,Decomposition,Memory,device_grid>::template ghost_put_<op,prp...>(dec,
2393	ig_box,
2394	eg_box,
2395	loc_ig_box,
2396	loc_eg_box,
2397	gdb_ext,
2398	loc_grid,
2399	g_id_to_internal_ghost_box);
2400	}
2401
2402
2403	/! \brief Copy the give grid into this grid*
2404	*
2405	* It copy the first grid into the given grid (No ghost)
2406	*
2407	* \warning the Decomposition must be ensured to be the same, otherwise crashes can happen, if you want to copy the grid independently from the decomposition please use the operator equal
2408	*
2409	* \param g Grid to copy
2410	* \param use_memcpy use memcpy function if possible
2411	*
2412	* \return itself
2413	*
2414	*/
2415	grid_dist_id<dim,St,T,Decomposition,Memory,device_grid> & copy(grid_dist_id<dim,St,T,Decomposition,Memory,device_grid> & g, bool use_memcpy = true)
2416	{
2417	if (T::noPointers() == true && use_memcpy)
2418	{
2419	for (size_t i = `0` ; i < this->getN_loc_grid() ; i++)
2420	{
2421	auto & gs_src = this->get_loc_grid(i).getGrid();
2422
2423	long int start = gs_src.LinId(gdb_ext.get(i).Dbox.getKP1());
2424	long int stop = gs_src.LinId(gdb_ext.get(i).Dbox.getKP2());
2425
2426	if (stop < start) {continue;}
2427
2428	void * dst = static_cast<void >(static_cast<char* >(this->get_loc_grid(i).getPointer()) + startsizeof(T));
2429	void * src = static_cast<void >(static_cast<char* >(g.get_loc_grid(i).getPointer()) + startsizeof(T));
2430
2431	memcpy(dst,src,sizeof(T) * (stop + `1` - start));
2432	}
2433	}
2434	else
2435	{
2436	grid_key_dx<dim> cnt[`1`];
2437	cnt[`0`].zero();
2438
2439	for (size_t i = `0` ; i < this->getN_loc_grid() ; i++)
2440	{
2441	auto & dst = this->get_loc_grid(i);
2442	auto & src = g.get_loc_grid(i);
2443
2444	auto it = this->get_loc_grid_iterator_stencil(i,cnt);
2445
2446	while (it.isNext())
2447	{
2448	// center point
2449	auto Cp = it.template getStencil<`0`>();
2450
2451	dst.insert_o(Cp) = src.get_o(Cp);
2452
2453	++it;
2454	}
2455	}
2456	}
2457
2458	return *this;
2459	}
2460
2461	/! \brief Copy the give grid into this grid*
2462	*
2463	* It copy the first grid into the given grid (No ghost)
2464	*
2465	* \warning the Decomposition must be ensured to be the same, otherwise crashes can happen, if you want to copy the grid independently from the decomposition please use the operator equal
2466	*
2467	* \param g Grid to copy
2468	* \param use_memcpy use memcpy function if possible
2469	*
2470	* \return itself
2471	*
2472	*/
2473	grid_dist_id<dim,St,T,Decomposition,Memory,device_grid> & copy_sparse(grid_dist_id<dim,St,T,Decomposition,Memory,device_grid> & g, bool use_memcpy = true)
2474	{
2475	grid_key_dx<dim> cnt[`1`];
2476	cnt[`0`].zero();
2477
2478	for (size_t i = `0` ; i < this->getN_loc_grid() ; i++)
2479	{
2480	auto & dst = this->get_loc_grid(i);
2481	auto & src = g.get_loc_grid(i);
2482
2483	dst = src;
2484	}
2485	return *this;
2486	}
2487
2488	/! \brief Get the spacing on each dimension*
2489	*
2490	* \return the spacing of the grid on each dimension as a point
2491	*
2492	*/
2493	Point<dim,St> getSpacing()
2494	{
2495	return cd_sm.getCellBox().getP2();
2496	}
2497
2498	/! \brief Convert a g_dist_key_dx into a global key*
2499	*
2500	* \see grid_dist_key_dx
2501	* \see grid_dist_iterator
2502	*
2503	* \param k grid_dist_key_dx point (in general returned by the iterators)
2504	*
2505	* \return the global position in the grid
2506	*
2507	*/
2508	inline grid_key_dx<dim> getGKey(const grid_dist_key_dx<dim> & k)
2509	{
2510	#ifdef SE_CLASS2
2511	check_valid(this,`8`);
2512	#endif
2513	// Get the sub-domain id
2514	size_t sub_id = k.getSub();
2515
2516	grid_key_dx<dim> k_glob = k.getKey();
2517
2518	// shift
2519	k_glob = k_glob + gdb_ext.get(sub_id).origin;
2520
2521	return k_glob;
2522	}
2523
2524	/! \brief Add the computation cost on the decomposition using a resolution function*
2525	*
2526	*
2527	* \param md Model to use
2528	* \param ts It is an optional parameter approximately should be the number of ghost get between two
2529	* rebalancing at first decomposition this number can be ignored (default = 1) because not used
2530	*
2531	*/
2532	template <typename Model>inline void addComputationCosts(Model md=Model(), size_t ts = `1`)
2533	{
2534	CellDecomposer_sm<dim, St, shift<dim,St>> cdsm;
2535
2536	Decomposition & dec = getDecomposition();
2537	auto & dist = getDecomposition().getDistribution();
2538
2539	cdsm.setDimensions(dec.getDomain(), dec.getDistGrid().getSize(), `0`);
2540
2541	// Invert the id to positional
2542
2543	Point<dim,St> p;
2544	for (size_t i = `0`; i < dist.getNOwnerSubSubDomains() ; i++)
2545	{
2546	dist.getSubSubDomainPos(i,p);
2547	dec.setSubSubDomainComputationCost(dist.getOwnerSubSubDomain(i) , `1` + md.resolution(p));
2548	}
2549
2550	dec.computeCommunicationAndMigrationCosts(ts);
2551
2552	dist.setDistTol(md.distributionTol());
2553	}
2554
2555	/! \brief apply a convolution using the stencil N*
2556	*
2557	*
2558	*/
2559	template<unsigned int prop_src, unsigned int prop_dst, unsigned int stencil_size, unsigned int N, typename lambda_f, typename ... ArgsT >
2560	void conv(int (& stencil)[N][dim], grid_key_dx<`3`> start, grid_key_dx<`3`> stop , lambda_f func, ArgsT ... args)
2561	{
2562	for (int i = `0` ; i < loc_grid.size() ; i++)
2563	{
2564	Box<dim,long int> inte;
2565
2566	Box<dim,long int> base;
2567	for (int j = `0` ; j < dim ; j++)
2568	{
2569	base.setLow(j,(long int)start.get(j) - (long int)gdb_ext.get(i).origin.get(j));
2570	base.setHigh(j,(long int)stop.get(j) - (long int)gdb_ext.get(i).origin.get(j));
2571	}
2572
2573	Box<dim,long int> dom = gdb_ext.get(i).Dbox;
2574
2575	bool overlap = dom.Intersect(base,inte);
2576
2577	if (overlap == true)
2578	{
2579	loc_grid.get(i).template conv<prop_src,prop_dst,stencil_size>(stencil,inte.getKP1(),inte.getKP2(),func,args...);
2580	}
2581	}
2582	}
2583
2584	/! \brief apply a convolution using the stencil N*
2585	*
2586	*
2587	*/
2588	template<unsigned int prop_src, unsigned int prop_dst, unsigned int stencil_size, typename lambda_f, typename ... ArgsT >
2589	void conv_cross(grid_key_dx<`3`> start, grid_key_dx<`3`> stop , lambda_f func, ArgsT ... args)
2590	{
2591	for (int i = `0` ; i < loc_grid.size() ; i++)
2592	{
2593	Box<dim,long int> inte;
2594
2595	Box<dim,long int> base;
2596	for (int j = `0` ; j < dim ; j++)
2597	{
2598	base.setLow(j,(long int)start.get(j) - (long int)gdb_ext.get(i).origin.get(j));
2599	base.setHigh(j,(long int)stop.get(j) - (long int)gdb_ext.get(i).origin.get(j));
2600	}
2601
2602	Box<dim,long int> dom = gdb_ext.get(i).Dbox;
2603
2604	bool overlap = dom.Intersect(base,inte);
2605
2606	if (overlap == true)
2607	{
2608	loc_grid.get(i).template conv_cross<prop_src,prop_dst,stencil_size>(inte.getKP1(),inte.getKP2(),func,args...);
2609	}
2610	}
2611	}
2612
2613	/! \brief apply a convolution using the stencil N*
2614	*
2615	*
2616	*/
2617	template<unsigned int stencil_size, typename v_type, typename lambda_f, typename ... ArgsT >
2618	void conv_cross_ids(grid_key_dx<`3`> start, grid_key_dx<`3`> stop , lambda_f func, ArgsT ... args)
2619	{
2620	for (int i = `0` ; i < loc_grid.size() ; i++)
2621	{
2622	Box<dim,long int> inte;
2623
2624	Box<dim,long int> base;
2625	for (int j = `0` ; j < dim ; j++)
2626	{
2627	base.setLow(j,(long int)start.get(j) - (long int)gdb_ext.get(i).origin.get(j));
2628	base.setHigh(j,(long int)stop.get(j) - (long int)gdb_ext.get(i).origin.get(j));
2629	}
2630
2631	Box<dim,long int> dom = gdb_ext.get(i).Dbox;
2632
2633	bool overlap = dom.Intersect(base,inte);
2634
2635	if (overlap == true)
2636	{
2637	loc_grid.get(i).template conv_cross_ids<stencil_size,v_type>(inte.getKP1(),inte.getKP2(),func,args...);
2638	}
2639	}
2640	}
2641
2642	/! \brief apply a convolution using the stencil N*
2643	*
2644	*
2645	*/
2646	template<unsigned int prop_src1, unsigned int prop_src2, unsigned int prop_dst1, unsigned int prop_dst2, unsigned int stencil_size, unsigned int N, typename lambda_f, typename ... ArgsT >
2647	void conv2(int (& stencil)[N][dim], grid_key_dx<dim> start, grid_key_dx<dim> stop , lambda_f func, ArgsT ... args)
2648	{
2649	for (int i = `0` ; i < loc_grid.size() ; i++)
2650	{
2651	Box<dim,long int> inte;
2652
2653	Box<dim,long int> base;
2654	for (int j = `0` ; j < dim ; j++)
2655	{
2656	base.setLow(j,(long int)start.get(j) - (long int)gdb_ext.get(i).origin.get(j));
2657	base.setHigh(j,(long int)stop.get(j) - (long int)gdb_ext.get(i).origin.get(j));
2658	}
2659
2660	Box<dim,long int> dom = gdb_ext.get(i).Dbox;
2661
2662	bool overlap = dom.Intersect(base,inte);
2663
2664	if (overlap == true)
2665	{
2666	loc_grid.get(i).template conv2<prop_src1,prop_src2,prop_dst1,prop_dst2,stencil_size>(stencil,inte.getKP1(),inte.getKP2(),func,create_vcluster().rank(),args...);
2667	}
2668	}
2669	}
2670
2671	/! \brief apply a convolution using the stencil N*
2672	*
2673	*
2674	*/
2675	template<unsigned int prop_src1, unsigned int prop_src2, unsigned int prop_dst1, unsigned int prop_dst2, unsigned int stencil_size, typename lambda_f, typename ... ArgsT >
2676	void conv2(grid_key_dx<dim> start, grid_key_dx<dim> stop , lambda_f func, ArgsT ... args)
2677	{
2678	for (int i = `0` ; i < loc_grid.size() ; i++)
2679	{
2680	Box<dim,long int> inte;
2681
2682	Box<dim,long int> base;
2683	for (int j = `0` ; j < dim ; j++)
2684	{
2685	base.setLow(j,(long int)start.get(j) - (long int)gdb_ext.get(i).origin.get(j));
2686	base.setHigh(j,(long int)stop.get(j) - (long int)gdb_ext.get(i).origin.get(j));
2687	}
2688
2689	Box<dim,long int> dom = gdb_ext.get(i).Dbox;
2690
2691	bool overlap = dom.Intersect(base,inte);
2692
2693	if (overlap == true)
2694	{
2695	loc_grid.get(i).template conv2<prop_src1,prop_src2,prop_dst1,prop_dst2,stencil_size>(inte.getKP1(),inte.getKP2(),func,args...);
2696	}
2697	}
2698	}
2699
2700	template<typename NNtype>
2701	void findNeighbours()
2702	{
2703	for (int i = `0` ; i < loc_grid.size() ; i++)
2704	{
2705	loc_grid.get(i).findNeighbours();
2706	}
2707	}
2708
2709	/! \brief apply a convolution using the stencil N*
2710	*
2711	*
2712	*/
2713	template<unsigned int prop_src1, unsigned int prop_src2, unsigned int prop_dst1, unsigned int prop_dst2, unsigned int stencil_size, typename lambda_f, typename ... ArgsT >
2714	void conv_cross2(grid_key_dx<`3`> start, grid_key_dx<`3`> stop , lambda_f func, ArgsT ... args)
2715	{
2716	for (int i = `0` ; i < loc_grid.size() ; i++)
2717	{
2718	Box<dim,long int> inte;
2719
2720	Box<dim,long int> base;
2721	for (int j = `0` ; j < dim ; j++)
2722	{
2723	base.setLow(j,(long int)start.get(j) - (long int)gdb_ext.get(i).origin.get(j));
2724	base.setHigh(j,(long int)stop.get(j) - (long int)gdb_ext.get(i).origin.get(j));
2725	}
2726
2727	Box<dim,long int> dom = gdb_ext.get(i).Dbox;
2728
2729	bool overlap = dom.Intersect(base,inte);
2730
2731	if (overlap == true)
2732	{
2733	loc_grid.get(i).template conv_cross2<prop_src1,prop_src2,prop_dst1,prop_dst2,stencil_size>(inte.getKP1(),inte.getKP2(),func,args...);
2734	}
2735	}
2736	}
2737
2738	/! \brief Write the distributed grid information*
2739	*
2740	* * grid_X.vtk Output each local grids for each local processor X
2741	* * internal_ghost_X.vtk Internal ghost boxes in grid units for the local processor X
2742	*
2743	* \param output directory where to put the files + prefix
2744	* \param opt options
2745	*
2746	* \return true if the write operation succeed
2747	*
2748	*/
2749	bool write(std::string output, size_t opt = VTK_WRITER \| FORMAT_BINARY )
2750	{
2751	#ifdef SE_CLASS2
2752	check_valid(this,`8`);
2753	#endif
2754
2755	file_type ft = file_type::ASCII;
2756
2757	if (opt & FORMAT_BINARY)
2758	ft = file_type::BINARY;
2759
2760	// Create a writer and write
2761	VTKWriter<boost::mpl::pair<device_grid,float>,VECTOR_GRIDS> vtk_g;
2762	for (size_t i = `0` ; i < loc_grid.size() ; i++)
2763	{
2764	Point<dim,St> offset = getOffset(i);
2765
2766	if (opt & PRINT_GHOST)
2767	{vtk_g.add(loc_grid.get(i),offset,cd_sm.getCellBox().getP2(),gdb_ext.get(i).GDbox);}
2768	else
2769	{vtk_g.add(loc_grid.get(i),offset,cd_sm.getCellBox().getP2(),gdb_ext.get(i).Dbox);}
2770	}
2771	vtk_g.write(output + "_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk", prp_names, "grids", ft);
2772
2773	return true;
2774	}
2775
2776	/! \brief Write all grids indigually*
2777	*
2778	* \param output files
2779	*
2780	*/
2781	bool write_debug(std::string output)
2782	{
2783	for (int i = `0` ; i < getN_loc_grid() ; i++)
2784	{
2785	Point<dim,St> sp;
2786	Point<dim,St> offset;
2787
2788	for (int j = `0` ; j < dim ; j++)
2789	{sp.get(j) = this->spacing(j);}
2790
2791	offset = gdb_ext.get(i).origin;
2792
2793	get_loc_grid(i).write_debug(output + "_" + std::to_string(i) + "_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk",sp,offset);
2794	}
2795
2796	return true;
2797	}
2798
2799	/! \brief Write the distributed grid information*
2800	*
2801	* * grid_X.vtk Output each local grids for each local processor X
2802	* * internal_ghost_X.vtk Internal ghost boxes in grid units for the local processor X
2803	*
2804	* \param output directory where to put the files + prefix
2805	* \param i frame number
2806	* \param opt options
2807	*
2808	* \return true id the write succeed
2809	*
2810	*/
2811	bool write_frame(std::string output, size_t i, size_t opt = VTK_WRITER \| FORMAT_ASCII)
2812	{
2813	#ifdef SE_CLASS2
2814	check_valid(this,`8`);
2815	#endif
2816	file_type ft = file_type::ASCII;
2817
2818	if (opt & FORMAT_BINARY)
2819	ft = file_type::BINARY;
2820
2821	// Create a writer and write
2822	VTKWriter<boost::mpl::pair<device_grid,float>,VECTOR_GRIDS> vtk_g;
2823	for (size_t i = `0` ; i < loc_grid.size() ; i++)
2824	{
2825	Point<dim,St> offset = getOffset(i);
2826	vtk_g.add(loc_grid.get(i),offset,cd_sm.getCellBox().getP2(),gdb_ext.get(i).Dbox);
2827	}
2828	vtk_g.write(output + "_" + std::to_string(v_cl.getProcessUnitID()) + "_" + std::to_string(i) + ".vtk",prp_names,"grids",ft);
2829
2830	return true;
2831	}
2832
2833
2834
2835	/! \brief Get the i sub-domain grid*
2836	*
2837	* \param i sub-domain
2838	*
2839	* \return local grid
2840	*
2841	*/
2842	device_grid & get_loc_grid(size_t i)
2843	{
2844	return loc_grid.get(i);
2845	}
2846
2847	/! \brief Get the i sub-domain grid*
2848	*
2849	* \param i sub-domain
2850	*
2851	* \return local grid
2852	*
2853	*/
2854	grid_key_dx_iterator_sub<dim,no_stencil> get_loc_grid_iterator(size_t i)
2855	{
2856	return grid_key_dx_iterator_sub<dim,no_stencil>(loc_grid.get(i).getGrid(),
2857	gdb_ext.get(i).Dbox.getKP1(),
2858	gdb_ext.get(i).Dbox.getKP2());
2859	}
2860
2861	/! \brief Get the i sub-domain grid*
2862	*
2863	* \param i sub-domain
2864	*
2865	* \return local grid
2866	*
2867	*/
2868	template<unsigned int Np>
2869	grid_key_dx_iterator_sub<dim,stencil_offset_compute<dim,Np>,typename device_grid::linearizer_type>
2870	get_loc_grid_iterator_stencil(size_t i,const grid_key_dx<dim> (& stencil_pnt)[Np])
2871	{
2872	return grid_key_dx_iterator_sub<dim,stencil_offset_compute<dim,Np>,typename device_grid::linearizer_type>(loc_grid.get(i).getGrid(),
2873	gdb_ext.get(i).Dbox.getKP1(),
2874	gdb_ext.get(i).Dbox.getKP2(),
2875	stencil_pnt);
2876	}
2877
2878	/! \brief Return the number of local grid*
2879	*
2880	* \return the number of local grid
2881	*
2882	*/
2883	size_t getN_loc_grid()
2884	{
2885	return loc_grid.size();
2886	}
2887
2888
2889	/! \brief It return the id of structure in the allocation list*
2890	*
2891	* \see print_alloc and SE_CLASS2
2892	*
2893	* \return the id
2894	*
2895	*/
2896	long int who()
2897	{
2898	#ifdef SE_CLASS2
2899	return check_whoami(this,`8`);
2900	#else
2901	return -`1`;
2902	#endif
2903	}
2904
2905	/! \brief It print the internal ghost boxes and external ghost boxes in global unit*
2906	*
2907	*
2908	*/
2909	void debugPrint()
2910	{
2911	size_t tot_volume = `0`;
2912
2913	std::cout << "-------- External Ghost boxes ---------- " << std::endl;
2914
2915	for (size_t i = `0` ; i < eg_box.size() ; i++)
2916	{
2917	std::cout << "Processor: " << eg_box.get(i).prc << " Boxes:" << std::endl;
2918
2919	for (size_t j = `0`; j < eg_box.get(i).bid.size() ; j++)
2920	{
2921	std::cout << " Box: " << eg_box.get(i).bid.get(j).g_e_box.toString() << " Id: " << eg_box.get(i).bid.get(j).g_id << std::endl;
2922	tot_volume += eg_box.get(i).bid.get(j).g_e_box.getVolumeKey();
2923	}
2924	}
2925
2926	std::cout << "TOT volume external ghost " << tot_volume << std::endl;
2927
2928	std::cout << "-------- Internal Ghost boxes ---------- " << std::endl;
2929
2930	tot_volume = `0`;
2931	for (size_t i = `0` ; i < ig_box.size() ; i++)
2932	{
2933	std::cout << "Processor: " << ig_box.get(i).prc << " Boxes:" << std::endl;
2934
2935	for (size_t j = `0` ; j < ig_box.get(i).bid.size() ; j++)
2936	{
2937	std::cout << " Box: " << ig_box.get(i).bid.get(j).box.toString() << " Id: " << ig_box.get(i).bid.get(j).g_id << std::endl;
2938	tot_volume += ig_box.get(i).bid.get(j).box.getVolumeKey();
2939	}
2940	}
2941
2942	std::cout << "TOT volume internal ghost " << tot_volume << std::endl;
2943	}
2944
2945	/! \brief Set the properties names*
2946	*
2947	* It is useful to specify name for the properties in vtk writers
2948	*
2949	* \param names set of properties names
2950	*
2951	*/
2952	void setPropNames(const openfpm::vector<std::string> & names)
2953	{
2954	prp_names = names;
2955	}
2956
2957	/! \brief It delete all the points*
2958	*
2959	* This function on dense does nothing in case of dense grid but in case of
2960	* sparse_grid it kills all the points
2961	*
2962	*/
2963	void clear()
2964	{
2965	for (size_t i = `0` ; i < loc_grid.size() ; i++)
2966	{loc_grid.get(i).clear();}
2967	}
2968
2969	/! \brief construct link between levels*
2970	*
2971	* \praram grid_up grid level up
2972	* \param grid_dw grid level down
2973	*
2974	*/
2975	void construct_link(self & grid_up, self & grid_dw)
2976	{
2977	for (int i = `0` ; i < loc_grid.size() ; i++)
2978	{
2979	loc_grid.get(i).construct_link(grid_up.get_loc_grid(i),grid_dw.get_loc_grid(i),v_cl.getmgpuContext());
2980	}
2981	}
2982
2983	/! \brief construct link between current and the level down*
2984	*
2985	*
2986	* \param grid_dw grid level down
2987	*
2988	*/
2989	void construct_link_dw(self & grid_dw, openfpm::vector<offset_mv<dim>> & mvof)
2990	{
2991	for (int i = `0` ; i < loc_grid.size() ; i++)
2992	{
2993	Point<dim,int> p_dw;
2994	for(int j = `0` ; j < dim ; j++)
2995	{p_dw.get(j) = mvof.get(i).dw.get(j);}
2996
2997	loc_grid.get(i).construct_link_dw(grid_dw.get_loc_grid(i),gdb_ext.get(i).Dbox,p_dw,v_cl.getmgpuContext());
2998	}
2999	}
3000
3001	/! \brief construct link between current and the level up*
3002	*
3003	*
3004	* \param grid_dw grid level down
3005	*
3006	*/
3007	void construct_link_up(self & grid_up, openfpm::vector<offset_mv<dim>> & mvof)
3008	{
3009	for (int i = `0` ; i < loc_grid.size() ; i++)
3010	{
3011	Point<dim,int> p_up;
3012	for(int j = `0` ; j < dim ; j++)
3013	{p_up.get(j) = mvof.get(i).up.get(j);}
3014
3015	loc_grid.get(i).construct_link_up(grid_up.get_loc_grid(i),gdb_ext.get(i).Dbox,p_up,v_cl.getmgpuContext());
3016	}
3017	}
3018
3019	/! \brief construct link between current and the level up*
3020	*
3021	*
3022	* \param grid_dw grid level down
3023	*
3024	*/
3025	template<typename stencil_type>
3026	void tagBoundaries()
3027	{
3028	for (int i = `0` ; i < loc_grid.size() ; i++)
3029	{
3030	// we limit to the domain subset for tagging
3031
3032	Box_check<dim,unsigned int> chk(gdb_ext.get(i).Dbox);
3033
3034
3035	loc_grid.get(i).template tagBoundaries<stencil_type>(v_cl.getmgpuContext(),chk);
3036	}
3037	}
3038
3039	/! \brief It move all the grid parts that do not belong to the local processor to the respective processor*
3040	*
3041	*/
3042	void map(size_t opt = `0`)
3043	{
3044	// Save the background values
3045	T bv;
3046
3047	copy_aggregate_dual<decltype(loc_grid.get(`0`).getBackgroundValue()),
3048	T> ca(loc_grid.get(`0`).getBackgroundValue(),bv);
3049
3050	boost::mpl::for_each_ref<boost::mpl::range_c<int,`0`,T::max_prop>>(ca);
3051
3052	if (!(opt & NO_GDB_EXT_SWITCH))
3053	{
3054	gdb_ext_old = gdb_ext;
3055	loc_grid_old = loc_grid;
3056
3057	InitializeStructures(g_sz,bx_def,gint,bx_def.size() != `0`);
3058	}
3059
3060	getGlobalGridsInfo(gdb_ext_global);
3061
3062	this->map_(dec,cd_sm,loc_grid,loc_grid_old,gdb_ext,gdb_ext_old,gdb_ext_global);
3063
3064	loc_grid_old.clear();
3065	loc_grid_old.shrink_to_fit();
3066	gdb_ext_old.clear();
3067
3068	// reset ghost structure to recalculate
3069	reset_ghost_structures();
3070
3071	// Reset the background values
3072	setBackgroundValue(bv);
3073	}
3074
3075	/! \brief Save the grid state on HDF5*
3076	*
3077	* \param filename output filename
3078	*
3079	*/
3080	inline void save(const std::string & filename) const
3081	{
3082	HDF5_writer<GRID_DIST> h5s;
3083
3084	h5s.save(filename,loc_grid,gdb_ext);
3085	}
3086
3087	/! \brief Reload the grid from HDF5 file*
3088	*
3089	* \param filename output filename
3090	*
3091	*/
3092	inline void load(const std::string & filename)
3093	{
3094	HDF5_reader<GRID_DIST> h5l;
3095
3096	h5l.load<device_grid>(filename,loc_grid_old,gdb_ext_old);
3097
3098	if (v_cl.size() != `1`)
3099	{
3100	// Map the distributed grid
3101	map(NO_GDB_EXT_SWITCH);
3102	}
3103	else
3104	{
3105	for (int i = `0` ; i < gdb_ext_old.size() ; i++)
3106	{
3107	auto & lg = loc_grid_old.get(i);
3108	auto it_src = lg.getIterator(gdb_ext_old.get(i).Dbox.getKP1(),gdb_ext_old.get(i).Dbox.getKP2());
3109	auto & dg = loc_grid.get(`0`);
3110
3111	grid_key_dx<dim> orig;
3112	for (int j = `0` ; j < dim ; j++)
3113	{
3114	orig.set_d(j,gdb_ext_old.get(i).origin.get(j) + gdb_ext_old.get(i).Dbox.getKP1().get(j));
3115	}
3116
3117	while (it_src.isNext())
3118	{
3119	auto key = it_src.get();
3120	auto key_dst = key + orig;
3121
3122	dg.get_o(key_dst) = lg.get_o(key);
3123
3124	++it_src;
3125	}
3126	}
3127	}
3128	}
3129
3130	/! \brief This is a meta-function return which type of sub iterator a grid produce*
3131	*
3132	* \return the type of the sub-grid iterator
3133	*
3134	*/
3135	template <typename stencil = no_stencil>
3136	static grid_dist_iterator_sub<dim,device_grid> type_of_subiterator()
3137	{
3138	return grid_key_dx_iterator_sub<dim, stencil>();
3139	}
3140
3141	/! \brief Get the internal local ghost box*
3142	*
3143	* \return the internal local ghost box
3144	*
3145	*/
3146	const openfpm::vector<i_lbox_grid<dim>> & get_loc_ig_box()
3147	{
3148	return this->loc_ig_box;
3149	}
3150
3151	/! \brief Get the internal ghost box*
3152	*
3153	* \return the internal local ghost box
3154	*
3155	*/
3156	const openfpm::vector<i_lbox_grid<dim>> & get_ig_box()
3157	{
3158	return this->ig_box;
3159	}
3160
3161	void print_stats()
3162	{
3163	std::cout << "-- REPORT --" << std::endl;
3164	#ifdef ENABLE_GRID_DIST_ID_PERF_STATS
3165	std::cout << "Processor: " << v_cl.rank() << " Time spent in packing data: " << tot_pack << std::endl;
3166	std::cout << "Processor: " << v_cl.rank() << " Time spent in sending and receving data: " << tot_sendrecv << std::endl;
3167	std::cout << "Processor: " << v_cl.rank() << " Time spent in merging: " << tot_merge << std::endl;
3168	std::cout << "Processor: " << v_cl.rank() << " Time spent in local merging: " << tot_loc_merge << std::endl;
3169	#else
3170
3171	std::cout << "Enable ENABLE_GRID_DIST_ID_PERF_STATS if you want to activate this feature" << std::endl;
3172
3173	#endif
3174	}
3175
3176	void clear_stats()
3177	{
3178	#ifdef ENABLE_GRID_DIST_ID_PERF_STATS
3179	tot_pack = `0`;
3180	tot_sendrecv = `0`;
3181	tot_merge = `0`;
3182	#else
3183
3184	std::cout << "Enable ENABLE_GRID_DIST_ID_PERF_STATS if you want to activate this feature" << std::endl;
3185
3186	#endif
3187	}
3188
3189	#ifdef __NVCC__
3190
3191	/! \brief Set the number inserts each GPU thread do*
3192	*
3193	* \param n_ins number of insert per thread
3194	*
3195	*/
3196	void setNumberOfInsertPerThread(size_t n_ins)
3197	{
3198	gpu_n_insert_thread = n_ins;
3199	}
3200
3201	template<typename func_t,typename it_t, typename ... args_t>
3202	void iterateGridGPU(it_t & it, args_t ... args)
3203	{
3204	// setGPUInsertBuffer must be called in anycase even with 0 points to insert
3205	// the loop "it.isNextGrid()" does not guarantee to call it for all local grids
3206	for (size_t i = `0` ; i < loc_grid.size() ; i++)
3207	{loc_grid.get(i).setGPUInsertBuffer(`0ul`,`1ul`);}
3208
3209	while(it.isNextGrid())
3210	{
3211	Box<dim,size_t> b = it.getGridBox();
3212
3213	size_t i = it.getGridId();
3214
3215	auto ite = loc_grid.get(i).getGridGPUIterator(b.getKP1int(),b.getKP2int());
3216
3217	loc_grid.get(i).setGPUInsertBuffer(ite.nblocks(),ite.nthrs());
3218	loc_grid.get(i).initializeGPUInsertBuffer();
3219
3220	ite_gpu_dist<dim> itd = ite;
3221
3222	for (int j = `0` ; j < dim ; j++)
3223	{
3224	itd.origin.set_d(j,gdb_ext.get(i).origin.get(j));
3225	itd.start_base.set_d(j,`0`);
3226	}
3227
3228	CUDA_LAUNCH((grid_apply_functor),ite,loc_grid.get(i).toKernel(),itd,func_t(),args...);
3229
3230	it.nextGrid();
3231	}
3232	}
3233
3234	/! \brief Remove the points*
3235	*
3236	* \param box Remove all the points in the box
3237	*
3238	*/
3239	void removePoints(Box<dim,size_t> & box)
3240	{
3241	Box<dim,long int> box_i = box;
3242
3243	for (size_t i = `0` ; i < loc_grid.size() ; i++)
3244	{
3245	Box<dim,long int> bx = gdb_ext.get(i).Dbox + gdb_ext.get(i).origin;
3246
3247	Box<dim,long int> bout;
3248	bool inte = bx.Intersect(box,bout);
3249	bout -= gdb_ext.get(i).origin;
3250
3251	if (inte == true)
3252	{
3253	loc_grid.get(i).copyRemoveReset();
3254	loc_grid.get(i).remove(bout);
3255	loc_grid.get(i).removePoints(v_cl.getmgpuContext());
3256	}
3257	}
3258	}
3259
3260	/! \brief Move the memory from the device to host memory*
3261	*
3262	*/
3263	template<unsigned int ... prp> void deviceToHost()
3264	{
3265	for (size_t i = `0` ; i < loc_grid.size() ; i++)
3266	{
3267	loc_grid.get(i).template deviceToHost<prp ...>();
3268	}
3269	}
3270
3271	/! \brief Move the memory from the device to host memory*
3272	*
3273	*/
3274	template<unsigned int ... prp> void hostToDevice()
3275	{
3276	for (size_t i = `0` ; i < loc_grid.size() ; i++)
3277	{
3278	loc_grid.get(i).template hostToDevice<prp ...>();
3279	}
3280	}
3281
3282	#endif
3283
3284
3285	//! Define friend classes
3286	//\cond
3287	friend grid_dist_id<dim,St,T,typename Decomposition::extended_type,Memory,device_grid>;
3288	//\endcond
3289	};
3290
3291
3292	template<unsigned int dim, typename St, typename T, typename Memory = HeapMemory, typename Decomposition = CartDecomposition<dim,St> >
3293	using sgrid_dist_id = grid_dist_id<dim,St,T,Decomposition,Memory,sgrid_cpu<dim,T,Memory>>;
3294
3295	template<unsigned int dim, typename St, typename T, typename Memory = HeapMemory, typename Decomposition = CartDecomposition<dim,St>>
3296	using sgrid_dist_soa = grid_dist_id<dim,St,T,Decomposition,Memory,sgrid_soa<dim,T,Memory>>;
3297
3298	template<unsigned int dim, typename St, typename T, typename devg, typename Memory = HeapMemory, typename Decomposition = CartDecomposition<dim,St>>
3299	using grid_dist_id_devg = grid_dist_id<dim,St,T,Decomposition,Memory,devg>;
3300
3301	#ifdef __NVCC__
3302	template<unsigned int dim, typename St, typename T, typename Memory = CudaMemory, typename Decomposition = CartDecomposition<dim,St,CudaMemory,memory_traits_inte> >
3303	using sgrid_dist_id_gpu = grid_dist_id<dim,St,T,Decomposition,Memory,SparseGridGpu<dim,T>>;
3304
3305	template<unsigned int dim, typename St, typename T, typename Memory = CudaMemory, typename Decomposition = CartDecomposition<dim,St,CudaMemory,memory_traits_inte> >
3306	using sgrid_dist_sid_gpu = grid_dist_id<dim,St,T,Decomposition,Memory,SparseGridGpu<dim,T,default_edge<dim>::type::value,default_edge<dim>::tb::value,int>>;
3307	#endif
3308
3309	#endif
3310

Browse the source code of openfpm/src/Grid/grid_dist_id.hpp