Commit 2a805e31 authored by incardon's avatar incardon

Adding Draw particles functionalities

parent dd0e4d4d
/*
* DrawParticlesGrid.hpp
*
* Created on: Jan 3, 2017
* Author: i-bird
*/
#ifndef OPENFPM_NUMERICS_SRC_DRAW_DRAWPARTICLESGRID_HPP_
#define OPENFPM_NUMERICS_SRC_DRAW_DRAWPARTICLESGRID_HPP_
/*! \brief this class draw particles on subset of grid-like position
*
\verbatim
Point B
| (0.6,0.6)
+ + + + + | + + + +
+---------------+
+ | + + + + | + + + +
| |
+ | + + + + | + + + +
| |
+ | + + + + | + + + +
| |
+ | + + + + | + + + +
| |
+ | + + + + | + + + +
| |
+ | + + + + | + + + +
| |
+ | + + + + | + + + +
+---------------+
+ + + + + ^ + + + +
(-1.0,-1.0) |
|
(Point A)
\endverbatim
*
*
* Suppose we have a grid 9x9 from (-1.0,-1.0 to 0.6,0.6)
*
* Defined a Box (-0.9,-0.9 to -0.1,0.5)
*
* This class will return the points
*
* (-0.8 -0.8)
* (-0.6 -0.8)
* (-0.4 -0.8)
* ...
* (-0.1,-0.8) Point A
* ...
* (-0.1,0.5) Point B
*
*/
template<unsigned int dim, typename T>
class DrawParticlesGrid
{
/*! \brief Draw Particles
*
* \param sp grid spacing
*
*/
DrawParticlesGrid(T (& sp)[dim], Box<dim,T> & domain)
{
}
PointIterator getPointIterator()
{
}
};
#endif /* OPENFPM_NUMERICS_SRC_DRAW_DRAWPARTICLESGRID_HPP_ */
/*
* DrawParticles_unit_tests.hpp
*
* Created on: Jan 3, 2017
* Author: i-bird
*/
#ifndef OPENFPM_NUMERICS_SRC_DRAW_DRAWPARTICLES_UNIT_TESTS_HPP_
#define OPENFPM_NUMERICS_SRC_DRAW_DRAWPARTICLES_UNIT_TESTS_HPP_
#include "PointIterator.hpp"
BOOST_AUTO_TEST_SUITE( draw_particles )
BOOST_AUTO_TEST_CASE(point_iterator)
{
size_t sz[] = {30,17,28};
Box<3,double> domain({-1.2,0.5,-0.6},{1.0,3.1,1.312});
Box<3,double> sub_domain({-0.2,0.8,0.2},{1.0,1.1,1.0});
// Boundary conditions
size_t bc[3]={NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
// ghost, big enough to contain the interaction radius
Ghost<3,float> ghost(0.01);
vector_dist<3,double,aggregate<double>> vd(0,domain,bc,ghost);
PointIterator<3,double,CartDecomposition<3,double>> p(vd.getDecomposition(),sz,domain,sub_domain);
size_t cnt = 0;
bool good = true;
while (p.isNext())
{
vd.add();
vd.getLastPos()[0] = p.get().get(0);
vd.getLastPos()[1] = p.get().get(1);
vd.getLastPos()[2] = p.get().get(2);
good &= sub_domain.isInside(p.get());
cnt++;
++p;
}
Vcluster & v_cl = create_vcluster();
v_cl.sum(cnt);
BOOST_REQUIRE_EQUAL(cnt,sz[0]*sz[1]*sz[2]);
BOOST_REQUIRE_EQUAL(good,true);
}
BOOST_AUTO_TEST_SUITE_END()
#endif /* OPENFPM_NUMERICS_SRC_DRAW_DRAWPARTICLES_UNIT_TESTS_HPP_ */
/*
* PointIterator.hpp
*
* Created on: Jan 3, 2017
* Author: i-bird
*/
#ifndef OPENFPM_NUMERICS_SRC_DRAW_POINTITERATOR_HPP_
#define OPENFPM_NUMERICS_SRC_DRAW_POINTITERATOR_HPP_
/*! \brief this class draw particles on subset of grid-like position
*
\verbatim
Point B
| (0.6,0.6)
+ + + + + | + + + +
+---------------+
+ | + + + + | + + + +
| |
+ | + + + + | + + + +
| |
+ | + + + + | + + + +
| |
+ | + + + + | + + + +
| |
+ | + + + + | + + + +
| |
+ | + + + + | + + + +
| |
+ | + + + + | + + + +
+---------------+
+ + + + + ^ + + + +
(-1.0,-1.0) |
|
(Point A)
\endverbatim
*
*
* Suppose we have a grid 9x9 from (-1.0,-1.0 to 0.6,0.6)
*
* Defined a Box (-0.9,-0.9 to -0.1,0.5)
*
* This class will return the points
*
* (-0.8 -0.8)
* (-0.6 -0.8)
* (-0.4 -0.8)
* ...
* (-0.1,-0.8) Point A
* ...
* (-0.1,0.5) Point B
*
*/
template<unsigned int dim, typename T, typename Decomposition>
class PointIterator: protected grid_dist_id_iterator_dec<Decomposition>
{
//! Actual point
Point<dim,T> ap;
//! sub_domain
Box<dim,T> sub_domain;
//! Spacing
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])
{
for (size_t i = 0 ; i < dim ; i++)
sp[i] = (dom.getHigh(i) - dom.getLow(i)) / gs[i];
Point<dim,T> pstart = sub_dom.getP1();
grid_key_dx<dim> pkey;
for (size_t i = 0 ; i < dim ; i++)
pkey.set_d(i,std::ceil( (sub_dom.getLow(i) - dom.getLow(i)) / sp[i]));
return pkey;
}
static grid_key_dx<dim> getStop(size_t (& gs)[dim], Box<dim,T> & dom, Box<dim,T> & sub_dom, T (& sp)[dim])
{
for (size_t i = 0 ; i < dim ; i++)
sp[i] = (dom.getHigh(i) - dom.getLow(i)) / gs[i];
Point<dim,T> pstop = sub_dom.getP2();
grid_key_dx<dim> pkey;
for (size_t i = 0 ; i < dim ; i++)
pkey.set_d(i,std::ceil( (sub_dom.getLow(i) - dom.getLow(i)) / sp[i]));
return pkey;
}
public:
/*! \brief Draw Particles
*
* \param sp grid spacing
*
*/
PointIterator( Decomposition & dec, size_t (& sz)[dim], 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)
{
// Calculate the spacing
for (size_t i = 0 ; i < dim ; i++)
sp[i] = domain.getLow(i) / sz[i];
}
/*! \Return the actual point
*
*
*/
Point<dim,T> & get()
{
return ap;
}
/*! \brief Next point
*
* \return itself
*
*/
PointIterator & operator++()
{
grid_dist_id_iterator_dec<Decomposition>::operator++();
grid_key_dx<dim> key = grid_dist_id_iterator_dec<Decomposition>::get();
for (size_t i = 0 ; i < dim ; i++)
ap.get(i) = key.get(i) * sp[i] + sub_domain.getLow(i);
return *this;
}
bool isNext()
{
return grid_dist_id_iterator_dec<Decomposition>::isNext();
}
};
#endif /* OPENFPM_NUMERICS_SRC_DRAW_POINTITERATOR_HPP_ */
...@@ -15,10 +15,11 @@ ...@@ -15,10 +15,11 @@
#define HAS_VAL 1 #define HAS_VAL 1
#define HAS_POS_VAL 2 #define HAS_POS_VAL 2
//! Evaluate the constant field function
template <unsigned int, typename T> template <unsigned int, typename T>
struct has_val struct has_val
{ {
//! evaluate
static float call_val() static float call_val()
{ {
std::cerr << "Error the type " << demangle(typeid(T).name()) << "interpreted as constant field, does not have a function val() or val_pos(), please see the numeric examples in Finite Differences for more information\n"; std::cerr << "Error the type " << demangle(typeid(T).name()) << "interpreted as constant field, does not have a function val() or val_pos(), please see the numeric examples in Finite Differences for more information\n";
...@@ -26,16 +27,18 @@ struct has_val ...@@ -26,16 +27,18 @@ struct has_val
} }
}; };
//! Evaluate the constant field function
template <typename T> template <typename T>
struct has_val<HAS_VAL,T> struct has_val<HAS_VAL,T>
{ {
//! evaluate
static decltype(T::val()) call_val() static decltype(T::val()) call_val()
{ {
return T::val(); return T::val();
} }
}; };
//! Multiplication expression
template<typename v_expr> template<typename v_expr>
struct const_mul_functor_value struct const_mul_functor_value
{ {
...@@ -48,12 +51,13 @@ struct const_mul_functor_value ...@@ -48,12 +51,13 @@ struct const_mul_functor_value
//! sum functor //! sum functor
std::unordered_map<long int,typename last::stype> & cols; std::unordered_map<long int,typename last::stype> & cols;
//! grid size
const grid_sm<last::dims,void> & gs; const grid_sm<last::dims,void> & gs;
// grid mapping //! grid mapping
const grid_dist_id<last::dims,typename last::stype,scalar<size_t>,typename last::b_grid::decomposition::extended_type> & g_map; const grid_dist_id<last::dims,typename last::stype,scalar<size_t>,typename last::b_grid::decomposition::extended_type> & g_map;
// grid position //! grid position
grid_dist_key_dx<last::dims> & kmap; grid_dist_key_dx<last::dims> & kmap;
//! coefficent //! coefficent
...@@ -63,9 +67,21 @@ struct const_mul_functor_value ...@@ -63,9 +67,21 @@ struct const_mul_functor_value
typename last::stype (& spacing)[last::dims]; typename last::stype (& spacing)[last::dims];
/*! \brief constructor /*! \brief constructor
*
* \param g_map mapping grid
* \param kmap grid point (row) where we evaluate the non-zero colums
* \param gs grid size
* \param spacing grid spacing
* \param col non zero colums
* \param coefficent
* *
*/ */
const_mul_functor_value(const grid_dist_id<last::dims,typename last::stype,scalar<size_t>,typename last::b_grid::decomposition::extended_type> & g_map, grid_dist_key_dx<last::dims> & kmap, const grid_sm<last::dims,void> & gs, typename last::stype (& spacing)[last::dims], std::unordered_map<long int,typename last::stype> & cols, typename last::stype coeff) const_mul_functor_value(const grid_dist_id<last::dims,typename last::stype,scalar<size_t>,typename last::b_grid::decomposition::extended_type> & g_map,
grid_dist_key_dx<last::dims> & kmap,
const grid_sm<last::dims,void> & gs,
typename last::stype (& spacing)[last::dims],
std::unordered_map<long int,typename last::stype> & cols,
typename last::stype coeff)
:cols(cols),gs(gs),g_map(g_map),kmap(kmap),coeff(coeff),spacing(spacing) :cols(cols),gs(gs),g_map(g_map),kmap(kmap),coeff(coeff),spacing(spacing)
{}; {};
...@@ -80,6 +96,11 @@ struct const_mul_functor_value ...@@ -80,6 +96,11 @@ struct const_mul_functor_value
coeff *= has_val<is_const_field<cfield>::value * 1,cfield>::call_val(); coeff *= has_val<is_const_field<cfield>::value * 1,cfield>::call_val();
} }
/*! \brief Get the coefficent
*
* \return the coefficent
*
*/
typename last::stype getCoeff() typename last::stype getCoeff()
{ {
return coeff; return coeff;
......
...@@ -8,27 +8,34 @@ ...@@ -8,27 +8,34 @@
#ifndef OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATOR_ASSIGN_HPP_ #ifndef OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATOR_ASSIGN_HPP_
#define OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATOR_ASSIGN_HPP_ #define OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATOR_ASSIGN_HPP_
//! Construct a vector expression from a type T that is already an expression
//! it does nothing
template<typename T> template<typename T>
struct construct_expression struct construct_expression
{ {
//! It return the expression itself
static inline const T & construct(const T & e) static inline const T & construct(const T & e)
{ {
return e; return e;
} }
}; };
//! construct a vector expression from a double
template<> template<>
struct construct_expression<double> struct construct_expression<double>
{ {
//! It return an expression from a double
static inline vector_dist_expression<0,double> construct(double e) static inline vector_dist_expression<0,double> construct(double e)
{ {
return vector_dist_expression<0,double>(e); return vector_dist_expression<0,double>(e);
} }
}; };
//! construct a vector expression from a float
template<> template<>
struct construct_expression<float> struct construct_expression<float>
{ {
//! It return an expression from a float
static inline vector_dist_expression<0,float> construct(const float e) static inline vector_dist_expression<0,float> construct(const float e)
{ {
return vector_dist_expression<0,float>(e); return vector_dist_expression<0,float>(e);
......
...@@ -53,6 +53,7 @@ struct set_zero ...@@ -53,6 +53,7 @@ struct set_zero
} }
}; };
//! Create a point with all compunent set to zero
template<unsigned int dim, typename T> template<unsigned int dim, typename T>
struct set_zero<Point<dim,T>> struct set_zero<Point<dim,T>>
{ {
...@@ -276,7 +277,7 @@ public: ...@@ -276,7 +277,7 @@ public:
* \param o1 expression * \param o1 expression
* \param cl Cell-list * \param cl Cell-list
* \param ker kernel function * \param ker kernel function
* \param vector_orig vector containing the particle positions * \param vd vector containing the particle positions
* *
*/ */
vector_dist_expression_op(const exp1 & o1, NN & cl, Kernel & ker, const vector_orig & vd) vector_dist_expression_op(const exp1 & o1, NN & cl, Kernel & ker, const vector_orig & vd)
......
...@@ -51,9 +51,6 @@ template <typename vector,typename Kernel, typename NN_type> bool check_values_a ...@@ -51,9 +51,6 @@ template <typename vector,typename Kernel, typename NN_type> bool check_values_a
base2 += vd.template getProp<C>(p); base2 += vd.template getProp<C>(p);
if (base1 != base2)
std::cout << "CAZZO " << base1 << " " << base2 << std::endl;
ret &= base1 == base2; ret &= base1 == base2;
++it2; ++it2;
......
...@@ -23,19 +23,20 @@ template <unsigned int dim, typename T> inline vector_dist_expression<16384,Poin ...@@ -23,19 +23,20 @@ template <unsigned int dim, typename T> inline vector_dist_expression<16384,Poin
} }
/*! \brief Main class that encapsulate a vector properties /*! \brief This class represent a constant parameter in a vector expression
* *
* \tparam prp property involved * \tparam point type of point
* \tparam vector involved
* *
*/ */
template<typename point> template<typename point>
class vector_dist_expression<16384,point> class vector_dist_expression<16384,point>
{ {
//! constant point stored
point p; point p;
public: public:
//! vector expression from a constant point
vector_dist_expression(point p) vector_dist_expression(point p)
:p(p) :p(p)
{} {}
...@@ -50,7 +51,9 @@ public: ...@@ -50,7 +51,9 @@ public:
/*! \brief Evaluate the expression /*! \brief Evaluate the expression
* *
* \param key where to evaluate the expression * \param k where to evaluate the expression (ignored)
*
* \return the point stored
* *
*/ */
inline point value(const vect_dist_key_dx & k) const inline point value(const vect_dist_key_dx & k) const
......
...@@ -126,29 +126,33 @@ public: ...@@ -126,29 +126,33 @@ public:
#include "Vector/Vector.hpp" #include "Vector/Vector.hpp"
//! stub when library compiled without eigen
template<typename T> template<typename T>
class umfpack_solver class umfpack_solver
{ {
public: public:
//! stub solve
template<typename impl> static Vector<T> solve(const SparseMatrix<T,impl> & A, const Vector<T> & b) template<typename impl> static Vector<T> solve(const SparseMatrix<T,impl> & A, const Vector<T> & b)
{ {
std::cerr << __FILE__ << ":" << __LINE__ << " Error Umfpack only support double precision" << "/n"; std::cerr << __FILE__ << ":" << __LINE__ << " Error Umfpack only support double precision" << "/n";
} }
//! stub solve
void best_solve() void best_solve()
{ {
std::cerr << __FILE__ << ":" << __LINE__ << " Error Umfpack only support double precision" << "/n"; std::cerr << __FILE__ << ":" << __LINE__ << " Error Umfpack only support double precision" << "/n";
} }
}; };
//! stub when library compiled without eigen
template<> template<>
class umfpack_solver<double> class umfpack_solver<double>
{ {
public: public:
//! stub solve
template<unsigned int impl, typename id_type> static Vector<double> solve(SparseMatrix<double,id_type,impl> & A, const Vector<double> & b, size_t opt = UMFPACK_NONE) template<unsigned int impl, typename id_type> static Vector<double> solve(SparseMatrix<double,id_type,impl> & A, const Vector<double> & b, size_t opt = UMFPACK_NONE)
{ {
std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use umfpack you must compile OpenFPM with linear algebra support" << "/n"; std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use umfpack you must compile OpenFPM with linear algebra support" << "/n";
...@@ -158,6 +162,7 @@ public: ...@@ -158,6 +162,7 @@ public:
return x; return x;
} }
//! stub solve
void best_solve() void best_solve()
{ {
std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use umfpack you must compile OpenFPM with linear algebra support" << "/n"; std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use umfpack you must compile OpenFPM with linear algebra support" << "/n";
......
...@@ -13,3 +13,4 @@ ...@@ -13,3 +13,4 @@
#include "util/util_num_unit_tests.hpp" #include "util/util_num_unit_tests.hpp"
#include "PSE/Kernels_unit_tests.hpp" #include "PSE/Kernels_unit_tests.hpp"
#include "Operators/Vector/vector_dist_operators_unit_tests.hpp" #include "Operators/Vector/vector_dist_operators_unit_tests.hpp"
#include "Draw/DrawParticles_unit_tests.hpp"
...@@ -18,8 +18,8 @@ struct is_const_field: std::false_type {}; ...@@ -18,8 +18,8 @@ struct is_const_field: std::false_type {};
/*! \brief is_constant check if a type define a constant field /*! \brief is_constant check if a type define a constant field
* *
* ### Example * ### Example
* \snippet util_num.hpp Constant fields struct definition * \snippet util_num_unit_tests.hpp Constant fields struct definition
* \snippet util_num.hpp Usage of is_const_field * \snippet util_num_unit_tests.hpp Usage of is_const_field