diff --git a/Makefile.am b/Makefile.am
index 7c618d3e326a45c8f6719691f58bde9a5a1c4c26..bad0185690207b20c536ceace77718261a0f6833 100755
--- a/Makefile.am
+++ b/Makefile.am
@@ -1,3 +1,7 @@
SUBDIRS = src
-bin_PROGRAMS =
\ No newline at end of file
+bin_PROGRAMS =
+
+test:
+ cd src && make test
+
diff --git a/src/Makefile.am b/src/Makefile.am
index 6a3c9c8c2bd412793e0398dde67f9810648f1115..802c23f261255c7172f16e112cf6b93025b64336 100755
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -10,3 +10,7 @@ nobase_include_HEADERS = PSE/Kernels.hpp PSE/Kernels_test_util.hpp
.cu.o :
$(NVCC) $(NVCCFLAGS) -o $@ -c $<
+
+test: numerics
+ source $(HOME)/openfpm_vars && cd .. && mpirun -np 2 ./src/numerics && mpirun -np 3 ./src/numerics && mpirun -np 4 ./src/numerics
+
diff --git a/src/Operators/Vector/vector_dist_operators.hpp b/src/Operators/Vector/vector_dist_operators.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..64918f017764f97734e95a99206b37a511dbd575
--- /dev/null
+++ b/src/Operators/Vector/vector_dist_operators.hpp
@@ -0,0 +1,1023 @@
+/*
+ * vector_dist_operators.hpp
+ *
+ * Created on: Jun 11, 2016
+ * Author: i-bird
+ */
+
+#ifndef OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_HPP_
+#define OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_HPP_
+
+
+#define VECT_SUM 1
+#define VECT_SUB 2
+#define VECT_MUL 3
+#define VECT_DIV 4
+#define VECT_NORM 5
+#define VECT_NORM2 6
+#define VECT_APPLYKER 7
+
+/*! \brief is_expression check if a type is simple a type or is just an encapsulation of an expression
+ *
+ * return true if T::is_expression is a valid expression
+ *
+ */
+
+template<typename ObjType, typename Sfinae = void>
+struct is_expression: std::false_type {};
+
+template<typename ObjType>
+struct is_expression<ObjType, typename Void<typename ObjType::is_expression>::type> : std::true_type
+{};
+
+template<typename exp, bool is_exp = is_expression<exp>::value>
+struct apply_kernel_rtype
+{
+ typedef typename exp::orig_type rtype;
+};
+
+template<typename exp>
+struct apply_kernel_rtype<exp,false>
+{
+ typedef exp rtype;
+};
+
+
+/*! \brief Apply the kernel to particle differently that is a number or is an expression
+ *
+ *
+ */
+template<typename T, typename vector, typename exp,typename NN_type, typename Kernel, typename rtype, bool is_exp=is_expression<T>::value>
+struct apply_kernel_is_number_or_expression
+{
+ inline static rtype apply(const vector & vd, NN_type & cl, const exp & v_exp, const vect_dist_key_dx & key, Kernel & lker)
+ {
+ // accumulator
+ rtype pse = static_cast<typename rtype::coord_type>(0.0);
+
+ // position of particle p
+ Point<vector::dims,typename vector::stype> p = vd.getPos(key);
+
+ // property of the particle x
+ rtype prp_p = v_exp.value(key);
+
+ // Get the neighborhood of the particle
+ auto NN = cl.template getNNIterator<NO_CHECK>(cl.getCell(p));
+ while(NN.isNext())
+ {
+ auto nnp = NN.get();
+
+ // Calculate contribution given by the kernel value at position p,
+ // given by the Near particle, exclude itself
+ if (nnp != key.getKey())
+ {
+ // property of the particle x
+ rtype prp_q = v_exp.value(nnp);
+
+ // position of the particle q
+ Point<vector::dims,typename vector::stype> q = vd.getPos(nnp);
+
+ for (size_t i = 0 ; i < T::nvals ; i++)
+ {
+ float val = lker.value(p,q,prp_p.value(i),prp_q.value(i));
+
+ // W(x-y)
+ pse.value(i) += val;
+ }
+ }
+
+ // Next particle
+ ++NN;
+ }
+
+ return pse;
+ }
+};
+
+template<typename T, typename vector, typename exp,typename NN_type, typename Kernel, typename rtype>
+struct apply_kernel_is_number_or_expression<T,vector,exp,NN_type,Kernel,rtype,false>
+{
+
+ /*! \brief Apply the kernel
+ *
+ * \param cl Neighborhood of particles
+ * \param v_exp vector expression
+ * \param key particle id
+ * \param lker kernel function
+ *
+ */
+ inline static rtype apply(const vector & vd, NN_type & cl, const exp & v_exp, const vect_dist_key_dx & key, const Kernel & lker)
+ {
+ // accumulator
+ rtype pse = 0;
+
+ // Get f(x) at the position of the particle
+ rtype & prp_p = v_exp.value(key);
+
+ // position of particle p
+ auto & p = vd.getPos(key);
+
+ // Get the neighborhood of the particle
+ auto NN = cl.template getNNIterator<NO_CHECK>(cl.getCell(p));
+ while(NN.isNext())
+ {
+ auto nnp = NN.get();
+
+ // Calculate contribution given by the kernel value at position p,
+ // given by the Near particle, exclude itself
+ if (nnp != key.getKey())
+ {
+ // Get f(x) at the position of the particle
+ rtype & prp_q = vd.getPos(nnp);
+
+ // position of particle q
+ auto & q = v_exp.value_pos(key);
+
+ // W(x-y)
+ auto ker = lker.value(p,q,prp_p,prp_q);
+
+ pse += ker;
+ }
+
+ // Next particle
+ ++NN;
+ }
+
+ return pse;
+ }
+};
+
+/*! \brief has_init check if a type has defined a
+ * method called init
+ *
+ *
+ * return true if T::init() is a valid expression (function pointers)
+ * and produce a defined type
+ *
+ */
+
+template<typename ObjType, typename Sfinae = void>
+struct has_init: std::false_type {};
+
+template<typename ObjType>
+struct has_init<ObjType, typename Void<typename ObjType::has_init>::type> : std::true_type
+{};
+
+/*! \brief Call the init function if needed
+ *
+ * \param r_exp expression
+ *
+ */
+template <typename T, bool has_init = has_init<T>::value >
+struct call_init_if_needed
+{
+ static inline void call(T & r_exp)
+ {
+ r_exp.init();
+ }
+};
+
+template <typename T>
+struct call_init_if_needed<T,false>
+{
+ static inline void call(T & r_exp)
+ {
+ }
+};
+
+
+
+/*! \brief Unknown operation specialization
+ *
+ * \tparam exp1 expression1
+ * \tparam exp2 expression2
+ *
+ */
+template <typename exp1, typename exp2, unsigned int op>
+class vector_dist_expression_op
+{
+
+};
+
+/*! \brief Sum operation
+ *
+ * \tparam exp1 expression1
+ * \tparam exp2 expression2
+ *
+ */
+template <typename exp1, typename exp2>
+class vector_dist_expression_op<exp1,exp2,VECT_SUM>
+{
+ const exp1 o1;
+ const exp2 o2;
+
+public:
+
+ inline vector_dist_expression_op(const exp1 & o1, const exp2 & o2)
+ :o1(o1),o2(o2)
+ {}
+
+ /*! \brief Evaluate the expression
+ *
+ * \param key where to evaluate the expression
+ *
+ */
+ template<typename r_type=typename std::remove_reference<decltype(o1.value(vect_dist_key_dx(0)) + o2.value(vect_dist_key_dx(0)))>::type > inline r_type value(const vect_dist_key_dx & key) const
+ {
+ return o1.value(key) + o2.value(key);
+ }
+
+};
+
+/*! \brief Subtraction operation
+ *
+ * \tparam exp1 expression1
+ * \tparam exp2 expression2
+ *
+ */
+template <typename exp1, typename exp2>
+class vector_dist_expression_op<exp1,exp2,VECT_SUB>
+{
+ const exp1 o1;
+ const exp2 o2;
+
+public:
+
+ inline vector_dist_expression_op(const exp1 & o1, const exp2 & o2)
+ :o1(o1),o2(o2)
+ {}
+
+ /*! \brief Evaluate the expression
+ *
+ * \param key where to evaluate the expression
+ *
+ */
+ template<typename r_type=typename std::remove_reference<decltype(o1.value(vect_dist_key_dx(0)) - o2.value(vect_dist_key_dx(0)))>::type > inline r_type value(const vect_dist_key_dx & key) const
+ {
+ return o1.value(key) - o2.value(key);
+ }
+};
+
+/*! \brief Multiplication operation
+ *
+ * \tparam exp1 expression1
+ * \tparam exp2 expression2
+ *
+ */
+template <typename exp1, typename exp2>
+class vector_dist_expression_op<exp1,exp2,VECT_MUL>
+{
+ const exp1 o1;
+ const exp2 o2;
+
+public:
+
+ vector_dist_expression_op(const exp1 & o1, const exp2 & o2)
+ :o1(o1),o2(o2)
+ {}
+
+ /*! \brief Evaluate the expression
+ *
+ * \param key where to evaluate the expression
+ *
+ */
+ template<typename r_type=typename std::remove_reference<decltype(o1.value(vect_dist_key_dx(0)) * o2.value(vect_dist_key_dx(0)))>::type > inline r_type value(const vect_dist_key_dx & key) const
+ {
+ return o1.value(key) * o2.value(key);
+ }
+};
+
+/*! \brief Division operation
+ *
+ * \tparam exp1 expression1
+ * \tparam exp2 expression2
+ *
+ */
+template <typename exp1, typename exp2>
+class vector_dist_expression_op<exp1,exp2,VECT_DIV>
+{
+ const exp1 o1;
+ const exp2 o2;
+
+public:
+
+ vector_dist_expression_op(const exp1 & o1, const exp2 & o2)
+ :o1(o1),o2(o2)
+ {}
+
+ /*! \brief Evaluate the expression
+ *
+ * \param key where to evaluate the expression
+ *
+ */
+ template<typename r_type=typename std::remove_reference<decltype(o1.value(vect_dist_key_dx(0)) / o2.value(vect_dist_key_dx(0)))>::type > inline r_type value(const vect_dist_key_dx & key) const
+ {
+ return o1.value(key) / o2.value(key);
+ }
+};
+
+/*! \brief norm operation
+ *
+ * \tparam exp1 expression1
+ * \tparam exp2 expression2
+ *
+ */
+template <typename exp1>
+class vector_dist_expression_op<exp1,void,VECT_NORM>
+{
+ const exp1 o1;
+
+public:
+
+ vector_dist_expression_op(const exp1 & o1)
+ :o1(o1)
+ {}
+
+ /*! \brief Evaluate the expression
+ *
+ * \param key where to evaluate the expression
+ *
+ */
+ template<typename r_type=typename std::remove_reference<decltype(norm(o1.value(vect_dist_key_dx(0))))>::type > inline r_type value(const vect_dist_key_dx & key) const
+ {
+ return norm(o1.value(key));
+ }
+};
+
+/*! \brief Apply kernel operation
+ *
+ * \tparam exp1 expression1
+ * \tparam NN list
+ *
+ */
+template <typename exp1,typename vector_type>
+class vector_dist_expression_op<exp1,vector_type,VECT_APPLYKER>
+{
+ typedef typename boost::mpl::at<vector_type,boost::mpl::int_<0>>::type NN;
+ typedef typename boost::mpl::at<vector_type,boost::mpl::int_<1>>::type Kernel;
+ typedef typename boost::mpl::at<vector_type,boost::mpl::int_<2>>::type vector_orig;
+
+ const exp1 o1;
+ NN & cl;
+ Kernel & ker;
+ const vector_orig & vd;
+
+public:
+
+ vector_dist_expression_op(const exp1 & o1, NN & cl, Kernel & ker, const vector_orig & vd)
+ :o1(o1),cl(cl),ker(ker),vd(vd)
+ {}
+
+ /*! \brief Evaluate the expression
+ *
+ * \param key where to evaluate the expression
+ *
+ */
+ template<typename r_type=typename apply_kernel_rtype<decltype(o1.value(vect_dist_key_dx(0)))>::rtype > inline r_type value(const vect_dist_key_dx & key) const
+ {
+ typedef typename apply_kernel_rtype<decltype(o1.value(vect_dist_key_dx(0)))>::rtype rtype;
+
+ return apply_kernel_is_number_or_expression<decltype(o1.value(key)),vector_orig,exp1,NN,Kernel,rtype>::apply(vd,cl,o1,key,ker);
+ }
+};
+
+/*! \brief Main class that encapsulate a vector properties
+ *
+ * \tparam prp property involved
+ * \tparam vector involved
+ *
+ */
+template<unsigned int prp, typename vector>
+class vector_dist_expression
+{
+ vector & v;
+
+public:
+
+ vector_dist_expression(vector & v)
+ :v(v)
+ {}
+
+ /*! \brief Evaluate the expression
+ *
+ * \param key where to evaluate the expression
+ *
+ */
+ inline auto value(const vect_dist_key_dx & k) const -> decltype(v.template getProp<prp>(k))
+ {
+ return v.template getProp<prp>(k);
+ }
+
+
+ /*! \brief Fill the vector property with the evaluated expression
+ *
+ * \param v_exp expression to evaluate
+ *
+ */
+ template<typename exp1, typename exp2, unsigned int op> vector & operator=(const vector_dist_expression_op<exp1,exp2,op> & v_exp)
+ {
+ auto it = v.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ auto exp = v_exp.value(key);
+ call_init_if_needed<decltype(exp)>::call(exp);
+ v.template getProp<prp>(key) = v_exp.value(key);
+
+ ++it;
+ }
+
+ return v;
+ }
+
+ /*! \brief Fill the vector property with the double
+ *
+ * \param d value to fill
+ *
+ */
+ vector & operator=(double d)
+ {
+ auto it = v.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ v.template getProp<prp>(key) = d;
+
+ ++it;
+ }
+
+ return v;
+ }
+};
+
+/*! \brief Main class that encapsulate a double constant
+ *
+ * \param prp no meaning
+ *
+ */
+template<unsigned int prp>
+class vector_dist_expression<prp,double>
+{
+ double d;
+
+public:
+
+ inline vector_dist_expression(double & d)
+ :d(d)
+ {}
+
+ /*! \brief Evaluate the expression
+ *
+ * It just return the velue set in the constructor
+ *
+ */
+ inline double value(const vect_dist_key_dx & k) const
+ {
+ return d;
+ }
+};
+
+
+/* \brief sum two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<unsigned int p1, unsigned int p2, typename v1, typename v2>
+inline vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression<p2,v2>,VECT_SUM>
+operator+(const vector_dist_expression<p1,v1> & va, const vector_dist_expression<p2,v2> & vb)
+{
+ vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression<p2,v2>,VECT_SUM> exp_sum(va,vb);
+
+ return exp_sum;
+}
+
+/* \brief sum two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<typename exp1 , typename exp2, unsigned int op1, unsigned int prp1, typename v1>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<prp1,v1>,VECT_SUM>
+operator+(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression<prp1,v1> & vb)
+{
+ vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<prp1,v1>,VECT_SUM> exp_sum(va,vb);
+
+ return exp_sum;
+}
+
+/* \brief sum two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<typename exp1 , typename exp2, unsigned int op1, unsigned int prp1, typename v1>
+inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression_op<exp1,exp2,op1>,VECT_SUM>
+operator+(const vector_dist_expression<prp1,v1> & va, const vector_dist_expression_op<exp1,exp2,op1> & vb)
+{
+ vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression_op<exp1,exp2,op1>,VECT_SUM> exp_sum(va,vb);
+
+ return exp_sum;
+}
+
+/* \brief sum two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<typename exp1 , typename exp2, unsigned int op1, typename exp3 , typename exp4, unsigned int op2>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,VECT_SUM>
+operator+(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression_op<exp3,exp4,op2> & vb)
+{
+ vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,VECT_SUM> exp_sum(va,vb);
+
+ return exp_sum;
+}
+
+/* \brief sum two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<unsigned int prp1 , typename v1>
+inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,VECT_SUM>
+operator+(const vector_dist_expression<prp1,v1> & va, double d)
+{
+ vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,VECT_SUM> exp_sum(va,vector_dist_expression<0,double>(d));
+
+ return exp_sum;
+}
+
+/* \brief sum two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<unsigned int prp1 , typename v1>
+inline vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<prp1,v1>,VECT_SUM>
+operator+(double d, const vector_dist_expression<prp1,v1> & vb)
+{
+ vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<prp1,v1>,VECT_SUM> exp_sum(vector_dist_expression<0,double>(d),vb);
+
+ return exp_sum;
+}
+
+/* \brief sum two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<typename exp1 , typename exp2, unsigned int op1>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,VECT_SUM>
+operator+(const vector_dist_expression_op<exp1,exp2,op1> & va, double d)
+{
+ vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,VECT_SUM> exp_sum(va,vector_dist_expression<0,double>(d));
+
+ return exp_sum;
+}
+
+
+/* \brief subtract two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<unsigned int p1, unsigned int p2, typename v1, typename v2>
+inline vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression<p2,v2>,VECT_SUB>
+operator-(const vector_dist_expression<p1,v1> & va, const vector_dist_expression<p2,v2> & vb)
+{
+ vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression<p2,v2>,VECT_SUB> exp_sum(va,vb);
+
+ return exp_sum;
+}
+
+
+/* \brief subtract two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<typename exp1, typename exp2, unsigned int op1, unsigned int p2, typename v2>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<p2,v2>,VECT_SUB>
+operator-(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression<p2,v2> & vb)
+{
+ vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<p2,v2>,VECT_SUB> exp_sum(va,vb);
+
+ return exp_sum;
+}
+
+/* \brief subtract two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<typename exp1, typename exp2, unsigned int op1, unsigned int p2, typename v2>
+inline vector_dist_expression_op<vector_dist_expression<p2,v2>,vector_dist_expression_op<exp1,exp2,op1>,VECT_SUB>
+operator-(const vector_dist_expression<p2,v2> & va, const vector_dist_expression_op<exp1,exp2,op1> & vb)
+{
+ vector_dist_expression_op<vector_dist_expression<p2,v2>, vector_dist_expression_op<exp1,exp2,op1>,VECT_SUB> exp_sum(va,vb);
+
+ return exp_sum;
+}
+
+/* \brief subtract two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<typename exp1, typename exp2, unsigned int op1, typename exp3, typename exp4, unsigned int op2>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,VECT_SUB>
+operator-(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression_op<exp3,exp4,op2> & vb)
+{
+ vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,VECT_SUB> exp_sum(va,vb);
+
+ return exp_sum;
+}
+
+/* \brief subtract two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<unsigned int prp1, typename v1>
+inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,VECT_SUB>
+operator-(const vector_dist_expression<prp1,v1> & va, double d)
+{
+ vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,VECT_SUB> exp_sum(va,vector_dist_expression<0,double>(d));
+
+ return exp_sum;
+}
+
+/* \brief subtract two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<unsigned int prp1, typename v1>
+inline vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<prp1,v1>,VECT_SUB>
+operator-(double d, const vector_dist_expression<prp1,v1> & vb)
+{
+ vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<prp1,v1>,VECT_SUB> exp_sum(vector_dist_expression<0,double>(d),vb);
+
+ return exp_sum;
+}
+
+/* \brief Multiply two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<unsigned int p2, typename v2>
+inline vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<p2,v2>,VECT_MUL>
+operator*(double d, const vector_dist_expression<p2,v2> & vb)
+{
+ vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<p2,v2>,VECT_MUL> exp_sum(vector_dist_expression<0,double>(d),vb);
+
+ return exp_sum;
+}
+
+/* \brief Multiply two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<unsigned int p2, typename v2>
+inline vector_dist_expression_op<vector_dist_expression<p2,v2>,vector_dist_expression<0,double>,VECT_MUL>
+operator*(const vector_dist_expression<p2,v2> & va, double d)
+{
+ vector_dist_expression_op<vector_dist_expression<p2,v2>,vector_dist_expression<0,double>,VECT_MUL> exp_sum(va,vector_dist_expression<0,double>(d));
+
+ return exp_sum;
+}
+
+/* \brief Multiply two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<unsigned int p1, typename v1,unsigned int p2, typename v2>
+inline vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression<p2,v2>,VECT_MUL>
+operator*(const vector_dist_expression<p1,v1> & va, const vector_dist_expression<p2,v2> & vb)
+{
+ vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression<p2,v2>,VECT_MUL> exp_sum(va,vb);
+
+ return exp_sum;
+}
+
+/* \brief Multiply two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<unsigned int p1, typename v1, typename exp1, typename exp2, unsigned int op1>
+inline vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression_op<exp1,exp2,op1>,VECT_MUL>
+operator*(const vector_dist_expression<p1,v1> & va, const vector_dist_expression_op<exp1,exp2,op1> & vb)
+{
+ vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression_op<exp1,exp2,op1>,VECT_MUL> exp_sum(va,vb);
+
+ return exp_sum;
+}
+
+/* \brief Multiply two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<unsigned int p1, typename v1, typename exp1, typename exp2, unsigned int op1>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<p1,v1>,VECT_MUL>
+operator*(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression<p1,v1> & vb)
+{
+ vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<p1,v1>,VECT_MUL> exp_sum(va,vb);
+
+ return exp_sum;
+}
+
+/* \brief Multiply two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<typename exp1, typename exp2, unsigned int op1, typename exp3 , typename exp4, unsigned int op2>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,VECT_MUL>
+operator*(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression_op<exp3,exp4,op2> & vb)
+{
+ vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,VECT_MUL> exp_sum(va,vb);
+
+ return exp_sum;
+}
+
+/* \brief Divide two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<typename exp1, typename exp2, unsigned int op1>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,VECT_DIV>
+operator/(const vector_dist_expression_op<exp1,exp2,op1> & va, double d)
+{
+ vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,VECT_DIV> exp_sum(va,vector_dist_expression<0,double>(d));
+
+ return exp_sum;
+}
+
+
+/* \brief Divide two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<typename exp1, typename exp2, unsigned int op1>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,VECT_DIV>
+operator/(double d, const vector_dist_expression_op<exp1,exp2,op1> & va)
+{
+ vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,VECT_DIV> exp_sum(vector_dist_expression<0,double>(d),va);
+
+ return exp_sum;
+}
+
+/* \brief Divide two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<unsigned int prp1, typename v1>
+inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,VECT_DIV>
+operator/(const vector_dist_expression<prp1,v1> & va, double d)
+{
+ vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,VECT_DIV> exp_sum(va,vector_dist_expression<0,double>(d));
+
+ return exp_sum;
+}
+
+/* \brief Divide two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<unsigned int prp1, typename v1>
+inline vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<prp1,v1>,VECT_DIV>
+operator/(double d, const vector_dist_expression<prp1,v1> & va)
+{
+ vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<prp1,v1>,VECT_DIV> exp_sum(vector_dist_expression<0,double>(d),va);
+
+ return exp_sum;
+}
+
+/* \brief Divide two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<unsigned int prp1, typename v1, unsigned int prp2, typename v2>
+inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<prp2,v2>,VECT_DIV>
+operator/(const vector_dist_expression<prp1,v1> & va, const vector_dist_expression<prp2,v2> & vb)
+{
+ vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<prp2,v2>,VECT_DIV> exp_sum(va,vb);
+
+ return exp_sum;
+}
+
+/* \brief Divide two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<unsigned int prp1, typename v1, typename exp1,typename exp2, unsigned int op1>
+inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression_op<exp1,exp2,op1>,VECT_DIV>
+operator/(const vector_dist_expression<prp1,v1> & va, const vector_dist_expression_op<exp1,exp2,op1> & vb)
+{
+ vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression_op<exp1,exp2,op1>,VECT_DIV> exp_sum(va,vb);
+
+ return exp_sum;
+}
+
+/* \brief Divide two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<unsigned int prp1, typename v1, typename exp1,typename exp2, unsigned int op1>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<prp1,v1>,VECT_DIV>
+operator/(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression<prp1,v1> & vb)
+{
+ vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<prp1,v1>,VECT_DIV> exp_sum(va,vb);
+
+ return exp_sum;
+}
+
+/* \brief Divide two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<typename exp1,typename exp2, unsigned int op1, typename exp3, typename exp4, unsigned int op2>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,VECT_DIV>
+operator/(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression_op<exp3,exp4,op2> & vb)
+{
+ vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,VECT_DIV> exp_sum(va,vb);
+
+ return exp_sum;
+}
+/////////// NORM operator ///////////////////////
+
+/* \brief Divide two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<typename exp1, typename exp2, unsigned int op1>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,void,VECT_NORM>
+norm(const vector_dist_expression_op<exp1,exp2,op1> & va)
+{
+ vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,void,VECT_NORM> exp_sum(va);
+
+ return exp_sum;
+}
+
+
+/* \brief Divide two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<typename exp1, typename exp2, unsigned int op1>
+inline vector_dist_expression_op<vector_dist_expression<0,double>,void,VECT_NORM>
+norm(double d)
+{
+ vector_dist_expression_op<vector_dist_expression<0,double>,void,VECT_NORM> exp_sum(vector_dist_expression<0,double>(d));
+
+ return exp_sum;
+}
+
+/* \brief Divide two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<unsigned int prp1, typename v1>
+inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,void,VECT_NORM>
+norm(const vector_dist_expression<prp1,v1> & va)
+{
+ vector_dist_expression_op<vector_dist_expression<prp1,v1>,void,VECT_NORM> exp_sum(va);
+
+ return exp_sum;
+}
+
+
+///////////////////////////////////////////// Apply kernel operator ////
+////////////////////////////////////////////////////////////////////////
+
+/* \brief Divide two distributed vector expression
+ *
+ * \param va vector expression one
+ * \param vb vector expression two
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<typename exp1, typename exp2, unsigned int op1, typename NN, typename Kernel, typename vector_type>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,boost::mpl::vector<NN,Kernel,vector_type>,VECT_APPLYKER>
+applyKernel(const vector_dist_expression_op<exp1,exp2,op1> & va, vector_type & vd, NN & cl, Kernel & ker)
+{
+ vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,boost::mpl::vector<NN,Kernel,vector_type>,VECT_APPLYKER> exp_sum(va,cl,ker,vd);
+
+ return exp_sum;
+}
+
+
+
+#endif /* OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_HPP_ */
diff --git a/src/Operators/Vector/vector_dist_operators_unit_tests.hpp b/src/Operators/Vector/vector_dist_operators_unit_tests.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ef8f9b24dd64a3a08e499649753eb7d869f8f965
--- /dev/null
+++ b/src/Operators/Vector/vector_dist_operators_unit_tests.hpp
@@ -0,0 +1,751 @@
+/*
+ * vector_dist_operators_unit_tests.hpp
+ *
+ * Created on: Jun 11, 2016
+ * Author: i-bird
+ */
+
+#ifndef OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_UNIT_TESTS_HPP_
+#define OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_UNIT_TESTS_HPP_
+
+constexpr int A = 0;
+constexpr int B = 1;
+constexpr int C = 2;
+
+constexpr int PA = 3;
+constexpr int PB = 4;
+constexpr int PC = 5;
+
+//////////////////// Here we define all the function to checl the operators
+
+template <unsigned int prp, typename vector> bool check_values(const vector & v,float a)
+{
+ bool ret = true;
+ auto it = v.getDomainIterator();
+
+ while (it.isNext())
+ {
+ ret &= v.template getProp<prp>(it.get()) == a;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+template <typename vector> bool check_values_complex_expr(const vector & vd)
+{
+ bool ret = true;
+ auto it = vd.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ float base1 = vd.template getProp<B>(key) + 2.0 + vd.template getProp<B>(key) - 2.0*vd.template getProp<C>(key) / 5.0;
+ float base2 = vd.template getProp<A>(key);
+
+ ret &= base1 == base2;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_sum(const vector & vd, double d)
+{
+ bool ret = true;
+ auto it = vd.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ rtype base1 = vd.template getProp<B>(key) + d;
+ rtype base2 = vd.template getProp<A>(key);
+
+ ret &= base1 == base2;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_sum(const vector & vd1, const vector & vd2)
+{
+ bool ret = true;
+ auto it = vd1.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ rtype base1 = vd1.template getProp<B>(key) + vd2.template getProp<C>(key);
+ rtype base2 = vd1.template getProp<A>(key);
+
+ ret &= base1 == base2;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_sum_3(const vector & vd1)
+{
+ bool ret = true;
+ auto it = vd1.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ rtype base1 = vd1.template getProp<B>(key) + vd1.template getProp<C>(key) + vd1.template getProp<B>(key);
+ rtype base2 = vd1.template getProp<A>(key);
+
+ ret &= base1 == base2;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_sum_4(const vector & vd1)
+{
+ bool ret = true;
+ auto it = vd1.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ rtype base1 = vd1.template getProp<B>(key) + vd1.template getProp<C>(key) + vd1.template getProp<B>(key) + vd1.template getProp<C>(key);
+ rtype base2 = vd1.template getProp<A>(key);
+
+ ret &= base1 == base2;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_sub(const vector & vd, double d)
+{
+ bool ret = true;
+ auto it = vd.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ rtype base1 = vd.template getProp<B>(key) - d;
+ rtype base2 = vd.template getProp<A>(key);
+
+ ret &= base1 == base2;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_sub(double d, const vector & vd)
+{
+ bool ret = true;
+ auto it = vd.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ rtype base1 = d - vd.template getProp<B>(key);
+ rtype base2 = vd.template getProp<A>(key);
+
+ ret &= base1 == base2;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_sub(const vector & vd1, const vector & vd2)
+{
+ bool ret = true;
+ auto it = vd1.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ rtype base1 = vd1.template getProp<C>(key) - vd2.template getProp<B>(key);
+ rtype base2 = vd1.template getProp<A>(key);
+
+ ret &= base1 == base2;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_sub_31(const vector & vd1)
+{
+ bool ret = true;
+ auto it = vd1.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ rtype base1 = vd1.template getProp<B>(key) - (vd1.template getProp<C>(key) + vd1.template getProp<B>(key));
+ rtype base2 = vd1.template getProp<A>(key);
+
+ ret &= base1 == base2;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_sub_32(const vector & vd1)
+{
+ bool ret = true;
+ auto it = vd1.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ rtype base1 = (vd1.template getProp<C>(key) + vd1.template getProp<B>(key)) - vd1.template getProp<B>(key);
+ rtype base2 = vd1.template getProp<A>(key);
+
+ ret &= base1 == base2;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_sub_4(const vector & vd1)
+{
+ bool ret = true;
+ auto it = vd1.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ rtype base1 = (vd1.template getProp<C>(key) + vd1.template getProp<B>(key)) - (vd1.template getProp<C>(key) + vd1.template getProp<B>(key));
+ rtype base2 = vd1.template getProp<A>(key);
+
+ ret &= base1 == base2;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_mul(const vector & vd, double d)
+{
+ bool ret = true;
+ auto it = vd.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ rtype base1 = vd.template getProp<B>(key) * d;
+ rtype base2 = vd.template getProp<A>(key);
+
+ ret &= base1 == base2;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_mul(const vector & vd1, const vector & vd2)
+{
+ bool ret = true;
+ auto it = vd1.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ rtype base1 = vd1.template getProp<C>(key) * vd2.template getProp<B>(key);
+ rtype base2 = vd1.template getProp<A>(key);
+
+ ret &= base1 == base2;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_mul_3(const vector & vd1)
+{
+ bool ret = true;
+ auto it = vd1.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ rtype base1 = vd1.template getProp<B>(key) * (vd1.template getProp<B>(key) + vd1.template getProp<C>(key));
+ rtype base2 = vd1.template getProp<A>(key);
+
+ ret &= base1 == base2;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_mul_4(const vector & vd1)
+{
+ bool ret = true;
+ auto it = vd1.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ rtype base1 = (vd1.template getProp<B>(key) + vd1.template getProp<C>(key)) * (vd1.template getProp<B>(key) + vd1.template getProp<C>(key));
+ rtype base2 = vd1.template getProp<A>(key);
+
+ ret &= base1 == base2;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_div(const vector & vd, double d)
+{
+ bool ret = true;
+ auto it = vd.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ rtype base1 = vd.template getProp<B>(key) / d;
+ rtype base2 = vd.template getProp<A>(key);
+
+ ret &= base1 == base2;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_div(double d, const vector & vd)
+{
+ bool ret = true;
+ auto it = vd.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ rtype base1 = d / vd.template getProp<B>(key);
+ rtype base2 = vd.template getProp<A>(key);
+
+ ret &= base1 == base2;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_div(const vector & vd1, const vector & vd2)
+{
+ bool ret = true;
+ auto it = vd1.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ rtype base1 = vd1.template getProp<C>(key) / vd2.template getProp<B>(key);
+ rtype base2 = vd1.template getProp<A>(key);
+
+ ret &= base1 == base2;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_div_31(const vector & vd1)
+{
+ bool ret = true;
+ auto it = vd1.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ rtype base1 = vd1.template getProp<B>(key) / (vd1.template getProp<B>(key) + vd1.template getProp<C>(key));
+ rtype base2 = vd1.template getProp<A>(key);
+
+ ret &= base1 == base2;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_div_32(const vector & vd1)
+{
+ bool ret = true;
+ auto it = vd1.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ rtype base1 = (vd1.template getProp<C>(key) + vd1.template getProp<B>(key)) / vd1.template getProp<B>(key);
+ rtype base2 = vd1.template getProp<A>(key);
+
+ ret &= base1 == base2;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_div_4(const vector & vd1)
+{
+ bool ret = true;
+ auto it = vd1.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ rtype base1 = (vd1.template getProp<B>(key) + vd1.template getProp<C>(key)) / (vd1.template getProp<B>(key) + vd1.template getProp<C>(key));
+ rtype base2 = vd1.template getProp<A>(key);
+
+ ret &= base1 == base2;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+template <typename vector> bool check_values_scal_norm_dist(vector & vd)
+{
+ bool ret = true;
+ auto it = vd.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ float base1 = vd.template getProp<PB>(key) * vd.template getProp<PC>(key) + norm(vd.template getProp<PC>(key) + vd.template getProp<PB>(key)) + distance(vd.template getProp<PC>(key),vd.template getProp<PB>(key));
+ float base2 = vd.template getProp<A>(key);
+
+ ret &= base1 == base2;
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+
+ return ret;
+}
+
+template <typename vector> void fill_values(const vector & v)
+{
+ auto it = v.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto p = it.get();
+
+ v.template getProp<A>(p) = p.getKey()+1;
+ v.template getProp<B>(p) = 2.0*p.getKey()+1;
+ v.template getProp<C>(p) = 3.0*p.getKey()+1;
+
+ for (size_t k = 0 ; k < 3 ; k++)
+ {
+ v.template getProp<PA>(p)[k] = p.getKey()+1+k;
+ v.template getProp<PB>(p)[k] = 2.0*p.getKey()+1+k;
+ v.template getProp<PC>(p)[k] = 3.0*p.getKey()+1+k;
+ }
+
+ ++it;
+ }
+}
+
+//! Exponential kernel
+struct exp_kernel
+{
+ float var;
+
+ exp_kernel(float var)
+ :var(var)
+ {}
+
+ inline float value(Point<3,float> & p, Point<3,float> q,float pA,float pB)
+ {
+ float dist = norm(p-q);
+
+ return (pA - pB) * exp(dist * dist / var);
+ }
+};
+
+////////////////////////////////////////////////////////////////////////////
+
+BOOST_AUTO_TEST_SUITE( vector_dist_test )
+
+
+BOOST_AUTO_TEST_CASE( vector_dist_operators_test )
+{
+ if (create_vcluster().getProcessingUnits() > 3)
+ return;
+
+ Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
+
+ // Boundary conditions
+ size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+
+ // ghost
+ Ghost<3,float> ghost(0.05);
+
+ // vector type
+ typedef vector_dist<3,float,aggregate<float,float,float,VectorS<3,float>,VectorS<3,float>,VectorS<3,float>>> vtype;
+
+ vector_dist<3,float,aggregate<float,float,float,VectorS<3,float>,VectorS<3,float>,VectorS<3,float>>> vd(100,box,bc,ghost);
+
+ vd.getV<A>() = 1.0;
+ vd.getV<B>() = 2.0f;
+ vd.getV<C>() = 3.0;
+
+ check_values<A>(vd,1.0);
+ check_values<B>(vd,2.0);
+ check_values<C>(vd,3.0);
+
+ fill_values(vd);
+
+ vd.getV<A>() = vd.getV<B>() + 2.0 + vd.getV<B>() - 2.0*vd.getV<C>() / 5.0;
+ check_values_complex_expr(vd);
+
+ // Various combination of 2 operator
+
+ vd.getV<A>() = vd.getV<B>() + 2.0;
+ check_values_sum<float,vtype,A,B,C>(vd,2.0);
+ vd.getV<A>() = 2.0 + vd.getV<B>();
+ check_values_sum<float,vtype,A,B,C>(vd,2.0);
+ vd.getV<A>() = vd.getV<C>() + vd.getV<B>();
+ check_values_sum<float,vtype,A,B,C>(vd,vd);
+
+ vd.getV<A>() = vd.getV<B>() - 2.0;
+ check_values_sub<float,vtype,A,B,C>(vd,2.0);
+ vd.getV<A>() = 2.0 - vd.getV<B>();
+ check_values_sub<float,vtype,A,B,C>(2.0,vd);
+ vd.getV<A>() = vd.getV<C>() - vd.getV<B>();
+ check_values_sub<float,vtype,A,B,C>(vd,vd);
+
+ vd.getV<A>() = vd.getV<B>() * 2.0;
+ check_values_mul<float,vtype,A,B,C>(vd,2.0);
+ vd.getV<A>() = 2.0 * vd.getV<B>();
+ check_values_mul<float,vtype,A,B,C>(vd,2.0);
+ vd.getV<A>() = vd.getV<C>() * vd.getV<B>();
+ check_values_mul<float,vtype,A,B,C>(vd,vd);
+
+ vd.getV<A>() = vd.getV<B>() / 2.0;
+ check_values_div<float,vtype,A,B,C>(vd,2.0);
+ vd.getV<A>() = 2.0 / vd.getV<B>();
+ check_values_div<float,vtype,A,B,C>(2.0,vd);
+ vd.getV<A>() = vd.getV<C>() / vd.getV<B>();
+ check_values_div<float,vtype,A,B,C>(vd,vd);
+
+ // Variuos combination 3 operator
+
+ vd.getV<A>() = vd.getV<B>() + (vd.getV<C>() + vd.getV<B>());
+ check_values_sum_3<float,vtype,A,B,C>(vd);
+ vd.getV<A>() = (vd.getV<C>() + vd.getV<B>()) + vd.getV<B>();
+ check_values_sum_3<float,vtype,A,B,C>(vd);
+ vd.getV<A>() = (vd.getV<C>() + vd.getV<B>()) + (vd.getV<C>() + vd.getV<B>());
+ check_values_sum_4<float,vtype,A,B,C>(vd);
+
+ vd.getV<A>() = vd.getV<B>() - (vd.getV<C>() + vd.getV<B>());
+ check_values_sub_31<float,vtype,A,B,C>(vd);
+ vd.getV<A>() = (vd.getV<C>() + vd.getV<B>()) - vd.getV<B>();
+ check_values_sub_32<float,vtype,A,B,C>(vd);
+ vd.getV<A>() = (vd.getV<C>() + vd.getV<B>()) - (vd.getV<C>() + vd.getV<B>());
+ check_values_sub_4<float,vtype,A,B,C>(vd);
+
+ vd.getV<A>() = vd.getV<B>() * (vd.getV<C>() + vd.getV<B>());
+ check_values_mul_3<float,vtype,A,B,C>(vd);
+ vd.getV<A>() = (vd.getV<C>() + vd.getV<B>()) * vd.getV<B>();
+ check_values_mul_3<float,vtype,A,B,C>(vd);
+ vd.getV<A>() = (vd.getV<C>() + vd.getV<B>()) * (vd.getV<C>() + vd.getV<B>());
+ check_values_mul_4<float,vtype,A,B,C>(vd);
+
+ vd.getV<A>() = vd.getV<B>() / (vd.getV<C>() + vd.getV<B>());
+ check_values_div_31<float,vtype,A,B,C>(vd);
+ vd.getV<A>() = (vd.getV<C>() + vd.getV<B>()) / vd.getV<B>();
+ check_values_div_32<float,vtype,A,B,C>(vd);
+ vd.getV<A>() = (vd.getV<C>() + vd.getV<B>()) / (vd.getV<C>() + vd.getV<B>());
+ check_values_div_4<float,vtype,A,B,C>(vd);
+
+ // We try with vectors
+
+ // Various combination of 2 operator
+
+ vd.getV<PA>() = vd.getV<PB>() + 2.0;
+ check_values_sum<VectorS<3,float>,vtype,PA,PB,PC>(vd,2.0);
+ vd.getV<PA>() = 2.0 + vd.getV<PB>();
+ check_values_sum<VectorS<3,float>,vtype,PA,PB,PC>(vd,2.0);
+ vd.getV<PA>() = vd.getV<PC>() + vd.getV<PB>();
+ check_values_sum<VectorS<3,float>,vtype,PA,PB,PC>(vd,vd);
+
+ vd.getV<PA>() = vd.getV<PB>() - 2.0;
+ check_values_sub<VectorS<3,float>,vtype,PA,PB,PC>(vd,2.0);
+ vd.getV<PA>() = 2.0 - vd.getV<PB>();
+ check_values_sub<VectorS<3,float>,vtype,PA,PB,PC>(2.0,vd);
+ vd.getV<PA>() = vd.getV<PC>() - vd.getV<PB>();
+ check_values_sub<VectorS<3,float>,vtype,PA,PB,PC>(vd,vd);
+
+ vd.getV<PA>() = vd.getV<PB>() * 2.0;
+ check_values_mul<VectorS<3,float>,vtype,PA,PB,PC>(vd,2.0);
+ vd.getV<PA>() = 2.0 * vd.getV<PB>();
+ check_values_mul<VectorS<3,float>,vtype,PA,PB,PC>(vd,2.0);
+ vd.getV<PA>() = vd.getV<PC>() * vd.getV<PB>();
+ check_values_mul<VectorS<3,float>,vtype,PA,PB,PC>(vd,vd);
+
+ vd.getV<PA>() = vd.getV<PB>() / 2.0;
+ check_values_div<VectorS<3,float>,vtype,PA,PB,PC>(vd,2.0);
+ vd.getV<PA>() = 2.0 / vd.getV<PB>();
+ check_values_div<VectorS<3,float>,vtype,PA,PB,PC>(2.0,vd);
+ vd.getV<PA>() = vd.getV<PC>() / vd.getV<PB>();
+ check_values_div<VectorS<3,float>,vtype,PA,PB,PC>(vd,vd);
+
+ // Variuos combination 3 operator
+
+ vd.getV<PA>() = vd.getV<PB>() + (vd.getV<PC>() + vd.getV<PB>());
+ check_values_sum_3<VectorS<3,float>,vtype,PA,PB,PC>(vd);
+ vd.getV<PA>() = (vd.getV<PC>() + vd.getV<PB>()) + vd.getV<PB>();
+ check_values_sum_3<VectorS<3,float>,vtype,PA,PB,PC>(vd);
+ vd.getV<PA>() = (vd.getV<PC>() + vd.getV<PB>()) + (vd.getV<PC>() + vd.getV<PB>());
+ check_values_sum_4<VectorS<3,float>,vtype,PA,PB,PC>(vd);
+
+ vd.getV<PA>() = vd.getV<PB>() - (vd.getV<PC>() + vd.getV<PB>());
+ check_values_sub_31<VectorS<3,float>,vtype,PA,PB,PC>(vd);
+ vd.getV<PA>() = (vd.getV<PC>() + vd.getV<PB>()) - vd.getV<PB>();
+ check_values_sub_32<VectorS<3,float>,vtype,PA,PB,PC>(vd);
+ vd.getV<PA>() = (vd.getV<PC>() + vd.getV<PB>()) - (vd.getV<PC>() + vd.getV<PB>());
+ check_values_sub_4<VectorS<3,float>,vtype,PA,PB,PC>(vd);
+
+ vd.getV<PA>() = vd.getV<PB>() * (vd.getV<PC>() + vd.getV<PB>());
+ check_values_mul_3<VectorS<3,float>,vtype,PA,PB,PC>(vd);
+ vd.getV<PA>() = (vd.getV<PC>() + vd.getV<PB>()) * vd.getV<PB>();
+ check_values_mul_3<VectorS<3,float>,vtype,PA,PB,PC>(vd);
+ vd.getV<PA>() = (vd.getV<PC>() + vd.getV<PB>()) * (vd.getV<PC>() + vd.getV<PB>());
+ check_values_mul_4<VectorS<3,float>,vtype,PA,PB,PC>(vd);
+ vd.getV<A>() = vd.getV<PB>() * (vd.getV<PC>() + vd.getV<PB>());
+ check_values_mul_3<float,vtype,A,PB,PC>(vd);
+ vd.getV<A>() = (vd.getV<PC>() + vd.getV<PB>()) * vd.getV<PB>();
+ check_values_mul_3<float,vtype,A,PB,PC>(vd);
+ vd.getV<A>() = (vd.getV<PC>() + vd.getV<PB>()) * (vd.getV<PC>() + vd.getV<PB>());
+ check_values_mul_4<float,vtype,A,PB,PC>(vd);
+
+ vd.getV<PA>() = vd.getV<PB>() / (vd.getV<PC>() + vd.getV<PB>());
+ check_values_div_31<VectorS<3,float>,vtype,PA,PB,PC>(vd);
+ vd.getV<PA>() = (vd.getV<PC>() + vd.getV<PB>()) / vd.getV<PB>();
+ check_values_div_32<VectorS<3,float>,vtype,PA,PB,PC>(vd);
+ vd.getV<PA>() = (vd.getV<PC>() + vd.getV<PB>()) / (vd.getV<PC>() + vd.getV<PB>());
+ check_values_div_4<VectorS<3,float>,vtype,PA,PB,PC>(vd);
+
+ // normalization function
+
+ vd.getV<A>() = vd.getV<PB>() * vd.getV<PC>() + norm(vd.getV<PC>() + vd.getV<PB>()) + distance(vd.getV<PC>(),vd.getV<PB>());
+ check_values_scal_norm_dist(vd);
+
+ // we apply an exponential kernel to calculate something
+
+ auto cl = vd.getCellList(0.01);
+ exp_kernel ker(0.2);
+
+ vd.getV<A>() = applyKernel(vd.getV<PC>() * vd.getV<PB>() + norm(vd.getV<PB>()),vd,cl,ker) + vd.getV<C>();
+ check_values_apply_kernel(vd);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+
+
+#endif /* OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_UNIT_TESTS_HPP_ */
diff --git a/src/main.cpp b/src/main.cpp
index 5273f673c29eed21fd51e5c101c16daddff2e759..db22b925a5b46a32bd87ac7794db2334b4117e8b 100755
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -12,3 +12,4 @@
#include "FiniteDifference/eq_unit_test_3d.hpp"
#include "util/util_num_unit_tests.hpp"
#include "PSE/Kernels_unit_tests.hpp"
+#include "Operators/Vector/vector_dist_operators_unit_tests.hpp"