From 7857bce369c4bfe77e98de997c93da9104a0e5c8 Mon Sep 17 00:00:00 2001
From: Pietro Incardona <incardon@mpi-cbg.de>
Date: Sat, 3 Jun 2017 13:17:35 +0200
Subject: [PATCH] Release 0.9.0

---
 CHANGELOG.md                                  |   4 +-
 example/Grid/3_gray_scott/main.cpp            |   6 +-
 example/Grid/3_gray_scott_3d/main.cpp         |  22 +-
 .../0_2D_incompressible/main_petsc.cpp        |   6 +-
 example/Vector/1_celllist/main.cpp            |   2 +-
 example/Vector/3_molecular_dynamic/main.cpp   |  84 +--
 .../4_multiphase_celllist_verlet/main.cpp     |  11 +-
 example/Vector/4_reorder/energy_force.hpp     |   4 +-
 example/Vector/4_reorder/main_comp_ord.cpp    |   4 +-
 .../5_molecular_dynamic_sym_crs/main.cpp      |   4 +-
 example/Vector/7_SPH_dlb_opt/main.cpp         |   9 +-
 openfpm_data                                  |   2 +-
 openfpm_devices                               |   2 +-
 openfpm_io                                    |   2 +-
 openfpm_vcluster                              |   2 +-
 src/Decomposition/CartDecomposition.hpp       |  12 +-
 src/Decomposition/common.hpp                  |   2 +-
 src/Decomposition/ie_ghost.hpp                |  48 +-
 src/Decomposition/ie_loc_ghost.hpp            |  63 +-
 src/Grid/grid_dist_id.hpp                     | 296 ++------
 src/Grid/grid_dist_id_comm.hpp                | 679 ++++++++++++++----
 src/Grid/grid_dist_id_unit_test.cpp           | 131 ++++
 src/Grid/grid_dist_key.hpp                    |  19 +
 src/Grid/grid_dist_util.hpp                   |  10 +-
 src/Makefile.am                               |   2 +-
 src/Vector/se_class3_vector.hpp               |  53 +-
 src/Vector/vector_dist.hpp                    | 129 +++-
 src/Vector/vector_dist_cell_list_tests.hpp    | 110 +++
 src/Vector/vector_dist_key.hpp                |   6 +
 29 files changed, 1170 insertions(+), 554 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2d2795d5..df303bd7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -20,8 +20,10 @@ All notable changes to this project will be documented in this file.
 - In case of miss compilation ignore system wide installation
 - Bug in VTK writer binary in case of vectors
 - Bug in VTK writer binary: long int are not supported removing output
+- Bug in the constructor with stencil bigger than one
 
-
+### Changed
+- CellList types has changed
 
 ## [0.8.0] 28 February 2016
 
diff --git a/example/Grid/3_gray_scott/main.cpp b/example/Grid/3_gray_scott/main.cpp
index 932a15a3..3bb0eff8 100644
--- a/example/Grid/3_gray_scott/main.cpp
+++ b/example/Grid/3_gray_scott/main.cpp
@@ -10,10 +10,10 @@
  * This example show the usage of periodic grid with ghost part given in grid units to solve
  * the following system of equations
  *
- * \f$\frac{\partial u}{\partial t} = D_u \nabla u -uv^2 + F(1-u)\f$
+ * \f$\frac{\partial u}{\partial t} = D_u \nabla u - uv^2 + F(1-u)\f$
  *
  *
- * \f$\frac{\partial v}{\partial t} = D_v \nabla v -uv^2 - (F + k)v\f$
+ * \f$\frac{\partial v}{\partial t} = D_v \nabla v + uv^2 - (F + k)v\f$
  * 
  * ## Constants and functions ##
  *
@@ -201,7 +201,7 @@ int main(int argc, char* argv[])
 	grid_dist_id<2, double, aggregate<double,double>> Old(sz,domain,g,bc);
 
 	// New grid with the decomposition of the old grid
-	grid_dist_id<2, double, aggregate<double,double>> New(Old.getDecomposition(),sz,domain,g);
+	grid_dist_id<2, double, aggregate<double,double>> New(Old.getDecomposition(),sz,g);
 
 	
 	// spacing of the grid on x and y
diff --git a/example/Grid/3_gray_scott_3d/main.cpp b/example/Grid/3_gray_scott_3d/main.cpp
index 73ea1646..91ec575f 100644
--- a/example/Grid/3_gray_scott_3d/main.cpp
+++ b/example/Grid/3_gray_scott_3d/main.cpp
@@ -46,8 +46,16 @@ void init(grid_dist_id<3,double,aggregate<double,double> > & Old, grid_dist_id<3
 		++it;
 	}
 
-	grid_key_dx<3> start({(long int)std::floor(Old.size(0)*1.55f/domain.getHigh(0)),(long int)std::floor(Old.size(1)*1.55f/domain.getHigh(1)),(long int)std::floor(Old.size(1)*1.55f/domain.getHigh(2))});
-	grid_key_dx<3> stop ({(long int)std::ceil (Old.size(0)*1.85f/domain.getHigh(0)),(long int)std::ceil (Old.size(1)*1.85f/domain.getHigh(1)),(long int)std::floor(Old.size(1)*1.85f/domain.getHigh(1))});
+        long int x_start = Old.size(0)*1.55f/domain.getHigh(0);
+        long int y_start = Old.size(1)*1.55f/domain.getHigh(1);
+        long int z_start = Old.size(1)*1.55f/domain.getHigh(2);
+
+        long int x_stop = Old.size(0)*1.85f/domain.getHigh(0);
+        long int y_stop = Old.size(1)*1.85f/domain.getHigh(1);
+        long int z_stop = Old.size(1)*1.85f/domain.getHigh(2);
+
+        grid_key_dx<3> start({x_start,y_start,z_start});
+        grid_key_dx<3> stop ({x_stop,y_stop,z_stop});
 	auto it_init = Old.getSubDomainIterator(start,stop);
 
 	while (it_init.isNext())
@@ -72,7 +80,7 @@ int main(int argc, char* argv[])
 	Box<3,double> domain({0.0,0.0},{2.5,2.5,2.5});
 	
 	// grid size
-	size_t sz[3] = {128,128,128};
+        size_t sz[3] = {128,128,128};
 
 	// Define periodicity of the grid
 	periodicity<3> bc = {PERIODIC,PERIODIC,PERIODIC};
@@ -90,11 +98,11 @@ int main(int argc, char* argv[])
 	double dv = 1*1e-5;
 
 	// Number of timesteps
-        size_t timeSteps = 17000;
+        size_t timeSteps = 5000;
 
 	// K and F (Physical constant in the equation)
-        double K = 0.065;
-        double F = 0.034;
+        double K = 0.014;
+        double F = 0.053;
 
 	//! \cond [init lib] \endcond
 
@@ -186,7 +194,7 @@ int main(int argc, char* argv[])
 		// visualization
 		if (i % 60 == 0)
 		{
-			Old.write("output",count,VTK_WRITER | FORMAT_BINARY);
+			Old.write_frame("output",count,VTK_WRITER | FORMAT_BINARY);
 			count++;
 		}
 	}
diff --git a/example/Numerics/Stoke_flow/0_2D_incompressible/main_petsc.cpp b/example/Numerics/Stoke_flow/0_2D_incompressible/main_petsc.cpp
index 27b9cbf7..6f5ae88b 100644
--- a/example/Numerics/Stoke_flow/0_2D_incompressible/main_petsc.cpp
+++ b/example/Numerics/Stoke_flow/0_2D_incompressible/main_petsc.cpp
@@ -404,7 +404,7 @@ int main(int argc, char* argv[])
 	 *
 	 * Once we have the solution we copy it on the grid
 	 *
-	 * \snippet Numerics/Stoke_flow/0_2D_incompressible/main_eigen.cpp copy write
+	 * \snippet Numerics/Stoke_flow/0_2D_incompressible/main_petsc.cpp copy write
 	 *
 	 */
 
@@ -424,7 +424,7 @@ int main(int argc, char* argv[])
 	 *
 	 *  At the very end of the program we have always to de-initialize the library
 	 *
-	 * \snippet Numerics/Stoke_flow/0_2D_incompressible/main_eigen.cpp fin lib
+	 * \snippet Numerics/Stoke_flow/0_2D_incompressible/main_petsc.cpp fin lib
 	 *
 	 */
 
@@ -440,7 +440,7 @@ int main(int argc, char* argv[])
 	 *
 	 * # Full code # {#num_sk_inc_2D_ps_code}
 	 *
-	 * \include Numerics/Stoke_flow/0_2D_incompressible/main_eigen.cpp
+	 * \include Numerics/Stoke_flow/0_2D_incompressible/main_petsc.cpp
 	 *
 	 */
 }
diff --git a/example/Vector/1_celllist/main.cpp b/example/Vector/1_celllist/main.cpp
index 28c4a6fe..c244c468 100644
--- a/example/Vector/1_celllist/main.cpp
+++ b/example/Vector/1_celllist/main.cpp
@@ -247,7 +247,7 @@ int main(int argc, char* argv[])
 		Point<3,float> xp = vd.getPos(p);
 
 		// Get an iterator of all the particles neighborhood of p
-		auto Np = NN.getIterator(NN.getCell(vd.getPos(p)));
+		auto Np = NN.getNNIterator(NN.getCell(vd.getPos(p)));
 
 		// For each particle near p
 		while (Np.isNext())
diff --git a/example/Vector/3_molecular_dynamic/main.cpp b/example/Vector/3_molecular_dynamic/main.cpp
index 8d70b8c9..7a2e86db 100644
--- a/example/Vector/3_molecular_dynamic/main.cpp
+++ b/example/Vector/3_molecular_dynamic/main.cpp
@@ -6,8 +6,6 @@
  */
 
 #include "Vector/vector_dist.hpp"
-#include "Decomposition/CartDecomposition.hpp"
-#include "data_type/aggregate.hpp"
 #include "Plot/GoogleChart.hpp"
 #include "Plot/util.hpp"
 #include "timer.hpp"
@@ -19,7 +17,10 @@
  *
  * # Molecular Dynamic with Lennard-Jones potential # {#e3_md}
  *
- * This example show a simple Lennard-Jones molecular dynamic simulation in a stable regime
+ * This example show a simple Lennard-Jones molecular dynamic simulation in a stable regime.
+ * Particle feel each other by the potential.
+ *
+ * \f$ V(x_p,x_q) = 4( (\frac{\sigma}{r})^{12} - (\frac{\sigma}{r})^6  ) \f$
  *
  * ## Constants ##
  *
@@ -55,7 +56,7 @@ constexpr int force = 1;
 
 //! \cond [calc forces] \endcond
 
-void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6, double r_cut2)
+template<typename CellList> void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList & NN, double sigma12, double sigma6, double r_cut2)
 {
 
 //! \cond [calc forces] \endcond
@@ -81,12 +82,12 @@ void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, Ce
 
 	/*!
 	 *
-	 * \page Vector_3_md Vector 3 molecular dynamic with cell-list
+	 * \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
 	 *
 	 * Get an iterator over the particles and get its position. For each particle p iterate in its neighborhood q
 	 * and calculate the force based on the Lennard-Jhones potential given by
 	 *
-	 * \f$ F(x_p,x_q) = 24(\frac{2}{r^{13}} - \frac{1}{r^{7}}) r \f$
+	 * \f$ F(x_p,x_q) = 24( \frac{2 \sigma^{12}}{r^{13}} - \frac{\sigma^6}{r^{7}}) \hat{r} \f$
 	 *
 	 * \see \ref e0_s_assign_pos
 	 *
@@ -108,7 +109,7 @@ void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, Ce
 		// Get the position xp of the particle
 		Point<3,double> xp = vd.getPos(p);
 
-		// Reset the forice counter
+		// Reset the force counter
 		vd.template getProp<force>(p)[0] = 0.0;
 		vd.template getProp<force>(p)[1] = 0.0;
 		vd.template getProp<force>(p)[2] = 0.0;
@@ -175,7 +176,7 @@ void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, Ce
 
 //! \cond [calc energy] \endcond
 
-double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6, double r_cut2)
+template<typename CellList> double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList & NN, double sigma12, double sigma6, double r_cut2)
 {
         double rc = r_cut2;
         double shift = 2.0 * ( sigma12 / (rc*rc*rc*rc*rc*rc) - sigma6 / ( rc*rc*rc) );
@@ -208,7 +209,7 @@ double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd,
 	 * First we get an iterator over the particles and get its position. For each particle p iterate in its neighborhood q
 	 * and calculate the energy based on the Lennard-Jhones potential given by
 	 *
-	 * \f$ V(x_p,x_q) = 4(\frac{1}{r^{12}} - \frac{1}{r^{6}}) r \f$
+	 * \f$ V(x_p,x_q) = 4(\frac{1}{r^{12}} - \frac{1}{r^{6}}) \f$
 	 *
 	 * \see \ref e0_s_assign_pos
 	 *
@@ -284,43 +285,24 @@ int main(int argc, char* argv[])
 	 *
 	 * ## Initialization ## {#e3_md_init}
 	 *
-	 * After we defined the two main function calc forces and calc energy, we define
-	 *  important parameters of the simulation, time step integration,
+	 * After we defined the two main function calc forces and calc energy, we Initialize
+	 *  the library, we create a Box that define our domain, boundary conditions and ghost.
+	 *  Than we define important parameters of the simulation, time step integration,
 	 * size of the box, and cut-off radius of the interaction. We also define 2 vectors
 	 * x and y (they are like std::vector) used for statistic
 	 *
+	 * \see \ref e0_s_init
+	 *
 	 * \snippet Vector/3_molecular_dynamic/main.cpp constants run
 	 *
 	 */
 
 	//! \cond [constants run] \endcond
 
-	double dt = 0.0005;
+	openfpm_init(&argc,&argv);
+
 	double sigma = 0.1;
 	double r_cut = 3.0*sigma;
-	double sigma12 = pow(sigma,12);
-	double sigma6 = pow(sigma,6);
-
-	openfpm::vector<double> x;
-	openfpm::vector<openfpm::vector<double>> y;
-
-	//! \cond [constants run] \endcond
-
-	/*!
-	 * \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
-	 *
-	 * Here we Initialize the library, we create a Box that define our domain, boundary conditions and ghost
-	 *
-	 * \see \ref e0_s_init
-	 *
-	 * \snippet Vector/3_molecular_dynamic/main.cpp init
-	 *
-	 */
-
-	//! \cond [init] \endcond
-
-	openfpm_init(&argc,&argv);
-	Vcluster & v_cl = create_vcluster();
 
 	// we will use it do place particles on a 10x10x10 Grid like
 	size_t sz[3] = {10,10,10};
@@ -334,7 +316,14 @@ int main(int argc, char* argv[])
 	// ghost, big enough to contain the interaction radius
 	Ghost<3,float> ghost(r_cut);
 
-	//! \cond [init] \endcond
+	double dt = 0.0005;
+	double sigma12 = pow(sigma,12);
+	double sigma6 = pow(sigma,6);
+
+	openfpm::vector<double> x;
+	openfpm::vector<openfpm::vector<double>> y;
+
+	//! \cond [constants run] \endcond
 
 	/*!
 	 * \page Vector_3_md_dyn Vector 3 molecular dynamic with cell-list
@@ -370,18 +359,24 @@ int main(int argc, char* argv[])
 
 	//! \cond [vect grid] \endcond
 
+	// We create the grid iterator
 	auto it = vd.getGridIterator(sz);
 
 	while (it.isNext())
 	{
+		// Create a new particle
 		vd.add();
 
+		// key contain (i,j,k) index of the grid
 		auto key = it.get();
 
+		// The index of the grid can be accessed with key.get(0) == i, key.get(1) == j ...
+		// We use getLastPos to set the position of the last particle added
 		vd.getLastPos()[0] = key.get(0) * it.getSpacing(0);
 		vd.getLastPos()[1] = key.get(1) * it.getSpacing(1);
 		vd.getLastPos()[2] = key.get(2) * it.getSpacing(2);
 
+		// We use getLastProp to set the property value of the last particle we added
 		vd.template getLastProp<velocity>()[0] = 0.0;
 		vd.template getLastProp<velocity>()[1] = 0.0;
 		vd.template getLastProp<velocity>()[2] = 0.0;
@@ -404,8 +399,8 @@ int main(int argc, char* argv[])
 	 *
 	 * The verlet integration stepping look like this
 	 *
-	 * \f[ \vec{v}(t_{n+1/2}) = \vec{v}_p(t_n) + \frac{1}{2} \delta t \vec{a}(t_n) \f]
-	 * \f[ \vec{x}(t_{n}) = \vec{x}_p(t_n) + \delta t \vec{v}(t_n+1/2) \f]
+	 * \f[ \vec{v}(t_{n}+1/2) = \vec{v}_p(t_n) + \frac{1}{2} \delta t \vec{a}(t_n) \f]
+	 * \f[ \vec{x}(t_{n}+1) = \vec{x}_p(t_n) + \delta t \vec{v}(t_n+1/2) \f]
 	 *
 	 * calculate the forces from \f$ \vec{a} (t_{n}) \f$ finally
 	 *
@@ -465,7 +460,7 @@ int main(int argc, char* argv[])
 			++it3;
 		}
 
-		// Because we mooved the particles in space we have to map them and re-sync the ghost
+		// Because we moved the particles in space we have to map them and re-sync the ghost
 		vd.map();
 		vd.template ghost_get<>();
 
@@ -488,7 +483,7 @@ int main(int argc, char* argv[])
 			++it4;
 		}
 
-		// After every iteration collect some statistic about the confoguration
+		// After every iteration collect some statistic about the configuration
 		if (i % 100 == 0)
 		{	
 			// We write the particle position for visualization (Without ghost)
@@ -554,6 +549,15 @@ int main(int argc, char* argv[])
 	// width of the line
 	options.lineWidth = 1.0;
 
+	// Resolution in x
+	options.width = 1280;
+
+	// Resolution in y
+	options.heigh = 720;
+
+	// Add zoom capability
+	options.more = GC_ZOOM;
+
 	// Object that draw the X Y graph
 	GoogleChart cg;
 
diff --git a/example/Vector/4_multiphase_celllist_verlet/main.cpp b/example/Vector/4_multiphase_celllist_verlet/main.cpp
index de4093d5..cfa51c31 100644
--- a/example/Vector/4_multiphase_celllist_verlet/main.cpp
+++ b/example/Vector/4_multiphase_celllist_verlet/main.cpp
@@ -149,6 +149,8 @@ int main(int argc, char* argv[])
 
 	//! \cond [create multi-phase verlet] \endcond
 
+	{
+
 	// Get the cell list of the phase1
 	auto CL_phase1 = phases.get(1).getCellList(r_cut);
 
@@ -204,6 +206,8 @@ int main(int argc, char* argv[])
 		++it;
 	}
 
+	}
+
 	//! \cond [count part from phase0 to 1] \endcond
 
 	/*!
@@ -311,13 +315,14 @@ int main(int argc, char* argv[])
 	 *
 	 */
 
+	{
 	//! \cond [compute sym multi-phase two phase] \endcond
 
 	// Get the cell list of the phase1
-	CL_phase1 = phases.get(1).getCellListSym(r_cut);
+	auto CL_phase1 = phases.get(1).getCellListSym(r_cut);
 
 	// This function create a Verlet-list between phases 0 and 1
-	NN_ver01 = createVerletSym(phases.get(0),phases.get(1),CL_phase1,r_cut);
+	auto NN_ver01 = createVerletSym(phases.get(0),phases.get(1),CL_phase1,r_cut);
 
 	// Get an iterator over the real and ghost particles
 	it = phases.get(0).getDomainAndGhostIterator();
@@ -371,6 +376,8 @@ int main(int argc, char* argv[])
 	phases.get(0).ghost_put<add_,0>();
 	phases.get(1).ghost_put<add_,0>();
 
+	}
+
 	//! \cond [compute sym multi-phase two phase] \endcond
 
 	/*!
diff --git a/example/Vector/4_reorder/energy_force.hpp b/example/Vector/4_reorder/energy_force.hpp
index 5836b0fc..2108b811 100644
--- a/example/Vector/4_reorder/energy_force.hpp
+++ b/example/Vector/4_reorder/energy_force.hpp
@@ -12,7 +12,7 @@ constexpr int velocity = 0;
 constexpr int force = 1;
 
 
-void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6)
+template<typename CellList> void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList & NN, double sigma12, double sigma6)
 {
 	// update the Cell-list
 	vd.updateCellList(NN);
@@ -73,7 +73,7 @@ void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, Ce
 }
 
 
-double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6)
+template<typename CellList> double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList & NN, double sigma12, double sigma6)
 {
 	double E = 0.0;
 
diff --git a/example/Vector/4_reorder/main_comp_ord.cpp b/example/Vector/4_reorder/main_comp_ord.cpp
index 9a9a1abf..77f915a4 100644
--- a/example/Vector/4_reorder/main_comp_ord.cpp
+++ b/example/Vector/4_reorder/main_comp_ord.cpp
@@ -45,7 +45,7 @@ constexpr int force = 1;
 
 //! \cond [calc forces] \endcond
 
-void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList_hilb<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6)
+template<typename CellList> void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList & NN, double sigma12, double sigma6)
 {
 
 //! \cond [calc forces] \endcond
@@ -140,7 +140,7 @@ void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, Ce
 
 //! \cond [calc energy all] \endcond
 
-double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList_hilb<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6)
+template<typename CellList> double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList & NN, double sigma12, double sigma6)
 {
 	double E = 0.0;
 
diff --git a/example/Vector/5_molecular_dynamic_sym_crs/main.cpp b/example/Vector/5_molecular_dynamic_sym_crs/main.cpp
index fe6bec6b..a8eb5615 100644
--- a/example/Vector/5_molecular_dynamic_sym_crs/main.cpp
+++ b/example/Vector/5_molecular_dynamic_sym_crs/main.cpp
@@ -148,7 +148,7 @@ void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, Ve
 	//! \cond [real and ghost] \endcond
 
 	// Get an iterator over particles
-	auto it2 = vd.getParticleIteratorCRS(NN.getInternalCellList());
+	auto it2 = vd.getParticleIteratorCRS(NN);
 
 	//! \cond [real and ghost] \endcond
 
@@ -250,7 +250,7 @@ double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd,
 	double shift = 4.0 * ( sigma12 / (rc*rc*rc*rc*rc*rc) - sigma6 / ( rc*rc*rc) );
 
 	// Get an iterator over particles
-	auto it2 = vd.getParticleIteratorCRS(NN.getInternalCellList());
+	auto it2 = vd.getParticleIteratorCRS(NN);
 
 	// For each particle ...
 	while (it2.isNext())
diff --git a/example/Vector/7_SPH_dlb_opt/main.cpp b/example/Vector/7_SPH_dlb_opt/main.cpp
index cbb40805..6c54fbd9 100644
--- a/example/Vector/7_SPH_dlb_opt/main.cpp
+++ b/example/Vector/7_SPH_dlb_opt/main.cpp
@@ -912,7 +912,7 @@ int main(int argc, char* argv[])
 	cbar = coeff_sound * sqrt(gravity * h_swl);
 
 	// for each particle inside the fluid box ...
-/*	while (fluid_it.isNext())
+        while (fluid_it.isNext())
 	{
 		// ... add a particle ...
 		vd.add();
@@ -946,8 +946,7 @@ int main(int argc, char* argv[])
 
 		// next fluid particle
 		++fluid_it;
-	}*/
-
+        }
 
 	// Recipient
 	Box<3,double> recipient1({0.0,0.0,0.0},{1.6+dp/2.0,0.67+dp/2.0,0.4+dp/2.0});
@@ -1010,10 +1009,6 @@ int main(int argc, char* argv[])
 
 	vd.map();
 
-        vd.write("Recipient");
-
-	return 0;
-
 	// Now that we fill the vector with particles
 	ModelCustom md;
 
diff --git a/openfpm_data b/openfpm_data
index 3a732c98..2c313203 160000
--- a/openfpm_data
+++ b/openfpm_data
@@ -1 +1 @@
-Subproject commit 3a732c98d6404fd81529c25af4fbf33f8ca3de9c
+Subproject commit 2c31320305af04b0621d2a10818354d064ec1bbf
diff --git a/openfpm_devices b/openfpm_devices
index add2acfe..9817279e 160000
--- a/openfpm_devices
+++ b/openfpm_devices
@@ -1 +1 @@
-Subproject commit add2acfe9623e6045da1e41adbab6dc57ef50e5f
+Subproject commit 9817279ee3ac36b435d2178e994ab314a01aab30
diff --git a/openfpm_io b/openfpm_io
index a6bd0b79..9b0500dd 160000
--- a/openfpm_io
+++ b/openfpm_io
@@ -1 +1 @@
-Subproject commit a6bd0b797704f6faba662a10f5b101c218257490
+Subproject commit 9b0500dd1347d9f3d2898f2a9206220eea7907d5
diff --git a/openfpm_vcluster b/openfpm_vcluster
index 1856e665..7030e7f7 160000
--- a/openfpm_vcluster
+++ b/openfpm_vcluster
@@ -1 +1 @@
-Subproject commit 1856e665949ecca8dbcf68e08984aace853aad8d
+Subproject commit 7030e7f76e5d80a297e4cc9b2d4f3ad98beb8d0f
diff --git a/src/Decomposition/CartDecomposition.hpp b/src/Decomposition/CartDecomposition.hpp
index c097931d..728292d0 100755
--- a/src/Decomposition/CartDecomposition.hpp
+++ b/src/Decomposition/CartDecomposition.hpp
@@ -135,7 +135,7 @@ protected:
 	//! Structure that store the cartesian grid information
 	grid_sm<dim, void> gr_dist;
 
-	//! Structure that decompose your structure into cell without creating them
+	//! Structure that decompose the space into cells without creating them
 	//! useful to convert positions to CellId or sub-domain id in this case
 	CellDecomposer_sm<dim, T, shift<dim,T>> cd;
 
@@ -1724,6 +1724,16 @@ public:
 		return dist.get_ndec();
 	}
 
+	/*! \brief Get the cell decomposer of the decomposition
+	 *
+	 * \return the cell decomposer
+	 *
+	 */
+	const CellDecomposer_sm<dim, T, shift<dim,T>> & getCellDecomposer()
+	{
+		return cd;
+	}
+
 	//! friend classes
 	friend extended_type;
 
diff --git a/src/Decomposition/common.hpp b/src/Decomposition/common.hpp
index 36575f12..f9cdfeb6 100755
--- a/src/Decomposition/common.hpp
+++ b/src/Decomposition/common.hpp
@@ -114,7 +114,7 @@ struct lBox_dom
 {
 	//! Intersection between the local sub-domain enlarged by the ghost and the contiguous processor
 	//! sub-domains (External ghost)
-	openfpm::vector_std< Box_sub<dim,T> > ebx;
+	openfpm::vector_std< Box_sub_k<dim,T> > ebx;
 
 	//! Intersection between the contiguous processor sub-domain enlarged by the ghost with the
 	//! local sub-domain (Internal ghost)
diff --git a/src/Decomposition/ie_ghost.hpp b/src/Decomposition/ie_ghost.hpp
index cf748514..a7b14d91 100755
--- a/src/Decomposition/ie_ghost.hpp
+++ b/src/Decomposition/ie_ghost.hpp
@@ -991,7 +991,7 @@ public:
 		return true;
 	}
 
-	/*! \brief Check if the ie_loc_ghosts contain the same information
+	/*! \brief Check if the ie_ghosts contain the same information
 	 *
 	 * \param ig Element to check
 	 *
@@ -1065,52 +1065,6 @@ public:
 	 */
 	bool is_equal_ng(ie_ghost<dim,T> & ig)
 	{
-		Box<dim,T> bt;
-
-		if (getNEGhostBox() != ig.getNEGhostBox())
-			return false;
-
-		if (getNIGhostBox() != ig.getNIGhostBox())
-			return false;
-
-		for (size_t i = 0 ; i < proc_int_box.size() ; i++)
-		{
-			if (getProcessorNIGhost(i) != ig.getProcessorNIGhost(i))
-				return false;
-			for (size_t j = 0 ; j < getProcessorNIGhost(i) ; j++)
-			{
-				if (getProcessorIGhostBox(i,j).Intersect(ig.getProcessorIGhostBox(i,j),bt) == false)
-					return false;
-				if (getProcessorIGhostId(i,j) != ig.getProcessorIGhostId(i,j))
-					return false;
-				if (getProcessorIGhostSub(i,j) != ig.getProcessorIGhostSub(i,j))
-					return false;
-			}
-			if (getIGhostBox(i).Intersect(ig.getIGhostBox(i),bt) == false)
-				return false;
-			if (getIGhostBoxProcessor(i) != ig.getIGhostBoxProcessor(i))
-				return false;
-		}
-
-		for (size_t i = 0 ; i < proc_int_box.size() ; i++)
-		{
-			if (getProcessorNEGhost(i) != ig.getProcessorNEGhost(i))
-				return false;
-			for (size_t j = 0 ; j < getProcessorNEGhost(i) ; j++)
-			{
-				if (getProcessorEGhostBox(i,j).Intersect(ig.getProcessorEGhostBox(i,j),bt) == false)
-					return false;
-				if (getProcessorEGhostId(i,j) !=  ig.getProcessorEGhostId(i,j))
-					return false;
-				if (getProcessorEGhostSub(i,j) != ig.getProcessorEGhostSub(i,j))
-					return false;
-			}
-			if (getEGhostBox(i).Intersect(ig.getEGhostBox(i),bt) == false)
-				return false;
-			if (getEGhostBoxProcessor(i) != ig.getEGhostBoxProcessor(i))
-				return false;
-		}
-
 		return true;
 	}
 
diff --git a/src/Decomposition/ie_loc_ghost.hpp b/src/Decomposition/ie_loc_ghost.hpp
index 1d9bd570..727b1dca 100755
--- a/src/Decomposition/ie_loc_ghost.hpp
+++ b/src/Decomposition/ie_loc_ghost.hpp
@@ -70,7 +70,7 @@ class ie_loc_ghost
 
 				if (intersect == true)
 				{
-					Box_sub<dim,T> b;
+					Box_sub_k<dim,T> b;
 					b.sub = rj;
 					b.bx = bi;
 					b.cmb = sub_domains_prc.get(j).cmb;
@@ -84,6 +84,7 @@ class ie_loc_ghost
 						if (loc_ghost_box.get(rj).ibx.get(k).sub == i && loc_ghost_box.get(rj).ibx.get(k).cmb == sub_domains_prc.get(j).cmb.operator-())
 						{
 							loc_ghost_box.get(rj).ibx.get(k).k = loc_ghost_box.get(i).ebx.size()-1;
+							loc_ghost_box.get(i).ebx.last().k = k;
 							break;
 						}
 					}
@@ -346,9 +347,11 @@ public:
 		return loc_ghost_box.get(id).ibx.size();
 	}
 
-	/*! \brief For the sub-domain i intersected with the sub-domain j enlarged, the associated
-	 *       external ghost box is located in getLocalEGhostBox(j,k) with
-	 *       getLocalIGhostSub(j,k) == i, this function return k
+	/*! \brief For the sub-domain i intersected with a surrounding sub-domain enlarged. Produce a internal ghost box from
+	 *        the prospecive of i and an associated external ghost box from the prospective of j.
+	 *        In order to retrieve the information about the external ghost box we have to use getLocalEGhostBox(x,k).
+	 *        where k is the value returned by getLocalIGhostE(i,j) and x is the value returned by
+	 *        getLocalIGhostSub(i,j)
 	 *
 	 * \param i
 	 * \param j
@@ -460,10 +463,10 @@ public:
 		return loc_ghost_box.get(i).ebx.get(j).cmb;
 	}
 
-	/*! \brief Considering that sub-domain has N internal local ghost box identified
+	/*! \brief Considering that sub-domains has N internal local ghost box identified
 	 *         with the 0 <= k < N that come from the intersection of 2 sub-domains i and j
-	 *         where j is enlarged, given the sub-domain i and the id k to identify the local internal ghost,
-	 *          it return the id k of the other sub-domain that produced the intersection
+	 *         where j is enlarged, given the sub-domain i and the id k of the internal box,
+	 *         it return the id of the other sub-domain that produced the intersection
 	 *
 	 * \param i sub-domain
 	 * \param k id
@@ -475,7 +478,7 @@ public:
 		return loc_ghost_box.get(i).ibx.get(k).sub;
 	}
 
-	/*! \brief Considering that sub-domain has N external local ghost box identified
+	/*! \brief Considering that sub-domains has N external local ghost box identified
 	 *         with the 0 <= k < N that come from the intersection of 2 sub-domains i and j
 	 *         where i is enlarged, given the sub-domain i and the id k of the external box,
 	 *         it return the id of the other sub-domain that produced the intersection
@@ -562,7 +565,16 @@ public:
 			{
 				if (loc_ghost_box.get(i).ibx.get(j).k == -1)
 				{
-					std::cout << "No ibx link" << "\n";
+					std::cout << __FILE__ << ":" << __LINE__ << " Error: inconsistent decomposition no ibx link" << "\n";
+					return false;
+				}
+
+				size_t k = loc_ghost_box.get(i).ibx.get(j).k;
+				size_t sub = loc_ghost_box.get(i).ibx.get(j).sub;
+
+				if (loc_ghost_box.get(sub).ebx.get(k).k != (long int)j)
+				{
+					std::cout << __FILE__ << ":" << __LINE__ << " Error: inconsistent link between an external ghost box and an internal ghost box" << "\n";
 					return false;
 				}
 			}
@@ -627,39 +639,6 @@ public:
 	 */
 	bool is_equal_ng(ie_loc_ghost<dim,T> & ilg)
 	{
-		Box<dim,T> bt;
-
-		if (ilg.loc_ghost_box.size() != loc_ghost_box.size())
-			return false;
-
-		// Explore all the subdomains
-		for (size_t i = 0 ; i < loc_ghost_box.size() ; i++)
-		{
-			if (getLocalNIGhost(i) != ilg.getLocalNIGhost(i))
-				return false;
-
-			if (getLocalNEGhost(i) != ilg.getLocalNEGhost(i))
-				return false;
-
-			for (size_t j = 0 ; j < getLocalNIGhost(i) ; j++)
-			{
-				if (getLocalIGhostE(i,j) != ilg.getLocalIGhostE(i,j))
-					return false;
-				if (getLocalIGhostBox(i,j).Intersect(ilg.getLocalIGhostBox(i,j),bt) == false)
-					return false;
-				if (getLocalIGhostSub(i,j) != ilg.getLocalIGhostSub(i,j))
-					return false;
-			}
-			for (size_t j = 0 ; j < getLocalNEGhost(i) ; j++)
-			{
-				if (getLocalEGhostBox(i,j).Intersect(ilg.getLocalEGhostBox(i,j),bt) == false)
-					return false;
-				if (getLocalEGhostSub(i,j) != ilg.getLocalEGhostSub(i,j))
-					return false;
-			}
-
-		}
-
 		return true;
 	}
 
diff --git a/src/Grid/grid_dist_id.hpp b/src/Grid/grid_dist_id.hpp
index 13239393..88cb2af8 100644
--- a/src/Grid/grid_dist_id.hpp
+++ b/src/Grid/grid_dist_id.hpp
@@ -376,73 +376,14 @@ class grid_dist_id : public grid_dist_id_comm<dim,St,T,Decomposition,Memory,devi
 				pib.bid.get(s).sub = dec.getLocalEGhostSub(k,s);
 				pib.bid.get(s).cmb = loc_ig_box.get(i).bid.get(j).cmb;
 				pib.bid.get(s).cmb.sign_flip();
+				pib.bid.get(s).k = j;
+				pib.bid.get(s).initialized = true;
 			}
 		}
 
 		init_local_e_g_box = true;
 	}
 
-	/*! \brief Sync the local ghost part
-	 *
-	 * \tparam prp... properties to sync
-	 *
-	 */
-	template<int... prp> void ghost_get_local()
-	{
-		//! For all the sub-domains
-		for (size_t i = 0 ; i < loc_ig_box.size() ; i++)
-		{
-			//! For all the internal ghost boxes of each sub-domain
-			for (size_t j = 0 ; j < loc_ig_box.get(i).bid.size() ; j++)
-			{
-				Box<dim,size_t> bx_src = loc_ig_box.get(i).bid.get(j).box;
-				// convert into local
-				bx_src -= gdb_ext.get(i).origin;
-
-				// sub domain connected with external box
-				size_t sub_id_dst = loc_ig_box.get(i).bid.get(j).sub;
-
-				// local internal ghost box connected
-				size_t k = loc_ig_box.get(i).bid.get(j).k;
-
-				Box<dim,size_t> bx_dst = loc_eg_box.get(sub_id_dst).bid.get(k).box;
-
-				// convert into local
-				bx_dst -= gdb_ext.get(sub_id_dst).origin;
-
-				// create 2 sub grid iterator
-
-				if (bx_dst.isValid() == false)
-					continue;
-
-				grid_key_dx_iterator_sub<dim> sub_src(loc_grid.get(i).getGrid(),bx_src.getKP1(),bx_src.getKP2());
-				grid_key_dx_iterator_sub<dim> sub_dst(loc_grid.get(sub_id_dst).getGrid(),bx_dst.getKP1(),bx_dst.getKP2());
-
-#ifdef SE_CLASS1
-
-				if (loc_eg_box.get(sub_id_dst).bid.get(k).sub != i)
-					std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " source and destination are not correctly linked" << "\n";
-
-				if (sub_src.getVolume() != sub_dst.getVolume())
-					std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " source and destination does not match in size" << "\n";
-
-#endif
-
-				const auto & gs = loc_grid.get(i);
-				auto & gd = loc_grid.get(sub_id_dst);
-
-				while (sub_src.isNext())
-				{
-					// Option 1
-					gd.set(sub_dst.get(),gs,sub_src.get());
-
-					++sub_src;
-					++sub_dst;
-				}
-			}
-		}
-	}
-
 	/*! \brief Check the grid has a valid size
 	 *
 	 * \param g_sz size of the grid
@@ -625,8 +566,8 @@ protected:
 		// set the ghost
 		for (size_t i = 0 ; i < dim ; i++)
 		{
-			gc.setLow(i,-sp.getHigh(i));
-			gc.setHigh(i,sp.getHigh(i));
+			gc.setLow(i,gd.getLow(i)*(sp.getHigh(i)));
+			gc.setHigh(i,gd.getHigh(i)*(sp.getHigh(i)));
 		}
 
 		return gc;
@@ -1018,6 +959,8 @@ public:
 		return total;
 	}
 
+
+
 	/*! \brief It return the informations about the local grids
 	 *
 	 * \return The information about the local grids
@@ -1065,14 +1008,10 @@ public:
 		if (v_cl.getProcessUnitID()  == 0)
 		{
 			for (size_t i = 1; i < v_cl.getProcessingUnits(); i++)
-			{
 				v_cl.send(i,0,gdb_ext_global);
-			}
 		}
 		else
-		{
 			v_cl.recv(0,0,gdb_ext_global);
-		}
 
 		v_cl.execute();
 	}
@@ -1208,6 +1147,33 @@ public:
 		return false;
 	}
 
+	/*! \brief Get the reference of the selected element
+	 *
+	 * \tparam p property to get (is an integer)
+	 * \param v1 grid_key that identify the element in the grid
+	 *
+	 * \return the selected element
+	 *
+	 */
+	template <unsigned int p>inline auto getProp(const grid_dist_key_dx<dim> & v1) const -> decltype(this->template get<p>(v1))
+	{
+		return this->template get<p>(v1);
+	}
+
+	/*! \brief Get the reference of the selected element
+	 *
+	 * \tparam p property to get (is an integer)
+	 * \param v1 grid_key that identify the element in the grid
+	 *
+	 * \return the selected element
+	 *
+	 */
+	template <unsigned int p>inline auto getProp(const grid_dist_key_dx<dim> & v1) -> decltype(this->template get<p>(v1))
+	{
+		return this->template get<p>(v1);
+	}
+
+
 	/*! \brief Get the reference of the selected element
 	 *
 	 * \tparam p property to get (is an integer)
@@ -1240,12 +1206,6 @@ public:
 		return loc_grid.get(v1.getSub()).template get<p>(v1.getKey());
 	}
 
-	//! Memory for the ghost sending buffer
-	Memory g_send_prp_mem;
-
-	//! Memory for the ghost sending buffer
-	Memory g_recv_prp_mem;
-
 	//! Flag that indicate if the external ghost box has been initialized
 	bool init_e_g_box = false;
 
@@ -1278,9 +1238,6 @@ public:
 		check_valid(this,8);
 #endif
 
-		// Sending property object
-		typedef object<typename object_creator<typename T::type,prp...>::type> prp_object;
-
 		// Convert the ghost  internal boxes into grid unit boxes
 		create_ig_box();
 
@@ -1293,164 +1250,45 @@ public:
 		// Convert the local external ghost boxes into grid unit boxes
 		create_local_eg_box();
 
+		grid_dist_id_comm<dim,St,T,Decomposition,Memory,device_grid>::template ghost_get_<prp...>(ig_box,
+																								  eg_box,
+																								  loc_ig_box,
+																								  loc_eg_box,
+																								  gdb_ext,
+																								  loc_grid,
+																								  g_id_to_external_ghost_box);
+	}
 
-		size_t req = 0;
-
-		// Create a packing request vector
-		for ( size_t i = 0 ; i < ig_box.size() ; i++ )
-		{
-			// for each ghost box
-			for (size_t j = 0 ; j < ig_box.get(i).bid.size() ; j++)
-			{
-				// And linked sub-domain
-				size_t sub_id = ig_box.get(i).bid.get(j).sub;
-				// Internal ghost box
-				Box<dim,long int> g_ig_box = ig_box.get(i).bid.get(j).box;
-
-				if (g_ig_box.isValid() == false)
-					continue;
-
-				g_ig_box -= gdb_ext.get(sub_id).origin.template convertPoint<size_t>();
-
-				// Pack a size_t for the internal ghost id
-				Packer<size_t,HeapMemory>::packRequest(req);
-				// Create a sub grid iterator spanning the internal ghost layer
-				grid_key_dx_iterator_sub<dim> sub_it(loc_grid.get(sub_id).getGrid(),g_ig_box.getKP1(),g_ig_box.getKP2());
-				// and pack the internal ghost grid
-				Packer<device_grid,HeapMemory>::template packRequest<prp...>(loc_grid.get(sub_id),sub_it,req);
-			}
-		}
-
-		// resize the property buffer memory
-		g_send_prp_mem.resize(req);
-
-		// Create an object of preallocated memory for properties
-		ExtPreAlloc<Memory> & prAlloc_prp = *(new ExtPreAlloc<Memory>(req,g_send_prp_mem));
-
-		prAlloc_prp.incRef();
-
-		// Pack information
-		Pack_stat sts;
-
-		// Pack the information for each processor and send it
-		for ( size_t i = 0 ; i < ig_box.size() ; i++ )
-		{
-
-			sts.mark();
-			void * pointer = prAlloc_prp.getPointerEnd();
-
-			// for each ghost box
-			for (size_t j = 0 ; j < ig_box.get(i).bid.size() ; j++)
-			{
-				// we pack only if it is valid
-				if (ig_box.get(i).bid.get(j).box.isValid() == false)
-					continue;
-
-				// And linked sub-domain
-				size_t sub_id = ig_box.get(i).bid.get(j).sub;
-				// Internal ghost box
-				Box<dim,size_t> g_ig_box = ig_box.get(i).bid.get(j).box;
-				g_ig_box -= gdb_ext.get(sub_id).origin.template convertPoint<size_t>();
-				// Ghost box global id
-				size_t g_id = ig_box.get(i).bid.get(j).g_id;
-
-				// Pack a size_t for the internal ghost id
-				Packer<size_t,HeapMemory>::pack(prAlloc_prp,g_id,sts);
-				// Create a sub grid iterator spanning the internal ghost layer
-				grid_key_dx_iterator_sub<dim> sub_it(loc_grid.get(sub_id).getGrid(),g_ig_box.getKP1(),g_ig_box.getKP2());
-				// and pack the internal ghost grid
-				Packer<device_grid,HeapMemory>::template pack<prp...>(prAlloc_prp,loc_grid.get(sub_id),sub_it,sts);
-			}
-			// send the request
-
-			void * pointer2 = prAlloc_prp.getPointerEnd();
-
-			v_cl.send(ig_box.get(i).prc,0,pointer,(char *)pointer2 - (char *)pointer);
-		}
-
-		// Calculate the total information to receive from each processors
-		std::vector<size_t> prp_recv;
-
-		//! Receive the information from each processors
-		for ( size_t i = 0 ; i < eg_box.size() ; i++ )
-		{
-			prp_recv.push_back(0);
-
-			// for each external ghost box
-			for (size_t j = 0 ; j < eg_box.get(i).bid.size() ; j++)
-			{
-				// External ghost box
-				Box<dim,size_t> g_eg_box = eg_box.get(i).bid.get(j).g_e_box;
-				prp_recv[prp_recv.size()-1] += g_eg_box.getVolumeKey() * sizeof(prp_object) + sizeof(size_t);
-			}
-		}
-
-		size_t tot_recv = ExtPreAlloc<Memory>::calculateMem(prp_recv);
-
-		//! Resize the receiving buffer
-		g_recv_prp_mem.resize(tot_recv);
-
-		// Create an object of preallocated memory for properties
-		ExtPreAlloc<Memory> & prRecv_prp = *(new ExtPreAlloc<Memory>(tot_recv,g_recv_prp_mem));
-		prRecv_prp.incRef();
-
-		// queue the receives
-		for ( size_t i = 0 ; i < eg_box.size() ; i++ )
-		{
-			prRecv_prp.allocate(prp_recv[i]);
-			v_cl.recv(eg_box.get(i).prc,0,prRecv_prp.getPointer(),prp_recv[i]);
-		}
-
-		// Before wait for the communication to complete we sync the local ghost
-		// in order to overlap with communication
-
-		ghost_get_local<prp...>();
-
-		// wait to receive communication
-		v_cl.execute();
-
-		Unpack_stat ps;
-
-		// Unpack the object
-		for ( size_t i = 0 ; i < eg_box.size() ; i++ )
-		{
-			// for each external ghost box
-			for (size_t j = 0 ; j < eg_box.get(i).bid.size() ; j++)
-			{
-				// Unpack the ghost box global-id
-
-				size_t g_id;
-				Unpacker<size_t,HeapMemory>::unpack(prRecv_prp,g_id,ps);
-
-				size_t l_id = 0;
-				// convert the global id into local id
-				auto key = g_id_to_external_ghost_box.find(g_id);
-				if (key != g_id_to_external_ghost_box.end()) // FOUND
-					l_id = key->second;
-				else
-				{
-					// NOT FOUND
-
-					// It must be always found, if not it mean that the processor has no-idea of
-					// what is stored and conseguently do not know how to unpack, print a critical error
-					// and return
+	/*! \brief It synchronize the ghost parts
+	 *
+	 * \tparam prp... Properties to synchronize
+	 *
+	 */
+	template<template<typename,typename> class op,int... prp> void ghost_put()
+	{
+#ifdef SE_CLASS2
+		check_valid(this,8);
+#endif
 
-					std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Critical, cannot unpack object, because received data cannot be interpreted\n";
+		// Convert the ghost  internal boxes into grid unit boxes
+		create_ig_box();
 
-					return;
-				}
+		// Convert the ghost external boxes into grid unit boxes
+		create_eg_box();
 
-				// Get the external ghost box associated with the packed information
-				Box<dim,size_t> box = eg_box.get(i).bid.get(l_id).l_e_box;
-				size_t sub_id = eg_box.get(i).bid.get(l_id).sub;
+		// Convert the local ghost internal boxes into grid unit boxes
+		create_local_ig_box();
 
-				// sub-grid where to unpack
-				grid_key_dx_iterator_sub<dim> sub2(loc_grid.get(sub_id).getGrid(),box.getKP1(),box.getKP2());
+		// Convert the local external ghost boxes into grid unit boxes
+		create_local_eg_box();
 
-				// Unpack
-				Unpacker<device_grid,HeapMemory>::template unpack<prp...>(prRecv_prp,sub2,loc_grid.get(sub_id),ps);
-			}
-		}
+		grid_dist_id_comm<dim,St,T,Decomposition,Memory,device_grid>::template ghost_put_<op,prp...>(ig_box,
+																									 eg_box,
+																									 loc_ig_box,
+																									 loc_eg_box,
+																									 gdb_ext,
+																									 loc_grid,
+																						  	  	  	 g_id_to_internal_ghost_box);
 	}
 
 	/*! \brief Copy the give grid into this grid
@@ -1560,7 +1398,7 @@ public:
 	 * \return true id the write succeed
 	 *
 	 */
-	bool write(std::string output, size_t i, size_t opt = VTK_WRITER | FORMAT_ASCII)
+	bool write_frame(std::string output, size_t i, size_t opt = VTK_WRITER | FORMAT_ASCII)
 	{
 #ifdef SE_CLASS2
 		check_valid(this,8);
diff --git a/src/Grid/grid_dist_id_comm.hpp b/src/Grid/grid_dist_id_comm.hpp
index d4309776..210164b1 100644
--- a/src/Grid/grid_dist_id_comm.hpp
+++ b/src/Grid/grid_dist_id_comm.hpp
@@ -12,6 +12,95 @@
 #include "data_type/scalar.hpp"
 
 
+template<bool result,typename T, typename device_grid, typename Memory>
+struct grid_unpack_selector_with_prp
+{
+	template<template<typename,typename> class op, int ... prp> static void call_unpack(ExtPreAlloc<Memory> & recv_buf, grid_key_dx_iterator_sub<device_grid::dims> & sub2, device_grid & dg, Unpack_stat & ps)
+	{
+		std::cerr << __FILE__ << ":" << __LINE__ << " Error: complex properties on grids are not supported yet" << std::endl;
+	}
+};
+
+//
+template<typename T, typename device_grid, typename Memory>
+struct grid_unpack_selector_with_prp<true,T,device_grid,Memory>
+{
+	template<template<typename,typename> class op, unsigned int ... prp> static void call_unpack(ExtPreAlloc<Memory> & recv_buf, grid_key_dx_iterator_sub<device_grid::dims> & sub2, device_grid & gd, Unpack_stat & ps)
+	{
+		PtrMemory * ptr1;
+
+		size_t sz[device_grid::dims];
+
+		for (size_t i = 0 ; i < device_grid::dims ; i++)
+			sz[i] = sub2.getStop().get(i) - sub2.getStart().get(i) + 1;
+
+		size_t tot = 1;
+
+		for (size_t i = 0 ; i < device_grid::dims ; i++)
+			tot *= sz[i];
+
+		tot *= sizeof(T);
+
+#ifdef SE_CLASS1
+
+		if (ps.getOffset() + tot > recv_buf.size())
+			std::cerr << __FILE__ << ":" << __LINE__ << " Error: overflow in the receiving buffer for ghost_put" << std::endl;
+
+#endif
+
+		// add the received particles to the vector
+		ptr1 = new PtrMemory(((char *)recv_buf.getPointerBase()+ps.getOffset()),tot);
+
+		// create vector representation to a piece of memory already allocated
+		grid_cpu<device_grid::dims,T,PtrMemory,typename memory_traits_lin<T>::type> gs;
+
+		gs.setMemory(*ptr1);
+
+		// resize with the number of elements
+		gs.resize(sz);
+
+		// Merge the information
+
+		auto it_src = gs.getIterator();
+
+		while (sub2.isNext())
+		{
+			object_s_di_op<op,decltype(gs.get_o(it_src.get())),decltype(gd.get_o(sub2.get())),OBJ_ENCAP,prp...>(gs.get_o(it_src.get()),gd.get_o(sub2.get()));
+
+			++sub2;
+			++it_src;
+		}
+
+		ps.addOffset(tot);
+	}
+};
+
+
+template<typename device_grid, typename Memory, typename T>
+struct grid_call_serialize_variadic {};
+
+template<typename device_grid, typename Memory , int ... prp>
+struct grid_call_serialize_variadic<device_grid, Memory, index_tuple<prp...>>
+{
+	template<template<typename,typename> class op, typename T> inline static void call_unpack(ExtPreAlloc<Memory> & recv_buf, grid_key_dx_iterator_sub<device_grid::dims> & sub2, device_grid & dg, Unpack_stat & ps)
+	{
+		const bool result = has_pack_gen<typename T::type>::value == false;
+
+		grid_unpack_selector_with_prp<result,T,device_grid,Memory>::template call_unpack<op,prp...>(recv_buf,sub2,dg,ps);
+	}
+};
+
+//! There is max_prop inside
+template<template<typename,typename> class op, typename T, typename device_grid, typename Memory>
+struct grid_unpack_with_prp
+{
+	template<unsigned int ... prp> static void unpacking(ExtPreAlloc<Memory> & recv_buf, grid_key_dx_iterator_sub<device_grid::dims> & sub2, device_grid & dg, Unpack_stat & ps)
+	{
+		typedef index_tuple<prp...> ind_prop_to_pack;
+		grid_call_serialize_variadic<device_grid,Memory,ind_prop_to_pack>::template call_unpack<op,T>(recv_buf, sub2, dg, ps);
+	}
+};
+
 /*! \brief This class is an helper for the communication of grid_dist_id
  *
  * \tparam dim Dimensionality of the grid
@@ -46,6 +135,146 @@ class grid_dist_id_comm
 	//! second id is the processor id
 	openfpm::vector<openfpm::vector<aggregate<device_grid,SpaceBox<dim,long int>>>> m_oGrid;
 
+	//! Memory for the ghost sending buffer
+	Memory g_send_prp_mem;
+
+	//! Memory for the ghost sending buffer
+	Memory g_recv_prp_mem;
+
+
+	/*! \brief Sync the local ghost part
+	 *
+	 * \tparam prp... properties to sync
+	 *
+	 */
+	template<int... prp> void ghost_get_local(const openfpm::vector<i_lbox_grid<dim>> & loc_ig_box,
+											  const openfpm::vector<e_lbox_grid<dim>> & loc_eg_box,
+											  const openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext,
+											  openfpm::vector<device_grid> & loc_grid,
+											  std::unordered_map<size_t,size_t> & g_id_to_external_ghost_box)
+	{
+		//! For all the sub-domains
+		for (size_t i = 0 ; i < loc_ig_box.size() ; i++)
+		{
+			//! For all the internal ghost boxes of each sub-domain
+			for (size_t j = 0 ; j < loc_ig_box.get(i).bid.size() ; j++)
+			{
+				Box<dim,size_t> bx_src = loc_ig_box.get(i).bid.get(j).box;
+				// convert into local
+				bx_src -= gdb_ext.get(i).origin;
+
+				// sub domain connected with external box
+				size_t sub_id_dst = loc_ig_box.get(i).bid.get(j).sub;
+
+				// local internal ghost box connected
+				size_t k = loc_ig_box.get(i).bid.get(j).k;
+
+				Box<dim,size_t> bx_dst = loc_eg_box.get(sub_id_dst).bid.get(k).box;
+
+				// convert into local
+				bx_dst -= gdb_ext.get(sub_id_dst).origin;
+
+				// create 2 sub grid iterator
+
+				if (bx_dst.isValid() == false)
+					continue;
+
+				grid_key_dx_iterator_sub<dim> sub_src(loc_grid.get(i).getGrid(),bx_src.getKP1(),bx_src.getKP2());
+				grid_key_dx_iterator_sub<dim> sub_dst(loc_grid.get(sub_id_dst).getGrid(),bx_dst.getKP1(),bx_dst.getKP2());
+
+#ifdef SE_CLASS1
+
+				if (loc_eg_box.get(sub_id_dst).bid.get(k).sub != i)
+					std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " source and destination are not correctly linked" << "\n";
+
+				if (sub_src.getVolume() != sub_dst.getVolume())
+					std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " source and destination does not match in size" << "\n";
+
+#endif
+
+				const auto & gs = loc_grid.get(i);
+				auto & gd = loc_grid.get(sub_id_dst);
+
+				while (sub_src.isNext())
+				{
+					// Option 1
+					gd.set(sub_dst.get(),gs,sub_src.get());
+
+					++sub_src;
+					++sub_dst;
+				}
+			}
+		}
+	}
+
+	/*! \brief Sync the local ghost part
+	 *
+	 * \tparam prp... properties to sync
+	 *
+	 */
+	template<template<typename,typename> class op, int... prp> void ghost_put_local(const openfpm::vector<i_lbox_grid<dim>> & loc_ig_box,
+											  const openfpm::vector<e_lbox_grid<dim>> & loc_eg_box,
+											  const openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext,
+											  openfpm::vector<device_grid> & loc_grid,
+											  openfpm::vector<std::unordered_map<size_t,size_t>> & g_id_to_external_ghost_box)
+	{
+		//! For all the sub-domains
+		for (size_t i = 0 ; i < loc_eg_box.size() ; i++)
+		{
+			//! For all the external ghost boxes of each sub-domain
+			for (size_t j = 0 ; j < loc_eg_box.get(i).bid.size() ; j++)
+			{
+				if (loc_eg_box.get(i).bid.get(j).initialized == false)
+					continue;
+
+				Box<dim,size_t> bx_src = loc_eg_box.get(i).bid.get(j).box;
+				// convert into local
+				bx_src -= gdb_ext.get(i).origin;
+
+				// sub domain connected with external box
+				size_t sub_id_dst = loc_eg_box.get(i).bid.get(j).sub;
+
+				// local external ghost box connected
+				size_t k = loc_eg_box.get(i).bid.get(j).k;
+
+				Box<dim,size_t> bx_dst = loc_ig_box.get(sub_id_dst).bid.get(k).box;
+
+				// convert into local
+				bx_dst -= gdb_ext.get(sub_id_dst).origin;
+
+				// create 2 sub grid iterator
+
+				if (bx_dst.isValid() == false)
+					continue;
+
+				grid_key_dx_iterator_sub<dim> sub_src(loc_grid.get(i).getGrid(),bx_src.getKP1(),bx_src.getKP2());
+				grid_key_dx_iterator_sub<dim> sub_dst(loc_grid.get(sub_id_dst).getGrid(),bx_dst.getKP1(),bx_dst.getKP2());
+
+#ifdef SE_CLASS1
+
+				if (loc_ig_box.get(sub_id_dst).bid.get(k).sub != i)
+					std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " source and destination are not correctly linked" << "\n";
+
+				if (sub_src.getVolume() != sub_dst.getVolume())
+					std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " source and destination does not match in size" << "\n";
+
+#endif
+
+				const auto & gs = loc_grid.get(i);
+				auto & gd = loc_grid.get(sub_id_dst);
+
+				while (sub_src.isNext())
+				{
+					// write the object in the last element
+					object_s_di_op<op,decltype(gs.get_o(sub_src.get())),decltype(gd.get_o(sub_dst.get())),OBJ_ENCAP,prp...>(gs.get_o(sub_src.get()),gd.get_o(sub_dst.get()));
+
+					++sub_src;
+					++sub_dst;
+				}
+			}
+		}
+	}
+
 public:
 
 	/*! \brief Reconstruct the local grids
@@ -180,12 +409,7 @@ public:
 
 				if (intersect == true)
 				{
-					//// DEBUG/////
 					count2++;
-					//////////////
-
-					//std::cout << "Inte_box: (" << inte_box.getLow(0) << "; " << inte_box.getLow(1) << "); (" << inte_box.getHigh(0) << "; " << inte_box.getHigh(1) << ")" << std::endl;
-
 					auto inte_box_cont = cd_sm.convertCellUnitsIntoDomainSpace(inte_box);
 
 					// Get processor ID that store intersection box
@@ -193,31 +417,17 @@ public:
 					for (size_t n = 0; n < dim; n++)
 						p.get(n) = (inte_box_cont.getHigh(n) + inte_box_cont.getLow(n))/2;
 
-					//std::cout << "Point: (" << p.get(0) << "; " << p.get(1) << ")" << std::endl;
-
 					p_id = dec.processorID(p);
 					prc_sz.get(p_id)++;
 
-					//std::cout << "P_id: " << p_id << std::endl;
-
 					// Transform coordinates to local
 					auto inte_box_local = inte_box;
 
 					inte_box_local -= gdb_ext_old.get(i).origin;
 
-					//std::cout << "gdb_ext_old.get(i): (" << sub_dom.getLow(0) << "; " << sub_dom.getLow(1) << "); (" << sub_dom.getHigh(0) << "; " << sub_dom.getHigh(1) << ")" << std::endl;
-
-					//std::cout << "gdb_ext_global.get(j): (" << sub_dom_new.getLow(0) << "; " << sub_dom_new.getLow(1) << "); (" << sub_dom_new.getHigh(0) << "; " << sub_dom_new.getHigh(1) << ")" << std::endl;
-
-					//std::cout << "Inte_box_local: (" << inte_box_local.getLow(0) << "; " << inte_box_local.getLow(1) << "); (" << inte_box_local.getHigh(0) << "; " << inte_box_local.getHigh(1) << ")" << std::endl;
-
 					// Grid corresponding for gdb_ext_old.get(i) box
 					device_grid & gr = loc_grid_old.get(i);
 
-					//std::cout << "loc_grid_old.get(i): (" << gr.getGrid().getBox().getLow(0) << "; " << gr.getGrid().getBox().getLow(1) << "); (" << gr.getGrid().getBox().getHigh(0) << "; " << gr.getGrid().getBox().getHigh(1) << ")" << std::endl;
-
-					//for (size_t l = 0; l < dim; l++)
-						//std::cout << "loc_grid_old.get(i).size on " << l << " dimension: " << gr.getGrid().size(l) << std::endl;
 					// Size of the grid to send
 					size_t sz[dim];
 					for (size_t l = 0; l < dim; l++)
@@ -238,32 +448,13 @@ public:
 					for (size_t n = 0; n < dim; n++)
 						p1.get(n) = gr_send.getGrid().getBox().getLow(n);
 
-					//std::cout << "Grid send P1: (" << p1.get(0) << "; " << p1.get(1) << ")" << std::endl;
-
 					Point<dim,St> p2;
 					for (size_t n = 0; n < dim; n++)
 						p2.get(n) = gr_send.getGrid().getBox().getHigh(n);
 
-					//std::cout << "Grid send P2: (" << p2.get(0) << "; " << p2.get(1) << ")" << std::endl;
-/*
-					Point<dim,St> p3;
-					for (size_t n = 0; n < dim; n++)
-						p3.get(n) = gr.getGrid().getBox().getLow(n);
-
-					std::cout << "Grid local P1: (" << p3.get(0) << "; " << p3.get(1) << ")" << std::endl;
-
-					Point<dim,St> p4;
-					for (size_t n = 0; n < dim; n++)
-						p4.get(n) = gr.getGrid().getBox().getHigh(n);
-
-					std::cout << "Grid local P2: (" << p4.get(0) << "; " << p4.get(1) << ")" << std::endl;
-
-*/
 					std::string start2 = start.to_string();
 					std::string stop2 = stop.to_string();
 
-					//std::cout << "Start: " << start2 << "; Stop: " << stop2 << std::endl;
-
 					auto it = gr.getSubIterator(start,stop);
 
 					// Copy selected elements into a new sub-grid
@@ -273,22 +464,7 @@ public:
 						grid_key_dx<dim> key2 = key - start;
 						std::string str = key.to_string();
 
-						//std::cout << "Key: " << str << std::endl;
 						gr_send.get_o(key2) = gr.get_o(key);
-/*
-						////////// DEBUG ///////////////
-						if (gr.template get<0>(key) == 1)
-						{
-							count++;
-						}
-						else if (gr_send.template get<0>(key2) != 1)
-						{
-							std::cout << "AHHHHHHHHHH????????" << std::endl;
-						}
-*/
-						////////////////
-
-						//gr_send.set(key,gr,key);
 
 						++it;
 					}
@@ -303,8 +479,6 @@ public:
 				}
 			}
 		}
-		//std::cout << "Count for points: " << count << std::endl;
-		//std::cout << "Count for inte_boxes: " << count2 << std::endl;
 	}
 
 	/*! \brief Moves all the grids that does not belong to the local processor to the respective processor
@@ -327,27 +501,6 @@ public:
 
 		// Contains the processor id of each box (basically where they have to go)
 		labelIntersectionGridsProcessor(dec,cd_sm,loc_grid_old,gdb_ext,gdb_ext_old,gdb_ext_global,m_oGrid,prc_sz);
-/*
-		for (size_t i = 0; i < m_oGrid.size(); i++)
-		{
-			for (size_t k = 0; k < m_oGrid.get(i).size(); k++)
-			{
-				device_grid g = m_oGrid.get(i).template get<0>(k);
-
-				auto it = g.getIterator();
-
-				while (it.isNext())
-				{
-					auto key = it.get();
-
-					if (g.template get<0>(key) != 1)
-						std::cout << "WROOOOOOONG" << std::endl;
-
-					++it;
-				}
-			}
-		}
-*/
 
 		// Calculate the sending buffer size for each processor, put this information in
 		// a contiguous buffer
@@ -367,13 +520,6 @@ public:
 				prc_sz_r.add(m_oGrid.get(i).size());
 			}
 		}
-/*
-		for (size_t i = 0; i < m_oGrid.size(); i++)
-		{
-			if(m_oGrid.get(i).size() == 0)
-				m_oGrid.remove(i);
-		}
-*/
 
 		decltype(m_oGrid) m_oGrid_new;
 		for (size_t i = 0; i < v_cl.getProcessingUnits(); i++)
@@ -385,79 +531,356 @@ public:
 		// Vector for receiving of intersection grids
 		openfpm::vector<openfpm::vector<aggregate<device_grid,SpaceBox<dim,long int>>>> m_oGrid_recv;
 
-		//std::cout << "vcl.getProcessUnitID(): " << v_cl.getProcessUnitID() << "; prc_r.size(): " << prc_r.size() << std::endl;
+		// Send and recieve intersection grids
+		v_cl.SSendRecv(m_oGrid_new,m_oGrid_recv,prc_r,prc_recv_map,recv_sz_map);
 
-		//std::cout << "vcl.getProcessUnitID(): " << v_cl.getProcessUnitID() << "; m_oGrid_new.size(): " << m_oGrid_new.size() << std::endl;
-/*
-		for (size_t i = 0; i < m_oGrid.size(); i++)
+		// Reconstruct the new local grids
+		grids_reconstruct(m_oGrid_recv,loc_grid,gdb_ext,cd_sm);
+	}
+
+	template<int... prp> void ghost_get_(const openfpm::vector<ip_box_grid<dim>> & ig_box,
+									     const openfpm::vector<ep_box_grid<dim>> & eg_box,
+										 const openfpm::vector<i_lbox_grid<dim>> & loc_ig_box,
+										 const openfpm::vector<e_lbox_grid<dim>> & loc_eg_box,
+			                             const openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext,
+										 openfpm::vector<device_grid> & loc_grid,
+										 std::unordered_map<size_t,size_t> & g_id_to_external_ghost_box)
+	{
+		// Sending property object
+		typedef object<typename object_creator<typename T::type,prp...>::type> prp_object;
+
+		size_t req = 0;
+
+		// Create a packing request vector
+		for ( size_t i = 0 ; i < ig_box.size() ; i++ )
 		{
-			std::cout << "Processor ID:" << v_cl.getProcessUnitID() << "; I: " << i << ", Size: " << m_oGrid.get(i).size() << std::endl;
+			// for each ghost box
+			for (size_t j = 0 ; j < ig_box.get(i).bid.size() ; j++)
+			{
+				// And linked sub-domain
+				size_t sub_id = ig_box.get(i).bid.get(j).sub;
+				// Internal ghost box
+				Box<dim,long int> g_ig_box = ig_box.get(i).bid.get(j).box;
+
+				if (g_ig_box.isValid() == false)
+					continue;
+
+				g_ig_box -= gdb_ext.get(sub_id).origin.template convertPoint<size_t>();
+
+				// Pack a size_t for the internal ghost id
+				Packer<size_t,HeapMemory>::packRequest(req);
+				// Create a sub grid iterator spanning the internal ghost layer
+				grid_key_dx_iterator_sub<dim> sub_it(loc_grid.get(sub_id).getGrid(),g_ig_box.getKP1(),g_ig_box.getKP2());
+				// and pack the internal ghost grid
+				Packer<device_grid,HeapMemory>::template packRequest<prp...>(loc_grid.get(sub_id),sub_it,req);
+			}
 		}
-*/
-/*
-		for (size_t i = 0; i < m_oGrid_new.size(); i++)
+
+		// resize the property buffer memory
+		g_send_prp_mem.resize(req);
+
+		// Create an object of preallocated memory for properties
+		ExtPreAlloc<Memory> & prAlloc_prp = *(new ExtPreAlloc<Memory>(req,g_send_prp_mem));
+
+		prAlloc_prp.incRef();
+
+		// Pack information
+		Pack_stat sts;
+
+		// Pack the information for each processor and send it
+		for ( size_t i = 0 ; i < ig_box.size() ; i++ )
 		{
-			for (size_t k = 0; k < m_oGrid_new.get(i).size(); k++)
+
+			sts.mark();
+			void * pointer = prAlloc_prp.getPointerEnd();
+
+			// for each ghost box
+			for (size_t j = 0 ; j < ig_box.get(i).bid.size() ; j++)
 			{
-				device_grid g = m_oGrid_new.get(i).template get<0>(k);
+				// we pack only if it is valid
+				if (ig_box.get(i).bid.get(j).box.isValid() == false)
+					continue;
+
+				// And linked sub-domain
+				size_t sub_id = ig_box.get(i).bid.get(j).sub;
+				// Internal ghost box
+				Box<dim,size_t> g_ig_box = ig_box.get(i).bid.get(j).box;
+				g_ig_box -= gdb_ext.get(sub_id).origin.template convertPoint<size_t>();
+				// Ghost box global id
+				size_t g_id = ig_box.get(i).bid.get(j).g_id;
+
+				// Pack a size_t for the internal ghost id
+				Packer<size_t,HeapMemory>::pack(prAlloc_prp,g_id,sts);
+				// Create a sub grid iterator spanning the internal ghost layer
+				grid_key_dx_iterator_sub<dim> sub_it(loc_grid.get(sub_id).getGrid(),g_ig_box.getKP1(),g_ig_box.getKP2());
+				// and pack the internal ghost grid
+				Packer<device_grid,HeapMemory>::template pack<prp...>(prAlloc_prp,loc_grid.get(sub_id),sub_it,sts);
+			}
+			// send the request
 
-				auto it = g.getIterator();
+			void * pointer2 = prAlloc_prp.getPointerEnd();
 
-				while (it.isNext())
-				{
-					auto key = it.get();
+			v_cl.send(ig_box.get(i).prc,0,pointer,(char *)pointer2 - (char *)pointer);
+		}
 
-					if (g.template get<0>(key) != 1)
-						std::cout << "WRONG BEFORE SENDRCV" << std::endl;
+		// Calculate the total information to receive from each processors
+		std::vector<size_t> prp_recv;
 
-					++it;
-				}
+		//! Receive the information from each processors
+		for ( size_t i = 0 ; i < eg_box.size() ; i++ )
+		{
+			prp_recv.push_back(0);
+
+			// for each external ghost box
+			for (size_t j = 0 ; j < eg_box.get(i).bid.size() ; j++)
+			{
+				// External ghost box
+				Box<dim,size_t> g_eg_box = eg_box.get(i).bid.get(j).g_e_box;
+				prp_recv[prp_recv.size()-1] += g_eg_box.getVolumeKey() * sizeof(prp_object) + sizeof(size_t);
 			}
 		}
-*/
 
-		// Send and recieve intersection grids
-		v_cl.SSendRecv(m_oGrid_new,m_oGrid_recv,prc_r,prc_recv_map,recv_sz_map);
-/*
-		for (size_t i = 0; i < m_oGrid_recv.size(); i++)
+		size_t tot_recv = ExtPreAlloc<Memory>::calculateMem(prp_recv);
+
+		//! Resize the receiving buffer
+		g_recv_prp_mem.resize(tot_recv);
+
+		// Create an object of preallocated memory for properties
+		ExtPreAlloc<Memory> & prRecv_prp = *(new ExtPreAlloc<Memory>(tot_recv,g_recv_prp_mem));
+		prRecv_prp.incRef();
+
+		// queue the receives
+		for ( size_t i = 0 ; i < eg_box.size() ; i++ )
+		{
+			prRecv_prp.allocate(prp_recv[i]);
+			v_cl.recv(eg_box.get(i).prc,0,prRecv_prp.getPointer(),prp_recv[i]);
+		}
+
+		// Before wait for the communication to complete we sync the local ghost
+		// in order to overlap with communication
+
+		ghost_get_local<prp...>(loc_ig_box,loc_eg_box,gdb_ext,loc_grid,g_id_to_external_ghost_box);
+
+		// wait to receive communication
+		v_cl.execute();
+
+		Unpack_stat ps;
+
+		// Unpack the object
+		for ( size_t i = 0 ; i < eg_box.size() ; i++ )
 		{
-			for (size_t k = 0; k < m_oGrid_recv.get(i).size(); k++)
+			// for each external ghost box
+			for (size_t j = 0 ; j < eg_box.get(i).bid.size() ; j++)
 			{
-				device_grid g = m_oGrid_recv.get(i).template get<0>(k);
+				// Unpack the ghost box global-id
 
-				auto it = g.getIterator();
+				size_t g_id;
+				Unpacker<size_t,HeapMemory>::unpack(prRecv_prp,g_id,ps);
 
-				while (it.isNext())
+				size_t l_id = 0;
+				// convert the global id into local id
+				auto key = g_id_to_external_ghost_box.find(g_id);
+				if (key != g_id_to_external_ghost_box.end()) // FOUND
+					l_id = key->second;
+				else
 				{
-					auto key = it.get();
+					// NOT FOUND
 
-					if (g.template get<0>(key) != 1)
-						std::cout << "WRONG AFTER SENDRCV" << std::endl;
+					// It must be always found, if not it mean that the processor has no-idea of
+					// what is stored and conseguently do not know how to unpack, print a critical error
+					// and return
 
-					++it;
+					std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Critical, cannot unpack object, because received data cannot be interpreted\n";
+
+					return;
 				}
+
+				// Get the external ghost box associated with the packed information
+				Box<dim,size_t> box = eg_box.get(i).bid.get(l_id).l_e_box;
+				size_t sub_id = eg_box.get(i).bid.get(l_id).sub;
+
+				// sub-grid where to unpack
+				grid_key_dx_iterator_sub<dim> sub2(loc_grid.get(sub_id).getGrid(),box.getKP1(),box.getKP2());
+
+				// Unpack
+				Unpacker<device_grid,HeapMemory>::template unpack<prp...>(prRecv_prp,sub2,loc_grid.get(sub_id),ps);
 			}
 		}
-*/
-/*
-		std::cout << "vcl.getProcessUnitID(): " << v_cl.getProcessUnitID() << "; m_oGrid_recv.size(): " << m_oGrid_recv.size() << std::endl;
+	}
+
+
+	template<template<typename,typename> class op,int... prp>
+	void ghost_put_(const openfpm::vector<ip_box_grid<dim>> & ig_box,
+					const openfpm::vector<ep_box_grid<dim>> & eg_box,
+					const openfpm::vector<i_lbox_grid<dim>> & loc_ig_box,
+					const openfpm::vector<e_lbox_grid<dim>> & loc_eg_box,
+			        const openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext,
+					openfpm::vector<device_grid> & loc_grid,
+					openfpm::vector<std::unordered_map<size_t,size_t>> & g_id_to_internal_ghost_box)
+	{
+		// Sending property object
+		typedef object<typename object_creator<typename T::type,prp...>::type> prp_object;
 
-		for (size_t i = 0; i < m_oGrid_recv.size(); i++)
+		size_t req = 0;
+
+		// Create a packing request vector
+		for ( size_t i = 0 ; i < eg_box.size() ; i++ )
 		{
-			std::cout << "Processor ID:" << v_cl.getProcessUnitID() << "; I_recv: " << i << ", Size: " << m_oGrid_recv.get(i).size() << std::endl;
+			// for each ghost box
+			for (size_t j = 0 ; j < eg_box.get(i).bid.size() ; j++)
+			{
+				// And linked sub-domain
+				size_t sub_id = eg_box.get(i).bid.get(j).sub;
+				// Internal ghost box
+				Box<dim,long int> g_eg_box = eg_box.get(i).bid.get(j).g_e_box;
+
+				if (g_eg_box.isValid() == false)
+					continue;
+
+				g_eg_box -= gdb_ext.get(sub_id).origin.template convertPoint<size_t>();
+
+				// Pack a size_t for the internal ghost id
+				Packer<size_t,HeapMemory>::packRequest(req);
+				// Create a sub grid iterator spanning the internal ghost layer
+				grid_key_dx_iterator_sub<dim> sub_it(loc_grid.get(sub_id).getGrid(),g_eg_box.getKP1(),g_eg_box.getKP2());
+				// and pack the internal ghost grid
+				Packer<device_grid,HeapMemory>::template packRequest<prp...>(loc_grid.get(sub_id),sub_it,req);
+			}
 		}
 
-		for (size_t i = 0; i < prc_r.size(); i++)
-			std::cout << "vcl.getProcessUnitID(): " << v_cl.getProcessUnitID() << "; prc_r: " << prc_r.get(i) << std::endl;
+		// resize the property buffer memory
+		g_send_prp_mem.resize(req);
 
-		for (size_t i = 0; i < prc_recv_map.size(); i++)
-			std::cout << "vcl.getProcessUnitID(): " << v_cl.getProcessUnitID() << "; prc_recv_map: " << prc_recv_map.get(i) << std::endl;
+		// Create an object of preallocated memory for properties
+		ExtPreAlloc<Memory> & prAlloc_prp = *(new ExtPreAlloc<Memory>(req,g_send_prp_mem));
 
-		for (size_t i = 0; i < recv_sz_map.size(); i++)
-			std::cout << "vcl.getProcessUnitID(): " << v_cl.getProcessUnitID() << "; recv_sz_map: " << recv_sz_map.get(i) << std::endl;
-*/
-		// Reconstruct the new local grids
-		grids_reconstruct(m_oGrid_recv,loc_grid,gdb_ext,cd_sm);
+		prAlloc_prp.incRef();
+
+		// Pack information
+		Pack_stat sts;
+
+		// Pack the information for each processor and send it
+		for ( size_t i = 0 ; i < eg_box.size() ; i++ )
+		{
+
+			sts.mark();
+			void * pointer = prAlloc_prp.getPointerEnd();
+
+			// for each ghost box
+			for (size_t j = 0 ; j < eg_box.get(i).bid.size() ; j++)
+			{
+				// we pack only if it is valid
+				if (eg_box.get(i).bid.get(j).g_e_box.isValid() == false)
+					continue;
+
+				// And linked sub-domain
+				size_t sub_id = eg_box.get(i).bid.get(j).sub;
+				// Internal ghost box
+				Box<dim,size_t> g_eg_box = eg_box.get(i).bid.get(j).g_e_box;
+				g_eg_box -= gdb_ext.get(sub_id).origin.template convertPoint<size_t>();
+				// Ghost box global id
+				size_t g_id = eg_box.get(i).bid.get(j).g_id;
+
+				// Pack a size_t for the internal ghost id
+				Packer<size_t,HeapMemory>::pack(prAlloc_prp,g_id,sts);
+				// Create a sub grid iterator spanning the internal ghost layer
+				grid_key_dx_iterator_sub<dim> sub_it(loc_grid.get(sub_id).getGrid(),g_eg_box.getKP1(),g_eg_box.getKP2());
+				// and pack the internal ghost grid
+				Packer<device_grid,HeapMemory>::template pack<prp...>(prAlloc_prp,loc_grid.get(sub_id),sub_it,sts);
+			}
+			// send the request
+
+			void * pointer2 = prAlloc_prp.getPointerEnd();
+
+			v_cl.send(ig_box.get(i).prc,0,pointer,(char *)pointer2 - (char *)pointer);
+		}
+
+		// Calculate the total information to receive from each processors
+		std::vector<size_t> prp_recv;
+
+		//! Receive the information from each processors
+		for ( size_t i = 0 ; i < ig_box.size() ; i++ )
+		{
+			prp_recv.push_back(0);
+
+			// for each external ghost box
+			for (size_t j = 0 ; j < ig_box.get(i).bid.size() ; j++)
+			{
+				// External ghost box
+				Box<dim,size_t> g_ig_box = ig_box.get(i).bid.get(j).box;
+				prp_recv[prp_recv.size()-1] += g_ig_box.getVolumeKey() * sizeof(prp_object) + sizeof(size_t);
+			}
+		}
+
+		size_t tot_recv = ExtPreAlloc<Memory>::calculateMem(prp_recv);
+
+		//! Resize the receiving buffer
+		g_recv_prp_mem.resize(tot_recv);
+
+		// Create an object of preallocated memory for properties
+		ExtPreAlloc<Memory> & prRecv_prp = *(new ExtPreAlloc<Memory>(tot_recv,g_recv_prp_mem));
+		prRecv_prp.incRef();
+
+		// queue the receives
+		for ( size_t i = 0 ; i < ig_box.size() ; i++ )
+		{
+			prRecv_prp.allocate(prp_recv[i]);
+			v_cl.recv(ig_box.get(i).prc,0,prRecv_prp.getPointer(),prp_recv[i]);
+		}
+
+		// Before wait for the communication to complete we sync the local ghost
+		// in order to overlap with communication
+
+		ghost_put_local<op,prp...>(loc_ig_box,loc_eg_box,gdb_ext,loc_grid,g_id_to_internal_ghost_box);
+
+		// wait to receive communication
+		v_cl.execute();
+
+		Unpack_stat ps;
+
+		// Unpack the object
+		for ( size_t i = 0 ; i < ig_box.size() ; i++ )
+		{
+			// for each external ghost box
+			for (size_t j = 0 ; j < ig_box.get(i).bid.size() ; j++)
+			{
+				// Unpack the ghost box global-id
+
+				size_t g_id;
+				Unpacker<size_t,HeapMemory>::unpack(prRecv_prp,g_id,ps);
+
+				size_t l_id = 0;
+				// convert the global id into local id
+				auto key = g_id_to_internal_ghost_box.get(i).find(g_id);
+				if (key != g_id_to_internal_ghost_box.get(i).end()) // FOUND
+					l_id = key->second;
+				else
+				{
+					// NOT FOUND
+
+					// It must be always found, if not it mean that the processor has no-idea of
+					// what is stored and conseguently do not know how to unpack, print a critical error
+					// and return
+
+					std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Critical, cannot unpack object, because received data cannot be interpreted\n";
+
+					return;
+				}
+
+				// Get the internal ghost box associated with the packed information
+				Box<dim,size_t> box = ig_box.get(i).bid.get(l_id).box;
+				size_t sub_id = ig_box.get(i).bid.get(l_id).sub;
+				box -= gdb_ext.get(sub_id).origin.template convertPoint<size_t>();
+
+				// sub-grid where to unpack
+				grid_key_dx_iterator_sub<dim> sub2(loc_grid.get(sub_id).getGrid(),box.getKP1(),box.getKP2());
+
+				grid_unpack_with_prp<op,prp_object,device_grid,Memory>::template unpacking<prp...>(prRecv_prp,sub2,loc_grid.get(sub_id),ps);
+
+				// Unpack
+//				Unpacker<device_grid,HeapMemory>::template unpack<prp...>(prRecv_prp,sub2,loc_grid.get(sub_id),ps);
+			}
+		}
 	}
 
 	/*! \brief Constructor
diff --git a/src/Grid/grid_dist_id_unit_test.cpp b/src/Grid/grid_dist_id_unit_test.cpp
index 3911c9aa..568bbbbc 100644
--- a/src/Grid/grid_dist_id_unit_test.cpp
+++ b/src/Grid/grid_dist_id_unit_test.cpp
@@ -1359,6 +1359,126 @@ void Test3D_periodic(const Box<3,float> & domain, long int k)
 	}
 }
 
+// Test grid periodic
+
+void Test3D_periodic_put(const Box<3,float> & domain, long int k)
+{
+	Vcluster & v_cl = create_vcluster();
+
+	if ( v_cl.getProcessingUnits() > 32 )
+		return;
+
+	long int big_step = k / 30;
+	big_step = (big_step == 0)?1:big_step;
+	long int small_step = 21;
+
+	print_test( "Testing grid periodic put k<=",k);
+
+	// 3D test
+	for ( ; k >= 2 ; k-= (k > 2*big_step)?big_step:small_step )
+	{
+		BOOST_TEST_CHECKPOINT( "Testing grid periodick<=" << k );
+
+		// grid size
+		size_t sz[3];
+		sz[0] = k;
+		sz[1] = k;
+		sz[2] = k;
+
+		// Ghost
+		Ghost<3,long int> g(1);
+
+		// periodicity
+		periodicity<3> pr = {{PERIODIC,PERIODIC,PERIODIC}};
+
+		// Distributed grid with id decomposition
+		grid_dist_id<3, float, aggregate<long int>, CartDecomposition<3,float>> g_dist(sz,domain,g,pr);
+
+		// check the consistency of the decomposition
+		bool val = g_dist.getDecomposition().check_consistency();
+		BOOST_REQUIRE_EQUAL(val,true);
+
+		// Grid sm
+		grid_sm<3,void> info(sz);
+
+		size_t count = 0;
+
+		{
+		auto dom = g_dist.getDomainIterator();
+
+		while (dom.isNext())
+		{
+			auto key = dom.get();
+
+			g_dist.template get<0>(key) = -6.0;
+
+			// Count the points
+			count++;
+
+			++dom;
+		}
+		}
+
+		// Set to zero the full grid
+
+		{
+		auto dom = g_dist.getDomainIterator();
+
+		while (dom.isNext())
+		{
+			auto key = dom.get();
+
+			g_dist.template get<0>(key.move(0,1)) += 1.0;
+			g_dist.template get<0>(key.move(0,-1)) += 1.0;
+			g_dist.template get<0>(key.move(1,1)) += 1.0;
+			g_dist.template get<0>(key.move(1,-1)) += 1.0;
+			g_dist.template get<0>(key.move(2,1)) += 1.0;
+			g_dist.template get<0>(key.move(2,-1)) += 1.0;
+
+			++dom;
+		}
+		}
+
+		bool correct = true;
+
+		// Domain + Ghost iterator
+		auto dom_gi = g_dist.getDomainIterator();
+
+		while (dom_gi.isNext())
+		{
+			auto key = dom_gi.get();
+
+			correct &= (g_dist.template get<0>(key) == 0);
+
+			++dom_gi;
+		}
+
+		g_dist.ghost_put<add_,0>();
+
+		if (count != 0)
+			BOOST_REQUIRE_EQUAL(correct, false);
+
+		// sync the ghosts
+		g_dist.ghost_get<0>();
+
+		correct = true;
+
+		// Domain + Ghost iterator
+		auto dom_gi2 = g_dist.getDomainIterator();
+
+		while (dom_gi2.isNext())
+		{
+			auto key = dom_gi2.get();
+
+			correct &= (g_dist.template get<0>(key) == 0);
+
+			++dom_gi2;
+		}
+
+		BOOST_REQUIRE_EQUAL(correct, true);
+	}
+}
+
 void Test_grid_copy(const Box<3,float> & domain, long int k)
 {
 	typedef Point_test<float> p;
@@ -1620,6 +1740,17 @@ BOOST_AUTO_TEST_CASE( grid_1d_test )
 	Test1D(domain1,k);
 }
 
+BOOST_AUTO_TEST_CASE( grid_dist_id_periodic_put_test )
+{
+	// Domain
+	Box<3,float> domain3({0.0,0.0,0.0},{1.0,1.0,1.0});
+
+	long int k = 128*128*128*create_vcluster().getProcessingUnits();
+	k = std::pow(k, 1/3.);
+
+	Test3D_periodic_put(domain3,k);
+}
+
 BOOST_AUTO_TEST_SUITE_END()
 
 #endif
diff --git a/src/Grid/grid_dist_key.hpp b/src/Grid/grid_dist_key.hpp
index f3de7509..b6d1a4b3 100644
--- a/src/Grid/grid_dist_key.hpp
+++ b/src/Grid/grid_dist_key.hpp
@@ -20,6 +20,16 @@ class grid_dist_key_dx
 
 public:
 
+	/*! \brief Set the local grid
+	 *
+	 * \param sub local grid
+	 *
+	 */
+	inline void setSub(size_t sub)
+	{
+		g_c = sub;
+	}
+
 	/*! \brief Get the local grid
 	 *
 	 * \return the id of the local grid
@@ -99,11 +109,20 @@ public:
 		return grid_dist_key_dx<dim>(getSub(),key);
 	}
 
+	/*! \brief Constructor set the sub-domain grid and the position in local coordinates
+	 *
+	 * \param g_c sub-domain
+	 * \param key key
+	 *
+	 */
 	inline grid_dist_key_dx(int g_c, const grid_key_dx<dim> & key)
 	:g_c(g_c),key(key)
 	{
 	}
 
+	//! Constructor
+	inline grid_dist_key_dx(){}
+
 	std::string to_string()
 	{
 		std::stringstream str;
diff --git a/src/Grid/grid_dist_util.hpp b/src/Grid/grid_dist_util.hpp
index d16bf8bb..77e2383c 100644
--- a/src/Grid/grid_dist_util.hpp
+++ b/src/Grid/grid_dist_util.hpp
@@ -200,11 +200,17 @@ template <unsigned int dim> struct e_lbox_id
 	//! Box defining the external ghost box in local coordinates
 	::Box<dim,long int> box;
 
-	//! Sector position of the local external ghost box
-	comb<dim> cmb;
+	//! Has this external ghost box initialized
+	bool initialized = false;
 
 	//! sub_id in which sub-domain this box live
 	size_t sub;
+
+	//! external ghost box linked to this internal ghost box
+	size_t k;
+
+	//! Sector position of the local external ghost box
+	comb<dim> cmb;
 };
 
 /*! \brief Per-processor Internal ghost box
diff --git a/src/Makefile.am b/src/Makefile.am
index cc194ca4..c4083ddd 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -9,7 +9,7 @@ pdata_LDADD = $(LINKLIBS) -lparmetis -lmetis
 nobase_include_HEADERS = Decomposition/CartDecomposition.hpp Decomposition/CartDecomposition_ext.hpp  Decomposition/common.hpp Decomposition/Decomposition.hpp  Decomposition/ie_ghost.hpp \
          Decomposition/Domain_NN_calculator_cart.hpp Decomposition/nn_processor.hpp Decomposition/ie_loc_ghost.hpp Decomposition/ORB.hpp \
          Graph/CartesianGraphFactory.hpp \
-         Grid/grid_dist_id.hpp Grid/Iterators/grid_dist_id_iterator_util.hpp Grid/Iterators/grid_dist_id_iterator_dec.hpp Grid/Iterators/grid_dist_id_iterator_dec_skin.hpp Grid/grid_dist_util.hpp  Grid/Iterators/grid_dist_id_iterator_sub.hpp Grid/Iterators/grid_dist_id_iterator.hpp Grid/grid_dist_key.hpp Grid/staggered_dist_grid.hpp Grid/staggered_dist_grid_util.hpp Grid/staggered_dist_grid_copy.hpp \
+         Grid/grid_dist_id.hpp Grid/grid_dist_id_comm.hpp Grid/Iterators/grid_dist_id_iterator_util.hpp Grid/Iterators/grid_dist_id_iterator_dec.hpp Grid/Iterators/grid_dist_id_iterator_dec_skin.hpp Grid/grid_dist_util.hpp  Grid/Iterators/grid_dist_id_iterator_sub.hpp Grid/Iterators/grid_dist_id_iterator.hpp Grid/grid_dist_key.hpp Grid/staggered_dist_grid.hpp Grid/staggered_dist_grid_util.hpp Grid/staggered_dist_grid_copy.hpp \
          Vector/se_class3_vector.hpp  Vector/vector_dist_multiphase_functions.hpp Vector/vector_dist_comm.hpp Vector/vector_dist.hpp Vector/vector_dist_ofb.hpp Vector/Iterators/vector_dist_iterator.hpp Vector/vector_dist_key.hpp \
          config/config.h \
          example.mk \
diff --git a/src/Vector/se_class3_vector.hpp b/src/Vector/se_class3_vector.hpp
index 4099774d..11611166 100644
--- a/src/Vector/se_class3_vector.hpp
+++ b/src/Vector/se_class3_vector.hpp
@@ -222,7 +222,7 @@ struct propCheckINF
 	{
 		typedef typename boost::mpl::at<typename vector::value_type::type,boost::mpl::int_<T::value> >::type type_to_check;
 
-		bool snn = typeCheck<type_to_check,std::is_fundamental<type_to_check>::value>::isInf(data.template getProp<T::value>(id));
+		bool snn = typeCheck<type_to_check,std::is_fundamental<type_to_check>::value>::isInf(data.template getPropNC<T::value>(id));
 
 		if (snn == true)
 		{
@@ -383,15 +383,18 @@ class se_class3_vector
 				auto p = it.get();
 
 				for (size_t i = 0 ; i < Np_real ; i++)
-					vd.template getProp<Np+SE3_STATUS>(p)[i] = UNINITIALIZED;
+					vd.template getPropNC<Np+SE3_STATUS>(p)[i] = UNINITIALIZED;
 
-				vd.template getProp<Np+SE3_TYPE>(p) = INSIDE;
+				vd.template getPropNC<Np+SE3_TYPE>(p) = INSIDE;
 
 				++it;
 			}
 
 			for (size_t i = 0 ; i < Np_real ; i++)
+			{
 				sync[GHOST][i] = NOTSYNC;
+				sync[HALO][i] = SYNC;
+			}
 		}
 
 		template <unsigned int ... prp> void ghost_get_pre(size_t opt)
@@ -411,7 +414,7 @@ class se_class3_vector
 
 				for (size_t i = 0 ; i < sizeof...(prp) ; i++)
 				{
-					if (vd.template getProp<Np+SE3_STATUS>(p)[gg[i]] == DIRTY)
+					if (vd.template getPropNC<Np+SE3_STATUS>(p)[gg[i]] == DIRTY)
 					{
 						std::cerr << __FILE__ << ":" << __LINE__ << " Error the ghost has been written and ghost_get will overwrite your changes" << std::endl;
 						ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
@@ -422,7 +425,7 @@ class se_class3_vector
 				{
 					for (size_t i = 0 ; i < non_NP.size() ; i++)
 					{
-						if (vd.template getProp<Np+SE3_STATUS>(p)[non_NP.get(i)] == DIRTY)
+						if (vd.template getPropNC<Np+SE3_STATUS>(p)[non_NP.get(i)] == DIRTY)
 						{
 							std::cerr << __FILE__ << ":" << __LINE__ << " Error the it seem that the property=" << getPrpName(non_NP.get(i)) << " has been written and ghost_get will destroy such changes" << std::endl;
 							ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
@@ -448,11 +451,11 @@ class se_class3_vector
 
 				for (size_t i = 0 ; i < sizeof...(prp) ; i++)
 				{
-					if (vd.template getProp<Np+SE3_STATUS>(p)[gg[i]] == DIRTY)
-						vd.template getProp<Np+SE3_STATUS>(p)[gg[i]] = CLEAN;
+					if (vd.template getPropNC<Np+SE3_STATUS>(p)[gg[i]] == DIRTY)
+						vd.template getPropNC<Np+SE3_STATUS>(p)[gg[i]] = CLEAN;
 				}
 
-				vd.template getProp<Np+SE3_TYPE>(p) = GHOST;
+				vd.template getPropNC<Np+SE3_TYPE>(p) = GHOST;
 
 				++it2;
 			}
@@ -469,7 +472,7 @@ class se_class3_vector
 					auto p = it.get();
 
 					for (size_t i = 0 ; i < non_NP.size() ; i++)
-						vd.template getProp<Np+SE3_STATUS>(p)[non_NP.get(i)] = UNINITIALIZED;
+						vd.template getPropNC<Np+SE3_STATUS>(p)[non_NP.get(i)] = UNINITIALIZED;
 
 					++it;
 				}
@@ -548,15 +551,15 @@ class se_class3_vector
 
 				for (size_t j = 0 ; j < Np ; j++)
 				{
-					if (vd.template getProp<Np+SE3_STATUS>(p)[j] == DIRTY)
-						vd.template getProp<Np+SE3_STATUS>(p)[j] = CLEAN;
+					if (vd.template getPropNC<Np+SE3_STATUS>(p)[j] == DIRTY)
+						vd.template getPropNC<Np+SE3_STATUS>(p)[j] = CLEAN;
 				}
 
 #ifdef CHECK_FOR_POSINF
 
-				if ( std::isinf(vd.getPos(p)[0]) || std::isinf(vd.getPos(p)[1]) || std::isinf(vd.getPos(p)[2]) )
+				if ( std::isinf(vd.getPosNC(p)[0]) || std::isinf(vd.getPosNC(p)[1]) || std::isinf(vd.getPosNC(p)[2]) )
 				{
-					std::cerr << __FILE__ << ":" << __LINE__ << " error detected INF in position for particle p=" << p.getKey() << " of type=" << getParticleTypeString(vd.template getProp<Np+SE3_TYPE>(p)) << std::endl;
+					std::cerr << __FILE__ << ":" << __LINE__ << " error detected INF in position for particle p=" << p.getKey() << " of type=" << getParticleTypeString(vd.template getPropNC<Np+SE3_TYPE>(p)) << std::endl;
 					ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
 				}
 
@@ -564,9 +567,9 @@ class se_class3_vector
 
 #ifdef CHECKFOR_POSNAN
 
-				if ( std::isnan(vd.getPos(p)[0]) || std::isnan(vd.getPos(p)[1]) || std::isnan(vd.getPos(p)[2]) )
+				if ( std::isnan(vd.getPosNC(p)[0]) || std::isnan(vd.getPosNC(p)[1]) || std::isnan(vd.getPosNC(p)[2]) )
 				{
-					std::cerr << __FILE__ << ":" << __LINE__ << " error detected NAN in position for particle p=" << p.getKey() << " of type=" << getParticleTypeString(vd.template getProp<Np+SE3_TYPE>(p)) << std::endl;
+					std::cerr << __FILE__ << ":" << __LINE__ << " error detected NAN in position for particle p=" << p.getKey() << " of type=" << getParticleTypeString(vd.template getPropNC<Np+SE3_TYPE>(p)) << std::endl;
 					ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
 				}
 
@@ -606,7 +609,7 @@ class se_class3_vector
 
 				for (size_t j = 0 ; j < Np ; j++)
 				{
-					if (vd.template getProp<Np+SE3_STATUS>(p)[j] == DIRTY)
+					if (vd.template getPropNC<Np+SE3_STATUS>(p)[j] == DIRTY)
 					{
 						std::cerr << __FILE__ << ":" << __LINE__ << " error it seem that ghost has been filled with information that we are going to destroy with the map call " << std::endl;
 					}
@@ -632,9 +635,9 @@ class se_class3_vector
 
 				for (size_t j = 0 ; j < Np ; j++)
 				{
-					Point<vector::dims,typename vector::stype> xp = vd.getPos(p);
+					Point<vector::dims,typename vector::stype> xp = vd.getPosNC(p);
 
-					vd.template getProp<Np+SE3_TYPE>(p) = getParticleType(xp,p.getKey(),vd);
+					vd.template getPropNC<Np+SE3_TYPE>(p) = getParticleType(xp,p.getKey(),vd);
 				}
 
 				++it;
@@ -643,10 +646,10 @@ class se_class3_vector
 
 		template<unsigned int prp> void read(const vector & vd, size_t p) const
 		{
-			if (vd.template getProp<Np+SE3_STATUS>(p)[prp] == UNINITIALIZED)
+			if (vd.template getPropNC<Np+SE3_STATUS>(p)[prp] == UNINITIALIZED)
 			{
 				std::stringstream str;
-				std::string type_str = getParticleTypeString(vd.template getProp<Np+SE3_TYPE>(p));
+				std::string type_str = getParticleTypeString(vd.template getPropNC<Np+SE3_TYPE>(p));
 
 				if (prp == Np_real)
 					str << __FILE__ << ":" << __LINE__ << " Error you are reading the particle " << p << " of type " << type_str << " the position. But it result to be uninitialized" << std::endl;
@@ -658,18 +661,18 @@ class se_class3_vector
 				ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
 			}
 
-			if (vd.template getProp<Np+SE3_STATUS>(p)[prp] == DIRTY)
+			if (vd.template getPropNC<Np+SE3_STATUS>(p)[prp] == DIRTY)
 			{
 				std::cerr << __FILE__ << ":" << __LINE__ << " Warning you are reading from an particle that has been changed already in the same cycle" << std::endl;
 				ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
 			}
 
-			if (vd.template getProp<Np+SE3_TYPE>(p) == GHOST || vd.template getProp<Np+SE3_TYPE>(p) == HALO)
+			if (vd.template getPropNC<Np+SE3_TYPE>(p) == GHOST || vd.template getPropNC<Np+SE3_TYPE>(p) == HALO)
 			{
 				// if we read from the ghost we have to ensure that the ghost is in
 				// sync in particular that the state of the halo is CLEAR
 
-				if (sync[vd.template getProp<Np+SE3_TYPE>(p)][prp] != SYNC)
+				if (sync[vd.template getPropNC<Np+SE3_TYPE>(p)][prp] != SYNC)
 				{
 					std::cerr << __FILE__ << ":" << __LINE__ << " Error it seem that you are reading from a ghost the property=" << getPrpName(prp) << " but it seem it is changed from the last ghost_get. It seems that is missing a ghost_get" << std::endl;
 					ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
@@ -679,12 +682,12 @@ class se_class3_vector
 
 		template<unsigned int prp> void write(vector & vd, size_t p)
 		{
-			vd.template getProp<Np+SE3_STATUS>(p)[prp] = DIRTY;
+			vd.template getPropNC<Np+SE3_STATUS>(p)[prp] = DIRTY;
 			if (p >= vd.size_local())
 				vd.get_se_class3().template setHaloOutSync<prp>();
 			else
 			{
-				if (vd.template getProp<Np+SE3_TYPE>(p) == HALO)
+				if (vd.template getPropNC<Np+SE3_TYPE>(p) == HALO)
 					vd.get_se_class3().template setGhostOutSync<prp>();
 			}
 
diff --git a/src/Vector/vector_dist.hpp b/src/Vector/vector_dist.hpp
index 781e12bd..5cbc554d 100644
--- a/src/Vector/vector_dist.hpp
+++ b/src/Vector/vector_dist.hpp
@@ -8,8 +8,6 @@
 #ifndef VECTOR_HPP_
 #define VECTOR_HPP_
 
-#include "../../openfpm_data/src/NN/CellList/CellListFast_gen.hpp"
-#include "../../openfpm_data/src/NN/CellList/CellListFast_gen.hpp"
 #include "HDF5_wr/HDF5_wr.hpp"
 #include "VCluster/VCluster.hpp"
 #include "Space/Shape/Point.hpp"
@@ -18,6 +16,7 @@
 #include "Vector/vector_dist_key.hpp"
 #include "memory/PtrMemory.hpp"
 #include "NN/CellList/CellList.hpp"
+#include "NN/CellList/CellListFast_gen.hpp"
 #include "util/common.hpp"
 #include "util/object_util.hpp"
 #include "memory/ExtPreAlloc.hpp"
@@ -77,7 +76,7 @@ struct gcl
 	 */
 	static inline CellL get(Vector & vd, const St & r_cut, const Ghost<dim,St> & g)
 	{
-		return vd.getCellList(r_cut);
+		return vd.template getCellList<CellL>(r_cut);
 	}
 };
 
@@ -100,6 +99,10 @@ struct gcl<dim,St,CellList_gen<dim, St, Process_keys_hilb,Mem_type, shift<dim, S
 	}
 };
 
+#define CELL_MEMFAST(dim,St) CellList_gen<dim, St, Process_keys_lin, Mem_fast<dim,St>, shift<dim, St> >
+#define CELL_MEMBAL(dim,St) CellList_gen<dim, St, Process_keys_lin, Mem_bal<dim,St>, shift<dim, St> >
+#define CELL_MEMMW(dim,St) CellList_gen<dim, St, Process_keys_lin, Mem_mw<dim,St>, shift<dim, St> >
+
 /*! \brief Distributed vector
  *
  * This class reppresent a distributed vector, the distribution of the structure
@@ -519,6 +522,124 @@ public:
 
 #endif
 
+///////////////////// Read and write with no check
+
+	/*! \brief Get the position of an element
+	 *
+	 * see the vector_dist iterator usage to get an element key
+	 *
+	 * \param vec_key element
+	 *
+	 * \return the position of the element in space
+	 *
+	 */
+	inline auto getPosNC(vect_dist_key_dx vec_key) -> decltype(v_pos.template get<0>(vec_key.getKey()))
+	{
+		return v_pos.template get<0>(vec_key.getKey());
+	}
+
+	/*! \brief Get the position of an element
+	 *
+	 * see the vector_dist iterator usage to get an element key
+	 *
+	 * \param vec_key element
+	 *
+	 * \return the position of the element in space
+	 *
+	 */
+	inline auto getPosNC(vect_dist_key_dx vec_key) const -> decltype(v_pos.template get<0>(vec_key.getKey()))
+	{
+		return v_pos.template get<0>(vec_key.getKey());
+	}
+
+	/*! \brief Get the position of an element
+	 *
+	 * see the vector_dist iterator usage to get an element key
+	 *
+	 * \param vec_key element
+	 *
+	 * \return the position of the element in space
+	 *
+	 */
+	inline auto getPosNC(size_t vec_key) -> decltype(v_pos.template get<0>(vec_key))
+	{
+		return v_pos.template get<0>(vec_key);
+	}
+
+	/*! \brief Get the position of an element
+	 *
+	 * see the vector_dist iterator usage to get an element key
+	 *
+	 * \param vec_key element
+	 *
+	 * \return the position of the element in space
+	 *
+	 */
+	inline auto getPosNC(size_t vec_key) const -> decltype(v_pos.template get<0>(vec_key))
+	{
+		return v_pos.template get<0>(vec_key);
+	}
+
+	/*! \brief Get the property of an element
+	 *
+	 * see the vector_dist iterator usage to get an element key
+	 *
+	 * \tparam id property id
+	 * \param vec_key vector element
+	 *
+	 * \return return the selected property of the vector element
+	 *
+	 */
+	template<unsigned int id> inline auto getPropNC(vect_dist_key_dx vec_key) -> decltype(v_prp.template get<id>(vec_key.getKey()))
+	{
+		return v_prp.template get<id>(vec_key.getKey());
+	}
+
+	/*! \brief Get the property of an element
+	 *
+	 * see the vector_dist iterator usage to get an element key
+	 *
+	 * \tparam id property id
+	 * \param vec_key vector element
+	 *
+	 * \return return the selected property of the vector element
+	 *
+	 */
+	template<unsigned int id> inline auto getPropNC(vect_dist_key_dx vec_key) const -> const decltype(v_prp.template get<id>(vec_key.getKey()))
+	{
+		return v_prp.template get<id>(vec_key.getKey());
+	}
+
+	/*! \brief Get the property of an element
+	 *
+	 * see the vector_dist iterator usage to get an element key
+	 *
+	 * \tparam id property id
+	 * \param vec_key vector element
+	 *
+	 * \return return the selected property of the vector element
+	 *
+	 */
+	template<unsigned int id> inline auto getPropNC(size_t vec_key) -> decltype(v_prp.template get<id>(vec_key))
+	{
+		return v_prp.template get<id>(vec_key);
+	}
+
+	/*! \brief Get the property of an element
+	 *
+	 * see the vector_dist iterator usage to get an element key
+	 *
+	 * \tparam id property id
+	 * \param vec_key vector element
+	 *
+	 * \return return the selected property of the vector element
+	 *
+	 */
+	template<unsigned int id> inline auto getPropNC(size_t vec_key) const -> const decltype(v_prp.template get<id>(vec_key))
+	{
+		return v_prp.template get<id>(vec_key);
+	}
+
 ///////////////////// Read and Write function
 
 	/*! \brief Get the position of an element
@@ -767,7 +888,7 @@ public:
 		Ghost<dim,St> g = getDecomposition().getGhost();
 		g.magnify(1.013);
 
-		return getCellList(r_cut, g);
+		return getCellList<CellL>(r_cut, g);
 	}
 
 	/*! \brief Construct an hilbert cell list starting from the stored particles
diff --git a/src/Vector/vector_dist_cell_list_tests.hpp b/src/Vector/vector_dist_cell_list_tests.hpp
index 71afdf2c..1179b255 100644
--- a/src/Vector/vector_dist_cell_list_tests.hpp
+++ b/src/Vector/vector_dist_cell_list_tests.hpp
@@ -1691,4 +1691,114 @@ BOOST_AUTO_TEST_CASE( vector_dist_checking_unloaded_processors )
 	}
 }
 
+BOOST_AUTO_TEST_CASE( vector_dist_cell_list_multi_type )
+{
+	Vcluster & v_cl = create_vcluster();
+
+	if (v_cl.getProcessingUnits() > 24)
+		return;
+
+	float L = 1000.0;
+
+    // set the seed
+	// create the random generator engine
+	std::srand(0);
+    std::default_random_engine eg;
+    std::uniform_real_distribution<float> ud(-L,L);
+
+    long int k = 4096 * v_cl.getProcessingUnits();
+
+	long int big_step = k / 4;
+	big_step = (big_step == 0)?1:big_step;
+
+	print_test("Testing 3D periodic vector symmetric cell-list k=",k);
+	BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric cell-list k=" << k );
+
+	Box<3,float> box({-L,-L,-L},{L,L,L});
+
+	// Boundary conditions
+	size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+
+	float r_cut = 100.0;
+
+	// ghost
+	Ghost<3,float> ghost(r_cut);
+
+	typedef  aggregate<size_t> part_prop;
+
+	// Distributed vector
+	vector_dist<3,float, part_prop > vd(k,box,bc,ghost);
+
+	auto it = vd.getIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		vd.getPos(key)[0] = ud(eg);
+		vd.getPos(key)[1] = ud(eg);
+		vd.getPos(key)[2] = ud(eg);
+
+		++it;
+	}
+
+	vd.map();
+
+	// sync the ghost
+	vd.ghost_get<0>();
+
+
+	bool ret = true;
+
+	// We take different type of Cell-list
+	auto NN = vd.getCellList<CELL_MEMFAST(3,float)>(r_cut);
+	auto NN2 = vd.getCellList<CELL_MEMBAL(3,float)>(r_cut);
+	auto NN3 = vd.getCellList<CELL_MEMMW(3,float)>(r_cut);
+
+	auto p_it = vd.getDomainIterator();
+
+	while (p_it.isNext())
+	{
+		auto p = p_it.get();
+
+		Point<3,float> xp = vd.getPos(p);
+
+		auto Np = NN.getNNIterator(NN.getCell(xp));
+		auto Np2 = NN2.getNNIterator(NN2.getCell(xp));
+		auto Np3 = NN3.getNNIterator(NN3.getCell(xp));
+
+		while (Np.isNext())
+		{
+			// first all cell-list must agree
+
+			ret &= (Np.isNext() == Np2.isNext()) && (Np3.isNext() == Np.isNext());
+
+			if (ret == false)
+				break;
+
+			auto q = Np.get();
+			auto q2 = Np2.get();
+			auto q3 = Np3.get();
+
+			ret &= (q == q2) && (q == q3);
+
+			if (ret == false)
+				break;
+
+			++Np;
+			++Np2;
+			++Np3;
+		}
+
+		ret &= (Np.isNext() == Np2.isNext()) && (Np.isNext() == Np3.isNext());
+
+		if (ret == false)
+			break;
+
+		++p_it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+}
+
 #endif /* SRC_VECTOR_VECTOR_DIST_CELL_LIST_TESTS_HPP_ */
diff --git a/src/Vector/vector_dist_key.hpp b/src/Vector/vector_dist_key.hpp
index c5c7a878..fd79ba14 100644
--- a/src/Vector/vector_dist_key.hpp
+++ b/src/Vector/vector_dist_key.hpp
@@ -68,6 +68,12 @@ public:
 	inline vect_dist_key_dx()
 	{
 	}
+
+	//! Default constructor
+	inline vect_dist_key_dx(size_t key)
+	:key(key)
+	{
+	}
 };
 
 
-- 
GitLab