Improving the documentation of numeric

src/FiniteDifference/Average.hpp

@@ -42,18 +42,28 @@ class Avg
 
 	/*! \brief Calculate the position where the derivative is calculated
 	 *
-	 * In case on non staggered case this function just return pos, in case of staggered,
-	 *  it calculate where the operator is calculated on a staggered grid
+	 * In case of non staggered case this function just return a null grid_key, in case of staggered,
+	 *  it calculate how the operator shift the calculation in the cell
+	 *
+	 * \param position where we are calculating the derivative
+	 * \param gs Grid info
+	 * \param s_pos staggered position of the properties
 	 *
 	 */
-	inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs)
+	inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, const openfpm::vector<comb<Sys_eqs::dims>> & s_pos = openfpm::vector<comb<Sys_eqs::dims>>())
 	{
 		std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
 	}
 };
 
-/*! \brief Average on direction i
+/*! \brief Central average scheme on direction i
+ *
+ * \verbatim
+ *
+ *  +0.5     +0.5
+ *   *---+---*
  *
+ * \endverbatim
  *
  */
 template<unsigned int d, typename arg, typename Sys_eqs>
@@ -61,8 +71,20 @@ class Avg<d,arg,Sys_eqs,CENTRAL>
 {
 	public:
 
-	/*! \brief fill the row
+	/*! \brief Calculate which colums of the Matrix are non zero
+	 *
+	 * stub_or_real it is just for change the argument type when testing, in normal
+	 * conditions it is a distributed map
+	 *
+	 * \param g_map It is the map explained in FDScheme
+	 * \param k_map position where the average is calculated
+	 * \param gs Grid info
+	 * \param cols non-zero colums calculated by the function
+	 * \param coeff coefficent (constant in front of the derivative)
 	 *
+	 * ### Example
+	 *
+	 * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
 	 *
 	 */
 	inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
@@ -87,43 +109,44 @@ class Avg<d,arg,Sys_eqs,CENTRAL>
 	}
 
 
-	/*! \brief Calculate the position where the derivative is calculated
+	/*! \brief Calculate the position where the average is calculated
 	 *
-	 * In case on non staggered case this function just return pos, in case of staggered,
-	 *  it calculate where the operator is calculated on a staggered grid
+	 * It follow the same concept of central derivative
 	 *
-	 *  \param pos from the position
-	 *  \param fld Field we are deriving, if not provided the function just return pos
-	 *  \param s_pos position of the properties in the staggered grid
+	 * \param position where we are calculating the derivative
+	 * \param gs Grid info
+	 * \param s_pos staggered position of the properties
 	 *
 	 */
-	inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, long int fld = -1, const openfpm::vector<comb<Sys_eqs::dims>> & s_pos = openfpm::vector<comb<Sys_eqs::dims>>())
+	inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, const comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
 	{
+		auto arg_pos = arg::position(pos,gs,s_pos);
 		if (is_grid_staggered<Sys_eqs>::value())
 		{
-			if (fld == -1)
-				return pos;
-
-			if (s_pos[fld][d] == 1)
+			if (arg_pos.get(d) == -1)
 			{
-				grid_key_dx<Sys_eqs::dims> ret = pos;
-				ret.set_d(d,0);
-				return pos;
+				arg_pos.set_d(d,0);
+				return arg_pos;
 			}
 			else
 			{
-				grid_key_dx<Sys_eqs::dims> ret = pos;
-				ret.set_d(d,1);
-				return pos;
+				arg_pos.set_d(d,-1);
+				return arg_pos;
 			}
 		}
 
-		return pos;
+		return arg_pos;
 	}
 };
 
-/*! \brief Average FORWARD on direction i
+/*! \brief FORWARD average on direction i
+ *
+ * \verbatim
+ *
+ *  +0.5    0.5
+ *    +------*
  *
+ * \endverbatim
  *
  */
 template<unsigned int d, typename arg, typename Sys_eqs>
@@ -131,8 +154,20 @@ class Avg<d,arg,Sys_eqs,FORWARD>
 {
 	public:
 
-	/*! \brief fill the row
+	/*! \brief Calculate which colums of the Matrix are non zero
 	 *
+	 * stub_or_real it is just for change the argument type when testing, in normal
+	 * conditions it is a distributed map
+	 *
+	 * \param g_map It is the map explained in FDScheme
+	 * \param k_map position where the average is calculated
+	 * \param gs Grid info
+	 * \param cols non-zero colums calculated by the function
+	 * \param coeff coefficent (constant in front of the derivative)
+	 *
+	 * ### Example
+	 *
+	 * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
 	 *
 	 */
 	inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
@@ -148,24 +183,30 @@ class Avg<d,arg,Sys_eqs,FORWARD>
 	}
 
 
-	/*! \brief Calculate the position where the derivative is calculated
+	/*! \brief Calculate the position in the cell where the average is calculated
 	 *
-	 * In case on non staggered case this function just return pos, in case of staggered,
-	 *  it calculate where the operator is calculated on a staggered grid
+	 * In case of non staggered case this function just return a null grid_key, in case of staggered,
+	 * the FORWARD scheme return the position of the staggered property
 	 *
-	 *  \param pos from the position
-	 *  \param fld Field we are deriving, if not provided the function just return pos
-	 *  \param s_pos position of the properties in the staggered grid
+	 * \param position where we are calculating the derivative
+	 * \param gs Grid info
+	 * \param s_pos staggered position of the properties
 	 *
 	 */
-	inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, long int fld = -1, const openfpm::vector<comb<Sys_eqs::dims>> & s_pos = openfpm::vector<comb<Sys_eqs::dims>>())
+	inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, const comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
 	{
-		return pos;
+		return arg::position(pos,gs,s_pos);
 	}
 };
 
-/*! \brief Average BACKWARD on direction i
+/*! \brief First order BACKWARD derivative on direction i
+ *
+ * \verbatim
  *
+ *  +0.5    0.5
+ *    *------+
+ *
+ * \endverbatim
  *
  */
 template<unsigned int d, typename arg, typename Sys_eqs>
@@ -173,8 +214,20 @@ class Avg<d,arg,Sys_eqs,BACKWARD>
 {
 	public:
 
-	/*! \brief fill the row
+	/*! \brief Calculate which colums of the Matrix are non zero
+	 *
+	 * stub_or_real it is just for change the argument type when testing, in normal
+	 * conditions it is a distributed map
+	 *
+	 * \param g_map It is the map explained in FDScheme
+	 * \param k_map position where the average is calculated
+	 * \param gs Grid info
+	 * \param cols non-zero colums calculated by the function
+	 * \param coeff coefficent (constant in front of the derivative)
 	 *
+	 * ### Example
+	 *
+	 * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
 	 *
 	 */
 	inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims], std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
@@ -189,19 +242,19 @@ class Avg<d,arg,Sys_eqs,BACKWARD>
 	}
 
 
-	/*! \brief Calculate the position where the derivative is calculated
+	/*! \brief Calculate the position in the cell where the average is calculated
 	 *
-	 * In case on non staggered case this function just return pos, in case of staggered,
-	 *  it calculate where the operator is calculated on a staggered grid
+	 * In case of non staggered case this function just return a null grid_key, in case of staggered,
+	 * the BACKWARD scheme return the position of the staggered property
 	 *
-	 *  \param pos from the position
-	 *  \param fld Field we are deriving, if not provided the function just return pos
-	 *  \param s_pos position of the properties in the staggered grid
+	 * \param position where we are calculating the derivative
+	 * \param gs Grid info
+	 * \param s_pos staggered position of the properties
 	 *
 	 */
-	inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, long int fld = -1, const openfpm::vector<comb<Sys_eqs::dims>> & s_pos = openfpm::vector<comb<Sys_eqs::dims>>())
+	inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, const comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
 	{
-		return pos;
+		return arg::position(pos,gs,s_pos);
 	}
 };
 

src/FiniteDifference/FDScheme.hpp

@@ -14,10 +14,11 @@
 #include "eq.hpp"
 #include "data_type/scalar.hpp"
 #include "NN/CellList/CellDecomposer.hpp"
+#include "Grid/staggered_dist_grid_util.hpp"
 
 /*! \brief Finite Differences
  *
- * This class is able to discreatize on a Matrix any operator. In order to create a consistent
+ * 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
  * points. Each processor also add an offset equal to the number of local
@@ -67,7 +68,53 @@
  *
  * \endverbatim
  *
- * \tparam dim Dimensionality of the finite differences scheme
+ * \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
  *
  */
 
@@ -111,12 +158,52 @@ private:
 	// Grid points that has each processor
 	openfpm::vector<size_t> pnt;
 
+	// Staggered position for each property
+	comb<Sys_eqs::dims> s_pos[Sys_eqs::nvar];
+
 	// 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
 	size_t s_pnt;
 
+	struct key_and_eq
+	{
+		grid_key_dx<Sys_eqs::dims> key;
+		size_t eq;
+	};
+
+	/*! \brief From the row Matrix position to the spatial position
+	 *
+	 * \param row Matrix
+	 *
+	 * \return spatial position
+	 *
+	 */
+	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());
+				ke.key -= pd.getKP1();
+				return ke;
+			}
+
+			++it;
+		}
+		std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " the row does not map to any position" << "\n";
+
+		return ke;
+	}
+
 	/* \brief calculate the mapping grid size with padding
 	 *
 	 * \param gs original grid size
@@ -165,7 +252,10 @@ private:
 		for (size_t i = 0 ; i < nz_rows.size() ; i++)
 		{
 			if (nz_rows.get(i) == false)
-				std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix row " << i <<  " is not filled\n";
+			{
+				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";
+			}
 		}
 
 		// all the colums must have a non zero element
@@ -176,21 +266,42 @@ private:
 		}
 	}
 
-	/* \brief Convert discrete ghost into continous ghost
+public:
+
+	/*! \brief set the staggered position for each property
 	 *
-	 * \param Ghost
+	 * \param vector containing the staggered position for each property
 	 *
 	 */
-	inline Ghost<Sys_eqs::dims,typename Sys_eqs::stype> convert_dg_cg(const Ghost<Sys_eqs::dims,typename Sys_eqs::stype> & g)
+	void setStagPos(comb<Sys_eqs::dims> (& sp)[Sys_eqs::nvar])
 	{
-		return g_map_type::convert_ghost(g,g_map.getCellDecomposer());
+		for (size_t i = 0 ; i < Sys_eqs::nvar ; i++)
+			s_pos[i] = sp[i];
 	}
 
-public:
+	/*! \brief compute the staggered position for each property
+	 *
+	 * This is compute from the value_type stored by Sys_eqs::b_grid::value_type
+	 *
+	 * ### Example
+	 *
+	 * \snippet eq_unit_test.hpp Compute staggered properties
+	 *
+	 */
+	void computeStag()
+	{
+		typedef typename Sys_eqs::b_grid::value_type prp_type;
 
-	/*! \brief Get the grid padding
+		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
 	 *
-	 * \return the grid padding
+	 * \return the padding specified
 	 *
 	 */
 	const Padding<Sys_eqs::dims> & getPadding()
@@ -199,6 +310,8 @@ public:
 	}
 
 	/*! \brief Return the map between the grid index position and the position in the distributed vector
+	 *
+	 * It is the map explained in the intro of the FDScheme
 	 *
 	 * \return the map
 	 *
@@ -210,13 +323,17 @@ public:
 
 	/*! \brief Constructor
 	 *
-	 * \param pd Padding
+	 * \param pd Padding, how many points out of boundary are present
+	 * \param domain extension of the domain
 	 * \param gs grid infos where Finite differences work
+	 * \param dec Decomposition of the domain
 	 *
 	 */
-	FDScheme(Padding<Sys_eqs::dims> & pd, const Box<Sys_eqs::dims,typename Sys_eqs::stype> & domain, const grid_sm<Sys_eqs::dims,void> & gs, const typename Sys_eqs::b_grid::decomposition & dec, Vcluster & v_cl)
+	FDScheme(Padding<Sys_eqs::dims> & pd, const Box<Sys_eqs::dims,typename Sys_eqs::stype> & domain, const grid_sm<Sys_eqs::dims,void> & gs, const typename Sys_eqs::b_grid::decomposition & dec)
 	:pd(pd),gs(gs),g_map(dec,padding(gs.getSize(),pd),domain,Ghost<Sys_eqs::dims,typename Sys_eqs::stype>(1)),row(0),row_b(0)
 	{
+		Vcluster & v_cl = dec.getVC();
+
 		// Calculate the size of the local domain
 		size_t sz = g_map.getLocalDomainSize();
 
@@ -337,17 +454,22 @@ public:
 			for ( auto it = cols.begin(); it != cols.end(); ++it )
 			{
 				trpl.add();
-				trpl.last().row() = Sys_eqs::nvar * gs.LinId(gkey) + id;
+				trpl.last().row() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
 				trpl.last().col() = it->first;
 				trpl.last().value() = it->second;
 
-				std::cout << "(" << trpl.last().row() << "," << trpl.last().col() << "," << trpl.last().value() << ")" << "\n";
+//				std::cout << "(" << trpl.last().row() << "," << trpl.last().col() << "," << trpl.last().value() << ")" << "\n";
 			}
 
-			b.get(Sys_eqs::nvar * gs.LinId(gkey) + id) = num;
+			b.get(g_map.template get<0>(key)*Sys_eqs::nvar + id) = num;
 
 			cols.clear();
-			std::cout << "\n";
+//			std::cout << "\n";
+
+			// if SE_CLASS1 is defined check the position
+#ifdef SE_CLASS1
+			T::position(gkey,gs,s_pos);
+#endif
 
 			++row;
 			++row_b;
@@ -359,7 +481,7 @@ public:
 
 	/*! \brief produce the Matrix
 	 *
-	 *  \tparam Syst_eq System of equations, or a single equation
+	 *  \return the Sparse matrix produced
 	 *
 	 */
 	typename Sys_eqs::SparseMatrix_type & getA()
@@ -378,7 +500,7 @@ public:
 
 	/*! \brief produce the B vector
 	 *
-	 *  \tparam Syst_eq System of equations, or a single equation
+	 *  \return the vector produced
 	 *
 	 */
 	typename Sys_eqs::Vector_type & getB()

src/FiniteDifference/FDScheme_unit_tests.hpp

@@ -121,6 +121,8 @@ BOOST_AUTO_TEST_SUITE( fd_test )
 
 BOOST_AUTO_TEST_CASE( der_central_non_periodic)
 {
+	//! [Usage of stencil derivative]
+
 	// grid size
 	size_t sz[2]={16,16};
 
@@ -155,6 +157,8 @@ BOOST_AUTO_TEST_CASE( der_central_non_periodic)
 	BOOST_REQUIRE_EQUAL(cols_y[17+16],1/spacing[1]/2.0);
 	BOOST_REQUIRE_EQUAL(cols_y[17-16],-1/spacing[1]/2.0);
 
+	//! [Usage of stencil derivative]
+
 	// filled colums
 
 	std::unordered_map<long int,float> cols_xx;
@@ -178,19 +182,19 @@ BOOST_AUTO_TEST_CASE( der_central_non_periodic)
 	BOOST_REQUIRE_EQUAL(cols_xx[34],-2/spacing[0]/spacing[0]/2.0/2.0);
 	BOOST_REQUIRE_EQUAL(cols_xx[36],1/spacing[0]/spacing[0]/2.0/2.0);
 
-	BOOST_REQUIRE_EQUAL(cols_xy[17],1/spacing[0]/spacing[1]);
-	BOOST_REQUIRE_EQUAL(cols_xy[19],-1/spacing[0]/spacing[1]);
-	BOOST_REQUIRE_EQUAL(cols_xy[49],-1/spacing[0]/spacing[1]);
-	BOOST_REQUIRE_EQUAL(cols_xy[51],1/spacing[0]/spacing[1]);
+	BOOST_REQUIRE_EQUAL(cols_xy[17],1/spacing[0]/spacing[1]/2.0/2.0);
+	BOOST_REQUIRE_EQUAL(cols_xy[19],-1/spacing[0]/spacing[1]/2.0/2.0);
+	BOOST_REQUIRE_EQUAL(cols_xy[49],-1/spacing[0]/spacing[1]/2.0/2.0);
+	BOOST_REQUIRE_EQUAL(cols_xy[51],1/spacing[0]/spacing[1]/2.0/2.0);
 
-	BOOST_REQUIRE_EQUAL(cols_yx[17],1/spacing[0]/spacing[1]);
-	BOOST_REQUIRE_EQUAL(cols_yx[19],-1/spacing[0]/spacing[1]);
-	BOOST_REQUIRE_EQUAL(cols_yx[49],-1/spacing[0]/spacing[1]);
-	BOOST_REQUIRE_EQUAL(cols_xy[51],1/spacing[0]/spacing[1]);
+	BOOST_REQUIRE_EQUAL(cols_yx[17],1/spacing[0]/spacing[1]/2.0/2.0);
+	BOOST_REQUIRE_EQUAL(cols_yx[19],-1/spacing[0]/spacing[1]/2.0/2.0);
+	BOOST_REQUIRE_EQUAL(cols_yx[49],-1/spacing[0]/spacing[1]/2.0/2.0);
+	BOOST_REQUIRE_EQUAL(cols_xy[51],1/spacing[0]/spacing[1]/2.0/2.0);
 
-	BOOST_REQUIRE_EQUAL(cols_yy[2],1/spacing[1]/spacing[1]);
-	BOOST_REQUIRE_EQUAL(cols_yy[34],-2/spacing[1]/spacing[1]);
-	BOOST_REQUIRE_EQUAL(cols_yy[66],1/spacing[1]/spacing[1]);
+	BOOST_REQUIRE_EQUAL(cols_yy[2],1/spacing[1]/spacing[1]/2.0/2.0);
+	BOOST_REQUIRE_EQUAL(cols_yy[34],-2/spacing[1]/spacing[1]/2.0/2.0);
+	BOOST_REQUIRE_EQUAL(cols_yy[66],1/spacing[1]/spacing[1]/2.0/2.0);
 
 	// Non periodic with one sided
 
@@ -462,6 +466,8 @@ BOOST_AUTO_TEST_CASE( avg_central_staggered_non_periodic)
 
 BOOST_AUTO_TEST_CASE( lap_periodic)
 {
+	//! [Laplacian usage]
+
 	// grid size
 	size_t sz[2]={16,16};
 
@@ -493,6 +499,8 @@ BOOST_AUTO_TEST_CASE( lap_periodic)
 
 	BOOST_REQUIRE_EQUAL(cols[17],-2/spacing[0]/spacing[0] - 2/spacing[1]/spacing[1]);
 
+	//! [Laplacian usage]
+
 	cols.clear();
 
 	Lap<Field<V,sys_pp>,sys_pp>::value(g_map,key1414,ginfo,spacing,cols,1);
@@ -512,6 +520,8 @@ BOOST_AUTO_TEST_CASE( lap_periodic)
 
 BOOST_AUTO_TEST_CASE( sum_periodic)
 {
+	//! [sum example]
+
 	// grid size
 	size_t sz[2]={16,16};
 
@@ -538,6 +548,8 @@ BOOST_AUTO_TEST_CASE( sum_periodic)
 
 	BOOST_REQUIRE_EQUAL(cols[17],2);
 
+	//! [sum example]
+
 	cols.clear();
 
 	sum<Field<V,sys_pp>, Field<V,sys_pp> , Field<V,sys_pp> ,sys_pp>::value(g_map,key11,ginfo,spacing,cols,1);
@@ -548,11 +560,15 @@ BOOST_AUTO_TEST_CASE( sum_periodic)
 
 	cols.clear();
 
+	//! [minus example]
+
 	sum<Field<V,sys_pp>, Field<V,sys_pp> , minus<Field<V,sys_pp>,sys_pp> ,sys_pp>::value(g_map,key11,ginfo,spacing,cols,1);
 
 	BOOST_REQUIRE_EQUAL(cols.size(),1ul);
 
 	BOOST_REQUIRE_EQUAL(cols[17],1ul);
+
+	//! [minus example]
 }
 
 //////////////// Position ////////////////////
@@ -560,73 +576,70 @@ BOOST_AUTO_TEST_CASE( sum_periodic)
 BOOST_AUTO_TEST_CASE( fd_test_use_staggered_position)
 {
 	// grid size
-/*	size_t sz[2]={16,16};
+	size_t sz[2]={16,16};
+
+	// grid_sm
+	grid_sm<2,void> ginfo(sz);
 
 	// Create a derivative row Matrix
-	grid_key_dx<2> key11(1,1);
 	grid_key_dx<2> key00(0,0);
+	grid_key_dx<2> key11(1,1);
 	grid_key_dx<2> key22(2,2);
 	grid_key_dx<2> key1515(15,15);
 
-	openfpm::vector<comb<2>> vx_c;
-	comb<2> cx = {{1,0}};
-	vx_c.add(cx);         // like v_x
-
-	openfpm::vector<comb<2>> vy_c;
-	comb<2> cy = {{0,1}};
-	vy_c.add(cy);         // like v_y
-
-	grid_key_dx<2> key_ret_x = D<x,Field<V,sys_nn>,sys_nn>::position(key11,0,vx_c);
-	grid_key_dx<2> key_ret_y = D<y,Field<V,sys_nn>,sys_nn>::position(key11,0,vx_c);
-
-	BOOST_REQUIRE_EQUAL(key_ret_x.get(0),1);
-	BOOST_REQUIRE_EQUAL(key_ret_x.get(1),1);
-
-	BOOST_REQUIRE_EQUAL(key_ret_y.get(0),0);
-	BOOST_REQUIRE_EQUAL(key_ret_y.get(1),0);
+	comb<2> vx_c[] = {{0,-1}};
+	comb<2> vy_c[] = {{-1,0}};
 
-	key_ret_x = D<x,Field<V,sys_nn>,sys_nn>::position(key11,0,vy_c);
-	key_ret_y = D<y,Field<V,sys_nn>,sys_nn>::position(key11,0,vy_c);
+	grid_key_dx<2> key_ret_vx_x = D<x,Field<V,syss_nn>,syss_nn>::position(key00,ginfo,vx_c);
+	grid_key_dx<2> key_ret_vx_y = D<y,Field<V,syss_nn>,syss_nn>::position(key00,ginfo,vx_c);
+	grid_key_dx<2> key_ret_vy_x = D<x,Field<V,syss_nn>,syss_nn>::position(key00,ginfo,vy_c);
+	grid_key_dx<2> key_ret_vy_y = D<y,Field<V,syss_nn>,syss_nn>::position(key00,ginfo,vy_c);