Fixing tests for numerics

configure.ac

@@ -122,13 +122,6 @@ if test x"$with_hdf5" = x"no"; then
     exit 207
 fi
 
-
-##########
-
-## Check for PETSC
-
-AX_LIB_PETSC()
-
 #########
 
 ###### Check for test coverage
@@ -314,6 +307,12 @@ AX_LAPACK([],[])
 
 AX_SUITESPARSE([],[])
 
+##########
+
+## Check for PETSC
+
+AX_LIB_PETSC()
+
 ###### Checking for EIGEN
 
 AX_EIGEN([],[])

m4/ax_petsc_lib.m4

@@ -105,7 +105,7 @@ AC_DEFUN([AX_LIB_PETSC], [
 			AX_OPENMP([CFLAGS="$OPENMP_CFLAGS"
 				   LDFLAGS="$OPENMP_LDFLAGS"],[])
                         CFLAGS="$CFLAGS -I$with_petsc/include $HDF5_INCLUDE $METIS_INCLUDE "
-                        LDFLAGS="$LDFLAGS -L$with_petsc/lib $HDF5_LDFLAGS  $HDF5_LIBS $METIS_LIB -lmetis "
+                        LDFLAGS="$LDFLAGS -L$with_petsc/lib $HDF5_LDFLAGS  $HDF5_LIBS $METIS_LIB -lmetis $SUITESPARSE_LIBS"
 			CC=$CXX
 
                         AC_LANG_SAVE

src/DMatrix/tests/EMatrix_unit_tests.cpp

+/*
+ * EMatrix_unit_tests.cpp
+ *
+ *  Created on: Feb 13, 2018
+ *      Author: i-bird
+ */
+
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+#include "DMatrix/EMatrix.hpp"
+#include "memory/HeapMemory.hpp"
+
+#ifdef HAVE_EIGEN
+
+BOOST_AUTO_TEST_SUITE (EMatrix_test)
+
+BOOST_AUTO_TEST_CASE( EMatrix_test_use)
+{
+	{
+	EMatrixXd em;
+
+	em.resize(8,5);
+
+	for (size_t i = 0 ; i < 8 ; i++)
+	{
+		for (size_t j = 0 ; j < 5 ; j++)
+		{em(i,j) = i*8+j;}
+	}
+
+	size_t pr = 0;
+	em.packRequest(pr);
+
+	// allocate the memory
+	HeapMemory pmem;
+	pmem.allocate(pr);
+	ExtPreAlloc<HeapMemory> & mem = *(new ExtPreAlloc<HeapMemory>(pr,pmem));
+	mem.incRef();
+
+	BOOST_REQUIRE_EQUAL(pr,8*5*sizeof(double) + 2*sizeof(size_t));
+
+	Pack_stat sts;
+	em.pack(mem,sts);
+
+	// Reset to zero
+	for (size_t i = 0 ; i < 8 ; i++)
+	{
+		for (size_t j = 0 ; j < 5 ; j++)
+		{em(i,j) = 0;}
+	}
+
+	Unpack_stat ps;
+	em.unpack(mem,ps);
+
+	for (size_t i = 0 ; i < 8 ; i++)
+	{
+		for (size_t j = 0 ; j < 5 ; j++)
+		{BOOST_REQUIRE_EQUAL(em(i,j),i*8+j);}
+	}
+	}
+
+
+	{
+	EMatrix3d em;
+
+	em.resize(3,3);
+
+	for (size_t i = 0 ; i < 3 ; i++)
+	{
+		for (size_t j = 0 ; j < 3 ; j++)
+		{em(i,j) = i*8+j;}
+	}
+
+	size_t pr = 0;
+	em.packRequest(pr);
+
+	// allocate the memory
+	HeapMemory pmem;
+	pmem.allocate(pr);
+	ExtPreAlloc<HeapMemory> & mem = *(new ExtPreAlloc<HeapMemory>(pr,pmem));
+	mem.incRef();
+
+	BOOST_REQUIRE_EQUAL(pr,3*3*sizeof(double) + 2*sizeof(size_t));
+
+	Pack_stat sts;
+	em.pack(mem,sts);
+
+	// Reset to zero
+	for (size_t i = 0 ; i < 3 ; i++)
+	{
+		for (size_t j = 0 ; j < 3 ; j++)
+		{em(i,j) = 0;}
+	}
+
+	Unpack_stat ps;
+	em.unpack(mem,ps);
+
+	for (size_t i = 0 ; i < 3 ; i++)
+	{
+		for (size_t j = 0 ; j < 3 ; j++)
+		{BOOST_REQUIRE_EQUAL(em(i,j),i*8+j);}
+	}
+	}
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+
+#endif

src/Draw/PointIteratorSkin.hpp

@@ -136,8 +136,15 @@ public:
 	 * \param sp grid spacing
 	 *
 	 */
-	PointIteratorSkin( Decomposition & dec, size_t (& sz)[dim], const Box<dim,T> & domain, const Box<dim,T> & sub_A, const Box<dim,T> & sub_B, size_t (& bc)[dim])
-	:grid_dist_id_iterator_dec_skin<Decomposition>(dec, sz, getAB(sz,domain,sub_A,sub_B,sp,RETURN_A), getAB(sz,domain,sub_A,sub_B,sp,RETURN_B), bc),domain(domain)
+	PointIteratorSkin( Decomposition & dec,
+			           size_t (& sz)[dim],
+					   const Box<dim,T> & domain,
+					   const Box<dim,T> & sub_A,
+					   const Box<dim,T> & sub_B,
+					   size_t (& bc)[dim])
+	:grid_dist_id_iterator_dec_skin<Decomposition>(dec, sz, getAB(sz,domain,sub_A,sub_B,sp,RETURN_A),
+	 getAB(sz,domain,sub_A,sub_B,sp,RETURN_B), bc),
+	 domain(domain)
 	{
 		sub_domainA.add(sub_A);
 		calculateAp();

src/FiniteDifference/FDScheme.hpp

@@ -441,6 +441,46 @@ private:
 		}
 	}
 
+	/*! \brief Impose an operator
+	 *
+	 * This function the RHS no matrix is imposed. This
+	 * function is usefull if the Matrix has been constructed and only
+	 * the right hand side b must be changed
+	 *
+	 * Ax = b
+	 *
+	 * \param num right hand side of the term (b term)
+	 * \param id Equation id in the system that we are imposing
+	 * \param it_d iterator that define where you want to impose
+	 *
+	 */
+	template<typename bop, typename iterator>
+	void impose_dit_b(bop num,
+					  long int id ,
+					  const iterator & it_d)
+	{
+
+		auto it = it_d;
+		grid_sm<Sys_eqs::dims,void> gs = g_map.getGridInfoVoid();
+
+		// iterate all the grid points
+		while (it.isNext())
+		{
+			// get the position
+			auto key = it.get();
+
+			b(g_map.template get<0>(key)*Sys_eqs::nvar + id) = num.get(key);
+
+			// if SE_CLASS1 is defined check the position
+#ifdef SE_CLASS1
+//			T::position(key,gs,s_pos);
+#endif
+
+			++row_b;
+			++it;
+		}
+	}
+
 	/*! \brief Impose an operator
 	 *
 	 * This function impose an operator on a particular grid region to produce the system
@@ -456,7 +496,7 @@ private:
 	 * \param it_d iterator that define where you want to impose
 	 *
 	 */
-	template<typename T, typename bop, typename iterator> void impose_dit(const T & op ,
+	template<typename T, typename bop, typename iterator> void impose_dit(const T & op,
 			                         bop num,
 									 long int id ,
 									 const iterator & it_d)
@@ -493,9 +533,8 @@ private:
 				trpl.last().value() = it->second;
 
 				if (trpl.last().row() == trpl.last().col())
-					is_diag = true;
+				{is_diag = true;}
 
-//				std::cout << "(" << trpl.last().row() << "," << trpl.last().col() << "," << trpl.last().value() << ")" << "\n";
 			}
 
 			// If does not have a diagonal entry put it to zero
@@ -812,6 +851,44 @@ public:
 
 	}
 
+	/*! \brief In case we want to impose a new b re-using FDScheme we have to call
+	 *         This function
+	 */
+	void new_b()
+	{row_b = 0;}
+
+	/*! \brief In case we want to impose a new A re-using FDScheme we have to call
+	 *         This function
+	 *
+	 */
+	void new_A()
+	{row = 0;}
+
+	/*! \brief Impose an operator
+	 *
+	 * This function impose an operator on a particular grid region to produce the system
+	 *
+	 * Ax = b
+	 *
+	 * ## Stokes equation 2D, lid driven cavity with one splipping wall
+	 * \snippet eq_unit_test.hpp Copy the solution to grid
+	 *
+	 * \param op Operator to impose (A term)
+	 * \param num right hand side of the term (b term)
+	 * \param id Equation id in the system that we are imposing
+	 * \param it_d iterator that define where you want to impose
+	 *
+	 */
+	template<unsigned int prp, typename b_term, typename iterator>
+	void impose_dit_b(b_term & b_t,
+					  const iterator & it_d,
+					  long int id = 0)
+	{
+		grid_b<b_term,prp> b(b_t);
+
+		impose_dit_b(b,id,it_d);
+	}
+
 	/*! \brief Impose an operator
 	 *
 	 * This function impose an operator on a particular grid region to produce the system

src/FiniteDifference/eq_unit_test_3d.cpp

@@ -229,6 +229,11 @@ template<typename solver_type,typename lid_nn_3d> void lid_driven_cavity_3d()
         std::string file1 = std::string("test/") + s + "lid_driven_cavity_3d_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + "_test_GCC6.vtk";
         std::string file2 = s + "lid_driven_cavity_3d_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk";
 
+    #elif __GNUC__ == 7
+
+        std::string file1 = std::string("test/") + s + "lid_driven_cavity_3d_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + "_test_GCC7.vtk";
+        std::string file2 = s + "lid_driven_cavity_3d_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk";
+
 	#else
 
         std::string file1 = std::string("test/") + s + "lid_driven_cavity_3d_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + "_test_GCC4.vtk";
@@ -254,7 +259,7 @@ template<typename solver_type,typename lid_nn_3d> void lid_driven_cavity_3d()
 
 BOOST_AUTO_TEST_CASE(lid_driven_cavity)
 {
-#ifdef HAVE_EIGEN
+#if defined(HAVE_EIGEN) && defined(HAVE_SUITESPARSE)
 	lid_driven_cavity_3d<umfpack_solver<double>,lid_nn_3d_eigen>();
 #endif
 #ifdef HAVE_PETSC

src/Makefile.am

@@ -2,8 +2,13 @@
 LINKLIBS = $(HDF5_LDFLAGS)  $(HDF5_LIBS) $(OPENMP_LDFLAGS) $(LIBHILBERT_LIB) $(PETSC_LIB) $(SUITESPARSE_LIBS) $(LAPACK_LIBS) $(BLAS_LIBS)  $(METIS_LIB) $(PARMETIS_LIB) $(DEFAULT_LIB) $(PTHREAD_LIBS) $(OPT_LIBS) $(BOOST_LDFLAGS) $(BOOST_PROGRAM_OPTIONS_LIB) $(BOOST_IOSTREAMS_LIB) $(LIBQUADMATH) $(OPENMP_LDFLAGS) $(LIBIFCORE)
 
 noinst_PROGRAMS = numerics
+<<<<<<< HEAD
 numerics_SOURCES = main.cpp Matrix/SparseMatrix_unit_tests.cpp interpolation/interpolation_unit_tests.cpp  Vector/Vector_unit_tests.cpp  Solvers/petsc_solver_unit_tests.cpp FiniteDifference/FDScheme_unit_tests.cpp FiniteDifference/eq_unit_test_3d.cpp FiniteDifference/eq_unit_test.cpp  Operators/Vector/vector_dist_operators_unit_tests.cpp ../../openfpm_vcluster/src/VCluster/VCluster.cpp ../../openfpm_devices/src/memory/HeapMemory.cpp ../../openfpm_devices/src/memory/PtrMemory.cpp ../../openfpm_devices/src/Memleak_check.cpp
 numerics_CXXFLAGS = -Wno-unknown-pragmas -Wno-int-in-bool-context $(HDF5_CPPFLAGS) $(OPENMP_CFLAGS) $(LIBHILBERT_INCLUDE) $(AM_CXXFLAGS) $(INCLUDES_PATH) $(BOOST_CPPFLAGS) $(SUITESPARSE_INCLUDE) $(METIS_INCLUDE) $(PARMETIS_INCLUDE) $(EIGEN_INCLUDE) $(PETSC_INCLUDE) -Wno-deprecated-declarations -Wno-unused-local-typedefs
+=======
+numerics_SOURCES = main.cpp Solvers/petsc_solver_unit_tests.cpp DMatrix/tests/EMatrix_unit_tests.cpp  ../../openfpm_vcluster/src/VCluster/VCluster.cpp ../../openfpm_devices/src/memory/HeapMemory.cpp ../../openfpm_devices/src/memory/PtrMemory.cpp ../../openfpm_devices/src/Memleak_check.cpp
+numerics_CXXFLAGS = $(HDF5_CPPFLAGS) $(OPENMP_CFLAGS) $(LIBHILBERT_INCLUDE) $(AM_CXXFLAGS) $(INCLUDES_PATH) $(BOOST_CPPFLAGS) $(SUITESPARSE_INCLUDE) $(METIS_INCLUDE) $(PARMETIS_INCLUDE) $(EIGEN_INCLUDE) $(PETSC_INCLUDE) -Wno-deprecated-declarations -Wno-unused-local-typedefs
+>>>>>>> 0b425037b139f9c2dbe989da5468eebb215bc623
 numerics_CFLAGS = $(CUDA_CFLAGS)
 numerics_LDADD = $(LINKLIBS) -lparmetis -lmetis
 nobase_include_HEADERS = Matrix/SparseMatrix.hpp Matrix/SparseMatrix_Eigen.hpp Matrix/SparseMatrix_petsc.hpp \
@@ -13,6 +18,7 @@ util/petsc_util.hpp util/linalgebra_lib.hpp util/util_num.hpp util/grid_dist_tes
 FiniteDifference/Average.hpp FiniteDifference/Derivative.hpp FiniteDifference/FD_util_include.hpp  FiniteDifference/eq.hpp FiniteDifference/FDScheme.hpp FiniteDifference/Laplacian.hpp FiniteDifference/mul.hpp FiniteDifference/sum.hpp FiniteDifference/util/common.hpp \
 PSE/Kernels.hpp PSE/Kernels_test_util.hpp Operators/Vector/vector_dist_operators_extensions.hpp Operators/Vector/vector_dist_operators.hpp Operators/Vector/vector_dist_operators_apply_kernel.hpp Operators/Vector/vector_dist_operators_functions.hpp Operators/Vector/vector_dist_operator_assign.hpp \
 Draw/DrawParticles.hpp Draw/PointIterator.hpp Draw/PointIteratorSkin.hpp \
+DMatrix/EMatrix.hpp \
 interpolation/interpolation.hpp interpolation/mp4_kernel.hpp interpolation/z_spline.hpp Solvers/petsc_solver_AMG_report.hpp
 
 .cu.o :

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/Operators/Vector/vector_dist_operators_functions.hpp

@@ -204,7 +204,6 @@ fun_name(double d, const vector_dist_expression_op<exp1,exp2,op1> & va)\
 
 CREATE_VDIST_ARG2_FUNC(pmul,pmul,VECT_PMUL)
 
-
 ////////// Special function reduce /////////////////////////
 
 
@@ -291,5 +290,17 @@ rsum(const vector_dist_expression<prp1,v1> & va, const vector_type & vd)
 	return exp_sum;
 }
 
+namespace openfpm
+{
+	/*! \brief General distance formula
+	 *
+	 *
+	 */
+	template <typename T, typename P> auto distance(T exp1, P exp2) -> decltype(norm(exp1 - exp2))
+	{
+		return norm(exp1 - exp2);
+	}
+}
+
 
 #endif /* OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_FUNCTIONS_HPP_ */

src/Operators/Vector/vector_dist_operators_unit_tests.cpp

@@ -940,7 +940,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_operators_test )
 
 	// normalization function
 
-	vA = vVB * vVC + norm(vVC + vVB) + distance(vVC,vVB);
+	vA = vVB * vVC + norm(vVC + vVB) + openfpm::distance(vVC,vVB);
 	check_values_scal_norm_dist(vd);
 
 	Point<3,float> p0({2.0,2.0,2.0});

src/PSE/Kernels_test_util.hpp

@@ -86,7 +86,7 @@ template<typename T, typename Kernel> void PSE_test(size_t Npart, size_t overlap
     size_t bc[1]={NON_PERIODIC};
 	Ghost<1,T> g(20*eps);
 
-	vector_dist<1,T, aggregate<T>, CartDecomposition<1,T> > vd(Npart,box,bc,g);
+	vector_dist<1,T, aggregate<T> > vd(Npart,box,bc,g);
 
 	auto it2 = vd.getIterator();
 

src/Solvers/umfpack_solver.hpp

@@ -14,7 +14,7 @@
 #define SOLVER_PRINT_RESIDUAL_NORM_INFINITY 1
 #define SOLVER_PRINT_DETERMINANT 2
 
-#ifdef HAVE_EIGEN
+#if defined(HAVE_EIGEN) && defined(HAVE_SUITESPARSE)
 
 /////// Compiled with EIGEN support
 

src/Vector/Vector_eigen.hpp

@@ -75,6 +75,16 @@ public:
 		row() = i;
 		value() = val;
 	}
+
+	/*! \brief Indicate that the structure has no pointer
+	 *
+	 * \return true
+	 *
+	 */
+	static inline bool noPointers()
+	{
+		return true;
+	}
 };
 
 template<typename T>

src/Vector/Vector_petsc.hpp

@@ -79,6 +79,16 @@ public:
 		rw() = i;
 		val() = val;
 	}
+
+	/*! \brief Indicate that the structure has no pointer
+	 *
+	 * \return true
+	 *
+	 */
+	static inline bool noPointers()
+	{
+		return true;
+	}
 };
 
 constexpr unsigned int row_id = 0;

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
@@ -434,24 +433,24 @@ struct inte_calc_impl
 		// calculate the position of the particle in the grid
 		// coordinated
 		for (size_t i = 0 ; i < vector::dims ; i++)
-			x0[i] = (p.get(i)-domain.getLow(i))*dx[i];
+		{x0[i] = (p.get(i)-domain.getLow(i))*dx[i];}
 
 		// convert into integer
 		for (size_t i = 0 ; i < vector::dims ; i++)
-			ip[i][0] = (int)x0[i];
+		{ip[i][0] = (int)x0[i];}
 
 		// convert the global grid position into local grid position
 		grid_key_dx<vector::dims> base;
 
 		for (size_t i = 0 ; i < vector::dims ; i++)
-			base.set_d(i,ip[i][0] - gd.getLocalGridsInfo().get(sub).origin.get(i) - (long int)kernel::np/2 + 1);
+		{base.set_d(i,ip[i][0] - gd.getLocalGridsInfo().get(sub).origin.get(i) - (long int)kernel::np/2 + 1);}
 
 		// convenient grid of points
 
 		for (size_t j = 0 ; j < kernel::np-1 ; j++)
 		{
 			for (size_t i = 0 ; i < vector::dims ; i++)
-				ip[i][j+1] = (int)ip[i][j]+1;
+			{ip[i][j+1] = (int)ip[i][j]+1;}
 		}
 
 		for (size_t i = 0 ; i < vector::dims ; i++)
@@ -460,13 +459,13 @@ struct inte_calc_impl
 		for (long int j = 0 ; j < kernel::np ; j++)
 		{
 			for (size_t i = 0 ; i < vector::dims ; i++)
-				x[i][j] = - xp[i] + typename vector::stype((long int)j - (long int)kernel::np/2 + 1);
+			{x[i][j] = - xp[i] + typename vector::stype((long int)j - (long int)kernel::np/2 + 1);}
 		}
 
 		for (size_t j = 0 ; j < kernel::np ; j++)
 		{
 			for (size_t i = 0 ; i < vector::dims ; i++)
-				a[i][j] = kernel::value(x[i][j],j);
+			{a[i][j] = kernel::value(x[i][j],j);}
 		}
 
 		calculate_aint<vector::dims,vector,kernel::np>::value(sz,a_int,a);
@@ -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;
 		}
 	}