Some small refactorization

src/FiniteDifference/FDScheme.hpp

@@ -18,13 +18,13 @@
 
 /*! \brief Finite Differences
  *
- * This class is able to discreatize 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 contiguos range on grid points without
- * holes. In order to ensure this, each processor produce a contiguos local labelling of its local
+ * 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
  * 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
- * labelling. An example is shown in the figures down, here we have
- * a grid 8x6 divided across two processor each processor label locally its grid points
+ * 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
  *
  * \verbatim
  *
@@ -266,6 +266,56 @@ private:
 		}
 	}
 
+	/*! \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;
+		}
+	}
+
 public:
 
 	/*! \brief set the staggered position for each property
@@ -387,6 +437,7 @@ public:
 	 * Ax = b
 	 *
 	 * ## Stokes equation, lid driven cavity with one splipping wall
+	 * \snippet eq_unit_test.hpp lid-driven cavity 2D
 	 *
 	 * \param op Operator to impose (A term)
 	 * \param num right hand side of the term (b term)
@@ -412,7 +463,8 @@ public:
 	 *
 	 * Ax = b
 	 *
-	 * ## Stokes equation, lid driven cavity with one splipping wall
+	 * ## 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)
@@ -513,6 +565,42 @@ public:
 
 		return b;
 	}
+
+	/*! \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 scheme Discretization scheme
+	 * \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();
+				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_map.getDecomposition().getGhost(),this->getPadding());
+				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);
+			}
+		}
+	}
 };
 
 #define EQS_FIELDS 0

src/FiniteDifference/eq_unit_test.hpp

@@ -39,18 +39,18 @@ struct lid_nn
 	typedef float stype;
 
 	// type of base grid, it is the distributed grid that will store the result
-	// Note the first property is a 2D vector (velocity), the second is a scalar
+	// Note the first property is a 2D vector (velocity), the second is a scalar (Pressure)
 	typedef grid_dist_id<2,float,aggregate<float[2],float>,CartDecomposition<2,float>> b_grid;
 
 	// type of SparseMatrix, for the linear system, this parameter is bounded by the solver
-	// that you are using
+	// that you are using, in case of umfpack it is the only possible choice
 	typedef SparseMatrix<double,int> SparseMatrix_type;
 
 	// type of Vector for the linear system, this parameter is bounded by the solver
-	// that you are using
+	// that you are using, in case of umfpack it is the only possible choice
 	typedef Vector<double> Vector_type;
 
-	// Define that the underline grid where we discretize the operators is staggered
+	// Define that the underline grid where we discretize the system of equation is staggered
 	static const int grid_type = STAGGERED_GRID;
 };
 
@@ -142,8 +142,15 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
 {
 	Vcluster & v_cl = *global_v_cluster;
 
+	if (v_cl.getProcessingUnits() > 3)
+		return;
+
 	//! [lid-driven cavity 2D]
 
+	// velocity in the grid is the property 0, pressure is the property 1
+	constexpr int velocity = 0;
+	constexpr int pressure = 1;
+
 	// Domain, a rectangle
 	Box<2,float> domain({0.0,0.0},{3.0,1.0});
 
@@ -168,7 +175,7 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
 	// Distributed grid that store the solution
 	grid_dist_id<2,float,aggregate<float[2],float>,CartDecomposition<2,float>> g_dist(szu,domain,g);
 
-	// Ghost stencil
+	// It is the maximum extension of the stencil
 	Ghost<2,long int> stencil_max(1);
 
 	// Finite difference scheme
@@ -219,11 +226,14 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
 
 	auto x = umfpack_solver<double>::solve(fd.getA(),fd.getB());
 
-	// Copy the solution to grid
-	x.copy<FDScheme<lid_nn>,decltype(g_dist),0,1>(fd,{0,0},{sz[0]-1,sz[1]-1},g_dist);
-
 	//! [lid-driven cavity 2D]
 
+	//! [Copy the solution to grid]
+
+	fd.copy<velocity,pressure>(x,{0,0},{sz[0]-1,sz[1]-1},g_dist);
+
+	//! [Copy the solution to grid]
+
 	g_dist.write("lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()));
 
 	if (v_cl.getProcessUnitID() == 0)

src/FiniteDifference/eq_unit_test_3d.hpp

@@ -132,6 +132,9 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
 {
 	Vcluster & v_cl = *global_v_cluster;
 
+	if (v_cl.getProcessingUnits() > 3)
+		return;
+
 	// Domain
 	Box<3,float> domain({0.0,0.0,0.0},{3.0,1.0,1.0});
 
@@ -231,7 +234,7 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
 	auto x = umfpack_solver<double>::solve(fd.getA(),fd.getB());
 
 	// Bring the solution to grid
-	x.copy<FDScheme<lid_nn_3d>,decltype(g_dist),velocity,pressure>(fd,{0,0},{sz[0]-1,sz[1]-1,sz[2]-1},g_dist);
+	fd.copy<velocity,pressure>(x,{0,0},{sz[0]-1,sz[1]-1,sz[2]-1},g_dist);
 
 	g_dist.write("lid_driven_cavity_3d_p" + std::to_string(v_cl.getProcessingUnits()));
 

src/Vector/Vector_eigen.hpp

@@ -88,53 +88,6 @@ class Vector<T,Eigen::Matrix<T, Eigen::Dynamic, 1>>
 	//size of each chunk
 	mutable openfpm::vector<size_t> sz;
 
-	/*! \brief Copy in a staggered grid
-	 *
-	 *
-	 */
-	template<typename scheme, typename Grid_dst ,unsigned int ... pos> void copy_staggered_to_staggered(scheme & sc, Grid_dst & g_dst)
-	{
-		// get the map
-		const auto & g_map = sc.getMap();
-
-		// check that g_dst is staggered
-		if (is_grid_staggered<Grid_dst>::value() == 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<typename scheme::Sys_eqs_typ,Grid_dst,typename std::remove_reference<decltype(*this)>::type> cp(key_dst,g_dst,*this,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 Here we collect the full matrix on master
 	 *
@@ -350,42 +303,6 @@ public:
 		return v;
 	}
 
-	/*! \brief Copy the vector into the grid
-	 *
-	 * ## Copy from the vector to the destination grid
-	 * \snippet eq_unit_tests.hpp
-	 *
-	 * \tparam scheme Discretization scheme
-	 * \tparam Grid_dst type of the target grid
-	 * \tparam pos target properties
-	 *
-	 * \param scheme Discretization scheme
-	 * \param start point
-	 * \param stop point
-	 * \param g_dst Destination grid
-	 *
-	 */
-	template<typename scheme, typename Grid_dst ,unsigned int ... pos> void copy(scheme & sc,const long int (& start)[scheme::Sys_eqs_typ::dims], const long int (& stop)[scheme::Sys_eqs_typ::dims], Grid_dst & g_dst)
-	{
-		if (is_grid_staggered<typename scheme::Sys_eqs_typ>::value())
-		{
-			if (g_dst.is_staggered() == true)
-				copy_staggered_to_staggered<scheme,Grid_dst,pos...>(sc,g_dst);
-			else
-			{
-				// Create a temporal staggered grid and copy the data there
-				auto & g_map = sc.getMap();
-				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_map.getDecomposition().getGhost(),sc.getPadding());
-				stg.setDefaultStagPosition();
-				copy_staggered_to_staggered<scheme,decltype(stg),pos...>(sc,stg);
-
-				// sync the ghost and interpolate to the normal grid
-				stg.template ghost_get<pos...>();
-				stg.template to_normal<Grid_dst,pos...>(g_dst,sc.getPadding(),start,stop);
-			}
-		}
-	}
-
 	/*! \brief Scatter the vector information to the other processors
 	 *
 	 * Eigen does not have a real parallel vector, so in order to work we have to scatter