dec_optimizer.hpp 20 KB
Newer Older
incardon's avatar
incardon committed
1 2 3
#ifndef DEC_OPTIMIZER_HPP
#define DEC_OPTIMIZER_HPP

Pietro Incardona's avatar
Pietro Incardona committed
4
#include "Grid/iterators/grid_key_dx_iterator_sub.hpp"
5
#include "Grid/iterators/grid_skin_iterator.hpp"
Pietro Incardona's avatar
Pietro Incardona committed
6

incardon's avatar
incardon committed
7 8
/*! \brief this class represent a wavefront of dimension dim
 *
Pietro Incardona's avatar
Pietro Incardona committed
9
 * \tparam dim Dimensionality of the wavefront (dimensionality of the space
incardon's avatar
incardon committed
10 11 12
 *                                       where it live so the wavefront
 *                                       is dim-1)
 *
Pietro Incardona's avatar
Pietro Incardona committed
13 14 15
 * Each wavefront is identified by one starting point and one stop point.
 * More or less a wavefront is just a box defined in the integer space
 *
incardon's avatar
incardon committed
16 17
 */
template <unsigned int dim>
Pietro Incardona's avatar
Pietro Incardona committed
18
class wavefront : public Box<dim,size_t>
incardon's avatar
incardon committed
19 20 21
{
public:

Pietro Incardona's avatar
Pietro Incardona committed
22
	//! start point is the property with id 0 (first property)
incardon's avatar
incardon committed
23 24
	static const int start = 0;

Pietro Incardona's avatar
Pietro Incardona committed
25 26
	//! stop point is the property with id 1 (second property)
	static const int stop = 1;
incardon's avatar
incardon committed
27 28 29 30 31
};

/*! \brief This class take a graph representing the space decomposition and produce a
 *         simplified version
 *
incardon's avatar
incardon committed
32 33 34
 * Given a Graph_CSR and a seed point produce an alternative decomposition in boxes with
 * less sub-domain. In the following we referee with sub-domain the boxes produced by this
 * algorithm and sub-sub-domain the sub-domain before reduction
incardon's avatar
incardon committed
35 36 37 38 39 40
 *
 */

template <unsigned int dim, typename Graph>
class dec_optimizer
{
Pietro Incardona's avatar
Pietro Incardona committed
41
	//! Contain information about the grid size
42
	grid_sm<dim,void> gh;
incardon's avatar
incardon committed
43 44 45

private:

46 47 48 49 50 51 52 53 54
	/*! \brief Expand one wavefront
	 *
	 * \param v_w wavefronts
	 * \param w_comb wavefront expansion combinations
	 * \param d direction of expansion
	 *
	 */
	void expand_one_wf(openfpm::vector<wavefront<dim>> & v_w, std::vector<comb<dim>> & w_comb , size_t d)
	{
55
		for (size_t j = 0 ; j < dim ; j++)
56 57 58 59 60 61 62 63 64
		{
			v_w.template get<wavefront<dim>::stop>(d)[j] = v_w.template get<wavefront<dim>::stop>(d)[j] + w_comb[d].c[j];
			v_w.template get<wavefront<dim>::start>(d)[j] = v_w.template get<wavefront<dim>::start>(d)[j] + w_comb[d].c[j];
		}
	}


	/*! \brief Adjust the other wavefronts
	 *
Pietro Incardona's avatar
Pietro Incardona committed
65 66 67
	 * \param v_w array of wavefronts
	 * \param hyp Hyper cube used to adjust the wavefront
	 * \param w_comb for each wavefront indicate their position (normal to the face of the wavefront)
68 69 70 71 72 73 74 75 76 77 78
	 * \param d direction
	 *
	 */
	void adjust_others_wf(openfpm::vector<wavefront<dim>> & v_w,  HyperCube<dim> & hyp, std::vector<comb<dim>> & w_comb, size_t d)
	{
		// expand the intersection of the wavefronts

		std::vector<comb<dim>> q_comb = SubHyperCube<dim,dim-1>::getCombinations_R(w_comb[d],dim-2);

		// Eliminate the w_comb[d] direction

79
		for (size_t k = 0 ; k < q_comb.size() ; k++)
80
		{
81
			for (size_t j = 0 ; j < dim ; j++)
82 83 84 85 86 87 88 89 90
			{
				if (w_comb[d].c[j] != 0)
				{
					q_comb[k].c[j] = 0;
				}
			}
		}

		// for all the combinations
91
		for (size_t j = 0 ; j < q_comb.size() ; j++)
92 93 94 95 96 97 98 99
		{
			size_t id = hyp.LinId(q_comb[j]);

			// get the combination of the direction d
			bool is_pos = hyp.isPositive(d);

			// is positive, modify the stop point or the starting point

100
			for (size_t s = 0 ; s < dim ; s++)
101 102 103 104 105 106 107 108 109
			{
				if (is_pos == true)
				{v_w.template get<wavefront<dim>::stop>(id)[s] = v_w.template get<wavefront<dim>::stop>(id)[s] + w_comb[d].c[s];}
				else
				{v_w.template get<wavefront<dim>::start>(id)[s] = v_w.template get<wavefront<dim>::start>(id)[s] + w_comb[d].c[s];}
			}
		}
	}

Pietro Incardona's avatar
Pietro Incardona committed
110
	/*! \brief Fill the wavefront position
incardon's avatar
incardon committed
111 112 113 114
	 *
	 * \tparam prp property to set
	 *
	 * \param graph we are processing
Pietro Incardona's avatar
Pietro Incardona committed
115
	 * \param v_w array of wavefronts
incardon's avatar
incardon committed
116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133
	 *
	 */
	template<unsigned int prp> void write_wavefront(Graph & graph,openfpm::vector<wavefront<dim>> & v_w)
	{
		// fill the wall domain with 0

		fill_domain<prp>(graph,gh.getBox(),0);

		// fill the wavefront

		for (int i = 0 ; i < v_w.size() ; i++)
		{
			Box<dim,size_t> box = wavefront<dim>::getBox(v_w.get(i));

			fill_domain<prp>(graph,box,1);
		}
	}

Pietro Incardona's avatar
Pietro Incardona committed
134
	/*! \brief Fill the domain
incardon's avatar
incardon committed
135 136 137 138
	 *
	 * \tparam p_sub property to set with the sub-domain id
	 *
	 * \param graph we are processing
Pietro Incardona's avatar
Pietro Incardona committed
139
	 * \param box Box to fill
incardon's avatar
incardon committed
140 141 142 143
	 * \param ids value to fill with
	 *
	 */

144
	template<unsigned int p_sub> void fill_domain(Graph & graph,const Box<dim,size_t> & box, long int ids)
incardon's avatar
incardon committed
145 146
	{
		// Create a subgrid iterator
incardon's avatar
incardon committed
147
		grid_key_dx_iterator_sub<dim,no_stencil,do_not_print_warning_on_adjustment<dim>> g_sub(gh,box.getKP1(),box.getKP2());
incardon's avatar
incardon committed
148 149 150 151 152 153 154 155 156 157 158

		// iterate through all grid points

		while (g_sub.isNext())
		{
			// get the actual key
			const grid_key_dx<dim> & gk = g_sub.get();

			// get the vertex and set the sub id

			graph.vertex(gh.LinId(gk)).template get<p_sub>() = ids;
incardon's avatar
incardon committed
159 160 161

			// next subdomain
			++g_sub;
incardon's avatar
incardon committed
162 163 164
		}
	}

Pietro Incardona's avatar
Pietro Incardona committed
165
	/*! \brief Add the boundary domain of id p_id to the queue
incardon's avatar
incardon committed
166
	 *
Pietro Incardona's avatar
Pietro Incardona committed
167 168
	 * \tparam p_sub property id where to store the sub-domain decomposition
	 * \tparam p_id property id where is stored the decomposition
incardon's avatar
incardon committed
169
	 *
Pietro Incardona's avatar
Pietro Incardona committed
170 171
	 * \param domains vector with sub-sub-domains still to process
	 * \param v_w array of wave-fronts
incardon's avatar
incardon committed
172
	 * \param graph we are processing
Pietro Incardona's avatar
Pietro Incardona committed
173 174 175
	 * \param w_comb wavefront combination, it is the normal vector to the wavefront
	 * \param pr_id processor id for which we are optimizing the decomposition
	 * \param bc boundary conditions
incardon's avatar
incardon committed
176 177
	 *
	 */
178
	template<unsigned int p_sub, unsigned int p_id> void add_to_queue(openfpm::vector<size_t> & domains, openfpm::vector<wavefront<dim>> & v_w, Graph & graph,  std::vector<comb<dim>> & w_comb, long int pr_id, const size_t(& bc)[dim])
incardon's avatar
incardon committed
179 180 181 182 183 184 185 186
	{
		// create a new queue
		openfpm::vector<size_t> domains_new;

		// push in the new queue, the old domains of the queue that are not assigned element

		for (size_t j = 0 ; j < domains.size() ; j++)
		{
incardon's avatar
incardon committed
187 188
			long int gs = graph.vertex(domains.get(j)).template get<p_sub>();
			if (gs < 0)
incardon's avatar
incardon committed
189 190 191 192 193 194 195
			{
				// not assigned push it

				domains_new.add(domains.get(j));
			}
		}

196 197
		// Create an Hyper-cube
		HyperCube<dim> hyp;
incardon's avatar
incardon committed
198

199
		for (size_t d = 0 ; d < v_w.size() ; d++)
incardon's avatar
incardon committed
200
		{
201 202
			expand_one_wf(v_w,w_comb,d);
			adjust_others_wf(v_w,hyp,w_comb,d);
incardon's avatar
incardon committed
203 204
		}

incardon's avatar
incardon committed
205
		// for each expanded wavefront create a sub-grid iterator and add the sub-domain
incardon's avatar
incardon committed
206

207
		for (size_t d = 0 ; d < v_w.size() ; d++)
incardon's avatar
incardon committed
208 209
		{
			// Create a sub-grid iterator
incardon's avatar
incardon committed
210
			grid_key_dx_iterator_sub_bc<dim,no_stencil,do_not_print_warning_on_adjustment<dim>> g_sub(gh,v_w.template get<wavefront<dim>::start>(d),v_w.template get<wavefront<dim>::stop>(d),bc);
incardon's avatar
incardon committed
211 212 213 214 215 216 217 218

			// iterate through all grid points

			while (g_sub.isNext())
			{
				// get the actual key
				const grid_key_dx<dim> & gk = g_sub.get();

incardon's avatar
incardon committed
219
				// get the vertex and if does not have a sub-id and is assigned ...
incardon's avatar
incardon committed
220 221
				long int pid = graph.vertex(gh.LinId(gk)).template get<p_sub>();

incardon's avatar
incardon committed
222 223 224
				// Get the processor id of the sub-sub-domain
				long int pp_id = graph.vertex(gh.LinId(gk)).template get<p_id>();

Pietro Incardona's avatar
Pietro Incardona committed
225
				// if the sub-sub-domain is not assigned
incardon's avatar
incardon committed
226
				if (pid < 0)
incardon's avatar
incardon committed
227
				{
Pietro Incardona's avatar
Pietro Incardona committed
228
					// ... and we are not processing the full graph
incardon's avatar
incardon committed
229 230
					if (pr_id != -1)
					{
Pietro Incardona's avatar
Pietro Incardona committed
231
						// ... and the processor id of the sub-sub-domain match the part we are processing, add to the queue
incardon's avatar
incardon committed
232 233

						if ( pr_id == pp_id)
incardon's avatar
incardon committed
234 235 236 237
							domains_new.add(gh.LinId(gk));
					}
					else
						domains_new.add(gh.LinId(gk));
incardon's avatar
incardon committed
238
				}
incardon's avatar
incardon committed
239 240

				++g_sub;
incardon's avatar
incardon committed
241 242 243 244
			}
		}

		// copy the new queue to the old one (it not copied, C++11 move semantic)
incardon's avatar
incardon committed
245
		domains.swap(domains_new);
incardon's avatar
incardon committed
246 247
	}

Pietro Incardona's avatar
Pietro Incardona committed
248
	/*! \brief Find the biggest hyper-cube
incardon's avatar
incardon committed
249 250
	 *
	 * starting from one initial sub-domain find the biggest hyper-cube
incardon's avatar
incardon committed
251
	 * output the box, and fill a list of neighborhood processor
incardon's avatar
incardon committed
252
	 *
253 254
	 * \tparam p_sub id of the property storing the sub-decomposition
	 * \tparam p_id id of the property containing the decomposition
incardon's avatar
incardon committed
255 256
	 *
	 * \param start_p initial domain
incardon's avatar
incardon committed
257 258 259 260
	 * \param graph representing the grid of sub-sub-domain
	 * \param box produced box
	 * \param v_w Wavefronts
	 * \param w_comb wavefronts directions (0,0,1) (0,0,-1) (0,1,0) (0,-1,0) ...
incardon's avatar
incardon committed
261 262 263 264 265 266 267 268 269 270 271 272 273 274 275
	 *
	 */
	template <unsigned int p_sub, unsigned int p_id> void expand_from_point(size_t start_p, Graph & graph, Box<dim,size_t> & box, openfpm::vector<wavefront<dim>> & v_w , std::vector<comb<dim>> & w_comb)
	{
		// We assume that Graph is the rapresentation of a cartesian graph
		// this mean that the direction d is at the child d

		// Get the number of wavefronts
		size_t n_wf = w_comb.size();

		// Create an Hyper-cube
		HyperCube<dim> hyp;

		// direction of expansion

incardon's avatar
incardon committed
276
		size_t domain_id = graph.vertex(start_p).template get<p_id>();
incardon's avatar
incardon committed
277 278 279 280 281 282 283 284 285 286 287
		bool can_expand = true;

		// while is possible to expand

		while (can_expand)
		{
			// reset can expand
			can_expand = false;

			// for each direction of expansion expand the wavefront

288
			for (size_t d = 0 ; d < n_wf ; d++)
incardon's avatar
incardon committed
289 290 291 292 293 294 295
			{
				// number of processed sub-domain
				size_t n_proc_sub = 0;

				// flag to indicate if the wavefront can expand
				bool w_can_expand = true;

incardon's avatar
incardon committed
296 297 298
				// Create an iterator of the expanded wavefront
				grid_key_dx<dim> start = grid_key_dx<dim>(v_w.template get<wavefront<dim>::start>(d)) + w_comb[d];
				grid_key_dx<dim> stop = grid_key_dx<dim>(v_w.template get<wavefront<dim>::stop>(d)) + w_comb[d];
incardon's avatar
incardon committed
299
				grid_key_dx_iterator_sub<dim,no_stencil,do_not_print_warning_on_adjustment<dim>> it(gh,start,stop);
incardon's avatar
incardon committed
300

incardon's avatar
incardon committed
301
				// for each sub-domain in the expanded wavefront
incardon's avatar
incardon committed
302 303 304
				while (it.isNext())
				{
					// get the wavefront sub-domain id
incardon's avatar
incardon committed
305
					size_t sub_w_e = gh.LinId(it.get());
incardon's avatar
incardon committed
306

incardon's avatar
incardon committed
307
					// we get the processor id of the neighborhood sub-domain on direction d
incardon's avatar
incardon committed
308
					// (expanded wavefront)
incardon's avatar
incardon committed
309
					size_t exp_p = graph.vertex(sub_w_e).template get<p_id>();
incardon's avatar
incardon committed
310

incardon's avatar
incardon committed
311 312 313
					// Check if already assigned
					long int ass = graph.vertex(sub_w_e).template get<p_sub>();

incardon's avatar
incardon committed
314
					// we check if it is the same processor id and is not assigned
incardon's avatar
incardon committed
315
					w_can_expand &= ((exp_p == domain_id) & (ass < 0));
incardon's avatar
incardon committed
316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333

					// next domain
					++it;

					// increase the number of processed sub-domain
					n_proc_sub++;
				}

				// if we did not processed sub-domain, we cannot expand
				w_can_expand &= (n_proc_sub != 0);

				// if you can expand one wavefront we did not reach the end
				can_expand |= w_can_expand;

				// if we can expand the wavefront expand it
				if (w_can_expand == true)
				{
					// expand the wavefront
334
					for (size_t j = 0 ; j < dim ; j++)
incardon's avatar
incardon committed
335 336 337 338 339 340 341
					{
						v_w.template get<wavefront<dim>::stop>(d)[j] = v_w.template get<wavefront<dim>::stop>(d)[j] + w_comb[d].c[j];
						v_w.template get<wavefront<dim>::start>(d)[j] = v_w.template get<wavefront<dim>::start>(d)[j] + w_comb[d].c[j];
					}

					// expand the intersection of the wavefronts

Pietro Incardona's avatar
Pietro Incardona committed
342 343 344
					if (dim >= 2)
					{
						std::vector<comb<dim>> q_comb = SubHyperCube<dim,dim-1>::getCombinations_R(w_comb[d],dim-2);
incardon's avatar
incardon committed
345

Pietro Incardona's avatar
Pietro Incardona committed
346
						// Eliminate the w_comb[d] direction
incardon's avatar
incardon committed
347

Pietro Incardona's avatar
Pietro Incardona committed
348
						for (size_t k = 0 ; k < q_comb.size() ; k++)
incardon's avatar
incardon committed
349
						{
Pietro Incardona's avatar
Pietro Incardona committed
350
							for (size_t j = 0 ; j < dim ; j++)
incardon's avatar
incardon committed
351
							{
Pietro Incardona's avatar
Pietro Incardona committed
352 353 354 355
								if (w_comb[d].c[j] != 0)
								{
									q_comb[k].c[j] = 0;
								}
incardon's avatar
incardon committed
356 357 358
							}
						}

Pietro Incardona's avatar
Pietro Incardona committed
359 360 361 362
						// for all the combinations
						for (size_t j = 0 ; j < q_comb.size() ; j++)
						{
							size_t id = hyp.LinId(q_comb[j]);
incardon's avatar
incardon committed
363

Pietro Incardona's avatar
Pietro Incardona committed
364
							// get the combination of the direction d
incardon's avatar
incardon committed
365

Pietro Incardona's avatar
Pietro Incardona committed
366
							bool is_pos = hyp.isPositive(d);
incardon's avatar
incardon committed
367

Pietro Incardona's avatar
Pietro Incardona committed
368
							// is positive, modify the stop point or the starting point
incardon's avatar
incardon committed
369

Pietro Incardona's avatar
Pietro Incardona committed
370 371 372 373 374 375 376
							for (size_t s = 0 ; s < dim ; s++)
							{
								if (is_pos == true)
								{v_w.template get<wavefront<dim>::stop>(id)[s] = v_w.template get<wavefront<dim>::stop>(id)[s] + w_comb[d].c[s];}
								else
								{v_w.template get<wavefront<dim>::start>(id)[s] = v_w.template get<wavefront<dim>::start>(id)[s] + w_comb[d].c[s];}
							}
incardon's avatar
incardon committed
377 378 379 380 381 382 383 384
						}
					}
				}
			}
		}

		// get back the hyper-cube produced

385
		for (size_t i = 0 ; i < dim ; i++)
incardon's avatar
incardon committed
386 387
		{
			// get the index of the wavefront direction
incardon's avatar
incardon committed
388 389
			size_t p_f = hyp.positiveFace(i);
			size_t n_f = hyp.negativeFace(i);
incardon's avatar
incardon committed
390

incardon's avatar
incardon committed
391 392 393 394 395 396
			// set the box
			box.setHigh(i,v_w.template get<wavefront<dim>::stop>(p_f)[i]);
			box.setLow(i,v_w.template get<wavefront<dim>::start>(n_f)[i]);
		}
	}

incardon's avatar
incardon committed
397
	/*! \brief Initialize the wavefronts
incardon's avatar
incardon committed
398
	 *
Pietro Incardona's avatar
Pietro Incardona committed
399 400
	 * \param start_p starting point for the wavefront set
	 * \param v_w Wavefront array
incardon's avatar
incardon committed
401 402 403 404 405 406
	 *
	 */
	void InitializeWavefront(grid_key_dx<dim> & start_p, openfpm::vector<wavefront<dim>> & v_w)
	{
		// Wavefront to initialize

407
		for (size_t i = 0 ; i < v_w.size() ; i++)
incardon's avatar
incardon committed
408
		{
409
			for (size_t j = 0 ; j < dim ; j++)
incardon's avatar
incardon committed
410
			{
incardon's avatar
incardon committed
411 412
				v_w.template get<wavefront<dim>::start>(i)[j] = start_p.get(j);
				v_w.template get<wavefront<dim>::stop>(i)[j] = start_p.get(j);
incardon's avatar
incardon committed
413 414 415 416
			}
		}
	}

incardon's avatar
incardon committed
417 418 419 420 421
	/*! \brief Get the first seed
	 *
	 * search in the graph for one sub-domain labelled with processor id
	 * to use as seed
	 *
Pietro Incardona's avatar
Pietro Incardona committed
422 423
	 * \tparam p_id property id containing the decomposition
	 * \tparam p_sub property id that will contain the sub-domain decomposition
incardon's avatar
incardon committed
424
	 *
Pietro Incardona's avatar
Pietro Incardona committed
425
	 * \param graph Graph
incardon's avatar
incardon committed
426 427
	 * \param id processor id
	 *
Pietro Incardona's avatar
Pietro Incardona committed
428 429
	 * \return a valid seed key
	 *
incardon's avatar
incardon committed
430
	 */
431
	template<unsigned int p_id, unsigned int p_sub> grid_key_dx<dim> search_seed(Graph & graph, long int id)
incardon's avatar
incardon committed
432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452
	{
		// if no processor is selected return the first point
		if (id < -1)
		{
			grid_key_dx<dim> key;
			key.zero();

			return key;
		}

		// Create a grid iterator
		grid_key_dx_iterator<dim> g_sub(gh);

		// iterate through all grid points

		while (g_sub.isNext())
		{
			// get the actual key
			const grid_key_dx<dim> & gk = g_sub.get();

			// if the subdomain has the id we are searching stop
453
			if ((long int)graph.vertex(gh.LinId(gk)).template get<p_id>() == id && graph.vertex(gh.LinId(gk)).template get<p_sub>() == -1)
incardon's avatar
incardon committed
454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475
			{
				return gk;
			}

			++g_sub;
		}

		// If not found return an invalid key
		grid_key_dx<dim> key;
		key.invalid();

		return key;
	}


	/*! \brief optimize the graph
	 *
	 * Starting from a domain (hyper-cubic), it create wavefront at the boundary and expand
	 * the boundary until the wavefronts cannot expand any more.
	 * To the domains inside the hyper-cube one sub-id is assigned. This procedure continue until
	 * all the domain of one p_id has a sub-id
	 *
Pietro Incardona's avatar
Pietro Incardona committed
476 477
	 * \tparam p_id property containing the decomposition
	 * \tparam p_sub property to fill with the sub-domain decomposition
incardon's avatar
incardon committed
478 479
	 *
	 * \param start_p seed point
incardon's avatar
incardon committed
480
	 * \param graph we are processing
Pietro Incardona's avatar
Pietro Incardona committed
481 482 483
	 * \param pr_id Processor id (if p_id == -1 the optimization is done for all the processors)
	 * \param lb list of sub-domain boxes produced by the algorithm
	 * \param box_nn_processor for each sub-domain it list all the neighborhood processors
484
	 * \param ghe Ghost extension in sub-sub-domain units in each direction
Pietro Incardona's avatar
Pietro Incardona committed
485 486 487
	 * \param init_sub_id when true p_sub property is initially set to -1 [default true]
	 * \param sub_id starting sub_id to enumerate them [default 0]
	 * \param bc boundary conditions
488 489
	 *
	 * \return last assigned sub-id
incardon's avatar
incardon committed
490 491
	 *
	 */
492
	template <unsigned int p_sub, unsigned int p_id> size_t optimize(grid_key_dx<dim> & start_p, Graph & graph, long int pr_id, openfpm::vector<Box<dim,size_t>> & lb, openfpm::vector< openfpm::vector<size_t> > & box_nn_processor , const Ghost<dim,long int> & ghe ,const size_t (& bc)[dim], bool init_sub_id = true, size_t sub_id = 0)
incardon's avatar
incardon committed
493 494 495 496
	{
		// queue
		openfpm::vector<size_t> v_q;

497 498 499
		// box list 2
		openfpm::vector< openfpm::vector<size_t> > box_nn_processor2;

incardon's avatar
incardon committed
500 501 502 503 504 505 506 507 508
		// Create an hyper-cube
		HyperCube<dim> hyp;

		// Get the wavefront combinations
		std::vector<comb<dim>> w_comb = hyp.getCombinations_R(dim-1);

		// wavefronts
		openfpm::vector<wavefront<dim>> v_w(w_comb.size());

incardon's avatar
incardon committed
509
		// fill the sub decomposition with negative number
510 511 512

		if (init_sub_id == true)
			fill_domain<p_sub>(graph,gh.getBox(),-1);
incardon's avatar
incardon committed
513 514 515 516 517 518 519 520 521

		// push the first domain
		v_q.add(gh.LinId(start_p));

		while (v_q.size() != 0)
		{
			// Box
			Box<dim,size_t> box;

incardon's avatar
incardon committed
522 523 524 525 526 527
			// Get the grid_key position from the linearized id
			start_p = gh.InvLinId(v_q.get(0));

			// Initialize the wavefronts from the domain start_p
			InitializeWavefront(start_p,v_w);

incardon's avatar
incardon committed
528
			// Create the biggest box containing the domain
incardon's avatar
incardon committed
529
			expand_from_point<p_sub,p_id>(v_q.get(0),graph,box,v_w,w_comb);
incardon's avatar
incardon committed
530 531 532

			// Add the created box to the list of boxes
			lb.add(box);
incardon's avatar
incardon committed
533 534 535 536

			// fill the domain
			fill_domain<p_sub>(graph,box,sub_id);

incardon's avatar
incardon committed
537
			// add the surrounding sub-domain to the queue
538
			add_to_queue<p_sub,p_id>(v_q,v_w,graph,w_comb,pr_id,bc);
incardon's avatar
incardon committed
539 540 541 542

			// increment the sub_id
			sub_id++;
		}
543 544 545 546

		return sub_id;
	}

547
	/*! \brief Construct the sub-domain processor list
Pietro Incardona's avatar
Pietro Incardona committed
548 549
	 *
	 * \tparam p_id property that contain the decomposition
550 551 552
	 *
	 * Each entry is a sub-domain, the list of numbers indicate the neighborhood processors
	 *
Pietro Incardona's avatar
Pietro Incardona committed
553 554 555 556 557 558
	 * \param graph graph to process
	 * \param box_nn_processor for each sub-domain it list all the neighborhood processors
	 * \param subs vector of sub-domains
	 * \param ghe ghost extensions
	 * \param bc boundary conditions
	 * \param pr_id processor that we are processing
559 560 561 562 563 564 565 566 567 568 569
	 *
	 */
	template<unsigned int p_id> void construct_box_nn_processor(Graph & graph, openfpm::vector< openfpm::vector<size_t> > & box_nn_processor, const openfpm::vector<Box<dim,size_t>> & subs, const Ghost<dim,long int> & ghe, const size_t (& bc)[dim], long int pr_id)
	{
		std::unordered_map<size_t,size_t> map;

		for (size_t i = 0 ; i < subs.size() ; i++)
		{
			map.clear();
			Box<dim,size_t> sub = subs.get(i);
			sub.enlarge(ghe);
Pietro Incardona's avatar
Pietro Incardona committed
570

571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593
			grid_skin_iterator_bc<dim> gsi(gh,subs.get(i),sub,bc);

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

				size_t pp_id = graph.vertex(gh.LinId(key)).template get<p_id>();
				if (pr_id != (long int)pp_id)
					map[pp_id] = pp_id;

				++gsi;
			}

			// Add the keys to box_nn_processors

			box_nn_processor.add();
			for ( auto it = map.begin(); it != map.end(); ++it )
			{
				box_nn_processor.last().add(it->first);
			}
		}
	}

594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615
public:

	/*! \brief Constructor
	 *
	 * \param g Graph to simplify
	 * \param sz size of the grid on each dimension
	 *
	 */

	dec_optimizer(Graph & g, const size_t (& sz)[dim])
	:gh(sz)
	{
		// The graph g is suppose to represent a cartesian grid
		// No check is performed on g
	}

	/*! \brief optimize the graph
	 *
	 * Starting from a sub-sub-domain, it create wavefronts at the boundary and expand
	 * the boundary until the wavefronts cannot expand any more, creating a sub-domain covering more sub-sub-domain.
	 * This procedure continue until all the domain is covered by a sub-domains
	 *
Pietro Incardona's avatar
Pietro Incardona committed
616 617
	 * \tparam p_id property containing the processor decomposition
	 * \tparam p_sub property to fill with the sub-domain decomposition
618 619 620
	 *
	 * \param start_p seed point
	 * \param graph we are processing
Pietro Incardona's avatar
Pietro Incardona committed
621 622
	 * \param ghe ghost size
	 * \param bc boundary conditions
623 624
	 *
	 */
625
	template <unsigned int p_sub, unsigned int p_id> void optimize(grid_key_dx<dim> & start_p, Graph & graph, const Ghost<dim,long int> & ghe , const size_t (& bc)[dim])
626 627 628 629 630 631 632 633
	{
		// temporal vector
		openfpm::vector<Box<dim,size_t>> tmp;

		// temporal vector
		openfpm::vector< openfpm::vector<size_t> > box_nn_processor;

		// optimize
634
		optimize<p_sub,p_id>(start_p,graph,-1,tmp, box_nn_processor,ghe,bc);
635 636 637 638 639 640 641 642
	}

	/*! \brief optimize the graph
	 *
	 * Starting from a sub-sub-domain, it create wavefronts at the boundary and expand
	 * the boundary until the wavefronts cannot expand any more, creating a sub-domain covering more sub-sub-domain.
	 * This procedure continue until all the sub-domain of the processor p_id are covered by a sub-domains
	 *
Pietro Incardona's avatar
Pietro Incardona committed
643 644
	 * \tparam p_id property containing the decomposition
	 * \tparam p_sub property to fill with the sub-domain decomposition
645 646
	 *
	 * \param graph we are processing
Pietro Incardona's avatar
Pietro Incardona committed
647 648 649 650
	 * \param pr_id Processor id (if p_id == -1 the optimization is done for all the processors)
	 * \param lb list of sub-domain boxes
	 * \param box_nn_processor for each sub-domain it list all the neighborhood processors
	 * \param ghe ghost size
651 652
	 *
	 */
653
	template <unsigned int p_sub, unsigned int p_id> void optimize(Graph & graph, long int pr_id, openfpm::vector<Box<dim,size_t>> & lb, openfpm::vector< openfpm::vector<size_t> > & box_nn_processor, const Ghost<dim,long int> & ghe, const size_t (& bc)[dim])
654
	{
655 656 657 658 659 660
		grid_key_dx<dim> key_seed;
		key_seed.zero();

		// if processor is -1 call optimize with -1 to do on all processors and exit
		if (pr_id == -1)
		{
661
			optimize<p_sub,p_id>(key_seed,graph,pr_id,lb,box_nn_processor,ghe,bc);
662 663 664 665

			// Construct box box_nn_processor from the constructed domain
			construct_box_nn_processor<p_id>(graph,box_nn_processor,lb,ghe,bc,pr_id);

666 667 668
			return;
		}

669 670 671 672 673
		size_t sub_id = 0;

		// fill the sub decomposition with negative number
		fill_domain<p_sub>(graph,gh.getBox(),-1);

674
		key_seed = search_seed<p_id,p_sub>(graph,pr_id);
675 676 677 678

		while (key_seed.isValid())
		{
			// optimize
679
			sub_id = optimize<p_sub,p_id>(key_seed,graph,pr_id,lb,box_nn_processor,ghe,bc,false,sub_id);
680 681 682 683

			// new seed
			key_seed = search_seed<p_id,p_sub>(graph,pr_id);
		}
684 685 686

		// Construct box box_nn_processor from the constructed domain
		construct_box_nn_processor<p_id>(graph,box_nn_processor,lb,ghe,bc,pr_id);
incardon's avatar
incardon committed
687 688 689 690
	}
};

#endif