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

4 5
#include "Grid/iterators/grid_key_dx_iterator_sub.hpp"

incardon's avatar
incardon committed
6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
/*! \brief this class represent a wavefront of dimension dim
 *
 * \dim Dimensionality of the wavefront (dimensionality of the space
 *                                       where it live so the wavefront
 *                                       is dim-1)
 *
 */

template <unsigned int dim>
class wavefront
{
public:

	typedef boost::fusion::vector<size_t[dim],size_t[dim]> type;

	typedef typename memory_traits_inte<type>::type memory_int;
	typedef typename memory_traits_lin<type>::type memory_lin;

	type data;

	static const int start = 0;
	static const int stop = 1;
	static const int max_prop = 2;

	/* \brief Get the key to the point 1
	 *
	 * \return the key to the point 1
	 *
	 */

	grid_key_dx<dim> getKP1()
	{
		// grid key to return
		grid_key_dx<dim> ret(boost::fusion::at_c<start>(data));

		return ret;
	}

	/* \brief Get the key to point 2
	 *
	 * \return the key to the point 2
	 *
	 */

	grid_key_dx<dim> getKP2()
	{
		// grid key to return
		grid_key_dx<dim> ret(boost::fusion::at_c<stop>(data));

		return ret;
	}

	/* \brief produce a box from an encapsulated object
	 *
	 * \param encap encapsulated object
	 *
	 */

incardon's avatar
incardon committed
64
	template<typename encap> static Box<dim,size_t> getBox(const encap && enc)
incardon's avatar
incardon committed
65 66 67 68 69
	{
		Box<dim,size_t> bx;

		// Create the object from the encapsulation

incardon's avatar
incardon committed
70
		getBox(enc,bx);
incardon's avatar
incardon committed
71 72 73 74 75 76 77 78 79 80

		return bx;
	}

	/* \brief produce a box from an encapsulated object
	 *
	 * \param encap encapsulated object
	 *
	 */

incardon's avatar
incardon committed
81
	template<typename encap> static void getBox(const encap & enc, Box<dim,size_t> & bx)
incardon's avatar
incardon committed
82 83 84 85 86
	{
		// Create the object from the encapsulation

		for (int i = 0 ; i < dim ; i++)
		{
incardon's avatar
incardon committed
87 88
			bx.setLow(i,enc.template get<wavefront::start>()[i]);
			bx.setHigh(i,enc.template get<wavefront::stop>()[i]);
incardon's avatar
incardon committed
89 90 91 92 93 94 95
		}
	}
};

/*! \brief This class take a graph representing the space decomposition and produce a
 *         simplified version
 *
incardon's avatar
incardon committed
96 97 98
 * 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
99 100 101 102 103 104 105 106
 *
 */

template <unsigned int dim, typename Graph>
class dec_optimizer
{
	// create a grid header for helping

107
	grid_sm<dim,void> gh;
incardon's avatar
incardon committed
108 109 110

private:

111 112 113 114 115
	/*! \brief Expand one wavefront
	 *
	 * \param v_w wavefronts
	 * \param w_comb wavefront expansion combinations
	 * \param d direction of expansion
116
	 * \param bc boundary condition
117 118 119 120
	 *
	 */
	void expand_one_wf(openfpm::vector<wavefront<dim>> & v_w, std::vector<comb<dim>> & w_comb , size_t d)
	{
121
		for (size_t j = 0 ; j < dim ; j++)
122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141
		{
			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
	 *
	 * \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

142
		for (size_t k = 0 ; k < q_comb.size() ; k++)
143
		{
144
			for (size_t j = 0 ; j < dim ; j++)
145 146 147 148 149 150 151 152 153
			{
				if (w_comb[d].c[j] != 0)
				{
					q_comb[k].c[j] = 0;
				}
			}
		}

		// for all the combinations
154
		for (size_t j = 0 ; j < q_comb.size() ; j++)
155 156 157 158 159 160 161 162 163
		{
			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

164
			for (size_t s = 0 ; s < dim ; s++)
165 166 167 168 169 170 171 172 173
			{
				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
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 209
	/* \brief Fill the wavefront position
	 *
	 * \tparam prp property to set
	 *
	 * \param graph we are processing
	 * \param Box to fill
	 * \param id value to fill with
	 *
	 */

	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);
		}
	}

	/* \brief Fill the domain
	 *
	 * \tparam p_sub property to set with the sub-domain id
	 *
	 * \param graph we are processing
	 * \param Box to fill
	 * \param ids value to fill with
	 *
	 */

210
	template<unsigned int p_sub> void fill_domain(Graph & graph,const Box<dim,size_t> & box, long int ids)
incardon's avatar
incardon committed
211 212
	{
		// Create a subgrid iterator
213
		grid_key_dx_iterator_sub<dim,do_not_print_warning_on_adjustment<dim>> g_sub(gh,box.getKP1(),box.getKP2());
incardon's avatar
incardon committed
214 215 216 217 218 219 220 221 222 223 224

		// 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
225 226 227

			// next subdomain
			++g_sub;
incardon's avatar
incardon committed
228 229 230
		}
	}

incardon's avatar
incardon committed
231
	/* \brief Add the boundary domain of id p_id to the queue
incardon's avatar
incardon committed
232 233 234 235 236 237
	 *
	 * \tparam i-property where is stored the decomposition
	 *
	 * \param domains vector with domains to process
	 * \param graph we are processing
	 * \param w_comb hyper-cube combinations
incardon's avatar
incardon committed
238
	 * \param p_id processor id
incardon's avatar
incardon committed
239
	 * \param box_nn_processor list of neighborhood processors for the box
240
	 * \param bc Boundary conditions
incardon's avatar
incardon committed
241 242
	 *
	 */
243
	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, openfpm::vector< openfpm::vector<size_t> > & box_nn_processor, const size_t(& bc)[dim])
incardon's avatar
incardon committed
244 245 246 247 248 249 250 251
	{
		// 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
252 253
			long int gs = graph.vertex(domains.get(j)).template get<p_sub>();
			if (gs < 0)
incardon's avatar
incardon committed
254 255 256 257 258 259 260
			{
				// not assigned push it

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

261 262
		// Create an Hyper-cube
		HyperCube<dim> hyp;
incardon's avatar
incardon committed
263

264
		for (size_t d = 0 ; d < v_w.size() ; d++)
incardon's avatar
incardon committed
265
		{
266 267
			expand_one_wf(v_w,w_comb,d);
			adjust_others_wf(v_w,hyp,w_comb,d);
incardon's avatar
incardon committed
268 269
		}

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

272
		for (size_t d = 0 ; d < v_w.size() ; d++)
incardon's avatar
incardon committed
273 274
		{
			// Create a sub-grid iterator
275
			grid_key_dx_iterator_sub_bc<dim,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
276 277 278 279 280 281 282 283

			// 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
284
				// get the vertex and if does not have a sub-id and is assigned ...
incardon's avatar
incardon committed
285 286
				long int pid = graph.vertex(gh.LinId(gk)).template get<p_sub>();

incardon's avatar
incardon committed
287 288 289 290 291 292 293 294 295 296
				// Get the processor id of the sub-sub-domain
				long int pp_id = graph.vertex(gh.LinId(gk)).template get<p_id>();

				if (pp_id != pr_id)
				{
					// this box is contiguous to other processors

					box_nn_processor.get(box_nn_processor.size()-1).add(pp_id);
				}

297
				// if the sub-sub-domain is not assigned
incardon's avatar
incardon committed
298
				if (pid < 0)
incardon's avatar
incardon committed
299
				{
300
					// ... and we are not processing the full graph
incardon's avatar
incardon committed
301 302
					if (pr_id != -1)
					{
303
						// ... and the processor id of the sub-sub-domain match the part we are processing, add to the queue
incardon's avatar
incardon committed
304 305

						if ( pr_id == pp_id)
incardon's avatar
incardon committed
306 307 308 309
							domains_new.add(gh.LinId(gk));
					}
					else
						domains_new.add(gh.LinId(gk));
incardon's avatar
incardon committed
310
				}
incardon's avatar
incardon committed
311 312

				++g_sub;
incardon's avatar
incardon committed
313 314 315
			}
		}

incardon's avatar
incardon committed
316 317 318 319 320
		// sort and make it unique the numbers in processor list
	    std::sort(box_nn_processor.get(box_nn_processor.size()-1).begin(), box_nn_processor.get(box_nn_processor.size()-1).end());
	    auto last = std::unique(box_nn_processor.get(box_nn_processor.size()-1).begin(), box_nn_processor.get(box_nn_processor.size()-1).end());
	    box_nn_processor.get(box_nn_processor.size()-1).erase(last, box_nn_processor.get(box_nn_processor.size()-1).end());

incardon's avatar
incardon committed
321
		// copy the new queue to the old one (it not copied, C++11 move semantic)
incardon's avatar
incardon committed
322
		domains.swap(domains_new);
incardon's avatar
incardon committed
323 324 325 326 327
	}

	/* \brief Find the biggest hyper-cube
	 *
	 * starting from one initial sub-domain find the biggest hyper-cube
incardon's avatar
incardon committed
328
	 * output the box, and fill a list of neighborhood processor
incardon's avatar
incardon committed
329
	 *
330 331
	 * \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
332 333
	 *
	 * \param start_p initial domain
incardon's avatar
incardon committed
334 335 336 337
	 * \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
338 339 340 341 342 343 344 345 346 347 348 349 350 351 352
	 *
	 */
	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
353
		size_t domain_id = graph.vertex(start_p).template get<p_id>();
incardon's avatar
incardon committed
354 355 356 357 358 359 360 361 362 363 364
		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

365
			for (size_t d = 0 ; d < n_wf ; d++)
incardon's avatar
incardon committed
366 367 368 369 370 371 372
			{
				// 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
373 374 375
				// 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];
376
				grid_key_dx_iterator_sub<dim,do_not_print_warning_on_adjustment<dim>> it(gh,start,stop);
incardon's avatar
incardon committed
377

incardon's avatar
incardon committed
378
				// for each sub-domain in the expanded wavefront
incardon's avatar
incardon committed
379 380 381
				while (it.isNext())
				{
					// get the wavefront sub-domain id
incardon's avatar
incardon committed
382
					size_t sub_w_e = gh.LinId(it.get());
incardon's avatar
incardon committed
383

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

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

incardon's avatar
incardon committed
391
					// we check if it is the same processor id and is not assigned
incardon's avatar
incardon committed
392
					w_can_expand &= ((exp_p == domain_id) & (ass < 0));
incardon's avatar
incardon committed
393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410

					// 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
411
					for (size_t j = 0 ; j < dim ; j++)
incardon's avatar
incardon committed
412 413 414 415 416 417 418 419 420 421 422
					{
						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

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

					// Eliminate the w_comb[d] direction

423
					for (size_t k = 0 ; k < q_comb.size() ; k++)
incardon's avatar
incardon committed
424
					{
425
						for (size_t j = 0 ; j < dim ; j++)
incardon's avatar
incardon committed
426 427 428 429 430 431 432 433 434
						{
							if (w_comb[d].c[j] != 0)
							{
								q_comb[k].c[j] = 0;
							}
						}
					}

					// for all the combinations
435
					for (size_t j = 0 ; j < q_comb.size() ; j++)
incardon's avatar
incardon committed
436 437 438 439 440 441 442 443 444
					{
						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

445
						for (size_t s = 0 ; s < dim ; s++)
incardon's avatar
incardon committed
446 447
						{
							if (is_pos == true)
incardon's avatar
incardon committed
448
							{v_w.template get<wavefront<dim>::stop>(id)[s] = v_w.template get<wavefront<dim>::stop>(id)[s] + w_comb[d].c[s];}
incardon's avatar
incardon committed
449
							else
incardon's avatar
incardon committed
450
							{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
451 452 453 454 455 456 457 458
						}
					}
				}
			}
		}

		// get back the hyper-cube produced

459
		for (size_t i = 0 ; i < dim ; i++)
incardon's avatar
incardon committed
460 461
		{
			// get the index of the wavefront direction
incardon's avatar
incardon committed
462 463
			size_t p_f = hyp.positiveFace(i);
			size_t n_f = hyp.negativeFace(i);
incardon's avatar
incardon committed
464

incardon's avatar
incardon committed
465 466 467 468 469 470
			// 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
471
	/*! \brief Initialize the wavefronts
incardon's avatar
incardon committed
472
	 *
incardon's avatar
incardon committed
473
	 * \param starting point of the wavefront set
incardon's avatar
incardon committed
474 475 476 477 478 479 480
	 * \param v_w Wavefront to initialize
	 *
	 */
	void InitializeWavefront(grid_key_dx<dim> & start_p, openfpm::vector<wavefront<dim>> & v_w)
	{
		// Wavefront to initialize

481
		for (size_t i = 0 ; i < v_w.size() ; i++)
incardon's avatar
incardon committed
482
		{
483
			for (size_t j = 0 ; j < dim ; j++)
incardon's avatar
incardon committed
484
			{
incardon's avatar
incardon committed
485 486
				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
487 488 489 490
			}
		}
	}

incardon's avatar
incardon committed
491 492 493 494 495 496 497 498 499 500 501 502
	/*! \brief Get the first seed
	 *
	 * search in the graph for one sub-domain labelled with processor id
	 * to use as seed
	 *
	 * \tparam p_id property in the graph storing the sub-domain id
	 *
	 * \param Graph graph
	 * \param id processor id
	 *
	 */

503
	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
504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524
	{
		// 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
525
			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
526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551
			{
				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
	 *
	 * \tparam j property containing the decomposition
	 * \tparam i property to fill with the sub-decomposition
	 *
	 * \param start_p seed point
incardon's avatar
incardon committed
552 553
	 * \param graph we are processing
	 * \param p_id Processor id (if p_id == -1 the optimization is done for all the processors)
554
	 * \param list of sub-domain boxes produced by the algorithm
incardon's avatar
incardon committed
555
	 * \param box_nn_processor for each box it list all the neighborhood processor
556
	 * \param bc Boundary condition
557 558 559 560
	 * \param init_sub_id when true p_sub property is initial set to -1 [default true]
	 * \param sub_id starting sub_id enumeration [default 0]
	 *
	 * \return last assigned sub-id
incardon's avatar
incardon committed
561 562
	 *
	 */
563
	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 size_t (& bc)[dim], bool init_sub_id = true, size_t sub_id = 0)
incardon's avatar
incardon committed
564 565 566 567 568 569 570 571 572 573 574 575 576
	{
		// queue
		openfpm::vector<size_t> v_q;

		// 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
577
		// fill the sub decomposition with negative number
578 579 580

		if (init_sub_id == true)
			fill_domain<p_sub>(graph,gh.getBox(),-1);
incardon's avatar
incardon committed
581 582 583 584 585 586 587 588 589

		// 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
590 591 592 593 594 595
			// 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
596 597 598
			// Create a list for the box
			box_nn_processor.add();

incardon's avatar
incardon committed
599
			// Create the biggest box containing the domain
incardon's avatar
incardon committed
600
			expand_from_point<p_sub,p_id>(v_q.get(0),graph,box,v_w,w_comb);
incardon's avatar
incardon committed
601 602 603

			// Add the created box to the list of boxes
			lb.add(box);
incardon's avatar
incardon committed
604 605 606 607

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

incardon's avatar
incardon committed
608
			// add the surrounding sub-domain to the queue
609
			add_to_queue<p_sub,p_id>(v_q,v_w,graph,w_comb,pr_id,box_nn_processor,bc);
incardon's avatar
incardon committed
610 611 612 613

			// increment the sub_id
			sub_id++;
		}
614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 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 669 670 671 672 673 674

		return sub_id;
	}

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
	 *
	 * \tparam j property containing the processor decomposition
	 * \tparam i property to fill with the sub-domain-decomposition id
	 *
	 * \param start_p seed point
	 * \param graph we are processing
	 *
	 */
	template <unsigned int p_sub, unsigned int p_id> void optimize(grid_key_dx<dim> & start_p, Graph & graph, const size_t (& bc)[dim])
	{
		// temporal vector
		openfpm::vector<Box<dim,size_t>> tmp;

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

		// optimize
		optimize<p_sub,p_id>(start_p,graph,-1,tmp, box_nn_processor,bc);
	}

	/*! \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
	 *
	 * \tparam j property containing the decomposition
	 * \tparam i property to fill with the sub-domain-decomposition id
	 *
	 * \param graph we are processing
	 * \param p_id Processor id (if p_id == -1 the optimization is done for all the processors)
	 * \param list of sub-domain boxes
	 *
	 */
	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 size_t (& bc)[dim])
	{
675 676 677 678 679 680 681 682 683 684
		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)
		{
			optimize<p_sub,p_id>(key_seed,graph,pr_id,lb,box_nn_processor,bc);
			return;
		}

685 686 687 688 689
		size_t sub_id = 0;

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

690
		key_seed = search_seed<p_id,p_sub>(graph,pr_id);
691 692 693 694 695 696 697 698 699

		while (key_seed.isValid())
		{
			// optimize
			sub_id = optimize<p_sub,p_id>(key_seed,graph,pr_id,lb,box_nn_processor,bc,false,sub_id);

			// new seed
			key_seed = search_seed<p_id,p_sub>(graph,pr_id);
		}
incardon's avatar
incardon committed
700 701 702 703
	}
};

#endif