diff --git a/src/FiniteDifference/Average.hpp b/src/FiniteDifference/Average.hpp
index 25288267eaad6a53e436e3f7c9322e26762a771c..2d578bd592c8540c75cdfb21dfa202d5d885feba 100644
--- a/src/FiniteDifference/Average.hpp
+++ b/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);
}
};
diff --git a/src/FiniteDifference/Derivative.hpp b/src/FiniteDifference/Derivative.hpp
index 3fd1710fc7ef81c528f786fc63e881f57778efdc..dfd295ba165b26e30306cd8a19c5dafae074a667 100644
--- a/src/FiniteDifference/Derivative.hpp
+++ b/src/FiniteDifference/Derivative.hpp
@@ -19,6 +19,7 @@
#include "FiniteDifference/util/common.hpp"
#include "util/util_num.hpp"
+
/*! \brief Derivative second order on h (spacing)
*
* \tparam d on which dimension derive
@@ -29,9 +30,16 @@
template<unsigned int d, typename Field, typename Sys_eqs, unsigned int impl=CENTRAL>
class D
{
- /*! \brief Create the row of the Matrix
+ /*! \brief Calculate which colums of the Matrix are non zero
+ *
+ * \param pos position where the derivative is calculated
+ * \param gs Grid info
+ * \param cols non-zero colums calculated by the function
+ * \param coeff coefficent (constant in front of the derivative)
*
- * \tparam ord
+ * ### Example
+ *
+ * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
*
*/
inline static void value(const grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
@@ -41,18 +49,28 @@ class D
/*! \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 comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
{
std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
}
};
-/*! \brief Derivative on direction i
+/*! \brief Second order central Derivative scheme on direction i
*
+ * \verbatim
+ *
+ * -1 +1
+ * *---+---*
+ *
+ * \endverbatim
*
*/
template<unsigned int d, typename arg, typename Sys_eqs>
@@ -60,8 +78,16 @@ class D<d,arg,Sys_eqs,CENTRAL>
{
public:
- /*! \brief fill the row
+ /*! \brief Calculate which colums of the Matrix are non zero
+ *
+ * \param pos position where the derivative 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)
@@ -88,42 +114,78 @@ class D<d,arg,Sys_eqs,CENTRAL>
/*! \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 on non staggered case this function just return a null grid_key, in case of staggered,
+ * it calculate how the operator shift in the cell
+ *
+ \verbatim
+
+ +--$--+
+ | |
+ # * #
+ | |
+ 0--$--+
+
+ # = velocity(y)
+ $ = velocity(x)
+ * = pressure
+
+ \endverbatim
+ *
+ * Consider this 2D staggered cell and a second order central derivative scheme, this lead to
*
- * \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
+ * \f$ \frac{\partial v_y}{\partial x} \f$ is calculated on position (*), so the function return the grid_key {0,0}
+ *
+ * \f$ \frac{\partial v_y}{\partial y} \f$ is calculated on position (0), so the function return the grid_key {-1,-1}
+ *
+ * \f$ \frac{\partial v_x}{\partial x} \f$ is calculated on position (0), so the function return the grid_key {-1,-1}
+ *
+ * \f$ \frac{\partial v_x}{\partial y} \f$ is calculated on position (*), so the function return the grid_key {0,0}
+ *
+ * \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 Derivative on direction i
+/*! \brief Second order one sided Derivative scheme on direction i
+ *
+ * \verbatim
+ *
+ * -1.5 2.0 -0.5
+ * +------*------*
*
+ * or
+ *
+ * -0.5 2.0 -1.5
+ * *------*------+
+ *
+ * in the bulk
+ *
+ * -1 +1
+ * *---+---*
+ *
+ * \endverbatim
*
*/
template<unsigned int d, typename arg, typename Sys_eqs>
@@ -131,15 +193,23 @@ class D<d,arg,Sys_eqs,CENTRAL_B_ONE_SIDE>
{
public:
- /*! \brief fill the row
+ /*! \brief Calculate which colums of the Matrix are non zero
+ *
+ * \param pos position where the derivative 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
*
*/
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)
{
#ifdef SE_CLASS1
if (Sys_eqs::boundary[d] == PERIODIC)
- std::cerr << __FILE__ << ":" << __LINE__ << " error, it make no sense use one sided derivation with periodic boundary\n";
+ std::cerr << __FILE__ << ":" << __LINE__ << " error, it make no sense use one sided derivation with periodic boundary, please use CENTRAL\n";
#endif
grid_key_dx<Sys_eqs::dims> pos = g_map.getGKey(kmap);
@@ -188,47 +258,64 @@ public:
/*! \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 on non staggered case this function just return a null grid_key, in case of staggered,
+ * it calculate how the operator shift in the cell
*
- */
- inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, long int fld = 0, const openfpm::vector<comb<Sys_eqs::dims>> & s_pos = openfpm::vector<comb<Sys_eqs::dims>>())
- {
- if (is_grid_staggered<Sys_eqs>::value())
- {
- if (fld == -1)
- return pos;
+ \verbatim
- if (s_pos[fld][d] == 1)
- {
- grid_key_dx<Sys_eqs::dims> ret = pos;
- ret.set_d(d,0);
- return pos;
- }
- else
- {
- grid_key_dx<Sys_eqs::dims> ret = pos;
- ret.set_d(d,1);
- return pos;
- }
- }
+ +--$--+
+ | |
+ # * #
+ | |
+ 0--$--+
- return pos;
+ # = velocity(y)
+ $ = velocity(x)
+ * = pressure
+
+ \endverbatim
+ *
+ * In the one side stencil the cell position depend if you are or not at the boundary.
+ * outside the boundary is simply the central scheme, at the boundary it is simply the
+ * staggered position of the property
+ *
+ * \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, const comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
+ {
+ return arg::position(pos,gs,s_pos);
}
};
-/*! \brief Derivative FORWARD on direction i
+/*! \brief First order FORWARD derivative on direction i
+ *
+ * \verbatim
+ *
+ * -1.0 1.0
+ * +------*
+ *
+ * \endverbatim
*
- *g
*/
template<unsigned int d, typename arg, typename Sys_eqs>
class D<d,arg,Sys_eqs,FORWARD>
{
public:
- /*! \brief fill the row
+ /*! \brief Calculate which colums of the Matrix are non zero
+ *
+ * \param pos position where the derivative 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)
@@ -246,31 +333,45 @@ class D<d,arg,Sys_eqs,FORWARD>
/*! \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,
+ * 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 Derivative BACKWARD on direction i
+/*! \brief First order BACKWARD derivative on direction i
+ *
+ * \verbatim
+ *
+ * -1.0 1.0
+ * *------+
+ *
+ * \endverbatim
*
- *g
*/
template<unsigned int d, typename arg, typename Sys_eqs>
class D<d,arg,Sys_eqs,BACKWARD>
{
public:
- /*! \brief fill the row
+ /*! \brief Calculate which colums of the Matrix are non zero
+ *
+ * \param pos position where the derivative 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)
@@ -288,17 +389,17 @@ class D<d,arg,Sys_eqs,BACKWARD>
/*! \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,
+ * 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);
}
};
diff --git a/src/FiniteDifference/FDScheme.hpp b/src/FiniteDifference/FDScheme.hpp
index 7eb338a932ab8c6d811506cfd6dec91aea264fe8..cae23f7b62a087b70132c165dae57debd8f80d43 100644
--- a/src/FiniteDifference/FDScheme.hpp
+++ b/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()
diff --git a/src/FiniteDifference/FDScheme_unit_tests.hpp b/src/FiniteDifference/FDScheme_unit_tests.hpp
index 212d855dcdb4cad825f4392c5b3e996691c21c13..9079bdc4ca0a30098cd7c271e17659417bfab4e3 100644
--- a/src/FiniteDifference/FDScheme_unit_tests.hpp
+++ b/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);
- BOOST_REQUIRE_EQUAL(key_ret_x.get(0),0);
- BOOST_REQUIRE_EQUAL(key_ret_x.get(1),0);
-
- BOOST_REQUIRE_EQUAL(key_ret_y.get(0),1);
- BOOST_REQUIRE_EQUAL(key_ret_y.get(1),1);
+ BOOST_REQUIRE_EQUAL(key_ret_vx_x.get(0),0);
+ BOOST_REQUIRE_EQUAL(key_ret_vx_x.get(1),0);
+ BOOST_REQUIRE_EQUAL(key_ret_vx_y.get(0),-1);
+ BOOST_REQUIRE_EQUAL(key_ret_vx_y.get(1),-1);
+ BOOST_REQUIRE_EQUAL(key_ret_vy_y.get(0),0);
+ BOOST_REQUIRE_EQUAL(key_ret_vy_y.get(1),0);
+ BOOST_REQUIRE_EQUAL(key_ret_vy_x.get(0),-1);
+ BOOST_REQUIRE_EQUAL(key_ret_vy_x.get(1),-1);
// Composed derivative
- grid_key_dx<2> key_ret_xx = D<x,D<x,Field<V,sys_nn>,sys_nn>,sys_nn>::position(key00,0);
- grid_key_dx<2> key_ret_xy = D<x,D<y,Field<V,sys_nn>,sys_nn>,sys_nn>::position(key00,0);
- grid_key_dx<2> key_ret_yx = D<y,D<x,Field<V,sys_nn>,sys_nn>,sys_nn>::position(key00,0);
- grid_key_dx<2> key_ret_yy = D<y,D<y,Field<V,sys_nn>,sys_nn>,sys_nn>::position(key00,0);
+ grid_key_dx<2> key_ret_xx = D<x,D<x,Field<V,syss_nn>,syss_nn>,syss_nn>::position(key00,ginfo,vx_c);
+ grid_key_dx<2> key_ret_xy = D<x,D<y,Field<V,syss_nn>,syss_nn>,syss_nn>::position(key00,ginfo,vx_c);
+ grid_key_dx<2> key_ret_yx = D<y,D<x,Field<V,syss_nn>,syss_nn>,syss_nn>::position(key00,ginfo,vx_c);
+ grid_key_dx<2> key_ret_yy = D<y,D<y,Field<V,syss_nn>,syss_nn>,syss_nn>::position(key00,ginfo,vx_c);
- BOOST_REQUIRE_EQUAL(key_ret_xx.get(0),1);
+ BOOST_REQUIRE_EQUAL(key_ret_xx.get(0),-1);
BOOST_REQUIRE_EQUAL(key_ret_xx.get(1),0);
- BOOST_REQUIRE_EQUAL(key_ret_xy.get(0),1);
- BOOST_REQUIRE_EQUAL(key_ret_xy.get(1),0);
-
BOOST_REQUIRE_EQUAL(key_ret_xy.get(0),0);
- BOOST_REQUIRE_EQUAL(key_ret_xy.get(1),1);
+ BOOST_REQUIRE_EQUAL(key_ret_xy.get(1),-1);
- BOOST_REQUIRE_EQUAL(key_ret_xy.get(0),0);
- BOOST_REQUIRE_EQUAL(key_ret_xy.get(1),1);
+ BOOST_REQUIRE_EQUAL(key_ret_yx.get(0),0);
+ BOOST_REQUIRE_EQUAL(key_ret_yx.get(1),-1);
- ////////////////// Non periodic with one sided
+ BOOST_REQUIRE_EQUAL(key_ret_yy.get(0),-1);
+ BOOST_REQUIRE_EQUAL(key_ret_yy.get(1),0);
- key_ret_x = D<x,Field<V,sys_nn>,sys_nn,CENTRAL_B_ONE_SIDE>::position(key00,0,vx_c);
- key_ret_y = D<y,Field<V,sys_nn>,sys_nn,CENTRAL_B_ONE_SIDE>::position(key00,0,vx_c);
+ ////////////////// Non periodic one sided
- BOOST_REQUIRE_EQUAL(key_ret_x.get(0),0);
- BOOST_REQUIRE_EQUAL(key_ret_x.get(1),1);
+ key_ret_vx_x = D<x,Field<V,syss_nn>,syss_nn,CENTRAL_B_ONE_SIDE>::position(key00,ginfo,vx_c);
+ key_ret_vx_y = D<y,Field<V,syss_nn>,syss_nn,CENTRAL_B_ONE_SIDE>::position(key00,ginfo,vx_c);
+ key_ret_vy_x = D<x,Field<V,syss_nn>,syss_nn,CENTRAL_B_ONE_SIDE>::position(key00,ginfo,vy_c);
+ key_ret_vy_y = D<y,Field<V,syss_nn>,syss_nn,CENTRAL_B_ONE_SIDE>::position(key00,ginfo,vy_c);
- BOOST_REQUIRE_EQUAL(key_ret_x,0);
- BOOST_REQUIRE_EQUAL(key_ret_y,0);
+ BOOST_REQUIRE_EQUAL(key_ret_vx_x.get(0),-1);
+ BOOST_REQUIRE_EQUAL(key_ret_vx_x.get(1),0);
+ BOOST_REQUIRE_EQUAL(key_ret_vx_y.get(0),-1);
+ BOOST_REQUIRE_EQUAL(key_ret_vx_y.get(1),0);
+ BOOST_REQUIRE_EQUAL(key_ret_vy_y.get(0),0);
+ BOOST_REQUIRE_EQUAL(key_ret_vy_y.get(1),-1);
+ BOOST_REQUIRE_EQUAL(key_ret_vy_x.get(0),0);
+ BOOST_REQUIRE_EQUAL(key_ret_vy_x.get(1),-1);
// Border left*/
-
-
}
BOOST_AUTO_TEST_SUITE_END()
diff --git a/src/FiniteDifference/Laplacian.hpp b/src/FiniteDifference/Laplacian.hpp
index 40af5b0a6a56707aa6a822bde0dbd4485abbff50..d40d91f048bc6545efd79947b63c449ab18e8922 100644
--- a/src/FiniteDifference/Laplacian.hpp
+++ b/src/FiniteDifference/Laplacian.hpp
@@ -18,9 +18,19 @@
template<typename arg, typename Sys_eqs, unsigned int impl=CENTRAL>
class Lap
{
- /*! \brief Create the row of the Matrix
+ /*! \brief Calculate which colums of the Matrix are non zero
*
- * \tparam ord
+ * stub_or_real it is just for change the argument type when testing, in normal
+ * conditions it is a distributed map
+ *
+ * \param pos position where the laplacian 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 Laplacian usage
*
*/
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, std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
@@ -28,10 +38,11 @@ class Lap
std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
}
- /*! \brief Calculate the position where the laplacian is calculated
+ /*! \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
+ * \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)
@@ -40,7 +51,19 @@ class Lap
}
};
-/*! \brief Derivative on direction i
+/*! \brief Laplacian second order approximation CENTRAL Scheme
+ *
+ * \verbatim
+
+ 1.0
+ *
+ | -4.0
+ 1.0 *---+---* 1.0
+ |
+ *
+ 1.0
+
+ * \endverbatim
*
*
*/
@@ -49,8 +72,20 @@ class Lap<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 pos position where the laplacian 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 Laplacian usage
*
*/
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)
@@ -75,13 +110,17 @@ public:
/*! \brief Calculate the position where the derivative is calculated
*
- * In case of 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 CENTRAL Laplacian scheme return the position of the staggered property
+ *
+ * \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 openfpm::vector<comb<Sys_eqs::dims>> & s_pos, long int & fld)
+ 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);
}
};
diff --git a/src/FiniteDifference/eq.hpp b/src/FiniteDifference/eq.hpp
index a01b5d30471566e4818bbcbe0047a3ae408eb935..6d85cc9e6dbff7d37666ac9ffad86b05df36221c 100644
--- a/src/FiniteDifference/eq.hpp
+++ b/src/FiniteDifference/eq.hpp
@@ -14,8 +14,8 @@
#define STAGGERED_GRID 1
#define NORMAL_GRID 0
-#define PERIODIC true
-#define NON_PERIODIC false
+//#define PERIODIC true
+//#define NON_PERIODIC false
#include "data_type/scalar.hpp"
#include "util/util_num.hpp"
@@ -94,10 +94,16 @@ public:
*/
static void value(const map_grid & 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)
{
- if (Sys_eqs::ord == EQS_FIELD)
- cols[g_map.template get<0>(kmap)*Sys_eqs::nvar + f] += coeff;
- else
- cols[g_map.template get<0>(kmap) + f * gs.size()] += coeff;
+ cols[g_map.template get<0>(kmap)*Sys_eqs::nvar + f] += coeff;
+ }
+
+ /*! \brief
+ *
+ *
+ */
+ 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 grid_key_dx<Sys_eqs::dims>(s_pos[f]);
}
};
@@ -108,10 +114,7 @@ class ConstField
inline size_t mat_factor(size_t nvar, size_t sz, const size_t ord)
{
- if (ord == EQS_FIELD)
- return nvar;
- else
- return sz;
+ return nvar;
}
#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_EQ_HPP_ */
diff --git a/src/FiniteDifference/eq_unit_test.hpp b/src/FiniteDifference/eq_unit_test.hpp
index 842ffdc272d26524316333737e4eb0ec7f0558a6..aa4ed9a6f5daec6838b2deb32d3b077963de4a27 100644
--- a/src/FiniteDifference/eq_unit_test.hpp
+++ b/src/FiniteDifference/eq_unit_test.hpp
@@ -21,40 +21,44 @@
BOOST_AUTO_TEST_SUITE( eq_test_suite )
-// Stokes flow
+//! [Definition of the system]
struct lid_nn
{
// dimensionaly of the equation (2D problem 3D problem ...)
static const unsigned int dims = 2;
- // number of fields in the system
+
+ // number of fields in the system v_x, v_y, P so a total of 3
static const unsigned int nvar = 3;
- static const unsigned int ord = EQS_FIELD;
- // boundary at X and Y
+ // boundary conditions PERIODIC OR NON_PERIODIC
static const bool boundary[];
- //
- static constexpr unsigned int num_cfields = 0;
-
// type of space float, double, ...
typedef float stype;
- // type of base grid
+ // 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
typedef grid_dist_id<2,float,aggregate<float[2],float>,CartDecomposition<2,float>> b_grid;
- // type of SparseMatrix
+ // type of SparseMatrix, for the linear system, this parameter is bounded by the solver
+ // that you are using
typedef SparseMatrix<double,int> SparseMatrix_type;
- // type of Vector
+ // type of Vector for the linear system, this parameter is bounded by the solver
+ // that you are using
typedef Vector<double> Vector_type;
- // Define the the underline grid is staggered
+ // Define that the underline grid where we discretize the operators is staggered
static const int grid_type = STAGGERED_GRID;
};
const bool lid_nn::boundary[] = {NON_PERIODIC,NON_PERIODIC};
+//! [Definition of the system]
+
+//! [Definition of the equation of the system in the bulk and at the boundary]
+
// Constant Field
struct eta
{
@@ -63,12 +67,12 @@ struct eta
static float val() {return 1.0;}
};
-// Model the equations
-
+// Convenient constants
constexpr unsigned int v[] = {0,1};
constexpr unsigned int P = 2;
constexpr unsigned int ic = 2;
+// Create field that we have v_x, v_y, P
typedef Field<v[x],lid_nn> v_x;
typedef Field<v[y],lid_nn> v_y;
typedef Field<P,lid_nn> Prs;
@@ -94,6 +98,30 @@ typedef D<y,v_y,lid_nn,FORWARD> dy_vy;
typedef sum<dx_vx,dy_vy,lid_nn> ic_eq;
+// Equation for boundary conditions
+
+/* Consider the staggered cell
+ *
+ \verbatim
+
+ +--$--+
+ | |
+ # * #
+ | |
+ 0--$--+
+
+ # = velocity(y)
+ $ = velocity(x)
+ * = pressure
+
+ \endverbatim
+ *
+ *
+ * If we want to impose v_y = 0 on 0 we have to interpolate between # of this cell
+ * and # of the previous cell on y, (Average) or Avg operator
+ *
+ */
+
// Directional Avg
typedef Avg<x,v_y,lid_nn> avg_vy;
typedef Avg<y,v_x,lid_nn> avg_vx;
@@ -105,57 +133,74 @@ typedef Avg<y,v_x,lid_nn,FORWARD> avg_vx_f;
#define EQ_2 1
#define EQ_3 2
-// Lid driven cavity, uncompressible fluid
+//! [Definition of the equation of the system in the bulk and at the boundary]
+
+// Lid driven cavity, incompressible fluid
BOOST_AUTO_TEST_CASE(lid_driven_cavity)
{
- // Domain
+ //! [lid-driven cavity 2D]
+
+ // Domain, a rectangle
Box<2,float> domain({0.0,0.0},{3.0,1.0});
- // Ghost
+ // Ghost (Not important in this case but required)
Ghost<2,float> g(0.01);
+ // Grid points on x=256 and y=64
long int sz[] = {256,64};
size_t szu[2];
szu[0] = (size_t)sz[0];
szu[1] = (size_t)sz[1];
+ // We need one more point on the left and down part of the domain
+ // This is given by the boundary conditions that we impose, the
+ // reason is mathematical in order to have a well defined system
+ // and cannot be discussed here
Padding<2> pd({1,1},{0,0});
- // Initialize the global VCluster
+ // Initialize openfpm
init_global_v_cluster(&boost::unit_test::framework::master_test_suite().argc,&boost::unit_test::framework::master_test_suite().argv);
- // Initialize openfpm
+ // Distributed grid that store the solution
grid_dist_id<2,float,aggregate<float[2],float>,CartDecomposition<2,float>> g_dist(szu,domain,g);
- // Distributed grid
- FDScheme<lid_nn> fd(pd,domain,g_dist.getGridInfo(),g_dist.getDecomposition(),g_dist.getVC());
-
- // start and end of the bulk
+ // Finite difference scheme
+ FDScheme<lid_nn> fd(pd,domain,g_dist.getGridInfo(),g_dist.getDecomposition());
+ // Here we impose the equation, we start from the incompressibility Eq imposed in the bulk with the
+ // exception of the first point {0,0} and than we set P = 0 in {0,0}, why we are doing this is again
+ // mathematical to have a well defined system, an intuitive explanation is that P and P + c are both
+ // solution for the incompressibility equation, this produce an ill-posed problem to make it well posed
+ // we set one point in this case {0,0} the pressure to a fixed constant for convenience P = 0
fd.impose(ic_eq(),0.0, EQ_3, {0,0},{sz[0]-2,sz[1]-2},true);
fd.impose(Prs(), 0.0, EQ_3, {0,0},{0,0});
+
+ // Here we impose the Eq1 and Eq2
fd.impose(vx_eq(),0.0, EQ_1, {1,0},{sz[0]-2,sz[1]-2});
fd.impose(vy_eq(),0.0, EQ_2, {0,1},{sz[0]-2,sz[1]-2});
- // v_x
- // R L
+ // v_x and v_y
+ // Imposing B1
fd.impose(v_x(),0.0, EQ_1, {0,0},{0,sz[1]-2});
- fd.impose(v_x(),0.0, EQ_1, {sz[0]-1,0},{sz[0]-1,sz[1]-2});
-
- // T B
- fd.impose(avg_vx_f(),0.0, EQ_1, {0,-1},{sz[0]-1,-1});
- fd.impose(avg_vx(),0.0, EQ_1, {0,sz[1]-1},{sz[0]-1,sz[1]-1});
-
- // v_y
- // R L
fd.impose(avg_vy_f(),0.0, EQ_2 , {-1,0},{-1,sz[1]-1});
+ // Imposing B2
+ fd.impose(v_x(),0.0, EQ_1, {sz[0]-1,0},{sz[0]-1,sz[1]-2});
fd.impose(avg_vy(),1.0, EQ_2, {sz[0]-1,0},{sz[0]-1,sz[1]-1});
- // T B
+ // Imposing B3
+ fd.impose(avg_vx_f(),0.0, EQ_1, {0,-1},{sz[0]-1,-1});
fd.impose(v_y(), 0.0, EQ_2, {0,0},{sz[0]-2,0});
+ // Imposing B4
+ fd.impose(avg_vx(),0.0, EQ_1, {0,sz[1]-1},{sz[0]-1,sz[1]-1});
fd.impose(v_y(), 0.0, EQ_2, {0,sz[1]-1},{sz[0]-2,sz[1]-1});
+ // When we pad the grid, there are points of the grid that are not
+ // touched by the previous condition. Mathematically this lead
+ // to have too many variables for the conditions that we are imposing.
+ // Here we are imposing variables that we do not touch to zero
+ //
+
// Padding pressure
fd.impose(Prs(), 0.0, EQ_3, {-1,-1},{sz[0]-1,-1});
fd.impose(Prs(), 0.0, EQ_3, {-1,sz[1]-1},{sz[0]-1,sz[1]-1});
@@ -168,9 +213,11 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
auto x = umfpack_solver<double>::solve(fd.getA(),fd.getB());
- // Bring the solution to grid
+ // 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]
+
g_dist.write("lid_driven_cavity");
}
diff --git a/src/FiniteDifference/mul.hpp b/src/FiniteDifference/mul.hpp
index fe8791a47eb0057c77dbc839351a1a570845ecd9..15ae6056eab90ef92d46ed0270e11808aefdba8e 100644
--- a/src/FiniteDifference/mul.hpp
+++ b/src/FiniteDifference/mul.hpp
@@ -88,7 +88,7 @@ struct const_mul_functor_value
/*! \brief It model an expression expr1 * expr2
*
- * \warning expr1 MUST be a constant expression
+ * \warning expr1 MUST be a constant expression only expr2 depend form variable, this requirement ensure linearity in the solving variable of the equations
*
* \tparam expr1
* \tparam expr2
@@ -105,9 +105,12 @@ struct mul
typedef typename boost::mpl::at<v_expr, boost::mpl::int_<v_sz::type::value - 1> >::type Sys_eqs;
- /*! \brief Create the row of the Matrix
+ /*! \brief Calculate which colums of the Matrix are non zero
*
- * \tparam ord
+ * \param pos position where the multiplication is calculated
+ * \param gs Grid info
+ * \param cols non-zero colums calculated by the function
+ * \param coeff coefficent (constant in front of the derivative)
*
*/
inline static void value(const grid_dist_id<Sys_eqs::dims,typename Sys_eqs::stype,scalar<size_t>,typename Sys_eqs::b_grid::decomposition> & 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)
@@ -123,15 +126,18 @@ struct mul
last_m::value(g_map,kmap,gs,spacing,cols,mfv.coeff);
}
- /*! \brief Calculate the position where the derivative is calculated
+ /*! \brief Calculate the position in the cell where the mul operator is performed
*
- * 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 just return the position of the staggered property in the last expression
+ *
+ * \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 comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
{
- std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
+ return boost::mpl::at<v_expr, boost::mpl::int_<v_sz::type::value - 2> >::type::position(pos,gs,s_pos);
}
};
diff --git a/src/FiniteDifference/sum.hpp b/src/FiniteDifference/sum.hpp
index e1b58af7495be1c2d7e9328e7f70410b5dc6670e..5e485e88f24f99c9b6c9025a896019ad1be1102c 100644
--- a/src/FiniteDifference/sum.hpp
+++ b/src/FiniteDifference/sum.hpp
@@ -51,14 +51,13 @@ struct sum_functor_value
:cols(cols),gs(gs),g_map(g_map),kmap(kmap),key(key),coeff(coeff),spacing(spacing)
{};
-
-
//! It call this function for every expression in the sum
template<typename T>
void operator()(T& t) const
{
boost::mpl::at<v_expr, boost::mpl::int_<T::value> >::type::value(g_map,kmap,gs,spacing,cols,coeff);
}
+
};
/*! \brief It model an expression expr1 + ... exprn
@@ -66,6 +65,10 @@ struct sum_functor_value
* \tparam expr.. two or more expression to be summed
* \tparam Sys_eqs stystem specification
*
+ * ## Example
+ *
+ * \snippet FDScheme_unit_tests.hpp sum example
+ *
*/
template<typename ... expr>
struct sum
@@ -78,9 +81,16 @@ struct sum
typedef typename boost::mpl::at<v_expr, boost::mpl::int_<v_sz::type::value - 1> >::type Sys_eqs;
- /*! \brief Create the row of the Matrix
+ /*! \brief Calculate which colums of the Matrix are non zero
*
- * \tparam ord
+ * \param pos position where the derivative 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 sum example
*
*/
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)
@@ -92,15 +102,19 @@ struct sum
boost::mpl::for_each_ref< boost::mpl::range_c<int,0,v_sz::type::value - 1> >(sm);
}
- /*! \brief Calculate the position where the derivative is calculated
+ /*! \brief Calculate the position in the cell where the sum operator is performed
*
- * 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 is done for the first element, the others are supposed to be in the same position
+ * it just return the position of the staggered property in the first expression
+ *
+ * \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 comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
{
- std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
+ return boost::mpl::at<v_expr, boost::mpl::int_<0> >::type::position(pos,gs,s_pos);
}
};
@@ -112,19 +126,20 @@ struct minus
*
* \tparam ord
*
+ * \snipper FDScheme_unit_tests.hpp minus example
+ *
*/
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)
{
arg::value(g_map,kmap,gs,spacing,cols,-coeff);
}
- /*! \brief Calculate the position where the derivative is calculated
+ /*! \brief Calculate the position where the minus 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 just return the position of the staggered property at first expression
*
*/
- 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 comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
{
std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
}
diff --git a/src/Makefile.am b/src/Makefile.am
index 7001f3d419507c1408a5d20f04fb5887ab5e2883..a417cea0a34341a360b50a9d16703bd7fdd78858 100755
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -3,7 +3,7 @@ LINKLIBS = $(SUITESPARSE_LIBS) $(LAPACK_LIBS) $(BLAS_LIBS) $(METIS_LIB) $(DEFA
noinst_PROGRAMS = numerics
numerics_SOURCES = main.cpp ../../openfpm_vcluster/src/VCluster.cpp ../../openfpm_devices/src/memory/HeapMemory.cpp ../../openfpm_devices/src/memory/PtrMemory.cpp ../../openfpm_devices/src/Memleak_check.cpp
-numerics_CXXFLAGS = $(INCLUDES_PATH) $(BOOST_CPPFLAGS) $(SUITESPARSE_INCLUDE) $(METIS_INCLUDE) $(EIGEN_INCLUDE) -Wno-deprecated-declarations
+numerics_CXXFLAGS = $(INCLUDES_PATH) $(BOOST_CPPFLAGS) $(SUITESPARSE_INCLUDE) $(METIS_INCLUDE) $(EIGEN_INCLUDE) -Wno-deprecated-declarations -Wunused-local-typedefs
numerics_CFLAGS = $(CUDA_CFLAGS)
numerics_LDADD = $(LINKLIBS) -lmetis
diff --git a/src/Vector/Vector_eigen_util.hpp b/src/Vector/Vector_eigen_util.hpp
index 823500090d48ab6142f87b73b383ebb114f17df9..da16b2a7868c1b429e6fdb7c4cf08b5664c79a42 100644
--- a/src/Vector/Vector_eigen_util.hpp
+++ b/src/Vector/Vector_eigen_util.hpp
@@ -22,10 +22,7 @@ struct copy_ele_sca_array
{
template<typename Grid> inline static void copy(Grid & grid_dst, const grid_dist_key_dx<Eqs_sys::dims> & key, const Ev & x,size_t lin_id, size_t base_id, size_t gs_size)
{
- if (Eqs_sys::ord == EQS_FIELD)
- grid_dst.template get<T::value>(key) = x(lin_id * Eqs_sys::nvar + base_id);
- else
- grid_dst.template get<T::value>(key) = x(base_id * gs_size + lin_id);
+ grid_dst.template get<T::value>(key) = x(lin_id * Eqs_sys::nvar + base_id);
}
};
@@ -45,10 +42,7 @@ struct copy_ele_sca_array<copy_type,T,Ev,Eqs_sys,1>
{
for (size_t i = 0 ; i < std::extent<copy_type>::value ; i++)
{
- if (Eqs_sys::ord == EQS_FIELD)
- grid_dst.template get<T::value>(key)[i] = x(lin_id * Eqs_sys::nvar + base_id + i);
- else
- grid_dst.template get<T::value>(key)[i] = x(base_id * gs_size + lin_id);
+ grid_dst.template get<T::value>(key)[i] = x(lin_id * Eqs_sys::nvar + base_id + i);
}
}
};
@@ -70,14 +64,10 @@ struct interp_ele_sca_array
for (size_t i = 0 ; i < interp_pos.get(0).size() ; i++)
{
- size_t gs_size = g_map.getGridInfoVoid().size();
auto key_m = key_src.move(interp_pos.get(0)[i]);
size_t lin_id = g_map.template get<0>(key_m);
- if (scheme::Sys_eqs_typ::ord == EQS_FIELD)
- grid_dst.template get<T::value>(key) += x(lin_id * scheme::Sys_eqs_typ::nvar + base_id);
- else
- grid_dst.template get<T::value>(key) += x(base_id * gs_size + lin_id);
+ grid_dst.template get<T::value>(key) += x(lin_id * scheme::Sys_eqs_typ::nvar + base_id);
division += 1.0;
}
@@ -106,14 +96,10 @@ struct interp_ele_sca_array<copy_type,T,Ev,scheme,1>
division = 0.0;
for (size_t i = 0 ; i < interp_pos.get(j).size() ; i++)
{
- size_t gs_size = g_map.getGridInfoVoid().size();
auto key_m = key_src.move(interp_pos.get(j)[i]);
size_t lin_id = g_map.template get<0>(key_m);
- if (scheme::Sys_eqs_typ::ord == EQS_FIELD)
- grid_dst.template get<T::value>(key)[j] += x(lin_id * scheme::Sys_eqs_typ::nvar + base_id);
- else
- grid_dst.template get<T::value>(key)[j] += x(base_id * gs_size + lin_id);
+ grid_dst.template get<T::value>(key)[j] += x(lin_id * scheme::Sys_eqs_typ::nvar + base_id);
division += 1.0;
}