CartDecomposition.hpp 51.5 KB
Newer Older
incardon's avatar
incardon committed
1 2 3
/*
 * CartDecomposition.hpp
 *
4
 *  Created on: Oct 07, 2015
5
 *      Author: Pietro Incardona, Antonio Leo
incardon's avatar
incardon committed
6 7 8 9 10 11
 */

#ifndef CARTDECOMPOSITION_HPP
#define CARTDECOMPOSITION_HPP

#include "config.h"
12
#include <cmath>
incardon's avatar
incardon committed
13
#include "VCluster/VCluster.hpp"
14
#include "Graph/CartesianGraphFactory.hpp"
incardon's avatar
incardon committed
15
#include "Decomposition.hpp"
incardon's avatar
incardon committed
16
#include "Vector/map_vector.hpp"
incardon's avatar
incardon committed
17 18 19 20 21
#include <vector>
#include <initializer_list>
#include "SubdomainGraphNodes.hpp"
#include "dec_optimizer.hpp"
#include "Space/Shape/Box.hpp"
incardon's avatar
incardon committed
22
#include "Space/Shape/Point.hpp"
incardon's avatar
incardon committed
23
#include "NN/CellList/CellDecomposer.hpp"
incardon's avatar
incardon committed
24 25
#include <unordered_map>
#include "NN/CellList/CellList.hpp"
incardon's avatar
incardon committed
26
#include "Space/Ghost.hpp"
27 28
#include "common.hpp"
#include "ie_loc_ghost.hpp"
29 30
#include "ie_ghost.hpp"
#include "nn_processor.hpp"
31
#include "GraphMLWriter/GraphMLWriter.hpp"
32 33 34 35
#include "Distribution/ParMetisDistribution.hpp"
#include "Distribution/DistParMetisDistribution.hpp"
#include "Distribution/MetisDistribution.hpp"
#include "DLB/DLB.hpp"
36
#include "util/se_util.hpp"
37
#include "util/mathutil.hpp"
Pietro Incardona's avatar
Pietro Incardona committed
38
#include "CartDecomposition_ext.hpp"
39 40
#include "data_type/aggregate.hpp"
#include "Domain_NN_calculator_cart.hpp"
incardon's avatar
incardon committed
41

42 43
#define CARTDEC_ERROR 2000lu

incardon's avatar
incardon committed
44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
/*! \brief It spread the sub-sub-domain on a regular cartesian grid of size dim
 *
 * \warning this function only guarantee that the division on each direction is
 *          2^n with some n and does not guarantee that the number of
 *          sub-sub-domain is preserved
 *
 * \param div number of division on each direction as output
 * \param n_sub number of sub-domain
 * \param dim_r dimension reduction
 *
 */
template<unsigned int dim> static void nsub_to_div2(size_t (& div)[dim], size_t n_sub, size_t dim_r)
{
	for (size_t i = 0; i < dim; i++)
	{
		if (i < dim_r)
		{div[i] = openfpm::math::round_big_2(pow(n_sub, 1.0 / dim_r));}
		else
		{div[i] = 1;}
	}
}

/*! \brief It spread the sub-sub-domain on a regular cartesian grid of size dim
 *
 * \warning this function only guarantee that the division on each direction is
 *          2^n with some n and does not guarantee that the number of
 *          sub-sub-domain is preserved
 *
 * \param div number of division on each direction as output
 * \param n_sub number of sub-domain
 * \param dim_r dimension reduction
 *
 */
template<unsigned int dim> static void nsub_to_div(size_t (& div)[dim], size_t n_sub, size_t dim_r)
{
	for (size_t i = 0; i < dim; i++)
	{
		if (i < dim_r)
		{div[i] = std::floor(pow(n_sub, 1.0 / dim_r));}
		else
		{div[i] = 1;}
	}
}

88 89
#define COMPUTE_SKIN_SUB 1

incardon's avatar
incardon committed
90
/**
Pietro Incardona's avatar
Pietro Incardona committed
91
 * \brief This class decompose a space into sub-sub-domains and distribute them across processors
incardon's avatar
incardon committed
92 93 94 95
 *
 * \tparam dim is the dimensionality of the physical domain we are going to decompose.
 * \tparam T type of the space we decompose, Real, Integer, Complex ...
 * \tparam Memory Memory factory used to allocate memory
96
 * \tparam Distribution type of distribution, can be ParMetisDistribution or MetisDistribution
incardon's avatar
incardon committed
97
 *
98
 * Given an N-dimensional space, this class decompose the space into a Cartesian grid of small
incardon's avatar
incardon committed
99 100 101 102
 * sub-sub-domain. To each sub-sub-domain is assigned an id that identify at which processor is
 * assigned (in general the union of all the sub-sub-domain assigned to a processor is
 * simply connected space), a second step merge several sub-sub-domain with same id into bigger region
 *  sub-domain. Each sub-domain has an extended space called ghost part
103 104 105 106
 *
 * Assuming that VCluster.getProcessUnitID(), equivalent to the MPI processor rank, return the processor local
 * processor id, we define
 *
107 108
 * * local processor: processor rank
 * * local sub-domain: sub-domain given to the local processor
109 110 111 112 113 114 115 116 117 118 119 120
 * * external ghost box: (or ghost box) are the boxes that compose the ghost space of the processor, or the
 *   boxes produced expanding every local sub-domain by the ghost extension and intersecting with the sub-domain
 *   of the other processors
 * * Near processors are the processors adjacent to the local processor, where with adjacent we mean all the processor
 *   that has a non-zero intersection with the ghost part of the local processor, or all the processors that
 *   produce non-zero external boxes with the local processor, or all the processor that should communicate
 *   in case of ghost data synchronization
 * * internal ghost box: is the part of ghost of the near processor that intersect the space of the
 *       processor, or the boxes produced expanding the sub-domain of the near processors with the local sub-domain
 * * Near processor sub-domain: is a sub-domain that live in the a near (or contiguous) processor
 * * Near processor list: the list of all the near processor of the local processor (each processor has a list
 *                        of the near processor)
Pietro Incardona's avatar
Pietro Incardona committed
121
 * * Local ghosts internal or external are all the ghosts that does not involve inter-processor communications
122 123
 *
 * \see calculateGhostBoxes() for a visualization of internal and external ghost boxes
incardon's avatar
incardon committed
124
 *
125
 * ### Create a Cartesian decomposition object on a Box space, distribute, calculate internal and external ghost boxes
Pietro Incardona's avatar
Pietro Incardona committed
126
 *
127 128
 * \snippet CartDecomposition_unit_test.hpp Create CartDecomposition
 *
incardon's avatar
incardon committed
129 130
 */

Pietro Incardona's avatar
Pietro Incardona committed
131
template<unsigned int dim, typename T, typename Memory, typename Distribution>
132
class CartDecomposition: public ie_loc_ghost<dim, T>, public nn_prcs<dim, T>, public ie_ghost<dim, T>, public domain_nn_calculator_cart<dim>
incardon's avatar
incardon committed
133 134
{
public:
incardon's avatar
incardon committed
135

incardon's avatar
incardon committed
136 137 138 139
	//! Type of the domain we are going to decompose
	typedef T domain_type;

	//! It simplify to access the SpaceBox element
140
	typedef SpaceBox<dim, T> Box;
incardon's avatar
incardon committed
141

Pietro Incardona's avatar
Pietro Incardona committed
142 143 144 145 146 147 148
	//! This class is base of itself
	typedef CartDecomposition<dim,T,Memory,Distribution> base_type;

	//! This class admit a class defined on an extended domain
	typedef CartDecomposition_ext<dim,T,Memory,Distribution> extended_type;

protected:
incardon's avatar
incardon committed
149

incardon's avatar
incardon committed
150 151 152
	//! Indicate the communication weight has been set
	bool commCostSet = false;

incardon's avatar
incardon committed
153
	//! This is the key type to access  data_s, for example in the case of vector
incardon's avatar
incardon committed
154
	//! acc_key is size_t
Pietro Incardona's avatar
Pietro Incardona committed
155 156 157 158 159 160
	typedef typename openfpm::vector<SpaceBox<dim, T>,
			Memory,
			typename memory_traits_lin<SpaceBox<dim, T>>::type,
			memory_traits_lin,
			openfpm::vector_grow_policy_default,
			openfpm::vect_isel<SpaceBox<dim, T>>::value>::access_key acc_key;
incardon's avatar
incardon committed
161 162

	//! the set of all local sub-domain as vector
163
	openfpm::vector<SpaceBox<dim, T>> sub_domains;
incardon's avatar
incardon committed
164

165 166
	//! the remote set of all sub-domains as vector of 'sub_domains' vectors
	mutable openfpm::vector<Box_map<dim, T>> sub_domains_global;
167

incardon's avatar
incardon committed
168
	//! for each sub-domain, contain the list of the neighborhood processors
incardon's avatar
incardon committed
169 170
	openfpm::vector<openfpm::vector<long unsigned int> > box_nn_processor;

incardon's avatar
incardon committed
171
	//! Structure that contain for each sub-sub-domain box the processor id
incardon's avatar
incardon committed
172
	//! exist for efficient global communication
173
	CellList<dim,T,Mem_fast<>,shift<dim,T>> fine_s;
incardon's avatar
incardon committed
174

incardon's avatar
incardon committed
175
	//! Structure that store the cartesian grid information
176
	grid_sm<dim, void> gr;
incardon's avatar
incardon committed
177

178 179 180
	//! Structure that store the cartesian grid information
	grid_sm<dim, void> gr_dist;

incardon's avatar
incardon committed
181
	//! Structure that decompose the space into cells without creating them
incardon's avatar
incardon committed
182
	//! useful to convert positions to CellId or sub-domain id in this case
183
	CellDecomposer_sm<dim, T, shift<dim,T>> cd;
incardon's avatar
incardon committed
184 185

	//! rectangular domain to decompose
186
	::Box<dim,T> domain;
incardon's avatar
incardon committed
187 188 189 190

	//! Box Spacing
	T spacing[dim];

191 192 193 194
	//! Magnification factor between distribution and
	//! decomposition
	size_t magn[dim];

incardon's avatar
incardon committed
195 196 197
	//! Runtime virtual cluster machine
	Vcluster & v_cl;

198
	//! Create distribution
199 200
	Distribution dist;

201
	//! Processor bounding box
202
	::Box<dim,T> bbox;
incardon's avatar
incardon committed
203

204
	//! reference counter of the object in case is shared between object
205 206
	long int ref_cnt;

207
	//! ghost info
208 209
	Ghost<dim,T> ghost;

210
	//! Boundary condition info
211
	size_t bc[dim];
incardon's avatar
incardon committed
212

213 214 215 216 217 218
	//! Processor domain bounding box
	::Box<dim,size_t> proc_box;

	//! set of Boxes produced by the decomposition optimizer
	openfpm::vector<::Box<dim, size_t>> loc_box;

219 220 221
	/*! \brief It convert the box from the domain decomposition into sub-domain
	 *
	 * The decomposition box from the domain-decomposition contain the box in integer
222 223
	 * coordinates. This box is converted into a continuos box. It also adjust loc_box
	 * if the distribution grid and the decomposition grid are different.
224 225 226
	 *
	 * \param loc_box local box
	 *
227
	 * \return the corresponding sub-domain
228 229
	 *
	 */
230
	template<typename Memory_bx> SpaceBox<dim,T> convertDecBoxIntoSubDomain(encapc<1,::Box<dim,size_t>,Memory_bx> loc_box)
231 232 233 234 235 236 237 238
	{
		// A point with all coordinate to one
		size_t one[dim];
		for (size_t i = 0 ; i < dim ; i++)	{one[i] = 1;}

		SpaceBox<dim, size_t> sub_dc = loc_box;
		SpaceBox<dim, size_t> sub_dce = sub_dc;
		sub_dce.expand(one);
239 240 241 242 243 244 245 246 247
		sub_dce.mul(magn);

		// shrink by one
		for (size_t i = 0 ; i < dim ; i++)
		{
			loc_box.template get<Box::p1>()[i] = sub_dce.getLow(i);
			loc_box.template get<Box::p2>()[i] = sub_dce.getHigh(i) - 1;
		}

248 249 250 251 252 253 254 255 256 257 258 259 260
		SpaceBox<dim, T> sub_d(sub_dce);
		sub_d.mul(spacing);
		sub_d += domain.getP1();

		// we add the

		// Fixing sub-domains to cover all the domain

		// Fixing sub_d
		// if (loc_box) is at the boundary we have to ensure that the box span the full
		// domain (avoiding rounding off error)
		for (size_t i = 0; i < dim; i++)
		{
261
			if (sub_dc.getHigh(i) == gr.size(i) - 1)
incardon's avatar
incardon committed
262
			{sub_d.setHigh(i, domain.getHigh(i));}
263 264

			if (sub_dc.getLow(i) == 0)
incardon's avatar
incardon committed
265
			{sub_d.setLow(i,domain.getLow(i));}
266 267 268 269 270
		}

		return sub_d;
	}

271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301
	void collect_all_sub_domains(openfpm::vector<Box_map<dim,T>> & sub_domains_global)
	{
#ifdef SE_CLASS2
		check_valid(this,8);
#endif

		sub_domains_global.clear();
		openfpm::vector<Box_map<dim,T>> bm;

		for (size_t i = 0 ; i < sub_domains.size() ; i++)
		{
			Box_map<dim,T> tmp;
			tmp.box = ::SpaceBox<dim,T>(sub_domains.get(i));
			tmp.prc = v_cl.rank();

			bm.add(tmp);

		}

		v_cl.SGather(bm,sub_domains_global,0);

		size_t size = sub_domains_global.size();

		v_cl.max(size);
		v_cl.execute();

		sub_domains_global.resize(size);

		v_cl.Bcast(sub_domains_global,0);
		v_cl.execute();
	}
Pietro Incardona's avatar
Pietro Incardona committed
302 303 304

public:

305 306 307 308 309 310 311 312 313 314 315 316 317 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
	void initialize_fine_s(const ::Box<dim,T> & domain)
	{
		fine_s.clear();
		size_t div_g[dim];

		// We reduce the size of the cells by a factor 8 in 3d 4 in 2d
		for (size_t i = 0 ; i < dim ; i++)
		{div_g[i] = gr.size(i)/2;}

		fine_s.Initialize(domain,div_g);
	}

	void construct_fine_s()
	{
		collect_all_sub_domains(sub_domains_global);

		// now draw all sub-domains in fine-s

		for (size_t i = 0 ; i < sub_domains_global.size() ; i++)
		{

			// get the cells this box span
			const grid_key_dx<dim> p1 = fine_s.getCellGrid(sub_domains_global.get(i).box.getP1());
			const grid_key_dx<dim> p2 = fine_s.getCellGrid(sub_domains_global.get(i).box.getP2());

			// Get the grid and the sub-iterator
			auto & gi = fine_s.getGrid();
			grid_key_dx_iterator_sub<dim> g_sub(gi,p1,p2);

			// add the box-id to the cell list
			while (g_sub.isNext())
			{
				auto key = g_sub.get();
				fine_s.addCell(gi.LinId(key),i);
				++g_sub;
			}
		}
	}

344
	/*! \brief Constructor, it decompose and distribute the sub-domains across the processors
incardon's avatar
incardon committed
345
	 *
346
	 * \param v_cl Virtual cluster, used internally for communications
Pietro Incardona's avatar
Pietro Incardona committed
347
	 * \param bc boundary conditions
348
	 * \param opt option (one option is to construct)
349
	 *
incardon's avatar
incardon committed
350
	 */
351
	void createSubdomains(Vcluster & v_cl, const size_t (& bc)[dim], size_t opt = 0)
incardon's avatar
incardon committed
352
	{
353 354
		int p_id = v_cl.getProcessUnitID();

incardon's avatar
incardon committed
355 356 357
		// Calculate the total number of box and and the spacing
		// on each direction
		// Get the box containing the domain
358
		SpaceBox<dim, T> bs = domain.getBox();
incardon's avatar
incardon committed
359

360
		for (unsigned int i = 0; i < dim; i++)
incardon's avatar
incardon committed
361 362
		{
			// Calculate the spacing
incardon's avatar
incardon committed
363
			spacing[i] = (bs.getHigh(i) - bs.getLow(i)) / gr.size(i);
incardon's avatar
incardon committed
364 365
		}

incardon's avatar
incardon committed
366
		// fill the structure that store the processor id for each sub-domain
367
		initialize_fine_s(domain);
incardon's avatar
incardon committed
368

incardon's avatar
incardon committed
369 370
		// Optimize the decomposition creating bigger spaces
		// And reducing Ghost over-stress
371
		dec_optimizer<dim, Graph_CSR<nm_v, nm_e>> d_o(dist.getGraph(), gr_dist.getSize());
incardon's avatar
incardon committed
372

373 374 375 376 377 378 379 380 381 382
		// Ghost
		Ghost<dim,long int> ghe;

		// Set the ghost
		for (size_t i = 0 ; i < dim ; i++)
		{
			ghe.setLow(i,static_cast<long int>(ghost.getLow(i)/spacing[i]) - 1);
			ghe.setHigh(i,static_cast<long int>(ghost.getHigh(i)/spacing[i]) + 1);
		}

incardon's avatar
incardon committed
383
		// optimize the decomposition
384
		d_o.template optimize<nm_v::sub_id, nm_v::proc_id>(dist.getGraph(), p_id, loc_box, box_nn_processor,ghe,bc);
incardon's avatar
incardon committed
385

386
		// Initialize
incardon's avatar
incardon committed
387
		if (loc_box.size() > 0)
388
		{
389
			bbox = convertDecBoxIntoSubDomain(loc_box.get(0));
390
			proc_box = loc_box.get(0);
391
			sub_domains.add(bbox);
392
		}
incardon's avatar
incardon committed
393 394 395 396 397 398
		else
		{
			// invalidate all the boxes
			for (size_t i = 0 ; i < dim ; i++)
			{
				proc_box.setLow(i,0.0);
399
				proc_box.setHigh(i,0);
incardon's avatar
incardon committed
400 401

				bbox.setLow(i,0.0);
402
				bbox.setHigh(i,0);
incardon's avatar
incardon committed
403 404
			}
		}
405

incardon's avatar
incardon committed
406
		// convert into sub-domain
407
		for (size_t s = 1; s < loc_box.size(); s++)
incardon's avatar
incardon committed
408
		{
409
			SpaceBox<dim,T> sub_d = convertDecBoxIntoSubDomain(loc_box.get(s));
410

incardon's avatar
incardon committed
411 412
			// add the sub-domain
			sub_domains.add(sub_d);
incardon's avatar
incardon committed
413 414 415

			// Calculate the bound box
			bbox.enclose(sub_d);
416
			proc_box.enclose(loc_box.get(s));
incardon's avatar
incardon committed
417
		}
incardon's avatar
incardon committed
418

419
		nn_prcs<dim,T>::create(box_nn_processor, sub_domains);
420
		nn_prcs<dim,T>::applyBC(domain,ghost,bc);
421

incardon's avatar
incardon committed
422
		// fill fine_s structure
incardon's avatar
incardon committed
423 424 425
		// fine_s structure contain the processor id for each sub-sub-domain
		// with sub-sub-domain we mean the sub-domain decomposition before
		// running dec_optimizer (before merging sub-domains)
426

427
		///////////////////////////////// TODO //////////////////////////////////////////
incardon's avatar
incardon committed
428

429 430 431 432 433
		construct_fine_s();

		/////////////////////////////////////////////////////////////////////////////////

/*		grid_key_dx_iterator<dim> git(gr);
434 435 436 437 438 439 440

		while (git.isNext())
		{
			auto key = git.get();
			grid_key_dx<dim> key2;

			for (size_t i = 0 ; i < dim ; i++)
441
			{key2.set_d(i,key.get(i) / magn[i]);}
442 443 444 445

			size_t lin = gr_dist.LinId(key2);
			size_t lin2 = gr.LinId(key);

446 447
			// Here we draw the fine_s in the cell-list

448 449 450
			fine_s.get(lin2) = dist.getGraph().template vertex_p<nm_v::proc_id>(lin);

			++git;
451
		}*/
452

453 454 455 456 457 458 459 460 461 462
		Initialize_geo_cell_lists();
	}

	/*! \brief Initialize geo_cell lists
	 *
	 *
	 *
	 */
	void Initialize_geo_cell_lists()
	{
463 464 465
		// Get the processor bounding Box
		::Box<dim,T> bound = getProcessorBounds();

466 467
		// Check if the box is valid
		if (bound.isValidN() == true)
incardon's avatar
incardon committed
468 469 470
		{
			// Not necessary, but I prefer
			bound.enlarge(ghost);
471

incardon's avatar
incardon committed
472 473 474
			// calculate the sub-divisions
			size_t div[dim];
			for (size_t i = 0; i < dim; i++)
incardon's avatar
incardon committed
475
			{div[i] = (size_t) ((bound.getHigh(i) - bound.getLow(i)) / cd.getCellBox().getP2()[i]);}
476

incardon's avatar
incardon committed
477 478 479 480
			// Initialize the geo_cell structure
			ie_ghost<dim,T>::Initialize_geo_cell(bound,div);

			// Initialize shift vectors
incardon's avatar
incardon committed
481
			ie_ghost<dim,T>::generateShiftVectors(domain,bc);
incardon's avatar
incardon committed
482
		}
incardon's avatar
incardon committed
483 484
	}

485 486 487 488 489 490
	/*! \brief Calculate communication and migration costs
	 *
	 * \param ts how many timesteps have passed since last calculation, used to approximate the cost
	 */
	void computeCommunicationAndMigrationCosts(size_t ts)
	{
491
		float migration = 0;
492 493

		SpaceBox<dim, T> cellBox = cd.getCellBox();
494 495
		float b_s = static_cast<float>(cellBox.getHigh(0));
		float gh_s = static_cast<float>(ghost.getHigh(0));
496 497 498 499 500

		// compute the gh_area for 2 dim case
		float gh_v = (gh_s * b_s);

		// multiply for sub-sub-domain side for each domain
501
		for (size_t i = 2; i < dim; i++)
502 503
		{
			/* coverity[dead_error_line] */
504
			gh_v *= b_s;
505
		}
506 507 508 509 510 511 512 513 514

		size_t norm = (size_t) (1.0 / gh_v);

		migration = pow(b_s, dim);

		size_t prev = 0;

		for (size_t i = 0; i < dist.getNSubSubDomains(); i++)
		{
incardon's avatar
incardon committed
515
			dist.setMigrationCost(i, norm * migration /* * dist.getSubSubDomainComputationCost(i)*/ );
516 517 518

			for (size_t s = 0; s < dist.getNSubSubDomainNeighbors(i); s++)
			{
incardon's avatar
incardon committed
519 520 521
				// We have to remove dist.getSubSubDomainComputationCost(i) otherwise the graph is
				// not directed
				dist.setCommunicationCost(i, s, 1 /** dist.getSubSubDomainComputationCost(i)*/  *  ts);
522 523 524
			}
			prev += dist.getNSubSubDomainNeighbors(i);
		}
incardon's avatar
incardon committed
525 526

		commCostSet = true;
527
	}
incardon's avatar
incardon committed
528

529
	/*! \brief Create the sub-domain that decompose your domain
incardon's avatar
incardon committed
530 531 532 533 534
	 *
	 */
	void CreateSubspaces()
	{
		// Create a grid where each point is a space
535
		grid_sm<dim, void> g(div);
incardon's avatar
incardon committed
536 537 538 539 540 541 542 543 544 545 546

		// create a grid_key_dx iterator
		grid_key_dx_iterator<dim> gk_it(g);

		// Divide the space into subspaces
		while (gk_it.isNext())
		{
			//! iterate through all subspaces
			grid_key_dx<dim> key = gk_it.get();

			//! Create a new subspace
547
			SpaceBox<dim, T> tmp;
incardon's avatar
incardon committed
548 549

			//! fill with the Margin of the box
550
			for (int i = 0; i < dim; i++)
incardon's avatar
incardon committed
551
			{
552 553
				tmp.setHigh(i, (key.get(i) + 1) * spacing[i]);
				tmp.setLow(i, key.get(i) * spacing[i]);
incardon's avatar
incardon committed
554 555 556 557 558
			}

			//! add the space box
			sub_domains.add(tmp);

559
			// Next sub-domain
incardon's avatar
incardon committed
560 561 562 563
			++gk_it;
		}
	}

564

incardon's avatar
incardon committed
565
	/*! \brief It calculate the internal ghost boxes
566 567 568 569 570 571 572 573
	 *
	 * Example: Processor 10 calculate
	 * B8_0 B9_0 B9_1 and B5_0
	 *
	 *
	 *
	 \verbatim

Pietro Incardona's avatar
Pietro Incardona committed
574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597
	+----------------------------------------------------+
	|                                                    |
	|                 Processor 8                        |
	|                 Sub+domain 0                       +-----------------------------------+
	|                                                    |                                   |
	|                                                    |                                   |
	++--------------+---+---------------------------+----+        Processor 9                |
	 |              |   |     B8_0                  |    |        Subdomain 0                |
	 |              +------------------------------------+                                   |
	 |              |   |                           |    |                                   |
	 |              |   |                           |B9_0|                                   |
	 |              | B |    Local processor        |    |                                   |
	 | Processor 5  | 5 |    Subdomain 0            |    |                                   |
	 | Subdomain 0  | _ |                           +----------------------------------------+
	 |              | 0 |                           |    |                                   |
	 |              |   |                           |    |                                   |
	 |              |   |                           |    |        Processor 9                |
	 |              |   |                           |B9_1|        Subdomain 1                |
	 |              |   |                           |    |                                   |
	 |              |   |                           |    |                                   |
	 |              |   |                           |    |                                   |
	 +--------------+---+---------------------------+----+                                   |
														 |                                   |
														 +-----------------------------------+
598 599 600 601 602 603 604


 \endverbatim

       and also
       G8_0 G9_0 G9_1 G5_0 (External ghost boxes)

Pietro Incardona's avatar
Pietro Incardona committed
605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630
\verbatim

		  +----------------------------------------------------+
		  |                 Processor 8                        |
		  |                 Subdomain 0                        +-----------------------------------+
		  |                                                    |                                   |
		  |           +---------------------------------------------+                              |
		  |           |         G8_0                           |    |                              |
	+-----+---------------+------------------------------------+    |   Processor 9                |
	|                 |   |                                    |    |   Subdomain 0                |
	|                 |   |                                    |G9_0|                              |
	|                 |   |                                    |    |                              |
	|                 |   |                                    |    |                              |
	|                 |   |        Local processor             |    |                              |
	|  Processor 5    |   |        Sub+domain 0                |    |                              |
	|  Subdomain 0    |   |                                    +-----------------------------------+
	|                 |   |                                    |    |                              |
	|                 | G |                                    |    |                              |
	|                 | 5 |                                    |    |   Processor 9                |
	|                 | | |                                    |    |   Subdomain 1                |
	|                 | 0 |                                    |G9_1|                              |
	|                 |   |                                    |    |                              |
	|                 |   |                                    |    |                              |
	+---------------------+------------------------------------+    |                              |
					  |                                        |    |                              |
					  +----------------------------------------+----+------------------------------+
631 632 633 634

	 \endverbatim

	 *
incardon's avatar
incardon committed
635
	 * ghost margins for each dimensions (p1 negative part) (p2 positive part)
636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668
	 *
	 *
	 \verbatim

	 	 	 	 	 ^ p2[1]
	 	 	 	 	 |
	 	 	 	 	 |
	 	 	 	+----+----+
	 	 	 	|         |
	 	 	 	|         |
	 p1[0]<-----+         +----> p2[0]
	 	 	 	|         |
	 	 	 	|         |
	 	 	 	+----+----+
	 	 	 	 	 |
	 	 	 	 	 v  p1[1]

	 \endverbatim

	 *
	 *
	 */
	void calculateGhostBoxes()
	{
		// Intersect all the local sub-domains with the sub-domains of the contiguous processors

		// create the internal structures that store ghost information
		ie_ghost<dim, T>::create_box_nn_processor_ext(v_cl, ghost, sub_domains, box_nn_processor, *this);
		ie_ghost<dim, T>::create_box_nn_processor_int(v_cl, ghost, sub_domains, box_nn_processor, *this);

		ie_loc_ghost<dim,T>::create(sub_domains,domain,ghost,bc);
	}

669 670 671 672
	template<typename T2> inline size_t processorID_impl(T2 & p) const
	{
		// Get the number of elements in the cell

incardon's avatar
incardon committed
673
		size_t e = -1;
674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700
		size_t cl = fine_s.getCell(p);
		size_t n_ele = fine_s.getNelements(cl);

		for (size_t i = 0 ; i < n_ele ; i++)
		{
			e = fine_s.get(cl,i);

			if (sub_domains_global.get(e).box.isInsideNP(p) == true)
			{
				break;
			}
		}

#ifdef SE_CLASS1

		if (n_ele == 0)
		{
			std::cout << __FILE__ << ":" << __LINE__ << " I cannot detect in which processor this particle go" << std::endl;
			return -1;
		}

#endif

		return sub_domains_global.get(e).prc;
	}


incardon's avatar
incardon committed
701 702
public:

Pietro Incardona's avatar
Pietro Incardona committed
703
	//! Space dimensions
704 705
	static constexpr int dims = dim;

Pietro Incardona's avatar
Pietro Incardona committed
706
	//! Space type
707 708
	typedef T stype;

incardon's avatar
incardon committed
709 710 711 712 713 714 715 716 717 718 719 720 721 722
	//! Increment the reference counter
	void incRef()
	{ref_cnt++;}

	//! Decrement the reference counter
	void decRef()
	{ref_cnt--;}

	//! Return the reference counter
	long int ref()
	{
		return ref_cnt;
	}

incardon's avatar
incardon committed
723 724
	/*! \brief Cartesian decomposition constructor
	 *
725
	 * \param v_cl Virtual cluster, used internally to handle or pipeline communication
incardon's avatar
incardon committed
726 727
	 *
	 */
728 729
	CartDecomposition(Vcluster & v_cl)
	:nn_prcs<dim, T>(v_cl), v_cl(v_cl), dist(v_cl),ref_cnt(0)
incardon's avatar
incardon committed
730 731 732 733
	{
		// Reset the box to zero
		bbox.zero();
	}
incardon's avatar
incardon committed
734

735 736 737 738 739
	/*! \brief Cartesian decomposition copy constructor
	 *
     * \param cart object to copy
	 *
	 */
incardon's avatar
incardon committed
740
	CartDecomposition(const CartDecomposition<dim,T,Memory,Distribution> & cart)
741
	:nn_prcs<dim,T>(cart.v_cl),v_cl(cart.v_cl),dist(v_cl),ref_cnt(0)
742 743 744 745 746 747 748 749 750
	{
		this->operator=(cart);
	}

	/*! \brief Cartesian decomposition copy constructor
	 *
     * \param cart object to copy
	 *
	 */
incardon's avatar
incardon committed
751
	CartDecomposition(CartDecomposition<dim,T,Memory,Distribution> && cart)
752
	:nn_prcs<dim,T>(cart.v_cl),v_cl(cart.v_cl),dist(v_cl),ref_cnt(0)
753 754 755 756
	{
		this->operator=(cart);
	}

incardon's avatar
incardon committed
757 758
	//! Cartesian decomposition destructor
	~CartDecomposition()
759 760
	{
	}
incardon's avatar
incardon committed
761

762 763 764 765 766 767 768 769 770 771 772 773 774 775
	/*! \brief class to select the returned id by ghost_processorID
	 *
	 */
	class box_id
	{
	public:
		/*! \brief Return the box id
		 *
		 * \param p structure containing the id informations
		 * \param b_id box_id
		 *
		 * \return box id
		 *
		 */
776
		inline static size_t id(p_box<dim, T> & p, size_t b_id)
777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795
		{
			return b_id;
		}
	};

	/*! \brief class to select the returned id by ghost_processorID
	 *
	 */
	class processor_id
	{
	public:
		/*! \brief Return the processor id
		 *
		 * \param p structure containing the id informations
		 * \param b_id box_id
		 *
		 * \return processor id
		 *
		 */
796
		inline static size_t id(p_box<dim, T> & p, size_t b_id)
797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815
		{
			return p.proc;
		}
	};

	/*! \brief class to select the returned id by ghost_processorID
	 *
	 */
	class lc_processor_id
	{
	public:
		/*! \brief Return the near processor id
		 *
		 * \param p structure containing the id informations
		 * \param b_id box_id
		 *
		 * \return local processor id
		 *
		 */
816
		inline static size_t id(p_box<dim, T> & p, size_t b_id)
817 818 819 820 821
		{
			return p.lc_proc;
		}
	};

822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841
	/*! \brief class to select the returned id by ghost_processorID
	 *
	 */
	class shift_id
	{
	public:
		/*! \brief Return the shift id
		 *
		 * \param p structure containing the id informations
		 * \param b_id box_id
		 *
		 * \return shift_id id
		 *
		 */
		inline static size_t id(p_box<dim,T> & p, size_t b_id)
		{
			return p.shift_id;
		}
	};

842 843
	/*! \brief Apply boundary condition to the point
	 *
Pietro Incardona's avatar
Pietro Incardona committed
844 845 846 847 848
	 * If the particle go out to the right, bring back the particle on the left
	 * in case of periodic, nothing in case of non periodic
	 *
	 * \param pt Point to apply the boundary condition. (it's coordinated are changed according the
	 *        the explanation before)
849 850 851 852 853 854 855
	 *
	 */
	void applyPointBC(float (& pt)[dim]) const
	{
		for (size_t i = 0 ; i < dim ; i++)
		{
			if (bc[i] == PERIODIC)
856
			{pt[i] = openfpm::math::periodic_l(pt[i],domain.getHigh(i),domain.getLow(i));}
857 858 859 860 861
		}
	}

	/*! \brief Apply boundary condition to the point
	 *
Pietro Incardona's avatar
Pietro Incardona committed
862 863 864 865 866
	 * If the particle go out to the right, bring back the particle on the left
	 * in case of periodic, nothing in case of non periodic
	 *
	 * \param pt Point to apply the boundary conditions.(it's coordinated are changed according the
	 *        the explanation before)
867 868 869 870 871 872 873
	 *
	 */
	void applyPointBC(Point<dim,T> & pt) const
	{
		for (size_t i = 0 ; i < dim ; i++)
		{
			if (bc[i] == PERIODIC)
874
			{pt.get(i) = openfpm::math::periodic_l(pt.get(i),domain.getHigh(i),domain.getLow(i));}
875 876 877 878 879
		}
	}

	/*! \brief Apply boundary condition to the point
	 *
Pietro Incardona's avatar
Pietro Incardona committed
880 881 882 883 884
	 * If the particle go out to the right, bring back the particle on the left
	 * in case of periodic, nothing in case of non periodic
	 *
	 * \param pt encapsulated point object (it's coordinated are changed according the
	 *        the explanation before)
885 886 887 888 889 890 891
	 *
	 */
	template<typename Mem> void applyPointBC(encapc<1,Point<dim,T>,Mem> && pt) const
	{
		for (size_t i = 0 ; i < dim ; i++)
		{
			if (bc[i] == PERIODIC)
892
			{pt.template get<0>()[i] = openfpm::math::periodic_l(pt.template get<0>()[i],domain.getHigh(i),domain.getLow(i));}
893 894 895
		}
	}

896 897 898 899 900 901 902
	/*! \brief It create another object that contain the same decomposition information but with different ghost boxes
	 *
	 * \param g ghost
	 *
	 * \return a duplicated decomposition with different ghost boxes
	 *
	 */
incardon's avatar
incardon committed
903
	CartDecomposition<dim,T,Memory,Distribution> duplicate(const Ghost<dim,T> & g) const
904
	{
incardon's avatar
incardon committed
905
		CartDecomposition<dim,T,Memory,Distribution> cart(v_cl);
906 907 908 909 910 911 912 913

		cart.box_nn_processor = box_nn_processor;
		cart.sub_domains = sub_domains;
		cart.fine_s = fine_s;

		cart.gr = gr;
		cart.cd = cd;
		cart.domain = domain;
914 915
		for (size_t i = 0 ; i < dim ; i++)
		{cart.spacing[i] = spacing[i];};
916 917 918 919

		cart.bbox = bbox;
		cart.ghost = g;

920 921
		cart.dist = dist;

922 923 924
		for (size_t i = 0 ; i < dim ; i++)
			cart.bc[i] = bc[i];

925 926
		(static_cast<nn_prcs<dim,T> &>(cart)).create(box_nn_processor, sub_domains);
		(static_cast<nn_prcs<dim,T> &>(cart)).applyBC(domain,ghost,bc);
927

928 929
		cart.Initialize_geo_cell_lists();
		cart.calculateGhostBoxes();
930 931 932 933 934

		return cart;
	}

	/*! \brief It create another object that contain the same information and act in the same way
935
	 *
Pietro Incardona's avatar
Pietro Incardona committed
936
	 * \return a duplicated CartDecomposition object
937 938
	 *
	 */
incardon's avatar
incardon committed
939
	CartDecomposition<dim,T,Memory,Distribution> duplicate() const
940
	{
incardon's avatar
incardon committed
941
		CartDecomposition<dim,T,Memory,Distribution> cart(v_cl);
942 943 944 945 946 947 948 949 950

		(static_cast<ie_loc_ghost<dim,T>*>(&cart))->operator=(static_cast<ie_loc_ghost<dim,T>>(*this));
		(static_cast<nn_prcs<dim,T>*>(&cart))->operator=(static_cast<nn_prcs<dim,T>>(*this));
		(static_cast<ie_ghost<dim,T>*>(&cart))->operator=(static_cast<ie_ghost<dim,T>>(*this));

		cart.sub_domains = sub_domains;
		cart.box_nn_processor = box_nn_processor;
		cart.fine_s = fine_s;
		cart.gr = gr;
incardon's avatar
incardon committed
951 952 953
		cart.gr_dist = gr_dist;
		cart.dist = dist;
		cart.commCostSet = commCostSet;
954 955
		cart.cd = cd;
		cart.domain = domain;
incardon's avatar
incardon committed
956
		cart.sub_domains_global = sub_domains_global;
957 958
		for (size_t i = 0 ; i < dim ; i++)
		{cart.spacing[i] = spacing[i];};
959 960 961

		cart.ghost = ghost;

962 963
		cart.bbox = bbox;

964 965 966
		for (size_t i = 0 ; i < dim ; i++)
			cart.bc[i] = this->bc[i];

967 968 969 970 971 972 973
		return cart;
	}

	/*! \brief Copy the element
	 *
	 * \param cart element to copy
	 *
Pietro Incardona's avatar
Pietro Incardona committed
974 975
	 * \return itself
	 *
976
	 */
977
	CartDecomposition<dim,T,Memory, Distribution> & operator=(const CartDecomposition & cart)
978 979 980 981 982 983 984 985 986
	{
		static_cast<ie_loc_ghost<dim,T>*>(this)->operator=(static_cast<ie_loc_ghost<dim,T>>(cart));
		static_cast<nn_prcs<dim,T>*>(this)->operator=(static_cast<nn_prcs<dim,T>>(cart));
		static_cast<ie_ghost<dim,T>*>(this)->operator=(static_cast<ie_ghost<dim,T>>(cart));

		sub_domains = cart.sub_domains;
		box_nn_processor = cart.box_nn_processor;
		fine_s = cart.fine_s;
		gr = cart.gr;
incardon's avatar
incardon committed
987 988 989
		gr_dist = cart.gr_dist;
		dist = cart.dist;
		commCostSet = cart.commCostSet;
990 991
		cd = cart.cd;
		domain = cart.domain;
incardon's avatar
incardon committed
992
		sub_domains_global = cart.sub_domains_global;
993 994

		for (size_t i = 0 ; i < dim ; i++)
incardon's avatar
incardon committed
995 996 997 998
		{
			spacing[i] = cart.spacing[i];
			magn[i] = cart.magn[i];
		};
999 1000 1001

		ghost = cart.ghost;

1002 1003
		bbox = cart.bbox;

1004 1005 1006
		for (size_t i = 0 ; i < dim ; i++)
			bc[i] = cart.bc[i];

1007 1008 1009 1010 1011 1012 1013
		return *this;
	}

	/*! \brief Copy the element, move semantic
	 *
	 * \param cart element to copy
	 *
Pietro Incardona's avatar
Pietro Incardona committed
1014 1015
	 * \return itself
	 *
1016
	 */
incardon's avatar
incardon committed
1017
	CartDecomposition<dim,T,Memory,Distribution> & operator=(CartDecomposition && cart)
1018
	{
Pietro Incardona's avatar
Pietro Incardona committed
1019 1020 1021
		static_cast<ie_loc_ghost<dim,T>*>(this)->operator=(static_cast<ie_loc_ghost<dim,T>>(cart));
		static_cast<nn_prcs<dim,T>*>(this)->operator=(static_cast<nn_prcs<dim,T>>(cart));
		static_cast<ie_ghost<dim,T>*>(this)->operator=(static_cast<ie_ghost<dim,T>>(cart));
1022 1023 1024 1025 1026

		sub_domains.swap(cart.sub_domains);
		box_nn_processor.swap(cart.box_nn_processor);
		fine_s.swap(cart.fine_s);
		gr = cart.gr;
incardon's avatar
incardon committed
1027 1028 1029
		gr_dist = cart.gr_dist;
		dist = cart.dist;
		commCostSet = cart.commCostSet;
1030 1031
		cd = cart.cd;
		domain = cart.domain;
incardon's avatar
incardon committed
1032 1033
		sub_domains_global.swap(cart.sub_domains_global);

1034
		for (size_t i = 0 ; i < dim ; i++)
incardon's avatar
incardon committed
1035 1036 1037 1038
		{
			spacing[i] = cart.spacing[i];
			magn[i] = cart.magn[i];
		};
1039 1040 1041

		ghost = cart.ghost;

1042
		bbox = cart.bbox;
1043

1044
		for (size_t i = 0 ; i < dim ; i++)
1045 1046 1047
			bc[i] = cart.bc[i];

		return *this;
1048 1049
	}

1050 1051 1052 1053 1054 1055
	/*! \brief The default grid size
	 *
	 *  The default grid is always an isotropic grid that adapt with the number of processors,
	 *  it define in how many cell it will be divided the space for a particular required minimum
	 *  number of sub-domain
	 *
Pietro Incardona's avatar
Pietro Incardona committed
1056 1057 1058 1059
	 * \param n_sub number of subdomains per processors
	 *
	 * \return grid dimension (it is one number because on the other dimensions is the same)
	 *
1060 1061 1062 1063 1064
	 */
	static size_t getDefaultGrid(size_t n_sub)
	{
		// Calculate the number of sub-sub-domain on
		// each dimension
1065
		return openfpm::math::round_big_2(pow(n_sub, 1.0 / dim));
1066 1067
	}

1068
	/*! \brief Given a point return in which processor the particle should go
Pietro Incardona's avatar
Pietro Incardona committed
1069 1070
	 *
	 * \param p point
incardon's avatar
incardon committed
1071 1072 1073 1074
	 *
	 * \return processorID
	 *
	 */
incardon's avatar
incardon committed
1075
	template<typename Mem> size_t inline processorID(const encapc<1, Point<dim,T>, Mem> & p) const
incardon's avatar
incardon committed
1076
	{
1077
		return processorID_impl(p);
incardon's avatar
incardon committed
1078 1079
	}

1080
	/*! \brief Given a point return in which processor the particle should go
Pietro Incardona's avatar
Pietro Incardona committed
1081 1082
	 *
	 * \param p point
incardon's avatar
incardon committed
1083
	 *
1084
	 * \return processorID
incardon's avatar
incardon committed
1085 1086
	 *
	 */
1087
	size_t inline processorID(const Point<dim,T> &p) const
incardon's avatar
incardon committed
1088
	{
1089
		return processorID_impl(p);
incardon's avatar
incardon committed
1090 1091
	}

1092
	/*! \brief Given a point return in which processor the particle should go
Pietro Incardona's avatar
Pietro Incardona committed
1093 1094
	 *
	 * \param p point
incardon's avatar
incardon committed
1095 1096 1097 1098
	 *
	 * \return processorID
	 *
	 */
1099
	size_t inline processorID(const T (&p)[dim]) const
incardon's avatar
incardon committed
1100
	{
1101
		return processorID_impl(p);
incardon's avatar
incardon committed
1102
	}
incardon's avatar
incardon committed
1103

Pietro Incardona's avatar
Pietro Incardona committed
1104
	/*! \brief Given a point return in which processor the point/particle should go
incardon's avatar
incardon committed
1105
	 *
1106 1107
	 * Boundary conditions are considered
	 *
Pietro Incardona's avatar
Pietro Incardona committed
1108 1109
	 * \param p point
	 *
1110
	 * \return processorID
incardon's avatar
incardon committed
1111 1112
	 *
	 */
1113
	template<typename Mem> size_t inline processorIDBC(encapc<1, Point<dim,T>, Mem> p)
incardon's avatar
incardon committed
1114
	{
1115 1116 1117
		Point<dim,T> pt = p;
		applyPointBC(pt);

1118

1119
		return processorID_impl(pt);
incardon's avatar
incardon committed
1120 1121
	}

1122
	/*! \brief Given a point return in which processor the particle should go
1123 1124
	 *
	 * Boundary conditions are considered
incardon's avatar
incardon committed
1125
	 *
Pietro Incardona's avatar
Pietro Incardona committed
1126 1127
	 * \param p point
	 *
incardon's avatar
incardon committed
1128 1129 1130
	 * \return processorID
	 *
	 */
1131
	size_t inline processorIDBC(const Point<dim,T> &p) const
1132 1133 1134
	{
		Point<dim,T> pt = p;
		applyPointBC(pt);
incardon's avatar
incardon committed
1135

1136 1137
		// Get the number of elements in the cell

1138
		return processorID_impl(pt);
1139 1140 1141 1142 1143 1144
	}

	/*! \brief Given a point return in which processor the particle should go
	 *
	 * Boundary consition are considered
	 *
Pietro Incardona's avatar
Pietro Incardona committed
1145 1146
	 * \param p point position
	 *
1147 1148 1149
	 * \return processorID
	 *
	 */
1150
	size_t inline processorIDBC(const T (&p)[dim]) const
incardon's avatar
incardon committed
1151
	{
1152 1153 1154
		Point<dim,T> pt = p;
		applyPointBC(pt);

1155
		return processorID_impl(pt);
incardon's avatar
incardon committed
1156 1157
	}

1158 1159 1160 1161 1162 1163 1164
	/*! \brief Get the periodicity on i dimension
	 *
	 * \param i dimension
	 *
	 * \return the periodicity in direction i
	 *
	 */
1165
	inline size_t periodicity(size_t i) const
1166 1167 1168 1169
	{
		return bc[i];
	}

1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180
	/*! \brief Get the periodicity
	 *
	 *
	 * \return the periodicity
	 *
	 */
	inline const size_t (& periodicity() const) [dim]
	{
		return bc;
	}

1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197
	/*! \brief Calculate magnification
	 *
	 * \param gm distribution grid
	 *
	 */
	void calculate_magn(const grid_sm<dim,void> & gm)
	{
		if (gm.size() == 0)
		{
			for (size_t i = 0 ; i < dim ; i++)
				magn[i] = 1;
		}
		else
		{
			for (size_t i = 0 ; i < dim ; i++)
			{
				if (gr.size(i) % gm.size(i) != 0)
incardon's avatar
incardon committed
1198
					std::cerr << __FILE__ << ":" << __LINE__ << ".Error the decomposition grid specified as gr.size(" << i << ")=" << gr.size(i) << " is not multiple of the distribution grid gm.size(" << i << ")=" << gm.size(i) << std::endl;
1199 1200 1201 1202 1203 1204

				magn[i] = gr.size(i) / gm.size(i);
			}
		}
	}

incardon's avatar
incardon committed
1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275
	/*! \brief Set the best parameters for the decomposition
	 *
	 * It based on number of processors and dimensionality find a "good" parameter setting
	 *
	 * \param domain_ domain to decompose
	 * \param bc boundary conditions
	 * \param ghost Ghost size
	 * \param sec_dist Distribution grid. The distribution grid help in reducing the underlying
	 *                 distribution problem simplifying decomposition problem. This is done in order to
	 *                 reduce the load/balancing dynamic load balancing problem
	 *
	 * \param dec_gran number of sub-sub-domain for each processor
	 *
	 */
	void setGoodParameters(::Box<dim,T> domain_,
						   const size_t (& bc)[dim],
						   const Ghost<dim,T> & ghost,
						   size_t dec_gran,
						   const grid_sm<dim,void> & sec_dist = grid_sm<dim,void>())
	{
		size_t div[dim];

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

		// Calculate the maximum number (before merging) of sub-domain on
		// each dimension

		nsub_to_div2(div,n_sub,dim);

/*		for (size_t i = 0; i < dim; i++)
		{
			div[i] = openfpm::math::round_big_2(pow(n_sub, 1.0 / dim));
		}*/

		if (dim > 3)
		{
			long int dim_r = dim-1;
			do
			{
				// Check for adjustment
				size_t tot_size = 1;
				for (size_t i = 0 ; i < dim ; i++)
				{tot_size *= div[i];}

				// the granularity is too coarse increase the divisions
				if (tot_size / n_proc > 0.75*dec_gran )
				{break;}

				nsub_to_div(div,n_sub,dim_r);

				dim_r--;
			} while(dim_r > 0);
		}

		setParameters(div,domain_,bc,ghost,sec_dist);
	}

	/*! \brief return the parameters of the decomposition
	 *
	 * \param div_ number of divisions in each dimension
	 *
	 */
	void getParameters(size_t (& div_)[dim])
	{
		for (size_t i = 0 ; i < dim ; i++)
		{div_[i] = this->gr.size(i);}
	}
1276

incardon's avatar
incardon committed
1277 1278
	/*! \brief Set the parameter of the decomposition
	 *
1279
	 * \param div_ storing into how many sub-sub-domains to decompose on each dimension
1280
	 * \param domain_ domain to decompose
Pietro Incardona's avatar
Pietro Incardona committed
1281
	 * \param bc boundary conditions
1282
	 * \param ghost Ghost size
1283 1284 1285
	 * \param sec_dist Distribution grid. The distribution grid help in reducing the underlying
	 *                 distribution problem simplifying decomposition problem. This is done in order to
	 *                 reduce the load/balancing dynamic load balancing problem
incardon's avatar
incardon committed
1286 1287
	 *
	 */
incardon's avatar
incardon committed
1288 1289 1290 1291 1292
	void setParameters(const size_t (& div_)[dim],
					   ::Box<dim,T> domain_,
						const size_t (& bc)[dim],
						const Ghost<dim,T> & ghost,
						const grid_sm<dim,void> & sec_dist = grid_sm<dim,void>())
incardon's avatar
incardon committed
1293
	{
1294 1295 1296 1297
		// set the boundary conditions
		for (size_t i = 0 ; i < dim ; i++)
			this->bc[i] = bc[i];

incardon's avatar
incardon committed
1298 1299
		// set the ghost
		this->ghost = ghost;
incardon's avatar
incardon committed
1300

1301
		// Set the decomposition parameters
incardon's avatar
incardon committed
1302
		gr.setDimensions(div_);
incardon's avatar
incardon committed
1303
		domain = domain_;
1304
		cd.setDimensions(domain, div_, 0);
incardon's avatar
incardon committed
1305

1306 1307 1308 1309 1310 1311 1312 1313 1314 1315 1316 1317
		// We we have a secondary grid costruct a reduced graph
		if (sec_dist.size(0) != 0)
		{
			calculate_magn(sec_dist);
			gr_dist.setDimensions(sec_dist.getSize());
		}
		else
		{
			calculate_magn(sec_dist);
			gr_dist = gr;
		}

1318
		// init distribution
1319
		dist.createCartGraph(gr_dist, domain);
1320 1321 1322

	}

Pietro Incardona's avatar
Pietro Incardona committed
1323 1324 1325 1326
	/*! \brief Delete the decomposition and reset the data-structure
	 *
	 *
	 */
1327 1328 1329 1330 1331
	void reset()
	{
		sub_domains.clear();
		box_nn_processor.clear();
		fine_s.clear();
incardon's avatar
incardon committed
1332
		loc_box.clear();
1333 1334 1335 1336 1337
		nn_prcs<dim, T>::reset();
		ie_ghost<dim, T>::reset();
		ie_loc_ghost<dim, T>::reset();
	}

1338 1339 1340 1341 1342
	/*! \brief Start decomposition
	 *
	 */
	void decompose()
	{
1343 1344
		reset();

incardon's avatar
incardon committed
1345
		if (commCostSet == false)
incardon's avatar
incardon committed
1346
		{computeCommunicationAndMigrationCosts(1);}
1347 1348

		dist.decompose();
incardon's avatar
incardon committed
1349

1350
		createSubdomains(v_cl,bc);
1351 1352

		calculateGhostBoxes();
incardon's avatar
incardon committed
1353 1354

		domain_nn_calculator_cart<dim>::reset();
1355
		domain_nn_calculator_cart<dim>::setParameters(proc_box);
incardon's avatar
incardon committed
1356 1357
	}

1358
	/*! \brief Refine the decomposition, available only for ParMetis distribution, for Metis it is a null call
Pietro Incardona's avatar
Pietro Incardona committed
1359 1360
	 *
	 * \param ts number of time step from the previous load balancing
1361 1362
	 *
	 */
incardon's avatar
incardon committed
1363
	void refine(size_t ts)
1364
	{
1365 1366
		reset();

incardon's avatar
incardon committed
1367
		if (commCostSet == false)
incardon's avatar
incardon committed
1368
		{computeCommunicationAndMigrationCosts(ts);}
1369 1370

		dist.refine();
1371 1372 1373 1374

		createSubdomains(v_cl,bc);

		calculateGhostBoxes();
incardon's avatar
incardon committed
1375 1376

		domain_nn_calculator_cart<dim>::reset();
1377
		domain_nn_calculator_cart<dim>::setParameters(proc_box);
1378 1379
	}

incardon's avatar
incardon committed
1380 1381 1382 1383 1384 1385 1386 1387 1388 1389
	/*! \brief Refine the decomposition, available only for ParMetis distribution, for Metis it is a null call
	 *
	 * \param ts number of time step from the previous load balancing
	 *
	 */
	void redecompose(size_t ts)
	{
		reset();

		if (commCostSet == false)
incardon's avatar
incardon committed
1390
		{computeCommunicationAndMigrationCosts(ts);}
incardon's avatar
incardon committed
1391 1392 1393 1394 1395 1396 1397 1398

		dist.redecompose();

		createSubdomains(v_cl,bc);

		calculateGhostBoxes();

		domain_nn_calculator_cart<dim>::reset();
1399
		domain_nn_calculator_cart<dim>::setParameters(proc_box);
incardon's avatar
incardon committed
1400 1401
	}

1402
	/*! \brief Refine the decomposition, available only for ParMetis distribution, for Metis it is a null call
Pietro Incardona's avatar
Pietro Incardona committed
1403 1404
	 *
	 * \param dlb Dynamic load balancing object
1405
	 *
1406
	 * \return true if the re-balance has been executed, false otherwise
1407
	 */
incardon's avatar
incardon committed
1408
	bool refine(DLB & dlb)
1409
	{
1410 1411 1412 1413 1414 1415 1416
		// if the DLB heuristic to use is the "Unbalance Threshold" get unbalance percentage
		if (dlb.getHeurisitc() == DLB::Heuristic::UNBALANCE_THRLD)
		{
			float unbalance = dist.getUnbalance();
			dlb.setUnbalance(unbalance);
			if (v_cl.getProcessUnitID() == 0)
			{
1417
				std::cout << std::setprecision(3) << unbalance << "\n";
1418 1419 1420
			}
		}

1421 1422
		if (dlb.rebalanceNeeded())
		{
incardon's avatar
incardon committed
1423
			refine(dlb.getNTimeStepSinceDLB());
1424

1425
			return true;
1426
		}
1427 1428 1429
		return false;
	}

Pietro Incardona's avatar
Pietro Incardona committed
1430 1431
//	size_t n_step = 0;

1432 1433 1434 1435 1436 1437 1438 1439 1440 1441 1442 1443 1444 1445 1446 1447
	/*! \brief Get the current un-balance value
	 *
	 * \return the un-balance percentage value
	 */
	float getUnbalance()
	{
		return dist.getUnbalance();
	}

	/*! \brief Compute the processor load counting the total weights of its vertices
	 *
	 * \return the current processor load
	 */
	size_t getProcessorLoad()
	{
		return dist.getProcessorLoad();
1448 1449 1450 1451 1452 1453 1454 1455
	}

	/*! \brief function that return the position of the cell in the space
	 *
	 * \param id vertex id
	 * \param pos vector that will contain x, y, z
	 *
	 */
1456
	inline void getSubSubDomainPosition(size_t id, T (&pos)[dim])
1457
	{
1458
		dist.getSubSubDomainPosition(id, pos);
1459 1460
	}

1461
	//TODO fix in Parmetis distribution to get only the right amount of vertices
1462 1463
	/*! \brief Get the number of sub-sub-domains in this sub-graph
	 *
1464
	 * \return number of sub-sub-domains in this sub-graph
1465 1466 1467 1468 1469 1470
	 */
	size_t getNSubSubDomains()
	{
		return dist.getNSubSubDomains();
	}

Pietro Incardona's avatar
Pietro Incardona committed
1471
	/*! \brief Function that set the computational cost for a of a sub-sub domain
1472 1473
	 *
	 * \param id vertex id
Pietro Incardona's avatar
Pietro Incardona committed
1474
	 * \param weight compotational cost
1475 1476 1477 1478
	 *
	 */
	inline void setSubSubDomainComputationCost(size_t id, size_t weight)
	{
1479
		dist.setComputationCost(id, weight);
1480 1481
	}

Pietro Incardona's avatar
Pietro Incardona committed
1482
	/*! \brief function that return the computation cost of the sub-sub-domain id
1483
	 *
Pietro Incardona's avatar
Pietro Incardona committed
1484 1485 1486
	 * \param id sub-sub-domain id
	 *
	 * \return the computational cost
1487 1488 1489 1490
	 *
	 */
	inline size_t getSubSubDomainComputationCost(size_t id)
	{
incardon's avatar
incardon committed
1491
		return dist.getSubSubDomainComputationCost(id);
1492 1493 1494 1495 1496 1497 1498 1499 1500 1501 1502
	}

	/*! \brief Operator to access the size of the sub-graph
	 *
	 * \return the size of the subgraph
	 */
	size_t subSize()
	{
		return dist.subSize();
	}

1503
	/*! \brief Get the number of local sub-domains
incardon's avatar