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

#include <vector>
5
#include <unordered_map>
incardon's avatar
incardon committed
6
#include "Grid/map_grid.hpp"
incardon's avatar
incardon committed
7
#include "VCluster/VCluster.hpp"
incardon's avatar
incardon committed
8
#include "Space/SpaceBox.hpp"
9
#include "util/mathutil.hpp"
incardon's avatar
incardon committed
10 11 12
#include "Iterators/grid_dist_id_iterator_dec.hpp"
#include "Iterators/grid_dist_id_iterator.hpp"
#include "Iterators/grid_dist_id_iterator_sub.hpp"
incardon's avatar
incardon committed
13
#include "grid_dist_key.hpp"
14
#include "NN/CellList/CellDecomposer.hpp"
15 16
#include "util/object_util.hpp"
#include "memory/ExtPreAlloc.hpp"
17
#include "VTKWriter/VTKWriter.hpp"
18 19
#include "Packer_Unpacker/Packer.hpp"
#include "Packer_Unpacker/Unpacker.hpp"
20
#include "Decomposition/CartDecomposition.hpp"
21
#include "data_type/aggregate.hpp"
22
#include "hdf5.h"
Yaroslav's avatar
Yaroslav committed
23
#include "grid_dist_id_comm.hpp"
incardon's avatar
incardon committed
24
#include "HDF5_wr/HDF5_wr.hpp"
incardon's avatar
incardon committed
25
#include "SparseGrid/SparseGrid.hpp"
26 27 28 29
#ifdef __NVCC__
#include "SparseGridGpu/SparseGridGpu.hpp"
#include "cuda/grid_dist_id_kernels.cuh"
#endif
incardon's avatar
incardon committed
30

incardon's avatar
incardon committed
31

32
//! Internal ghost box sent to construct external ghost box into the other processors
Pietro Incardona's avatar
Pietro Incardona committed
33 34 35
template<unsigned int dim>
struct Box_fix
{
36
	//! Box in global unit
Pietro Incardona's avatar
Pietro Incardona committed
37
	Box<dim,size_t> bx;
38
	//! In which sector live the box
39
	comb<dim> cmb;
40
	//! Global id of the internal ghost box
Pietro Incardona's avatar
Pietro Incardona committed
41
	size_t g_id;
42
	//! from which sub-domain this internal ghost box is generated (or with which sub-domain is overlapping)
43
	size_t r_sub;
Pietro Incardona's avatar
Pietro Incardona committed
44 45
};

incardon's avatar
incardon committed
46
#define NO_GDB_EXT_SWITCH 0x1000
incardon's avatar
incardon committed
47 48 49

/*! \brief This is a distributed grid
 *
50 51
 * Implementation of a distributed grid the decomposition is geometrical, grid
 * is splitted across several processor
incardon's avatar
incardon committed
52 53
 *
 * \param dim Dimensionality of the grid
54 55
 * \param St Type of space where the grid is living
 * \param T object the grid is storing
incardon's avatar
incardon committed
56 57 58 59
 * \param Decomposition Class that decompose the grid for example CartDecomposition
 * \param Mem Is the allocator
 * \param device type of base structure is going to store the data
 *
60
 * ### Create a distributed grid and access it
Pietro Incardona's avatar
Pietro Incardona committed
61
 * \snippet grid_dist_id_unit_test.cpp Create and access a distributed grid
62
 * ### Synchronize the ghosts and check the information
Pietro Incardona's avatar
Pietro Incardona committed
63
 * \snippet grid_dist_id_unit_test.cpp Synchronize the ghost and check the information
64
 * ### Create and access a distributed grid for complex structures
Pietro Incardona's avatar
Pietro Incardona committed
65
 * \snippet grid_dist_id_unit_test.cpp Create and access a distributed grid complex
incardon's avatar
incardon committed
66
 * ### Synchronize a distributed grid for complex structures
Pietro Incardona's avatar
Pietro Incardona committed
67
 * \snippet grid_dist_id_unit_test.cpp Synchronized distributed grid complex
incardon's avatar
incardon committed
68
 * ### Usage of a grid dist iterator sub
incardon's avatar
incardon committed
69
 * \snippet grid_dist_id_iterators_unit_tests.hpp Usage of a sub_grid iterator
incardon's avatar
incardon committed
70
 * ### Construct two grid with the same decomposition
Pietro Incardona's avatar
Pietro Incardona committed
71
 * \snippet grid_dist_id_unit_test.cpp Construct two grid with the same decomposition
72
 *
incardon's avatar
incardon committed
73
 */
74 75 76 77 78 79
template<unsigned int dim,
		 typename St,
		 typename T,
		 typename Decomposition = CartDecomposition<dim,St>,
		 typename Memory=HeapMemory ,
		 typename device_grid=grid_cpu<dim,T> >
Yaroslav's avatar
Yaroslav committed
80
class grid_dist_id : public grid_dist_id_comm<dim,St,T,Decomposition,Memory,device_grid>
incardon's avatar
incardon committed
81
{
82
	//! Domain
83 84
	Box<dim,St> domain;

85
	//! Ghost expansion
86
	Ghost<dim,St> ghost;
incardon's avatar
incardon committed
87

incardon's avatar
incardon committed
88 89 90
	//! Ghost expansion
	Ghost<dim,long int> ghost_int;

incardon's avatar
incardon committed
91
	//! Local grids
92
	mutable openfpm::vector<device_grid> loc_grid;
incardon's avatar
incardon committed
93

94 95
	//! Old local grids
	mutable openfpm::vector<device_grid> loc_grid_old;
incardon's avatar
incardon committed
96 97

	//! Space Decomposition
98
	Decomposition dec;
incardon's avatar
incardon committed
99

incardon's avatar
incardon committed
100 101 102 103 104 105 106 107 108 109 110
	//! gdb_ext markers
	//! In the case where the grid is defined everywhere
	//! gdb_ext_marker is useless and so is empty
	//! in the case we have a grid defined on a smaller set
	//! of boxes gbd_ext_markers indicate the division across
	//! subdomains. For example Sub-domain 0 produce 2 grid
	//! Sub-domain 1 produce 3 grid Sub-domain 2 produce 2 grid
	//! Sub-domain 3 produce 1 grid
	//! gdb_ext_markers contain 0,2,5,7,8
	openfpm::vector<size_t> gdb_ext_markers;

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

114 115 116 117 118
	//! Global gdb_ext
	mutable openfpm::vector<GBoxes<device_grid::dims>> gdb_ext_global;

	//! Extension of each old grid (old): Domain and ghost + domain
	openfpm::vector<GBoxes<device_grid::dims>> gdb_ext_old;
Yaroslav's avatar
Yaroslav committed
119

incardon's avatar
incardon committed
120 121 122
	//! Size of the grid on each dimension
	size_t g_sz[dim];

123
	//! Structure that divide the space into cells
Pietro Incardona's avatar
Pietro Incardona committed
124
	CellDecomposer_sm<dim,St,shift<dim,St>> cd_sm;
incardon's avatar
incardon committed
125

126 127 128
	//! size of the insert pool on gpu
	size_t gpu_insert_pool_size;

129
	//! Communicator class
130
	Vcluster<> & v_cl;
incardon's avatar
incardon committed
131

132 133 134
	//! properties names
	openfpm::vector<std::string> prp_names;

135
	//! It map a global ghost id (g_id) to the external ghost box information
Pietro Incardona's avatar
Pietro Incardona committed
136
	//! It is unique across all the near processor
137 138
	std::unordered_map<size_t,size_t> g_id_to_external_ghost_box;

incardon's avatar
incardon committed
139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208
	/*! Link a received external ghost box to the linked eg_box.
	 * When the grid is defined everywhere for each received external ghost box
	 * exists one eg_box linked to it that contain the information
	 * on how to transfer the information to the associated sub-domain grid.
	 * Unfortunately when we specify where the grid is defined, a received external
	 * ghost box can be linked to multiple sub-domain grids (one sub-domain can have multiple
	 * sub grids).
	 * So in standard situation (grid defined everywhere)
	 * a received external ghost box is linked to a single
	 * eg_box entry and eb_gid_list play mainly no role.
	 * (play no role but must be filled ghost_get expect this structure to be
	 * filled consistently, it will be clear later how to do it in this case).
	 * When the grid is not defined everywhere a received ghost box can be linked
	 * to multiple external ghost boxes. (Like in figure)
	 *
	 * \verbatim
  +--------------------------------------------+------
  |        Sub-domain                          | Another sub-domain
  |          +------+           +---------+    |
  |          |      |           |         |    |
  |          |      |           |         |    |
  |          |      |           |         |    |
  |          |  3   |           |    4    |    |
  |   empty  |      |    empty  |         |    |
  |          |      |           |         |    |
  |          |      |           |         |    |
  |          |      |           |         |    |                         1
  |          |      |           |         |    |
+-+-----+----+------+-----------+---------+----+-----+-----   Processor bound
        |***##########*********#############****|****|
        |                                            |                   0
        |                                            |
        |                                            |
        |                  9                         |
        |                                            |
        |                                            |
        |                                            |
        +--------------------------------------------+

	 * \endverbatim
	 *
	 * As we can see here the grid number 9 on processo 0 has an internal ghost box
	 * The internal ghost-box is sent to processor 1 and is a received external
	 * ghost box. This external ghost box is partially shared in two separated grids.
	 * It is important to note that 3 and 4 are grid defined externally and are
	 * not defined by the sub-domain border. It is important also to note that the sub-domain
	 * granularity in processor 1 define the granularity of the internal ghost box in
	 * processor 0 and consequently every external ghost box in processor 1 is linked
	 *  uniquely with one internal ghost box in processor 0. On the other hand if we have
	 *  a secondary granularity define by external boxes like 3 and 4 this is not anymore
	 *  true and one internal ghost box in 0 can be linked with multiple grids.
	 *  The granularity of the space division is different from the granularity of where
	 *  the grid is defined. Space decomposition exist independently from the data-structure
	 *  and can be shared across multiple data-structure this mean that cannot be redefined
	 *  based on where is the grid definitions.
	 * The internal ghost box could be redefined in order to respect the granularity.
	 * We do not do this for 3 main reason.
	 *
	 * 1) The definition box must be communicated across processors.
	 * 2) An interprocessor global-id link must be established with lower sub-domain
	 *    granularty
	 * 3) Despite the points * are not linked, but must be anyway sent
	 *    to processor 1, this mean that make not too much sense to increase
	 *    the granularity in advance on processor 0, but it is better receive
	 *    the information an than solve the lower granularity locally
	 *    on processor 1
	 *
	 */
	openfpm::vector<e_box_multi<dim>> eb_gid_list;

Pietro Incardona's avatar
Pietro Incardona committed
209 210 211 212
	//! It map a global ghost id (g_id) to the internal ghost box information
	//! (is unique for processor), it is not unique across all the near processor
	openfpm::vector<std::unordered_map<size_t,size_t>> g_id_to_internal_ghost_box;

213
	//! Receiving size
214 215
	openfpm::vector<size_t> recv_sz;

216
	//! Receiving buffer for particles ghost get
217 218
	openfpm::vector<HeapMemory> recv_mem_gg;

219
	//! Grid informations object
incardon's avatar
incardon committed
220 221
	grid_sm<dim,T> ginfo;

222
	//! Grid informations object without type
incardon's avatar
incardon committed
223 224
	grid_sm<dim,void> ginfo_v;

incardon's avatar
incardon committed
225 226 227 228 229 230
	//! Set of boxes that define where the grid is defined
	openfpm::vector<Box<dim,long int>> bx_def;

	//! Indicate if we have to use bx_def to define the grid
	bool use_bx_def = false;

231 232 233 234 235 236
	//! Indicate if the local internal ghost box has been initialized
	bool init_local_i_g_box = false;

	//! Indicate if the local external ghost box has been initialized
	bool init_local_e_g_box = false;

incardon's avatar
incardon committed
237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260
	//! Flag that indicate if the external ghost box has been initialized
	bool init_e_g_box = false;

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

	//! Flag that indicate if the internal and external ghost box has been fixed
	bool init_fix_ie_g_box = false;

	//! Internal ghost boxes in grid units
	openfpm::vector<ip_box_grid<dim>> ig_box;

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

	//! Local internal ghost boxes in grid units
	openfpm::vector<i_lbox_grid<dim>> loc_ig_box;

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

	//! Number of sub-sub-domain for each processor
	size_t v_sub_unit_factor = 64;

261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288
	/*! \brief Call-back to allocate buffer to receive incoming objects (external ghost boxes)
	 *
	 * \param msg_i message size required to receive from i
	 * \param total_msg message size to receive from all the processors
	 * \param total_p the total number of processor want to communicate with you
	 * \param i processor id
	 * \param ri request id (it is an id that goes from 0 to total_p, and is unique
	 *           every time message_alloc is called)
	 * \param ptr void pointer parameter for additional data to pass to the call-back
	 *
	 */
	static void * msg_alloc_external_box(size_t msg_i ,size_t total_msg, size_t total_p, size_t i, size_t ri, void * ptr)
	{
		grid_dist_id<dim,St,T,Decomposition,Memory,device_grid> * g = static_cast<grid_dist_id<dim,St,T,Decomposition,Memory,device_grid> *>(ptr);

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

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

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

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

incardon's avatar
incardon committed
289 290 291 292 293 294 295
	/*! \brief this function is for optimization of the ghost size
	 *
	 * Because the decomposition work in continuum and discrete ghost is
	 *  converted in continuum, in some case continuum ghost because of
	 *  rounding-off error can produce ghost bigger than the discrete selected
	 *   one. This function adjust for this round-off error
	 *
incardon's avatar
incardon committed
296
	 * \param sub_id sub-domain id
incardon's avatar
incardon committed
297 298 299 300
	 * \param sub_domain_other the other sub-domain
	 * \param ib internal ghost box to adjust
	 *
	 */
incardon's avatar
incardon committed
301
	void set_for_adjustment(size_t sub_id,
incardon's avatar
incardon committed
302 303 304 305 306
							const Box<dim,St> & sub_domain_other,
							const comb<dim> & cmb,
							Box<dim,long int> & ib,
							Ghost<dim,long int> & g)
	{
307
		if (g.isInvalidGhost() == true || use_bx_def == true)
incardon's avatar
incardon committed
308 309
		{return;}

incardon's avatar
incardon committed
310 311 312 313 314 315 316 317
		Box<dim,long int> sub_domain;

		if (use_bx_def == false)
		{
			sub_domain = gdb_ext.get(sub_id).Dbox;
			sub_domain += gdb_ext.get(sub_id).origin;
		}

incardon's avatar
incardon committed
318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343
		// Convert from SpaceBox<dim,St> to SpaceBox<dim,long int>
		Box<dim,long int> sub_domain_other_exp = cd_sm.convertDomainSpaceIntoGridUnits(sub_domain_other,dec.periodicity());

		// translate sub_domain_other based on cmb
		for (size_t i = 0 ; i < dim ; i++)
		{
			if (cmb.c[i] == 1)
			{
				sub_domain_other_exp.setLow(i,sub_domain_other_exp.getLow(i) - ginfo.size(i));
				sub_domain_other_exp.setHigh(i,sub_domain_other_exp.getHigh(i) - ginfo.size(i));
			}
			else if (cmb.c[i] == -1)
			{
				sub_domain_other_exp.setLow(i,sub_domain_other_exp.getLow(i) + ginfo.size(i));
				sub_domain_other_exp.setHigh(i,sub_domain_other_exp.getHigh(i) + ginfo.size(i));
			}
		}

		sub_domain_other_exp.enlarge(g);
		if (sub_domain_other_exp.Intersect(sub_domain,ib) == false)
		{
			for (size_t i = 0 ; i < dim ; i++)
			{ib.setHigh(i,ib.getLow(i) - 1);}
		}
	}

344 345


346
	/*! \brief Create per-processor internal ghost boxes list in grid units and g_id_to_external_ghost_box
347 348 349 350
	 *
	 */
	void create_ig_box()
	{
incardon's avatar
incardon committed
351 352 353
		// temporal vector used for computation
		openfpm::vector_std<result_box<dim>> ibv;

Pietro Incardona's avatar
Pietro Incardona committed
354 355
		if (init_i_g_box == true)	return;

356 357 358
		// Get the grid info
		auto g = cd_sm.getGrid();

Pietro Incardona's avatar
Pietro Incardona committed
359
		g_id_to_internal_ghost_box.resize(dec.getNNProcessors());
360 361 362 363 364 365 366 367 368 369 370

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

			pib.prc = dec.IDtoProc(i);
			for (size_t j = 0 ; j < dec.getProcessorNIGhost(i) ; j++)
			{
				// Get the internal ghost boxes and transform into grid units
incardon's avatar
incardon committed
371
				::Box<dim,St> ib_dom = dec.getProcessorIGhostBox(i,j);
372
				::Box<dim,long int> ib = cd_sm.convertDomainSpaceIntoGridUnits(ib_dom,dec.periodicity());
incardon's avatar
incardon committed
373

374 375 376 377 378
				size_t sub_id = dec.getProcessorIGhostSub(i,j);
				size_t r_sub = dec.getProcessorIGhostSSub(i,j);

				auto & n_box = dec.getNearSubdomains(dec.IDtoProc(i));

incardon's avatar
incardon committed
379
				set_for_adjustment(sub_id,
380 381 382
						           n_box.get(r_sub),dec.getProcessorIGhostPos(i,j),
						           ib,ghost_int);

incardon's avatar
incardon committed
383 384 385 386 387
				// Here we intersect the internal ghost box with the definition boxes
				// this operation make sense when the grid is not defined in the full
				// domain and we have to intersect the internal ghost box with all the box
				// that define where the grid is defined
				bx_intersect<dim>(bx_def,use_bx_def,ib,ibv);
incardon's avatar
incardon committed
388

incardon's avatar
incardon committed
389 390 391 392 393
				for (size_t k = 0 ; k < ibv.size() ; k++)
				{
					// Check if ib is valid if not it mean that the internal ghost does not contain information so skip it
					if (ibv.get(k).bx.isValid() == false)
					{continue;}
incardon's avatar
incardon committed
394

incardon's avatar
incardon committed
395 396
					// save the box and the sub-domain id (it is calculated as the linearization of P1)
					::Box<dim,size_t> cvt = ibv.get(k).bx;
incardon's avatar
incardon committed
397

incardon's avatar
incardon committed
398 399 400
					i_box_id<dim> bid_t;
					bid_t.box = cvt;
					bid_t.g_id = dec.getProcessorIGhostId(i,j) | (k) << 52;
incardon's avatar
incardon committed
401

incardon's avatar
incardon committed
402 403 404 405
					bid_t.sub = convert_to_gdb_ext(dec.getProcessorIGhostSub(i,j),
							                       ibv.get(k).id,
												   gdb_ext,
												   gdb_ext_markers);
406

incardon's avatar
incardon committed
407 408 409
					bid_t.cmb = dec.getProcessorIGhostPos(i,j);
					bid_t.r_sub = dec.getProcessorIGhostSSub(i,j);
					pib.bid.add(bid_t);
Pietro Incardona's avatar
Pietro Incardona committed
410

incardon's avatar
incardon committed
411 412
					g_id_to_internal_ghost_box.get(i)[bid_t.g_id | (k) << 52 ] = pib.bid.size()-1;
				}
413 414 415 416 417 418
			}
		}

		init_i_g_box = true;
	}

incardon's avatar
incardon committed
419

420 421 422 423 424 425 426 427 428 429
	/*! \brief Create per-processor internal ghost box list in grid units
	 *
	 */
	void create_eg_box()
	{
		// Get the grid info
		auto g = cd_sm.getGrid();

		if (init_e_g_box == true)	return;

430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452
		// Here we collect all the calculated internal ghost box in the sector different from 0 that this processor has

		openfpm::vector<size_t> prc;
		openfpm::vector<size_t> prc_recv;
		openfpm::vector<size_t> sz_recv;
		openfpm::vector<openfpm::vector<Box_fix<dim>>> box_int_send(dec.getNNProcessors());
		openfpm::vector<openfpm::vector<Box_fix<dim>>> box_int_recv;

		for(size_t i = 0 ; i < dec.getNNProcessors() ; i++)
		{
			for (size_t j = 0 ; j < ig_box.get(i).bid.size() ; j++)
			{
				box_int_send.get(i).add();
				box_int_send.get(i).last().bx = ig_box.get(i).bid.get(j).box;
				box_int_send.get(i).last().g_id = ig_box.get(i).bid.get(j).g_id;
				box_int_send.get(i).last().r_sub = ig_box.get(i).bid.get(j).r_sub;
				box_int_send.get(i).last().cmb = ig_box.get(i).bid.get(j).cmb;
			}
			prc.add(dec.IDtoProc(i));
		}

		v_cl.SSendRecv(box_int_send,box_int_recv,prc,prc_recv,sz_recv);

Pietro Incardona's avatar
Pietro Incardona committed
453
		eg_box.resize(dec.getNNProcessors());
454

Pietro Incardona's avatar
Pietro Incardona committed
455
		for (size_t i = 0 ; i < eg_box.size() ; i++)
incardon's avatar
incardon committed
456
		{eg_box.get(i).prc = dec.IDtoProc(i);}
457 458 459

		for (size_t i = 0 ; i < box_int_recv.size() ; i++)
		{
incardon's avatar
incardon committed
460
			size_t pib_cnt = 0;
461
			size_t p_id = dec.ProctoID(prc_recv.get(i));
Pietro Incardona's avatar
Pietro Incardona committed
462
			auto&& pib = eg_box.get(p_id);
463 464
			pib.prc = prc_recv.get(i);

incardon's avatar
incardon committed
465 466 467
			eg_box.get(p_id).recv_pnt = 0;
			eg_box.get(p_id).n_r_box = box_int_recv.get(i).size();

468 469 470
			// For each received internal ghost box
			for (size_t j = 0 ; j < box_int_recv.get(i).size() ; j++)
			{
incardon's avatar
incardon committed
471 472 473
				size_t vol_recv = box_int_recv.get(i).get(j).bx.getVolumeKey();

				eg_box.get(p_id).recv_pnt += vol_recv;
474 475
				size_t send_list_id = box_int_recv.get(i).get(j).r_sub;

incardon's avatar
incardon committed
476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492
				if (use_bx_def == true)
				{
					// First we understand if the internal ghost box sent intersect
					// some local extended sub-domain.

					// eb_gid_list, for an explanation check the declaration
					eb_gid_list.add();
					eb_gid_list.last().e_id = p_id;

					// Now we have to check if a received external ghost box intersect one
					// or more sub-grids
					for (size_t k = 0 ; k < gdb_ext.size() ; k++)
					{
						Box<dim,long int> bx = gdb_ext.get(k).GDbox;
						bx += gdb_ext.get(k).origin;

						Box<dim,long int> output;
493
						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);
incardon's avatar
incardon committed
494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545

						// it intersect one sub-grid
						if (bx.Intersect(flp_i,output))
						{
							// link

							size_t g_id = box_int_recv.get(i).get(j).g_id;
							add_eg_box<dim>(k,box_int_recv.get(i).get(j).cmb,output,
									g_id,
									gdb_ext.get(k).origin,
									box_int_recv.get(i).get(j).bx.getP1(),
									pib.bid);

							eb_gid_list.last().eb_list.add(pib.bid.size() - 1);

							g_id_to_external_ghost_box[g_id] = eb_gid_list.size() - 1;
						}
					}

					// now we check if exist a full match across the full intersected
					// ghost parts

					bool no_match = true;
					for (size_t k = 0 ; k < eb_gid_list.last().eb_list.size() ; k++)
					{
						size_t eb_id = eb_gid_list.last().eb_list.get(k);

						if (pib.bid.get(eb_id).g_e_box == box_int_recv.get(i).get(j).bx)
						{
							// full match found

							eb_gid_list.last().full_match = eb_id;
							no_match = false;

							break;
						}
					}

					// This is the case where a full match has not been found. In this case we
					// generate an additional gdb_ext and local grid with only external ghost

					if (no_match == true)
					{
						// Create a grid with the same size of the external ghost

						size_t sz[dim];
						for (size_t s = 0 ; s < dim ; s++)
						{sz[s] = box_int_recv.get(i).get(j).bx.getHigh(s) - box_int_recv.get(i).get(j).bx.getLow(s) + 1;}

						// Add an unlinked gdb_ext
						// An unlinked gdb_ext is an empty domain with only a external ghost
						// part
546
						Box<dim,long int> output = flip_box(box_int_recv.get(i).get(j).bx,box_int_recv.get(i).get(j).cmb,ginfo);
incardon's avatar
incardon committed
547 548 549 550 551

						GBoxes<dim> tmp;
						tmp.GDbox = box_int_recv.get(i).get(j).bx;
						tmp.GDbox -= tmp.GDbox.getP1();
						tmp.origin = output.getP1();
552
						for (size_t s = 0 ; s < dim ; s++)
incardon's avatar
incardon committed
553 554
						{
							// we set an invalid box, there is no-domain
555 556
							tmp.Dbox.setLow(s,0);
							tmp.Dbox.setHigh(s,-1);
incardon's avatar
incardon committed
557
						}
incardon's avatar
incardon committed
558
						tmp.k = (size_t)-1;
incardon's avatar
incardon committed
559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605
						gdb_ext.add(tmp);

						// create the local grid

						loc_grid.add();
						loc_grid.last().resize(sz);

						// Add an external ghost box

						size_t g_id = box_int_recv.get(i).get(j).g_id;
						add_eg_box<dim>(gdb_ext.size()-1,box_int_recv.get(i).get(j).cmb,output,
								g_id,
								gdb_ext.get(gdb_ext.size()-1).origin,
								box_int_recv.get(i).get(j).bx.getP1(),
								pib.bid);

						// now we map the received ghost box to the information of the
						// external ghost box created
						eb_gid_list.last().full_match = pib.bid.size() - 1;
						eb_gid_list.last().eb_list.add(pib.bid.size() - 1);
						g_id_to_external_ghost_box[g_id] = eb_gid_list.size() - 1;
					}

					// now we create lr_e_box

					size_t fm = eb_gid_list.last().full_match;
					size_t sub_id = pib.bid.get(fm).sub;

					for ( ; pib_cnt < pib.bid.size() ; pib_cnt++)
					{
						pib.bid.get(pib_cnt).lr_e_box -= gdb_ext.get(sub_id).origin;
					}

				}
				else
				{
					// Get the list of the sent sub-domains
					// and recover the id of the sub-domain from
					// the sent list
					const openfpm::vector<size_t> & s_sub = dec.getSentSubdomains(p_id);

					size_t sub_id = s_sub.get(send_list_id);

					e_box_id<dim> bid_t;
					bid_t.sub = sub_id;
					bid_t.cmb = box_int_recv.get(i).get(j).cmb;
					bid_t.cmb.sign_flip();
606
					::Box<dim,long int> ib = flip_box(box_int_recv.get(i).get(j).bx,box_int_recv.get(i).get(j).cmb,ginfo);
incardon's avatar
incardon committed
607 608 609 610 611 612 613 614 615 616 617 618 619 620 621
					bid_t.g_e_box = ib;
					bid_t.g_id = box_int_recv.get(i).get(j).g_id;
					// Translate in local coordinate
					Box<dim,long int> tb = ib;
					tb -= gdb_ext.get(sub_id).origin;
					bid_t.l_e_box = tb;

					pib.bid.add(bid_t);
					eb_gid_list.add();
					eb_gid_list.last().eb_list.add(pib.bid.size()-1);
					eb_gid_list.last().e_id = p_id;
					eb_gid_list.last().full_match = pib.bid.size()-1;

					g_id_to_external_ghost_box[bid_t.g_id] = eb_gid_list.size()-1;
				}
622 623 624
			}
		}

Pietro Incardona's avatar
Pietro Incardona committed
625
		init_e_g_box = true;
626 627
	}

incardon's avatar
incardon committed
628 629 630 631 632
	/*! \brief Create local internal ghost box in grid units
	 *
	 */
	void create_local_ig_box()
	{
incardon's avatar
incardon committed
633 634
		openfpm::vector_std<result_box<dim>> ibv;

incardon's avatar
incardon committed
635 636 637 638 639
		// Get the grid info
		auto g = cd_sm.getGrid();

		if (init_local_i_g_box == true)	return;

640
		// Get the number of sub-domains
Pietro Incardona's avatar
Pietro Incardona committed
641
		for (size_t i = 0 ; i < dec.getNSubDomain() ; i++)
incardon's avatar
incardon committed
642 643 644 645 646 647
		{
			loc_ig_box.add();
			auto&& pib = loc_ig_box.last();

			for (size_t j = 0 ; j < dec.getLocalNIGhost(i) ; j++)
			{
incardon's avatar
incardon committed
648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663
				if (use_bx_def == true)
				{
					// Get the internal ghost boxes and transform into grid units
					::Box<dim,St> ib_dom = dec.getLocalIGhostBox(i,j);
					::Box<dim,long int> ib = cd_sm.convertDomainSpaceIntoGridUnits(ib_dom,dec.periodicity());

					// Here we intersect the internal ghost box with the definition boxes
					// this operation make sense when the grid is not defined in the full
					// domain and we have to intersect the internal ghost box with all the box
					// that define where the grid is defined
					bx_intersect<dim>(bx_def,use_bx_def,ib,ibv);

					for (size_t k = 0 ; k < ibv.size() ; k++)
					{
						// Check if ib is valid if not it mean that the internal ghost does not contain information so skip it
						if (ibv.get(k).bx.isValid() == false)
664
						{continue;}
incardon's avatar
incardon committed
665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689

						pib.bid.add();
						pib.bid.last().box = ibv.get(k).bx;

						pib.bid.last().sub_gdb_ext = convert_to_gdb_ext(i,
																ibv.get(k).id,
																gdb_ext,
																gdb_ext_markers);

						pib.bid.last().sub = dec.getLocalIGhostSub(i,j);

						// It will be filled later
						pib.bid.last().cmb = dec.getLocalIGhostPos(i,j);
					}
				}
				else
				{
					// Get the internal ghost boxes and transform into grid units
					::Box<dim,St> ib_dom = dec.getLocalIGhostBox(i,j);
					::Box<dim,long int> ib = cd_sm.convertDomainSpaceIntoGridUnits(ib_dom,dec.periodicity());

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

690 691 692
					size_t sub_id = i;
					size_t r_sub = dec.getLocalIGhostSub(i,j);

incardon's avatar
incardon committed
693
					set_for_adjustment(sub_id,dec.getSubDomain(r_sub),
694 695 696 697 698 699
							           dec.getLocalIGhostPos(i,j),ib,ghost_int);

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

incardon's avatar
incardon committed
700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720
					pib.bid.add();
					pib.bid.last().box = ib;
					pib.bid.last().sub = dec.getLocalIGhostSub(i,j);
					pib.bid.last().k.add(dec.getLocalIGhostE(i,j));
					pib.bid.last().cmb = dec.getLocalIGhostPos(i,j);
					pib.bid.last().sub_gdb_ext = i;
				}
			}

			if (use_bx_def == true)
			{
				// unfortunately boxes that define where the grid is located can generate
				// additional internal ghost boxes

				for (size_t j = gdb_ext_markers.get(i) ; j < gdb_ext_markers.get(i+1) ; j++)
				{
					// intersect within box in the save sub-domain

					for (size_t k = gdb_ext_markers.get(i) ; k < gdb_ext_markers.get(i+1) ; k++)
					{
						if (j == k)	{continue;}
incardon's avatar
incardon committed
721

incardon's avatar
incardon committed
722 723 724 725 726
						// extend k and calculate the internal ghost box
						Box<dim,long int> bx_e =  gdb_ext.get(k).GDbox;
						bx_e += gdb_ext.get(k).origin;
						Box<dim,long int> bx = gdb_ext.get(j).Dbox;
						bx += gdb_ext.get(j).origin;
incardon's avatar
incardon committed
727

incardon's avatar
incardon committed
728 729 730 731
						Box<dim,long int> output;
						if (bx.Intersect(bx_e, output) == true)
						{
							pib.bid.add();
incardon's avatar
incardon committed
732

incardon's avatar
incardon committed
733
							pib.bid.last().box = output;
incardon's avatar
incardon committed
734

incardon's avatar
incardon committed
735 736
							pib.bid.last().sub_gdb_ext = j;
							pib.bid.last().sub = i;
incardon's avatar
incardon committed
737

incardon's avatar
incardon committed
738 739 740 741 742 743
//							if (use_bx_def == true)
//							{pib.bid.last().k = -1;}
//							else
//							{pib.bid.last().k = dec.getLocalIGhostE(i,j);}
							// these ghost always in the quadrant zero
							pib.bid.last().cmb.zero();
incardon's avatar
incardon committed
744

incardon's avatar
incardon committed
745 746 747
						}
					}
				}
incardon's avatar
incardon committed
748
			}
incardon's avatar
incardon committed
749

incardon's avatar
incardon committed
750 751
		}

incardon's avatar
incardon committed
752

incardon's avatar
incardon committed
753 754 755
		init_local_i_g_box = true;
	}

756
	/*! \brief Create per-processor external ghost boxes list in grid units
incardon's avatar
incardon committed
757 758 759 760 761 762 763 764 765
	 *
	 */
	void create_local_eg_box()
	{
		// Get the grid info
		auto g = cd_sm.getGrid();

		if (init_local_e_g_box == true)	return;

Pietro Incardona's avatar
Pietro Incardona committed
766
		loc_eg_box.resize(dec.getNSubDomain());
767 768 769 770 771 772

		// Get the number of sub-domain
		for (size_t i = 0 ; i < dec.getNSubDomain() ; i++)
		{
			for (size_t j = 0 ; j < loc_ig_box.get(i).bid.size() ; j++)
			{
incardon's avatar
incardon committed
773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791
				long int volume_linked = 0;

				size_t le_sub = loc_ig_box.get(i).bid.get(j).sub;
				auto & pib = loc_eg_box.get(le_sub);

				if (use_bx_def == true)
				{

					// We check if an external local ghost box intersect one
					// or more sub-grids
					for (size_t k = 0 ; k < gdb_ext.size() ; k++)
					{
						Box<dim,long int> bx = gdb_ext.get(k).Dbox;
						bx += gdb_ext.get(k).origin;

						Box<dim,long int> gbx = gdb_ext.get(k).GDbox;
						gbx += gdb_ext.get(k).origin;

						Box<dim,long int> output;
792
						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);
incardon's avatar
incardon committed
793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811

						bool intersect_domain = bx.Intersect(flp_i,output);
						bool intersect_gdomain = gbx.Intersect(flp_i,output);

						// it intersect one sub-grid
						if (intersect_domain == false && intersect_gdomain == true)
						{
							// fill the link variable
							loc_ig_box.get(i).bid.get(j).k.add(pib.bid.size());
							size_t s = loc_ig_box.get(i).bid.get(j).k.last();

							comb<dim> cmb = loc_ig_box.get(i).bid.get(j).cmb;
							cmb.sign_flip();

							add_loc_eg_box(le_sub,
										   loc_ig_box.get(i).bid.get(j).sub_gdb_ext,
										   j,
										   k,
										   pib.bid,
812
										   output,
incardon's avatar
incardon committed
813 814 815
										   cmb);


816
							volume_linked += pib.bid.last().ebox.getVolumeKey();
incardon's avatar
incardon committed
817 818 819 820 821 822 823 824 825 826 827
						}
					}
				}
				else
				{
					size_t k = loc_ig_box.get(i).bid.get(j).sub;
					auto & pib = loc_eg_box.get(k);

					size_t s = loc_ig_box.get(i).bid.get(j).k.get(0);
					pib.bid.resize(dec.getLocalNEGhost(k));

828
					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);
incardon's avatar
incardon committed
829 830 831 832 833 834 835
					pib.bid.get(s).sub = dec.getLocalEGhostSub(k,s);
					pib.bid.get(s).cmb = loc_ig_box.get(i).bid.get(j).cmb;
					pib.bid.get(s).cmb.sign_flip();
					pib.bid.get(s).k = j;
					pib.bid.get(s).initialized = true;
					pib.bid.get(s).sub_gdb_ext = k;
				}
836 837 838
			}
		}

Pietro Incardona's avatar
Pietro Incardona committed
839
		init_local_e_g_box = true;
incardon's avatar
incardon committed
840 841
	}

incardon's avatar
incardon committed
842

843
	/*! \brief Check the grid has a valid size
844 845
	 *
	 * \param g_sz size of the grid
incardon's avatar
incardon committed
846 847 848 849 850 851 852
	 *
	 */
	inline void check_size(const size_t (& g_sz)[dim])
	{
		for (size_t i = 0 ; i < dim ; i++)
		{
			if (g_sz[i] < 2)
incardon's avatar
incardon committed
853
				std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " distributed grids with size smaller than 2 are not supported\n";
incardon's avatar
incardon committed
854 855
		}
	}
incardon's avatar
incardon committed
856

incardon's avatar
incardon committed
857
	/*! \brief Create the grids on memory
incardon's avatar
incardon committed
858 859 860
	 *
	 * \param bx_def Where the grid is defined
	 * \param use_bx_def use the array that define where the grid is defined
incardon's avatar
incardon committed
861 862
	 *
	 */
incardon's avatar
incardon committed
863 864 865
	void Create(openfpm::vector<Box<dim,long int>> & bx_def,
			    const Ghost<dim,long int> & g,
			    bool use_bx_def)
incardon's avatar
incardon committed
866
	{
867
		// create gdb
incardon's avatar
incardon committed
868 869 870
		create_gdb_ext<dim,Decomposition>(gdb_ext,gdb_ext_markers,dec,cd_sm,bx_def,g,use_bx_def);

		size_t n_grid = gdb_ext.size();
871

incardon's avatar
incardon committed
872
		// create local grids for each hyper-cube
873
		loc_grid.resize(n_grid);
incardon's avatar
incardon committed
874 875 876 877 878 879 880

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

		// Allocate the grids
		for (size_t i = 0 ; i < n_grid ; i++)
		{
881
			SpaceBox<dim,long int> sp_tg = gdb_ext.get(i).GDbox;
incardon's avatar
incardon committed
882 883

			// Get the size of the local grid
884 885 886 887
			// The boxes indicate the extension of the index the size
			// is this extension +1
			// for example a 1D box (interval) from 0 to 3 in one dimension have
			// the points 0,1,2,3 = so a total of 4 points
888
			for (size_t j = 0 ; j < dim ; j++)
incardon's avatar
incardon committed
889
			{l_res[j] = (sp_tg.getHigh(j) >= 0)?(sp_tg.getHigh(j)+1):0;}
incardon's avatar
incardon committed
890 891 892 893 894 895 896

			// Set the dimensions of the local grid
			loc_grid.get(i).resize(l_res);
		}
	}


897 898 899 900 901 902
    /*! \brief Initialize the Cell decomposer of the grid enforcing perfect overlap of the cells
	 *
	 * \param cd_old the CellDecomposer we are trying to mach
	 * \param ext extension of the domain
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
903
	inline void InitializeCellDecomposer(const CellDecomposer_sm<dim,St,shift<dim,St>> & cd_old, const Box<dim,size_t> & ext)
904 905 906 907 908
	{
		// Initialize the cell decomposer
		cd_sm.setDimensions(cd_old,ext);
	}

909
    /*! \brief Initialize the Cell decomposer of the grid
incardon's avatar
incardon committed
910
	 *
911 912
	 * \param g_sz Size of the grid
	 * \param bc boundary conditions
incardon's avatar
incardon committed
913 914
	 *
	 */
915
	inline void InitializeCellDecomposer(const size_t (& g_sz)[dim], const size_t (& bc)[dim])
incardon's avatar
incardon committed
916
	{
incardon's avatar
incardon committed
917
		// check that the grid has valid size
incardon's avatar
incardon committed
918 919
		check_size(g_sz);

920
		// get the size of the cell decomposer
incardon's avatar
incardon committed
921
		size_t c_g[dim];
922
		getCellDecomposerPar<dim>(c_g,g_sz,bc);
incardon's avatar
incardon committed
923 924 925

		// Initialize the cell decomposer
		cd_sm.setDimensions(domain,c_g,0);
incardon's avatar
incardon committed
926
	}
incardon's avatar
incardon committed
927

incardon's avatar
incardon committed
928 929 930
	/*! \brief Initialize the grid
	 *
	 * \param g_sz Global size of the grid
931
	 * \param bc boundary conditions
incardon's avatar
incardon committed
932 933
	 *
	 */
934
	inline void InitializeDecomposition(const size_t (& g_sz)[dim], const size_t (& bc)[dim])
incardon's avatar
incardon committed
935
	{
incardon's avatar
incardon committed
936
		// fill the global size of the grid
incardon's avatar
incardon committed
937
		for (size_t i = 0 ; i < dim ; i++)	{this->g_sz[i] = g_sz[i];}
incardon's avatar
incardon committed
938 939 940 941

		// Get the number of processor and calculate the number of sub-domain
		// for decomposition
		size_t n_proc = v_cl.getProcessingUnits();
incardon's avatar
incardon committed
942
		size_t n_sub = n_proc * v_sub_unit_factor;
incardon's avatar
incardon committed
943 944 945 946

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

		// Create the sub-domains
951
		dec.setParameters(div,domain,bc,ghost);
Pietro Incardona's avatar
Pietro Incardona committed
952
		dec.decompose();
incardon's avatar
incardon committed
953 954
	}

incardon's avatar
incardon committed
955 956 957 958 959 960
	/*! \brief Initialize the grid
	 *
	 * \param g_sz Global size of the grid
	 *
	 */
	inline void InitializeStructures(const size_t (& g_sz)[dim])
incardon's avatar
incardon committed
961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981
	{
		// an empty
		openfpm::vector<Box<dim,long int>> empty;

		// Ghost zero
		Ghost<dim,long int> zero;

		InitializeStructures(g_sz,empty,zero,false);
	}

	/*! \brief Initialize the grid
	 *
	 * \param g_sz Global size of the grid
	 * \param g ghost extension of the grid in integer unit
	 * \param bx set of boxes that define where is defined the grid
	 *
	 */
	inline void InitializeStructures(const size_t (& g_sz)[dim],
			                         openfpm::vector<Box<dim,long int>> & bx,
									 const Ghost<dim,long int> & g,
									 bool use_bx_def)
incardon's avatar
incardon committed
982 983 984 985 986
	{
		// fill the global size of the grid
		for (size_t i = 0 ; i < dim ; i++)	{this->g_sz[i] = g_sz[i];}

		// Create local grid
incardon's avatar
incardon committed
987
		Create(bx,g,use_bx_def);
incardon's avatar
incardon committed
988 989
	}

incardon's avatar
incardon committed
990 991 992
	// Ghost as integer
	Ghost<dim,long int> gint = Ghost<dim,long int>(0);

993
protected:
994 995

	/*! \brief Given a local sub-domain i with a local grid Domain + ghost return the part of the local grid that is domain
996 997
	 *
	 * \param i sub-domain
998 999 1000 1001 1002 1003 1004 1005 1006
	 *
	 * \return the Box defining the domain in the local grid
	 *
	 */
	Box<dim,size_t> getDomain(size_t i)
	{
		return gdb_ext.get(i).Dbox;
	}

1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027
	/*! \brief Convert a ghost from grid point units into continus space
	 *
	 * \param gd Ghost in continuous space
	 * \param cd_sm CellDecomposer of the grid
	 *
	 * \return the ghost in continuous unit
	 *
	 */
	static inline Ghost<dim,float> convert_ghost(const Ghost<dim,long int> & gd, const CellDecomposer_sm<dim,St,shift<dim,St>> & cd_sm)
	{
		Ghost<dim,float> gc;

		// get the grid spacing
		Box<dim,St> sp = cd_sm.getCellBox();

		// enlarge 0.001 of the spacing
		sp.magnify_fix_P1(1.1);

		// set the ghost
		for (size_t i = 0 ; i < dim ; i++)
		{
incardon's avatar
incardon committed
1028 1029
			gc.setLow(i,gd.getLow(i)*(sp.getHigh(i)));
			gc.setHigh(i,gd.getHigh(i)*(sp.getHigh(i)));
1030 1031 1032 1033 1034
		}

		return gc;
	}

incardon's avatar
incardon committed
1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061
	/*! \brief Set the minimum number of sub-domain per processor
	 *
	 * \param n_sub
	 *
	 */
	void setDecompositionGranularity(size_t n_sub)
	{
		this->v_sub_unit_factor = n_sub;
	}

	void reset_ghost_structures()
	{
		g_id_to_internal_ghost_box.clear();
		ig_box.clear();
		init_i_g_box = false;

		eg_box.clear();
		eb_gid_list.clear();
		init_e_g_box = false;

		init_local_i_g_box = false;
		loc_ig_box.clear();

		init_local_e_g_box = false;
		loc_eg_box.clear();
	}

incardon's avatar
incardon committed
1062 1063
public:

1064
	//! Which kind of grid the structure store
incardon's avatar
incardon committed
1065 1066
	typedef device_grid d_grid;

1067
	//! Decomposition used
incardon's avatar
incardon committed
1068 1069
	typedef Decomposition decomposition;

1070
	//! value_type
1071 1072
	typedef T value_type;

1073
	//! Type of space
1074 1075
	typedef St stype;

1076
	//! Type of Memory
1077 1078
	typedef Memory memory_type;

1079
	//! Type of device grid
1080 1081
	typedef device_grid device_grid_type;

1082
	//! Number of dimensions
1083 1084 1085 1086
	static const unsigned int dims = dim;

	/*! \brief Get the domain where the grid is defined
	 *
1087
	 * \return the domain of the grid
1088 1089 1090 1091 1092 1093 1094
	 *
	 */
	inline const Box<dim,St> getDomain() const
	{
		return domain;
	}

incardon's avatar
incardon committed
1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106
	/*! \brief Get the point where it start the origin of the grid of the sub-domain i
	 *
	 * \param i sub-domain
	 *
	 * \return the point
	 *
	 */
	Point<dim,St> getOffset(size_t i)
	{
		return pmul(Point<dim,St>(gdb_ext.get(i).origin), cd_sm.getCellBox().getP2()) + getDomain().getP1();
	}

1107
    /*! \brief Get the spacing of the grid in direction i
1108 1109
     *
     * \param i dimension
1110 1111 1112 1113 1114 1115 1116 1117 1118
     *
     * \return the spacing
     *
     */
    inline St spacing(size_t i) const
    {
    	return cd_sm.getCellBox().getHigh(i);
    }

1119 1120 1121 1122 1123 1124 1125 1126 1127 1128
	/*! \brief Return the total number of points in the grid
	 *
	 * \return number of points
	 *
	 */
	size_t size() const
	{
		return ginfo_v.size();
	}

incardon's avatar
incardon committed
1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140
	/*! \brief set the background value
	 *
	 * You can use this function make sense in case of sparse in case of dense
	 * it does nothing
	 *
	 */
	void setBackgroundValue(T & bv)
	{
		for (size_t i = 0 ; i < loc_grid.size() ; i++)
		{meta_copy<T>::meta_copy_(bv,loc_grid.get(i).getBackgroundValue());}
	}

1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153
	/*! \brief set the background value
	 *
	 * You can use this function make sense in case of sparse in case of dense
	 * it does nothing
	 *
	 */
	template<unsigned int p>
	void setBackgroundValue(const typename boost::mpl::at<typename T::type,boost::mpl::int_<p>>::type & bv)
	{
		for (size_t i = 0 ; i < loc_grid.size() ; i++)
		{loc_grid.get(i).template setBackgroundValue<p>(bv);}
	}

incardon's avatar
incardon committed
1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171
	/*! \brief Return the local total number of points inserted in the grid
	 *
	 * in case of dense grid it return the number of local points, in case of
	 * sparse it return the number of inserted points
	 *
	 * \return number of points
	 *
	 */
	size_t size_local_inserted() const
	{
		size_t lins = 0;

		for (size_t i = 0 ; i < loc_grid.size() ; i++)
		{lins += loc_grid.get(i).size_inserted();}

		return lins;
	}

Pietro Incardona's avatar
Pietro Incardona committed
1172 1173 1174 1175 1176 1177 1178 1179 1180
	/*! \brief Return the total number of points in the grid
	 *
	 * \param i direction
	 *
	 * \return number of points on direction i
	 *
	 */
	size_t size(size_t i) const
	{
1181
		return ginfo_v.size(i);
Pietro Incardona's avatar
Pietro Incardona committed
1182 1183
	}

incardon's avatar
incardon committed
1184 1185 1186 1187 1188 1189 1190 1191 1192
	/*! \brief Default Copy constructor on this class make no sense and is unsafe, this definition disable it
	 *
	 * \param g grid to copy
	 *
	 */
	grid_dist_id(const grid_dist_id<dim,St,T,Decomposition,Memory,device_grid> & g)
	:grid_dist_id_comm<dim,St,T,Decomposition,Memory,device_grid>(g),
	 domain(g.domain),
	 ghost(g.ghost),
1193
	 ghost_int(g.ghost_int),
incardon's avatar
incardon committed
1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218
	 loc_grid(g.loc_grid),
	 loc_grid_old(g.loc_grid_old),
	 dec(g.dec),
	 gdb_ext_markers(g.gdb_ext_markers),
	 gdb_ext(g.gdb_ext),
	 gdb_ext_global(g.gdb_ext_global),
	 gdb_ext_old(g.gdb_ext_old),
	 cd_sm(g.cd_sm),
	 v_cl(g.v_cl),
	 prp_names(g.prp_names),
	 g_id_to_external_ghost_box(g.g_id_to_external_ghost_box),
	 g_id_to_internal_ghost_box(g.g_id_to_internal_ghost_box),
	 ginfo(g.ginfo),
	 ginfo_v(g.ginfo_v),
	 init_local_i_g_box(g.init_local_i_g_box),
	 init_local_e_g_box(g.init_local_e_g_box)
	{
#ifdef SE_CLASS2
		check_new(this,8,GRID_DIST_EVENT,4);
#endif

		for (size_t i = 0 ; i < dim ; i++)
		{g_sz[i] = g.g_sz[i];}
	}

1219 1220 1221 1222 1223 1224
	/*! \brief This constructor is special, it construct an expanded grid that perfectly overlap with the previous
	 *
	 * The key-word here is "perfectly overlap". Using the default constructor you could create
	 * something similar, but because of rounding-off error it can happen that it is not perfectly overlapping
	 *
	 * \param g previous grid
1225
	 * \param gh Ghost part in grid units
1226 1227 1228
	 * \param ext extension of the grid (must be positive on every direction)
	 *
	 */
incardon's avatar
incardon committed
1229 1230 1231 1232 1233
	template<typename H>
	grid_dist_id(const grid_dist_id<dim,St,H,typename Decomposition::base_type,Memory,grid_cpu<dim,H>> & g,
			     const Ghost<dim,long int> & gh,
				 Box<dim,size_t> ext)
	:ghost_int(gh),dec(create_vcluster()),v_cl(create_vcluster())
1234 1235 1236 1237 1238
	{
#ifdef SE_CLASS2
		check_new(this,8,GRID_DIST_EVENT,4);
#endif

1239 1240 1241 1242 1243 1244 1245 1246 1247 1248
		size_t ext_dim[dim];
		for (size_t i = 0 ; i < dim ; i++) {ext_dim[i] = g.getGridInfoVoid().size(i) + ext.getKP1().get(i) + ext.getKP2().get(i);}

		// Set the grid info of the extended grid
		ginfo.setDimensions(ext_dim);
		ginfo_v.setDimensions(ext_dim);

		InitializeCellDecomposer(g.getCellDecomposer(),ext);

		ghost = convert_ghost(gh,cd_sm);
1249 1250 1251 1252 1253

		// Extend the grid by the extension part and calculate the domain

		for (size_t i = 0 ; i < dim ; i++)
		{
1254
			g_sz[i] = g.size(i) + ext.getLow(i) + ext.getHigh(i);
1255

1256 1257 1258 1259 1260 1261 1262 1263 1264 1265
			if (g.getDecomposition().periodicity(i) == NON_PERIODIC)
			{
				this->domain.setLow(i,g.getDomain().getLow(i) - ext.getLow(i) * g.spacing(i) - g.spacing(i) / 2.0);
				this->domain.setHigh(i,g.getDomain().getHigh(i) + ext.getHigh(i) * g.spacing(i) + g.spacing(i) / 2.0);
			}
			else
			{
				this->domain.setLow(i,g.getDomain().getLow(i) - ext.getLow(i) * g.spacing(i));
				this->domain.setHigh(i,g.getDomain().getHigh(i) + ext.getHigh(i) * g.spacing(i));
			}
1266 1267
		}

1268
		dec.setParameters(g.getDecomposition(),ghost,this->domain);
Pietro Incardona's avatar
Pietro Incardona committed
1269

incardon's avatar
incardon committed
1270 1271 1272 1273
		// an empty
		openfpm::vector<Box<dim,long int>> empty;

		InitializeStructures(g.getGridInfoVoid().getSize(),empty,gh,false);
1274
	}
incardon's avatar
incardon committed
1275

1276 1277 1278 1279 1280 1281 1282
    /*! 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
     *
     * \param dec Decomposition
     * \param g_sz grid size on each dimension
     * \param ghost Ghost part
     *
     */
incardon's avatar
incardon committed
1283 1284 1285 1286 1287
    grid_dist_id(const Decomposition & dec,
    		     const size_t (& g_sz)[dim],
				 const Ghost<dim,St> & ghost)
    :domain(dec.getDomain()),ghost(ghost),ghost_int(INVALID_GHOST),dec(dec),v_cl(create_vcluster()),
	 ginfo(g_sz),ginfo_v(g_sz)
Pietro Incardona's avatar
Pietro Incardona committed
1288
	{
1289 1290 1291
#ifdef SE_CLASS2
		check_new(this,8,GRID_DIST_EVENT,4);
#endif
1292

Pietro Incardona's avatar
Pietro Incardona committed
1293
		InitializeCellDecomposer(g_sz,dec.periodicity());
incardon's avatar
incardon committed
1294 1295 1296

		this->dec = dec.duplicate(ghost);

Pietro Incardona's avatar
Pietro Incardona committed
1297 1298 1299
		InitializeStructures(g_sz);
	}

1300 1301 1302 1303 1304 1305 1306
    /*! 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
     *
     * \param dec Decomposition
     * \param g_sz grid size on each dimension
     * \param ghost Ghost part
     *
     */
incardon's avatar
incardon committed
1307 1308 1309 1310
    grid_dist_id(Decomposition && dec, const size_t (& g_sz)[dim],
    		     const Ghost<dim,St> & ghost)
    :domain(dec.getDomain()),ghost(ghost),dec(dec),ginfo(g_sz),
	 ginfo_v(g_sz),v_cl(create_vcluster()),ghost_int(INVALID_GHOST)
incardon's avatar
incardon committed
1311
	{
1312 1313 1314
#ifdef SE_CLASS2
		check_new(this,8,GRID_DIST_EVENT,4);
#endif
1315

Pietro Incardona's avatar
Pietro Incardona committed
1316
		InitializeCellDecomposer(g_sz,dec.periodicity());
incardon's avatar
incardon committed
1317 1318 1319

		this->dec = dec.duplicate(ghost);

incardon's avatar
incardon committed
1320
		InitializeStructures(g_sz);
incardon's avatar
incardon committed
1321 1322
	}

1323
    /*! 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
1324
     *
1325 1326
     * \param dec Decomposition
     * \param g_sz grid size on each dimension
1327
     * \param g Ghost part (given in grid units)
1328 1329
     *
     * \warning In very rare case the ghost part can be one point bigger than the one specified
1330 1331
     *
     */
incardon's avatar
incardon committed
1332 1333 1334 1335
	grid_dist_id(const Decomposition & dec, const size_t (& g_sz)[dim],
			     const Ghost<dim,long int> & g)
	:domain(dec.getDomain()),ghost_int(g),dec(create_vcluster()),v_cl(create_vcluster()),
	 ginfo(g_sz),ginfo_v(g_sz)
incardon's avatar
incardon committed
1336
	{
1337 1338 1339
#ifdef SE_CLASS2
		check_new(this,8,GRID_DIST_EVENT,4);
#endif
1340

Pietro Incardona's avatar
Pietro Incardona committed
1341
		InitializeCellDecomposer(g_sz,dec.periodicity());
1342 1343

		ghost = convert_ghost(g,cd_sm);
1344
		this->dec = dec.duplicate(ghost);
1345

incardon's avatar
incardon committed
1346 1347 1348
		// an empty
		openfpm::vector<Box<dim,long int>> empty;

1349
		// Initialize structures
incardon's avatar
incardon committed
1350
		InitializeStructures(g_sz,empty,g,false);
incardon's avatar
incardon committed
1351
	}
incardon's avatar
incardon committed
1352

1353 1354 1355 1356
    /*! 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
     *
     * \param dec Decomposition
     * \param g_sz grid size on each dimension
1357
     * \param g Ghost part (given in grid units)
1358 1359 1360 1361
     *
     * \warning In very rare case the ghost part can be one point bigger than the one specified
     *
     */
incardon's avatar
incardon committed
1362 1363 1364 1365
	grid_dist_id(Decomposition && dec, const size_t (& g_sz)[dim],
			     const Ghost<dim,long int> & g)
	:domain(dec.getDomain()),dec(dec),v_cl(create_vcluster()),ginfo(g_sz),
	 ginfo_v(g_sz),ghost_int(g)
incardon's avatar
incardon committed
1366
	{
1367 1368 1369
#ifdef SE_CLASS2
		check_new(this,8,GRID_DIST_EVENT,4);
#endif
Pietro Incardona's avatar
Pietro Incardona committed
1370
		InitializeCellDecomposer(g_sz,dec.periodicity());
incardon's avatar
incardon committed
1371

Pietro Incardona's avatar
Pietro Incardona committed
1372
		ghost = convert_ghost(g,cd_sm);
incardon's avatar
incardon committed
1373 1374 1375 1376
		this->dec = dec.duplicate(ghost);

		// an empty
		openfpm::vector<Box<dim,long int>> empty;
incardon's avatar
incardon committed
1377 1378

		// Initialize structures
incardon's avatar
incardon committed
1379
		InitializeStructures(g_sz,empty,g,false);
incardon's avatar
incardon committed
1380 1381
	}

1382 1383 1384 1385
    /*! It construct a grid of a specified size, defined on a specified Box space, and having a specified ghost size
     *
     * \param g_sz grid size on each dimension
     * \param domain Box that contain the grid
1386
     * \param g Ghost part (given in grid units)
1387 1388 1389 1390
     *
     * \warning In very rare case the ghost part can be one point bigger than the one specified
     *
     */
incardon's avatar
incardon committed
1391 1392
	grid_dist_id(const size_t (& g_sz)[dim],const Box<dim,St> & domain, const Ghost<dim,St> & g, size_t opt = 0)
	:grid_dist_id(g_sz,domain,g,create_non_periodic<dim>(),opt)
1393 1394 1395 1396 1397 1398 1399
	{
	}

    /*! It construct a grid of a specified size, defined on a specified Box space, having a specified ghost size and periodicity
     *
     * \param g_sz grid size on each dimension
     * \param domain Box that contain the grid
1400
     * \param g Ghost part of the domain (given in grid units)
1401 1402 1403 1404
     *
     * \warning In very rare case the ghost part can be one point bigger than the one specified
     *
     */
incardon's avatar
incardon committed
1405 1406
	grid_dist_id(const size_t (& g_sz)[dim],const Box<dim,St> & domain, const Ghost<dim,long int> & g, size_t opt = 0)
	:grid_dist_id(g_sz,domain,g,create_non_periodic<dim>(),opt)
1407 1408 1409 1410 1411 1412 1413
	{
	}

    /*! It construct a grid of a specified size, defined on a specified Box space, having a specified ghost size, and specified periodicity
     *
     * \param g_sz grid size on each dimension
     * \param domain Box that contain the grid
1414
     * \param g Ghost part (given in grid units)
1415 1416 1417 1418 1419
     * \param p Boundary conditions
     *
     * \warning In very rare case the ghost part can be one point bigger than the one specified
     *
     */
incardon's avatar
incardon committed
1420 1421
	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)
	:domain(domain),ghost(g),ghost_int(INVALID_GHOST),dec(create_vcluster()),v_cl(create_vcluster()),ginfo(g_sz),ginfo_v(g_sz)
incardon's avatar
incardon committed
1422
	{
Pietro Incardona's avatar
Pietro Incardona committed
1423 1424 1425
#ifdef SE_CLASS2
		check_new(this,8,GRID_DIST_EVENT,4);
#endif
incardon's avatar
incardon committed
1426

incardon's avatar
incardon committed
1427 1428 1429
		if (opt >> 32 != 0)
		{this->setDecompositionGranularity(opt >> 32);}

1430 1431
		InitializeCellDecomposer(g_sz,p.bc);
		InitializeDecomposition(g_sz, p.bc);
incardon's avatar
incardon committed
1432 1433 1434
		InitializeStructures(g_sz);
	}

incardon's avatar
incardon committed
1435

1436
    /*! It construct a grid of a specified size, defined on a specified Box space, having a specified ghost size and periodicity
1437 1438 1439
     *
     * \param g_sz grid size on each dimension
     * \param domain Box that contain the grid
1440 1441
     * \param g Ghost part of the domain (given in grid units)
     * \param p periodicity
1442 1443 1444 1445
     *
     * \warning In very rare case the ghost part can be one point bigger than the one specified
     *
     */
incardon's avatar
incardon committed
1446 1447
	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)
	:domain(domain),ghost_int(g),dec(create_vcluster()),v_cl(create_vcluster()),ginfo(g_sz),ginfo_v(g_sz)
incardon's avatar
incardon committed
1448
	{
1449 1450 1451
#ifdef SE_CLASS2
		check_new(this,8,GRID_DIST_EVENT,4);
#endif
incardon's avatar
incardon committed
1452 1453 1454 1455

		if (opt >> 32 != 0)
		{this->setDecompositionGranularity(opt >> 32);}

1456
		InitializeCellDecomposer(g_sz,p.bc);
incardon's avatar
incardon committed
1457

Pietro Incardona's avatar
Pietro Incardona committed
1458
		ghost = convert_ghost(g,cd_sm);
1459

1460
		InitializeDecomposition(g_sz,p.bc);
incardon's avatar
incardon committed
1461 1462 1463 1464

		// an empty
		openfpm::vector<Box<dim,long int>> empty;

incardon's avatar
incardon committed
1465
		// Initialize structures
incardon's avatar
incardon committed
1466 1467 1468 1469 1470 1471 1472 1473 1474 1475 1476 1477 1478 1479 1480 1481 1482 1483 1484 1485 1486 1487 1488 1489 1490 1491 1492 1493 1494 1495 1496 1497 1498 1499 1500 1501
		InitializeStructures(g_sz,empty,g,false);
	}

	/*! \brief It construct a grid on the full domain restricted
	 *         to the set of boxes specified
	 *
	 * In particular the grid is defined in the space equal to the
	 *  domain intersected the boxes defined by bx
	 *
	 * \param g_sz grid size on each dimension
	 * \param domain where the grid is constructed
	 * \param g ghost size
	 * \param p periodicity of the grid
	 * \param bx set of boxes where the grid is defined
	 *
	 *
	 */
	grid_dist_id(const size_t (& g_sz)[dim],
			     const Box<dim,St> & domain,
				 const Ghost<dim,long int> & g,
				 const periodicity<dim> & p,
				 openfpm::vector<Box<dim,long int>> & bx_def)
	:domain(domain),dec(create_vcluster()),v_cl(create_vcluster()),ginfo(g_sz),ginfo_v(g_sz),gint(g)
	{
#ifdef SE_CLASS2
		check_new(this,8,GRID_DIST_EVENT,4);
#endif

		InitializeCellDecomposer(g_sz,p.bc);

		ghost = convert_ghost(g,cd_sm);

		InitializeDecomposition(g_sz, p.bc);
		InitializeStructures(g_sz,bx_def,g,true);
		this->bx_def = bx_def;
		this->use_bx_def = true;
incardon's avatar
incardon committed
1502 1503
	}

incardon's avatar
incardon committed
1504 1505 1506 1507 1508
	/*! \brief Get an object containing the grid informations
	 *
	 * \return an information object about this grid
	 *
	 */
1509
	const grid_sm<dim,T> & getGridInfo() const
incardon's avatar
incardon committed
1510
	{
1511 1512 1513
#ifdef SE_CLASS2
		check_valid(this,8);
#endif
incardon's avatar
incardon committed
1514 1515 1516 1517 1518 1519 1520 1521
		return ginfo;
	}

	/*! \brief Get an object containing the grid informations without type
	 *
	 * \return an information object about this grid
	 *
	 */
1522
	const grid_sm<dim,void> & getGridInfoVoid() const
incardon's avatar
incardon committed
1523
	{
1524 1525 1526
#ifdef SE_CLASS2
		check_valid(this,8);
#endif
incardon's avatar
incardon committed
1527 1528 1529
		return ginfo_v;
	}

1530
	/*! \brief Get the object that store the information about the decomposition
incardon's avatar
incardon committed
1531 1532 1533 1534 1535 1536
	 *
	 * \return the decomposition object
	 *
	 */
	Decomposition & getDecomposition()
	{
1537 1538 1539
#ifdef SE_CLASS2
		check_valid(this,8);
#endif
incardon's avatar
incardon committed
1540 1541 1542
		return dec;
	}

Pietro Incardona's avatar
Pietro Incardona committed
1543 1544 1545 1546 1547 1548 1549 1550 1551 1552 1553 1554 1555
	/*! \brief Get the object that store the information about the decomposition
	 *
	 * \return the decomposition object
	 *
	 */
	const Decomposition & getDecomposition() const
	{
#ifdef SE_CLASS2
		check_valid(this,8);
#endif
		return dec;
	}

incardon's avatar
incardon committed
1556 1557 1558 1559 1560
	/*! \brief Return the cell decomposer
	 *
	 * \return the cell decomposer
	 *
	 */
1561
	const CellDecomposer_sm<dim,St,shift<dim,St>> & getCellDecomposer() const
incardon's avatar
incardon committed
1562
	{
1563 1564 1565
#ifdef SE_CLASS2
		check_valid(this,8);
#endif
incardon's avatar
incardon committed
1566 1567 1568 1569
		return cd_sm;
	}

	/*! \brief Check that the global grid key is inside the grid domain
1570 1571
	 *
	 * \param gk point to check
incardon's avatar
incardon committed
1572
	 *
incardon's avatar
incardon committed
1573
	 * \return true if is inside
incardon's avatar
incardon committed
1574
	 *
incardon's avatar
incardon committed
1575 1576 1577
	 */
	bool