FDScheme.hpp 25.2 KB
Newer Older
incardon's avatar
incardon committed
1 2 3 4 5 6 7 8 9 10
/*
 * FiniteDifferences.hpp
 *
 *  Created on: Sep 17, 2015
 *      Author: i-bird
 */

#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_HPP_
#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_HPP_

11
#include "../Matrix/SparseMatrix.hpp"
incardon's avatar
incardon committed
12
#include "Grid/grid_dist_id.hpp"
13
#include "Grid/Iterators/grid_dist_id_iterator_sub.hpp"
incardon's avatar
incardon committed
14
#include "eq.hpp"
15
#include "NN/CellList/CellDecomposer.hpp"
16
#include "Grid/staggered_dist_grid_util.hpp"
17 18 19
#include "Grid/grid_dist_id.hpp"
#include "Vector/Vector_util.hpp"
#include "Grid/staggered_dist_grid.hpp"
incardon's avatar
incardon committed
20 21

/*! \brief Finite Differences
22
 *
23 24 25
 * This class is able to discretize on a Matrix any system of equations producing a linear system of type \f$Ax=b\f$. In order to create a consistent
 * Matrix it is required that each processor must contain a contiguous range on grid points without
 * holes. In order to ensure this, each processor produce a contiguous local labeling of its local
26 27
 * points. Each processor also add an offset equal to the number of local
 * points of the processors with id smaller than him, to produce a global and non overlapping
28 29
 * labeling. An example is shown in the figures down, here we have
 * a grid 8x6 divided across four processors each processor label locally its grid points
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 64 65 66 67 68 69 70 71
 *
 * \verbatim
 *
+--------------------------+
| 1   2   3   4| 1  2  3  4|
|              |           |
| 5   6   7   8| 5  6  7  8|
|              |           |
| 9  10  11  12| 9 10 11 12|
+--------------------------+
|13  14  15| 13 14 15 16 17|
|          |               |
|16  17  18| 18 19 20 21 22|
|          |               |
|19  20  21| 23 24 25 26 27|
+--------------------------+

 *
 *
 * \endverbatim
 *
 * To the local relabelling is added an offset to make the local id global and non overlapping
 *
 *
 * \verbatim
 *
+--------------------------+
| 1   2   3   4|23 24 25 26|
|              |           |
| 5   6   7   8|27 28 29 30|
|              |           |
| 9  10  12  13|31 32 33 34|
+--------------------------+
|14  15  16| 35 36 37 38 39|
|          |               |
|17  18  19| 40 41 42 43 44|
|          |               |
|20  21  22| 45 46 47 48 49|
+--------------------------+
 *
 *
 * \endverbatim
incardon's avatar
incardon committed
72
 *
73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
 * \tparam Sys_eqs Definition of the system of equations
 *
 * # Examples
 *
 * ## Solve lid-driven cavity 2D for incompressible fluid (inertia=0 --> Re=0)
 *
 * In this case the system of equation to solve is
 *
 * \f$
\left\{
\begin{array}{c}
\eta\nabla v_x + \partial_x P = 0 \quad Eq1 \\
\eta\nabla v_y + \partial_y P = 0 \quad Eq2 \\
\partial_x v_x + \partial_y v_y = 0 \quad Eq3
\end{array}
\right.  \f$

and  boundary conditions

 * \f$
\left\{
\begin{array}{c}
v_x = 0, v_y = 0 \quad x = 0 \quad B1\\
v_x = 0, v_y = 1.0 \quad x = L \quad B2\\
v_x = 0, v_y = 0 \quad y = 0 \quad B3\\
v_x = 0, v_y = 0 \quad y = L \quad B4\\
\end{array}
\right.  \f$

 *
 * with \f$v_x\f$ and \f$v_y\f$ the velocity in x and y and \f$P\f$ Pressure
 *
 * In order to solve such system first we define the general properties of the system
 *
 *	\snippet eq_unit_test.hpp Definition of the system
 *
 * ## Define the equations of the system
 *
 * \snippet eq_unit_test.hpp Definition of the equation of the system in the bulk and at the boundary
 *
 * ## Define the domain and impose the equations
 *
 * \snippet eq_unit_test.hpp lid-driven cavity 2D
 *
 * # 3D
 *
 * A 3D case is given in the examples
incardon's avatar
incardon committed
120 121 122 123 124 125
 *
 */

template<typename Sys_eqs>
class FDScheme
{
126 127
public:

128
	//! Distributed grid map
incardon's avatar
incardon committed
129
	typedef grid_dist_id<Sys_eqs::dims,typename Sys_eqs::stype,aggregate<size_t>,typename Sys_eqs::b_grid::decomposition::extended_type> g_map_type;
130

131
	//! Type that specify the properties of the system of equations
132 133 134 135
	typedef Sys_eqs Sys_eqs_typ;

private:

incardon's avatar
incardon committed
136
	//! Encapsulation of the b term as constant
137 138
	struct constant_b
	{
incardon's avatar
incardon committed
139
		//! scalar
140 141
		typename Sys_eqs::stype scal;

incardon's avatar
incardon committed
142 143 144 145 146
		/*! \brief Constrictor from a scalar
		 *
		 * \param scal scalar
		 *
		 */
147 148 149 150 151
		constant_b(typename Sys_eqs::stype scal)
		{
			this->scal = scal;
		}

incardon's avatar
incardon committed
152 153 154 155 156 157 158 159 160
		/*! \brief Get the b term on a grid point
		 *
		 * \note It does not matter the grid point it is a scalar
		 *
		 * \param  key grid position (unused because it is a constant)
		 *
		 * \return the scalar
		 *
		 */
161 162 163 164 165 166
		typename Sys_eqs::stype get(grid_dist_key_dx<Sys_eqs::dims> & key)
		{
			return scal;
		}
	};

incardon's avatar
incardon committed
167
	//! Encapsulation of the b term as grid
168 169 170
	template<typename grid, unsigned int prp>
	struct grid_b
	{
incardon's avatar
incardon committed
171
		//! b term fo the grid
172 173
		grid & gr;

incardon's avatar
incardon committed
174 175 176 177 178
		/*! \brief gr grid that encapsulate
		 *
		 * \param gr grid
		 *
		 */
179 180 181 182
		grid_b(grid & gr)
		:gr(gr)
		{}

incardon's avatar
incardon committed
183 184 185 186 187 188 189
		/*! \brief Get the value of the b term on a grid point
		 *
		 * \param key grid point
		 *
		 * \return the value
		 *
		 */
190 191 192 193 194 195
		typename Sys_eqs::stype get(grid_dist_key_dx<Sys_eqs::dims> & key)
		{
			return gr.template get<prp>(key);
		}
	};

196
	//! Padding
197 198
	Padding<Sys_eqs::dims> pd;

199
	//! Sparse matrix triplet type
200 201
	typedef typename Sys_eqs::SparseMatrix_type::triplet_type triplet;

202
	//! Vector b
203
	typename Sys_eqs::Vector_type b;
204

205
	//! Domain Grid informations
206 207
	const grid_sm<Sys_eqs::dims,void> & gs;

208
	//! Get the grid spacing
209
	typename Sys_eqs::stype spacing[Sys_eqs::dims];
210

211
	//! mapping grid
212
	g_map_type g_map;
213

214
	//! row of the matrix
215 216
	size_t row;

217
	//! row on b
218 219
	size_t row_b;

220
	//! Grid points that has each processor
221 222
	openfpm::vector<size_t> pnt;

223
	//! Staggered position for each property
224 225
	comb<Sys_eqs::dims> s_pos[Sys_eqs::nvar];

226 227 228 229
	//! Each point in the grid has a global id, to decompose correctly the Matrix each processor contain a
	//! contiguos range of global id, example processor 0 can have from 0 to 234 and processor 1 from 235 to 512
	//! no processors can have holes in the sequence, this number indicate where the sequence start for this
	//! processor
230 231
	size_t s_pnt;

232 233 234
	/*! \brief Equation id + space position
	 *
	 */
235 236
	struct key_and_eq
	{
237
		//! space position
238
		grid_key_dx<Sys_eqs::dims> key;
239 240

		//! equation id
241 242 243 244 245
		size_t eq;
	};

	/*! \brief From the row Matrix position to the spatial position
	 *
246
	 * \param row Matrix row
247
	 *
248
	 * \return spatial position + equation id
249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273
	 *
	 */
	inline key_and_eq from_row_to_key(size_t row)
	{
		key_and_eq ke;
		auto it = g_map.getDomainIterator();

		while (it.isNext())
		{
			size_t row_low = g_map.template get<0>(it.get());

			if (row >= row_low * Sys_eqs::nvar && row < row_low * Sys_eqs::nvar + Sys_eqs::nvar)
			{
				ke.eq = row - row_low * Sys_eqs::nvar;
				ke.key = g_map.getGKey(it.get());
				return ke;
			}

			++it;
		}
		std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " the row does not map to any position" << "\n";

		return ke;
	}

274
	/*! \brief calculate the mapping grid size with padding
275
	 *
276 277
	 * \param sz original grid size
	 * \param pd padding
278 279 280 281 282 283 284 285 286 287 288 289 290 291
	 *
	 * \return padded grid size
	 *
	 */
	inline const std::vector<size_t> padding( const size_t (& sz)[Sys_eqs::dims], Padding<Sys_eqs::dims> & pd)
	{
		std::vector<size_t> g_sz_pad(Sys_eqs::dims);

		for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
			g_sz_pad[i] = sz[i] + pd.getLow(i) + pd.getHigh(i);

		return g_sz_pad;
	}

292 293 294 295 296
	/*! \brief Check if the Matrix is consistent
	 *
	 */
	void consistency()
	{
297 298
		openfpm::vector<triplet> & trpl = A.getMatrixTriplets();

299 300 301 302 303
		// A and B must have the same rows
		if (row != row_b)
			std::cerr << "Error " << __FILE__ << ":" << __LINE__ << "the term B and the Matrix A for Ax=B must contain the same number of rows\n";

		// Indicate all the non zero rows
304
		openfpm::vector<unsigned char> nz_rows;
305 306 307
		nz_rows.resize(row_b);

		for (size_t i = 0 ; i < trpl.size() ; i++)
308
			nz_rows.get(trpl.get(i).row() - s_pnt*Sys_eqs::nvar) = true;
309 310

		// Indicate all the non zero colums
311
		// This check can be done only on single processor
312

incardon's avatar
incardon committed
313
		Vcluster<> & v_cl = create_vcluster();
314
		if (v_cl.getProcessingUnits() == 1)
315
		{
316 317
			openfpm::vector<unsigned> nz_cols;
			nz_cols.resize(row_b);
318

319 320 321 322 323
			for (size_t i = 0 ; i < trpl.size() ; i++)
				nz_cols.get(trpl.get(i).col()) = true;

			// all the rows must have a non zero element
			for (size_t i = 0 ; i < nz_rows.size() ; i++)
324
			{
325 326 327 328 329
				if (nz_rows.get(i) == false)
				{
					key_and_eq ke = from_row_to_key(i);
					std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix row " << i <<  " is not filled, position " << ke.key.to_string() << " equation: " << ke.eq << "\n";
				}
330
			}
331

332 333 334 335 336 337
			// all the colums must have a non zero element
			for (size_t i = 0 ; i < nz_cols.size() ; i++)
			{
				if (nz_cols.get(i) == false)
					std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix colum " << i << " is not filled\n";
			}
338 339 340
		}
	}

341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390
	/*! \brief Copy a given solution vector in a staggered grid
	 *
	 * \tparam Vct Vector type
	 * \tparam Grid_dst target grid
	 * \tparam pos set of properties
	 *
	 * \param v Vector
	 * \param g_dst target staggered grid
	 *
	 */
	template<typename Vct, typename Grid_dst ,unsigned int ... pos> void copy_staggered(Vct & v, Grid_dst & g_dst)
	{
		// check that g_dst is staggered
		if (g_dst.is_staggered() == false)
			std::cerr << __FILE__ << ":" << __LINE__ << " The destination grid must be staggered " << std::endl;

#ifdef SE_CLASS1

		if (g_map.getLocalDomainSize() != g_dst.getLocalDomainSize())
			std::cerr << __FILE__ << ":" << __LINE__ << " The staggered and destination grid in size does not match " << std::endl;
#endif

		// sub-grid iterator over the grid map
		auto g_map_it = g_map.getDomainIterator();

		// Iterator over the destination grid
		auto g_dst_it = g_dst.getDomainIterator();

		while (g_map_it.isNext() == true)
		{
			typedef typename to_boost_vmpl<pos...>::type vid;
			typedef boost::mpl::size<vid> v_size;

			auto key_src = g_map_it.get();
			size_t lin_id = g_map.template get<0>(key_src);

			// destination point
			auto key_dst = g_dst_it.get();

			// Transform this id into an id for the Eigen vector

			copy_ele<Sys_eqs_typ, Grid_dst,Vct> cp(key_dst,g_dst,v,lin_id,g_map.size());

			boost::mpl::for_each_ref<boost::mpl::range_c<int,0,v_size::value>>(cp);

			++g_map_it;
			++g_dst_it;
		}
	}

391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406
	/*! \brief Copy a given solution vector in a normal grid
	 *
	 * \tparam Vct Vector type
	 * \tparam Grid_dst target grid
	 * \tparam pos set of property
	 *
	 * \param v Vector
	 * \param g_dst target normal grid
	 *
	 */
	template<typename Vct, typename Grid_dst ,unsigned int ... pos> void copy_normal(Vct & v, Grid_dst & g_dst)
	{
		// check that g_dst is staggered
		if (g_dst.is_staggered() == true)
			std::cerr << __FILE__ << ":" << __LINE__ << " The destination grid must be normal " << std::endl;

incardon's avatar
incardon committed
407 408
		grid_key_dx<Grid_dst::dims> start;
		grid_key_dx<Grid_dst::dims> stop;
409

incardon's avatar
incardon committed
410 411 412 413 414
		for (size_t i = 0 ; i < Grid_dst::dims ; i++)
		{
			start.set_d(i,pd.getLow(i));
			stop.set_d(i,g_map.size(i) - pd.getHigh(i));
		}
415 416

		// sub-grid iterator over the grid map
incardon's avatar
incardon committed
417
		auto g_map_it = g_map.getSubDomainIterator(start,stop);
418 419 420 421

		// Iterator over the destination grid
		auto g_dst_it = g_dst.getDomainIterator();

incardon's avatar
incardon committed
422
		while (g_dst_it.isNext() == true)
423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458
		{
			typedef typename to_boost_vmpl<pos...>::type vid;
			typedef boost::mpl::size<vid> v_size;

			auto key_src = g_map_it.get();
			size_t lin_id = g_map.template get<0>(key_src);

			// destination point
			auto key_dst = g_dst_it.get();

			// Transform this id into an id for the vector

			copy_ele<Sys_eqs_typ, Grid_dst,Vct> cp(key_dst,g_dst,v,lin_id,g_map.size());

			boost::mpl::for_each_ref<boost::mpl::range_c<int,0,v_size::value>>(cp);

			++g_map_it;
			++g_dst_it;
		}
	}

	/*! \brief Impose an operator
	 *
	 * This function impose an operator on a particular grid region to produce the system
	 *
	 * Ax = b
	 *
	 * ## Stokes equation 2D, lid driven cavity with one splipping wall
	 * \snippet eq_unit_test.hpp Copy the solution to grid
	 *
	 * \param op Operator to impose (A term)
	 * \param num right hand side of the term (b term)
	 * \param id Equation id in the system that we are imposing
	 * \param it_d iterator that define where you want to impose
	 *
	 */
incardon's avatar
incardon committed
459 460 461 462 463 464 465 466 467 468
	template<typename T, typename bop, typename iterator> void impose_dit(const T & op ,
			                         bop num,
									 long int id ,
									 const iterator & it_d)
	{
		openfpm::vector<triplet> & trpl = A.getMatrixTriplets();

		auto it = it_d;
		grid_sm<Sys_eqs::dims,void> gs = g_map.getGridInfoVoid();

incardon's avatar
incardon committed
469
		std::unordered_map<long int,typename Sys_eqs::stype> cols;
incardon's avatar
incardon committed
470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540

		// iterate all the grid points
		while (it.isNext())
		{
			// get the position
			auto key = it.get();

			// Add padding
			for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
				key.getKeyRef().set_d(i,key.getKeyRef().get(i) + pd.getLow(i));

			// Calculate the non-zero colums
			T::value(g_map,key,gs,spacing,cols,1.0);

			// indicate if the diagonal has been set
			bool is_diag = false;

			// create the triplet
			for ( auto it = cols.begin(); it != cols.end(); ++it )
			{
				trpl.add();
				trpl.last().row() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
				trpl.last().col() = it->first;
				trpl.last().value() = it->second;

				if (trpl.last().row() == trpl.last().col())
					is_diag = true;

//				std::cout << "(" << trpl.last().row() << "," << trpl.last().col() << "," << trpl.last().value() << ")" << "\n";
			}

			// If does not have a diagonal entry put it to zero
			if (is_diag == false)
			{
				trpl.add();
				trpl.last().row() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
				trpl.last().col() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
				trpl.last().value() = 0.0;
			}

			b(g_map.template get<0>(key)*Sys_eqs::nvar + id) = num.get(key);

			cols.clear();

			// if SE_CLASS1 is defined check the position
#ifdef SE_CLASS1
//			T::position(key,gs,s_pos);
#endif

			++row;
			++row_b;
			++it;
		}
	}

	/*! \brief Impose an operator
	 *
	 * This function impose an operator on a particular grid region to produce the system
	 *
	 * Ax = b
	 *
	 * ## Stokes equation 2D, lid driven cavity with one splipping wall
	 * \snippet eq_unit_test.hpp Copy the solution to grid
	 *
	 * \param op Operator to impose (A term)
	 * \param num right hand side of the term (b term)
	 * \param id Equation id in the system that we are imposing
	 * \param it_d iterator that define where you want to impose
	 *
	 */
	template<typename T, typename bop, typename iterator> void impose_git(const T & op ,
541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606
			                         bop num,
									 long int id ,
									 const iterator & it_d)
	{
		openfpm::vector<triplet> & trpl = A.getMatrixTriplets();

		auto it = it_d;
		grid_sm<Sys_eqs::dims,void> gs = g_map.getGridInfoVoid();

		std::unordered_map<long int,float> cols;

		// iterate all the grid points
		while (it.isNext())
		{
			// get the position
			auto key = it.get();

			// Calculate the non-zero colums
			T::value(g_map,key,gs,spacing,cols,1.0);

			// indicate if the diagonal has been set
			bool is_diag = false;

			// create the triplet
			for ( auto it = cols.begin(); it != cols.end(); ++it )
			{
				trpl.add();
				trpl.last().row() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
				trpl.last().col() = it->first;
				trpl.last().value() = it->second;

				if (trpl.last().row() == trpl.last().col())
					is_diag = true;

//				std::cout << "(" << trpl.last().row() << "," << trpl.last().col() << "," << trpl.last().value() << ")" << "\n";
			}

			// If does not have a diagonal entry put it to zero
			if (is_diag == false)
			{
				trpl.add();
				trpl.last().row() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
				trpl.last().col() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
				trpl.last().value() = 0.0;
			}

			b(g_map.template get<0>(key)*Sys_eqs::nvar + id) = num.get(key);

			cols.clear();

			// if SE_CLASS1 is defined check the position
#ifdef SE_CLASS1
//			T::position(key,gs,s_pos);
#endif

			++row;
			++row_b;
			++it;
		}
	}

	/*! \brief Construct the gmap structure
	 *
	 */
	void construct_gmap()
	{
incardon's avatar
incardon committed
607
		Vcluster<> & v_cl = create_vcluster();
608 609 610 611 612 613 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

		// Calculate the size of the local domain
		size_t sz = g_map.getLocalDomainSize();

		// Get the total size of the local grids on each processors
		v_cl.allGather(sz,pnt);
		v_cl.execute();
		s_pnt = 0;

		// calculate the starting point for this processor
		for (size_t i = 0 ; i < v_cl.getProcessUnitID() ; i++)
			s_pnt += pnt.get(i);

		// resize b if needed
		b.resize(Sys_eqs::nvar * g_map.size(),Sys_eqs::nvar * sz);

		// Calculate the starting point

		// Counter
		size_t cnt = 0;

		// Create the re-mapping grid
		auto it = g_map.getDomainIterator();

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

			g_map.template get<0>(key) = cnt + s_pnt;

			++cnt;
			++it;
		}

		// sync the ghost
		g_map.template ghost_get<0>();
	}

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
	/*! \initialize the object FDScheme
	 *
	 * \param domain simulation domain
	 *
	 */
	void Initialize(const Box<Sys_eqs::dims,typename Sys_eqs::stype> & domain)
	{
		construct_gmap();

		// Create a CellDecomposer and calculate the spacing

		size_t sz_g[Sys_eqs::dims];
		for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
		{
			if (Sys_eqs::boundary[i] == NON_PERIODIC)
				sz_g[i] = gs.getSize()[i] - 1;
			else
				sz_g[i] = gs.getSize()[i];
		}

		CellDecomposer_sm<Sys_eqs::dims,typename Sys_eqs::stype> cd(domain,sz_g,0);

		for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
			spacing[i] = cd.getCellBox().getHigh(i);
	}

672 673 674
public:

	/*! \brief set the staggered position for each property
incardon's avatar
incardon committed
675
	 *
676
	 * \param sp vector containing the staggered position for each property
incardon's avatar
incardon committed
677 678
	 *
	 */
679
	void setStagPos(comb<Sys_eqs::dims> (& sp)[Sys_eqs::nvar])
incardon's avatar
incardon committed
680
	{
681 682
		for (size_t i = 0 ; i < Sys_eqs::nvar ; i++)
			s_pos[i] = sp[i];
incardon's avatar
incardon committed
683 684
	}

685 686 687
	/*! \brief compute the staggered position for each property
	 *
	 * This is compute from the value_type stored by Sys_eqs::b_grid::value_type
688
	 * the position of the staggered properties
689 690 691 692 693 694
	 *
	 *
	 */
	void computeStag()
	{
		typedef typename Sys_eqs::b_grid::value_type prp_type;
incardon's avatar
incardon committed
695

696 697 698 699 700 701 702 703
		openfpm::vector<comb<Sys_eqs::dims>> c_prp[prp_type::max_prop];

		stag_set_position<Sys_eqs::dims,prp_type> ssp(c_prp);

		boost::mpl::for_each_ref< boost::mpl::range_c<int,0,prp_type::max_prop> >(ssp);
	}

	/*! \brief Get the specified padding
704
	 *
705
	 * \return the padding specified
706 707 708 709 710 711 712 713
	 *
	 */
	const Padding<Sys_eqs::dims> & getPadding()
	{
		return pd;
	}

	/*! \brief Return the map between the grid index position and the position in the distributed vector
714 715
	 *
	 * It is the map explained in the intro of the FDScheme
716 717 718 719 720 721 722 723 724
	 *
	 * \return the map
	 *
	 */
	const g_map_type & getMap()
	{
		return g_map;
	}

725
	/*! \brief Constructor
726 727
	 *
	 * The periodicity is given by the grid b_g
728
	 *
729
	 * \param stencil maximum extension of the stencil on each directions
730
	 * \param domain the domain
731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749
	 * \param b_g object grid that will store the solution
	 *
	 */
	FDScheme(const Ghost<Sys_eqs::dims,long int> & stencil,
			 const Box<Sys_eqs::dims,typename Sys_eqs::stype> & domain,
			 const typename Sys_eqs::b_grid & b_g)
	:pd({0,0,0},{0,0,0}),gs(b_g.getGridInfoVoid()),g_map(b_g,stencil,pd),row(0),row_b(0)
	{
		Initialize(domain);
	}

	/*! \brief Constructor
	 *
	 * The periodicity is given by the grid b_g
	 *
	 * \param pd Padding, how many points out of boundary are present
	 * \param stencil maximum extension of the stencil on each directions
	 * \param domain the domain
	 * \param b_g object grid that will store the solution
750 751
	 *
	 */
752 753 754 755
	FDScheme(Padding<Sys_eqs::dims> & pd,
			 const Ghost<Sys_eqs::dims,long int> & stencil,
			 const Box<Sys_eqs::dims,typename Sys_eqs::stype> & domain,
			 const typename Sys_eqs::b_grid & b_g)
756
	:pd(pd),gs(b_g.getGridInfoVoid()),g_map(b_g,stencil,pd),row(0),row_b(0)
757
	{
758
		Initialize(domain);
759 760
	}

incardon's avatar
incardon committed
761 762
	/*! \brief Impose an operator
	 *
763 764 765 766 767
	 * This function impose an operator on a box region to produce the system
	 *
	 * Ax = b
	 *
	 * ## Stokes equation, lid driven cavity with one splipping wall
768
	 * \snippet eq_unit_test.hpp lid-driven cavity 2D
incardon's avatar
incardon committed
769
	 *
770 771 772 773 774
	 * \param op Operator to impose (A term)
	 * \param num right hand side of the term (b term)
	 * \param id Equation id in the system that we are imposing
	 * \param start starting point of the box
	 * \param stop stop point of the box
775
	 * \param skip_first skip the first point
incardon's avatar
incardon committed
776 777
	 *
	 */
778 779 780 781 782 783
	template<typename T> void impose(const T & op,
			                         typename Sys_eqs::stype num,
									 long int id,
									 const long int (& start)[Sys_eqs::dims],
									 const long int (& stop)[Sys_eqs::dims],
									 bool skip_first = false)
incardon's avatar
incardon committed
784
	{
785 786
		grid_key_dx<Sys_eqs::dims> start_k;
		grid_key_dx<Sys_eqs::dims> stop_k;
787

788 789 790 791 792 793 794 795 796
        bool increment = false;
        if (skip_first == true)
        {
                start_k = grid_key_dx<Sys_eqs::dims>(start);
                stop_k = grid_key_dx<Sys_eqs::dims>(start);

                auto it = g_map.getSubDomainIterator(start_k,stop_k);

                if (it.isNext() == true)
Pietro Incardona's avatar
Pietro Incardona committed
797
                        increment = true;
798 799 800 801 802 803 804 805 806 807 808
        }

        // add padding to start and stop
        start_k = grid_key_dx<Sys_eqs::dims>(start);
        stop_k = grid_key_dx<Sys_eqs::dims>(stop);

        auto it = g_map.getSubDomainIterator(start_k,stop_k);

        if (increment == true)
                ++it;

809 810
        constant_b b(num);

incardon's avatar
incardon committed
811
        impose_git(op,b,id,it);
812 813 814 815 816 817 818 819 820

	}

	/*! \brief Impose an operator
	 *
	 * This function impose an operator on a particular grid region to produce the system
	 *
	 * Ax = b
	 *
821 822
	 * ## Stokes equation 2D, lid driven cavity with one splipping wall
	 * \snippet eq_unit_test.hpp Copy the solution to grid
823 824 825 826 827 828 829
	 *
	 * \param op Operator to impose (A term)
	 * \param num right hand side of the term (b term)
	 * \param id Equation id in the system that we are imposing
	 * \param it_d iterator that define where you want to impose
	 *
	 */
incardon's avatar
incardon committed
830
	template<typename T> void impose_dit(const T & op ,
831 832 833
									 typename Sys_eqs::stype num,
									 long int id ,
									 grid_dist_iterator_sub<Sys_eqs::dims,typename g_map_type::d_grid> it_d)
834
	{
835
		constant_b b(num);
836

incardon's avatar
incardon committed
837
		impose_dit(op,b,id,it_d);
838
	}
839

840 841 842 843 844 845 846 847 848 849 850 851 852 853 854
	/*! \brief Impose an operator
	 *
	 * This function impose an operator on a particular grid region to produce the system
	 *
	 * Ax = b
	 *
	 * ## Stokes equation 2D, lid driven cavity with one splipping wall
	 * \snippet eq_unit_test.hpp Copy the solution to grid
	 *
	 * \param op Operator to impose (A term)
	 * \param num right hand side of the term (b term)
	 * \param id Equation id in the system that we are imposing
	 * \param it_d iterator that define where you want to impose
	 *
	 */
incardon's avatar
incardon committed
855
	template<unsigned int prp, typename T, typename b_term, typename iterator> void impose_dit(const T & op ,
856
									 b_term & b_t,
857 858
									 const iterator & it_d,
									 long int id = 0)
859 860
	{
		grid_b<b_term,prp> b(b_t);
incardon's avatar
incardon committed
861

incardon's avatar
incardon committed
862
		impose_dit(op,b,id,it_d);
incardon's avatar
incardon committed
863
	}
incardon's avatar
incardon committed
864

865
	//! type of the sparse matrix
866 867
	typename Sys_eqs::SparseMatrix_type A;

incardon's avatar
incardon committed
868 869
	/*! \brief produce the Matrix
	 *
870
	 *  \return the Sparse matrix produced
incardon's avatar
incardon committed
871 872
	 *
	 */
873 874 875 876 877
	typename Sys_eqs::SparseMatrix_type & getA()
	{
#ifdef SE_CLASS1
		consistency();
#endif
878 879 880
		A.resize(g_map.size()*Sys_eqs::nvar,g_map.size()*Sys_eqs::nvar,
				 g_map.getLocalDomainSize()*Sys_eqs::nvar,
				 g_map.getLocalDomainSize()*Sys_eqs::nvar);
881 882 883 884 885 886 887

		return A;

	}

	/*! \brief produce the B vector
	 *
888
	 *  \return the vector produced
889 890 891
	 *
	 */
	typename Sys_eqs::Vector_type & getB()
incardon's avatar
incardon committed
892
	{
893 894 895
#ifdef SE_CLASS1
		consistency();
#endif
incardon's avatar
incardon committed
896

897
		return b;
incardon's avatar
incardon committed
898
	}
899 900 901 902 903 904 905 906 907 908

	/*! \brief Copy the vector into the grid
	 *
	 * ## Copy the solution into the grid
	 * \snippet eq_unit_test.hpp Copy the solution to grid
	 *
	 * \tparam scheme Discretization scheme
	 * \tparam Grid_dst type of the target grid
	 * \tparam pos target properties
	 *
909
	 * \param v Vector that contain the solution of the system
910 911 912 913 914 915 916 917 918 919 920 921 922 923 924
	 * \param start point
	 * \param stop point
	 * \param g_dst Destination grid
	 *
	 */
	template<unsigned int ... pos, typename Vct, typename Grid_dst> void copy(Vct & v,const long int (& start)[Sys_eqs_typ::dims], const long int (& stop)[Sys_eqs_typ::dims], Grid_dst & g_dst)
	{
		if (is_grid_staggered<Sys_eqs>::value())
		{
			if (g_dst.is_staggered() == true)
				copy_staggered<Vct,Grid_dst,pos...>(v,g_dst);
			else
			{
				// Create a temporal staggered grid and copy the data there
				auto & g_map = this->getMap();
925 926 927 928 929 930 931 932 933 934 935

				// Convert the ghost in grid units

				Ghost<Grid_dst::dims,long int> g_int;
				for (size_t i = 0 ; i < Grid_dst::dims ; i++)
				{
					g_int.setLow(i,g_map.getDecomposition().getGhost().getLow(i) / g_map.spacing(i));
					g_int.setHigh(i,g_map.getDecomposition().getGhost().getHigh(i) / g_map.spacing(i));
				}

				staggered_grid_dist<Grid_dst::dims,typename Grid_dst::stype,typename Grid_dst::value_type,typename Grid_dst::decomposition::extended_type, typename Grid_dst::memory_type, typename Grid_dst::device_grid_type> stg(g_dst,g_int,this->getPadding());
936 937 938 939 940 941 942 943
				stg.setDefaultStagPosition();
				copy_staggered<Vct,decltype(stg),pos...>(v,stg);

				// sync the ghost and interpolate to the normal grid
				stg.template ghost_get<pos...>();
				stg.template to_normal<Grid_dst,pos...>(g_dst,this->getPadding(),start,stop);
			}
		}
944 945 946 947
		else
		{
			copy_normal<Vct,Grid_dst,pos...>(v,g_dst);
		}
948
	}
949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975

	/*! \brief Copy the vector into the grid
	 *
	 * ## Copy the solution into the grid
	 * \snippet eq_unit_test.hpp Copy the solution to grid
	 *
	 * \tparam scheme Discretization scheme
	 * \tparam Grid_dst type of the target grid
	 * \tparam pos target properties
	 *
	 * \param v Vector that contain the solution of the system
	 * \param g_dst Destination grid
	 *
	 */
	template<unsigned int ... pos, typename Vct, typename Grid_dst> void copy(Vct & v, Grid_dst & g_dst)
	{
		long int start[Sys_eqs_typ::dims];
		long int stop[Sys_eqs_typ::dims];

		for (size_t i = 0 ; i < Sys_eqs_typ::dims ; i++)
		{
			start[i] = 0;
			stop[i] = g_dst.size(i);
		}

		this->copy<pos...>(v,start,stop,g_dst);
	}
incardon's avatar
incardon committed
976 977 978 979 980 981 982
};

#define EQS_FIELDS 0
#define EQS_SPACE 1


#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_HPP_ */