diff --git a/Makefile.am b/Makefile.am
new file mode 100755
index 0000000000000000000000000000000000000000..7c618d3e326a45c8f6719691f58bde9a5a1c4c26
--- /dev/null
+++ b/Makefile.am
@@ -0,0 +1,3 @@
+SUBDIRS = src
+
+bin_PROGRAMS = 
\ No newline at end of file
diff --git a/configure.ac b/configure.ac
new file mode 100755
index 0000000000000000000000000000000000000000..fb719e480b548e5be486efc3caa43c783b9a7362
--- /dev/null
+++ b/configure.ac
@@ -0,0 +1,204 @@
+#                                               -*- Autoconf -*-
+# Process this file with autoconf to produce a configure script.
+
+
+AC_PREREQ(2.59)
+AC_INIT(FULL-PACKAGE-NAME, VERSION, BUG-REPORT-ADDRESS)
+AC_CANONICAL_SYSTEM
+AC_CONFIG_SRCDIR([src/main.cpp])
+AM_INIT_AUTOMAKE
+AC_CONFIG_HEADER([src/config.h])
+m4_ifdef([MYSQL_FOUND],,[m4_include([m4/ax_lib_mysql.m4])])
+m4_ifdef([AX_CHECK_COMPILER_FLAGS],,[m4_include([m4/ax_check_compiler_flags.m4])])
+m4_ifdef([ACX_PTHREAD],,[m4_include([m4/acx_pthread.m4])])
+m4_ifdef([AX_CHECK_CL],,[m4_include([m4/ax_opencl.m4])])
+m4_ifdef([AX_BOOST_BASE],,[m4_include([m4/ax_boost_base.m4])])
+m4_ifdef([AX_BOOST_PROGRAM_OPTIONS],,[m4_include([m4/ax_boost_program_options.m4])])
+m4_ifdef([AX_BOOST_THREAD],,[m4_include([m4/ax_boost_thread.m4])])
+m4_ifdef([ACX_MPI],,[m4_include([m4/acx_mpi.m4])])
+m4_ifdef([AX_OPENMP],,[m4_include([m4/ax_openmp.m4])])
+m4_ifdef([AX_GCC_X86_CPUID],,[m4_include([m4/ax_gcc_x86_cpuid.m4])])
+m4_ifdef([AX_GCC_ARCHFLAG],,[m4_include([m4/ax_gcc_archflag.m4])])
+m4_ifdef([AX_CUDA],,[m4_include([m4/ax_cuda.m4])])
+
+CXXFLAGS+=" --std=c++11 "
+NVCCFLAGS=" "
+INCLUDES_PATH=" "
+
+# Checks for programs.
+AC_PROG_CXX
+
+# Checks g++ flags
+
+AC_CANONICAL_HOST
+
+# Check target architetture
+
+AX_GCC_ARCHFLAG([], [CXXFLAGS="$CXXFLAGS $ax_cv_gcc_archflag"], [])
+
+###### Check for debug compilation
+
+AC_MSG_CHECKING(whether to build with debug information)
+debuger=no
+AC_ARG_ENABLE(debug,
+	AC_HELP_STRING(
+	    [--enable-debug],
+	    [enable debug data generation (def=no)]
+	),
+  	debuger="$enableval"
+)
+
+
+
+AC_MSG_RESULT($debuger)
+if test x"$debuger" = x"yes"; then
+	AC_DEFINE([DEBUG_MODE],[],[Debug])
+	AC_DEFINE([DEBUG],[],[Debug])
+	CXXFLAGS="$CXXFLAGS -g3  -Wall -O0 "
+	NVCCFLAGS+="$NVCCFLAGS -g -O0 "
+else
+	CXXFLAGS="$CXXFLAGS -Wall -O3 -g3 "
+	NVCCFLAGS+="$NVCCFLAGS -O3 "
+fi
+
+## Check for memcheck
+
+AC_MSG_CHECKING(whether to build with memcheck capabilities)
+debuger=no
+AC_ARG_ENABLE(memcheck,
+        AC_HELP_STRING(
+            [--enable-memcheck],
+            [enable memory check (def=no)]
+        ),
+        memcheck="$enableval"
+)
+
+AC_MSG_RESULT($memcheck)
+
+if test x"$memcheck" = x"yes"; then
+        AC_DEFINE([MEMLEAK_CHECK],[],[Memory check, corruption and leak])
+fi
+
+#########
+
+####### include OpenFPM_devices include path
+
+INCLUDES_PATH+="-I. -I../../metis_install/include -I../../OpenFPM_IO/src -I../../OpenFPM_data/src -I../../OpenFPM_devices/src -I../../OpenFPM_vcluster/src/"
+
+####### Checking for GPU support
+
+AX_CUDA
+
+AC_MSG_CHECKING(whether to build with GPU support)
+gpu_support=no
+AC_ARG_ENABLE(gpu,
+	AC_HELP_STRING(
+	    [--enable-gpu],
+	    [enable gpu support]
+	),
+  	gpu_support="$enableval"
+)
+
+
+
+AC_MSG_RESULT($gpu_support)
+if test x"$gpu_support" = x"yes"; then
+	AC_DEFINE([GPU],[],[GPU support])
+fi
+
+
+# Set this conditional if cuda is wanted
+
+AM_CONDITIONAL(BUILDCUDA, test ! x$NVCC = x"no")
+
+###########################
+
+no_avx=no
+no_sse42=no
+no_sse41=no
+no_sse3=no
+no_sse2=no
+no_sse=no
+no_mmx=no
+AX_CHECK_COMPILER_FLAGS([-msse4.2],[CXXFLAGS="$CXXFLAGS -mavx"],[no_avx=yes])
+AX_CHECK_COMPILER_FLAGS([-msse4.2],[CXXFLAGS="$CXXFLAGS -msse4.2"],[no_sse42=yes])
+AX_CHECK_COMPILER_FLAGS([-msse4.1],[CXXFLAGS="$CXXFLAGS -msse4.1"],[no_sse41=yes])
+AX_CHECK_COMPILER_FLAGS([-msse3],[CXXFLAGS="$CXXFLAGS -msse3"],[no_sse3=yes])
+AX_CHECK_COMPILER_FLAGS([-msse2],[CXXFLAGS="$CXXFLAGS -msse2"],[no_sse2=yes])
+AX_CHECK_COMPILER_FLAGS([-msse],[CXXFLAGS="$CXXFLAGS -msse"],[no_sse=yes])
+AX_CHECK_COMPILER_FLAGS([-mmmx],[CXXFLAGS="$CXXFLAGS -mmmx"],[no_mmx=yes])
+AX_CHECK_COMPILER_FLAGS([-Wno-unused-but-set-variable],[CXXFLAGS="$CXXFLAGS -Wno-unused-but-set-variable"],[])
+
+AC_SUBST(NVCCFLAGS)
+AC_SUBST(INCLUDES_PATH)
+
+# Checks for typedefs, structures, and compiler characteristics.
+
+# Checks for library functions.
+
+AC_CONFIG_FILES([Makefile src/Makefile])
+AC_OUTPUT
+echo ""
+echo   "***********************************"
+echo   "*                                 *"
+arch_str="${ax_cv_gcc_archflag#-march=#-mtune=}"
+arch_str="${arch_str#-mtune=}"
+n_arch_str=${#arch_str}
+for (( X=0; X<23-n_arch_str; X++ ))
+do
+   arch_str="$arch_str "
+done
+echo "*    arch: $arch_str*"
+if [ test x"$no_sse42" = x"no" ]; then
+    echo "*    sse4.2: yes                  *"
+else
+    echo "*    sse4.2: no                   *"
+fi
+if [ test x"$no_sse41" = x"no" ]; then
+    echo "*    sse4.1: yes                  *"
+else
+    echo "*    sse4.1: no                   *"
+fi
+if [ test x"$no_sse3" = x"no" ]; then
+    echo "*    sse3: yes                    *"
+else
+    echo "*    sse3: no                     *"
+fi
+if [ test x"$no_sse2" = x"no" ]; then
+    echo "*    sse2: yes                    *"
+else
+    echo "*    sse2: no                     *"
+fi
+if [ test x"$no_sse" = x"no" ]; then
+    echo "*    sse: yes                     *"
+else
+    echo "*    sse: no                      *"
+fi
+if [ test x"$no_mmx" = x"no" ]; then
+  echo "*    mmx: yes                     *"
+else
+  echo "*    mmx: no                      *"
+fi
+if [ test x"$profiler" = x"yes" ]; then
+  echo "*    profiler: yes                *"
+else
+  echo "*    profiler: no                 *"
+fi
+if [ test x"$debuger" = x"yes" ]; then
+  echo "*    debug: yes                   *"
+else
+  echo "*    debug: no                    *"
+fi
+if [ test x"$no_64" = x"no" ]; then
+  echo "*    64 bit: yes                  *"
+else
+  echo "*    64 bit: no                   *"
+fi
+if [ test x"$gpu_support" = x"no" ]; then
+  echo "*    gpu: no                      *"
+else
+  echo "*    gpu: yes                     *"
+fi
+echo   "*                                 *"
+echo   "***********************************"
+
diff --git a/src/dec_optimizer.hpp b/src/dec_optimizer.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..3eeba6451ce7868e7c9fa9daad0b945bfefa9d85
--- /dev/null
+++ b/src/dec_optimizer.hpp
@@ -0,0 +1,475 @@
+#ifndef DEC_OPTIMIZER_HPP
+#define DEC_OPTIMIZER_HPP
+
+/*! \brief this class represent a wavefront of dimension dim
+ *
+ * \dim Dimensionality of the wavefront (dimensionality of the space
+ *                                       where it live so the wavefront
+ *                                       is dim-1)
+ *
+ */
+
+template <unsigned int dim>
+class wavefront
+{
+public:
+
+	typedef boost::fusion::vector<size_t[dim],size_t[dim]> type;
+
+	typedef typename memory_traits_inte<type>::type memory_int;
+	typedef typename memory_traits_lin<type>::type memory_lin;
+
+	type data;
+
+	static const int start = 0;
+	static const int stop = 1;
+	static const int max_prop = 2;
+
+	/* \brief Get the key to the point 1
+	 *
+	 * \return the key to the point 1
+	 *
+	 */
+
+	grid_key_dx<dim> getKP1()
+	{
+		// grid key to return
+		grid_key_dx<dim> ret(boost::fusion::at_c<start>(data));
+
+		return ret;
+	}
+
+	/* \brief Get the key to point 2
+	 *
+	 * \return the key to the point 2
+	 *
+	 */
+
+	grid_key_dx<dim> getKP2()
+	{
+		// grid key to return
+		grid_key_dx<dim> ret(boost::fusion::at_c<stop>(data));
+
+		return ret;
+	}
+
+	/* \brief produce a box from an encapsulated object
+	 *
+	 * \param encap encapsulated object
+	 *
+	 */
+
+	template<typename encap> static Box<dim,size_t> getBox(encap & enc)
+	{
+		Box<dim,size_t> bx;
+
+		// Create the object from the encapsulation
+
+		for (int i = 0 ; i < dim ; i++)
+		{
+			bx.setHigh(i,enc.template get<wavefront::start>()[i]);
+			bx.setLow(i,enc.template get<wavefront::stop>()[i]);
+		}
+
+		return bx;
+	}
+
+	/* \brief produce a box from an encapsulated object
+	 *
+	 * \param encap encapsulated object
+	 *
+	 */
+
+	template<typename encap> static void getBox(encap & enc, Box<dim,size_t> & bx)
+	{
+		// Create the object from the encapsulation
+
+		for (int i = 0 ; i < dim ; i++)
+		{
+			bx.setHigh(i,enc.template get<wavefront::start>()[i]);
+			bx.setLow(i,enc.template get<wavefront::stop>()[i]);
+		}
+	}
+};
+
+/*! \brief This class take a graph representing the space decomposition and produce a
+ *         simplified version
+ *
+ * Given a Graph_CSR and a seed point produce an alternative decomposition to reduce
+ * the ghost over-stress
+ *
+ */
+
+template <unsigned int dim, typename Graph>
+class dec_optimizer
+{
+	// create a grid header for helping
+
+	grid<dim,void> gh;
+
+private:
+
+	/* \brief Fill the wavefront position
+	 *
+	 * \tparam prp property to set
+	 *
+	 * \param graph we are processing
+	 * \param Box to fill
+	 * \param id value to fill with
+	 *
+	 */
+
+	template<unsigned int prp> void write_wavefront(Graph & graph,openfpm::vector<wavefront<dim>> & v_w)
+	{
+		// fill the wall domain with 0
+
+		fill_domain<prp>(graph,gh.getBox(),0);
+
+		// fill the wavefront
+
+		for (int i = 0 ; i < v_w.size() ; i++)
+		{
+			Box<dim,size_t> box = wavefront<dim>::getBox(v_w.get(i));
+
+			fill_domain<prp>(graph,box,1);
+		}
+	}
+
+	/* \brief Fill the domain
+	 *
+	 * \tparam p_sub property to set with the sub-domain id
+	 *
+	 * \param graph we are processing
+	 * \param Box to fill
+	 * \param ids value to fill with
+	 *
+	 */
+
+	template<unsigned int p_sub> void fill_domain(Graph & graph, Box<dim,size_t> & box, size_t ids)
+	{
+		// Create a subgrid iterator
+		grid_key_dx_iterator_sub<dim> g_sub(gh,box.getKP1(),box.getKP2());
+
+		// iterate through all grid points
+
+		while (g_sub.isNext())
+		{
+			// get the actual key
+			const grid_key_dx<dim> & gk = g_sub.get();
+
+			// get the vertex and set the sub id
+
+			graph.vertex(gh.LinId(gk)).template get<p_sub>() = ids;
+		}
+	}
+
+	/* \brief Add the boundary domain to the queue
+	 *
+	 * \tparam i-property where is stored the decomposition
+	 *
+	 * \param domains vector with domains to process
+	 * \param graph we are processing
+	 * \param w_comb hyper-cube combinations
+	 *
+	 */
+	template<unsigned int p_sub> void add_to_queue(openfpm::vector<size_t> & domains, openfpm::vector<wavefront<dim>> & v_w, Graph & graph,  std::vector<comb<dim>> & w_comb)
+	{
+		// create a new queue
+		openfpm::vector<size_t> domains_new;
+
+		// push in the new queue, the old domains of the queue that are not assigned element
+
+		for (size_t j = 0 ; j < domains.size() ; j++)
+		{
+			if (graph.vertex(domains.get(j)).template get<p_sub>() < 0)
+			{
+				// not assigned push it
+
+				domains_new.add(domains.get(j));
+			}
+		}
+
+		// take the wavefront expand on direction d of one
+
+		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];
+			}
+		}
+
+		// for each expanded wavefront create a sub-grid iterator and add the subbomain
+
+		for (int d = 0 ; d < v_w.size() ; d++)
+		{
+			// Create a sub-grid iterator
+			grid_key_dx_iterator_sub<dim> g_sub(gh,v_w.template get<wavefront<dim>::start>(d),v_w.template get<wavefront<dim>::stop>(d));
+
+			// iterate through all grid points
+
+			while (g_sub.isNext())
+			{
+				// get the actual key
+				const grid_key_dx<dim> & gk = g_sub.get();
+
+				// get the vertex and check if is not assigned add to the queue
+
+				if (graph.vertex(gh.LinId(gk)).template get<p_sub>() < 0)
+				{
+					// add
+
+					domains_new.add(gh.LinId(gk));
+				}
+			}
+		}
+
+		// copy the new queue to the old one (it not copied, C++11 move semantic)
+		domains = domains_new;
+	}
+
+	/* \brief Find the biggest hyper-cube
+	 *
+	 * starting from one initial sub-domain find the biggest hyper-cube
+	 *
+	 * \tparam j id of the property storing the sub-decomposition
+	 * \tparam i id of the property containing the decomposition
+	 *
+	 * \param start_p initial domain
+	 * \param graph representing the grid
+	 *
+	 */
+	template <unsigned int p_sub, unsigned int p_id> void expand_from_point(size_t start_p, Graph & graph, Box<dim,size_t> & box, openfpm::vector<wavefront<dim>> & v_w , std::vector<comb<dim>> & w_comb)
+	{
+		// We assume that Graph is the rapresentation of a cartesian graph
+		// this mean that the direction d is at the child d
+
+		// Get the number of wavefronts
+		size_t n_wf = w_comb.size();
+
+		// Create an Hyper-cube
+		HyperCube<dim> hyp;
+
+		// direction of expansion
+
+		size_t domain_id = graph.vertex(gh.LinId(start_p)).template get<p_id>();
+		bool can_expand = true;
+
+		// while is possible to expand
+
+		while (can_expand)
+		{
+			// reset can expand
+			can_expand = false;
+
+			// for each direction of expansion expand the wavefront
+
+			for (int d = 0 ; d < n_wf ; d++)
+			{
+				// number of processed sub-domain
+				size_t n_proc_sub = 0;
+
+				// flag to indicate if the wavefront can expand
+				bool w_can_expand = true;
+
+				// Create an iterator of the wavefront
+				grid_key_dx_iterator_sub<dim> it(gh,v_w.template get<wavefront<dim>::start>(d),v_w.template get<wavefront<dim>::stop>(d));
+
+				// for each subdomain
+				while (it.isNext())
+				{
+					// get the wavefront sub-domain id
+					size_t sub_w = gh.LinId(it.get());
+
+					// we get the processor id of the neighborhood sub-domain on direction d
+					size_t exp_p = graph.vertex(graph.getChild(sub_w,d)).template get<p_id>();
+
+					// we check if it is the same processor id
+					w_can_expand &= exp_p != domain_id;
+
+					// next domain
+					++it;
+
+					// increase the number of processed sub-domain
+					n_proc_sub++;
+				}
+
+				// if we did not processed sub-domain, we cannot expand
+				w_can_expand &= (n_proc_sub != 0);
+
+				// if you can expand one wavefront we did not reach the end
+				can_expand |= w_can_expand;
+
+				// if we can expand the wavefront expand it
+				if (w_can_expand == true)
+				{
+					// expand the wavefront
+					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];
+
+						std::cout << "stop: " << v_w.template get<wavefront<dim>::stop>(d)[j] << "\n";
+						std::cout << "start: " << v_w.template get<wavefront<dim>::start>(d)[j] << "\n";
+					}
+
+					// 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)[j] = v_w.template get<wavefront<dim>::stop>(d)[j] + w_comb[id].c[j];}
+							else
+							{v_w.template get<wavefront<dim>::start>(id)[j] = v_w.template get<wavefront<dim>::start>(d)[j] + w_comb[id].c[j];}
+						}
+					}
+				}
+			}
+
+			// Debug output the wavefront graph for debug
+
+			// duplicate the graph
+
+			Graph g_debug = graph.duplicate();
+			write_wavefront<nm_part_v::id>(g_debug,v_w);
+
+			VTKWriter<Graph> vtk(g_debug);
+			vtk.template write<1>("vtk_debug.vtk");
+
+			////////////////////////////////////
+		}
+
+		// get back the hyper-cube produced
+
+		for (int i = 0 ; i < dim ; i++)
+		{
+			// get the index of the wavefront direction
+			int w = 0;
+
+			for (int j = 0 ; j < dim ; j++)
+			{
+				if (w_comb[i].c[j] == 1)
+				{
+					w = j;
+					break;
+				}
+			}
+
+			// set the box
+			box.setHigh(i,v_w.template get<wavefront<dim>::stop>(i)[w]);
+			box.setLow(i,v_w.template get<wavefront<dim>::start>(i)[w]);
+		}
+	}
+
+public:
+
+	/*! \brief Constructor
+	 *
+	 * \param g Graph to simplify
+	 * \param sz size of the grid on each dimension
+	 *
+	 */
+
+	dec_optimizer(Graph & g, std::vector<size_t> sz)
+	:gh(sz)
+	{
+		// The graph g is suppose to represent a cartesian grid
+		// No check is performed on g
+	}
+
+	/*! \brief optimize the graph
+	 *
+	 * Starting from a domain (hyper-cubic), it create wavefront at the boundary and expand
+	 * the boundary until the wavefronts cannot expand any more.
+	 * To the domains inside the hyper-cube one sub-id is assigned. This procedure continue until
+	 * all the domain of one id has a sub-id
+	 *
+	 * \tparam j property containing the decomposition
+	 * \tparam i property to fill with the sub-decomposition
+	 *
+	 * \param Seed point
+	 * \param graph we are processing
+	 *
+	 */
+
+	template <unsigned int p_sub, unsigned int p_id> void optimize(grid_key_dx<dim> & start_p, Graph & graph)
+	{
+		// sub-domain id
+		size_t sub_id =  0;
+
+		// queue
+		openfpm::vector<size_t> v_q;
+
+		// Create an hyper-cube
+		HyperCube<dim> hyp;
+
+		// Get the wavefront combinations
+		std::vector<comb<dim>> w_comb = hyp.getCombinations_R(dim-1);
+
+		// wavefronts
+		openfpm::vector<wavefront<dim>> v_w(w_comb.size());
+
+		// Initialize the wavefronts from the domain start_p
+
+		for (int i = 0 ; i < v_w.size() ; i++)
+		{
+			for (int j = 0 ; j < dim ; j++)
+			{
+				v_w.template get<wavefront<dim>::start>(i)[j] = start_p.get(start_p.get(j));
+				v_w.template get<wavefront<dim>::stop>(i)[j] = start_p.get(start_p.get(j));
+			}
+		}
+
+		// push the first domain
+
+		v_q.add(gh.LinId(start_p));
+
+		while (v_q.size() != 0)
+		{
+			// Box
+			Box<dim,size_t> box;
+
+			// Create the biggest box containing the domain
+			expand_from_point<p_sub,p_id>(gh.LinId(start_p),graph,box,v_w,w_comb);
+
+			// fill the domain
+			fill_domain<p_sub>(graph,box,sub_id);
+
+			// create the queue
+			add_to_queue<p_sub>(v_q,v_w,graph,w_comb);
+
+			// increment the sub_id
+			sub_id++;
+		}
+	}
+};
+
+#endif
diff --git a/src/dec_optimizer_unit_test.hpp b/src/dec_optimizer_unit_test.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..7ab5aa2099f8a8e6d54a4a25c07de512e464e387
--- /dev/null
+++ b/src/dec_optimizer_unit_test.hpp
@@ -0,0 +1,71 @@
+/*
+ * dec_optimize.hpp
+ *
+ *  Created on: Jan 16, 2015
+ *      Author: i-bird
+ */
+
+#ifndef DEC_OPTIMIZE_HPP_
+#define DEC_OPTIMIZE_HPP_
+
+
+#include "Graph/CartesianGraphFactory.hpp"
+#include "map_graph.hpp"
+#include "metis_util.hpp"
+#include "dec_optimizer.hpp"
+
+#undef GS_SIZE
+#define GS_SIZE 8
+
+BOOST_AUTO_TEST_SUITE( dec_optimizer_test )
+
+BOOST_AUTO_TEST_CASE( dec_optimizer_test_use)
+{
+	CartesianGraphFactory<3,Graph_CSR<nm_v,nm_e>> g_factory;
+	CartesianGraphFactory<3,Graph_CSR<nm_part_v,nm_part_e>> g_factory_part;
+
+	// Cartesian grid
+	std::vector<size_t> sz;
+	sz.push_back(4);
+	sz.push_back(4);
+	sz.push_back(1);
+
+	// Box
+	Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
+
+	// Graph to decompose
+
+	Graph_CSR<nm_v,nm_e> g = g_factory.construct<nm_e::communication,float,2,0,1,2>(sz,box);
+
+	// Processor graph
+
+	Graph_CSR<nm_part_v,nm_part_e> gp = g_factory_part.construct<NO_EDGE,float,2>(sz,box);
+
+	// Convert the graph to metis
+
+	Metis<Graph_CSR<nm_v,nm_e>> met(g,16);
+
+	// decompose
+
+	met.decompose<nm_part_v::id>(gp);
+	met.decompose<nm_v::id>();
+
+	// optimize
+
+	dec_optimizer<3,Graph_CSR<nm_v,nm_e>> d_o(g,sz);
+
+	grid_key_dx<3> keyZero(0,0,0);
+	d_o.optimize<nm_v::sub_id,nm_v::id>(keyZero,g);
+
+	// Write the VTK file
+
+	VTKWriter<Graph_CSR<nm_part_v,nm_part_e>> vtk(gp);
+	vtk.write("vtk_partition.vtk");
+	VTKWriter<Graph_CSR<nm_v,nm_e>> vtk2(g);
+	vtk2.write("vtk_partition2.vtk");
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+
+#endif /* DEC_OPTIMIZE_HPP_ */
diff --git a/src/gargabe.hpp b/src/gargabe.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4b610b42fb754ac55e3362607858d2c20d893f64
--- /dev/null
+++ b/src/gargabe.hpp
@@ -0,0 +1,105 @@
+/*
+ * gargabe.hpp
+ *
+ *  Created on: Jan 13, 2015
+ *      Author: i-bird
+ */
+
+#ifndef GARGABE_HPP_
+#define GARGABE_HPP_
+
+
+
+	template <unsigned int j, unsigned int i, typename Graph> void optimize(size_t start_p, Graph & graph)
+	{
+		// We assume that Graph is the rapresentation of a cartesian graph
+		// this mean that the direction d is at the child d
+
+		// Create an Hyper-cube
+
+		HyperCube<dim> hyp;
+
+		// Get the number of wavefronts
+
+		size_t n_wf = hyp.getNumberOfElements_R(0);
+
+		// Get the number of intersecting wavefront
+
+
+
+		// Get the number of sub-dimensional common wavefront
+		// basically are a list of all the subdomain common to two or more
+
+		// Create n_wf wavefront queue
+
+		openfpm::vector<wavefront> v_w;
+		v.reserve(n_wf);
+
+		// direction of expansion
+
+		size_t domain_id = 0;
+		int exp_dir = 0;
+		bool can_expand = true;
+
+		// while is possible to expand
+
+		while (can_expand)
+		{
+			// for each direction of expansion expand the wavefront
+
+			for (int d = 0 ; d < n_wf ; d++)
+			{
+				// get the wavefront at direction d
+
+				openfpm::vector<size_t> & wf_d = v_w.get<wavefront::domains>(d);
+
+				// flag to indicate if the wavefront can expand
+
+				bool w_can_expand = true;
+
+				// for each subdomain
+
+				for (size_t sub = 0 ; sub < wf_d.size() ; sub++)
+				{
+					// check if the adjacent domain in direction d exist
+					// and is of the same id
+
+					// get the starting subdomain
+					size_t sub_w = wf_d.get<0>(sub);
+
+					// we get the processor id of the neighborhood sub-domain on direction d
+					size_t exp_p = graph.getChild(sub_w,d).get<j>();
+
+					// we check if it is the same processor id
+					if (exp_p != domain_id)
+					{
+						w_can_expand = false;
+					}
+				}
+
+				// if we can expand the wavefront expand it
+				if (w_can_expand == true)
+				{
+					// for each subdomain
+					for (size_t sub = 0 ; sub < wf_d.size() ; sub++)
+					{
+						// update the position of the wavefront
+						wf_d.get<0>(sub) = wf_d.get<0>(sub) + gh.stride(d);
+					}
+
+					// here we add sub-domains to all the other queues
+					// get the face of the hyper-cube
+
+					SubHyperCube<dim,dim-1> sub_hyp = hyp.getSubHyperCube(d);
+
+					std::vector<comb<dim>> q_comb = sub_hyp.getCombinations_R(dim-2);
+				}
+			}
+		}
+
+		// For each point in the Hyper-cube check if we can move the wave front
+
+
+	}
+
+#endif /* GARGABE_HPP_ */
diff --git a/src/hypercube_unit_test.hpp b/src/hypercube_unit_test.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ca63a875de2563cd56c0262f04b06c28b3678368
--- /dev/null
+++ b/src/hypercube_unit_test.hpp
@@ -0,0 +1,384 @@
+#ifndef HYPERCUBE_UNIT_TEST_HPP
+#define HYPERCUBE_UNIT_TEST_HPP
+
+/*! \brief Check if the linearization is correct
+ *
+ * \param cmb vector of combination to checl
+ *
+ */
+
+template<unsigned int dim> void check_lin(std::vector<comb<dim>> cmb)
+{
+	// Check the linearization
+	for (int i = 0 ; i < cmb.size() ; i++)
+	{
+		BOOST_REQUIRE_EQUAL(HyperCube<dim>::LinId(cmb[i]),i);
+	}
+}
+
+/*! \brief Check if the combinations are dinstict
+ *
+ * \param combs combinations to check
+ *
+ */
+
+template<unsigned int dim> bool isDinstict(std::vector<comb<dim>> combs)
+{
+	// Check if the combination are dinstinct
+
+	for (int i = 0 ; i < combs.size() ; i++)
+	{
+		for (int j = i+1 ; j < combs.size() ; j++)
+		{
+			if (combs[i] == combs[j])
+				return false;
+		}
+
+	}
+
+	return true;
+}
+
+/*! \brief isSubdecomposition check if combs are elements that exist in c as sub-elements
+ *
+ * \param combs elements to check
+ * \param element that must contain all the combs
+ *
+ */
+
+template<unsigned int dim> bool isSubdecomposition(std::vector<comb<dim>> combs, comb<dim> c)
+{
+	for (int i = 0 ; i < combs.size() ; i++)
+	{
+		if (c.isSub(combs[i]) == false)
+			return false;
+	}
+
+	return true;
+}
+
+template<unsigned int dim> bool isValid(std::vector<comb<dim>> combs)
+{
+	// Check if the combinations are valid
+
+	for (int i = 0 ; i < combs.size() ; i++)
+	{
+		if (combs[i].isValid() == false)
+			return false;
+	}
+
+	return true;
+}
+
+BOOST_AUTO_TEST_SUITE( Hyper_cube )
+
+BOOST_AUTO_TEST_CASE( Hyper_cube_use)
+{
+	// Test Hyper-Cube functionality
+
+	size_t ele[8];
+
+	// Point
+//	ele[0] = HyperCube<0>::getNumberOfElements_R(0);
+
+	// Get combination for each dimensions
+//	std::vector<comb<0>> v_c0 = HyperCube<0>::getCombinations_R(0);
+
+	// Check
+//	BOOST_REQUIRE_EQUAL(ele[0],1);
+
+	// Line
+	// Number of vertex
+	ele[0] = HyperCube<1>::getNumberOfElements_R(0);
+	// Number of edge
+	ele[1] = HyperCube<1>::getNumberOfElements_R(1);
+
+	// Get combination for each dimensions
+	std::vector<comb<1>> v_c1_0 = HyperCube<1>::getCombinations_R(0);
+
+	// Check
+	BOOST_REQUIRE_EQUAL(ele[0],2);
+	BOOST_REQUIRE_EQUAL(ele[1],1);
+
+	// Fill the expected vector
+
+	comb<1> c[] = {1,-1};
+
+	// Check the linearization
+	check_lin(v_c1_0);
+
+	boost_check_array(&c[0],&v_c1_0[0],2);
+
+	// Square
+	// Number of vertex
+	ele[0] = HyperCube<2>::getNumberOfElements_R(0);
+	// Number of edge
+	ele[1] = HyperCube<2>::getNumberOfElements_R(1);
+	// Number of faces
+	ele[2] = HyperCube<2>::getNumberOfElements_R(2);
+
+	// Get combination for each dimensions
+	std::vector<comb<2>> v_c2_0 = HyperCube<2>::getCombinations_R(0);
+	std::vector<comb<2>> v_c2_1 = HyperCube<2>::getCombinations_R(1);
+
+	// Check
+	BOOST_REQUIRE_EQUAL(ele[0],4);
+	BOOST_REQUIRE_EQUAL(ele[1],4);
+	BOOST_REQUIRE_EQUAL(ele[2],1);
+
+	// Check combination
+
+	comb<2> c2_0[] = {{1,1},{1,-1},{-1,1},{-1,-1}};
+	comb<2> c2_1[] = {{1,0},{-1,0},{0,1},{0,-1}};
+	check_lin(v_c2_0);
+	check_lin(v_c2_1);
+
+	boost_check_array(&c2_0[0],&v_c2_0[0],4);
+	boost_check_array(&c2_1[0],&v_c2_1[0],4);
+
+	// Cube
+	// Number of vertex
+	ele[0] = HyperCube<3>::getNumberOfElements_R(0);
+	// Number of edge
+	ele[1] = HyperCube<3>::getNumberOfElements_R(1);
+	// Number of faces
+	ele[2] = HyperCube<3>::getNumberOfElements_R(2);
+	// Number of Cubes
+	ele[3] = HyperCube<3>::getNumberOfElements_R(3);
+
+	// Get combination for each dimensions
+	std::vector<comb<3>> v_c3_0 = HyperCube<3>::getCombinations_R(0);
+	std::vector<comb<3>> v_c3_1 = HyperCube<3>::getCombinations_R(1);
+	std::vector<comb<3>> v_c3_2 = HyperCube<3>::getCombinations_R(2);
+
+	// Check
+	BOOST_REQUIRE_EQUAL(ele[0],8);
+	BOOST_REQUIRE_EQUAL(ele[1],12);
+	BOOST_REQUIRE_EQUAL(ele[2],6);
+	BOOST_REQUIRE_EQUAL(ele[3],1);
+
+	// Check combination
+
+	comb<3> c3_0[] = {{1,1,1},{1,1,-1},{1,-1,1},{1,-1,-1},{-1,1,1},{-1,1,-1},{-1,-1,1},{-1,-1,-1}};
+	comb<3> c3_1[] = {{1,1,0},{1,-1,0},{-1,1,0},{-1,-1,0},{1,0,1},{1,0,-1},{-1,0,1},{-1,0,-1},{0,1,1},{0,1,-1},{0,-1,1},{0,-1,-1}};
+	comb<3> c3_2[] = {{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}};
+	check_lin(v_c3_0);
+	check_lin(v_c3_1);
+	check_lin(v_c3_2);
+
+	boost_check_array(&c3_0[0],&v_c3_0[0],8);
+	boost_check_array(&c3_1[0],&v_c3_1[0],12);
+	boost_check_array(&c3_2[0],&v_c3_2[0],6);
+
+	// Tesseract
+	// Number of vertex
+	ele[0] = HyperCube<4>::getNumberOfElements_R(0);
+	// Number of edge
+	ele[1] = HyperCube<4>::getNumberOfElements_R(1);
+	// Number of faces
+	ele[2] = HyperCube<4>::getNumberOfElements_R(2);
+	// Number of Cubes
+	ele[3] = HyperCube<4>::getNumberOfElements_R(3);
+	// Number of Tessaract
+	ele[4] = HyperCube<4>::getNumberOfElements_R(4);
+
+	// Get combination for each dimensions
+	std::vector<comb<4>> v_c4_0 = HyperCube<4>::getCombinations_R(0);
+	std::vector<comb<4>> v_c4_1 = HyperCube<4>::getCombinations_R(1);
+	std::vector<comb<4>> v_c4_2 = HyperCube<4>::getCombinations_R(2);
+	std::vector<comb<4>> v_c4_3 = HyperCube<4>::getCombinations_R(3);
+	check_lin(v_c4_0);
+	check_lin(v_c4_1);
+	check_lin(v_c4_2);
+	check_lin(v_c4_3);
+
+	// Check
+	BOOST_REQUIRE_EQUAL(ele[0],16);
+	BOOST_REQUIRE_EQUAL(ele[1],32);
+	BOOST_REQUIRE_EQUAL(ele[2],24);
+	BOOST_REQUIRE_EQUAL(ele[3],8);
+	BOOST_REQUIRE_EQUAL(ele[4],1);
+
+	// Penteract
+	// Number of vertex
+	ele[0] = HyperCube<5>::getNumberOfElements_R(0);
+	// Number of edge
+	ele[1] = HyperCube<5>::getNumberOfElements_R(1);
+	// Number of faces
+	ele[2] = HyperCube<5>::getNumberOfElements_R(2);
+	// Number of Cubes
+	ele[3] = HyperCube<5>::getNumberOfElements_R(3);
+	// Number of Tessaract
+	ele[4] = HyperCube<5>::getNumberOfElements_R(4);
+	// Number of Penteract
+	ele[5] = HyperCube<5>::getNumberOfElements_R(5);
+
+	// Get combination for each dimensions
+	std::vector<comb<5>> v_c5_0 = HyperCube<5>::getCombinations_R(0);
+	std::vector<comb<5>> v_c5_1 = HyperCube<5>::getCombinations_R(1);
+	std::vector<comb<5>> v_c5_2 = HyperCube<5>::getCombinations_R(2);
+	std::vector<comb<5>> v_c5_3 = HyperCube<5>::getCombinations_R(3);
+	std::vector<comb<5>> v_c5_4 = HyperCube<5>::getCombinations_R(4);
+	check_lin(v_c5_0);
+	check_lin(v_c5_1);
+	check_lin(v_c5_2);
+	check_lin(v_c5_3);
+	check_lin(v_c5_4);
+
+	// Check
+	BOOST_REQUIRE_EQUAL(ele[0],32);
+	BOOST_REQUIRE_EQUAL(ele[1],80);
+	BOOST_REQUIRE_EQUAL(ele[2],80);
+	BOOST_REQUIRE_EQUAL(ele[3],40);
+	BOOST_REQUIRE_EQUAL(ele[4],10);
+	BOOST_REQUIRE_EQUAL(ele[5],1);
+
+
+	// Test SubHypercube 2D and 3D
+
+	std::vector<comb<2>> sc2_0 = HyperCube<2>::getCombinations_R(0);
+
+	for (size_t i = 0 ; i < sc2_0.size() ; i++)
+	{
+		// Expecting one element equal to c2[i]
+		std::vector<comb<2>> combs = SubHyperCube<2,0>::getCombinations_R(sc2_0[i],0);
+		BOOST_REQUIRE_EQUAL(combs.size(),1);
+		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
+		BOOST_REQUIRE_EQUAL(isValid(combs),true);
+		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc2_0[i]),true);
+	}
+
+	std::vector<comb<2>> sc2_1 = HyperCube<2>::getCombinations_R(1);
+
+	for (size_t i = 0 ; i < sc2_1.size() ; i++)
+	{
+		// Expecting two elements, valid, distinct,  sub-decomposition of c2[i]
+		std::vector<comb<2>> combs = SubHyperCube<2,1>::getCombinations_R(sc2_1[i],0);
+		BOOST_REQUIRE_EQUAL(combs.size(),2);
+		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
+		BOOST_REQUIRE_EQUAL(isValid(combs),true);
+		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc2_1[i]),true);
+
+		// Expecting one element, valid, distinct,  sub-decomposition of c2[i]
+		combs = SubHyperCube<2,1>::getCombinations_R(sc2_1[i],1);
+		BOOST_REQUIRE_EQUAL(combs.size(),1);
+		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
+		BOOST_REQUIRE_EQUAL(isValid(combs),true);
+		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc2_1[i]),true);
+	}
+
+	std::vector<comb<2>> sc2_2 = HyperCube<2>::getCombinations_R(2);
+
+	for (size_t i = 0 ; i < sc2_2.size() ; i++)
+	{
+		// Expecting two elements, valid, distinct,  sub-decomposition of sc2_2[i]
+		std::vector<comb<2>> combs = SubHyperCube<2,2>::getCombinations_R(sc2_2[i],0);
+		BOOST_REQUIRE_EQUAL(combs.size(),4);
+		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
+		BOOST_REQUIRE_EQUAL(isValid(combs),true);
+		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc2_2[i]),true);
+
+		// Expecting one element, valid, distinct,  sub-decomposition of c2[i]
+		combs = SubHyperCube<2,2>::getCombinations_R(sc2_2[i],1);
+		BOOST_REQUIRE_EQUAL(combs.size(),4);
+		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
+		BOOST_REQUIRE_EQUAL(isValid(combs),true);
+		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc2_2[i]),true);
+
+		// Expecting one element, valid, distinct,  sub-decomposition of c2[i]
+		combs = SubHyperCube<2,2>::getCombinations_R(sc2_2[i],2);
+		BOOST_REQUIRE_EQUAL(combs.size(),1);
+		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
+		BOOST_REQUIRE_EQUAL(isValid(combs),true);
+		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc2_2[i]),true);
+	}
+
+	////////////// 3D ////////////////
+
+	std::vector<comb<3>> sc3_0 = HyperCube<3>::getCombinations_R(0);
+
+	for (size_t i = 0 ; i < sc3_0.size() ; i++)
+	{
+		// Expecting one element equal to sc3[i]
+		std::vector<comb<3>> combs = SubHyperCube<3,0>::getCombinations_R(sc3_0[i],0);
+		BOOST_REQUIRE_EQUAL(combs.size(),1);
+		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
+		BOOST_REQUIRE_EQUAL(isValid(combs),true);
+		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc3_0[i]),true);
+	}
+
+	std::vector<comb<3>> sc3_1 = HyperCube<3>::getCombinations_R(1);
+
+	for (size_t i = 0 ; i < sc3_1.size() ; i++)
+	{
+		// Expecting two elements, valid, distinct,  sub-decomposition of sc3[i]
+		std::vector<comb<3>> combs = SubHyperCube<3,1>::getCombinations_R(sc3_1[i],0);
+		BOOST_REQUIRE_EQUAL(combs.size(),2);
+		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
+		BOOST_REQUIRE_EQUAL(isValid(combs),true);
+		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc3_1[i]),true);
+
+		// Expecting one element, valid, distinct,  sub-decomposition of c2[i]
+		combs = SubHyperCube<3,1>::getCombinations_R(sc3_1[i],1);
+		BOOST_REQUIRE_EQUAL(combs.size(),1);
+		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
+		BOOST_REQUIRE_EQUAL(isValid(combs),true);
+		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc3_1[i]),true);
+	}
+
+	std::vector<comb<3>> sc3_2 = HyperCube<3>::getCombinations_R(2);
+
+	for (size_t i = 0 ; i < sc3_2.size() ; i++)
+	{
+		// Expecting two elements, valid, distinct,  sub-decomposition of sc3_2[i]
+		std::vector<comb<3>> combs = SubHyperCube<3,2>::getCombinations_R(sc3_2[i],0);
+		BOOST_REQUIRE_EQUAL(combs.size(),4);
+		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
+		BOOST_REQUIRE_EQUAL(isValid(combs),true);
+		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc3_2[i]),true);
+
+		// Expecting one element, valid, distinct,  sub-decomposition of c2[i]
+		combs = SubHyperCube<3,2>::getCombinations_R(sc3_2[i],1);
+		BOOST_REQUIRE_EQUAL(combs.size(),4);
+		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
+		BOOST_REQUIRE_EQUAL(isValid(combs),true);
+		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc3_2[i]),true);
+
+		// Expecting one element, valid, distinct,  sub-decomposition of c2[i]
+		combs = SubHyperCube<3,2>::getCombinations_R(sc3_2[i],2);
+		BOOST_REQUIRE_EQUAL(combs.size(),1);
+		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
+		BOOST_REQUIRE_EQUAL(isValid(combs),true);
+		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc3_2[i]),true);
+	}
+
+/*	comb<5> c0;
+	for (int i = 0 ; i < 5 ; i++)
+		c0.c[i] = 1;
+
+	std::vector<comb<5>> combs_0_0 = SubHyperCube<5,0>::getCombinations_R(c0,0);
+	std::vector<comb<5>> combs_1_0 = SubHyperCube<5,1>::getCombinations_R(c,0);
+	std::vector<comb<5>> combs_1_1 = SubHyperCube<5,1>::getCombinations_R(c,1);
+	std::vector<comb<5>> combs_2_0 = SubHyperCube<5,2>::getCombinations_R(c,0);
+	std::vector<comb<5>> combs_2_1 = SubHyperCube<5,2>::getCombinations_R(c,1);
+	std::vector<comb<5>> combs_2_2 = SubHyperCube<5,2>::getCombinations_R(c,2);
+	std::vector<comb<5>> combs_3_0 = SubHyperCube<5,3>::getCombinations_R(c,0);
+	std::vector<comb<5>> combs_3_1 = SubHyperCube<5,3>::getCombinations_R(c,1);
+	std::vector<comb<5>> combs_3_2 = SubHyperCube<5,3>::getCombinations_R(c,2);
+	std::vector<comb<5>> combs_3_3 = SubHyperCube<5,3>::getCombinations_R(c,3);
+	std::vector<comb<5>> combs_4_0 = SubHyperCube<5,4>::getCombinations_R(c,0);
+	std::vector<comb<5>> combs_4_1 = SubHyperCube<5,4>::getCombinations_R(c,1);
+	std::vector<comb<5>> combs_4_2 = SubHyperCube<5,4>::getCombinations_R(c,2);
+	std::vector<comb<5>> combs_4_3 = SubHyperCube<5,4>::getCombinations_R(c,3);
+	std::vector<comb<5>> combs_4_4 = SubHyperCube<5,4>::getCombinations_R(c,4);
+	std::vector<comb<5>> combs_5_0 = SubHyperCube<5,5>::getCombinations_R(c,0);
+	std::vector<comb<5>> combs_5_1 = SubHyperCube<5,5>::getCombinations_R(c,1);
+	std::vector<comb<5>> combs_5_2 = SubHyperCube<5,5>::getCombinations_R(c,2);
+	std::vector<comb<5>> combs_5_3 = SubHyperCube<5,5>::getCombinations_R(c,3);
+	std::vector<comb<5>> combs_5_4 = SubHyperCube<5,5>::getCombinations_R(c,4);
+	std::vector<comb<5>> combs_5_5 = SubHyperCube<5,5>::getCombinations_R(c,5);*/
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+#endif
diff --git a/src/main.cpp b/src/main.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..70c9befaad1c1ab13979c7e98453fafe9e5e8fd3
--- /dev/null
+++ b/src/main.cpp
@@ -0,0 +1,21 @@
+#include <iostream>
+#include "Graph/CartesianGraphFactory.hpp"
+
+#define BOOST_DISABLE_ASSERTS
+
+
+#define BOOST_TEST_MODULE "C++ test module for OpenFPM_pdata project"
+#include <boost/test/included/unit_test.hpp>
+
+#include <grid_dist.hpp>
+#include "Point_test.hpp"
+#include "Decomposition/CartDecomposition.hpp"
+#include "memory/HeapMemory.hpp"
+#include "Space/Shape/Box.hpp"
+#include "util.hpp"
+
+#include "hypercube_unit_test.hpp"
+#include "CartesianGraphFactory_unit_test.hpp"
+#include "metis_util_unit_test.hpp"
+#include "dec_optimizer_unit_test.hpp"
+
diff --git a/src/metis_util.hpp b/src/metis_util.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0ea26f6d478a05a01c3abfa1f1870d6fffe22480
--- /dev/null
+++ b/src/metis_util.hpp
@@ -0,0 +1,306 @@
+/*
+ * metis_util.hpp
+ *
+ *  Created on: Nov 21, 2014
+ *      Author: Pietro Incardona
+ */
+
+#ifndef METIS_UTIL_HPP
+#define METIS_UTIL_HPP
+
+#include <iostream>
+#include "metis.h"
+#include "VTKWriter.hpp"
+
+/*! \brief Metis graph structure
+ *
+ * Metis graph structure
+ *
+ */
+struct Metis_graph
+{
+	//! The number of vertices in the graph
+	idx_t * nvtxs;
+
+	//! number of balancing constrains
+	//! more practical, are the number of weights for each vertex
+	//! PS even we you specify vwgt == NULL ncon must be set at leat to
+	//! one
+	idx_t * ncon;
+
+	//! For each vertex it store the adjacency lost start for the vertex i
+	idx_t * xadj;
+
+	//! For each vertex it store a list of all neighborhood vertex
+	idx_t * adjncy;
+
+	//! Array that store the weight for each vertex
+	idx_t * vwgt;
+
+	//! Array of the vertex size, basically is the total communication amount
+	idx_t * vsize;
+	//! The weight of the edge
+	idx_t * adjwgt;
+	//! number of part to partition the graph
+	idx_t * nparts;
+	//! Desired weight for each partition (one for each constrain)
+	real_t * tpwgts;
+	//! For each partition load imbalance tollerated
+	real_t * ubvec;
+	//! Additional option for the graph partitioning
+	idx_t * options;
+	//! return the total comunication cost for each partition
+	idx_t * objval;
+	//! Is a output vector containing the partition for each vertex
+	idx_t * part;
+};
+
+//! Balance communication and computation
+#define BALANCE_CC 1
+//! Balance communication computation and memory
+#define BALANCE_CCM 2
+//! Balance computation and comunication and others
+#define BALANCE_CC_O(c) c+1
+
+/*! \brief Helper class to define Metis graph
+ *
+ * \tparam graph structure that store the graph
+ *
+ */
+
+template<typename Graph>
+class Metis
+{
+	// Graph in metis reppresentation
+	Metis_graph Mg;
+
+	// Original graph
+	Graph & g;
+
+	/*! \brief Construct Adjacency list
+	 *
+	 * \param g Graph
+	 *
+	 */
+	void constructAdjList(Graph & g)
+	{
+		// create xadj and adjlist
+
+		Mg.xadj = new idx_t[g.getNVertex()+1];
+		Mg.adjncy = new idx_t[g.getNEdge()];
+
+		//! Get a vertex iterator
+		auto it = g.getVertexIterator();
+
+		//! starting point in the adjacency list
+		size_t prev = 0;
+
+		// actual position
+		size_t id = 0;
+
+		// for each vertex calculate the position of the starting point in the adjacency list
+		while (it.isNext())
+		{
+			// Calculate the starting point in the adjacency list
+			Mg.xadj[id] = prev;
+
+			// Create the adjacency list
+			for (size_t s = 0 ; s < g.getNChilds(it) ; s++)
+			{Mg.adjncy[prev+s] = g.getChild(it,s);}
+
+			// update the position for the next vertex
+			prev += g.getNChilds(it);
+
+			++it;
+			id++;
+		}
+
+		// Fill the last point
+		Mg.xadj[id] = prev;
+	}
+
+public:
+
+	/*! \brief Constructor
+	 *
+	 * Construct a metis graph from Graph_CSR
+	 *
+	 * \param g Graph we want to convert to decompose
+	 * \param nc number of partitions
+	 *
+	 */
+	Metis(Graph & g, size_t nc)
+	:g(g)
+	{
+		// Get the number of vertex
+
+		Mg.nvtxs = new idx_t[1];
+		Mg.nvtxs[0] = g.getNVertex();
+
+		// Set the number of constrains
+
+		Mg.ncon = new idx_t[1];
+		Mg.ncon[0] = 1;
+
+		// construct the adjacency list
+		constructAdjList(g);
+
+		// Set to null the weight of the vertex
+
+		Mg.vwgt = NULL;
+
+		// Put the total communication size to NULL
+
+		Mg.vsize = NULL;
+
+		// Set to null the weight of the edge
+
+		Mg.adjwgt = NULL;
+
+		// Set the total number of partitions
+
+		Mg.nparts = new idx_t[1];
+		Mg.nparts[0] = nc;
+
+		// Set to null the desired weight for each partition (one for each constrain)
+
+		Mg.tpwgts = NULL;
+
+		//! Set to null the partition load imbalance tollerace
+
+		Mg.ubvec = NULL;
+
+		//! Set tp null additional option for the graph partitioning
+
+		Mg.options = NULL;
+
+		//! set the objective value
+		Mg.objval = new idx_t[1];
+
+		//! Is an output vector containing the partition for each vertex
+		Mg.part = new idx_t[g.getNVertex()];
+	}
+
+	/*! \brief destructor
+	 *
+	 * Destructor, It destroy all the memory allocated
+	 *
+	 */
+	~Metis()
+	{
+		// Deallocate the Mg structure
+		if (Mg.nvtxs != NULL)
+		{
+			delete [] Mg.nvtxs;
+		}
+
+		if (Mg.ncon != NULL)
+		{
+			delete [] Mg.ncon;
+		}
+
+		if (Mg.xadj != NULL)
+		{
+			delete [] Mg.xadj;
+		}
+
+		if (Mg.adjncy != NULL)
+		{
+			delete [] Mg.adjncy;
+		}
+		if (Mg.vwgt != NULL)
+		{
+			delete [] Mg.vwgt;
+		}
+		if (Mg.adjwgt != NULL)
+		{
+			delete [] Mg.adjwgt;
+		}
+		if (Mg.nparts != NULL)
+		{
+			delete [] Mg.nparts;
+		}
+		if (Mg.tpwgts != NULL)
+		{
+			delete [] Mg.tpwgts;
+		}
+		if (Mg.ubvec != NULL)
+		{
+			delete [] Mg.ubvec;
+		}
+		if (Mg.options != NULL)
+		{
+			delete [] Mg.options;
+		}
+		if (Mg.objval != NULL)
+		{
+			delete [] Mg.objval;
+		}
+		if (Mg.part != NULL)
+		{
+			delete [] Mg.part;
+		}
+	}
+
+	/*! \brief Decompose the graph
+	 *
+	 * \tparam i which property store the decomposition
+	 *
+	 */
+
+	template<unsigned int i>
+	void decompose()
+	{
+		// Decompose
+		METIS_PartGraphKway(Mg.nvtxs,Mg.ncon,Mg.xadj,Mg.adjncy,Mg.vwgt,Mg.vsize,Mg.adjwgt,
+				            Mg.nparts,Mg.tpwgts,Mg.ubvec,Mg.options,Mg.objval,Mg.part);
+
+		// vertex id
+
+		size_t id = 0;
+
+		// For each vertex store the processor that contain the data
+
+		auto it = g.getVertexIterator();
+
+		while (it.isNext())
+		{
+			g.vertex(it).template get<i>() = Mg.part[id];
+
+			++id;
+			++it;
+		}
+	}
+
+	/*! \brief Decompose the graph
+	 *
+	 * \tparam i which property store the decomposition
+	 *
+	 */
+
+	template<unsigned int i, typename Graph_part>
+	void decompose(Graph_part & gp)
+	{
+		// Decompose
+		METIS_PartGraphKway(Mg.nvtxs,Mg.ncon,Mg.xadj,Mg.adjncy,Mg.vwgt,Mg.vsize,Mg.adjwgt,
+				            Mg.nparts,Mg.tpwgts,Mg.ubvec,Mg.options,Mg.objval,Mg.part);
+
+		// vertex id
+
+		size_t id = 0;
+
+		// For each vertex store the processor that contain the data
+
+		auto it = gp.getVertexIterator();
+
+		while (it.isNext())
+		{
+			gp.vertex(it).template get<i>() = Mg.part[id];
+
+			++id;
+			++it;
+		}
+	}
+};
+
+#endif
diff --git a/src/metis_util_unit_test.hpp b/src/metis_util_unit_test.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..dc5442b67c1bbd4057b4fd53de68915c33714761
--- /dev/null
+++ b/src/metis_util_unit_test.hpp
@@ -0,0 +1,186 @@
+/*
+ * metis_util_unit_test.hpp
+ *
+ *  Created on: Dec 21, 2014
+ *      Author: i-bird
+ */
+
+#ifndef METIS_UTIL_UNIT_TEST_HPP_
+#define METIS_UTIL_UNIT_TEST_HPP_
+
+#include "Graph/CartesianGraphFactory.hpp"
+#include "map_graph.hpp"
+#include "metis_util.hpp"
+
+#undef GS_SIZE
+#define GS_SIZE 8
+
+/*!
+ *
+ * Test node
+ *
+ */
+
+struct nm_v
+{
+	//! The node contain 3 unsigned long integer for communication computation memory and id
+	typedef boost::fusion::vector<float,float,float,size_t,size_t,size_t,size_t,size_t> type;
+
+	typedef typename memory_traits_inte<type>::type memory_int;
+	typedef typename memory_traits_lin<type>::type memory_lin;
+
+	//! type of the positional field
+	typedef float s_type;
+
+	//! Attributes name
+	struct attributes
+	{
+		static const std::string name[];
+	};
+
+	//! The data
+	type data;
+
+	//! computation property id in boost::fusion::vector
+	static const unsigned int x = 0;
+	//! computation property id in boost::fusion::vector
+	static const unsigned int y = 1;
+	//! memory property id in boost::fusion::vector
+	static const unsigned int z = 2;
+	//! computation property id in boost::fusion::vector
+	static const unsigned int communication = 3;
+	//! computation property id in boost::fusion::vector
+	static const unsigned int computation = 4;
+	//! memory property id in boost::fusion::vector
+	static const unsigned int memory = 5;
+	//! memory property id in boost::fusion::vector
+	static const unsigned int id = 6;
+	//! memory property sub_id in boost::fusion::vector
+	static const unsigned int sub_id = 7;
+
+	//! total number of properties boost::fusion::vector
+	static const unsigned int max_prop = 8;
+};
+
+const std::string nm_v::attributes::name[] = {"x","y","z","communication","computation","memory","id","sub_id"};
+
+/*!
+ *
+ * Test node
+ *
+ */
+
+struct nm_e
+{
+	//! The node contain 3 unsigned long integer for comunication computation and memory
+	typedef boost::fusion::vector<size_t> type;
+
+	typedef typename memory_traits_inte<type>::type memory_int;
+	typedef typename memory_traits_lin<type>::type memory_lin;
+
+	//! Attributes name
+	struct attributes
+	{
+		static const std::string name[];
+	};
+
+	//! The data
+	type data;
+
+	//! computation property id in boost::fusion::vector
+	static const unsigned int communication = 0;
+	//! total number of properties boost::fusion::vector
+	static const unsigned int max_prop = 1;
+};
+
+const std::string nm_e::attributes::name[] = {"communication"};
+
+struct nm_part_v
+{
+	//! The node contain 3 unsigned long integer for comunication computation and memory
+	typedef boost::fusion::vector<size_t> type;
+
+	typedef typename memory_traits_inte<type>::type memory_int;
+	typedef typename memory_traits_lin<type>::type memory_lin;
+
+	typedef float s_type;
+
+	//! Attributes name
+	struct attributes
+	{
+		static const std::string name[];
+	};
+
+	//! The data
+
+	type data;
+
+	//! partition id in the boost::fusion::vector
+	static const unsigned int id = 0;
+
+	//! total number of properties
+	static const unsigned int max_prop = 1;
+};
+
+const std::string nm_part_v::attributes::name[] = {"id"};
+
+struct nm_part_e
+{
+	//! The node contain 3 unsigned long integer for comunication computation and memory
+	typedef boost::fusion::vector<> type;
+
+	typedef typename memory_traits_inte<type>::type memory_int;
+	typedef typename memory_traits_lin<type>::type memory_lin;
+
+	//! The data
+
+	type data;
+
+	//! total number of properties
+	static const unsigned int max_prop = 0;
+};
+
+BOOST_AUTO_TEST_SUITE( Metis_test )
+
+BOOST_AUTO_TEST_CASE( Metis_test_use)
+{
+	CartesianGraphFactory<3,Graph_CSR<nm_v,nm_e>> g_factory;
+	CartesianGraphFactory<3,Graph_CSR<nm_part_v,nm_part_e>> g_factory_part;
+
+	// Cartesian grid
+	std::vector<size_t> sz;
+	sz.push_back(GS_SIZE);
+	sz.push_back(GS_SIZE);
+	sz.push_back(GS_SIZE);
+
+	// Box
+	Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
+
+	// Graph to decompose
+
+	Graph_CSR<nm_v,nm_e> g = g_factory.construct<nm_e::communication,float,2,0,1,2>(sz,box);
+
+	// Processor graph
+
+	Graph_CSR<nm_part_v,nm_part_e> gp = g_factory_part.construct<NO_EDGE,float,2>(sz,box);
+
+	// Convert the graph to metis
+
+	Metis<Graph_CSR<nm_v,nm_e>> met(g,16);
+
+	// decompose
+
+	met.decompose<nm_part_v::id>(gp);
+	met.decompose<nm_v::id>();
+
+	// Write the VTK file
+
+	VTKWriter<Graph_CSR<nm_part_v,nm_part_e>> vtk(gp);
+	vtk.write("vtk_partition.vtk");
+	VTKWriter<Graph_CSR<nm_v,nm_e>> vtk2(g);
+	vtk2.write("vtk_partition2.vtk");
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+#endif