Commit 5e304636 authored by incardon's avatar incardon

Adding left type expression and right type expression

parent bc9efd0e
...@@ -12,10 +12,37 @@ ...@@ -12,10 +12,37 @@
#include "PointIteratorSkin.hpp" #include "PointIteratorSkin.hpp"
#include "Vector/vector_dist.hpp" #include "Vector/vector_dist.hpp"
/*! \brief A class to draw/create particles based on simple shaped
*
* ## Draw box example
*
* \snippet Draw/DrawParticles_unit_tests.hpp DrawBox_example
*
*/
class DrawParticles class DrawParticles
{ {
public: public:
/*! \brief Draw particles in a box B excluding the area of a second box A (B - A)
*
* The function return an iterator over particles defined on a virtual grid in the simulation domain.
*
* \param vd particles where we are creating the particles
* \param sz indicate the grid size of the virtual grid.
* \param domain Domain where the virtual grid is defined (Must match the domain of vd)
* \param sub_B box contained in domain that define where the particle iterator must iterate,
* particles are placed strictly inside this box
* \param sub_A box contained in the domain that define where the particle iterator should not
* iterate (excluding area)
*
* \note Suppose to have a simulation domain of 1.5 on x and we use sz = 16. Consider now
* to have particles with spacing 0.1 on x. if we define a sub_A that on extend from 0.65
* to 0.95 the first fluid particle
* is at 0.70 and the last is at 0.90
*
* \return an iterator to the selected particles
*
*/
template<unsigned int dim, typename T, typename aggr, typename Decomposition> static PointIteratorSkin<dim,T,Decomposition> template<unsigned int dim, typename T, typename aggr, typename Decomposition> static PointIteratorSkin<dim,T,Decomposition>
DrawSkin(vector_dist<dim,T,aggr,Decomposition> & vd, DrawSkin(vector_dist<dim,T,aggr,Decomposition> & vd,
size_t (& sz)[dim], size_t (& sz)[dim],
...@@ -28,9 +55,30 @@ public: ...@@ -28,9 +55,30 @@ public:
for (size_t i = 0 ; i < dim ; i++) for (size_t i = 0 ; i < dim ; i++)
bc[i] = NON_PERIODIC; bc[i] = NON_PERIODIC;
return PointIteratorSkin<dim,T,Decomposition>(vd.getDecomposition(),sz,domain,sub_A, sub_B, bc); return PointIteratorSkin<dim,T,Decomposition>(vd.getDecomposition(),sz,vd.getDecomposition().getDomain(),sub_A, sub_B, bc);
} }
/*! \brief Draw particles in a box B excluding the areas of an array of boxes A_n
*
* The function return an iterator over particles defined on a virtual grid in the simulation domain.
*
* \param vd particles where we are creating the particles
* \param sz indicate the grid size of the virtual grid.
* \param domain Domain where the virtual grid is defined (Must match the domain of vd)
* \param sub_B box contained in domain that define where the particle iterator must iterate,
* particles are placed strictly inside this box
* \param sub_A array of boxes contained in the domain that define where the particle iterator should not
* iterate (excluding areas)
*
* \note Suppose to have a simulation domain of 1.5 on x and we use sz = 16. Consider now
* to have particles with spacing 0.1 on x. if we define a sub_A that on extend from 0.65
* to 0.95 the first fluid particle
* is at 0.70 and the last is at 0.90
*
* \return an iterator to the selected particles
*
*/
template<unsigned int dim, typename T, typename aggr, typename Decomposition> static PointIteratorSkin<dim,T,Decomposition> template<unsigned int dim, typename T, typename aggr, typename Decomposition> static PointIteratorSkin<dim,T,Decomposition>
DrawSkin(vector_dist<dim,T,aggr,Decomposition> & vd, DrawSkin(vector_dist<dim,T,aggr,Decomposition> & vd,
size_t (& sz)[dim], size_t (& sz)[dim],
...@@ -43,7 +91,7 @@ public: ...@@ -43,7 +91,7 @@ public:
for (size_t i = 0 ; i < dim ; i++) for (size_t i = 0 ; i < dim ; i++)
bc[i] = NON_PERIODIC; bc[i] = NON_PERIODIC;
PointIteratorSkin<dim,T,Decomposition> it(vd.getDecomposition(),sz,domain,sub_A.get(0), sub_B, bc); PointIteratorSkin<dim,T,Decomposition> it(vd.getDecomposition(),sz,vd.getDecomposition().getDomain(),sub_A.get(0), sub_B, bc);
for (size_t i = 1 ; i < sub_A.size() ; i++) for (size_t i = 1 ; i < sub_A.size() ; i++)
it.addBoxA(Box<dim,T>(sub_A.get(i))); it.addBoxA(Box<dim,T>(sub_A.get(i)));
...@@ -51,13 +99,31 @@ public: ...@@ -51,13 +99,31 @@ public:
return it; return it;
} }
/*! \brief Draw particles in a box
*
* The function return an iterator over particles defined on a virtual grid in the simulation domain.
*
* \param vd particles where we are creating the particles
* \param sz indicate the grid size of the virtual grid.
* \param domain Domain where the virtual grid is defined (Must match the domain of vd)
* \param sub box contained in domain that define where the particle iterator must iterate,
* particles are placed strictly inside this box
*
* \note Suppose to have a simulation domain of 1.5 on x and we use sz = 16. Consider now
* to have particles with spacing 0.1 on x. if we define a sub box that on extend from 0.65
* to 0.95 the first fluid particle
* is at 0.70 and the last is at 0.90
*
* \return an iterator to the selected particles
*
*/
template<unsigned int dim, typename T, typename aggr, typename Decomposition> static PointIterator<dim,T,Decomposition> template<unsigned int dim, typename T, typename aggr, typename Decomposition> static PointIterator<dim,T,Decomposition>
DrawBox(vector_dist<dim,T,aggr,Decomposition> & vd, DrawBox(vector_dist<dim,T,aggr,Decomposition> & vd,
size_t (& sz)[dim], size_t (& sz)[dim],
Box<dim,T> & domain, Box<dim,T> & domain,
Box<dim,T> & sub) Box<dim,T> & sub)
{ {
return PointIterator<dim,T,Decomposition>(vd.getDecomposition(),sz,domain,sub); return PointIterator<dim,T,Decomposition>(vd.getDecomposition(),sz,vd.getDecomposition().getDomain(),sub);
} }
}; };
......
...@@ -16,6 +16,8 @@ BOOST_AUTO_TEST_SUITE( draw_particles ) ...@@ -16,6 +16,8 @@ BOOST_AUTO_TEST_SUITE( draw_particles )
BOOST_AUTO_TEST_CASE(point_iterator) BOOST_AUTO_TEST_CASE(point_iterator)
{ {
//! [DrawBox_example]
size_t sz[] = {23,27,20}; size_t sz[] = {23,27,20};
Box<3,double> domain({-1.2,0.5,-0.6},{1.0,3.1,1.3}); Box<3,double> domain({-1.2,0.5,-0.6},{1.0,3.1,1.3});
...@@ -52,6 +54,8 @@ BOOST_AUTO_TEST_CASE(point_iterator) ...@@ -52,6 +54,8 @@ BOOST_AUTO_TEST_CASE(point_iterator)
++p; ++p;
} }
//! [DrawBox_example]
Vcluster & v_cl = create_vcluster(); Vcluster & v_cl = create_vcluster();
v_cl.sum(cnt); v_cl.sum(cnt);
......
...@@ -69,7 +69,7 @@ class PointIterator: protected grid_dist_id_iterator_dec<Decomposition> ...@@ -69,7 +69,7 @@ class PointIterator: protected grid_dist_id_iterator_dec<Decomposition>
//! Spacing //! Spacing
T sp[dim]; T sp[dim];
static grid_key_dx<dim> getStart(size_t (& gs)[dim], Box<dim,T> & dom, Box<dim,T> & sub_dom, T (& sp)[dim]) static grid_key_dx<dim> getStart(size_t (& gs)[dim], const Box<dim,T> & dom, Box<dim,T> & sub_dom, T (& sp)[dim])
{ {
for (size_t i = 0 ; i < dim ; i++) for (size_t i = 0 ; i < dim ; i++)
sp[i] = (dom.getHigh(i) - dom.getLow(i)) / (gs[i] -1); sp[i] = (dom.getHigh(i) - dom.getLow(i)) / (gs[i] -1);
...@@ -82,7 +82,7 @@ class PointIterator: protected grid_dist_id_iterator_dec<Decomposition> ...@@ -82,7 +82,7 @@ class PointIterator: protected grid_dist_id_iterator_dec<Decomposition>
return pkey; return pkey;
} }
static grid_key_dx<dim> getStop(size_t (& gs)[dim], Box<dim,T> & dom, Box<dim,T> & sub_dom, T (& sp)[dim]) static grid_key_dx<dim> getStop(size_t (& gs)[dim], const Box<dim,T> & dom, Box<dim,T> & sub_dom, T (& sp)[dim])
{ {
for (size_t i = 0 ; i < dim ; i++) for (size_t i = 0 ; i < dim ; i++)
sp[i] = (dom.getHigh(i) - dom.getLow(i)) / (gs[i] - 1); sp[i] = (dom.getHigh(i) - dom.getLow(i)) / (gs[i] - 1);
...@@ -97,6 +97,9 @@ class PointIterator: protected grid_dist_id_iterator_dec<Decomposition> ...@@ -97,6 +97,9 @@ class PointIterator: protected grid_dist_id_iterator_dec<Decomposition>
void calculateAp() void calculateAp()
{ {
if (grid_dist_id_iterator_dec<Decomposition>::isNext() == false)
return;
grid_key_dx<dim> key = grid_dist_id_iterator_dec<Decomposition>::get(); grid_key_dx<dim> key = grid_dist_id_iterator_dec<Decomposition>::get();
for (size_t i = 0 ; i < dim ; i++) for (size_t i = 0 ; i < dim ; i++)
...@@ -110,7 +113,7 @@ public: ...@@ -110,7 +113,7 @@ public:
* \param sp grid spacing * \param sp grid spacing
* *
*/ */
PointIterator( Decomposition & dec, size_t (& sz)[dim], Box<dim,T> & domain, Box<dim,T> & sub_domain) PointIterator( Decomposition & dec, size_t (& sz)[dim], const Box<dim,T> & domain, Box<dim,T> & sub_domain)
:grid_dist_id_iterator_dec<Decomposition>(dec, sz, getStart(sz,domain,sub_domain,sp), getStop(sz,domain,sub_domain,sp)),sub_domain(sub_domain),domain(domain) :grid_dist_id_iterator_dec<Decomposition>(dec, sz, getStart(sz,domain,sub_domain,sp), getStop(sz,domain,sub_domain,sp)),sub_domain(sub_domain),domain(domain)
{ {
calculateAp(); calculateAp();
...@@ -133,6 +136,7 @@ public: ...@@ -133,6 +136,7 @@ public:
PointIterator & operator++() PointIterator & operator++()
{ {
grid_dist_id_iterator_dec<Decomposition>::operator++(); grid_dist_id_iterator_dec<Decomposition>::operator++();
calculateAp(); calculateAp();
return *this; return *this;
......
...@@ -101,6 +101,9 @@ class PointIteratorSkin: protected grid_dist_id_iterator_dec_skin<Decomposition> ...@@ -101,6 +101,9 @@ class PointIteratorSkin: protected grid_dist_id_iterator_dec_skin<Decomposition>
void calculateAp() void calculateAp()
{ {
if (grid_dist_id_iterator_dec_skin<Decomposition>::isNext() == false)
return;
grid_key_dx<dim> key = grid_dist_id_iterator_dec_skin<Decomposition>::get(); grid_key_dx<dim> key = grid_dist_id_iterator_dec_skin<Decomposition>::get();
for (size_t i = 0 ; i < dim ; i++) for (size_t i = 0 ; i < dim ; i++)
...@@ -133,7 +136,7 @@ public: ...@@ -133,7 +136,7 @@ public:
* \param sp grid spacing * \param sp grid spacing
* *
*/ */
PointIteratorSkin( Decomposition & dec, size_t (& sz)[dim], Box<dim,T> & domain, const Box<dim,T> & sub_A, const Box<dim,T> & sub_B, size_t (& bc)[dim]) PointIteratorSkin( Decomposition & dec, size_t (& sz)[dim], const Box<dim,T> & domain, const Box<dim,T> & sub_A, const Box<dim,T> & sub_B, size_t (& bc)[dim])
:grid_dist_id_iterator_dec_skin<Decomposition>(dec, sz, getAB(sz,domain,sub_A,sub_B,sp,RETURN_A), getAB(sz,domain,sub_A,sub_B,sp,RETURN_B), bc),domain(domain) :grid_dist_id_iterator_dec_skin<Decomposition>(dec, sz, getAB(sz,domain,sub_A,sub_B,sp,RETURN_A), getAB(sz,domain,sub_A,sub_B,sp,RETURN_B), bc),domain(domain)
{ {
sub_domainA.add(sub_A); sub_domainA.add(sub_A);
...@@ -142,6 +145,7 @@ public: ...@@ -142,6 +145,7 @@ public:
/*! \Return the actual point /*! \Return the actual point
* *
* \return the actual point
* *
*/ */
Point<dim,T> & get() Point<dim,T> & get()
...@@ -157,11 +161,13 @@ public: ...@@ -157,11 +161,13 @@ public:
PointIteratorSkin & operator++() PointIteratorSkin & operator++()
{ {
grid_dist_id_iterator_dec_skin<Decomposition>::operator++(); grid_dist_id_iterator_dec_skin<Decomposition>::operator++();
calculateAp(); calculateAp();
while (grid_dist_id_iterator_dec_skin<Decomposition>::isNext() && isValidPoint() == false) while (grid_dist_id_iterator_dec_skin<Decomposition>::isNext() && isValidPoint() == false)
{ {
grid_dist_id_iterator_dec_skin<Decomposition>::operator++(); grid_dist_id_iterator_dec_skin<Decomposition>::operator++();
calculateAp(); calculateAp();
} }
......
...@@ -52,14 +52,20 @@ class D ...@@ -52,14 +52,20 @@ class D
* In case of non staggered case this function just return a null grid_key, in case of staggered, * 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 * it calculate how the operator shift the calculation in the cell
* *
* \param position where we are calculating the derivative * \param pos position where we are calculating the derivative
* \param gs Grid info * \param gs Grid info
* \param s_pos staggered position of the properties * \param s_pos staggered position of the properties
* *
* \return where (in which cell) the derivative is calculated
*
*/ */
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]) 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"; std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
return pos;
} }
}; };
...@@ -80,8 +86,10 @@ class D<d,arg,Sys_eqs,CENTRAL> ...@@ -80,8 +86,10 @@ class D<d,arg,Sys_eqs,CENTRAL>
/*! \brief Calculate which colums of the Matrix are non zero /*! \brief Calculate which colums of the Matrix are non zero
* *
* \param pos position where the derivative is calculated * \param g_map mapping grid
* \param kmap position where the derivative is calculated
* \param gs Grid info * \param gs Grid info
* \param spacing grid spacing
* \param cols non-zero colums calculated by the function * \param cols non-zero colums calculated by the function
* \param coeff coefficent (constant in front of the derivative) * \param coeff coefficent (constant in front of the derivative)
* *
...@@ -90,7 +98,12 @@ class D<d,arg,Sys_eqs,CENTRAL> ...@@ -90,7 +98,12 @@ class D<d,arg,Sys_eqs,CENTRAL>
* \snippet FDScheme_unit_tests.hpp Usage of stencil derivative * \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::extended_type>::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) inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::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)
{ {
// if the system is staggered the CENTRAL derivative is equivalent to a forward derivative // if the system is staggered the CENTRAL derivative is equivalent to a forward derivative
if (is_grid_staggered<Sys_eqs>::value()) if (is_grid_staggered<Sys_eqs>::value())
...@@ -141,12 +154,16 @@ class D<d,arg,Sys_eqs,CENTRAL> ...@@ -141,12 +154,16 @@ class D<d,arg,Sys_eqs,CENTRAL>
* *
* \f$ \frac{\partial v_x}{\partial y} \f$ is calculated on position (*), so the function return the grid_key {0,0} * \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 pos position where we are calculating the derivative
* \param gs Grid info * \param gs Grid info
* \param s_pos staggered position of the properties * \param s_pos staggered position of the properties
* *
* \return where (in which cell grid) the derivative is calculated
*
*/ */
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]) 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); auto arg_pos = arg::position(pos,gs,s_pos);
if (is_grid_staggered<Sys_eqs>::value()) if (is_grid_staggered<Sys_eqs>::value())
...@@ -195,8 +212,10 @@ public: ...@@ -195,8 +212,10 @@ public:
/*! \brief Calculate which colums of the Matrix are non zero /*! \brief Calculate which colums of the Matrix are non zero
* *
* \param pos position where the derivative is calculated * \param g_map mapping grid points
* \param kmap position where the derivative is calculated
* \param gs Grid info * \param gs Grid info
* \param spacing of the grid
* \param cols non-zero colums calculated by the function * \param cols non-zero colums calculated by the function
* \param coeff coefficent (constant in front of the derivative) * \param coeff coefficent (constant in front of the derivative)
* *
...@@ -205,7 +224,12 @@ public: ...@@ -205,7 +224,12 @@ public:
* \snippet FDScheme_unit_tests.hpp Usage of stencil derivative * \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::extended_type>::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) static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::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 #ifdef SE_CLASS1
if (Sys_eqs::boundary[d] == PERIODIC) if (Sys_eqs::boundary[d] == PERIODIC)
...@@ -279,10 +303,12 @@ public: ...@@ -279,10 +303,12 @@ public:
* outside the boundary is simply the central scheme, at the boundary it is simply the * outside the boundary is simply the central scheme, at the boundary it is simply the
* staggered position of the property * staggered position of the property
* *
* \param position where we are calculating the derivative * \param pos position where we are calculating the derivative
* \param gs Grid info * \param gs Grid info
* \param s_pos staggered position of the properties * \param s_pos staggered position of the properties
* *
* \return where (in which cell grid) the derivative is calculated
*
*/ */
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]) 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])
{ {
...@@ -308,8 +334,10 @@ class D<d,arg,Sys_eqs,FORWARD> ...@@ -308,8 +334,10 @@ class D<d,arg,Sys_eqs,FORWARD>
/*! \brief Calculate which colums of the Matrix are non zero /*! \brief Calculate which colums of the Matrix are non zero
* *
* \param pos position where the derivative is calculated * \param g_map mapping grid
* \param kmap position where the derivative is calculated
* \param gs Grid info * \param gs Grid info
* \param spacing grid spacing
* \param cols non-zero colums calculated by the function * \param cols non-zero colums calculated by the function
* \param coeff coefficent (constant in front of the derivative) * \param coeff coefficent (constant in front of the derivative)
* *
...@@ -318,7 +346,12 @@ class D<d,arg,Sys_eqs,FORWARD> ...@@ -318,7 +346,12 @@ class D<d,arg,Sys_eqs,FORWARD>
* \snippet FDScheme_unit_tests.hpp Usage of stencil derivative * \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::extended_type>::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) inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::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)
{ {
long int old_val = kmap.getKeyRef().get(d); long int old_val = kmap.getKeyRef().get(d);
...@@ -336,12 +369,16 @@ class D<d,arg,Sys_eqs,FORWARD> ...@@ -336,12 +369,16 @@ class D<d,arg,Sys_eqs,FORWARD>
* In case of non staggered case this function just return a null grid_key, in case of staggered, * 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 * the FORWARD scheme return the position of the staggered property
* *
* \param position where we are calculating the derivative * \param pos position where we are calculating the derivative
* \param gs Grid info * \param gs Grid info
* \param s_pos staggered position of the properties * \param s_pos staggered position of the properties
* *
* \return where (in which cell grid) the derivative is calculated
*
*/ */
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]) 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); return arg::position(pos,gs,s_pos);
} }
...@@ -364,8 +401,10 @@ class D<d,arg,Sys_eqs,BACKWARD> ...@@ -364,8 +401,10 @@ class D<d,arg,Sys_eqs,BACKWARD>
/*! \brief Calculate which colums of the Matrix are non zero /*! \brief Calculate which colums of the Matrix are non zero
* *
* \param pos position where the derivative is calculated * \param g_map mapping grid
* \param kmap position where the derivative is calculated
* \param gs Grid info * \param gs Grid info
* \param spacing of the grid
* \param cols non-zero colums calculated by the function * \param cols non-zero colums calculated by the function
* \param coeff coefficent (constant in front of the derivative) * \param coeff coefficent (constant in front of the derivative)
* *
...@@ -374,7 +413,12 @@ class D<d,arg,Sys_eqs,BACKWARD> ...@@ -374,7 +413,12 @@ class D<d,arg,Sys_eqs,BACKWARD>
* \snippet FDScheme_unit_tests.hpp Usage of stencil derivative * \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::extended_type>::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) inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::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)
{ {
long int old_val = kmap.getKeyRef().get(d); long int old_val = kmap.getKeyRef().get(d);
...@@ -392,10 +436,12 @@ class D<d,arg,Sys_eqs,BACKWARD> ...@@ -392,10 +436,12 @@ class D<d,arg,Sys_eqs,BACKWARD>
* In case of non staggered case this function just return a null grid_key, in case of staggered, * 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 * the BACKWARD scheme return the position of the staggered property
* *
* \param position where we are calculating the derivative * \param pos position where we are calculating the derivative
* \param gs Grid info * \param gs Grid info
* \param s_pos staggered position of the properties * \param s_pos staggered position of the properties
* *
* \return where (in which cell grid) the derivative is calculated
*
*/ */
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]) 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])
{ {
......
...@@ -126,59 +126,67 @@ class FDScheme ...@@ -126,59 +126,67 @@ class FDScheme
{ {
public: public:
// Distributed grid map //! Distributed grid map
typedef grid_dist_id<Sys_eqs::dims,typename Sys_eqs::stype,scalar<size_t>,typename Sys_eqs::b_grid::decomposition::extended_type> g_map_type; typedef grid_dist_id<Sys_eqs::dims,typename Sys_eqs::stype,scalar<size_t>,typename Sys_eqs::b_grid::decomposition::extended_type> g_map_type;
//! Type that specify the properties of the system of equations
typedef Sys_eqs Sys_eqs_typ; typedef Sys_eqs Sys_eqs_typ;
private: private:
// Padding //! Padding
Padding<Sys_eqs::dims> pd; Padding<Sys_eqs::dims> pd;
//! Sparse matrix triplet type
typedef typename Sys_eqs::SparseMatrix_type::triplet_type triplet; typedef typename Sys_eqs::SparseMatrix_type::triplet_type triplet;
// Vector b //! Vector b
typename Sys_eqs::Vector_type b; typename Sys_eqs::Vector_type b;
// Domain Grid informations //! Domain Grid informations
const grid_sm<Sys_eqs::dims,void> & gs; const grid_sm<Sys_eqs::dims,void> & gs;
// Get the grid spacing //! Get the grid spacing
typename Sys_eqs::stype spacing[Sys_eqs::dims]; typename Sys_eqs::stype spacing[Sys_eqs::dims];
// mapping grid //! mapping grid
g_map_type g_map; g_map_type g_map;
// row of the matrix //! row of the matrix
size_t row; size_t row;
// row on b //! row on b
size_t row_b; size_t row_b;
// Grid points that has each processor //! Grid points that has each processor
openfpm::vector<size_t> pnt; openfpm::vector<size_t> pnt;
// Staggered position for each property //! Staggered position for each property
comb<Sys_eqs::dims> s_pos[Sys_eqs::nvar]; 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 //! 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 //! 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 //! no processors can have holes in the sequence, this number indicate where the sequence start for this
// processor //! processor
size_t s_pnt; size_t s_pnt;
/*! \brief Equation id + space position
*
*/
struct key_and_eq struct key_and_eq
{ {
//! space position
grid_key_dx<Sys_eqs::dims> key; grid_key_dx<Sys_eqs::dims> key;
//! equation id
size_t eq; size_t eq;
}; };
/*! \brief From the row Matrix position to the spatial position /*! \brief From the row Matrix position to the spatial position
* *
* \param row Matrix * \param row Matrix row
* *
* \return spatial position * \return spatial position + equation id
* *
*/ */
inline key_and_eq from_row_to_key(size_t row) inline key_and_eq from_row_to_key(size_t row)
...@@ -205,9 +213,10 @@ private: ...@@ -205,9 +213,10 @@ private:
return ke; return ke;
} }
/* \brief calculate the mapping grid size with padding /*! \brief calculate the mapping grid size with padding
* *
* \param gs original grid size * \param sz original grid size
* \param pd padding
* *
* \return padded grid size * \return padded grid size
* *
...@@ -323,7 +332,7 @@ public: ...@@ -323,7 +332,7 @@ public:
/*! \brief set the staggered position for each property /*! \brief set the staggered position for each property
* *
* \param vector containing the staggered position for each property * \param sp vector containing the staggered position for each property
* *
*/ */
void setStagPos(comb<Sys_eqs::dims> (& sp)[Sys_eqs::nvar]) void setStagPos(comb<Sys_eqs::dims> (& sp)[Sys_eqs::nvar])
...@@ -335,10 +344,8 @@ public: ...@@ -335,10 +344,8 @@ public: