From 2a805e31cd0bc1e5fe89c8f6a03ba22b1691ad2d Mon Sep 17 00:00:00 2001
From: Pietro Incardona <incardon@mpi-cbg.de>
Date: Tue, 3 Jan 2017 22:19:06 +0100
Subject: [PATCH] Adding Draw particles functionalities
---
src/Draw/DrawParticlesGrid.hpp | 75 +++++++++
src/Draw/DrawParticles_unit_tests.hpp | 60 ++++++++
src/Draw/PointIterator.hpp | 144 ++++++++++++++++++
src/FiniteDifference/mul.hpp | 31 +++-
.../Vector/vector_dist_operator_assign.hpp | 7 +
.../vector_dist_operators_apply_kernel.hpp | 3 +-
...dist_operators_apply_kernel_unit_tests.hpp | 3 -
.../vector_dist_operators_extensions.hpp | 11 +-
src/Solvers/umfpack_solver.hpp | 7 +-
src/main.cpp | 1 +
src/util/util_num.hpp | 14 +-
src/util/util_num_unit_tests.hpp | 14 ++
12 files changed, 350 insertions(+), 20 deletions(-)
create mode 100644 src/Draw/DrawParticlesGrid.hpp
create mode 100644 src/Draw/DrawParticles_unit_tests.hpp
create mode 100644 src/Draw/PointIterator.hpp
diff --git a/src/Draw/DrawParticlesGrid.hpp b/src/Draw/DrawParticlesGrid.hpp
new file mode 100644
index 00000000..1767d2ee
--- /dev/null
+++ b/src/Draw/DrawParticlesGrid.hpp
@@ -0,0 +1,75 @@
+/*
+ * 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_ */
diff --git a/src/Draw/DrawParticles_unit_tests.hpp b/src/Draw/DrawParticles_unit_tests.hpp
new file mode 100644
index 00000000..e08d9c38
--- /dev/null
+++ b/src/Draw/DrawParticles_unit_tests.hpp
@@ -0,0 +1,60 @@
+/*
+ * 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_ */
diff --git a/src/Draw/PointIterator.hpp b/src/Draw/PointIterator.hpp
new file mode 100644
index 00000000..6a68625c
--- /dev/null
+++ b/src/Draw/PointIterator.hpp
@@ -0,0 +1,144 @@
+/*
+ * 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_ */
diff --git a/src/FiniteDifference/mul.hpp b/src/FiniteDifference/mul.hpp
index eab2e07e..759e2a64 100644
--- a/src/FiniteDifference/mul.hpp
+++ b/src/FiniteDifference/mul.hpp
@@ -15,10 +15,11 @@
#define HAS_VAL 1
#define HAS_POS_VAL 2
-
+//! Evaluate the constant field function
template <unsigned int, typename T>
struct has_val
{
+ //! evaluate
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";
@@ -26,16 +27,18 @@ struct has_val
}
};
+//! Evaluate the constant field function
template <typename T>
struct has_val<HAS_VAL,T>
{
+ //! evaluate
static decltype(T::val()) call_val()
{
return T::val();
}
};
-
+//! Multiplication expression
template<typename v_expr>
struct const_mul_functor_value
{
@@ -48,12 +51,13 @@ struct const_mul_functor_value
//! sum functor
std::unordered_map<long int,typename last::stype> & cols;
+ //! grid size
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;
- // grid position
+ //! grid position
grid_dist_key_dx<last::dims> & kmap;
//! coefficent
@@ -63,9 +67,21 @@ struct const_mul_functor_value
typename last::stype (& spacing)[last::dims];
/*! \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)
{};
@@ -80,6 +96,11 @@ struct const_mul_functor_value
coeff *= has_val<is_const_field<cfield>::value * 1,cfield>::call_val();
}
+ /*! \brief Get the coefficent
+ *
+ * \return the coefficent
+ *
+ */
typename last::stype getCoeff()
{
return coeff;
diff --git a/src/Operators/Vector/vector_dist_operator_assign.hpp b/src/Operators/Vector/vector_dist_operator_assign.hpp
index b977f4e2..44a3701d 100644
--- a/src/Operators/Vector/vector_dist_operator_assign.hpp
+++ b/src/Operators/Vector/vector_dist_operator_assign.hpp
@@ -8,27 +8,34 @@
#ifndef 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>
struct construct_expression
{
+ //! It return the expression itself
static inline const T & construct(const T & e)
{
return e;
}
};
+//! construct a vector expression from a double
template<>
struct construct_expression<double>
{
+ //! It return an expression from a double
static inline vector_dist_expression<0,double> construct(double e)
{
return vector_dist_expression<0,double>(e);
}
};
+//! construct a vector expression from a float
template<>
struct construct_expression<float>
{
+ //! It return an expression from a float
static inline vector_dist_expression<0,float> construct(const float e)
{
return vector_dist_expression<0,float>(e);
diff --git a/src/Operators/Vector/vector_dist_operators_apply_kernel.hpp b/src/Operators/Vector/vector_dist_operators_apply_kernel.hpp
index 447cfced..e2e1aead 100644
--- a/src/Operators/Vector/vector_dist_operators_apply_kernel.hpp
+++ b/src/Operators/Vector/vector_dist_operators_apply_kernel.hpp
@@ -53,6 +53,7 @@ struct set_zero
}
};
+//! Create a point with all compunent set to zero
template<unsigned int dim, typename T>
struct set_zero<Point<dim,T>>
{
@@ -276,7 +277,7 @@ public:
* \param o1 expression
* \param cl Cell-list
* \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)
diff --git a/src/Operators/Vector/vector_dist_operators_apply_kernel_unit_tests.hpp b/src/Operators/Vector/vector_dist_operators_apply_kernel_unit_tests.hpp
index 51079cac..28a4994f 100644
--- a/src/Operators/Vector/vector_dist_operators_apply_kernel_unit_tests.hpp
+++ b/src/Operators/Vector/vector_dist_operators_apply_kernel_unit_tests.hpp
@@ -51,9 +51,6 @@ template <typename vector,typename Kernel, typename NN_type> bool check_values_a
base2 += vd.template getProp<C>(p);
- if (base1 != base2)
- std::cout << "CAZZO " << base1 << " " << base2 << std::endl;
-
ret &= base1 == base2;
++it2;
diff --git a/src/Operators/Vector/vector_dist_operators_extensions.hpp b/src/Operators/Vector/vector_dist_operators_extensions.hpp
index a3df2298..adf208df 100644
--- a/src/Operators/Vector/vector_dist_operators_extensions.hpp
+++ b/src/Operators/Vector/vector_dist_operators_extensions.hpp
@@ -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 vector involved
+ * \tparam point type of point
*
*/
template<typename point>
class vector_dist_expression<16384,point>
{
+ //! constant point stored
point p;
public:
+ //! vector expression from a constant point
vector_dist_expression(point p)
:p(p)
{}
@@ -50,7 +51,9 @@ public:
/*! \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
diff --git a/src/Solvers/umfpack_solver.hpp b/src/Solvers/umfpack_solver.hpp
index f087871a..bdb85189 100644
--- a/src/Solvers/umfpack_solver.hpp
+++ b/src/Solvers/umfpack_solver.hpp
@@ -126,29 +126,33 @@ public:
#include "Vector/Vector.hpp"
+//! stub when library compiled without eigen
template<typename T>
class umfpack_solver
{
public:
+ //! stub solve
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";
}
+ //! stub solve
void best_solve()
{
std::cerr << __FILE__ << ":" << __LINE__ << " Error Umfpack only support double precision" << "/n";
}
};
-
+//! stub when library compiled without eigen
template<>
class umfpack_solver<double>
{
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)
{
std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use umfpack you must compile OpenFPM with linear algebra support" << "/n";
@@ -158,6 +162,7 @@ public:
return x;
}
+ //! stub solve
void best_solve()
{
std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use umfpack you must compile OpenFPM with linear algebra support" << "/n";
diff --git a/src/main.cpp b/src/main.cpp
index db22b925..30e4ef4b 100755
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -13,3 +13,4 @@
#include "util/util_num_unit_tests.hpp"
#include "PSE/Kernels_unit_tests.hpp"
#include "Operators/Vector/vector_dist_operators_unit_tests.hpp"
+#include "Draw/DrawParticles_unit_tests.hpp"
diff --git a/src/util/util_num.hpp b/src/util/util_num.hpp
index 5d953e47..2e3bbec9 100644
--- a/src/util/util_num.hpp
+++ b/src/util/util_num.hpp
@@ -18,8 +18,8 @@ struct is_const_field: std::false_type {};
/*! \brief is_constant check if a type define a constant field
*
* ### Example
- * \snippet util_num.hpp Constant fields struct definition
- * \snippet util_num.hpp Usage of is_const_field
+ * \snippet util_num_unit_tests.hpp Constant fields struct definition
+ * \snippet util_num_unit_tests.hpp Usage of is_const_field
*
* return true if T::constant_field exist, it does not matter which type is
*
@@ -37,8 +37,8 @@ struct has_val_pos: std::false_type {};
* internal structure with attributes
*
* ### Example
- * \snippet util.hpp Declaration of struct with attributes and without
- * \snippet util.hpp Check has_attributes
+ * \snippet util_test.hpp Declaration of struct with attributes and without
+ * \snippet util_test.hpp Check has_attributes
*
* return true if T::attributes::name[0] is a valid expression
* and produce a defined type
@@ -54,8 +54,7 @@ template<typename T, typename Sfinae = void>
struct is_testing: std::false_type {};
-/*! \brief has_attributes check if a type has defined an
- * internal structure with attributes
+/*! \brief is_testing check if a struct T has testing member defined
*
* ### Example
* \snippet util_num_unit_tests.hpp Is on test struct definition
@@ -83,12 +82,15 @@ struct is_testing<T, typename Void<typename T::testing >::type> : std::true_type
template<typename T, unsigned int dims, typename stype, typename decomposition, bool stub=is_testing<T>::value>
struct stub_or_real
{
+ //! switch type if we are on testing or not
typedef grid_dist_testing<dims> type;
};
+//! Case when we are not on testing
template<typename T, unsigned int dims, typename stype, typename decomposition>
struct stub_or_real<T,dims,stype,decomposition,false>
{
+ //! switch type if we are on testing or not
typedef grid_dist_id<dims,stype,scalar<size_t>,decomposition> type;
};
diff --git a/src/util/util_num_unit_tests.hpp b/src/util/util_num_unit_tests.hpp
index 0683b787..27a305a6 100644
--- a/src/util/util_num_unit_tests.hpp
+++ b/src/util/util_num_unit_tests.hpp
@@ -12,21 +12,30 @@
//! [Constant fields struct definition]
+//! define a constant field
struct anyname_field
{
+ //! define that is a constant field
typedef void const_field;
+ //! Evaluate the constant field (in this case always return 1.0)
float val() {return 1.0;}
};
+//! define a non-constant (in space) field
struct anyname_field_with_pos
{
+ //! It define that is a constant field
typedef void const_field;
+
+ //! It define that is not constant in space
typedef void with_position;
+ //! Evaluate the field in one point
float val_pos(grid_key_dx<2> & pos) {return 1.0;}
};
+//! stub field
struct no_field
{
};
@@ -35,15 +44,20 @@ struct no_field
//! [Is on test struct definition]
+//! define that we are on testing mode
struct on_test
{
+ //! specify testing mode
typedef void testing;
};
+//! Not on testing mode
struct not_on_test
{
};
+//! [Is on test struct definition]
+
BOOST_AUTO_TEST_SUITE( util_num_suite )
BOOST_AUTO_TEST_CASE( util_num )
--
GitLab