diff --git a/CHANGELOG.md b/CHANGELOG.md
index fa176d72caa6e413dfae2d4b3d235b46d689bf4d..c8ca27d71ea95a3649d1ca76b85fccd1d413ab33 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,6 +1,15 @@
 # Change Log
 All notable changes to this project will be documented in this file.
 
+## [0.9.0]
+
+### Added
+- Introduced getDomainIterator for Cell-list
+
+### Fixed
+- Installation of PETSC in case with MUMPS try without MUMPS
+- In case of miss compilation ignore system wide installation
+
 ## [0.8.0] 28 February 2016
 
 ### Added
diff --git a/openfpm_data b/openfpm_data
index 1152e735515f48cd9d6734f2d1a7b6d6dd0a4b91..dde50897665d36fc2723aabb0917291c9573c065 160000
--- a/openfpm_data
+++ b/openfpm_data
@@ -1 +1 @@
-Subproject commit 1152e735515f48cd9d6734f2d1a7b6d6dd0a4b91
+Subproject commit dde50897665d36fc2723aabb0917291c9573c065
diff --git a/src/Decomposition/CartDecomposition.hpp b/src/Decomposition/CartDecomposition.hpp
index cf3c91a1b194fc0115bbd9baa2ef267889e4bb0a..2ea88a1234e956a855dc043891e384de3e2d2e3f 100755
--- a/src/Decomposition/CartDecomposition.hpp
+++ b/src/Decomposition/CartDecomposition.hpp
@@ -1288,7 +1288,22 @@ public:
 		return domain_nn_calculator_cart<dim>::getDomainCells(shift,cell_shift,gs,proc_box,loc_box);
 	}
 
-	/*! \brief Get the anomalous cells
+	/*! \brief Get the CRS domain Cells
+	 *
+	 * It perform a linearization of the domain cells using the extension provided in gs
+	 *
+	 *
+	 * \param shift Cell padding
+	 * \param cell_shift where the domain cell start
+	 * \param gs grid extension
+	 *
+	 */
+	openfpm::vector<size_t> & getCRSDomainCells(grid_key_dx<dim> & shift, grid_key_dx<dim> & cell_shift, grid_sm<dim,void> & gs)
+	{
+		return domain_nn_calculator_cart<dim>::getCRSDomainCells(shift,cell_shift,gs,proc_box,loc_box);
+	}
+
+	/*! \brief Get the CRS anomalous cells
 	 *
 	 * This function include also a linearization of the indexes
 	 *
@@ -1301,9 +1316,9 @@ public:
 	 * \return the anomalous cells with neighborhood
 	 *
 	 */
-	openfpm::vector<subsub_lin<dim>> & getAnomDomainCells(grid_key_dx<dim> & shift, grid_key_dx<dim> & cell_shift, grid_sm<dim,void> & gs)
+	openfpm::vector<subsub_lin<dim>> & getCRSAnomDomainCells(grid_key_dx<dim> & shift, grid_key_dx<dim> & cell_shift, grid_sm<dim,void> & gs)
 	{
-		return domain_nn_calculator_cart<dim>::getAnomDomainCells(shift,cell_shift,gs,proc_box,loc_box);
+		return domain_nn_calculator_cart<dim>::getCRSAnomDomainCells(shift,cell_shift,gs,proc_box,loc_box);
 	}
 
 	/*! \brief Check if the particle is local considering boundary conditions
diff --git a/src/Decomposition/Distribution/SpaceDistribution.hpp b/src/Decomposition/Distribution/SpaceDistribution.hpp
index c7b2d9ac1197c1bf963a257bc32fd00897a0ffe4..346008d5acb726620f0e1ee5a3d827f7b064a39e 100644
--- a/src/Decomposition/Distribution/SpaceDistribution.hpp
+++ b/src/Decomposition/Distribution/SpaceDistribution.hpp
@@ -371,6 +371,18 @@ public:
 
 		return *this;
 	}
+
+	/*! \brief It return the decomposition id
+	 *
+	 * It just return 0
+	 *
+	 * \return 0
+	 *
+	 */
+	size_t get_ndec()
+	{
+		return 0;
+	}
 };
 
 
diff --git a/src/Decomposition/Domain_NN_calculator_cart.hpp b/src/Decomposition/Domain_NN_calculator_cart.hpp
index b4571fad5afe836c81673677b96ef18be3d6b8e2..e0889bc25ffe0f236c6d81530ca51ed6334db061 100644
--- a/src/Decomposition/Domain_NN_calculator_cart.hpp
+++ b/src/Decomposition/Domain_NN_calculator_cart.hpp
@@ -20,18 +20,24 @@ class domain_nn_calculator_cart
 	//! True if domain and anomalous domain cells are computed
 	bool are_domain_anom_computed;
 
-	//! Are linearized the domain cell
+	//! Are linearized the domain cells
+    bool are_dom_cells_lin;
+
+	//! Are linearized the CRS domain cell
     bool are_dom_lin;
 
-    //! are linearized the anomalous cells
+    //! are linearized the CRS anomalous cells
     bool are_anom_lin;
 
-	//! anomalous cell neighborhood
+	//! anomalous cell neighborhood for CRS
 	openfpm::vector<subsub<dim>> anom;
 
-	//! Set of normal domain cells
+	//! Set of normal domain cells for CRS
 	openfpm::vector<grid_key_dx<dim>> dom;
 
+	//! Set of domain cells
+	openfpm::vector<grid_key_dx<dim>> dom_cells;
+
 	//! Linearization is calculated out of a shift and grid dimension this is the shift
 	grid_key_dx<dim> shift_calc_dom;
 
@@ -45,12 +51,15 @@ class domain_nn_calculator_cart
 	grid_sm<dim,void> gs_calc_anom;
 
 
-	//! Set of normal domain cells linearized
+	//! Set of normal CRS domain cells linearized
 	openfpm::vector<size_t> dom_lin;
 
-	//! Set of normal domain cells linearized
+	//! Set of anomalous CRS domain cells linearized
 	openfpm::vector<subsub_lin<dim>> anom_lin;
 
+	//! Set of linearized domain cells
+	openfpm::vector<size_t> dom_cells_lin;
+
 	//! Processor box
 	Box<dim,long int> proc_box;
 
@@ -103,6 +112,7 @@ class domain_nn_calculator_cart
 	 */
 	void CalculateDomAndAnomCells(openfpm::vector<subsub<dim>> & sub_keys,
 			                      openfpm::vector<grid_key_dx<dim>> & dom_subsub,
+								  openfpm::vector<grid_key_dx<dim>> & dom_cells,
 								  const ::Box<dim,long int> & proc_box,
 								  grid_key_dx<dim> & shift,
 								  const openfpm::vector<::Box<dim, size_t>> & loc_box)
@@ -149,6 +159,8 @@ class domain_nn_calculator_cart
 					g.template get<0>(src).add(dst + shift);
 				}
 
+				dom_cells.add(key + shift);
+
 				++sub;
 			}
 		}
@@ -183,7 +195,7 @@ class domain_nn_calculator_cart
 public:
 
 	domain_nn_calculator_cart()
-	:are_domain_anom_computed(false),are_dom_lin(false),are_anom_lin(false)
+	:are_domain_anom_computed(false),are_dom_cells_lin(false),are_dom_lin(false),are_anom_lin(false)
 	{}
 
 	/*! \brief Get the domain Cells
@@ -200,7 +212,39 @@ public:
 	{
 		if (are_domain_anom_computed == false)
 		{
-			CalculateDomAndAnomCells(anom,dom,proc_box,shift,loc_box);
+			CalculateDomAndAnomCells(anom,dom,dom_cells,proc_box,shift,loc_box);
+			are_domain_anom_computed = true;
+		}
+
+		if (are_dom_cells_lin == false)
+		{
+			dom_cells_lin.clear();
+			shift_calc_dom = shift;
+			gs_calc_dom = gs;
+			for (size_t i = 0 ; i < dom_cells.size() ; i++)
+				dom_cells_lin.add(gs.LinId(dom_cells.get(i) - cell_shift));
+
+			are_dom_cells_lin = true;
+		}
+
+		return dom_cells_lin;
+	}
+
+	/*! \brief Get the domain Cells
+	 *
+	 * \param shift Shifting point
+	 * \param gs grid extension
+	 * \param proc_box processor bounding box
+	 * \param loc_box set of local sub-domains
+	 *
+	 * \return The set of domain cells
+	 *
+	 */
+	openfpm::vector<size_t> & getCRSDomainCells(grid_key_dx<dim> & shift, grid_key_dx<dim> & cell_shift, grid_sm<dim,void> & gs, Box<dim,size_t> & proc_box, openfpm::vector<::Box<dim, size_t>> & loc_box)
+	{
+		if (are_domain_anom_computed == false)
+		{
+			CalculateDomAndAnomCells(anom,dom,dom_cells,proc_box,shift,loc_box);
 			are_domain_anom_computed = true;
 		}
 
@@ -228,12 +272,12 @@ public:
 	 * \return The set of anomalous cells
 	 *
 	 */
-	openfpm::vector<subsub_lin<dim>> & getAnomDomainCells(grid_key_dx<dim> & shift, grid_key_dx<dim> & cell_shift, grid_sm<dim,void> & gs, Box<dim,size_t> & proc_box, openfpm::vector<::Box<dim, size_t>> & loc_box)
+	openfpm::vector<subsub_lin<dim>> & getCRSAnomDomainCells(grid_key_dx<dim> & shift, grid_key_dx<dim> & cell_shift, grid_sm<dim,void> & gs, Box<dim,size_t> & proc_box, openfpm::vector<::Box<dim, size_t>> & loc_box)
 	{
 		// if the neighborhood of each sub-sub-domains has not been calculated, calculate it
 		if (are_domain_anom_computed == false)
 		{
-			CalculateDomAndAnomCells(anom,dom,proc_box,shift,loc_box);
+			CalculateDomAndAnomCells(anom,dom,dom_cells,proc_box,shift,loc_box);
 			are_domain_anom_computed = true;
 		}
 
diff --git a/src/Vector/vector_dist.hpp b/src/Vector/vector_dist.hpp
index 4ddfc2cac16206b04ec639f68a6a656c3ddd079a..8a715c74c2988d8a8bb6825947604ecebf78ffa4 100644
--- a/src/Vector/vector_dist.hpp
+++ b/src/Vector/vector_dist.hpp
@@ -938,7 +938,9 @@ public:
 
 		grid_sm<dim,void> gs = NN.getInternalGrid();
 
-		ver.createVerletCrs(r_cut,g_m,v_pos,getDecomposition().getDomainCells(shift,cell_shift,gs),getDecomposition().getAnomDomainCells(shift,cell_shift,gs));
+		ver.createVerletCrs(r_cut,g_m,v_pos,
+				            getDecomposition().getCRSDomainCells(shift,cell_shift,gs),
+							getDecomposition().getCRSAnomDomainCells(shift,cell_shift,gs));
 
 		ver.set_ndec(getDecomposition().get_ndec());
 
@@ -1031,7 +1033,9 @@ public:
 
 				grid_sm<dim,void> gs = NN.getInternalGrid();
 
-				ver.updateCrs(getDecomposition().getDomain(),r_cut,v_pos,g_m,getDecomposition().getDomainCells(shift,cell_shift,gs),getDecomposition().getAnomDomainCells(shift,cell_shift,gs));
+				ver.updateCrs(getDecomposition().getDomain(),r_cut,v_pos,g_m,
+						      getDecomposition().getCRSDomainCells(shift,cell_shift,gs),
+							  getDecomposition().getCRSAnomDomainCells(shift,cell_shift,gs));
 			}
 			else
 			{
@@ -1061,60 +1065,6 @@ public:
 		}
 	}
 
-#if 0
-
-/*	void getVerletDeprecated(openfpm::vector<openfpm::vector<size_t>> & verlet, St r_cut)
-	{
-		// resize verlet to store the number of particles
-		verlet.resize(size_local());
-
-		// get the cell-list
-		auto cl = getCellList(r_cut);
-
-		// square of the cutting radius
-		St r_cut2 = r_cut * r_cut;
-
-		// iterate the particles
-		auto it_p = this->getDomainIterator();
-		while (it_p.isNext())
-		{
-			// key
-			vect_dist_key_dx key = it_p.get();
-
-			// Get the position of the particles
-			Point<dim, St> p = this->getPos(key);
-
-			// Clear the neighborhood of the particle
-			verlet.get(key.getKey()).clear();
-
-			// Get the neighborhood of the particle
-			auto NN = cl.template getNNIterator<NO_CHECK>(cl.getCell(p));
-			while (NN.isNext())
-			{
-				auto nnp = NN.get();
-
-				// p != q
-				if (nnp == key.getKey())
-				{
-					++NN;
-					continue;
-				}
-
-				Point<dim, St> q = this->getPos(nnp);
-
-				if (p.distance2(q) < r_cut2)
-					verlet.get(key.getKey()).add(nnp);
-
-				// Next particle
-				++NN;
-			}
-
-			// next particle
-			++it_p;
-		}
-	}*/
-
-#endif
 
 	/*! \brief Construct a cell list starting from the stored particles and reorder a vector according to the Hilberts curve
 	 *
@@ -1322,6 +1272,32 @@ public:
 		return vector_dist_iterator(g_m, v_pos.size());
 	}
 
+	/*! \brief Get an iterator that traverse the particles in the domain
+	 *
+	 * \return an iterator
+	 *
+	 */
+	template<typename CellList> ParticleIt_Cells<dim,CellList> getDomainIteratorCells(CellList & NN)
+	{
+#ifdef SE_CLASS3
+		se3.getIterator();
+#endif
+
+		// Shift
+		grid_key_dx<dim> cell_shift = NN.getShift();
+
+		// Shift
+		grid_key_dx<dim> shift = NN.getShift();
+
+		// Add padding
+		for (size_t i = 0 ; i < dim ; i++)
+			shift.set_d(i,shift.get(i) + NN.getPadding(i));
+
+		grid_sm<dim,void> gs = NN.getInternalGrid();
+
+		return ParticleIt_Cells<dim,CellList>(NN,getDecomposition().getDomainCells(shift,cell_shift,gs));
+	}
+
 	/*! \brief Get an iterator that traverse the particles in the domain
 	 *
 	 * \return an iterator
@@ -1768,7 +1744,9 @@ public:
 		grid_sm<dim,void> gs = NN.getInternalGrid();
 
 		// First we check that
-		return ParticleItCRS_Cells<dim,cli>(NN,getDecomposition().getDomainCells(shift,cell_shift,gs),getDecomposition().getAnomDomainCells(shift,cell_shift,gs),NN.getNNc_sym());
+		return ParticleItCRS_Cells<dim,cli>(NN,getDecomposition().getCRSDomainCells(shift,cell_shift,gs),
+				                               getDecomposition().getCRSAnomDomainCells(shift,cell_shift,gs),
+											   NN.getNNc_sym());
 	}
 
 	/*! \brief Return from which cell we have to start in case of CRS interation
diff --git a/src/Vector/vector_dist_NN_tests.hpp b/src/Vector/vector_dist_NN_tests.hpp
index a750bc95b6bfc5e915ef7e9f46268a9a8b96b5d0..c893c741e2dba2bf097e70dce8fd722cb232dae9 100644
--- a/src/Vector/vector_dist_NN_tests.hpp
+++ b/src/Vector/vector_dist_NN_tests.hpp
@@ -217,5 +217,90 @@ BOOST_AUTO_TEST_CASE( vector_dist_full_NN )
 	}
 }
 
+BOOST_AUTO_TEST_CASE( vector_dist_particle_iteration )
+{
+	Vcluster & v_cl = create_vcluster();
+
+	if (v_cl.getProcessingUnits() > 12)
+		return;
+
+    // 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);
+
+    long int k = 750 * v_cl.getProcessingUnits();
+
+	print_test("Testing 3D particle cell iterator=",k);
+	BOOST_TEST_CHECKPOINT( "Testing 3D full NN search k=" << k );
+
+	Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
+
+	// Boundary conditions
+	size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+
+	float r_cut = 0.01;
+
+	// ghost
+	Ghost<3,float> ghost(r_cut);
+
+	typedef  aggregate<float> part_prop;
+
+	// Distributed vector
+	vector_dist<3,float, part_prop > vd(k,box,bc,ghost,BIND_DEC_TO_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);
+
+			// Fill some properties randomly
+
+			vd.getProp<0>(key) = 0.0;
+
+			++it;
+		}
+
+	vd.map();
+
+	// sync the ghost
+	vd.ghost_get<0>();
+
+	openfpm::vector<size_t> ids;
+	ids.resize(vd.size_local());
+
+	auto NN = vd.getCellListSym(r_cut);
+
+	auto it_pcell = vd.getDomainIteratorCells(NN);
+
+	size_t count = 0;
+	while (it_pcell.isNext())
+	{
+		count++;
+
+		size_t id = it_pcell.get();
+		ids.get(id) = 1;
+
+		BOOST_REQUIRE(id < vd.size_local());
+
+		++it_pcell;
+	}
+
+	v_cl.sum(count);
+	v_cl.execute();
+
+	for (size_t i = 0 ; i < ids.size() ; i++)
+	{
+		std::cout << "OUTPUT: " << i << "    " << ids.get(i) <<  create_vcluster().getProcessUnitID() << std::endl;
+	}
+
+	BOOST_REQUIRE_EQUAL((long int)count,k);
+}
 
 #endif /* SRC_VECTOR_VECTOR_DIST_NN_TESTS_HPP_ */