diff --git a/Makefile.am b/Makefile.am
index 80f456e1b82788f91ba4f8aa28a659770b094e97..437d908cb4fb476ceb849aa3d0370e4dfe673d55 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -1,3 +1,26 @@
 SUBDIRS = src images openfpm_data openfpm_io openfpm_devices openfpm_vcluster openfpm_numerics
 
-bin_PROGRAMS = 
+bin_PROGRAMS =
+
+test_pdata:
+	cd src && make test
+
+test_data:
+	cd openfpm_data/src && make test
+
+test_devices:
+	cd openfpm_devices/src && make test
+
+test_vcluster:
+	cd openfpm_vcluster/src && make test
+
+test_io:
+	cd openfpm_io/src && make test
+
+test_numerics:
+	cd openfpm_numerics/src && make test
+
+test: test_devices test_data test_vcluster test_pdata test_io test_numerics
+
+.PHONY: test_pdata test_data test_devices test_vcluster test_io test_numerics
+
diff --git a/src/Grid/grid_dist_id.hpp b/src/Grid/grid_dist_id.hpp
index ca7f207f7cb052498246c066b885cd582c8275d8..358a5e591ff706ee75abbed1ce3c726e0136bac0 100644
--- a/src/Grid/grid_dist_id.hpp
+++ b/src/Grid/grid_dist_id.hpp
@@ -549,7 +549,7 @@ protected:
 	 */
 	Point<dim,St> getOffset(size_t i)
 	{
-		return Point<dim,St>(gdb_ext.get(i).origin) * cd_sm.getCellBox().getP2();
+		return pmul(Point<dim,St>(gdb_ext.get(i).origin), cd_sm.getCellBox().getP2());
 	}
 
 	/*! \brief Given a local sub-domain i with a local grid Domain + ghost return the part of the local grid that is domain
diff --git a/src/Makefile.am b/src/Makefile.am
index bf900784043faf7da206ac978de346a8ccebc311..9ed5530bd1f0bbfc126721f4cf7c8a8b76d1d05d 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -23,3 +23,6 @@ libofpm_pdata_a_CFLAGS =
 
 .cu.o :
 	$(NVCC) $(NVCCFLAGS) -o $@ -c $<
+
+test: pdata
+	source $(HOME)/openfpm_vars &&  cd .. && mpirun -np 3 ./src/pdata && mpirun -np 4 ./src/pdata
diff --git a/src/Vector/vector_dist_operators.hpp b/src/Vector/vector_dist_operators.hpp
index da74fadac3c5cc6d0625b65a6a9c804248ffb0e3..660779247ed836eb95aee4a0f5bace85e6050af2 100644
--- a/src/Vector/vector_dist_operators.hpp
+++ b/src/Vector/vector_dist_operators.hpp
@@ -8,17 +8,68 @@
 #ifndef SRC_VECTOR_VECTOR_DIST_OPERATORS_HPP_
 #define SRC_VECTOR_VECTOR_DIST_OPERATORS_HPP_
 
+/*! \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)
+	{
+	}
+};
+
+
 #define SUM 1
 #define SUB 2
 #define MUL 3
 #define DIV 4
 
+/*! \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,SUM>
 {
@@ -27,18 +78,27 @@ class vector_dist_expression_op<exp1,exp2,SUM>
 
 public:
 
-	typedef decltype(o1.value(vect_dist_key_dx(0))) r_type;
-
-	vector_dist_expression_op(const exp1 & o1, const exp2 & o2)
+	inline vector_dist_expression_op(const exp1 & o1, const exp2 & o2)
 	:o1(o1),o2(o2)
 	{}
 
-	template<typename r_type=typename std::remove_reference<decltype(o1.value(vect_dist_key_dx(0)))>::type > inline r_type value(const vect_dist_key_dx & key) const
+	/*! \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,SUB>
 {
@@ -47,17 +107,28 @@ class vector_dist_expression_op<exp1,exp2,SUB>
 
 public:
 
-	vector_dist_expression_op(const exp1 & o1, const exp2 & o2)
+	inline vector_dist_expression_op(const exp1 & o1, const exp2 & o2)
 	:o1(o1),o2(o2)
 	{}
 
-	template<typename r_type=typename std::remove_reference<decltype(o1.value(vect_dist_key_dx(0)))>::type > inline r_type value(const vect_dist_key_dx & key) const
+	/*! \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,MUL>
 {
@@ -70,13 +141,24 @@ public:
 	:o1(o1),o2(o2)
 	{}
 
-	template<typename r_type=typename std::remove_reference<decltype(o1.value(vect_dist_key_dx(0)))>::type > inline r_type value(const vect_dist_key_dx & key) const
+	/*! \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,DIV>
 {
@@ -89,17 +171,22 @@ public:
 	:o1(o1),o2(o2)
 	{}
 
-	template<typename r_type=typename std::remove_reference<decltype(o1.value(vect_dist_key_dx(0)))>::type > inline r_type value(const vect_dist_key_dx & key) const
+	/*! \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 Main class that encapsulate a vector operation expression
+/*! \brief Main class that encapsulate a vector properties
  *
- * \param prp property involved
- * \param exp expression type
+ * \tparam prp property involved
+ * \tparam vector involved
  *
  */
 template<unsigned int prp, typename vector>
@@ -113,12 +200,17 @@ public:
 	: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 evaluate the expression and store it in the vector
+	/*! \brief Fill the vector property with the evaluated expression
 	 *
 	 * \param v_exp expression to evaluate
 	 *
@@ -131,6 +223,8 @@ public:
 		{
 			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;
@@ -139,9 +233,9 @@ public:
 		return v;
 	}
 
-	/*! \brief evaluate the expression and store it in the vector
+	/*! \brief Fill the vector property with the double
 	 *
-	 * \param v_exp expression to evaluate
+	 * \param d value to fill
 	 *
 	 */
 	vector & operator=(double d)
@@ -161,10 +255,9 @@ public:
 	}
 };
 
-/*! \brief Main class that encapsulate a vector operation expression
+/*! \brief Main class that encapsulate a double constant
  *
- * \param prp property involved
- * \param exp expression type
+ * \param prp no meaning
  *
  */
 template<unsigned int prp>
@@ -178,6 +271,11 @@ public:
 	: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;
@@ -185,12 +283,12 @@ public:
 };
 
 
-/* \brief sum two distributed vectors
+/* \brief sum two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<unsigned int p1, unsigned int p2, typename v1, typename v2>
@@ -202,12 +300,12 @@ operator+(const vector_dist_expression<p1,v1> & va, const vector_dist_expression
 	return exp_sum;
 }
 
-/* \brief sum two distributed vectors
+/* \brief sum two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<typename exp1 , typename exp2, unsigned int op1, unsigned int prp1, typename v1>
@@ -219,12 +317,12 @@ operator+(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist
 	return exp_sum;
 }
 
-/* \brief sum two distributed vectors
+/* \brief sum two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<typename exp1 , typename exp2, unsigned int op1, unsigned int prp1, typename v1>
@@ -236,12 +334,12 @@ operator+(const vector_dist_expression<prp1,v1> & va, const vector_dist_expressi
 	return exp_sum;
 }
 
-/* \brief sum two distributed vectors
+/* \brief sum two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<typename exp1 , typename exp2, unsigned int op1, typename exp3 , typename exp4, unsigned int op2>
@@ -253,12 +351,12 @@ operator+(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist
 	return exp_sum;
 }
 
-/* \brief sum two distributed vectors
+/* \brief sum two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<unsigned int prp1 , typename v1>
@@ -270,12 +368,12 @@ operator+(const vector_dist_expression<prp1,v1> & va, double d)
 	return exp_sum;
 }
 
-/* \brief sum two distributed vectors
+/* \brief sum two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<unsigned int prp1 , typename v1>
@@ -287,12 +385,12 @@ operator+(double d, const vector_dist_expression<prp1,v1> & vb)
 	return exp_sum;
 }
 
-/* \brief sum two distributed vectors
+/* \brief sum two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<typename exp1 , typename exp2, unsigned int op1>
@@ -305,12 +403,12 @@ operator+(const vector_dist_expression_op<exp1,exp2,op1> & va, double d)
 }
 
 
-/* \brief sum two distributed vectors
+/* \brief subtract two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<unsigned int p1, unsigned int p2, typename v1, typename v2>
@@ -323,12 +421,12 @@ operator-(const vector_dist_expression<p1,v1> & va, const vector_dist_expression
 }
 
 
-/* \brief sum two distributed vectors
+/* \brief subtract two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<typename exp1, typename exp2, unsigned int op1, unsigned int p2, typename v2>
@@ -340,12 +438,12 @@ operator-(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist
 	return exp_sum;
 }
 
-/* \brief sum two distributed vectors
+/* \brief subtract two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<typename exp1, typename exp2, unsigned int op1, unsigned int p2, typename v2>
@@ -357,12 +455,12 @@ operator-(const vector_dist_expression<p2,v2> & va, const vector_dist_expression
 	return exp_sum;
 }
 
-/* \brief sum two distributed vectors
+/* \brief subtract two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<typename exp1, typename exp2, unsigned int op1, typename exp3, typename exp4, unsigned int op2>
@@ -374,12 +472,12 @@ operator-(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist
 	return exp_sum;
 }
 
-/* \brief sum two distributed vectors
+/* \brief subtract two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<unsigned int prp1, typename v1>
@@ -391,12 +489,12 @@ operator-(const vector_dist_expression<prp1,v1> & va, double d)
 	return exp_sum;
 }
 
-/* \brief sum two distributed vectors
+/* \brief subtract two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<unsigned int prp1, typename v1>
@@ -408,12 +506,12 @@ operator-(double d, const vector_dist_expression<prp1,v1> & vb)
 	return exp_sum;
 }
 
-/* \brief multiplication double with a vector expression
+/* \brief Multiply two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<unsigned int p2, typename v2>
@@ -425,12 +523,12 @@ operator*(double d, const vector_dist_expression<p2,v2> & vb)
 	return exp_sum;
 }
 
-/* \brief multiplication double with a vector expression
+/* \brief Multiply two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<unsigned int p2, typename v2>
@@ -442,12 +540,12 @@ operator*(const vector_dist_expression<p2,v2> & va, double d)
 	return exp_sum;
 }
 
-/* \brief multiplication double with a vector expression
+/* \brief Multiply two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<unsigned int p1, typename v1,unsigned int p2, typename v2>
@@ -459,12 +557,12 @@ operator*(const vector_dist_expression<p1,v1> & va, const vector_dist_expression
 	return exp_sum;
 }
 
-/* \brief multiplication double with a vector expression
+/* \brief Multiply two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<unsigned int p1, typename v1, typename exp1, typename exp2, unsigned int op1>
@@ -476,12 +574,12 @@ operator*(const vector_dist_expression<p1,v1> & va, const vector_dist_expression
 	return exp_sum;
 }
 
-/* \brief multiplication double with a vector expression
+/* \brief Multiply two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<unsigned int p1, typename v1, typename exp1, typename exp2, unsigned int op1>
@@ -493,12 +591,12 @@ operator*(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist
 	return exp_sum;
 }
 
-/* \brief multiplication double with a vector expression
+/* \brief Multiply two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<typename exp1, typename exp2, unsigned int op1, typename exp3 , typename exp4, unsigned int op2>
@@ -510,12 +608,12 @@ operator*(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist
 	return exp_sum;
 }
 
-/* \brief multiplication double with a vector expression
+/* \brief Divide two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<typename exp1, typename exp2, unsigned int op1>
@@ -528,12 +626,12 @@ operator/(const vector_dist_expression_op<exp1,exp2,op1> & va, double d)
 }
 
 
-/* \brief multiplication double with a vector expression
+/* \brief Divide two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<typename exp1, typename exp2, unsigned int op1>
@@ -545,12 +643,12 @@ operator/(double d, const vector_dist_expression_op<exp1,exp2,op1> & va)
 	return exp_sum;
 }
 
-/* \brief multiplication double with a vector expression
+/* \brief Divide two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<unsigned int prp1, typename v1>
@@ -562,29 +660,29 @@ operator/(const vector_dist_expression<prp1,v1> & va, double d)
 	return exp_sum;
 }
 
-/* \brief multiplication double with a vector expression
+/* \brief Divide two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \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>,DIV>
+inline vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<prp1,v1>,DIV>
 operator/(double d, const vector_dist_expression<prp1,v1> & va)
 {
-	vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,DIV> exp_sum(va,vector_dist_expression<0,double>(d));
+	vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<prp1,v1>,DIV> exp_sum(vector_dist_expression<0,double>(d),va);
 
 	return exp_sum;
 }
 
-/* \brief multiplication double with a vector expression
+/* \brief Divide two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<unsigned int prp1, typename v1, unsigned int prp2, typename v2>
@@ -596,12 +694,12 @@ operator/(const vector_dist_expression<prp1,v1> & va, const vector_dist_expressi
 	return exp_sum;
 }
 
-/* \brief multiplication double with a vector expression
+/* \brief Divide two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<unsigned int prp1, typename v1, typename exp1,typename exp2, unsigned int op1>
@@ -613,12 +711,12 @@ operator/(const vector_dist_expression<prp1,v1> & va, const vector_dist_expressi
 	return exp_sum;
 }
 
-/* \brief multiplication double with a vector expression
+/* \brief Divide two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<unsigned int prp1, typename v1, typename exp1,typename exp2, unsigned int op1>
@@ -630,12 +728,12 @@ operator/(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist
 	return exp_sum;
 }
 
-/* \brief multiplication double with a vector expression
+/* \brief Divide two distributed vector expression
  *
- * \param v1 vector1
- * \param v2 vector2
+ * \param va vector expression one
+ * \param vb vector expression two
  *
- * \return a grid_key_dx_expression that encapsulate the expression
+ * \return an object that encapsulate the expression
  *
  */
 template<typename exp1,typename exp2, unsigned int op1, typename exp3, typename exp4, unsigned int op2>
diff --git a/src/Vector/vector_dist_operators_unit_tests.hpp b/src/Vector/vector_dist_operators_unit_tests.hpp
index fbbed13fb76f80d862c30777e7a735452132a560..856f6ed41cd6f3f5a578a5d5f88b48019663261e 100644
--- a/src/Vector/vector_dist_operators_unit_tests.hpp
+++ b/src/Vector/vector_dist_operators_unit_tests.hpp
@@ -57,7 +57,7 @@ template <typename vector> bool check_values_complex_expr(const vector & vd)
 	return ret;
 }
 
-template <typename vector> bool check_values_sum(const vector & vd, double d)
+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();
@@ -66,8 +66,8 @@ template <typename vector> bool check_values_sum(const vector & vd, double d)
 	{
 		auto key = it.get();
 
-		float base1 = vd.template getProp<B>(key) + d;
-		float base2 = vd.template getProp<A>(key);
+		rtype base1 = vd.template getProp<B>(key) + d;
+		rtype base2 = vd.template getProp<A>(key);
 
 		ret &=  base1 == base2;
 
@@ -79,7 +79,7 @@ template <typename vector> bool check_values_sum(const vector & vd, double d)
 	return ret;
 }
 
-template <typename vector> bool check_values_sum(const vector & vd1, const vector & vd2)
+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();
@@ -88,8 +88,8 @@ template <typename vector> bool check_values_sum(const vector & vd1, const vecto
 	{
 		auto key = it.get();
 
-		float base1 = vd1.template getProp<B>(key) + vd2.template getProp<C>(key);
-		float base2 = vd1.template getProp<A>(key);
+		rtype base1 = vd1.template getProp<B>(key) + vd2.template getProp<C>(key);
+		rtype base2 = vd1.template getProp<A>(key);
 
 		ret &=  base1 == base2;
 
@@ -101,7 +101,7 @@ template <typename vector> bool check_values_sum(const vector & vd1, const vecto
 	return ret;
 }
 
-template <typename vector> bool check_values_sum_3(const vector & vd1)
+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();
@@ -110,8 +110,8 @@ template <typename vector> bool check_values_sum_3(const vector & vd1)
 	{
 		auto key = it.get();
 
-		float base1 = vd1.template getProp<B>(key) + vd1.template getProp<C>(key) + vd1.template getProp<B>(key);
-		float base2 = vd1.template getProp<A>(key);
+		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;
 
@@ -123,7 +123,7 @@ template <typename vector> bool check_values_sum_3(const vector & vd1)
 	return ret;
 }
 
-template <typename vector> bool check_values_sum_4(const vector & vd1)
+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();
@@ -132,8 +132,8 @@ template <typename vector> bool check_values_sum_4(const vector & vd1)
 	{
 		auto key = it.get();
 
-		float base1 = vd1.template getProp<B>(key) + vd1.template getProp<C>(key) + vd1.template getProp<B>(key) + vd1.template getProp<C>(key);
-		float base2 = vd1.template getProp<A>(key);
+		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;
 
@@ -145,7 +145,7 @@ template <typename vector> bool check_values_sum_4(const vector & vd1)
 	return ret;
 }
 
-template <typename vector> bool check_values_sub(const vector & vd, double d)
+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();
@@ -154,8 +154,8 @@ template <typename vector> bool check_values_sub(const vector & vd, double d)
 	{
 		auto key = it.get();
 
-		float base1 = vd.template getProp<B>(key) - d;
-		float base2 = vd.template getProp<A>(key);
+		rtype base1 = vd.template getProp<B>(key) - d;
+		rtype base2 = vd.template getProp<A>(key);
 
 		ret &=  base1 == base2;
 
@@ -167,7 +167,7 @@ template <typename vector> bool check_values_sub(const vector & vd, double d)
 	return ret;
 }
 
-template <typename vector> bool check_values_sub(double d, const vector & vd)
+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();
@@ -176,8 +176,30 @@ template <typename vector> bool check_values_sub(double d, const vector & vd)
 	{
 		auto key = it.get();
 
-		float base1 = d - vd.template getProp<B>(key);
-		float base2 = vd.template getProp<A>(key);
+		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;
 
@@ -189,7 +211,7 @@ template <typename vector> bool check_values_sub(double d, const vector & vd)
 	return ret;
 }
 
-template <typename vector> bool check_values_sub(const vector & vd1, const vector & vd2)
+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();
@@ -198,8 +220,8 @@ template <typename vector> bool check_values_sub(const vector & vd1, const vecto
 	{
 		auto key = it.get();
 
-		float base1 = vd1.template getProp<C>(key) - vd2.template getProp<B>(key);
-		float base2 = vd1.template getProp<A>(key);
+		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;
 
@@ -211,7 +233,53 @@ template <typename vector> bool check_values_sub(const vector & vd1, const vecto
 	return ret;
 }
 
-template <typename vector> bool check_values_mul(const vector & vd, double d)
+
+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();
@@ -220,8 +288,53 @@ template <typename vector> bool check_values_mul(const vector & vd, double d)
 	{
 		auto key = it.get();
 
-		float base1 = vd.template getProp<B>(key) * d;
-		float base2 = vd.template getProp<A>(key);
+		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;
 
@@ -234,7 +347,7 @@ template <typename vector> bool check_values_mul(const vector & vd, double d)
 }
 
 
-template <typename vector> bool check_values_mul(const vector & vd1, const vector & vd2)
+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();
@@ -243,8 +356,8 @@ template <typename vector> bool check_values_mul(const vector & vd1, const vecto
 	{
 		auto key = it.get();
 
-		float base1 = vd1.template getProp<C>(key) * vd2.template getProp<B>(key);
-		float base2 = vd1.template getProp<A>(key);
+		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;
 
@@ -257,7 +370,8 @@ template <typename vector> bool check_values_mul(const vector & vd1, const vecto
 }
 
 
-template <typename vector> bool check_values_div(const vector & vd, double d)
+
+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();
@@ -266,8 +380,8 @@ template <typename vector> bool check_values_div(const vector & vd, double d)
 	{
 		auto key = it.get();
 
-		float base1 = vd.template getProp<B>(key) / d;
-		float base2 = vd.template getProp<A>(key);
+		rtype base1 = vd.template getProp<B>(key) / d;
+		rtype base2 = vd.template getProp<A>(key);
 
 		ret &=  base1 == base2;
 
@@ -279,7 +393,7 @@ template <typename vector> bool check_values_div(const vector & vd, double d)
 	return ret;
 }
 
-template <typename vector> bool check_values_div(double d, const vector & vd)
+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();
@@ -288,8 +402,30 @@ template <typename vector> bool check_values_div(double d, const vector & vd)
 	{
 		auto key = it.get();
 
-		float base1 = d / vd.template getProp<B>(key);
-		float base2 = vd.template getProp<A>(key);
+		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;
 
@@ -301,7 +437,7 @@ template <typename vector> bool check_values_div(double d, const vector & vd)
 	return ret;
 }
 
-template <typename vector> bool check_values_div(const vector & vd1, const vector & vd2)
+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();
@@ -310,8 +446,8 @@ template <typename vector> bool check_values_div(const vector & vd1, const vecto
 	{
 		auto key = it.get();
 
-		float base1 = vd1.template getProp<C>(key) / vd2.template getProp<B>(key);
-		float base2 = vd1.template getProp<A>(key);
+		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;
 
@@ -323,6 +459,52 @@ template <typename vector> bool check_values_div(const vector & vd1, const vecto
 	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> void fill_values(const vector & v)
 {
 	auto it = v.getDomainIterator();
@@ -335,6 +517,13 @@ template <typename vector> void fill_values(const vector & v)
 		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) = p.getKey()+1+k;
+			v.template getProp<PB>(p) = 2.0*p.getKey()+1+k;
+			v.template getProp<PC>(p) = 3.0*p.getKey()+1+k;
+		}
+
 		++it;
 	}
 }
@@ -362,6 +551,9 @@ BOOST_AUTO_TEST_CASE( vector_dist_operators_test )
 	// 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;
@@ -380,53 +572,133 @@ BOOST_AUTO_TEST_CASE( vector_dist_operators_test )
 	// Various combination of 2 operator
 
 	vd.getV<A>() = vd.getV<B>() + 2.0;
-	check_values_sum(vd,2.0);
+	check_values_sum<float,vtype,A,B,C>(vd,2.0);
 	vd.getV<A>() = 2.0 + vd.getV<B>();
-	check_values_sum(vd,2.0);
+	check_values_sum<float,vtype,A,B,C>(vd,2.0);
 	vd.getV<A>() = vd.getV<C>() + vd.getV<B>();
-	check_values_sum(vd,vd);
+	check_values_sum<float,vtype,A,B,C>(vd,vd);
 
 	vd.getV<A>() = vd.getV<B>() - 2.0;
-	check_values_sub(vd,2.0);
+	check_values_sub<float,vtype,A,B,C>(vd,2.0);
 	vd.getV<A>() = 2.0 - vd.getV<B>();
-	check_values_sub(2.0,vd);
+	check_values_sub<float,vtype,A,B,C>(2.0,vd);
 	vd.getV<A>() = vd.getV<C>() - vd.getV<B>();
-	check_values_sub(vd,vd);
+	check_values_sub<float,vtype,A,B,C>(vd,vd);
 
 	vd.getV<A>() = vd.getV<B>() * 2.0;
-	check_values_mul(vd,2.0);
+	check_values_mul<float,vtype,A,B,C>(vd,2.0);
 	vd.getV<A>() = 2.0 * vd.getV<B>();
-	check_values_mul(vd,2.0);
+	check_values_mul<float,vtype,A,B,C>(vd,2.0);
 	vd.getV<A>() = vd.getV<C>() * vd.getV<B>();
-	check_values_mul(vd,vd);
+	check_values_mul<float,vtype,A,B,C>(vd,vd);
 
 	vd.getV<A>() = vd.getV<B>() / 2.0;
-	check_values_div(vd,2.0);
+	check_values_div<float,vtype,A,B,C>(vd,2.0);
 	vd.getV<A>() = 2.0 / vd.getV<B>();
-	check_values_div(vd,2.0);
+	check_values_div<float,vtype,A,B,C>(2.0,vd);
 	vd.getV<A>() = vd.getV<C>() / vd.getV<B>();
-	check_values_div(vd,vd);
+	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(vd);
+	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(vd);
+	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(vd);
+	check_values_sum_4<float,vtype,A,B,C>(vd);
 
-/*	vd.getV<A>() = vd.getV<B>() - (vd.getV<C>() + vd.getV<B>());
+	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>();
-	vd.getV<A>() = (vd.getV<C>() + vd.getV<B>()) / (vd.getV<C>() + 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
+
 
 }
 
diff --git a/src/Vector/vector_dist_unit_test.hpp b/src/Vector/vector_dist_unit_test.hpp
index a570c1b0500b410528b05e797f153547c0eabe22..7b09bbac89dd93425b26f619dbe7aa1665bb5a3e 100644
--- a/src/Vector/vector_dist_unit_test.hpp
+++ b/src/Vector/vector_dist_unit_test.hpp
@@ -140,7 +140,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost )
 	spacing = spacing / g_div;
 
 	// middle spacing
-	Point<2,float> m_spacing = spacing / 2;
+	Point<2,float> m_spacing = spacing / 2.0;
 
 	// set the ghost based on the radius cut off (make just a little bit smaller than the spacing)
 	Ghost<2,float> g(spacing.get(0) - spacing .get(0) * 0.0001);