Adding missing files

src/DMatrix/EMatrix.hpp

+/*
+ * EMatrix.hpp
+ *
+ *  Created on: Feb 12, 2018
+ *      Author: i-bird
+ */
+
+#ifndef OPENFPM_DATA_SRC_SPACE_EMATRIX_HPP_
+#define OPENFPM_DATA_SRC_SPACE_EMATRIX_HPP_
+
+#include "config.h"
+
+#ifdef HAVE_EIGEN
+#include <Eigen/Dense>
+#include "memory/ExtPreAlloc.hpp"
+#include "memory/HeapMemory.hpp"
+#include "Packer_Unpacker/Packer_util.hpp"
+#include "Packer_Unpacker/Packer.hpp"
+#include "Packer_Unpacker/Unpacker.hpp"
+
+/*! \brief it is just a wrapper for Eigen matrix compatible for openfpm serialization
+ *
+ *
+ */
+template<typename _Scalar, int _Rows, int _Cols>
+class EMatrix : public Eigen::Matrix<_Scalar,_Rows,_Cols>
+{
+public:
+
+	/*! \brief It calculate the number of byte required to serialize the object
+	 *
+	 * \tparam prp list of properties
+	 *
+	 * \param req reference to the total counter required to pack the information
+	 *
+	 */
+	template<int ... prp> inline void packRequest(size_t & req) const
+	{
+		// Memory required to serialize the Matrix
+		req += sizeof(_Scalar)*this->rows()*this->cols() + 2*sizeof(_Scalar);
+	}
+
+
+	/*! \brief pack a vector selecting the properties to pack
+	 *
+	 * \param mem preallocated memory where to pack the vector
+	 * \param sts pack-stat info
+	 *
+	 */
+	template<int ... prp> inline void pack(ExtPreAlloc<HeapMemory> & mem, Pack_stat & sts) const
+	{
+		//Pack the number of rows and colums
+		Packer<size_t, HeapMemory>::pack(mem,this->rows(),sts);
+		Packer<size_t, HeapMemory>::pack(mem,this->cols(),sts);
+
+		mem.allocate(sizeof(_Scalar)*this->rows()*this->cols());
+		void * dst = mem.getPointer();
+		void * src = (void *)this->data();
+
+		memcpy(dst,src,sizeof(_Scalar)*this->rows()*this->cols());
+		sts.incReq();
+	}
+
+	/*! \brief unpack a vector
+	 *
+	 * \param mem preallocated memory from where to unpack the vector
+	 * \param ps unpack-stat info
+	 */
+	template<int ... prp> inline void unpack(ExtPreAlloc<HeapMemory> & mem, Unpack_stat & ps)
+	{
+		size_t row;
+		size_t col;
+
+		//Pack the number of rows and colums
+		Unpacker<size_t, HeapMemory>::unpack(mem,row,ps);
+		Unpacker<size_t, HeapMemory>::unpack(mem,col,ps);
+
+		this->resize(row,col);
+
+		void * dst = this->data();
+		void * src = mem.getPointerOffset(ps.getOffset());
+
+		memcpy(dst,src,sizeof(_Scalar)*row*col);
+
+		ps.addOffset(sizeof(_Scalar)*row*col);
+	}
+
+};
+
+// Standard typedef from eigen
+typedef EMatrix< std::complex< double >, 2, 2 > 	EMatrix2cd;
+typedef EMatrix< std::complex< float >, 2, 2 > 	EMatrix2cf;
+typedef EMatrix< double, 2, 2 > 	EMatrix2d;
+typedef EMatrix< float, 2, 2 > 	EMatrix2f;
+typedef EMatrix< int, 2, 2 > 	EMatrix2i;
+typedef EMatrix< std::complex< double >, 2, Eigen::Dynamic > 	EMatrix2Xcd;
+typedef EMatrix< std::complex< float >, 2, Eigen::Dynamic > 	EMatrix2Xcf;
+typedef EMatrix< double, 2, Eigen::Dynamic > 	EMatrix2Xd;
+typedef EMatrix< float, 2, Eigen::Dynamic > 	EMatrix2Xf;
+typedef EMatrix< int, 2, Eigen::Dynamic > 	EMatrix2Xi;
+typedef EMatrix< std::complex< double >, 3, 3 > 	EMatrix3cd;
+typedef EMatrix< std::complex< float >, 3, 3 > 	EMatrix3cf;
+typedef EMatrix< double, 3, 3 > 	EMatrix3d;
+typedef EMatrix< float, 3, 3 > 	EMatrix3f;
+typedef EMatrix< int, 3, 3 > 	EMatrix3i;
+typedef EMatrix< std::complex< double >, 3, Eigen::Dynamic > 	EMatrix3Xcd;
+typedef EMatrix< std::complex< float >, 3, Eigen::Dynamic > 	EMatrix3Xcf;
+typedef EMatrix< double, 3, Eigen::Dynamic > 	EMatrix3Xd;
+typedef EMatrix< float, 3, Eigen::Dynamic > 	EMatrix3Xf;
+typedef EMatrix< int, 3, Eigen::Dynamic > 	EMatrix3Xi;
+typedef EMatrix< std::complex< double >, 4, 4 > 	EMatrix4cd;
+typedef EMatrix< std::complex< float >, 4, 4 > 	EMatrix4cf;
+typedef EMatrix< double, 4, 4 > 	EMatrix4d;
+typedef EMatrix< float, 4, 4 > 	EMatrix4f;
+typedef EMatrix< int, 4, 4 > 	EMatrix4i;
+typedef EMatrix< std::complex< double >, 4, Eigen::Dynamic > 	EMatrix4Xcd;
+typedef EMatrix< std::complex< float >, 4, Eigen::Dynamic > 	EMatrix4Xcf;
+typedef EMatrix< double, 4, Eigen::Dynamic > 	EMatrix4Xd;
+typedef EMatrix< float, 4, Eigen::Dynamic > 	EMatrix4Xf;
+typedef EMatrix< int, 4, Eigen::Dynamic > 	EMatrix4Xi;
+typedef EMatrix< std::complex< double >, Eigen::Dynamic, 2 > 	EMatrixX2cd;
+typedef EMatrix< std::complex< float >, Eigen::Dynamic, 2 > 	EMatrixX2cf;
+typedef EMatrix< double, Eigen::Dynamic, 2 > 	EMatrixX2d;
+typedef EMatrix< float, Eigen::Dynamic, 2 > 	EMatrixX2f;
+typedef EMatrix< int, Eigen::Dynamic, 2 > 	EMatrixX2i;
+typedef EMatrix< std::complex< double >, Eigen::Dynamic, 3 > 	EMatrixX3cd;
+typedef EMatrix< std::complex< float >, Eigen::Dynamic, 3 > 	EMatrixX3cf;
+typedef EMatrix< double, Eigen::Dynamic, 3 > 	EMatrixX3d;
+typedef EMatrix< float, Eigen::Dynamic, 3 > 	EMatrixX3f;
+typedef EMatrix< int, Eigen::Dynamic, 3 > 	EMatrixX3i;
+typedef EMatrix< std::complex< double >, Eigen::Dynamic, 4 > 	EMatrixX4cd;
+typedef EMatrix< std::complex< float >, Eigen::Dynamic, 4 > 	EMatrixX4cf;
+typedef EMatrix< double, Eigen::Dynamic, 4 > 	EMatrixX4d;
+typedef EMatrix< float, Eigen::Dynamic, 4 > 	EMatrixX4f;
+typedef EMatrix< int, Eigen::Dynamic, 4 > 	EMatrixX4i;
+typedef EMatrix< std::complex< double >, Eigen::Dynamic, Eigen::Dynamic > 	EMatrixXcd;
+typedef EMatrix< std::complex< float >, Eigen::Dynamic, Eigen::Dynamic > 	EMatrixXcf;
+typedef EMatrix< double, Eigen::Dynamic, Eigen::Dynamic > 	EMatrixXd;
+typedef EMatrix< float, Eigen::Dynamic, Eigen::Dynamic > 	EMatrixXf;
+typedef EMatrix< int, Eigen::Dynamic, Eigen::Dynamic > 	EMatrixXi;
+typedef EMatrix< std::complex< double >, 1, 2 > 	ERowVector2cd;
+typedef EMatrix< std::complex< float >, 1, 2 > 	ERowVector2cf;
+typedef EMatrix< double, 1, 2 > 	ERowVector2d;
+typedef EMatrix< float, 1, 2 > 	ERowVector2f;
+typedef EMatrix< int, 1, 2 > 	ERowVector2i;
+typedef EMatrix< std::complex< double >, 1, 3 > 	ERowVector3cd;
+typedef EMatrix< std::complex< float >, 1, 3 > 	ERowVector3cf;
+typedef EMatrix< double, 1, 3 > 	ERowVector3d;
+typedef EMatrix< float, 1, 3 > 	ERowVector3f;
+typedef EMatrix< int, 1, 3 > 	ERowVector3i;
+typedef EMatrix< std::complex< double >, 1, 4 > 	ERowVector4cd;
+typedef EMatrix< std::complex< float >, 1, 4 > 	ERowVector4cf;
+typedef EMatrix< double, 1, 4 > 	ERowVector4d;
+typedef EMatrix< float, 1, 4 > 	ERowVector4f;
+typedef EMatrix< int, 1, 4 > 	ERowVector4i;
+typedef EMatrix< std::complex< double >, 1, Eigen::Dynamic > 	ERowVectorXcd;
+typedef EMatrix< std::complex< float >, 1, Eigen::Dynamic > 	ERowVectorXcf;
+typedef EMatrix< double, 1, Eigen::Dynamic > 	ERowVectorXd;
+typedef EMatrix< float, 1, Eigen::Dynamic > 	ERowVectorXf;
+typedef EMatrix< int, 1, Eigen::Dynamic > 	ERowVectorXi;
+typedef EMatrix< std::complex< double >, 2, 1 > 	EVector2cd;
+typedef EMatrix< std::complex< float >, 2, 1 > 	EVector2cf;
+typedef EMatrix< double, 2, 1 > 	EVector2d;
+typedef EMatrix< float, 2, 1 > 	EVector2f;
+typedef EMatrix< int, 2, 1 > 	EVector2i;
+typedef EMatrix< std::complex< double >, 3, 1 > 	EVector3cd;
+typedef EMatrix< std::complex< float >, 3, 1 > 	EVector3cf;
+typedef EMatrix< double, 3, 1 > 	EVector3d;
+typedef EMatrix< float, 3, 1 > 	EVector3f;
+typedef EMatrix< int, 3, 1 > 	EVector3i;
+typedef EMatrix< std::complex< double >, 4, 1 > 	EVector4cd;
+typedef EMatrix< std::complex< float >, 4, 1 > 	EVector4cf;
+typedef EMatrix< double, 4, 1 > 	EVector4d;
+typedef EMatrix< float, 4, 1 > 	EVector4f;
+typedef EMatrix< int, 4, 1 > 	EVector4i;
+typedef EMatrix< std::complex< double >, Eigen::Dynamic, 1 > 	EVectorXcd;
+typedef EMatrix< std::complex< float >, Eigen::Dynamic, 1 > 	EVectorXcf;
+typedef EMatrix< double, Eigen::Dynamic, 1 > 	EVectorXd;
+typedef EMatrix< float, Eigen::Dynamic, 1 > 	EVectorXf;
+typedef EMatrix< int, Eigen::Dynamic, 1 > 	EVectorXi;
+
+#else
+
+// There is not EMatrix if we do not have Eigen
+
+#endif
+
+#endif /* OPENFPM_DATA_SRC_SPACE_EMATRIX_HPP_ */

src/Operators/Vector/vector_dist_operators_apply_kernel.hpp

@@ -8,6 +8,16 @@
 #ifndef OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_APPLY_KERNEL_HPP_
 #define OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_APPLY_KERNEL_HPP_
 
+//////// SET of small macro to make to define integral easy
+
+#define DEFINE_INTERACTION_3D(name) struct name \
+{\
+\
+	Point<3,double> value(const Point<3,double> & xp, const Point<3,double> xq)\
+	{
+
+#define END_INTERACTION }\
+						};
 
 /*! \brief is_expression check if a type is simple a type or is just an encapsulation of an expression
  *

src/interpolation/interpolation.hpp

@@ -11,6 +11,7 @@
 #include "NN/Mem_type/MemFast.hpp"
 #include "NN/CellList/CellList.hpp"
 #include "Grid/grid_dist_key.hpp"
+#include "Vector/vector_dist_key.hpp"
 
 #define INTERPOLATION_ERROR_OBJECT std::runtime_error("Runtime interpolation error");
 
@@ -407,8 +408,8 @@ struct inte_calc_impl
 	 *        kernel stencil for each local grid (sub-domain)
 	 *
 	 */
-	template<unsigned int prp_g, unsigned int prp_v, unsigned int m2p_or_p2m, unsigned int np_a_int, typename iterator, typename grid>
-																	 static inline void inte_calc(iterator & it,
+	template<unsigned int prp_g, unsigned int prp_v, unsigned int m2p_or_p2m, unsigned int np_a_int, typename grid>
+																	 static inline void inte_calc(const vect_dist_key_dx & key_p,
 																	 vector & vd,
 																	 const Box<vector::dims,typename vector::stype> & domain,
 																	 int (& ip)[vector::dims][kernel::np],
@@ -422,8 +423,6 @@ struct inte_calc_impl
 																	 const CellList<vector::dims,typename vector::stype,Mem_fast<>,shift<vector::dims,typename vector::stype>> & geo_cell,
 																	 openfpm::vector<agg_arr<openfpm::math::pow(kernel::np,vector::dims)>> & offsets)
 	{
-		auto key_p = it.get();
-
 		Point<vector::dims,typename vector::stype> p = vd.getPos(key_p);
 
 		// On which sub-domain we interpolate the particle
@@ -541,6 +540,18 @@ class interpolate
 	//! Type of the calculations
 	typedef typename vector::stype arr_type;
 
+	//! the offset for each sub-sub-domain
+	openfpm::vector<agg_arr<openfpm::math::pow(kernel::np,vector::dims)>> offsets;
+
+	//! kernel size
+	size_t sz[vector::dims];
+
+	//! grid spacing
+	typename vector::stype dx[vector::dims];
+
+	//! Simulation domain
+	Box<vector::dims,typename vector::stype> domain;
+
 	/*! \brief It calculate the interpolation stencil offsets
 	 *
 	 * \param offsets array where to store the linearized offset of the
@@ -624,6 +635,18 @@ public:
 				++g_sub;
 			}
 		}
+
+		for (size_t i = 0 ; i < vector::dims ; i++)
+		{sz[i] = kernel::np;}
+
+		calculate_the_offsets(offsets,sz);
+
+		// calculate spacing
+		for (size_t i = 0 ; i < vector::dims ; i++)
+		{dx[i] = 1.0/gd.spacing(i);}
+
+		// simulation domain
+		domain = vd.getDecomposition().getDomain();
 	};
 
 	/*! \brief Interpolate particles to mesh
@@ -651,23 +674,51 @@ public:
 
 #endif
 
-		Box<vector::dims,typename vector::stype> domain = vd.getDecomposition().getDomain();
+		// point position
+		typename vector::stype xp[vector::dims];
 
-		// grid spacing
-		typename vector::stype dx[vector::dims];
+		int ip[vector::dims][kernel::np];
+		typename vector::stype x[vector::dims][kernel::np];
+		typename vector::stype a[vector::dims][kernel::np];
 
-		for (size_t i = 0 ; i < vector::dims ; i++)
-			dx[i] = 1.0/gd.spacing(i);
+		typename vector::stype a_int[openfpm::math::pow(kernel::np,vector::dims)];
 
-		size_t sz[vector::dims];
+		auto it = vd.getDomainIterator();
 
-		for (size_t i = 0 ; i < vector::dims ; i++)
-			sz[i] = kernel::np;
+		while (it.isNext() == true)
+		{
+			auto key_p = it.get();
 
-		// Precalculate the offset for each sub-sub-domain
-		openfpm::vector<agg_arr<openfpm::math::pow(kernel::np,vector::dims)>> offsets;
+			inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_p2m,openfpm::math::pow(kernel::np,vector::dims)>(key_p,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
 
-		calculate_the_offsets(offsets,sz);
+			++it;
+		}
+	}
+
+	/*! \brief Interpolate mesh to particle
+	 *
+	 * Most of the time the particle set and the mesh are the same
+	 * as the one used in the constructor. They also can be different
+	 * as soon as they have the same decomposition
+	 *
+	 * \param gd grid or mesh
+	 * \param vd particle set
+	 *
+	 */
+	template<unsigned int prp_g, unsigned int prp_v> void m2p(grid & gd, vector & vd)
+	{
+#ifdef SE_CLASS1
+
+		if (!vd.getDecomposition().is_equal_ng(gd.getDecomposition()) )
+		{
+			std::cerr << __FILE__ << ":" << __LINE__ << " Error: the distribution of the vector of particles" <<
+					" and the grid is different. In order to interpolate the two data structure must have the" <<
+					" same decomposition" << std::endl;
+
+			ACTION_ON_ERROR(INTERPOLATION_ERROR_OBJECT)
+		}
+
+#endif
 
 		// point position
 		typename vector::stype xp[vector::dims];
@@ -676,19 +727,24 @@ public:
 		typename vector::stype x[vector::dims][kernel::np];
 		typename vector::stype a[vector::dims][kernel::np];
 
+		//		grid_cpu<vector::dims,aggregate<typename vector::stype>> a_int(sz);
+		//		a_int.setMemory();
 		typename vector::stype a_int[openfpm::math::pow(kernel::np,vector::dims)];
 
 		auto it = vd.getDomainIterator();
 
 		while (it.isNext() == true)
 		{
-			inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_p2m,openfpm::math::pow(kernel::np,vector::dims)>(it,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
+			auto key_p = it.get();
+
+			inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_m2p,openfpm::math::pow(kernel::np,vector::dims)>(key_p,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
 
 			++it;
 		}
 	}
 
-	/*! \brief Interpolate mesh to particle
+
+	/*! \brief Interpolate particles to mesh
 	 *
 	 * Most of the time the particle set and the mesh are the same
 	 * as the one used in the constructor. They also can be different
@@ -696,9 +752,10 @@ public:
 	 *
 	 * \param gd grid or mesh
 	 * \param vd particle set
+	 * \param p particle
 	 *
 	 */
-	template<unsigned int prp_g, unsigned int prp_v> void m2p(grid & gd, vector & vd)
+	template<unsigned int prp_v, unsigned int prp_g> inline void p2m(vector & vd, grid & gd,const vect_dist_key_dx & p)
 	{
 #ifdef SE_CLASS1
 
@@ -713,14 +770,6 @@ public:
 
 #endif
 
-		Box<vector::dims,typename vector::stype> domain = vd.getDecomposition().getDomain();
-
-		// grid spacing
-		typename vector::stype dx[vector::dims];
-
-		for (size_t i = 0 ; i < vector::dims ; i++)
-			dx[i] = 1.0/gd.spacing(i);
-
 		// point position
 		typename vector::stype xp[vector::dims];
 
@@ -728,28 +777,49 @@ public:
 		typename vector::stype x[vector::dims][kernel::np];
 		typename vector::stype a[vector::dims][kernel::np];
 
-		size_t sz[vector::dims];
+		typename vector::stype a_int[openfpm::math::pow(kernel::np,vector::dims)];
 
-		for (size_t i = 0 ; i < vector::dims ; i++)
-			sz[i] = kernel::np;
+		inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_p2m,openfpm::math::pow(kernel::np,vector::dims)>(p,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
+	}
 
-		// Precalculate the offset for each sub-sub-domain
-		openfpm::vector<agg_arr<openfpm::math::pow(kernel::np,vector::dims)>> offsets;
+	/*! \brief Interpolate mesh to particle
+	 *
+	 * Most of the time the particle set and the mesh are the same
+	 * as the one used in the constructor. They also can be different
+	 * as soon as they have the same decomposition
+	 *
+	 * \param gd grid or mesh
+	 * \param vd particle set
+	 * \param p particle
+	 *
+	 */
+	template<unsigned int prp_g, unsigned int prp_v> inline void m2p(grid & gd, vector & vd, const vect_dist_key_dx & p)
+	{
+#ifdef SE_CLASS1
 
-		calculate_the_offsets(offsets,sz);
+		if (!vd.getDecomposition().is_equal_ng(gd.getDecomposition()) )
+		{
+			std::cerr << __FILE__ << ":" << __LINE__ << " Error: the distribution of the vector of particles" <<
+					" and the grid is different. In order to interpolate the two data structure must have the" <<
+					" same decomposition" << std::endl;
+
+			ACTION_ON_ERROR(INTERPOLATION_ERROR_OBJECT)
+		}
+
+#endif
+
+		// point position
+		typename vector::stype xp[vector::dims];
+
+		int ip[vector::dims][kernel::np];
+		typename vector::stype x[vector::dims][kernel::np];
+		typename vector::stype a[vector::dims][kernel::np];
 
 		//		grid_cpu<vector::dims,aggregate<typename vector::stype>> a_int(sz);
 		//		a_int.setMemory();
 		typename vector::stype a_int[openfpm::math::pow(kernel::np,vector::dims)];
 
-		auto it = vd.getDomainIterator();
-
-		while (it.isNext() == true)
-		{
-			inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_m2p,openfpm::math::pow(kernel::np,vector::dims)>(it,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
-
-			++it;
-		}
+		inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_m2p,openfpm::math::pow(kernel::np,vector::dims)>(p,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
 	}
 
 };

src/interpolation/interpolation_unit_tests.hpp

@@ -82,8 +82,63 @@ template<typename vector, unsigned int mom_p> void momenta_vector(vector & vd,ty
 	}
 }
 
+template<unsigned int dim, typename T, typename grid, typename vector>
+void interp_test(grid & gd, vector & vd, bool single_particle)
+{
+	// Reset the grid
 
-BOOST_AUTO_TEST_CASE( interpolation_full_test_2D )
+	auto it2 = gd.getDomainGhostIterator();
+
+	while (it2.isNext())
+	{
+		auto key = it2.get();
+
+		gd.template get<0>(key) = 0.0;
+
+		++it2;
+	}
+
+	interpolate<vector,grid,mp4_kernel<float>> inte(vd,gd);
+
+	if (single_particle == false)
+	{inte.template p2m<0,0>(vd,gd);}
+	else
+	{
+		auto it = vd.getDomainIterator();
+
+		while (it.isNext())
+		{
+			auto p = it.get();
+
+			inte.template p2m<0,0>(vd,gd,p);
+
+			++it;
+		}
+	}
+
+	T mg[dim];
+	T mv[dim];
+
+	momenta_grid<grid,0>(gd,mg);
+	momenta_vector<vector,0>(vd,mv);
+
+	for (size_t i = 0 ; i < dim ; i++)
+	{BOOST_REQUIRE_CLOSE(mg[i],mv[i],0.001);}
+
+	momenta_grid<grid,1>(gd,mg);
+	momenta_vector<vector,1>(vd,mv);
+
+	for (size_t i = 0 ; i < dim ; i++)
+	{BOOST_REQUIRE_CLOSE(mg[i],mv[i],0.001);}
+
+	momenta_grid<grid,2>(gd,mg);
+	momenta_vector<vector,2>(vd,mv);
+
+	for (size_t i = 0 ; i < dim ; i++)
+	{BOOST_REQUIRE_CLOSE(mg[i],mv[i],0.001);}
+}
+
+BOOST_AUTO_TEST_CASE( interpolation_full_single_test_2D )
 {
 	Box<2,float> domain({0.0,0.0},{1.0,1.0});
 	size_t sz[2] = {64,64};
@@ -93,8 +148,7 @@ BOOST_AUTO_TEST_CASE( interpolation_full_test_2D )
 
 	size_t bc_v[2] = {PERIODIC,PERIODIC};
 
-	{
-	vector_dist<2,float,aggregate<float>> vd(4096,domain,bc_v,gv);
+	vector_dist<2,float,aggregate<float>> vd(65536,domain,bc_v,gv);
 	grid_dist_id<2,float,aggregate<float>> gd(vd.getDecomposition(),sz,gg);
 
 	// set one particle on vd
@@ -115,46 +169,52 @@ BOOST_AUTO_TEST_CASE( interpolation_full_test_2D )
 
 	vd.map();
 
-	// Reset the grid
+	interp_test<2,float>(gd,vd,true);
+}
 
-	auto it2 = gd.getDomainGhostIterator();
 
-	while (it2.isNext())
-	{
-		auto key = it2.get();
+BOOST_AUTO_TEST_CASE( interpolation_full_test_2D )
+{
+	Box<2,float> domain({0.0,0.0},{1.0,1.0});
+	size_t sz[2] = {64,64};
 
-		gd.template get<0>(key) = 0.0;
+	Ghost<2,long int> gg(2);
+	Ghost<2,float> gv(0.01);
 
-		++it2;
-	}
+	size_t bc_v[2] = {PERIODIC,PERIODIC};
 
-	interpolate<decltype(vd),decltype(gd),mp4_kernel<float>> inte(vd,gd);
+	{
+	vector_dist<2,float,aggregate<float>> vd(4096,domain,bc_v,gv);
+	grid_dist_id<2,float,aggregate<float>> gd(vd.getDecomposition(),sz,gg);
 
-	inte.p2m<0,0>(vd,gd);
+	// set one particle on vd
 
-	float mg[2];
-	float mv[2];
+	auto it = vd.getDomainIterator();
 
-	momenta_grid<decltype(gd),0>(gd,mg);
-	momenta_vector<decltype(vd),0>(vd,mv);
+	while (it.isNext())
+	{
+		auto p = it.get();
 
-	BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
-	BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
+		vd.getPos(p)[0] = (double)rand()/RAND_MAX;
+		vd.getPos(p)[1] = (double)rand()/RAND_MAX;
 
-	momenta_grid<decltype(gd),1>(gd,mg);
-	momenta_vector<decltype(vd),1>(vd,mv);
+		vd.getProp<0>(p) = 5.0;
 
-	BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
-	BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
+		++it;
+	}
 
-	momenta_grid<decltype(gd),2>(gd,mg);
-	momenta_vector<decltype(vd),2>(vd,mv);
+	vd.map();
 
-	BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
-	BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
+	// Reset the grid
+	interp_test<2,float>(gd,vd,false);
+
+	float mg[2];
+	float mv[2];
 
 	auto & v_cl = create_vcluster();
 
+	interpolate<decltype(vd),decltype(gd),mp4_kernel<float>> inte(vd,gd);
+
 	// We have to do a ghost get before interpolating m2p
 	// Before doing mesh to particle particle must be arranged
 	// into a grid like
@@ -248,8 +308,7 @@ BOOST_AUTO_TEST_CASE( interpolation_full_test_2D )
 	}
 }
 
-
-BOOST_AUTO_TEST_CASE( interpolation_full_test_3D )
+BOOST_AUTO_TEST_CASE( interpolation_full_single_test_3D )
 {
 	Box<3,double> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
 	size_t sz[3] = {64,64,64};
@@ -259,7 +318,6 @@ BOOST_AUTO_TEST_CASE( interpolation_full_test_3D )
 
 	size_t bc_v[3] = {PERIODIC,PERIODIC,PERIODIC};
 
-	{
 	vector_dist<3,double,aggregate<double>> vd(65536,domain,bc_v,gv);
 	grid_dist_id<3,double,aggregate<double>> gd(vd.getDecomposition(),sz,gg);
 
@@ -283,47 +341,51 @@ BOOST_AUTO_TEST_CASE( interpolation_full_test_3D )
 	vd.map();
 
 	// Reset the grid
+	interp_test<3,double>(gd,vd,true);
+}
 
-	auto it2 = gd.getDomainGhostIterator();
+BOOST_AUTO_TEST_CASE( interpolation_full_test_3D )
+{
+	Box<3,double> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
+	size_t sz[3] = {64,64,64};
 
-	while (it2.isNext())
-	{
-		auto key = it2.get();
+	Ghost<3,long int> gg(2);
+	Ghost<3,double> gv(0.01);
 
-		gd.template get<0>(key) = 0.0;
+	size_t bc_v[3] = {PERIODIC,PERIODIC,PERIODIC};
 
-		++it2;
-	}
+	{
+	vector_dist<3,double,aggregate<double>> vd(65536,domain,bc_v,gv);
+	grid_dist_id<3,double,aggregate<double>> gd(vd.getDecomposition(),sz,gg);
 
-	interpolate<decltype(vd),decltype(gd),mp4_kernel<double>> inte(vd,gd);
+	// set one particle on vd
 
-	inte.p2m<0,0>(vd,gd);
+	auto it = vd.getDomainIterator();
 
-	double mg[3];
-	double mv[3];
+	while (it.isNext())
+	{
+		auto p = it.get();
 
-	momenta_grid<decltype(gd),0>(gd,mg);
-	momenta_vector<decltype(vd),0>(vd,mv);
+		vd.getPos(p)[0] = (double)rand()/RAND_MAX;
+		vd.getPos(p)[1] = (double)rand()/RAND_MAX;
+		vd.getPos(p)[2] = (double)rand()/RAND_MAX;
 
-	BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
-	BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
-	BOOST_REQUIRE_CLOSE(mg[2],mv[2],0.001);
+		vd.getProp<0>(p) = 5.0;
 
-	momenta_grid<decltype(gd),1>(gd,mg);
-	momenta_vector<decltype(vd),1>(vd,mv);
+		++it;
+	}
 
-	BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
-	BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
-	BOOST_REQUIRE_CLOSE(mg[2],mv[2],0.001);
+	vd.map();
 
-	momenta_grid<decltype(gd),2>(gd,mg);
-	momenta_vector<decltype(vd),2>(vd,mv);
+	// Reset the grid
 
-	BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
-	BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
-	BOOST_REQUIRE_CLOSE(mg[2],mv[2],0.001);
+	// Reset the grid
+	interp_test<3,double>(gd,vd,false);
 
 	auto & v_cl = create_vcluster();
+	double mg[3];
+	double mv[3];
+	interpolate<decltype(vd),decltype(gd),mp4_kernel<double>> inte(vd,gd);
 
 	// We have to do a ghost get before interpolating m2p
 	// Before doing mesh to particle particle must be arranged