diff --git a/src/.deps/pdata-VCluster.Po b/src/.deps/pdata-VCluster.Po
index 5feb7880e5a35c697bddf5149f000cae7d6781e3..7d351b7406e69067fe92e63bdb157027349770f8 100644
--- a/src/.deps/pdata-VCluster.Po
+++ b/src/.deps/pdata-VCluster.Po
@@ -939,7 +939,7 @@ pdata-VCluster.o: ../../OpenFPM_vcluster/src/VCluster.cpp \
  /usr/include/boost/multi_array/view.hpp \
  /usr/include/boost/functional.hpp /usr/include/boost/call_traits.hpp \
  /usr/include/boost/detail/call_traits.hpp \
- ../../OpenFPM_data/src/ct_array.hpp \
+ ../../OpenFPM_data/src/util/ct_array.hpp \
  ../../OpenFPM_data/src/memory_array.hpp \
  ../../OpenFPM_data/src/memory.hpp ../../OpenFPM_data/src/memory.hpp \
  ../../OpenFPM_data/src/util/meta_copy.hpp \
@@ -955,14 +955,17 @@ pdata-VCluster.o: ../../OpenFPM_vcluster/src/VCluster.cpp \
  ../../OpenFPM_data/src/Grid/grid_key.hpp \
  ../../OpenFPM_data/src/Grid/comb.hpp \
  ../../OpenFPM_data/src/Grid/grid_key_expression.hpp \
- ../../OpenFPM_data/src/Space/Shape/Point.hpp \
  ../../OpenFPM_data/src/Grid/Encap.hpp \
+ ../../OpenFPM_data/src/Space/Shape/Point.hpp \
  ../../OpenFPM_data/src/Grid/grid_key.hpp \
  ../../OpenFPM_data/src/Grid/Encap.hpp \
  ../../OpenFPM_data/src/memory_array.hpp \
  ../../OpenFPM_devices/src/memory/HeapMemory.hpp \
  ../../OpenFPM_data/src/Vector/vect_isel.hpp \
  ../../OpenFPM_data/src/common.hpp \
+ ../../OpenFPM_data/src/util/object_s_di.hpp \
+ ../../OpenFPM_data/src/util/for_each_ref.hpp \
+ /usr/include/boost/fusion/include/size.hpp \
  ../../OpenFPM_data/src/Vector/map_vector_std.hpp \
  ../../OpenFPM_vcluster/src/MPI_IallreduceW.hpp
 
@@ -3024,7 +3027,7 @@ pdata-VCluster.o: ../../OpenFPM_vcluster/src/VCluster.cpp \
 
 /usr/include/boost/detail/call_traits.hpp:
 
-../../OpenFPM_data/src/ct_array.hpp:
+../../OpenFPM_data/src/util/ct_array.hpp:
 
 ../../OpenFPM_data/src/memory_array.hpp:
 
@@ -3058,10 +3061,10 @@ pdata-VCluster.o: ../../OpenFPM_vcluster/src/VCluster.cpp \
 
 ../../OpenFPM_data/src/Grid/grid_key_expression.hpp:
 
-../../OpenFPM_data/src/Space/Shape/Point.hpp:
-
 ../../OpenFPM_data/src/Grid/Encap.hpp:
 
+../../OpenFPM_data/src/Space/Shape/Point.hpp:
+
 ../../OpenFPM_data/src/Grid/grid_key.hpp:
 
 ../../OpenFPM_data/src/Grid/Encap.hpp:
@@ -3074,6 +3077,12 @@ pdata-VCluster.o: ../../OpenFPM_vcluster/src/VCluster.cpp \
 
 ../../OpenFPM_data/src/common.hpp:
 
+../../OpenFPM_data/src/util/object_s_di.hpp:
+
+../../OpenFPM_data/src/util/for_each_ref.hpp:
+
+/usr/include/boost/fusion/include/size.hpp:
+
 ../../OpenFPM_data/src/Vector/map_vector_std.hpp:
 
 ../../OpenFPM_vcluster/src/MPI_IallreduceW.hpp:
diff --git a/src/.deps/pdata-main.Po b/src/.deps/pdata-main.Po
index 8b6e781b90f945d9534d2e71515fc1f0a1695674..01321ba088d4669fc2e33835e7db436c005e0b06 100644
--- a/src/.deps/pdata-main.Po
+++ b/src/.deps/pdata-main.Po
@@ -901,7 +901,7 @@ pdata-main.o: main.cpp /usr/include/stdc-predef.h \
  /usr/include/boost/multi_array/view.hpp \
  /usr/include/boost/functional.hpp /usr/include/boost/call_traits.hpp \
  /usr/include/boost/detail/call_traits.hpp \
- ../../OpenFPM_data/src/ct_array.hpp \
+ ../../OpenFPM_data/src/util/ct_array.hpp \
  ../../OpenFPM_data/src/memory_array.hpp \
  ../../OpenFPM_data/src/memory.hpp ../../OpenFPM_data/src/memory.hpp \
  ../../OpenFPM_data/src/util/meta_copy.hpp \
@@ -920,14 +920,18 @@ pdata-main.o: main.cpp /usr/include/stdc-predef.h \
  ../../OpenFPM_data/src/Grid/grid_key.hpp \
  ../../OpenFPM_data/src/Grid/comb.hpp \
  ../../OpenFPM_data/src/Grid/grid_key_expression.hpp \
- ../../OpenFPM_data/src/Space/Shape/Point.hpp \
  ../../OpenFPM_data/src/Grid/Encap.hpp \
+ ../../OpenFPM_data/src/Space/Shape/Point.hpp \
  ../../OpenFPM_data/src/Grid/grid_key.hpp \
  ../../OpenFPM_data/src/Grid/Encap.hpp \
  ../../OpenFPM_data/src/memory_array.hpp \
  ../../OpenFPM_devices/src/memory/HeapMemory.hpp \
  ../../OpenFPM_data/src/Vector/vect_isel.hpp \
- ../../OpenFPM_data/src/common.hpp /home/i-bird/MPI/include/mpi.h \
+ ../../OpenFPM_data/src/common.hpp \
+ ../../OpenFPM_data/src/util/object_s_di.hpp \
+ ../../OpenFPM_data/src/util/for_each_ref.hpp \
+ /usr/include/boost/fusion/include/size.hpp \
+ /home/i-bird/MPI/include/mpi.h \
  /home/i-bird/MPI/include/mpi_portable_platform.h \
  /home/i-bird/MPI/include/openmpi/ompi/mpi/cxx/mpicxx.h \
  /home/i-bird/MPI/include/openmpi/ompi/mpi/cxx/constants.h \
@@ -1291,11 +1295,12 @@ pdata-main.o: main.cpp /usr/include/stdc-predef.h \
  ../../OpenFPM_data/src/Space/Shape/Box.hpp \
  ../../OpenFPM_data/src/Space/Ghost.hpp Grid/grid_dist_id_iterator.hpp \
  Grid/grid_dist_key.hpp ../../OpenFPM_data/src/Point_test.hpp \
- ../../OpenFPM_data/src/base_type.hpp Decomposition/CartDecomposition.hpp \
- Decomposition/Decomposition.hpp ../../OpenFPM_data/src/global_const.hpp \
- SubdomainGraphNodes.hpp metis_util.hpp \
- ../../metis_install/include/metis.h /usr/include/inttypes.h \
- ../../OpenFPM_IO/src/VTKWriter.hpp \
+ ../../OpenFPM_data/src/base_type.hpp \
+ ../../OpenFPM_data/src/Point_orig.hpp \
+ Decomposition/CartDecomposition.hpp Decomposition/Decomposition.hpp \
+ ../../OpenFPM_data/src/global_const.hpp SubdomainGraphNodes.hpp \
+ metis_util.hpp ../../metis_install/include/metis.h \
+ /usr/include/inttypes.h ../../OpenFPM_IO/src/VTKWriter.hpp \
  ../../OpenFPM_IO/src/VTKWriter_graph.hpp \
  ../../OpenFPM_IO/src/VTKWriter_vector_box.hpp \
  /usr/include/boost/math/special_functions/pow.hpp \
@@ -1362,13 +1367,14 @@ pdata-main.o: main.cpp /usr/include/stdc-predef.h \
  /usr/include/boost/config/abi_prefix.hpp \
  /usr/include/boost/config/abi_suffix.hpp dec_optimizer.hpp \
  ../../OpenFPM_data/src/NN/CellList/CellDecomposer.hpp \
+ ../../OpenFPM_data/src/Space/Matrix.hpp \
  /usr/include/c++/4.8.3/unordered_map \
  /usr/include/c++/4.8.3/bits/hashtable.h \
  /usr/include/c++/4.8.3/bits/hashtable_policy.h \
  /usr/include/c++/4.8.3/bits/unordered_map.h \
  ../../OpenFPM_data/src/NN/CellList/CellList.hpp \
- ../../OpenFPM_data/src/NN/CellList/CellListFast.hpp \
  ../../OpenFPM_data/src/NN/CellList/CellDecomposer.hpp \
+ ../../OpenFPM_data/src/NN/CellList/CellListFast.hpp \
  ../../OpenFPM_data/src/NN/CellList/CellNNIterator.hpp \
  ../../OpenFPM_data/src/NN/CellList/CellListBal.hpp \
  ../../OpenFPM_data/src/NN/CellList/CellListMem.hpp \
@@ -1387,9 +1393,8 @@ pdata-main.o: main.cpp /usr/include/stdc-predef.h \
  ../../OpenFPM_devices/src/memory/PtrMemory.hpp \
  ../../OpenFPM_data/src/util/object_util.hpp \
  ../../OpenFPM_data/src/util/object_creator.hpp \
- ../../OpenFPM_data/src/util/object_copy.hpp \
- ../../OpenFPM_data/src/util/for_each_ref.hpp \
- /usr/include/boost/fusion/include/size.hpp \
+ ../../OpenFPM_data/src/util/object_s_di.hpp \
+ ../../OpenFPM_data/src/util/object_si_d.hpp \
  ../../OpenFPM_devices/src/memory/ExtPreAlloc.hpp \
  Decomposition/CartDecomposition_unit_test.hpp \
  Decomposition/CartDecomposition.hpp
@@ -3372,7 +3377,7 @@ Graph/CartesianGraphFactory.hpp:
 
 /usr/include/boost/detail/call_traits.hpp:
 
-../../OpenFPM_data/src/ct_array.hpp:
+../../OpenFPM_data/src/util/ct_array.hpp:
 
 ../../OpenFPM_data/src/memory_array.hpp:
 
@@ -3414,10 +3419,10 @@ Graph/CartesianGraphFactory.hpp:
 
 ../../OpenFPM_data/src/Grid/grid_key_expression.hpp:
 
-../../OpenFPM_data/src/Space/Shape/Point.hpp:
-
 ../../OpenFPM_data/src/Grid/Encap.hpp:
 
+../../OpenFPM_data/src/Space/Shape/Point.hpp:
+
 ../../OpenFPM_data/src/Grid/grid_key.hpp:
 
 ../../OpenFPM_data/src/Grid/Encap.hpp:
@@ -3430,6 +3435,12 @@ Graph/CartesianGraphFactory.hpp:
 
 ../../OpenFPM_data/src/common.hpp:
 
+../../OpenFPM_data/src/util/object_s_di.hpp:
+
+../../OpenFPM_data/src/util/for_each_ref.hpp:
+
+/usr/include/boost/fusion/include/size.hpp:
+
 /home/i-bird/MPI/include/mpi.h:
 
 /home/i-bird/MPI/include/mpi_portable_platform.h:
@@ -4248,6 +4259,8 @@ Grid/grid_dist_key.hpp:
 
 ../../OpenFPM_data/src/base_type.hpp:
 
+../../OpenFPM_data/src/Point_orig.hpp:
+
 Decomposition/CartDecomposition.hpp:
 
 Decomposition/Decomposition.hpp:
@@ -4398,6 +4411,8 @@ dec_optimizer.hpp:
 
 ../../OpenFPM_data/src/NN/CellList/CellDecomposer.hpp:
 
+../../OpenFPM_data/src/Space/Matrix.hpp:
+
 /usr/include/c++/4.8.3/unordered_map:
 
 /usr/include/c++/4.8.3/bits/hashtable.h:
@@ -4408,10 +4423,10 @@ dec_optimizer.hpp:
 
 ../../OpenFPM_data/src/NN/CellList/CellList.hpp:
 
-../../OpenFPM_data/src/NN/CellList/CellListFast.hpp:
-
 ../../OpenFPM_data/src/NN/CellList/CellDecomposer.hpp:
 
+../../OpenFPM_data/src/NN/CellList/CellListFast.hpp:
+
 ../../OpenFPM_data/src/NN/CellList/CellNNIterator.hpp:
 
 ../../OpenFPM_data/src/NN/CellList/CellListBal.hpp:
@@ -4466,11 +4481,9 @@ Vector/vector_dist_key.hpp:
 
 ../../OpenFPM_data/src/util/object_creator.hpp:
 
-../../OpenFPM_data/src/util/object_copy.hpp:
+../../OpenFPM_data/src/util/object_s_di.hpp:
 
-../../OpenFPM_data/src/util/for_each_ref.hpp:
-
-/usr/include/boost/fusion/include/size.hpp:
+../../OpenFPM_data/src/util/object_si_d.hpp:
 
 ../../OpenFPM_devices/src/memory/ExtPreAlloc.hpp:
 
diff --git a/src/Decomposition/CartDecomposition.hpp b/src/Decomposition/CartDecomposition.hpp
index d6c197b495f6af3c6977ba246c68bb50e79ccf9b..0a4b56c5cf2acdc63c5cef2716dcd6a894a47d2c 100644
--- a/src/Decomposition/CartDecomposition.hpp
+++ b/src/Decomposition/CartDecomposition.hpp
@@ -41,7 +41,7 @@
  * \note sub-domain are the result of merging one or more sub-sub-domain (optimization)
  * \note Near processors are the processors adjacent to this processor
  * \note near processor sub-domain is a sub-domain that live in the a near (or contiguous) processor
- * \note external ghost box are the ghost space of the processors
+ * \note external ghost box (or ghost box) are the boxes that compose the ghost space of the processor
  * \note internal ghost box are the part of ghost of the near processor that intersect the space of the
  *       processor
  *
@@ -52,7 +52,7 @@ class CartDecomposition
 {
 	struct N_box
 	{
-		// id of the processor in the nn_processor list
+		// id of the processor in the nn_processor list (local processor id)
 		size_t id;
 
 		// Near processor sub-domains
@@ -276,6 +276,26 @@ private:
 
 			++it;
 		}
+
+		// Get the smallest sub-division on each direction
+		::Box<dim,T> unit = getSmallestSubdivision();
+		// Get the processor bounding Box
+		::Box<dim,T> bound = getProcessorBounds();
+
+		// calculate the sub-divisions (0.5 for rounding error)
+		size_t div[dim];
+		for (size_t i = 0 ; i < dim ; i++)
+			div[i] = (size_t)((bound.getHigh(i) - bound.getLow(i)) / unit.getHigh(i) + 0.5);
+
+		// Create shift
+		Point<dim,T> orig;
+
+		// p1 point of the Processor bound box is the shift
+		for (size_t i = 0 ; i < dim ; i++)
+			orig.get(i) = bound.getLow(i);
+
+		// Initialize the geo_cell structure
+		geo_cell.Initialize(domain,div,orig);
 	}
 
 	/*! \brief Create the subspaces that decompose your domain
@@ -343,13 +363,8 @@ private:
 		// cast the pointer
 		CartDecomposition<dim,T,device_l,Memory,Domain,data_s> * cd = static_cast< CartDecomposition<dim,T,device_l,Memory,Domain,data_s> *>(ptr);
 
-		if (cd->v_cl.getProcessUnitID() == 0)
-		{
-			std::cout << "Receiving from " << i << "       msg size: " << msg_i << "\n";
-		}
-
 		// Resize the memory
-		cd->nn_processor_subdomains[i].bx.resize(msg_i);
+		cd->nn_processor_subdomains[i].bx.resize(msg_i / sizeof(::Box<dim,T>) );
 
 		// Return the receive pointer
 		return cd->nn_processor_subdomains[i].bx.getPointer();
@@ -415,48 +430,145 @@ public:
 	~CartDecomposition()
 	{}
 
+	// It store all the boxes of the near processors in a linear array
+	struct p_box
+	{
+		//! Box that identify the intersection of the ghost of the near processor with the
+		//! processor sub-domain
+		::Box<dim,T> box;
+		//! local processor id
+		size_t lc_proc;
+		//! processor id
+		size_t proc;
+
+		/*! \brief Check if two p_box are the same
+		 *
+		 * \param pb box to check
+		 *
+		 */
+		bool operator==(const p_box & pb)
+		{
+			return pb.lc_proc == lc_proc;
+		}
+	};
+
 	openfpm::vector<size_t> ids;
 
+	/*! \brief class to select the returned id by ghost_processorID
+	 *
+	 */
+	class box_id
+	{
+	public:
+		/*! \brief Return the box id
+		 *
+		 * \param p structure containing the id informations
+		 * \param b_id box_id
+		 *
+		 * \return box id
+		 *
+		 */
+		inline static size_t id(p_box & p, size_t b_id)
+		{
+			return b_id;
+		}
+	};
+
+	/*! \brief class to select the returned id by ghost_processorID
+	 *
+	 */
+	class processor_id
+	{
+	public:
+		/*! \brief Return the processor id
+		 *
+		 * \param p structure containing the id informations
+		 * \param b_id box_id
+		 *
+		 * \return processor id
+		 *
+		 */
+		inline static size_t id(p_box & p, size_t b_id)
+		{
+			return p.proc;
+		}
+	};
+
+	/*! \brief class to select the returned id by ghost_processorID
+	 *
+	 */
+	class lc_processor_id
+	{
+	public:
+		/*! \brief Return the near processor id
+		 *
+		 * \param p structure containing the id informations
+		 * \param b_id box_id
+		 *
+		 * \return local processor id
+		 *
+		 */
+		inline static size_t id(p_box & p, size_t b_id)
+		{
+			return p.lc_proc;
+		}
+	};
+
+
+#define UNIQUE 1
+#define MULTIPLE 2
+
 	/*! \brief Given a position it return if the position belong to any neighborhood processor ghost
 	 * (Internal ghost)
 	 *
+	 * \tparam id type of if to get box_id processor_id lc_processor_id
 	 * \param p Particle position
+	 * \param opt intersection boxes of the same processor can overlap, so in general the function
+	 *        can produce more entry with the same processor, the UNIQUE option eliminate double entries
+	 *        (UNIQUE) is for particle data (MULTIPLE) is for grid data [default MULTIPLE]
 	 *
 	 * \param return the processor ids
 	 *
 	 */
-	inline const openfpm::vector<size_t> ghost_processorID(Point<dim,T> & p)
+	template <typename id> inline const openfpm::vector<size_t> ghost_processorID(Point<dim,T> & p, const int opt = MULTIPLE)
 	{
 		ids.clear();
 
-		// Check with geo-cell if a particle is inside one Cell caotaining boxes
+		// Check with geo-cell if a particle is inside one Cell containing boxes
 
-		auto cell_it = geo_cell.getCellIterator(p);
+		auto cell_it = geo_cell.getIterator(geo_cell.getCell(p));
 
 		// For each element in the cell, check if the point is inside the box
-		// if it is store the processor id
+		// if it is, store the processor id
 		while (cell_it.isNext())
 		{
 			size_t bid = cell_it.get();
 
 			if (vb_int.get(bid).box.isInside(p) == true)
 			{
-				ids.add(vb_int.get(bid).proc);
+				ids.add(id::id(vb_int.get(bid),bid));
 			}
+
+			++cell_it;
 		}
 
+		// Make the id unique
+		if (opt == UNIQUE)
+			ids.unique();
+
 		return ids;
 	}
 
 	/*! \brief Given a position it return if the position belong to any neighborhood processor ghost
 	 * (Internal ghost)
 	 *
+	 * \tparam id type of if to get box_id processor_id lc_processor_id
 	 * \param p Particle position
 	 *
 	 * \param return the processor ids
 	 *
 	 */
-	template<typename Mem> inline const openfpm::vector<size_t> ghost_processorID(const encapc<1,Point<dim,T>,Mem> & p)
+	template<typename id, typename Mem> inline const openfpm::vector<size_t> ghost_processorID(const encapc<1,Point<dim,T>,Mem> & p, const int opt = MULTIPLE)
 	{
 		ids.clear();
 
@@ -472,21 +584,23 @@ public:
 
 			if (vb_int.get(bid).box.isInside(p) == true)
 			{
-				ids.add(vb_int.get(bid).proc);
+				ids.add(id::id(vb_int.get(bid),bid));
 			}
+
+			++cell_it;
 		}
 
+		// Make the id unique
+		if (opt == UNIQUE)
+			ids.unique();
+
 		return ids;
 	}
 
-	// It store all the boxes of the near processors in a linear array
-	struct p_box
-	{
-		::Box<dim,T> box;
-		size_t proc;
-	};
+	// External ghost boxes for this processor, indicated with G8_0 G9_0 ...
+	openfpm::vector<p_box> vb_ext;
 
-	// Internal boxes for this processor domain, indicated with B8_0 B9_0 ..... in the figure
+	// Internal ghost boxes for this processor domain, indicated with B8_0 B9_0 ..... in the figure
 	// below as a linear vector
 	openfpm::vector<p_box> vb_int;
 
@@ -674,7 +788,16 @@ p1[0]<-----+         +----> p2[0]
 					bool intersect = sub_with_ghost.Intersect(::Box<dim,T>(p_box.get(b)),bi);
 
 					if (intersect == true)
+					{
+						struct p_box pb;
+
+						pb.box = bi;
+						pb.proc = p_id;
+						pb.lc_proc = ProctoID(p_id);
+
+						vb_ext.add(pb);
 						p_box_int.add(bi);
+					}
 				}
 			}
 
@@ -687,6 +810,9 @@ p1[0]<-----+         +----> p2[0]
 				// get the set of sub-domains of the contiguous processor p_id
 				openfpm::vector< ::Box<dim,T> > & nn_p_box = nn_processor_subdomains[p_id].bx;
 
+				// get the local processor id
+				size_t lc_proc = nn_processor_subdomains[p_id].id;
+
 				// near processor sub-domain intersections
 				openfpm::vector< ::Box<dim,T> > & p_box_int = box_nn_processor_int.get(i).get(j).nbx;
 
@@ -703,13 +829,18 @@ p1[0]<-----+         +----> p2[0]
 					n_sub.enlarge(ghost);
 
 					// Intersect with the local sub-domain
-
 					p_box b_int;
 					bool intersect = n_sub.Intersect(l_sub,b_int.box);
 
 					// store if it intersect
 					if (intersect == true)
 					{
+						// fill with the processor id
+						b_int.proc = p_id;
+
+						// fill the local processor id
+						b_int.lc_proc = lc_proc;
+
 						p_box_int.add(b_int.box);
 						vb_int.add(b_int);
 
@@ -784,7 +915,7 @@ p1[0]<-----+         +----> p2[0]
 	 *
 	 */
 
-	size_t inline processorID(T (&p)[dim])
+	size_t inline processorID(const T (&p)[dim]) const
 	{
 		return fine_s.get(cd.getCell(p));
 	}
@@ -973,9 +1104,9 @@ p1[0]<-----+         +----> p2[0]
 	 * \return true if it is local
 	 *
 	 */
-	template<typename Mem> bool isLocal(encapc<1, Point<dim,T>, Mem> p)
+	template<typename Mem> bool isLocal(const encapc<1, Point<dim,T>, Mem> p) const
 	{
-		return processorID<Mem>() == v_cl.getProcessUnitID();
+		return processorID<Mem>(p) == v_cl.getProcessUnitID();
 	}
 
 	/*! \brief Check if the particle is local
@@ -985,7 +1116,7 @@ p1[0]<-----+         +----> p2[0]
 	 * \return true if it is local
 	 *
 	 */
-	bool isLocal(T (&pos)[dim])
+	bool isLocal(const T (&pos)[dim]) const
 	{
 		return processorID(pos) == v_cl.getProcessUnitID();
 	}
@@ -1014,49 +1145,70 @@ p1[0]<-----+         +----> p2[0]
 		return geo_cell.getIterator(geo_cell.getCell(p));
 	}
 
-	/*! \brief if the point fall into the ghost of some near processor it return the processor number in which
-	 *  it fall
+
+	inline size_t getNNProcessors() const
+	{
+		return nn_processors.size();
+	}
+
+	/*! \brief Get the number of the calculated internal ghost boxes
 	 *
-	 * \param p Point
-	 * \return number of processors
+	 * \return the number of internal ghost boxes
+	 *
+	 */
+	inline size_t getNIGhostBox() const
+	{
+		return vb_int.size();
+	}
+
+	/*! \brief Give the internal ghost box id, it return at which processor it belong
+	 *
+	 * \return the processor id
 	 *
 	 */
-/*	inline size_t labelPointNp(Point<dim,T> & p)
+	inline ::Box<dim,T> getIGhostBox(size_t b_id) const
 	{
-		return geo_cell.getNelements(geo_cell.getCell(p));
-	}*/
+		return vb_int.get(b_id).box;
+	}
 
-	/*! \brief It return the label point cell
+	/*! \brief Give the processor id of the internal ghost box
 	 *
-	 * The labeling of a point p is regulated by a Cell list, give a point it give a cell-id
+	 * \return the processor id of the ghost box
 	 *
-	 * \param p Point
-	 * \return cell-id
+	 */
+	inline size_t getIGhostBoxProcessor(size_t b_id) const
+	{
+		return vb_int.get(b_id).proc;
+	}
+
+	/*! \brief Get the number of the calculated ghost boxes
+	 *
+	 * \return the number of internal ghost boxes
 	 *
 	 */
-/*	inline size_t labelPointCell(Point<dim,T> & p)
+	inline size_t getNGhostBox() const
 	{
-		return geo_cell.getCell(p);
-	}*/
+		return vb_ext.size();
+	}
 
-	/*! \brief get the number of near processors
+	/*! \brief Given the ghost box id, it return at which processor it belong
 	 *
-	 * \return the number of near processors
+	 * \return the processor id
 	 *
 	 */
-	inline size_t getNNProcessors()
+	inline ::Box<dim,T> getGhostBox(size_t b_id) const
 	{
-		return nn_processors.size();
+		return vb_ext.get(b_id).box;
 	}
 
-	/*! \brief Give the internal ghost box id, it return at which processor it belong
+	/*! \brief Give the processor id of the internal ghost box
 	 *
-	 * \return the number of near processors
+	 * \return the processor id of the ghost box
 	 *
 	 */
-	inline size_t getGhostBoxProcessor(size_t b_id)
+	inline size_t getGhostBoxProcessor(size_t b_id) const
 	{
-		return vb_int.get(b_id).proc;
+		return vb_ext.get(b_id).proc;
 	}
 
 	/*! \brief Convert the near processor ID to processor number
@@ -1070,6 +1222,18 @@ p1[0]<-----+         +----> p2[0]
 	{
 		return nn_processors.get(id);
 	}
+
+	/*! \brief Convert the processor id to local processor id
+	 *
+	 * \param processor id
+	 *
+	 * \return the local processor id
+	 *
+	 */
+	inline size_t ProctoID(size_t p)
+	{
+		return nn_processor_subdomains[p].id;
+	}
 };
 
 
diff --git a/src/Decomposition/CartDecomposition_unit_test.hpp b/src/Decomposition/CartDecomposition_unit_test.hpp
index 0fd1c9f388f22f7691d044d675e09939313ce1f5..eb2b351f948c7ac13bd3d536b7656b14a73d1209 100644
--- a/src/Decomposition/CartDecomposition_unit_test.hpp
+++ b/src/Decomposition/CartDecomposition_unit_test.hpp
@@ -39,7 +39,36 @@ BOOST_AUTO_TEST_CASE( CartDecomposition_test_use)
 	// create a ghost border
 	dec.calculateGhostBoxes(g);
 
-	//
+	// For each calculated ghost box
+	for (size_t i = 0 ; i < dec.getNIGhostBox() ; i++)
+	{
+		SpaceBox<3,float> b = dec.getIGhostBox(i);
+		size_t proc = dec.getIGhostBoxProcessor(i);
+
+		// sample one point inside the box
+		Point<3,float> p = b.rnd();
+
+		// Check that ghost_processorsID return that processor number
+		const openfpm::vector<size_t> & pr = dec.template ghost_processorID<CartDecomposition<3,float>::processor_id>(p);
+
+		bool found = false;
+
+		for (size_t j = 0; j < pr.size() ; j++)
+		{
+			if (pr.get(j) == proc)
+			{found = true; break;}
+		}
+
+		if (found == false)
+		{
+			int debug = 0;
+			debug++;
+
+			const openfpm::vector<size_t> pr2 = dec.template ghost_processorID<CartDecomposition<3,float>::processor_id>(p);
+		}
+
+		BOOST_REQUIRE_EQUAL(found,true);
+	}
 }
 
 BOOST_AUTO_TEST_SUITE_END()
diff --git a/src/Vector/vector_dist.hpp b/src/Vector/vector_dist.hpp
index cf8c479fa2ff0d2740e3eb5c74dcda3b00109eb7..5def7088c3a3237d7800484054024c8a54f43e3e 100644
--- a/src/Vector/vector_dist.hpp
+++ b/src/Vector/vector_dist.hpp
@@ -27,7 +27,9 @@
 #define PUT 2
 
 #define INTERNAL 0
+
 #define NO_POSITION 1
+#define WITH_POSITION 2
 
 /*! \brief Distributed vector
  *
@@ -38,6 +40,9 @@ class vector_dist
 {
 private:
 
+	// Ghost marker, all the particle with id > g_m are ghost all with g_m < are real particle
+	size_t g_m = 0;
+
 	// indicate from where the ghost particle start in the vector
 	size_t ghost_pointer;
 
@@ -63,7 +68,7 @@ public:
 
 	/*! \brief Constructor
 	 *
-	 * \param number of elements
+	 * \param Global number of elements
 	 *
 	 */
 	vector_dist(size_t np, Box box)
@@ -73,6 +78,9 @@ public:
 		v_pos = v_cl.template allocate<openfpm::vector<point>>(1);
 		v_prp = v_cl.template allocate<openfpm::vector<prop>>(1);
 
+		// convert to a local number of elements
+		np /= v_cl.getProcessingUnits();
+
 		// resize the position vector
 		v_pos.get(0).resize(np);
 
@@ -137,9 +145,9 @@ public:
 	 * \param vec_key vector element
 	 *
 	 */
-	template<unsigned int id> inline auto getProp(vect_dist_key_dx vec_key) -> decltype(v_prp.get(vec_key.v_c).template get<id>(vec_key.key))
+	template<unsigned int id> inline auto getProp(vect_dist_key_dx vec_key) -> decltype(v_prp.get(vec_key.getSub()).template get<id>(vec_key.getKey()))
 	{
-		return v_prp.get(vec_key.v_c).template get<id>(vec_key.key);
+		return v_prp.get(vec_key.getSub()).template get<id>(vec_key.getKey());
 	}
 
 	/*! \brief It store for each processor the position and properties vector of the particles
@@ -164,6 +172,9 @@ public:
 		dec.calculateGhostBoxes(g);
 	}
 
+	//! It map the processor id with the communication request into map procedure
+	openfpm::vector<size_t> p_map_req;
+
 	/*! \brief It communicate the particle to the respective processor
 	 *
 	 */
@@ -212,6 +223,9 @@ public:
 			++it;
 		}
 
+		// resize the map
+		p_map_req.resize(v_cl.getProcessingUnits());
+
 		// Create the sz and prc buffer
 
 		openfpm::vector<size_t> prc_sz_r;
@@ -221,6 +235,7 @@ public:
 		{
 			if (p_map.get(i) == 1)
 			{
+				p_map_req.get(i) = prc_r.size();
 				prc_r.add(i);
 				prc_sz_r.add(prc_sz.get(i));
 			}
@@ -233,8 +248,8 @@ public:
 		for (size_t i = 0 ;  i < prc_r.size() ; i++)
 		{
 			// Create the size required to store the particles position and properties to communicate
-			size_t s1 = openfpm::vector<point>::calculateMem(prc_sz_r.get(i),0);
-			size_t s2 = openfpm::vector<prop>::calculateMem(prc_sz_r.get(i),0);
+			size_t s1 = openfpm::vector<point,openfpm::device_cpu<point>,HeapMemory,openfpm::grow_policy_identity>::calculateMem(prc_sz_r.get(i),0);
+			size_t s2 = openfpm::vector<prop,openfpm::device_cpu<prop>,HeapMemory,openfpm::grow_policy_identity>::calculateMem(prc_sz_r.get(i),0);
 
 			// Preallocate the memory
 			size_t sz[2] = {s1,s2};
@@ -268,7 +283,7 @@ public:
 				continue;
 			}
 
-			lbl = (lbl > v_cl.getProcessUnitID())?lbl-1:lbl;
+			lbl = p_map_req.get(lbl);
 
 			pb.get(lbl).pos.set(prc_cnt.get(lbl),v_pos.get(up_v).get(key));
 			pb.get(lbl).prp.set(prc_cnt.get(lbl),v_prp.get(up_v).get(key));
@@ -295,7 +310,7 @@ public:
 		// Send and receive the particles
 
 		recv_cnt = 0;
-		v_cl.sendrecvMultipleMessagesPCX(prc_sz_r.size(),&p_map.get(0), &prc_sz_r.get(0), &prc_r.get(0) , &ptr.get(0) , vector_dist::message_alloc, this ,NEED_ALL_SIZE);
+		v_cl.sendrecvMultipleMessagesPCX(prc_sz_r.size(),&p_map.get(0), &prc_sz_r.get(0), &prc_r.get(0) , &ptr.get(0) , vector_dist::message_alloc_map, this ,NEED_ALL_SIZE);
 
 		// overwrite the outcoming particle with the incoming particle and resize the vectors
 
@@ -360,8 +375,6 @@ public:
 	// Memory for the ghost position sending buffer
 	Memory g_pos_mem;
 
-
-
 	/*! \brief It synchronize getting the ghost particles
 	 *
 	 * \prp Properties to get
@@ -369,11 +382,20 @@ public:
 	 * 		NO_RELABEL: If the particles does not move avoid to relabel and send particle position
 	 *
 	 */
-	template<int... prp> void ghost_get(size_t opt = NONE)
+	template<int... prp> void ghost_get(size_t opt = WITH_POSITION)
 	{
-		// Buffer that containe the number of elements for each processor
+		// Sending property object
+		typedef object<typename object_creator<typename prop::type,prp...>::type> prp_object;
+
+		// send vector for each processor
+		typedef  openfpm::vector<prp_object,openfpm::device_cpu<prp_object>,ExtPreAlloc<Memory>> send_vector;
+
+		// Buffer that contain the number of elements to send for each processor
 		ghost_prc_sz.clear();
 		ghost_prc_sz.resize(dec.getNNProcessors());
+		// Buffer that contain for each processor the id of the particle to send
+		opart.clear();
+		opart.resize(dec.getNNProcessors());
 
 		// Label the internal (assigned) particles
 		auto it = v_pos.get(INTERNAL).getIterator();
@@ -383,28 +405,22 @@ public:
 		{
 			auto key = it.get();
 
-			const openfpm::vector<size_t> & vp_id = dec.ghost_processorID(v_pos.get(INTERNAL).get(key));
+			const openfpm::vector<size_t> & vp_id = dec.template ghost_processorID<typename Decomposition::lc_processor_id>(v_pos.get(INTERNAL).get(key),UNIQUE);
 
 			for (size_t i = 0 ; i < vp_id.size() ; i++)
 			{
-				// Box id
-				size_t b_id = vp_id.get(i);
 				// processor id
-				size_t p_id = dec.getGhostBoxProcessor(b_id);
+				size_t p_id = vp_id.get(i);
 
 				// add particle to communicate
 				ghost_prc_sz.get(p_id)++;
 
-
 				opart.get(p_id).add(key);
 			}
 
 			++it;
 		}
 
-		// Sending property object
-		typedef object<typename object_creator<typename prop::type,prp...>::type> prp_object;
-
 		// Send buffer size in byte ( one buffer for all processors )
 		size_t size_byte_prp = 0;
 		size_t size_byte_pos = 0;
@@ -437,13 +453,10 @@ public:
 		// Create an object of preallocated memory for position
 		if (opt != NO_POSITION) prAlloc_pos = new ExtPreAlloc<Memory>(pap_pos,g_pos_mem);
 
-		// definition of the send vector for each processor
-		typedef  openfpm::vector<prp_object,openfpm::device_cpu<prp_object>,ExtPreAlloc<Memory>> send_vector;
-
 		// create a vector of send vector (ExtPreAlloc warrant that all the created vector are contiguous)
 		openfpm::vector<send_vector> g_send_prp;
 
-		// create a number of send buffer equal to the near processors
+		// create a number of send buffers equal to the near processors
 		g_send_prp.resize(ghost_prc_sz.size());
 		for (size_t i = 0 ; i < g_send_prp.size() ; i++)
 		{
@@ -465,7 +478,7 @@ public:
 				typedef encapc<1,prp_object,typename openfpm::vector<prp_object>::memory_t> encap_dst;
 
 				// Copy only the selected properties
-				object_copy<encap_src,encap_dst,ENCAP,prp...>(v_prp.get(INTERNAL).get(opart.get(i).get(j)),g_send_prp.get(i).get(j));
+				object_si_d<encap_src,encap_dst,ENCAP,prp...>(v_prp.get(INTERNAL).get(opart.get(i).get(j)),g_send_prp.get(i).get(j));
 			}
 		}
 
@@ -477,7 +490,7 @@ public:
 		openfpm::vector<send_pos_vector> g_pos_send;
 		if (opt != NO_POSITION)
 		{
-			// create a number of send buffer equal to the near processors
+			// create a number of send buffers equal to the near processors
 			g_pos_send.resize(ghost_prc_sz.size());
 			for (size_t i = 0 ; i < g_pos_send.size() ; i++)
 			{
@@ -506,11 +519,94 @@ public:
 			prc.add(dec.IDtoProc(i));
 		}
 
-		// Send receive the particles information
-//		v_cl.sendRecvMultipleMessageNBX();
+		// Send receive the particles properties information
+		v_cl.sendrecvMultipleMessagesNBX(prc,g_send_prp,msg_alloc_ghost_get,this);
+
+		// Mark the ghost part
+		g_m = v_prp.get(INTERNAL).size();
+
+		// Get the number of near processors
+		for (size_t i = 0 ; i < dec.getNNProcessors() ; i++)
+		{
+			// calculate the number of received elements
+			size_t n_ele = recv_sz.get(i) / sizeof(prp_object);
+
+			// add the received particles to the vector
+			PtrMemory * ptr1 = new PtrMemory(recv_mem_gg.get(i).getPointer(),recv_sz.get(i));
+
+			// create vector representation to a piece of memory already allocated
+			openfpm::vector<prp_object,openfpm::device_cpu<prp_object>,PtrMemory,openfpm::grow_policy_identity> v2;
+
+			v2.setMemory(*ptr1);
+
+			// resize with the number of elements
+			v2.resize(n_ele);
+
+			// Add the ghost particle
+			v_prp.get(INTERNAL).template add<prp_object,PtrMemory,openfpm::grow_policy_identity,prp...>(v2);
+		}
+
+		if (opt != NO_POSITION)
+		{
+			// Send receive the particles properties information
+			v_cl.sendrecvMultipleMessagesNBX(prc,g_pos_send,msg_alloc_ghost_get,this);
+
+			// Get the number of near processors
+			for (size_t i = 0 ; i < dec.getNNProcessors() ; i++)
+			{
+				// calculate the number of received elements
+				size_t n_ele = recv_sz.get(i) / sizeof(point);
+
+				// add the received particles to the vector
+				PtrMemory * ptr1 = new PtrMemory(recv_mem_gg.get(i).getPointer(),recv_sz.get(i));
+
+				// create vector representation to a piece of memory already allocated
+
+				openfpm::vector<point,openfpm::device_cpu<point>,PtrMemory,openfpm::grow_policy_identity> v2;
+
+				v2.setMemory(*ptr1);
+
+				// resize with the number of elements
+				v2.resize(n_ele);
+
+				// Add the ghost particle
+				v_pos.get(INTERNAL).template add<PtrMemory,openfpm::grow_policy_identity>(v2);
+			}
+		}
+	}
+
+	// Receiving size
+	openfpm::vector<size_t> recv_sz;
+
+	// Receiving buffer for particles ghost get
+	openfpm::vector<HeapMemory> recv_mem_gg;
+
+	/*! \brief Call-back to allocate buffer to receive incoming objects (particles)
+	 *
+	 * \param msg_i message size required to receive from i
+	 * \param total_msg message size to receive from all the processors
+	 * \param total_p the total number of processor want to communicate with you
+	 * \param i processor id
+	 * \param ri request id (it is an id that goes from 0 to total_p, and is unique
+	 *           every time message_alloc is called)
+	 * \param ptr void pointer parameter for additional data to pass to the call-back
+	 *
+	 */
+	static void * msg_alloc_ghost_get(size_t msg_i ,size_t total_msg, size_t total_p, size_t i, size_t ri, void * ptr)
+	{
+		vector_dist<point,prop,Box,Decomposition,Memory,with_id> * v = static_cast<vector_dist<point,prop,Box,Decomposition,Memory,with_id> *>(ptr);
+
+		v->recv_sz.resize(v->dec.getNNProcessors());
+		v->recv_mem_gg.resize(v->dec.getNNProcessors());
+
+		// Get the local processor id
+		size_t lc_id = v->dec.ProctoID(i);
 
+		// resize the receive buffer
+		v->recv_mem_gg.get(lc_id).resize(msg_i);
+		v->recv_sz.get(lc_id) = msg_i;
 
-		// add the received particles to the vector
+		return v->recv_mem_gg.get(lc_id).getPointer();
 	}
 
 	// Heap memory receiver
@@ -535,7 +631,7 @@ public:
 	 * \return the pointer where to store the message
 	 *
 	 */
-	static void * message_alloc(size_t msg_i ,size_t total_msg, size_t total_p, size_t i, size_t ri, void * ptr)
+	static void * message_alloc_map(size_t msg_i ,size_t total_msg, size_t total_p, size_t i, size_t ri, void * ptr)
 	{
 		// cast the pointer
 		vector_dist<point,prop,Box,Decomposition,Memory,with_id> * vd = static_cast<vector_dist<point,prop,Box,Decomposition,Memory,with_id> *>(ptr);
@@ -567,22 +663,42 @@ public:
 		return vector_dist_iterator<openfpm::vector<point>>(v_pos);
 	}
 
+	/*! \brief Get the iterator across the position of the ghost particles
+	 *
+	 * \return an iterator
+	 *
+	 */
+	vector_dist_iterator<openfpm::vector<point>> getGhostIterator()
+	{
+		return vector_dist_iterator<openfpm::vector<point>>(v_pos,g_m);
+	}
+
 	/*! \brief Get the iterator across the properties of the particles
 	 *
 	 * \return an iterator
 	 *
 	 */
-	vector_dist_iterator<openfpm::vector<point>> getPropIterator()
+	vector_dist_iterator<openfpm::vector<prop>> getPropIterator()
 	{
 		return vector_dist_iterator<openfpm::vector<prop>>(v_prp);
 	}
 
+	/*! \brief Get the iterator across the properties of the ghost particles
+	 *
+	 * \return an iterator
+	 *
+	 */
+	vector_dist_iterator<openfpm::vector<prop>> getGhostPropIterator()
+	{
+		return vector_dist_iterator<openfpm::vector<prop>>(v_prp,g_m);
+	}
+
 	/*! \brief Get the decomposition
 	 *
 	 * \return
 	 *
 	 */
-	Decomposition & getDecomposition()
+	const Decomposition & getDecomposition()
 	{
 		return dec;
 	}
diff --git a/src/Vector/vector_dist_iterator.hpp b/src/Vector/vector_dist_iterator.hpp
index 10863afe9ddf12e5df43a52892e9571767b7b24f..c3145b9438a5d522de5b5ae4db38880a9a696e99 100644
--- a/src/Vector/vector_dist_iterator.hpp
+++ b/src/Vector/vector_dist_iterator.hpp
@@ -27,11 +27,12 @@ class vector_dist_iterator
 
 	/*! \brief Constructor of the distributed grid
 	 *
-	 * \param gk std::vector of the local grid
+	 * \param gk the set of local vectors
+	 * \param offset iterator starting point
 	 *
 	 */
-	vector_dist_iterator(Vcluster_object_array<device_v> & gk)
-	:v_c(0),vList(gk),v_it(0)
+	vector_dist_iterator(Vcluster_object_array<device_v> & gk, size_t offset = 0)
+	:v_c(0),vList(gk),v_it(offset)
 	{
 	}
 
diff --git a/src/Vector/vector_dist_unit_test.hpp b/src/Vector/vector_dist_unit_test.hpp
index 951fe3af8026bca64c7bfb4e6247567f9daaa03a..ef41db84f3891180df1285491abeef85b18cf3a0 100644
--- a/src/Vector/vector_dist_unit_test.hpp
+++ b/src/Vector/vector_dist_unit_test.hpp
@@ -18,9 +18,11 @@ BOOST_AUTO_TEST_CASE( vector_dist_iterator_test_use )
 	typedef Point_test<float> p;
 	typedef Point<2,float> s;
 
+	Vcluster & v_cl = *global_v_cluster;
+
     // set the seed
 	// create the random generator engine
-	std::srand(global_v_cluster->getProcessUnitID());
+	std::srand(v_cl.getProcessUnitID());
     std::default_random_engine eg;
     std::uniform_real_distribution<float> ud(0.0f, 1.0f);
 
@@ -42,7 +44,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_iterator_test_use )
 	vd.map();
 
 	// Check if we have all the local particles
-
+	size_t cnt = 0;
 	auto & ct = vd.getDecomposition();
 	it = vd.getIterator();
 
@@ -50,10 +52,91 @@ BOOST_AUTO_TEST_CASE( vector_dist_iterator_test_use )
 	{
 		auto key = it.get();
 
-		// Check if local otherwise print the particle value
-
+		// Check if local
 		BOOST_REQUIRE_EQUAL(ct.isLocal(vd.template getPos<s::x>(key)),true);
 
+		cnt++;
+
+		++it;
+	}
+
+	//
+	v_cl.reduce(cnt);
+	v_cl.execute();
+	BOOST_REQUIRE_EQUAL(cnt,4096);
+}
+
+BOOST_AUTO_TEST_CASE( vector_dist_ghost )
+{
+	// Communication object
+	Vcluster & v_cl = *global_v_cluster;
+
+	typedef Point_test<float> p;
+	typedef Point<2,float> s;
+
+	Box<2,float> box({0.0,0.0},{1.0,1.0});
+	size_t g_div[]= {1000,1000};
+
+	// processor division on y direction
+	size_t point_div = g_div[1] / v_cl.getProcessingUnits();
+
+	// Create a grid info
+	grid_sm<2,void> g_info(g_div);
+
+	// Calculate the grid spacing
+	Point<2,float> spacing = box.getP2();
+	spacing = spacing / g_div;
+
+	// middle spacing
+	Point<2,float> m_spacing = spacing / 2;
+
+	// create a sub iterator
+	grid_key_dx<2> start(0,point_div * v_cl.getProcessUnitID());
+	grid_key_dx<2> stop(999,point_div * (v_cl.getProcessUnitID() + 1) - 1);
+	auto g_sub = g_info.getSubIterator(start,stop);
+
+	// Vector of particles
+	vector_dist<Point<2,float>, Point_test<float>, Box<2,float>, CartDecomposition<2,float> > vd(g_info.size(),box);
+
+	auto it = vd.getIterator();
+
+	while (it.isNext())
+	{
+		auto key_v = it.get();
+		auto key = g_sub.get();
+
+		// set the particle position
+
+		vd.template getPos<s::x>(key_v)[0] = key.get(0) * spacing[0] + m_spacing[0];
+		vd.template getPos<s::x>(key_v)[1] = key.get(1) * spacing[1] + m_spacing[1];
+
+		if (vd.template getPos<s::x>(key_v)[0] >= 1.0 || vd.template getPos<s::x>(key_v)[1] >= 1.0)
+		{
+			int debug = 0;
+			debug++;
+		}
+
+		++g_sub;
+		++it;
+	}
+
+	// redistribute the particles according to the decomposition
+	vd.map();
+
+	// Fill the scalar with the particle position
+	const auto & ct = vd.getDecomposition();
+	it = vd.getIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		// fill with the processor ID where these particle live
+		vd.template getProp<p::s>(key) = vd.getPos<s::x>(key)[0] + vd.getPos<s::x>(key)[1] * 16;
+		vd.template getProp<p::v>(key)[0] = v_cl.getProcessUnitID();
+		vd.template getProp<p::v>(key)[1] = v_cl.getProcessUnitID();
+		vd.template getProp<p::v>(key)[2] = v_cl.getProcessUnitID();
+
 		++it;
 	}
 
@@ -63,8 +146,58 @@ BOOST_AUTO_TEST_CASE( vector_dist_iterator_test_use )
 	vd.setGhost(g);
 
 	// do a ghost get
-
 	vd.template ghost_get<p::s,p::v>();
+
+	// Get the decomposition
+	const auto & dec = vd.getDecomposition();
+
+	// Get the ghost external boxes
+	openfpm::vector<size_t> vb(dec.getNGhostBox());
+
+	// Get the ghost iterator
+	auto g_it = vd.getGhostIterator();
+
+	// Check if the ghost particles contain the correct information
+	while (g_it.isNext())
+	{
+		auto key = g_it.get();
+
+		float x0 = vd.getPos<s::x>(key)[0];
+		float x1 = vd.getPos<s::x>(key)[1] * 16;
+		float scalar = vd.template getProp<p::s>(key);
+
+		// Check the received data
+		BOOST_REQUIRE_EQUAL(vd.getPos<s::x>(key)[0] + vd.getPos<s::x>(key)[1] * 16,vd.template getProp<p::s>(key));
+
+		bool is_in = false;
+		size_t b = 0;
+
+		// check if the received data is in one of the ghost boxes
+		for ( ; b < dec.getNGhostBox() ; b++)
+		{
+			if (dec.getGhostBox(b).isInside(vd.getPos<s::x>(key)) == true)
+			{is_in = true; break;}
+		}
+		BOOST_REQUIRE_EQUAL(is_in,true);
+
+		// Check that the particle come from the correct processor
+		BOOST_REQUIRE_EQUAL(vd.getProp<p::v>(key)[0],dec.getGhostBoxProcessor(b));
+
+		// Add
+		vb.get(b)++;
+
+		++g_it;
+	}
+
+    CellDecomposer_sm<2,float> cd(SpaceBox<2,float>(box),g_div,0);
+
+	for (size_t i = 0 ; i < vb.size() ; i++)
+	{
+		// Calculate how many particle should be in the box
+		size_t n_point = cd.getGridPoints(dec.getGhostBox(i)).getVolume();
+
+		BOOST_REQUIRE_EQUAL(n_point,vb.get(i));
+	}
 }
 
 BOOST_AUTO_TEST_SUITE_END()
diff --git a/src/dec_optimizer.hpp b/src/dec_optimizer.hpp
index c1998384b0f92da26b4841d6a7dcd8172f5e28c9..6b26f87a334735ce21ebd3bbef67021499ca3699 100644
--- a/src/dec_optimizer.hpp
+++ b/src/dec_optimizer.hpp
@@ -106,6 +106,68 @@ class dec_optimizer
 
 private:
 
+	/*! \brief Expand one wavefront
+	 *
+	 * \param v_w wavefronts
+	 * \param w_comb wavefront expansion combinations
+	 * \param d direction of expansion
+	 *
+	 */
+	void expand_one_wf(openfpm::vector<wavefront<dim>> & v_w, std::vector<comb<dim>> & w_comb , size_t d)
+	{
+		for (int j = 0 ; j < dim ; j++)
+		{
+			v_w.template get<wavefront<dim>::stop>(d)[j] = v_w.template get<wavefront<dim>::stop>(d)[j] + w_comb[d].c[j];
+			v_w.template get<wavefront<dim>::start>(d)[j] = v_w.template get<wavefront<dim>::start>(d)[j] + w_comb[d].c[j];
+		}
+	}
+
+
+	/*! \brief Adjust the other wavefronts
+	 *
+	 * \param d direction
+	 *
+	 */
+	void adjust_others_wf(openfpm::vector<wavefront<dim>> & v_w,  HyperCube<dim> & hyp, std::vector<comb<dim>> & w_comb, size_t d)
+	{
+		// expand the intersection of the wavefronts
+
+		std::vector<comb<dim>> q_comb = SubHyperCube<dim,dim-1>::getCombinations_R(w_comb[d],dim-2);
+
+		// Eliminate the w_comb[d] direction
+
+		for (int k = 0 ; k < q_comb.size() ; k++)
+		{
+			for (int j = 0 ; j < dim ; j++)
+			{
+				if (w_comb[d].c[j] != 0)
+				{
+					q_comb[k].c[j] = 0;
+				}
+			}
+		}
+
+		// for all the combinations
+		for (int j = 0 ; j < q_comb.size() ; j++)
+		{
+			size_t id = hyp.LinId(q_comb[j]);
+
+			// get the combination of the direction d
+
+			bool is_pos = hyp.isPositive(d);
+
+			// is positive, modify the stop point or the starting point
+
+			for (int s = 0 ; s < dim ; s++)
+			{
+				if (is_pos == true)
+				{v_w.template get<wavefront<dim>::stop>(id)[s] = v_w.template get<wavefront<dim>::stop>(id)[s] + w_comb[d].c[s];}
+				else
+				{v_w.template get<wavefront<dim>::start>(id)[s] = v_w.template get<wavefront<dim>::start>(id)[s] + w_comb[d].c[s];}
+			}
+		}
+	}
+
 	/* \brief Fill the wavefront position
 	 *
 	 * \tparam prp property to set
@@ -192,16 +254,13 @@ private:
 			}
 		}
 
-		// take the wavefront expand on direction d of one
+		// Create an Hyper-cube
+		HyperCube<dim> hyp;
 
 		for (int d = 0 ; d < v_w.size() ; d++)
 		{
-			// expand the wavefront
-			for (int j = 0 ; j < dim ; j++)
-			{
-				v_w.template get<wavefront<dim>::start>(d)[j] = v_w.template get<wavefront<dim>::start>(d)[j] + w_comb[d].c[j];
-				v_w.template get<wavefront<dim>::stop>(d)[j] = v_w.template get<wavefront<dim>::stop>(d)[j] + w_comb[d].c[j];
-			}
+			expand_one_wf(v_w,w_comb,d);
+			adjust_others_wf(v_w,hyp,w_comb,d);
 		}
 
 		// for each expanded wavefront create a sub-grid iterator and add the sub-domain