diff --git a/src/.deps/pdata-main.Po b/src/.deps/pdata-main.Po
index 01321ba088d4669fc2e33835e7db436c005e0b06..d2d126eb73dafc8573cf7700a17b3f27478b5b1d 100644
--- a/src/.deps/pdata-main.Po
+++ b/src/.deps/pdata-main.Po
@@ -1297,6 +1297,7 @@ pdata-main.o: main.cpp /usr/include/stdc-predef.h \
  Grid/grid_dist_key.hpp ../../OpenFPM_data/src/Point_test.hpp \
  ../../OpenFPM_data/src/base_type.hpp \
  ../../OpenFPM_data/src/Point_orig.hpp \
+ ../../OpenFPM_data/src/Grid/Encap.hpp \
  Decomposition/CartDecomposition.hpp Decomposition/Decomposition.hpp \
  ../../OpenFPM_data/src/global_const.hpp SubdomainGraphNodes.hpp \
  metis_util.hpp ../../metis_install/include/metis.h \
@@ -1396,6 +1397,8 @@ pdata-main.o: main.cpp /usr/include/stdc-predef.h \
  ../../OpenFPM_data/src/util/object_s_di.hpp \
  ../../OpenFPM_data/src/util/object_si_d.hpp \
  ../../OpenFPM_devices/src/memory/ExtPreAlloc.hpp \
+ ../../OpenFPM_IO/src/CSVWriter.hpp \
+ ../../OpenFPM_IO/src/csv_multiarray.hpp \
  Decomposition/CartDecomposition_unit_test.hpp \
  Decomposition/CartDecomposition.hpp
 
@@ -4261,6 +4264,8 @@ Grid/grid_dist_key.hpp:
 
 ../../OpenFPM_data/src/Point_orig.hpp:
 
+../../OpenFPM_data/src/Grid/Encap.hpp:
+
 Decomposition/CartDecomposition.hpp:
 
 Decomposition/Decomposition.hpp:
@@ -4487,6 +4492,10 @@ Vector/vector_dist_key.hpp:
 
 ../../OpenFPM_devices/src/memory/ExtPreAlloc.hpp:
 
+../../OpenFPM_IO/src/CSVWriter.hpp:
+
+../../OpenFPM_IO/src/csv_multiarray.hpp:
+
 Decomposition/CartDecomposition_unit_test.hpp:
 
 Decomposition/CartDecomposition.hpp:
diff --git a/src/Decomposition/CartDecomposition.hpp b/src/Decomposition/CartDecomposition.hpp
index 0a4b56c5cf2acdc63c5cef2716dcd6a894a47d2c..1ff5b8951af2126f8ddcfced9e5db0197ba720fe 100644
--- a/src/Decomposition/CartDecomposition.hpp
+++ b/src/Decomposition/CartDecomposition.hpp
@@ -253,14 +253,6 @@ private:
 			ss_box.contained(sub_d);
 		}
 
-		//++++++++++++++++++++++++++++++++++++++++ Debug output NN boxes
-		{
-		VTKWriter<openfpm::vector<::SpaceBox<dim,T>>,VECTOR_BOX> vtk_box1;
-		vtk_box1.add(sub_domains);
-		vtk_box1.write(std::string("loc_") + std::to_string(v_cl.getProcessUnitID()) + std::string(".vtk"));
-		}
-		//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
 		// fill fine_s structure
 		// fine_s structure contain the processor id for each sub-sub-domain
 		// with sub-sub-domain we mean the sub-domain decomposition before
@@ -720,41 +712,11 @@ p1[0]<-----+         +----> p2[0]
 			}
 		}
 
-		//++++++++++++++++++++++++++++++++++++++++ Debug output NN boxes
-		{
-		for (size_t b = 0 ; b < boxes.size() ; b++)
-		{
-			VTKWriter<openfpm::vector<::SpaceBox<dim,T>>,VECTOR_BOX> vtk_box1;
-			vtk_box1.add(boxes.get(b));
-			vtk_box1.write(std::string("Processor_") + std::to_string(v_cl.getProcessUnitID()) + "_" + std::to_string(nn_processors.get(b)) + std::string(".vtk"));
-		}
-		}
-
-		//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
 		// Intersect all the local sub-domains with the sub-domains of the contiguous processors
 
 		// Get the sub-domains of the near processors
 		v_cl.sendrecvMultipleMessagesNBX(nn_processors,boxes,CartDecomposition<dim,T,device_l,Memory,Domain,data_s>::message_alloc, this ,NEED_ALL_SIZE);
 
-		// ++++++++++++++++++++++++++++++++++++++++++ Check received boxes
-
-		{
-		VTKWriter<openfpm::vector<::Box<dim,T>>,VECTOR_BOX> vtk_box1;
-		for (size_t p = 0 ; p < nn_processors.size() ; p++)
-		{
-			size_t prc = nn_processors.get(p);
-
-			if (v_cl.getProcessUnitID() == 0)
-				std::cout << "Received from " << prc << "      n_boxes: " << nn_processor_subdomains[prc].bx.size() << "\n";
-
-			vtk_box1.add(nn_processor_subdomains[prc].bx);
-		}
-		vtk_box1.write(std::string("rb_Processor_") + std::to_string(v_cl.getProcessUnitID()) + std::string(".vtk"));
-		}
-
-		// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
 		box_nn_processor_int.resize(sub_domains.size());
 
 		// For each sub-domain
@@ -864,24 +826,6 @@ p1[0]<-----+         +----> p2[0]
 					}
 				}
 			}
-
-
-			// ++++++++++++++++++++++++++++++++++++++++ Debug +++++++++++++++++++++++++++++
-
-			{
-			VTKWriter<openfpm::vector<::Box<dim,T>>,VECTOR_BOX> vtk_box1;
-			for (size_t p = 0 ; p < box_nn_processor_int.size() ; p++)
-			{
-				for (size_t s = 0 ; s < box_nn_processor_int.get(p).size() ; s++)
-				{
-					vtk_box1.add(box_nn_processor_int.get(p).get(s).nbx);
-				}
-			}
-			vtk_box1.write(std::string("inte_Processor_") + std::to_string(v_cl.getProcessUnitID()) + std::string(".vtk"));
-			}
-
-			// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
 		}
 	}
 
@@ -1234,6 +1178,50 @@ p1[0]<-----+         +----> p2[0]
 	{
 		return nn_processor_subdomains[p].id;
 	}
+
+	/*! \brief Write the decomposition as VTK file
+	 *
+	 * The function generate several files
+	 *
+	 * 1) p_sub_X.vtk domain for the processor X as union of sub-domain
+	 * 2) sub_np_c_X.vtk sub-domain of the near processors contiguous to the processor X (Color encoded)
+	 * 3) sub_X_inte_g_np.vtk Intersection between the ghosts of the near processors and the processors X sub-domains (Color encoded)
+	 *
+	 * where X is the processor number
+	 *
+	 * \param output directory where to write the files
+	 *
+	 */
+	bool write(std::string output) const
+	{
+		//! p_sub_X.vtk domain for the processor X as union of sub-domain
+		VTKWriter<openfpm::vector<::SpaceBox<dim,T>>,VECTOR_BOX> vtk_box1;
+		vtk_box1.add(sub_domains);
+		vtk_box1.write(std::string("p_sub_") + std::to_string(v_cl.getProcessUnitID()) + std::string(".vtk"));
+
+		//! sub_np_c_X.vtk sub-domain of the near processors contiguous to the processor X (Color encoded)
+		VTKWriter<openfpm::vector<::Box<dim,T>>,VECTOR_BOX> vtk_box2;
+		for (size_t p = 0 ; p < nn_processors.size() ; p++)
+		{
+			size_t prc = nn_processors.get(p);
+			auto it = nn_processor_subdomains.find(prc);
+			if (it != nn_processor_subdomains.end())
+				vtk_box2.add(nn_processor_subdomains.at(prc).bx);
+		}
+		vtk_box2.write(std::string("sub_np_c_") + std::to_string(v_cl.getProcessUnitID()) + std::string(".vtk"));
+
+		//! sub_X_inte_g_np.vtk Intersection between the ghosts of the near processors and the processors X sub-domains (Color encoded)
+		VTKWriter<openfpm::vector<::Box<dim,T>>,VECTOR_BOX> vtk_box3;
+		for (size_t p = 0 ; p < box_nn_processor_int.size() ; p++)
+		{
+			for (size_t s = 0 ; s < box_nn_processor_int.get(p).size() ; s++)
+			{
+				auto & diocane = box_nn_processor_int.get(p).get(s).nbx;
+				vtk_box3.add(diocane);
+			}
+		}
+		vtk_box3.write(std::string("sub_") + std::to_string(v_cl.getProcessUnitID()) + std::string("_inte_g_np") + std::string(".vtk"));
+	}
 };
 
 
diff --git a/src/Graph/CartesianGraphFactory.hpp b/src/Graph/CartesianGraphFactory.hpp
index e8c2218d3095a6a881f85ba643bb951ddcd90b0d..9827982c1a41870d1c09dc4f4f09f04d776b947b 100644
--- a/src/Graph/CartesianGraphFactory.hpp
+++ b/src/Graph/CartesianGraphFactory.hpp
@@ -84,7 +84,7 @@ class Graph_constructor_impl
 {
 public:
 	//! Construct cartesian graph
-	static Graph construct(size_t (& sz)[dim], Box<dim,T> dom)
+	static Graph construct(const size_t (& sz)[dim], Box<dim,T> dom)
 	{
 		// Calculate the size of the hyper-cubes on each dimension
 
@@ -194,7 +194,7 @@ class Graph_constructor_impl<dim,Graph,NO_EDGE,T,dim_c,pos...>
 {
 public:
 	//! Construct cartesian graph
-	static Graph construct(size_t ( & sz)[dim], Box<dim,T> dom)
+	static Graph construct(const size_t ( & sz)[dim], Box<dim,T> dom)
 	{
 		// Calculate the size of the hyper-cubes on each dimension
 
@@ -369,7 +369,7 @@ public:
 	 *
 	 */
 	template <unsigned int se,typename T, unsigned int dim_c, int... pos>
-	static Graph construct(size_t (& sz)[dim], Box<dim,T> dom)
+	static Graph construct(const size_t (& sz)[dim], Box<dim,T> dom)
 	{
 		return Graph_constructor_impl<dim,Graph,se,T,dim_c,pos...>::construct(sz,dom);
 	}
diff --git a/src/Vector/vector_dist.hpp b/src/Vector/vector_dist.hpp
index 5def7088c3a3237d7800484054024c8a54f43e3e..e422e20f2c791b23be732d8efdb8434b6a356ba2 100644
--- a/src/Vector/vector_dist.hpp
+++ b/src/Vector/vector_dist.hpp
@@ -19,10 +19,12 @@
 #include "common.hpp"
 #include "util/object_util.hpp"
 #include "memory/ExtPreAlloc.hpp"
+#include "CSVWriter.hpp"
 
 #define NO_ID false
 #define ID true
 
+// Perform a ghost get or a ghost put
 #define GET	1
 #define PUT 2
 
@@ -31,6 +33,10 @@
 #define NO_POSITION 1
 #define WITH_POSITION 2
 
+// Write the particles with ghost
+#define NO_GHOST 0
+#define WITH_GHOST 2
+
 /*! \brief Distributed vector
  *
  */
@@ -314,6 +320,7 @@ public:
 
 		// overwrite the outcoming particle with the incoming particle and resize the vectors
 
+		size_t total_element = 0;
 		size_t o_p_id = 0;
 
 		for (size_t i = 0 ; i < v_proc.size() ; i++)
@@ -322,8 +329,13 @@ public:
 
 			size_t n_ele = v_proc.get(i) / (sizeof(point) + sizeof(prop));
 
-			PtrMemory * ptr1 = new PtrMemory(hp_recv.getPointer(),n_ele * sizeof(point));
-			PtrMemory * ptr2 = new PtrMemory((unsigned char *)hp_recv.getPointer() + n_ele * sizeof(point),n_ele * sizeof(prop));
+			// Pointer of the received positions for each near processor
+			void * ptr_pos = ((unsigned char *)hp_recv.getPointer()) + (total_element * (sizeof(point) + sizeof(prop)));
+			// Pointer of the received properties for each near processor
+			void * ptr_prp = ((unsigned char *)hp_recv.getPointer()) + (total_element * (sizeof(point) + sizeof(prop))) + n_ele * sizeof(point);
+
+			PtrMemory * ptr1 = new PtrMemory(ptr_pos,n_ele * sizeof(point));
+			PtrMemory * ptr2 = new PtrMemory(ptr_prp,n_ele * sizeof(prop));
 
 			// create vector representation to a piece of memory already allocated
 
@@ -352,6 +364,9 @@ public:
 				v_prp.get(0).add();
 				v_prp.get(0).set(v_prp.get(0).size()-1,vprp.get(j));
 			}
+
+			// increment the total number of element counter
+			total_element += n_ele;
 		}
 
 		// remove the hole (out-going particles) in the vector
@@ -525,8 +540,8 @@ public:
 		// 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++)
+		// Process the received data (recv_mem_gg != 0 if you have data)
+		for (size_t i = 0 ; i < dec.getNNProcessors() && recv_mem_gg.size() != 0 ; i++)
 		{
 			// calculate the number of received elements
 			size_t n_ele = recv_sz.get(i) / sizeof(prp_object);
@@ -543,7 +558,7 @@ public:
 			v2.resize(n_ele);
 
 			// Add the ghost particle
-			v_prp.get(INTERNAL).template add<prp_object,PtrMemory,openfpm::grow_policy_identity,prp...>(v2);
+			v_prp.get(INTERNAL).template add_prp<prp_object,PtrMemory,openfpm::grow_policy_identity,prp...>(v2);
 		}
 
 		if (opt != NO_POSITION)
@@ -551,8 +566,8 @@ public:
 			// 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++)
+			// Process the received data (recv_mem_gg != 0 if you have data)
+			for (size_t i = 0 ; i < dec.getNNProcessors() && recv_mem_gg.size() != 0 ; i++)
 			{
 				// calculate the number of received elements
 				size_t n_ele = recv_sz.get(i) / sizeof(point);
@@ -702,6 +717,30 @@ public:
 	{
 		return dec;
 	}
+
+	/*! \brief Output particle position and properties
+	 *
+	 * \param File output
+	 * \param opt NO_GHOST or WITH_GHOST
+	 *
+	 * \return if the file has been written correctly
+	 *
+	 */
+	inline bool write(std::string out, int opt = NO_GHOST)
+	{
+		if (hasEnding(out,".csv"))
+		{
+			// CSVWriter test
+			CSVWriter<openfpm::vector<point>, openfpm::vector<prop> > csv_writer;
+
+			std::string output = std::to_string(v_cl.getProcessUnitID()) + std::string("_") + out;
+
+			// Write the CSV
+			return csv_writer.write(output,v_pos.get(INTERNAL),v_prp.get(INTERNAL));
+		}
+
+		return false;
+	}
 };
 
 
diff --git a/src/Vector/vector_dist_iterator.hpp b/src/Vector/vector_dist_iterator.hpp
index c3145b9438a5d522de5b5ae4db38880a9a696e99..7ac4c34501cc158daf5ecce4805f0ef39cb38c93 100644
--- a/src/Vector/vector_dist_iterator.hpp
+++ b/src/Vector/vector_dist_iterator.hpp
@@ -34,6 +34,8 @@ class vector_dist_iterator
 	vector_dist_iterator(Vcluster_object_array<device_v> & gk, size_t offset = 0)
 	:v_c(0),vList(gk),v_it(offset)
 	{
+		if ( offset >= vList.get(0).size() )
+			v_c++;
 	}
 
 	// Destructor
diff --git a/src/Vector/vector_dist_unit_test.hpp b/src/Vector/vector_dist_unit_test.hpp
index ef41db84f3891180df1285491abeef85b18cf3a0..9276a8a72bdaab5ac376768f1f106895652c73c9 100644
--- a/src/Vector/vector_dist_unit_test.hpp
+++ b/src/Vector/vector_dist_unit_test.hpp
@@ -13,59 +13,6 @@
 
 BOOST_AUTO_TEST_SUITE( vector_dist_test )
 
-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(v_cl.getProcessUnitID());
-    std::default_random_engine eg;
-    std::uniform_real_distribution<float> ud(0.0f, 1.0f);
-
-	Box<2,float> box({0.0,0.0},{1.0,1.0});
-	vector_dist<Point<2,float>, Point_test<float>, Box<2,float>, CartDecomposition<2,float> > vd(4096,box);
-
-	auto it = vd.getIterator();
-
-	while (it.isNext())
-	{
-		auto key = it.get();
-
-		vd.template getPos<s::x>(key)[0] = ud(eg);
-		vd.template getPos<s::x>(key)[1] = ud(eg);
-
-		++it;
-	}
-
-	vd.map();
-
-	// Check if we have all the local particles
-	size_t cnt = 0;
-	auto & ct = vd.getDecomposition();
-	it = vd.getIterator();
-
-	while (it.isNext())
-	{
-		auto key = it.get();
-
-		// 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
@@ -75,7 +22,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost )
 	typedef Point<2,float> s;
 
 	Box<2,float> box({0.0,0.0},{1.0,1.0});
-	size_t g_div[]= {1000,1000};
+	size_t g_div[]= {16,16};
 
 	// processor division on y direction
 	size_t point_div = g_div[1] / v_cl.getProcessingUnits();
@@ -91,8 +38,8 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost )
 	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);
+	grid_key_dx<2> start(point_div * v_cl.getProcessUnitID(),0);
+	grid_key_dx<2> stop(point_div * (v_cl.getProcessUnitID() + 1) - 1,g_div[0]);
 	auto g_sub = g_info.getSubIterator(start,stop);
 
 	// Vector of particles
@@ -120,11 +67,18 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost )
 		++it;
 	}
 
+	// Debug write the particles
+	vd.write("Particles_before_map.csv");
+
 	// redistribute the particles according to the decomposition
 	vd.map();
 
+	// Debug write particles
+	vd.write("Particles_after_map.csv");
+
 	// Fill the scalar with the particle position
 	const auto & ct = vd.getDecomposition();
+
 	it = vd.getIterator();
 
 	while (it.isNext())
@@ -141,13 +95,19 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost )
 	}
 
 	// set the ghost based on the radius cut off
-	Ghost<2,float> g(0.01);
+	Ghost<2,float> g(spacing.get(0));
 
 	vd.setGhost(g);
 
+	//! Output the decomposition
+	ct.write(".");
+
 	// do a ghost get
 	vd.template ghost_get<p::s,p::v>();
 
+	// Debug write the particles with GHOST
+	vd.write("Particles_with_ghost.csv",WITH_GHOST);
+
 	// Get the decomposition
 	const auto & dec = vd.getDecomposition();
 
@@ -200,6 +160,59 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost )
 	}
 }
 
+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(v_cl.getProcessUnitID());
+    std::default_random_engine eg;
+    std::uniform_real_distribution<float> ud(0.0f, 1.0f);
+
+	Box<2,float> box({0.0,0.0},{1.0,1.0});
+	vector_dist<Point<2,float>, Point_test<float>, Box<2,float>, CartDecomposition<2,float> > vd(4096,box);
+
+	auto it = vd.getIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		vd.template getPos<s::x>(key)[0] = ud(eg);
+		vd.template getPos<s::x>(key)[1] = ud(eg);
+
+		++it;
+	}
+
+	vd.map();
+
+	// Check if we have all the local particles
+	size_t cnt = 0;
+	auto & ct = vd.getDecomposition();
+	it = vd.getIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		// 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_SUITE_END()
 
 #endif /* VECTOR_DIST_UNIT_TEST_HPP_ */
diff --git a/src/dec_optimizer.hpp b/src/dec_optimizer.hpp
index 6b26f87a334735ce21ebd3bbef67021499ca3699..5881558d5f59b2334b06f060e0219d559531ce9f 100644
--- a/src/dec_optimizer.hpp
+++ b/src/dec_optimizer.hpp
@@ -541,7 +541,7 @@ public:
 	 *
 	 */
 
-	dec_optimizer(Graph & g, size_t (& sz)[dim])
+	dec_optimizer(Graph & g, const size_t (& sz)[dim])
 	:gh(sz)
 	{
 		// The graph g is suppose to represent a cartesian grid
diff --git a/src/main.cpp b/src/main.cpp
index cb799f429a7c6701b823e7c86b009ab904f4a773..bdee86089aaed30fc6d36b449445bb34fbc4c7f0 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,5 +1,7 @@
 #include <iostream>
 #include "config.h"
+
+#define NO_WARNING
 #include "Graph/CartesianGraphFactory.hpp"
 
 #define BOOST_DISABLE_ASSERTS
@@ -8,13 +10,6 @@
 #define BOOST_TEST_MODULE "C++ test module for OpenFPM_pdata project"
 #include <boost/test/included/unit_test.hpp>
 
-/*struct MPIFixture {
-    MPIFixture() { MPI_Init(&boost::unit_test::framework::master_test_suite().argc,&boost::unit_test::framework::master_test_suite().argv);}
-    ~MPIFixture() { MPI_Finalize(); }
-};
-
-BOOST_GLOBAL_FIXTURE(MPIFixture);*/
-
 #include "Grid/grid_dist_id.hpp"
 #include "Point_test.hpp"
 #include "Decomposition/CartDecomposition.hpp"