diff --git a/src/Decomposition/CartDecomposition.cpp b/src/Decomposition/CartDecomposition.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5833b60dc05d89538435e1b868e89c8e0b3b1113
--- /dev/null
+++ b/src/Decomposition/CartDecomposition.cpp
@@ -0,0 +1,132 @@
+/*
+ * CartDecomposition.cpp
+ *
+ * Created on: Aug 15, 2014
+ * Author: Pietro Incardona
+ */
+
+#include "CartDecomposition.hpp"
+
+
+
+/*! \brief The the bulk part of the data set, or the data that does not depend
+ * from the ghosts layers
+ *
+ * The the bulk part of the data set, or the data that does not depend from the
+ * ghosts layers
+ *
+ */
+
+/*template<typename T> T CartDecomposition<T>::getBulk(T data)
+{
+ // for each element in data
+
+ for (size_t i = 0; i < data.size() ; i++)
+ {
+ if (localSpace.isInside())
+ }
+
+}
+
+template<typename T> T CartDecomposition<T>::getInternal()
+{
+
+}*/
+
+/*! \brief Check if is border or bulk
+ *
+ * \param neighboorhood define the neighboorhood of all the points
+ * \return true if border, false if bulk
+ *
+ */
+
+bool borderOrBulk(neighborhood & nb)
+{
+ device::grid<1,size_t> nbr = nb.next();
+
+ // check the neighborhood
+
+ // get neighborhood iterator
+
+ grid_key_dx_iterator<dim> iterator_nbr = nbr.getIterator();
+
+ while (iterator_nbr.hasNext())
+ {
+ grid_key_dx key_nbr = iterator_nbr.next();
+
+ // check if the neighboorhood is internal
+
+ if(subspace.isBound(data.template get<Point::x>(key_nbr)) == false)
+ {
+ // it is border
+
+ return true;
+
+ ret.bord.push_back(key);
+ break;
+ }
+ }
+
+ return false;
+}
+
+/*! \brief This function divide the data set into bulk, border, external and internal part
+ *
+ * \tparam dim dimensionality of the structure storing your data
+ * (example if they are in 3D grid, has to be 3)
+ * \tparam T type of object we are dividing
+ * \tparam device type of layout selected
+ * \param data 1-dimensional grid of point
+ * \param nb define the neighborhood of all the points
+ * \return a structure with the set of objects divided
+ *
+ */
+
+template<unsigned int dim, typename T, template<typename> class layout, typename Memory, template<unsigned int, typename> class Domain, template<typename, typename, typename> class data_s>
+dataDiv<T> CartDecomposition<dim,T,layout>::divide(device::grid<1,Point<dim,T>> & data, neighborhood & nb)
+{
+ //! allocate the 3 subset
+
+ dataDiv<T> ret;
+
+ ret.bord = new boost::shared_ptr<T>(new T());
+ ret.inte = new boost::shared_ptr<T>(new T());
+ ret.ext = new boost::shared_ptr<T>(new T());
+
+ //! get grid iterator
+
+ grid_key_dx_iterator<dim> iterator = data.getIterator();
+
+ //! we iterate trough all the set of objects
+
+ while (iterator.hasNext())
+ {
+ grid_key_dx<dim> key = iterator.next();
+
+ //! Check if the object is inside the subspace
+
+ if (subspace.isBound(data.template get<Point<3,T>::x>(key)))
+ {
+ //! Check if the neighborhood is inside the subspace
+
+ if (borderOrBulk(nb) == true)
+ {
+ // It is border
+
+ ret.bord.push_back(key);
+ }
+ else
+ {
+ // It is bulk
+
+ ret.bulk.push_back(key);
+ }
+ }
+ else
+ {
+ //! it is external
+
+ ret.ext.push_back(key);
+ }
+ }
+}
diff --git a/src/Decomposition/CartDecomposition.hpp b/src/Decomposition/CartDecomposition.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a3361c0d02036f056ef3f21bbe1bce00be795c75
--- /dev/null
+++ b/src/Decomposition/CartDecomposition.hpp
@@ -0,0 +1,512 @@
+/*
+ * CartDecomposition.hpp
+ *
+ * Created on: Aug 15, 2014
+ * Author: Pietro Incardona
+ */
+
+#ifndef CARTDECOMPOSITION_HPP
+#define CARTDECOMPOSITION_HPP
+
+#include "config.h"
+#include "Decomposition.hpp"
+#include "map_vector.hpp"
+#include <vector>
+#include "global_const.hpp"
+#include <initializer_list>
+#include "map_vector.hpp"
+#include "SubdomainGraphNodes.hpp"
+#include "metis_util.hpp"
+#include "dec_optimizer.hpp"
+#include "Space/Shape/Box.hpp"
+
+/**
+ * \brief This class decompose a space into subspaces
+ *
+ * This class decompose a space into regular hyper-cube subspaces, and give the possibilities to
+ * select one subspace
+ *
+ * \tparam dim is the dimensionality of the physical domain we are going to decompose.
+ * \tparam T type of the space we decompose, Real, Integer, Complex ...
+ * \tparam layout to use
+ * \tparam Memory Memory factory used to allocate memory
+ * \tparam Domain Structure that contain the information of your physical domain
+ * \tparam data type of structure that store the sub-domain decomposition can be an openfpm structure like
+ * vector, ...
+ *
+ * \note if PARALLEL_DECOMPOSITION macro is defined a parallel decomposition algorithm is used, basically
+ * each processor does not recompute the same decomposition
+ *
+ */
+
+template<unsigned int dim, typename T, template<typename> class device_l=openfpm::device_cpu, typename Memory=HeapMemory, template<unsigned int, typename> class Domain=Box, template<typename, typename, typename, typename> class data_s = openfpm::vector>
+class CartDecomposition
+{
+public:
+ //! Type of the domain we are going to decompose
+ typedef T domain_type;
+
+ //! It simplify to access the SpaceBox element
+ typedef SpaceBox<dim,T> Box;
+
+private:
+
+ //! This is the access_key to data_s, for example in the case of vector
+ //! acc_key is size_t
+ typedef typename data_s<SpaceBox<dim,T>,device_l<SpaceBox<dim,T>>,Memory,openfpm::vector_grow_policy_default>::access_key acc_key;
+
+ //! Subspace selected
+ //! access_key in case of grid is just the set of the index to access the grid
+ std::vector<acc_key> id_sub;
+
+ //! the margin of the sub-domain selected
+ SpaceBox<dim,T> sub_domain;
+
+ //! the set of all local sub-domain as vector
+ data_s<SpaceBox<dim,T>,device_l<SpaceBox<dim,T>>,Memory,openfpm::vector_grow_policy_default> sub_domains;
+
+ //! number of total sub-domain
+ size_t N_tot;
+
+ //! number of sub-domain on each dimension
+ size_t div[dim];
+
+ //! rectangular domain to decompose
+ Domain<dim,T> domain;
+
+ //! Box Spacing
+ T spacing[dim];
+
+ //! Runtime virtual cluster machine
+ Vcluster & v_cl;
+
+ /*! \brief Create internally the decomposition
+ *
+ * \param v_cl Virtual cluster, used internally to handle or pipeline communication
+ *
+ */
+ void CreateDecomposition(Vcluster & v_cl)
+ {
+ // Calculate the total number of box and and the spacing
+ // on each direction
+
+ N_tot = 1;
+
+ // Get the box containing the domain
+ SpaceBox<dim,T> bs = domain.getBox();
+
+ for (unsigned int i = 0; i < dim ; i++)
+ {
+ // Calculate the spacing
+
+ spacing[i] = (bs.getHigh(i) - bs.getLow(i)) / div[i];
+ N_tot *= div[i];
+ }
+
+ // Here we use METIS
+
+ // Create a cartesian grid graph
+ CartesianGraphFactory<3,Graph_CSR<nm_part_v,nm_part_e>> g_factory_part;
+
+ // Processor graph
+
+ Graph_CSR<nm_part_v,nm_part_e> gp = g_factory_part.construct<NO_EDGE,float,2>(div,domain);
+
+ // Get the number of processing units
+ size_t Np = v_cl.getProcessingUnits();
+
+ // Get the processor id
+ long int p_id = v_cl.getProcessUnitID();
+
+ // Convert the graph to metis
+ Metis<Graph_CSR<nm_part_v,nm_part_e>> met(gp,Np);
+
+ // decompose
+
+ met.decompose<nm_part_v::id>();
+
+ // Optimize the decomposition creating bigger spaces
+ // And reducing Ghost over-stress
+
+ dec_optimizer<3,Graph_CSR<nm_part_v,nm_part_e>> d_o(gp,div);
+
+ // set of Boxes produced by the decomposition optimizer
+ openfpm::vector<::Box<dim,size_t>> loc_box;
+
+ grid_key_dx<3> keyZero(0,0,0);
+ d_o.optimize<nm_part_v::sub_id,nm_part_v::id>(keyZero,gp,p_id,loc_box);
+
+ // convert into sub-domain
+ for (size_t s = 0 ; s < loc_box.size() ; s++)
+ {
+ SpaceBox<dim,T> sub_d(loc_box.get(s));
+
+ // re-scale with the spacing
+ sub_d.mul(spacing);
+
+ // add the sub-domain
+ sub_domains.add(sub_d);
+ }
+ }
+
+ /*! \brief Create the subspaces that decompose your domain
+ *
+ * Create the subspaces that decompose your domain
+ *
+ */
+
+ void CreateSubspaces()
+ {
+ // Create a grid where each point is a space
+
+ grid<3,void> g(div);
+
+ // create a grid_key_dx iterator
+
+ grid_key_dx_iterator<dim> gk_it(g);
+
+ // Divide the space into subspaces
+
+ while (gk_it.isNext())
+ {
+ //! iterate through all subspaces
+ grid_key_dx<dim> key = gk_it.get();
+
+ //! Create a new subspace
+
+ SpaceBox<dim,T> tmp;
+
+ //! fill with the Margin of the box
+
+ for (int i = 0 ; i < dim ; i++)
+ {
+ tmp.setHigh(i,(key.get(i)+1)*spacing[i]);
+ tmp.setLow(i,key.get(i)*spacing[i]);
+ }
+
+ //! add the space box
+
+ sub_domains.add(tmp);
+
+ // add the iterator
+
+ ++gk_it;
+ }
+ }
+
+public:
+
+ /*! \brief Cartesian decomposition copy constructor
+ *
+ * \param v_cl Virtual cluster, used internally to handle or pipeline communication
+ *
+ */
+ CartDecomposition(CartDecomposition<dim,T,device_l,Memory,Domain,data_s> && cd)
+ :sub_domain(cd.sub_domain),N_tot(cd.N_tot),domain(cd.domain),v_cl(cd.v_cl)
+ {
+ //! Subspace selected
+ //! access_key in case of grid is just the set of the index to access the grid
+ id_sub.swap(cd.id_sub);
+
+ //! the set of all local sub-domain as vector
+ sub_domains.swap(cd.sub_domains);
+
+ for (int i = 0 ; i < dim ; i++)
+ {
+ this->div[i] = div[dim];
+
+ //! Box Spacing
+ this->spacing[i] = spacing[i];
+ }
+ }
+
+ /*! \brief Cartesian decomposition constructor
+ *
+ * \param v_cl Virtual cluster, used internally to handle or pipeline communication
+ *
+ */
+ CartDecomposition(Vcluster & v_cl)
+ :id_sub(0),N_tot(0),v_cl(v_cl)
+ {}
+
+ /*! \brief Cartesian decomposition constructor, it divide the space in boxes
+ *
+ * \param dec is a vector that store how to divide on each dimension
+ * \param domain is the domain to divide
+ * \param v_cl are information of the cluster runtime machine
+ *
+ */
+ CartDecomposition(std::vector<size_t> dec, Domain<dim,T> domain, Vcluster & v_cl)
+ :id_sub(0),div(dec),domain(domain),v_cl(v_cl)
+ {
+ // Create the decomposition
+
+ CreateDecomposition(v_cl);
+ }
+
+ //! Cartesian decomposition destructor
+ ~CartDecomposition()
+ {}
+
+ /*! \brief Set the parameter of the decomposition
+ *
+ * \param div_ std::vector storing into how many domain to decompose on each dimension
+ * \param domain_ domain to decompose
+ *
+ */
+ void setParameters(std::vector<size_t> div_, Domain<dim,T> domain_)
+ {
+ // Set the decomposition parameters
+
+ div = div_;
+ domain = domain_;
+
+ //! Create the decomposition
+
+ CreateDecomposition(v_cl);
+ }
+
+ /*! \brief Set the parameter of the decomposition
+ *
+ * \param div_ std::vector storing into how many domain to decompose on each dimension
+ * \param domain_ domain to decompose
+ *
+ */
+ void setParameters(size_t div_[dim], Domain<dim,T> domain_)
+ {
+ // Set the decomposition parameters
+
+ for (int i = 0 ; i < dim ; i++)
+ div[i] = div_[i];
+
+ domain = domain_;
+
+ //! Create the decomposition
+
+ CreateDecomposition(v_cl);
+ }
+
+ /*! \brief Get the number of local local hyper-cubes or sub-domains
+ *
+ * \return the number of sub-domains
+ *
+ */
+ size_t getNLocalHyperCube()
+ {
+ return sub_domains.size();
+ }
+
+ /*! The the bulk part of the data set, or the data that
+ * does not depend from the ghosts layers
+ *
+ * \return the bulk of your data
+ *
+ */
+ T getBulk()
+ {
+
+ }
+
+ /*! \brief This function divide the data set into bulk, border, external and internal part
+ *
+ * \tparam dim dimensionality of the structure storing your data
+ * (example if they are in 3D grid, has to be 3)
+ * \tparam T type of object we are dividing
+ * \tparam device type of layout selected
+ * \param data 1-dimensional grid of point
+ * \param nb define the neighborhood of all the points
+ * \return a structure with the set of objects divided
+ *
+ */
+
+// dataDiv<T> CartDecomposition<dim,T,layout>::divide(layout::grid<1,Point<dim,T>> & data, neighborhood & nb);
+
+ /*! The the internal part of the data set, or the data that
+ * are inside the local space
+ *
+ * \return the internal part of your data
+ *
+ */
+ T getInternal()
+ {
+
+ }
+
+ /*! Get the internal part of the dataset, or the data that
+ * depend from the ghost layers
+ *
+ * \return the ghost part of your data
+ *
+ */
+
+ T getBorder()
+ {
+
+ }
+
+ /*! Get the external part of the dataset, or the data that
+ * are outside localSpace including ghost
+ *
+ * \return the external part of your data
+ *
+ */
+ T getExternal()
+ {
+
+ }
+
+ /*! \brief Get the number of one set of hyper-cube enclosing one particular
+ * subspace, the hyper-cube enclose your space, even if one box is enough
+ * can be more that one to increase occupancy
+ *
+ * In case of Cartesian decomposition it just return 1, each subspace
+ * has one hyper-cube, and occupancy 1
+ *
+ * \param id of the subspace
+ * \return the number of hyper-cube enclosing your space
+ *
+ */
+ size_t getNHyperCube(size_t id)
+ {
+ return 1;
+ }
+
+ /*! \brief Get the hyper-cube margins id_c has to be 0
+ *
+ * Get the hyper-cube margins id_c has to be 0, each subspace
+ * has one hyper-cube
+ *
+ * \param id of the subspace
+ * \param id_c
+ * \return The specified hyper-cube space
+ *
+ */
+ SpaceBox<dim,T> & getHyperCubeMargins(size_t id, size_t id_c)
+ {
+#ifdef DEBUG
+ // Check if this subspace exist
+ if (id >= N_tot)
+ {
+ std::cerr << "Error CartDecomposition: id > N_tot";
+ }
+ else if (id_c > 0)
+ {
+ // Each subspace is an hyper-cube so return error if id_c > 0
+ std::cerr << "Error CartDecomposition: id_c > 0";
+ }
+#endif
+
+ return sub_domains.get<Object>(id);
+ }
+
+ /*! \brief Get the total number of Hyper-cube
+ *
+ * Get the total number of Hyper-cube
+ *
+ * \return The total number of hyper-cube
+ *
+ */
+
+ size_t getNHyperCube()
+ {
+ return N_tot;
+ }
+
+ /*! \brief produce an hyper-cube approximation of the space decomposition
+ *
+ */
+
+ void hyperCube()
+ {
+ }
+
+ /*! \brief Select the local space
+ *
+ * Select the local space
+ *
+ * \param sub select the sub-space
+ *
+ */
+ void setSpace(size_t sub)
+ {
+ id_sub.push_back(sub);
+ }
+
+
+ /*! \brief Get the local grids
+ *
+ * Get the local grids
+ *
+ * \return the local grids
+ *
+ */
+
+ auto getLocalHyperCubes() -> decltype(sub_domains) &
+ {
+ return sub_domains;
+ }
+
+ /*! \brief Get the local hyper-cubes
+ *
+ * Get the local hyper-cubes
+ *
+ * \param lc is the id of the space
+ * \return the local hyper-cube
+ *
+ */
+
+ SpaceBox<dim,T> getLocalHyperCube(size_t lc)
+ {
+ // Create a space box
+ SpaceBox<dim,T> sp;
+
+ // fill the space box
+
+ for (size_t k = 0 ; k < dim ; k++)
+ {
+ // create the SpaceBox Low and High
+ sp.setLow(k,sub_domains.template get<Box::p1>(lc)[k]);
+ sp.setHigh(k,sub_domains.template get<Box::p2>(lc)[k]);
+ }
+
+ return sp;
+ }
+
+ /*! \brief Return the structure that store the physical domain
+ *
+ * Return the structure that store the physical domain
+ *
+ * \return The physical domain
+ *
+ */
+
+ Domain<dim,T> & getDomain()
+ {
+ return domain;
+ }
+
+ /*! \brief It return a graph that represent the domain decomposed into the cartesian grid
+ *
+ * It return a graph that represent the domain decomposed into the cartesian grid
+ *
+ */
+
+/* Graph<> createGraphModel()
+ {
+
+ }*/
+
+ /*! \brief It return a graph that represent the domain decomposed into the cartesian grid
+ *
+ * It return a graph that represent the domain decomposed into the cartesian grid
+ *
+ *
+ */
+/* Graph<> createLocalGraphMode()
+ {
+
+ }*/
+};
+
+
+#endif
diff --git a/src/Decomposition/CartDecomposition_unit_test.hpp b/src/Decomposition/CartDecomposition_unit_test.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..851cd64e0039e2ab9ab56487f231b5efa454f303
--- /dev/null
+++ b/src/Decomposition/CartDecomposition_unit_test.hpp
@@ -0,0 +1,40 @@
+#ifndef CARTDECOMPOSITION_UNIT_TEST_HPP
+#define CARTDECOMPOSITION_UNIT_TEST_HPP
+
+#include "CartDecomposition.hpp"
+#include "mathutil.hpp"
+
+BOOST_AUTO_TEST_SUITE( CartDecomposition_test )
+
+#define SUB_UNIT_FACTOR 64
+
+BOOST_AUTO_TEST_CASE( CartDecomposition_test_use)
+{
+ // Virtual cluster
+ Vcluster vcl(&boost::unit_test::framework::master_test_suite().argc,&boost::unit_test::framework::master_test_suite().argv);
+
+ CartDecomposition<3,float> dec(vcl);
+
+ // Physical domain
+ Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
+ size_t div[3];
+
+ // Get the number of processor and calculate the number of sub-domain
+ // for decomposition
+ size_t n_proc = vcl.getProcessingUnits();
+ size_t n_sub = n_proc * SUB_UNIT_FACTOR;
+
+ // Calculate the number of sub-domain on each dimension
+ for (int i = 0 ; i < 3 ; i++)
+ {div[i] = round_big_2(pow(n_sub,1.0/3));}
+
+ // Decompose
+ dec.setParameters(div,box);
+
+
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+
+#endif
diff --git a/src/Decomposition/Decomposition.cpp b/src/Decomposition/Decomposition.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/src/Decomposition/Decomposition.hpp b/src/Decomposition/Decomposition.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..35afa31f3ac1064c40c49208d0032826e60c920d
--- /dev/null
+++ b/src/Decomposition/Decomposition.hpp
@@ -0,0 +1,71 @@
+#include "Space/SpaceBox.hpp"
+
+#ifndef DECOMPOSITION_HPP_
+#define DECOMPOSITION_HPP_
+
+/**
+ *
+ * \brief class that store Internal part external and border part of a dataset
+ *
+ */
+
+template <typename T> class dataDiv
+{
+ //! Border part of the data
+ boost::shared_ptr<T> bord;
+ //! internal part of your data
+ boost::shared_ptr<T> inte;
+ //! external part of your data
+ boost::shared_ptr<T> ext;
+};
+
+/*! \brief This class define the domain decomposition interface
+ *
+ * This class define the domain decomposition interface, its main functionality
+ * is to divide a domain in several subspaces
+ *
+ * \tparam T structure that store the dataset
+ * \tparam S type of space is decomposing Real integer complex ...
+ *
+ */
+
+template<typename T, typename S>
+class Decomposition
+{
+ /*! \brief The the internal part of the data set, or the data that
+ * does not depend from the ghosts layers
+ *
+ * \return The internal part of the dataset
+ *
+ */
+ virtual T getInternal();
+
+
+ /*! Get the ghost part of the dataset
+ *
+ * \return The internal part of the dataset
+ *
+ */
+ virtual T getBorder();
+
+ /*! Get the external part of the dataset (outside the ghost)
+ *
+ * \return The external part of the dataset
+ *
+ */
+ virtual T getExternal();
+
+ //! divide the dataset from internal part and border
+ virtual dataDiv<T> divide();
+
+ //! Get the number of hyper-cube the space id is divided into
+ virtual size_t getNHyperCube(size_t id);
+
+ //! Get the hyper-cube margins
+ virtual std::vector<T> & getHyperCube(size_t id, size_t id_c);
+
+ //! destructor
+ virtual ~Decomposition(){}
+};
+
+#endif
diff --git a/src/Decomposition/DistModel.hpp b/src/Decomposition/DistModel.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5980dc8ac8a4731ebfb689ea09c02e91e92e4fd8
--- /dev/null
+++ b/src/Decomposition/DistModel.hpp
@@ -0,0 +1,30 @@
+#ifndef DIST_MODEL_HPP
+#define DIST_MODEL_HPP
+
+#include "metis.h"
+
+/*! \brief This class do graph partitioning
+ *
+ * This class do graph partitioning, it use METIS internaly
+ *
+ */
+
+template<typename Graph, typename algorithm>
+class GraphPartitioning
+{
+ //! Structure that store the graph
+ Graph & grp;
+
+ /*! Constructor
+ *
+ * It load the graph to partition
+ *
+ * \param g Graph to store
+ *
+ */
+ GraphPartitioning(Graph & g)
+ :grp(g)
+ {}
+};
+
+#endif
diff --git a/src/Grid/grid_dist_id.hpp b/src/Grid/grid_dist_id.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..12e53afeccdc99f815d073eb15b4de6cca70bff6
--- /dev/null
+++ b/src/Grid/grid_dist_id.hpp
@@ -0,0 +1,209 @@
+#ifndef COM_UNIT_HPP
+#define COM_UNIT_HPP
+
+#include <vector>
+#include "Grid/map_grid.hpp"
+#include "VCluster.hpp"
+#include "Space/SpaceBox.hpp"
+#include "mathutil.hpp"
+#include "grid_dist_id_iterator.hpp"
+#include "grid_dist_id_iterator_margin.hpp"
+#include "grid_dist_key.hpp"
+
+#define SUB_UNIT_FACTOR 64
+
+
+/*! \brief This is a distributed grid
+ *
+ * Implementation of a distributed grid with id decomposition. A distributed grid is a grid distributed
+ * across processors. The decomposition is performed on the id of the elements
+ *
+ * [Examples]
+ *
+ * on 1D where the id is from 1 to N
+ * processor k take M contiguous elements
+ *
+ * on 3D where (for example)
+ * processor k take M id-connected elements
+ *
+ * \param dim Dimensionality of the grid
+ * \param T type of grid
+ * \param Decomposition Class that decompose the grid for example CartDecomposition
+ * \param Mem Is the allocator
+ * \param device type of base structure is going to store the data
+ *
+ */
+
+template<unsigned int dim, typename T, typename Decomposition,typename Memory=HeapMemory , typename device_grid=grid_cpu<dim,T> >
+class grid_dist_id
+{
+ // Ghost expansion
+ Box<dim,size_t> ghost;
+
+ //! Local grids
+ Vcluster_object_array<device_grid> loc_grid;
+
+ //! Space Decomposition
+ Decomposition dec;
+
+ //! Size of the grid on each dimension
+ size_t g_sz[dim];
+
+ //! Communicator class
+
+ Vcluster & v_cl;
+
+ /*! \brief Get the grid size
+ *
+ * Get the grid size, given a domain, the resolution on it and another spaceBox
+ * it give the size on all directions of the local grid
+ *
+ * \param sp SpaceBox enclosing the local grid
+ * \param domain Space box enclosing the physical domain or part of it
+ * \param v_size grid size on this physical domain
+ *
+ * \return An std::vector representing the local grid on each dimension
+ *
+ */
+ std::vector<size_t> getGridSize(SpaceBox<dim,typename Decomposition::domain_type> & sp, Box<dim,typename Decomposition::domain_type> & domain, size_t (& v_size)[dim])
+ {
+ std::vector<size_t> tmp;
+ for (size_t d = 0 ; d < dim ; d++)
+ {
+ //! Get the grid size compared to the domain space and its resolution
+ typename Decomposition::domain_type dim_sz = (sp.getHigh(d) - sp.getLow(d)) / ((domain.getHigh(d) - domain.getLow(d)) / v_size[d]) + 0.5;
+
+ // push the size of the local grid
+ tmp.push_back(dim_sz);
+ }
+ return tmp;
+ }
+
+ /*! \brief Get the grid size
+ *
+ * Get the grid size, given a spaceBox
+ * it give the size on all directions of the local grid
+ *
+ * \param sp SpaceBox enclosing the local grid
+ * \param sz array to fill with the local grid size on each dimension
+ *
+ */
+ void getGridSize(SpaceBox<dim,size_t> & sp, size_t (& v_size)[dim])
+ {
+ for (size_t d = 0 ; d < dim ; d++)
+ {
+ // push the size of the local grid
+ v_size[d] = sp.getHigh(d) - sp.getLow(d);
+ }
+ }
+
+public:
+
+ //! constructor
+ grid_dist_id(Vcluster v_cl, Decomposition & dec, size_t (& g_sz)[dim], Box<dim,size_t> & ghost)
+ :ghost(ghost),loc_grid(NULL),v_cl(v_cl),dec(dec)
+ {
+ // fill the global size of the grid
+ for (int i = 0 ; i < dim ; i++) {this->g_sz[i] = g_sz[i];}
+
+ // Get the number of processor and calculate the number of sub-domain
+ // for decomposition
+ size_t n_proc = v_cl.getProcessingUnits();
+ size_t n_sub = n_proc * SUB_UNIT_FACTOR;
+
+ // Calculate the maximum number (before merging) of sub-domain on
+ // each dimension
+ size_t div[dim];
+ for (int i = 0 ; i < dim ; i++)
+ {div[i] = round_big_2(pow(n_sub,1.0/dim));}
+
+ // Create the sub-domains
+ dec.setParameters(div);
+ }
+
+ //! constructor
+ grid_dist_id(size_t (& g_sz)[dim])
+ :v_cl(*global_v_cluster),dec(Decomposition(v_cl))
+ {
+ // fill the global size of the grid
+ for (int i = 0 ; i < dim ; i++) {this->g_sz[i] = g_sz[i];}
+
+ // first compute a decomposition
+ Create();
+ }
+
+ /*! \brief Get the object that store the decomposition information
+ *
+ * \return the decomposition object
+ *
+ */
+
+ Decomposition & getDecomposition()
+ {
+ return dec;
+ }
+
+ /*! \brief Create the grid on memory
+ *
+ */
+
+ void Create()
+ {
+ // ! Create an hyper-cube approximation.
+ // ! In order to work on grid_dist the decomposition
+ // ! has to be a set of hyper-cube
+
+ dec.hyperCube();
+
+ // Get the number of local grid needed
+
+ size_t n_grid = dec.getNLocalHyperCube();
+
+ // create local grids for each hyper-cube
+
+ loc_grid = v_cl.allocate<device_grid>(n_grid);
+
+ // Size of the grid on each dimension
+ size_t l_res[dim];
+
+ // Allocate the grids
+
+ for (size_t i = 0 ; i < n_grid ; i++)
+ {
+ // Get the local hyper-cube
+
+ SpaceBox<dim,size_t> sp = dec.getLocalHyperCube(i);
+
+ // Calculate the local grid size
+
+ getGridSize(sp,l_res);
+
+ // Set the dimensions of the local grid
+
+ loc_grid.get(i).template resize<Memory>(l_res);
+ }
+ }
+
+ /*! \brief It return an iterator of the bulk part of the grid with a specified margin
+ *
+ * For margin we mean that every point is at least m points far from the border
+ *
+ * \param m margin
+ *
+ * \return An iterator to a grid with specified margins
+ *
+ */
+ grid_dist_iterator_margin<dim,device_grid> getBulkIterator(size_t margin)
+ {
+ grid_dist_iterator_margin<dim,device_grid> it(loc_grid,margin);
+
+ return it;
+ }
+
+ //! Destructor
+ ~grid_dist_id()
+ {
+ }
+};
+
+#endif
diff --git a/src/Grid/grid_dist_id_iterator.hpp b/src/Grid/grid_dist_id_iterator.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..3121c34cfa863cfa8f73d8fb0723bb3f1536cb50
--- /dev/null
+++ b/src/Grid/grid_dist_id_iterator.hpp
@@ -0,0 +1,111 @@
+/*
+ * grid_dist_id_iterator.hpp
+ *
+ * Created on: Feb 4, 2015
+ * Author: Pietro Incardona
+ */
+
+#ifndef GRID_DIST_ID_ITERATOR_HPP_
+#define GRID_DIST_ID_ITERATOR_HPP_
+
+#include "grid_dist_key.hpp"
+#include "Grid/grid.hpp"
+
+/*! \brief Distributed grid iterator
+ *
+ * Iterator across the local element of the distributed grid
+ *
+ */
+
+template<unsigned int dim, typename l_grid>
+class grid_dist_iterator
+{
+ //! grid list counter
+
+ size_t g_c;
+
+ //! List of the grids we are going to iterate
+ std::vector<l_grid> & gList;
+
+ //! Actual iterator
+ grid_key_dx_iterator<dim> a_it;
+
+ public:
+
+ /*! \brief Constructor of the distributed grid
+ *
+ * \param gk std::vector of the local grid
+ *
+ */
+ grid_dist_iterator(std::vector<l_grid> & gk)
+ :g_c(0),gList(gk)
+ {
+ // Initialize with the current iterator
+ // with the first grid
+
+ a_it = gList[0].getIterator();
+ }
+
+ // Destructor
+ ~grid_dist_iterator()
+ {}
+
+ /*! \brief Get the next element
+ *
+ * \return the next grid_key
+ *
+ */
+
+ grid_key_dx_iterator<dim> operator++()
+ {
+ a_it++;
+
+ // check if a_it is at the end
+
+ if (a_it.isEnd() == false)
+ return *this;
+ else
+ {
+ // switch to the new grid
+
+ g_c++;
+
+ // get the next grid iterator
+
+ a_it = a_it = gList[g_c].getIterator();
+
+ // increment to a valid point
+
+ a_it++;
+ }
+
+ return *this;
+ }
+
+ /*! \brief Check if there is the next element
+ *
+ * \return true if there is the next, false otherwise
+ *
+ */
+
+ bool isEnd()
+ {
+ // If there are no other grid stop
+
+ if (g_c >= gList.size())
+ return true;
+ }
+
+ /*! \brief Get the actual key
+ *
+ * \return the actual key
+ *
+ */
+ grid_dist_key_dx<dim> get()
+ {
+ return a_it;
+ }
+};
+
+
+#endif /* GRID_DIST_ID_ITERATOR_HPP_ */
diff --git a/src/Grid/grid_dist_id_iterator_margin.hpp b/src/Grid/grid_dist_id_iterator_margin.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..52192dd954d0a29cc0bbbba9a812e62e8aef9c91
--- /dev/null
+++ b/src/Grid/grid_dist_id_iterator_margin.hpp
@@ -0,0 +1,116 @@
+/*
+ * grid_dist_id_iterator_sub.hpp
+ *
+ * Created on: Feb 4, 2015
+ * Author: Pietro Incardona
+ */
+
+#ifndef GRID_DIST_ID_ITERATOR_SUB_HPP_
+#define GRID_DIST_ID_ITERATOR_SUB_HPP_
+
+#include "VCluster.hpp"
+
+/*! \brief Distributed grid iterator
+ *
+ * Iterator across the local element of the distributed grid
+ *
+ */
+
+template<unsigned int dim, typename device_grid>
+class grid_dist_iterator_margin
+{
+ //! grid list counter
+ size_t g_c;
+
+ //! List of the grids we are going to iterate
+ Vcluster_object_array<device_grid> & gList;
+
+ //! Actual iterator
+ grid_key_dx_iterator_sub<dim> * a_it;
+
+ //! margin of the grid iterator
+ size_t m;
+
+ public:
+
+ /*! \brief Constructor of the distributed grid
+ *
+ * \param gk std::vector of the local grid
+ *
+ */
+ grid_dist_iterator_margin(Vcluster_object_array<device_grid> & gk, size_t m)
+ :g_c(0),gList(gk),m(m)
+ {
+ // Initialize the current iterator
+ // with the first grid
+
+ a_it = new grid_key_dx_iterator_sub<dim>(gList[0].getSubIterator(m));
+ }
+
+ // Destructor
+ ~grid_dist_iterator_margin()
+ {
+ // delete the grid iterator
+
+ delete a_it;
+ }
+
+ /*! \brief Get the next element
+ *
+ * \return the next grid_key
+ *
+ */
+
+ grid_key_dx_iterator<dim> operator++()
+ {
+ a_it++;
+
+ // check if a_it is at the end
+
+ if (a_it.isEnd() == false)
+ return *this;
+ else
+ {
+ // switch to the new grid
+
+ g_c++;
+
+ // get the next grid iterator
+
+ a_it = a_it = gList[g_c].getIterator();
+
+ // increment to a valid point
+
+ a_it++;
+ }
+
+ return *this;
+ }
+
+ /*! \brief Check if there is the next element
+ *
+ * \return true if there is the next, false otherwise
+ *
+ */
+
+ bool isEnd()
+ {
+ // If there are no other grid stop
+
+ if (g_c >= gList.size())
+ return true;
+ }
+
+ /*! \brief Get the actual key
+ *
+ * \return the actual key
+ *
+ */
+ grid_dist_key_dx<dim> get()
+ {
+ return a_it;
+ }
+};
+
+
+#endif /* GRID_DIST_ID_ITERATOR_SUB_HPP_ */
diff --git a/src/Grid/grid_dist_id_unit_test.hpp b/src/Grid/grid_dist_id_unit_test.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e4ca239343054c9a97b67c9a9313cc24cd4b91b2
--- /dev/null
+++ b/src/Grid/grid_dist_id_unit_test.hpp
@@ -0,0 +1,109 @@
+#ifndef GRID_DIST_UNIT_TEST_HPP
+#define GRID_DIST_UNIT_TEST_HPP
+
+#include "grid_dist_id.hpp"
+#include "data_type/scalar.hpp"
+
+BOOST_AUTO_TEST_SUITE( grid_dist_id_test )
+
+template<typename iterator> void jacobi_iteration(iterator g_it, grid_dist_id<2, scalar<float>, CartDecomposition<2,size_t>> & g_dist)
+{
+ // scalar
+ typedef scalar<size_t> S;
+
+ // iterator
+
+ while(g_it.isNext())
+ {
+ // Jacobi update
+
+ auto pos = g_it.get();
+
+ g_dist.template get<S::ele>(pos) = (g_dist.template get<S::ele>(pos.move(0,1)) +
+ g_dist.template get<S::ele>(pos.move(0,-1)) +
+ g_dist.template get<S::ele>(pos.move(1,1)) +
+ g_dist.template get<S::ele>(pos.move(1,-1)) / 4.0);
+
+ ++g_it;
+ }
+}
+
+BOOST_AUTO_TEST_CASE( grid_dist_id_iterator_test_use)
+{
+ // grid size
+ size_t sz[2] = {1024,1024};
+
+ // Distributed grid with id decomposition
+
+ grid_dist_id<2, scalar<float>, CartDecomposition<2,size_t>> g_dist(sz);
+
+ // Create the grid on memory
+
+ g_dist.Create();
+
+ // get the Bulk iterator
+
+ auto bulk = g_dist.getBulkIterator(2);
+
+/* auto g_it = g_dist.getIteratorBulk();
+
+ auto g_it_halo = g_dist.getHalo();
+
+ // Let try to solve the poisson equation d2(u) = f with f = 1 and computation
+ // comunication overlap (100 Jacobi iteration)
+
+ for (int i = 0 ; i < 100 ; i++)
+ {
+ g_dist.ghost_get();
+
+ // Compute the bulk
+
+ jacobi_iteration(g_it);
+
+ g_dist.ghost_sync();
+
+ // Compute the halo
+
+ jacobi_iteration(g_it_halo);
+ }*/
+}
+
+BOOST_AUTO_TEST_CASE( grid_dist_id_poisson_test_use)
+{
+ // grid size
+ size_t sz[2] = {1024,1024};
+
+ // Distributed grid with id decomposition
+
+ grid_dist_id<2, scalar<float>, CartDecomposition<2,size_t>> g_dist(sz);
+
+ // Create the grid on memory
+
+ g_dist.Create();
+
+/* auto g_it = g_dist.getIteratorBulk();
+
+ auto g_it_halo = g_dist.getHalo();
+
+ // Let try to solve the poisson equation d2(u) = f with f = 1 and computation
+ // comunication overlap (100 Jacobi iteration)
+
+ for (int i = 0 ; i < 100 ; i++)
+ {
+ g_dist.ghost_get();
+
+ // Compute the bulk
+
+ jacobi_iteration(g_it);
+
+ g_dist.ghost_sync();
+
+ // Compute the halo
+
+ jacobi_iteration(g_it_halo);
+ }*/
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+#endif
diff --git a/src/Grid/grid_dist_key.hpp b/src/Grid/grid_dist_key.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..3a45a1e198b63220d80a302fc8530814bb80db24
--- /dev/null
+++ b/src/Grid/grid_dist_key.hpp
@@ -0,0 +1,22 @@
+#ifndef GRID_DIST_KEY_DX_HPP
+#define GRID_DIST_KEY_DX_HPP
+
+/*! \brief Grid key for a distributed grid
+ *
+ * Grid key for a distributed grid
+ *
+ */
+
+template<unsigned int dim>
+class grid_dist_key_dx
+{
+ //! grid list counter
+
+ size_t g_c;
+
+ //! Local grid iterator
+
+ grid_key_dx_iterator<dim> a_it;
+};
+
+#endif