diff --git a/openfpm_data.doc b/openfpm_data.doc
index 23ec9b0145beff468c5423d0f0d1218b095d2386..34bdaae0e920e3f2122ce87b00210cff4e2c171a 100644
--- a/openfpm_data.doc
+++ b/openfpm_data.doc
@@ -1407,7 +1407,7 @@ FORMULA_TRANSPARENT    = YES
 # The default value is: NO.
 # This tag requires that the tag GENERATE_HTML is set to YES.
 
-USE_MATHJAX            = NO
+USE_MATHJAX            = YES
 
 # When MathJax is enabled you can set the default output format to be used for
 # the MathJax output. See the MathJax site (see:
diff --git a/src/Grid/comb.hpp b/src/Grid/comb.hpp
index c6ee12eea95a152f99967c20d60acf288a587e1e..a62929a1ff7febca346e2eab9ceed452bc8886bc 100644
--- a/src/Grid/comb.hpp
+++ b/src/Grid/comb.hpp
@@ -250,6 +250,17 @@ struct comb
 		return c;
 	}
 
+	/*! \brief get the combination array pointer
+	 *
+	 * \return an array of char representing the combination
+	 *
+	 */
+
+	inline const char * getComb() const
+	{
+		return c;
+	}
+
 	/*! \brief get the index i of the combination
 	 *
 	 * NOTE: used on expression template
@@ -342,6 +353,27 @@ struct comb
 
 		return str.str();
 	}
+
+	/*! \brief Linearization
+	 *
+	 * From the combination produce a number
+	 *
+	 * \return the number
+	 *
+	 */
+	size_t inline lin() const
+	{
+		size_t ret = 0;
+		size_t accu = 1;
+
+		for (size_t i = 0 ; i < dim ; i++)
+		{
+			ret += (c[i] + 1) * accu;
+			accu *= 3;
+		}
+
+		return ret;
+	}
 };
 
 /*! brief specialization of comb in case of dim 0
@@ -464,6 +496,15 @@ struct comb<0>
 		return 0;
 	}
 
+	/* \brief produce an unique number from the combination
+	 *
+	 *
+	 *
+	 */
+	inline size_t lin()
+	{
+		return 0;
+	}
 };
 
 // create a specialization of std::vector<comb<0>>
diff --git a/src/Grid/grid_sm.hpp b/src/Grid/grid_sm.hpp
index 3b8aed216c5efbfec9ac0876f6f4a7ecf2c2a3d4..7e9d1ae2545d534e9b5a8e6792919d30936295ab 100755
--- a/src/Grid/grid_sm.hpp
+++ b/src/Grid/grid_sm.hpp
@@ -10,10 +10,12 @@
 #include "Space/Shape/Box.hpp"
 #include "grid_key.hpp"
 #include <iostream>
+#include "util/mathutil.hpp"
 
-// Box need the definition of grid_key_dx_r
-
+#define PERIODIC 1
+#define NON_PERIODIC 0
 
+// Box need the definition of grid_key_dx_r
 #define HARDWARE 1
 
 /*! \brief Class to check if the edge can be created or not
@@ -344,26 +346,44 @@ public:
 	 * \tparam check class that check the linearization, if this check fail the function return -1
 	 * \param gk grid_key_dx to linearize
 	 * \param sum_id shift on each dimension
+	 * \param bc boundary conditions
 	 *
 	 * \return The linearization of the gk key shifted by c, or -1 if the check fail
 	 */
 
-	template<typename check=NoCheck> inline mem_id LinId(const grid_key_dx<N> & gk, char sum_id[N]) const
+	template<typename check=NoCheck> inline mem_id LinId(const grid_key_dx<N> & gk, const char sum_id[N], const size_t (&bc)[N]) const
 	{
+		mem_id lid;
+
 		// Check the sum produce a valid key
 
-		if (check::valid(gk.k[0] + sum_id[0],sz[0]) == false)
-			return -1;
+		if (bc[0] == NON_PERIODIC)
+		{
+			if (check::valid(gk.k[0] + sum_id[0],sz[0]) == false)
+				return -1;
+
+			lid = gk.k[0] + sum_id[0];
+		}
+		else
+		{
+			lid = openfpm::math::positive_modulo(gk.k[0] + sum_id[0],sz[0]);
+		}
 
-		mem_id lid = gk.k[0] + sum_id[0];
 		for (mem_id i = 1 ; i < N ; i++)
 		{
 			// Check the sum produce a valid key
 
-			if (check::valid(gk.k[i] + sum_id[i],sz[i]) == false)
-				return -1;
+			if (bc[i] == NON_PERIODIC)
+			{
+				if (check::valid(gk.k[i] + sum_id[i],sz[i]) == false)
+					return -1;
 
-			lid += (gk.k[i] + sum_id[i]) * sz_s[i-1];
+				lid += (gk.k[i] + sum_id[i]) * sz_s[i-1];
+			}
+			else
+			{
+				lid += (openfpm::math::positive_modulo(gk.k[i] + sum_id[i],sz[i])) * sz_s[i-1];
+			}
 		}
 
 		return lid;
@@ -700,701 +720,7 @@ public:
 	template <unsigned int,typename> friend class grid_sm;
 };
 
-/**
- *
- * Grid key class iterator, iterate through the grid element
- *
- * \param dim dimensionality of the grid
- *
- * \note if you have a grid you can get this object from getIterator()
- *
- * ### Grid iterator declaration and usage
- * \snippet grid_unit_tests.hpp Grid iterator test usage
- *
- */
-
-template<unsigned int dim>
-class grid_key_dx_iterator
-{
-#ifdef DEBUG
-	// Actual status of the iterator, when the iterator is not initialized cannot be used
-	// and reinitialize must be called
-	bool initialized = false;
-#endif
-
-	grid_sm<dim,void> grid_base;
-
-	/*! \brief return the index i of the gk key
-	 *
-	 * \param i index to get
-	 *
-	 * \return index value
-	 *
-	 */
-
-	size_t get_gk(size_t i) const
-	{
-		return gk.get(i);
-	}
-
-protected:
-
-	grid_key_dx<dim> gk;
-
-public:
-
-	/*! \brief Default constructor
-	 *
-	 * \warning entremly unsafe
-	 * Before use the iterator you have call reinitialize
-	 *
-	 */
-	grid_key_dx_iterator()
-	{
-#ifdef DEBUG
-		initialized = false;
-#endif
-	}
-
-	/*! \brief Constructor from a grid_key_dx_iterator<dim>
-	 *
-	 * \param g_it grid_key_dx_iterator<dim>
-	 */
-	grid_key_dx_iterator(const grid_key_dx_iterator<dim> & g_it)
-	: grid_base(g_it.grid_base)
-	{
-		//! Initialize to 0 the index
-
-		for (size_t i = 0 ; i < dim ; i++)
-		{
-			gk.set_d(i,g_it.get_gk(i));
-		}
-
-#ifdef DEBUG
-		initialized = true;
-#endif
-	}
-
-	/*! \brief Constructor require a grid_sm<dim,T>
-	 *
-	 * \param g info of the grid on which iterate
-	 */
-	template<typename T> grid_key_dx_iterator(const grid_sm<dim,T> & g)
-	: grid_base(g)
-	{
-		reset();
-
-#ifdef DEBUG
-		initialized = true;
-#endif
-	}
-
-	/*! \brief Constructor from another grid_key_dx_iterator
-	 *
-	 * \param key_it grid_key_dx_iterator
-	 */
-	grid_key_dx_iterator<dim> operator=(const grid_key_dx_iterator<dim> & key_it)
-	{
-		grid_base = key_it.grid_base;
-
-		//! Initialize the index using key_it
-
-		for (size_t i = 0 ; i < dim ; i++)
-		{gk.set_d(i,key_it.get_gk(i));}
-
-		return *this;
-	}
-
-	/*! \brief Get the next element
-	 *
-	 * \return the next grid_key
-	 *
-	 */
-
-	grid_key_dx_iterator<dim> & operator++()
-	{
-		//! increment the first index
-
-		size_t id = gk.get(0);
-		gk.set_d(0,id+1);
-
-		//! check the overflow of all the index with exception of the last dimensionality
-
-		size_t i = 0;
-		for ( ; i < dim-1 ; i++)
-		{
-			size_t id = gk.get(i);
-			if (id >= grid_base.size(i))
-			{
-				// ! overflow, increment the next index
-
-				gk.set_d(i,0);
-				id = gk.get(i+1);
-				gk.set_d(i+1,id+1);
-			}
-			else
-			{
-				break;
-			}
-		}
-
-		return *this;
-	}
-
-	/*! \brief Set the dimension
-	 *
-	 * Set the dimension
-	 *
-	 * \param d is the dimension
-	 * \param sz set the counter to sz
-	 *
-	 */
-	void set(int d, size_t sz)
-	{
-		// set the counter dim to sz
-
-		gk.set_d(d,sz);
-	}
-
-	/*! \brief Check if there is the next element
-	 *
-	 * Check if there is the next element
-	 *
-	 * \return true if there is the next, false otherwise
-	 *
-	 */
-
-	bool isNext()
-	{
-		if (gk.get(dim-1) < (long int)grid_base.size(dim-1))
-		{
-			//! we did not reach the end of the grid
-
-			return true;
-		}
-
-		//! we reach the end of the grid
-		return false;
-	}
-
-	/*! \brief Get the actual key
-	 *
-	 * Get the actual key
-	 *
-	 * \return the actual key
-	 *
-	 */
-	const grid_key_dx<dim> & get()
-	{
-		return gk;
-	}
-
-	/*! \brief Reinitialize the grid_key_dx_iterator
-	 *
-	 * \param key form
-	 *
-	 */
-	void reinitialize(const grid_key_dx_iterator<dim> & key)
-	{
-		grid_base = key.grid_base;
-		reset();
-	}
-
-	/*! \brief Reset the iterator (it restart from the beginning)
-	 *
-	 */
-	void reset()
-	{
-		//! Initialize to 0 the index
-
-		for (size_t i = 0 ; i < dim ; i++)
-		{gk.set_d(i,0);}
-	}
-};
-
-
-/**
- *
- * Grid key class iterator, iterate through a starting linearized grid element
- * to a stop linearized grid element in particular if linearize is the function
- *  that linearize all the grid_key_dx, it create an iterator that pass through
- *  Linearize^(-1)(start) Linearize^(-1)(start+1) ....... Linearize^(-1)(stop)
- *
- * \param dim dimensionality of the grid
- *
- * Usage: In general you never create object directly, but you get it from a grid_cpu or grid_gpu with
- *        getIteratorLinStartStop()
- *
- */
-
-template<unsigned int dim>
-class grid_key_dx_iterator_sp : public grid_key_dx_iterator<dim>
-{
-	//! stop point
-	grid_key_dx<dim> gk_stop;
-
-public:
-
-	/*! \brief Constructor require a grid grid<dim,T>
-	 *
-	 * It construct an iterator from one index to another, in particular
-	 * if linearize is the function that linearize all the grid_key, it
-	 * create an iterator that pass through Linearize^(-1)(start)
-	 * Linearize^(-1)(start+1) ....... Linearize^(-1)(stop)
-	 *
-	 * For example for start (1,1) and stop (3,3) the point indicated with # are
-	 * explored by the iterator
-	 *
-	 * \verbatim
-	 *
-                      +-----+-----+-----+-----+-----+-----+ (6,5)
-                      |     |     |     |     |     |     |
-                      |     |     |     |     |     |     |
-                      |     |     |     |     |     |     |
-                      +-----+-----+-----+-----+-----+-----+
-                      |     |     |     |     |     |     |
-                      |     |     |     |     |     |     |
-                      |     |     |     |     |     |     |
-                      #-----#-----#-----#-----+-----+-----+
-                      |     |     |     |     |     |     |
-                      |     |     |     |     |     |     |
-                      |     |     |     |     |     |     |
-                      #-----#-----#-----#-----#-----#-----#
-                      |     |     |     |     |     |     |
-                      |     |     |     |     |     |     |
-                      +-----#-----#-----#-----#-----#-----#
-                      |     |     |     |     |     |     |
-                      |     |     |     |     |     |     |
-                      |     |     |     |     |     |     |
-                      +-----+-----+-----+-----+-----+-----+
-                    (0,0)
-	 *
-	 *
-	 * \endverbatim
-	 *
-	 *
-	 * \tparam T type of object that the grid store
-	 *
-	 * \param g Grid on which iterate
-	 * \param from starting point
-	 * \param to end point
-	 *
-	 */
-	template<typename T> grid_key_dx_iterator_sp(grid_sm<dim,T> & g, mem_id from, mem_id to)
-	:grid_key_dx_iterator<dim>(g)
-	 {
-		//! Convert to a grid_key
-		this->gk = g.InvLinId(from);
-
-		//! Convert to a grid_key
-		gk_stop = g.InvLinId(to);
-	 }
-
-	/*! \brief Check if there is the next element
-	 *
-	 * Check if there is the next element
-	 *
-	 * \return true if there is the next, false otherwise
-	 *
-	 */
-
-	bool isNext()
-	{
-		//! for all dimensions
-		for (int i = dim-1 ; i >= 0 ; i-- )
-		{
-			//! check if we still have points
-			if (this->gk.get(i) < gk_stop.get(i))
-				return true;
-			else if (this->gk.get(i) > gk_stop.get(i))
-				return false;
-		}
-
-		//! (Final point) we we still have one point
-		return true;
-	}
-};
-
-/* \brief grid_key_dx_iterator_sub can adjust the domain if the border go out-of-side
- *        in this case a warning is produced
- *
- * \tparam dim dimensionality
- *
- */
-template <unsigned int dim>
-class print_warning_on_adjustment
-{
-public:
-	inline static void pw1(size_t i, const grid_key_dx<dim> & gk_start) {std::cerr << "Warning: " << __FILE__ << ":" << __LINE__ << " start with index smaller than zero x[" << i << "]=" << gk_start.get(i) << "\n";}
-	inline static void pw2(size_t i, const grid_key_dx<dim> & gk_start) {std::cerr << "Warning: " << __FILE__ << ":" << __LINE__ << " Cropping start point x[" << i << "]=" << gk_start.get(i) << " to x[" << i << "]=0 \n"  ;}
-	inline static void pw3() {std::cerr << "Warning: " << __FILE__ << ":" << __LINE__ << " stop with smaller index than start \n";}
-	inline static void pw4(size_t i, const grid_key_dx<dim> & gk_stop, const grid_sm<dim,void> & grid_base) {std::cerr << "Warning: " << __FILE__ << ":" << __LINE__ << " stop index bigger than cell domain x[" << i << "]=" << gk_stop.get(i) << " >= " << grid_base.size(i) << "\n";}
-	inline static void pw5() {std::cerr << "Warning grid_key_dx_iterator_sub: " << __FILE__ << ":" << __LINE__ << " the starting point is not smaller or equal than than the stop point in all the coordinates" << "\n";}
-};
-
-/* \brief grid_key_dx_iterator_sub can adjust the domain if the border go out-of-side
- *        in this case a warning is produced
- *
- * \tparam dim dimensionality
- * \tparam gb type of grid
- *
- */
-template <unsigned int dim>
-class do_not_print_warning_on_adjustment
-{
-public:
-	inline static void pw1(size_t i, const grid_key_dx<dim> & gk_start) {}
-	inline static void pw2(size_t i, const grid_key_dx<dim> & gk_start) {}
-	inline static void pw3() {}
-	inline static void pw4(size_t i, const grid_key_dx<dim> & gk_stop, const grid_sm<dim,void> & grid_base) {}
-	inline static void pw5() {}
-};
-
-/**
- *
- * Grid key class iterator, iterate through a sub-grid defined by an hyper-cube
- *
- * \param dim dimensionality of the grid
- *
- * \note if you have a grid you can get this object from getIteratorSub()
- *
- * ### Sub grid iterator declaration and usage
- * \snippet grid_unit_tests.hpp Sub-grid iterator test usage
- *
- */
-
-template<unsigned int dim, typename warn>
-class grid_key_dx_iterator_sub : public grid_key_dx_iterator<dim>
-{
-#ifdef DEBUG
-	bool initialized = false;
-#endif
-
-	//! grid base where we are iterating
-	grid_sm<dim,void> grid_base;
-
-	//! start point
-	grid_key_dx<dim> gk_start;
-
-	//! stop point
-	grid_key_dx<dim> gk_stop;
-
-	/*! \brief Initialize gk
-	 *
-	 */
-	void Initialize()
-	{
-		// Check that start and stop are inside the domain otherwise crop them
-
-		for (size_t i = 0 ; i < dim ; i++)
-		{
-			// if start smaller than 0
-			if (gk_start.get(i) < 0)
-			{
-#ifdef DEBUG
-				warn::pw1(i,gk_start);
-#endif
-				if (gk_start.get(i) < gk_stop.get(i))
-				{
-#ifdef DEBUG
-					warn::pw2(i,gk_start);
-#endif
-					gk_start.set_d(i,0);
-				}
-				else
-				{
-#ifdef DEBUG
-					warn::pw3();
-#endif
-					// No points are available
-					gk_start.set_d(dim-1,gk_stop.get(dim-1)+1);
-					break;
-				}
-			}
-
-			// if stop bigger than the domain
-			if (gk_stop.get(i) >= (long int)grid_base.size(i))
-			{
-#ifdef DEBUG
-				warn::pw4(i,gk_stop,grid_base);
-#endif
-
-				if (gk_start.get(i) < (long int)grid_base.size(i))
-					gk_stop.set_d(i,grid_base.size(i)-1);
-				else
-				{
-					// No point are available
-					gk_start.set_d(dim-1,gk_stop.get(dim-1)+1);
-					break;
-				}
-			}
-		}
-
-		//! Initialize gk
-		for (unsigned int i = 0 ; i < dim ; i++)
-		{
-			this->gk.set_d(i,gk_start.get(i));
-		}
-
-#ifdef DEBUG
-		initialized = true;
-#endif
-	}
-
-public:
-
-	/*! \brief Default constructor
-	 *
-	 * \warning extremely unsafe
-	 * If you use this constructor before use the iterator you should call reinitialize first
-	 *
-	 */
-	grid_key_dx_iterator_sub()
-{}
-
-	/*! \brief Constructor from another grid_key_dx_iterator_sub
-	 *
-	 * It construct an iterator over an hyper-cube defined by start and stop,
-	 * \warning if start and stop are outside the domain defined by g the intersection
-	 * will be considered
-	 *
-	 * \param g_s_it grid_key_dx_iterator_sub
-	 *
-	 */
-	grid_key_dx_iterator_sub(const grid_key_dx_iterator_sub<dim> & g_s_it)
-	: grid_key_dx_iterator_sub(g_s_it.grid_base,g_s_it.gk_start,g_s_it.gk_stop)
-	{}
-
-
-	/*! \brief Constructor require a grid grid<dim,T>
-	 *
-	 * It construct an iterator over an hyper-cube defined by start and stop,
-	 * \warning if start and stop are outside the domain defined by g the intersection
-	 * will be considered
-	 *
-	 * \tparam T type of object that the grid store
-	 *
-	 * \param g Grid on which iterate
-	 * \param start starting point
-	 * \param stop end point
-	 *
-	 */
-	template<typename T> grid_key_dx_iterator_sub(const grid_sm<dim,T> & g, const grid_key_dx<dim> & start, const grid_key_dx<dim> & stop)
-	: grid_key_dx_iterator<dim>(g),grid_base(g),gk_start(start), gk_stop(stop)
-	{
-#ifdef DEBUG
-		//! If we are on debug check that the stop grid_key id bigger than the start
-		//! grid_key
-
-		for (unsigned int i = 0 ; i < dim ; i++)
-		{
-			if (gk_start.get(i) > gk_stop.get(i))
-			{
-				warn::pw5();
-			}
-		}
-#endif
-
-		Initialize();
-	}
-
-
-	/*! \brief Constructor require a grid grid<dim,T>
-	 *
-	 * It construct an iterator over an hyper-cube defined by start and stop,
-	 * \warning if start and stop are outside the domain defined by g the intersection
-	 * will be considered
-	 *
-	 * \tparam T type of object that the grid store
-	 *
-	 * \param g info of the grid where we are iterating
-	 * \param m Margin of the domain
-	 *
-	 */
-	template<typename T> grid_key_dx_iterator_sub(const grid_sm<dim,T> & g, const size_t m)
-			: grid_key_dx_iterator<dim>(g),grid_base(g)
-			  {
-		// Initialize the start and stop point
-		for (unsigned int i = 0 ; i < dim ; i++)
-		{
-			gk_start.set(i,m);
-			gk_stop.set(i,g.size(i)-m);
-		}
-
-		//
-		Initialize();
-			  }
-
-	/*! \brief Constructor require a grid grid<dim,T>
-	 *
-	 * It construct an iterator over an hyper-cube defined by start and stop,
-	 *
-	 * \tparam T type of object that the grid store
-	 *
-	 * \param g Grid on which iterate
-	 * \param start starting point
-	 * \param stop end point
-	 *
-	 */
-	template<typename T> grid_key_dx_iterator_sub(const grid_sm<dim,T> & g, const size_t (& start)[dim], const size_t (& stop)[dim])
-			: grid_key_dx_iterator<dim>(g),grid_base(g),gk_start(start), gk_stop(stop)
-			  {
-#ifdef DEBUG
-		//! If we are on debug check that the stop grid_key id bigger than the start
-		//! grid_key
-
-		for (unsigned int i = 0 ; i < dim ; i++)
-		{
-			if (start[i] > stop[i])
-			{
-				warn::pw5();
-			}
-		}
-#endif
-
-		Initialize();
-			  }
-
-	/*! \brief Get the next element
-	 *
-	 * Get the next element
-	 *
-	 * \return the next grid_key
-	 *
-	 */
-
-	grid_key_dx_iterator<dim> & operator++()
-			{
-#ifdef DEBUG
-		if (initialized == false)
-		{std::cerr << "Error: " << __FILE__ << __LINE__ << " using unitialized iterator" << "\n";}
-#endif
-
-		//! increment the first index
-
-		size_t id = this->gk.get(0);
-		this->gk.set_d(0,id+1);
-
-		//! check the overflow of all the index with exception of the last dimensionality
-
-		size_t i = 0;
-		for ( ; i < dim-1 ; i++)
-		{
-			size_t id = this->gk.get(i);
-			if ((long int)id > gk_stop.get(i))
-			{
-				// ! overflow, increment the next index
-
-				this->gk.set_d(i,gk_start.get(i));
-				id = this->gk.get(i+1);
-				this->gk.set_d(i+1,id+1);
-			}
-			else
-			{
-				break;
-			}
-		}
-
-		return *this;
-			}
-
-	/*! \brief Check if there is the next element
-	 *
-	 * Check if there is the next element
-	 *
-	 * \return true if there is the next, false otherwise
-	 *
-	 */
-
-	inline bool isNext()
-	{
-#ifdef DEBUG
-		if (initialized == false)
-		{std::cerr << "Error: " << __FILE__ << __LINE__ << " using unitialized iterator" << "\n";}
-#endif
 
-		if (this->gk.get(dim-1) <= gk_stop.get(dim-1))
-		{
-			//! we did not reach the end of the grid
-
-			return true;
-		}
-
-		//! we reach the end of the grid
-		return false;
-	}
-
-	/*! \brief Return the actual grid key iterator
-	 *
-	 */
-	inline grid_key_dx<dim> get()
-	{
-#ifdef DEBUG
-		if (initialized == false)
-		{std::cerr << "Error: " << __FILE__ << __LINE__ << " using unitialized iterator" << "\n";}
-#endif
-
-		return grid_key_dx_iterator<dim>::get();
-	}
-
-	/*! \brief Reinitialize the iterator
-	 *
-	 * it re-initialize the iterator with the passed grid_key_dx_iterator_sub
-	 * the actual position of the grid_key_dx_iterator_sub is ignored
-	 *
-	 * \param g_s_it grid_key_dx_iterator_sub
-	 *
-	 */
-
-	inline void reinitialize(const grid_key_dx_iterator_sub<dim> & g_s_it)
-	{
-		// Reinitialize the iterator
-
-		grid_key_dx_iterator<dim>::reinitialize(g_s_it);
-		grid_base = g_s_it.grid_base;
-		gk_start = g_s_it.gk_start;
-		gk_stop = g_s_it.gk_stop;
-
-
-#ifdef DEBUG
-		//! If we are on debug check that the stop grid_key id bigger than the start
-		//! grid_key
-
-		for (unsigned int i = 0 ; i < dim ; i++)
-		{
-			if (gk_start.get(i) > gk_stop.get(i))
-			{
-				warn::pw5();
-			}
-		}
-
-		initialized = true;
-#endif
-
-		Initialize();
-	}
-
-	/*! \brief Get the volume spanned by this sub-grid iterator
-	 *
-	 * \return the volume
-	 *
-	 */
-	inline size_t getVolume()
-	{
-		return Box<dim,long int>::getVolumeKey(gk_start.k, gk_stop.k);
-	}
-
-	/*! \brief Reset the iterator (it restart from the beginning)
-	 *
-	 */
-	inline void reset()
-	{
-		//! Initialize to 0 the index
-
-		for (size_t i = 0 ; i < dim ; i++)
-		{this->gk.set_d(i,gk_start.get(i));}
-	}
-};
 
 
 /*! \brief Emulate grid_key_dx with runtime dimensionality
diff --git a/src/Grid/grid_sm_unit_tests.hpp b/src/Grid/grid_sm_unit_tests.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a273c479bce3f315bd7a9489ebfadbf5f56ad969
--- /dev/null
+++ b/src/Grid/grid_sm_unit_tests.hpp
@@ -0,0 +1,174 @@
+/*
+ * grid_sm_test.hpp
+ *
+ *  Created on: Dec 14, 2015
+ *      Author: i-bird
+ */
+
+#ifndef OPENFPM_DATA_SRC_GRID_GRID_SM_UNIT_TESTS_HPP_
+#define OPENFPM_DATA_SRC_GRID_GRID_SM_UNIT_TESTS_HPP_
+
+#include "iterators/grid_key_dx_iterator_sub_bc.hpp"
+
+BOOST_AUTO_TEST_SUITE( grid_sm_test )
+
+
+BOOST_AUTO_TEST_CASE( grid_sm_linearization )
+{
+	const grid_key_dx<3> key1(1,2,3);
+	const grid_key_dx<3> zero(0,0,0);
+	const grid_key_dx<3> seven(7,7,7);
+	const comb<3> c({1,0,-1});
+	size_t sz[3] = {8,8,8};
+
+	grid_sm<3,int> gs(sz);
+	size_t bc[] = {NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
+
+	long int lin = gs.LinId<CheckExistence>(key1,c.getComb(),bc);
+	BOOST_REQUIRE_EQUAL(lin,146);
+	lin = gs.LinId<CheckExistence>(zero,c.getComb(),bc);
+	BOOST_REQUIRE_EQUAL(lin,-1);
+	lin = gs.LinId<CheckExistence>(seven,c.getComb(),bc);
+	BOOST_REQUIRE_EQUAL(lin,-1);
+
+	for (size_t i = 0 ; i < 3 ; i++)
+		bc[i] = PERIODIC;
+
+	lin = gs.LinId<CheckExistence>(key1,c.getComb(),bc);
+	BOOST_REQUIRE_EQUAL(lin,146);
+	lin = gs.LinId<CheckExistence>(zero,c.getComb(),bc);
+	BOOST_REQUIRE_EQUAL(lin,71);
+	lin = gs.LinId<CheckExistence>(seven,c.getComb(),bc);
+	BOOST_REQUIRE_EQUAL(lin,62);
+}
+
+
+BOOST_AUTO_TEST_CASE( grid_iterator_sub_p )
+{
+	const grid_key_dx<3> key1(4,4,4);
+	const grid_key_dx<3> key2(-1,-1,-1);
+	const grid_key_dx<3> key3(9,9,9);
+	size_t sz[3] = {8,8,8};
+
+	grid_sm<3,int> gs(sz);
+
+	grid_key_dx_iterator_sub_bc<3> it(gs,key2,key1,{PERIODIC,PERIODIC,PERIODIC});
+
+	size_t cnt = 0;
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		for (size_t i = 0 ; i < 3 ; i++)
+		{
+			BOOST_REQUIRE_EQUAL(key.get(i) >= (long int)0,true);
+			BOOST_REQUIRE_EQUAL(key.get(i) < (long int)sz[i],true);
+		}
+
+		cnt++;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(cnt,216ul);
+
+	grid_key_dx_iterator_sub_bc<3> it2(gs,key2,key3,{PERIODIC,PERIODIC,PERIODIC});
+
+	cnt = 0;
+
+	while (it2.isNext())
+	{
+		auto key = it2.get();
+
+		for (size_t i = 0 ; i < 3 ; i++)
+		{
+			BOOST_REQUIRE_EQUAL(key.get(i) >= (long int)0,true);
+			BOOST_REQUIRE_EQUAL(key.get(i) < (long int)sz[i],true);
+		}
+
+		cnt++;
+
+		++it2;
+	}
+
+	BOOST_REQUIRE_EQUAL(cnt,1331ul);
+
+	cnt = 0;
+
+	const grid_key_dx<3> key4(0,-1,0);
+	const grid_key_dx<3> key5(2,2,2);
+
+	grid_key_dx_iterator_sub_bc<3> it3(gs,key4,key5,{NON_PERIODIC,PERIODIC,NON_PERIODIC});
+
+	while (it3.isNext())
+	{
+		auto key = it3.get();
+
+		for (size_t i = 0 ; i < 3 ; i++)
+		{
+			BOOST_REQUIRE_EQUAL(key.get(i) >= (long int)0,true);
+			BOOST_REQUIRE_EQUAL(key.get(i) < (long int)sz[i],true);
+		}
+
+		cnt++;
+
+		++it3;
+	}
+
+	BOOST_REQUIRE_EQUAL(cnt,36ul);
+
+	// bc non periodic with out-of-bound
+
+	grid_key_dx_iterator_sub_bc<3,do_not_print_warning_on_adjustment<3>> it4(gs,key4,key5,{NON_PERIODIC,NON_PERIODIC,NON_PERIODIC});
+
+	cnt = 0;
+
+	while (it4.isNext())
+	{
+		auto key = it4.get();
+
+		for (size_t i = 0 ; i < 3 ; i++)
+		{
+			BOOST_REQUIRE_EQUAL(key.get(i) >= (long int)0,true);
+			BOOST_REQUIRE_EQUAL(key.get(i) < (long int)sz[i],true);
+		}
+
+		cnt++;
+
+		++it4;
+	}
+
+	BOOST_REQUIRE_EQUAL(cnt,27ul);
+
+	// case with no key
+
+	const grid_key_dx<3> key6(-1,-1,-1);
+	const grid_key_dx<3> key7(-1,-1,8);
+
+	grid_key_dx_iterator_sub_bc<3,do_not_print_warning_on_adjustment<3>> it5(gs,key6,key7,{NON_PERIODIC,NON_PERIODIC,NON_PERIODIC});
+
+	cnt = 0;
+
+	while (it5.isNext())
+	{
+		auto key = it5.get();
+
+		for (size_t i = 0 ; i < 3 ; i++)
+		{
+			BOOST_REQUIRE_EQUAL(key.get(i) >= 0,true);
+			BOOST_REQUIRE_EQUAL(key.get(i) < (long int)sz[i],true);
+		}
+
+		cnt++;
+
+		++it5;
+	}
+
+	BOOST_REQUIRE_EQUAL(cnt,0ul);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+
+#endif /* OPENFPM_DATA_SRC_GRID_GRID_SM_UNIT_TESTS_HPP_ */
diff --git a/src/Grid/grid_unit_tests.hpp b/src/Grid/grid_unit_tests.hpp
index b8414de7fba54ca591909c280165d3ffb7fb1f63..790c960cba16d2e41ab2b28c92c269decc13a18b 100644
--- a/src/Grid/grid_unit_tests.hpp
+++ b/src/Grid/grid_unit_tests.hpp
@@ -539,7 +539,7 @@ BOOST_AUTO_TEST_CASE( grid_iterator_test_use)
 		++g_it;
 	}
 
-	BOOST_REQUIRE_EQUAL(count, 16*16*16);
+	BOOST_REQUIRE_EQUAL(count, (size_t)16*16*16);
 	//! [Grid iterator test usage]
 	}
 
@@ -573,7 +573,7 @@ BOOST_AUTO_TEST_CASE( grid_iterator_test_use)
 		++g_it;
 	}
 
-	BOOST_REQUIRE_EQUAL(count, 14*14*14);
+	BOOST_REQUIRE_EQUAL(count, (size_t)14*14*14);
 
 	//! [Sub-grid iterator test usage]
 
@@ -618,7 +618,7 @@ BOOST_AUTO_TEST_CASE( grid_sub_iterator_test )
 		++g_it;
 	}
 
-	BOOST_REQUIRE_EQUAL(count, 14*14*14);
+	BOOST_REQUIRE_EQUAL(count, (size_t)14*14*14);
 
 	//! [Sub-grid iterator test usage]
 }
@@ -883,7 +883,7 @@ BOOST_AUTO_TEST_CASE( grid_iterator_sp_test )
 		++it;
 	}
 
-	BOOST_REQUIRE_EQUAL(count,2185);
+	BOOST_REQUIRE_EQUAL(count,2185ul);
 }
 
 /* \brief This is an ordinary test simple 3D with plain C array
diff --git a/src/Grid/iterators/grid_key_dx_iterator.hpp b/src/Grid/iterators/grid_key_dx_iterator.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f730b81d5ad99f73a54d59d5be9caa9d811e91a7
--- /dev/null
+++ b/src/Grid/iterators/grid_key_dx_iterator.hpp
@@ -0,0 +1,226 @@
+/*
+ * grid_key_dx_iterator.hpp
+ *
+ *  Created on: Dec 15, 2015
+ *      Author: i-bird
+ */
+
+#ifndef OPENFPM_DATA_SRC_GRID_ITERATORS_GRID_KEY_DX_ITERATOR_HPP_
+#define OPENFPM_DATA_SRC_GRID_ITERATORS_GRID_KEY_DX_ITERATOR_HPP_
+
+
+
+/**
+ *
+ * Grid key class iterator, iterate through the grid element
+ *
+ * \param dim dimensionality of the grid
+ *
+ * \note if you have a grid you can get this object from getIterator()
+ *
+ * ### Grid iterator declaration and usage
+ * \snippet grid_unit_tests.hpp Grid iterator test usage
+ *
+ */
+
+template<unsigned int dim>
+class grid_key_dx_iterator
+{
+#ifdef DEBUG
+	// Actual status of the iterator, when the iterator is not initialized cannot be used
+	// and reinitialize must be called
+	bool initialized = false;
+#endif
+
+	grid_sm<dim,void> grid_base;
+
+	/*! \brief return the index i of the gk key
+	 *
+	 * \param i index to get
+	 *
+	 * \return index value
+	 *
+	 */
+
+	size_t get_gk(size_t i) const
+	{
+		return gk.get(i);
+	}
+
+protected:
+
+	grid_key_dx<dim> gk;
+
+public:
+
+	/*! \brief Default constructor
+	 *
+	 * \warning entremly unsafe
+	 * Before use the iterator you have call reinitialize
+	 *
+	 */
+	grid_key_dx_iterator()
+	{
+#ifdef DEBUG
+		initialized = false;
+#endif
+	}
+
+	/*! \brief Constructor from a grid_key_dx_iterator<dim>
+	 *
+	 * \param g_it grid_key_dx_iterator<dim>
+	 */
+	grid_key_dx_iterator(const grid_key_dx_iterator<dim> & g_it)
+	: grid_base(g_it.grid_base)
+	{
+		//! Initialize to 0 the index
+
+		for (size_t i = 0 ; i < dim ; i++)
+		{
+			gk.set_d(i,g_it.get_gk(i));
+		}
+
+#ifdef DEBUG
+		initialized = true;
+#endif
+	}
+
+	/*! \brief Constructor require a grid_sm<dim,T>
+	 *
+	 * \param g info of the grid on which iterate
+	 */
+	template<typename T> grid_key_dx_iterator(const grid_sm<dim,T> & g)
+	: grid_base(g)
+	{
+		reset();
+
+#ifdef DEBUG
+		initialized = true;
+#endif
+	}
+
+	/*! \brief Constructor from another grid_key_dx_iterator
+	 *
+	 * \param key_it grid_key_dx_iterator
+	 */
+	grid_key_dx_iterator<dim> operator=(const grid_key_dx_iterator<dim> & key_it)
+	{
+		grid_base = key_it.grid_base;
+
+		//! Initialize the index using key_it
+
+		for (size_t i = 0 ; i < dim ; i++)
+		{gk.set_d(i,key_it.get_gk(i));}
+
+		return *this;
+	}
+
+	/*! \brief Get the next element
+	 *
+	 * \return the next grid_key
+	 *
+	 */
+
+	grid_key_dx_iterator<dim> & operator++()
+	{
+		//! increment the first index
+
+		size_t id = gk.get(0);
+		gk.set_d(0,id+1);
+
+		//! check the overflow of all the index with exception of the last dimensionality
+
+		size_t i = 0;
+		for ( ; i < dim-1 ; i++)
+		{
+			size_t id = gk.get(i);
+			if (id >= grid_base.size(i))
+			{
+				// ! overflow, increment the next index
+
+				gk.set_d(i,0);
+				id = gk.get(i+1);
+				gk.set_d(i+1,id+1);
+			}
+			else
+			{
+				break;
+			}
+		}
+
+		return *this;
+	}
+
+	/*! \brief Set the dimension
+	 *
+	 * Set the dimension
+	 *
+	 * \param d is the dimension
+	 * \param sz set the counter to sz
+	 *
+	 */
+	void set(int d, size_t sz)
+	{
+		// set the counter dim to sz
+
+		gk.set_d(d,sz);
+	}
+
+	/*! \brief Check if there is the next element
+	 *
+	 * Check if there is the next element
+	 *
+	 * \return true if there is the next, false otherwise
+	 *
+	 */
+
+	bool isNext()
+	{
+		if (gk.get(dim-1) < (long int)grid_base.size(dim-1))
+		{
+			//! we did not reach the end of the grid
+
+			return true;
+		}
+
+		//! we reach the end of the grid
+		return false;
+	}
+
+	/*! \brief Get the actual key
+	 *
+	 * Get the actual key
+	 *
+	 * \return the actual key
+	 *
+	 */
+	const grid_key_dx<dim> & get()
+	{
+		return gk;
+	}
+
+	/*! \brief Reinitialize the grid_key_dx_iterator
+	 *
+	 * \param key form
+	 *
+	 */
+	void reinitialize(const grid_key_dx_iterator<dim> & key)
+	{
+		grid_base = key.grid_base;
+		reset();
+	}
+
+	/*! \brief Reset the iterator (it restart from the beginning)
+	 *
+	 */
+	void reset()
+	{
+		//! Initialize to 0 the index
+
+		for (size_t i = 0 ; i < dim ; i++)
+		{gk.set_d(i,0);}
+	}
+};
+
+
+#endif /* OPENFPM_DATA_SRC_GRID_ITERATORS_GRID_KEY_DX_ITERATOR_HPP_ */
diff --git a/src/Grid/iterators/grid_key_dx_iterator_sp.hpp b/src/Grid/iterators/grid_key_dx_iterator_sp.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..6f75cac7376f1b5b6ece0672a524477760dc6f88
--- /dev/null
+++ b/src/Grid/iterators/grid_key_dx_iterator_sp.hpp
@@ -0,0 +1,116 @@
+/*
+ * grid_key_dx_iterator_sp.hpp
+ *
+ *  Created on: Dec 15, 2015
+ *      Author: i-bird
+ */
+
+#ifndef OPENFPM_DATA_SRC_GRID_ITERATORS_GRID_KEY_DX_ITERATOR_SP_HPP_
+#define OPENFPM_DATA_SRC_GRID_ITERATORS_GRID_KEY_DX_ITERATOR_SP_HPP_
+
+
+
+/**
+ *
+ * Grid key class iterator, iterate through a starting linearized grid element
+ * to a stop linearized grid element in particular if linearize is the function
+ *  that linearize all the grid_key_dx, it create an iterator that pass through
+ *  Linearize^(-1)(start) Linearize^(-1)(start+1) ....... Linearize^(-1)(stop)
+ *
+ * \param dim dimensionality of the grid
+ *
+ * Usage: In general you never create object directly, but you get it from a grid_cpu or grid_gpu with
+ *        getIteratorLinStartStop()
+ *
+ */
+
+template<unsigned int dim>
+class grid_key_dx_iterator_sp : public grid_key_dx_iterator<dim>
+{
+	//! stop point
+	grid_key_dx<dim> gk_stop;
+
+public:
+
+	/*! \brief Constructor require a grid grid<dim,T>
+	 *
+	 * It construct an iterator from one index to another, in particular
+	 * if linearize is the function that linearize all the grid_key, it
+	 * create an iterator that pass through Linearize^(-1)(start)
+	 * Linearize^(-1)(start+1) ....... Linearize^(-1)(stop)
+	 *
+	 * For example for start (1,1) and stop (3,3) the point indicated with # are
+	 * explored by the iterator
+	 *
+	 * \verbatim
+	 *
+                      +-----+-----+-----+-----+-----+-----+ (6,5)
+                      |     |     |     |     |     |     |
+                      |     |     |     |     |     |     |
+                      |     |     |     |     |     |     |
+                      +-----+-----+-----+-----+-----+-----+
+                      |     |     |     |     |     |     |
+                      |     |     |     |     |     |     |
+                      |     |     |     |     |     |     |
+                      #-----#-----#-----#-----+-----+-----+
+                      |     |     |     |     |     |     |
+                      |     |     |     |     |     |     |
+                      |     |     |     |     |     |     |
+                      #-----#-----#-----#-----#-----#-----#
+                      |     |     |     |     |     |     |
+                      |     |     |     |     |     |     |
+                      +-----#-----#-----#-----#-----#-----#
+                      |     |     |     |     |     |     |
+                      |     |     |     |     |     |     |
+                      |     |     |     |     |     |     |
+                      +-----+-----+-----+-----+-----+-----+
+                    (0,0)
+	 *
+	 *
+	 * \endverbatim
+	 *
+	 *
+	 * \tparam T type of object that the grid store
+	 *
+	 * \param g Grid on which iterate
+	 * \param from starting point
+	 * \param to end point
+	 *
+	 */
+	template<typename T> grid_key_dx_iterator_sp(grid_sm<dim,T> & g, mem_id from, mem_id to)
+	:grid_key_dx_iterator<dim>(g)
+	 {
+		//! Convert to a grid_key
+		this->gk = g.InvLinId(from);
+
+		//! Convert to a grid_key
+		gk_stop = g.InvLinId(to);
+	 }
+
+	/*! \brief Check if there is the next element
+	 *
+	 * Check if there is the next element
+	 *
+	 * \return true if there is the next, false otherwise
+	 *
+	 */
+
+	bool isNext()
+	{
+		//! for all dimensions
+		for (int i = dim-1 ; i >= 0 ; i-- )
+		{
+			//! check if we still have points
+			if (this->gk.get(i) < gk_stop.get(i))
+				return true;
+			else if (this->gk.get(i) > gk_stop.get(i))
+				return false;
+		}
+
+		//! (Final point) we we still have one point
+		return true;
+	}
+};
+
+
+#endif /* OPENFPM_DATA_SRC_GRID_ITERATORS_GRID_KEY_DX_ITERATOR_SP_HPP_ */
diff --git a/src/Grid/iterators/grid_key_dx_iterator_sub.hpp b/src/Grid/iterators/grid_key_dx_iterator_sub.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..b678bfa8c9725ff9d8fc322c4e576a332680c823
--- /dev/null
+++ b/src/Grid/iterators/grid_key_dx_iterator_sub.hpp
@@ -0,0 +1,405 @@
+/*
+ * grid_key_dx_iterator_sub.hpp
+ *
+ *  Created on: Dec 15, 2015
+ *      Author: i-bird
+ */
+
+#ifndef OPENFPM_DATA_SRC_GRID_ITERATORS_GRID_KEY_DX_ITERATOR_SUB_HPP_
+#define OPENFPM_DATA_SRC_GRID_ITERATORS_GRID_KEY_DX_ITERATOR_SUB_HPP_
+
+/* \brief grid_key_dx_iterator_sub can adjust the domain if the border go out-of-side
+ *        in this case a warning is produced
+ *
+ * \tparam dim dimensionality
+ *
+ */
+template <unsigned int dim>
+class print_warning_on_adjustment
+{
+public:
+	inline static void pw1(size_t i, const grid_key_dx<dim> & gk_start) {std::cerr << "Warning: " << __FILE__ << ":" << __LINE__ << " start with index smaller than zero x[" << i << "]=" << gk_start.get(i) << "\n";}
+	inline static void pw2(size_t i, const grid_key_dx<dim> & gk_start) {std::cerr << "Warning: " << __FILE__ << ":" << __LINE__ << " Cropping start point x[" << i << "]=" << gk_start.get(i) << " to x[" << i << "]=0 \n"  ;}
+	inline static void pw3() {std::cerr << "Warning: " << __FILE__ << ":" << __LINE__ << " stop with smaller index than start \n";}
+	inline static void pw4(size_t i, const grid_key_dx<dim> & gk_stop, const grid_sm<dim,void> & grid_base) {std::cerr << "Warning: " << __FILE__ << ":" << __LINE__ << " stop index bigger than cell domain x[" << i << "]=" << gk_stop.get(i) << " >= " << grid_base.size(i) << "\n";}
+	inline static void pw5() {std::cerr << "Warning grid_key_dx_iterator_sub: " << __FILE__ << ":" << __LINE__ << " the starting point is not smaller or equal than than the stop point in all the coordinates" << "\n";}
+};
+
+/* \brief grid_key_dx_iterator_sub can adjust the domain if the border go out-of-side
+ *        in this case a warning is produced
+ *
+ * \tparam dim dimensionality
+ * \tparam gb type of grid
+ *
+ */
+template <unsigned int dim>
+class do_not_print_warning_on_adjustment
+{
+public:
+	inline static void pw1(size_t i, const grid_key_dx<dim> & gk_start) {}
+	inline static void pw2(size_t i, const grid_key_dx<dim> & gk_start) {}
+	inline static void pw3() {}
+	inline static void pw4(size_t i, const grid_key_dx<dim> & gk_stop, const grid_sm<dim,void> & grid_base) {}
+	inline static void pw5() {}
+};
+
+
+/**
+ *
+ * Grid key class iterator, iterate through a sub-grid defined by an hyper-cube
+ *
+ * \param dim dimensionality of the grid
+ *
+ * \note if you have a grid you can get this object from getIteratorSub()
+ *
+ * ### Sub grid iterator declaration and usage
+ * \snippet grid_unit_tests.hpp Sub-grid iterator test usage
+ *
+ */
+
+template<unsigned int dim, typename warn>
+class grid_key_dx_iterator_sub : public grid_key_dx_iterator<dim>
+{
+#ifdef DEBUG
+	bool initialized = false;
+#endif
+
+	//! grid base where we are iterating
+	grid_sm<dim,void> grid_base;
+
+	//! start point
+	grid_key_dx<dim> gk_start;
+
+	//! stop point
+	grid_key_dx<dim> gk_stop;
+
+	/*! \brief Initialize gk
+	 *
+	 */
+	void Initialize()
+	{
+		// Check that start and stop are inside the domain otherwise crop them
+
+		for (size_t i = 0 ; i < dim ; i++)
+		{
+			// if start smaller than 0
+			if (gk_start.get(i) < 0)
+			{
+#ifdef DEBUG
+				warn::pw1(i,gk_start);
+#endif
+				if (gk_start.get(i) < gk_stop.get(i))
+				{
+#ifdef DEBUG
+					warn::pw2(i,gk_start);
+#endif
+					gk_start.set_d(i,0);
+				}
+				else
+				{
+#ifdef DEBUG
+					warn::pw3();
+#endif
+					// No points are available
+					gk_start.set_d(dim-1,gk_stop.get(dim-1)+1);
+					break;
+				}
+			}
+
+			// if stop bigger than the domain
+			if (gk_stop.get(i) >= (long int)grid_base.size(i))
+			{
+#ifdef DEBUG
+				warn::pw4(i,gk_stop,grid_base);
+#endif
+
+				if (gk_start.get(i) < (long int)grid_base.size(i))
+					gk_stop.set_d(i,grid_base.size(i)-1);
+				else
+				{
+					// No point are available
+					gk_start.set_d(dim-1,gk_stop.get(dim-1)+1);
+					break;
+				}
+			}
+		}
+
+		//! Initialize gk
+		for (unsigned int i = 0 ; i < dim ; i++)
+		{
+			this->gk.set_d(i,gk_start.get(i));
+		}
+
+#ifdef DEBUG
+		initialized = true;
+#endif
+	}
+
+public:
+
+	/*! \brief Default constructor
+	 *
+	 * \warning extremely unsafe
+	 * If you use this constructor before use the iterator you should call reinitialize first
+	 *
+	 */
+	grid_key_dx_iterator_sub()
+{}
+
+	/*! \brief Constructor from another grid_key_dx_iterator_sub
+	 *
+	 * It construct an iterator over an hyper-cube defined by start and stop,
+	 * \warning if start and stop are outside the domain defined by g the intersection
+	 * will be considered
+	 *
+	 * \param g_s_it grid_key_dx_iterator_sub
+	 *
+	 */
+	grid_key_dx_iterator_sub(const grid_key_dx_iterator_sub<dim> & g_s_it)
+	: grid_key_dx_iterator_sub(g_s_it.grid_base,g_s_it.gk_start,g_s_it.gk_stop)
+	{}
+
+
+	/*! \brief Constructor require a grid grid<dim,T>
+	 *
+	 * It construct an iterator over an hyper-cube defined by start and stop,
+	 * \warning if start and stop are outside the domain defined by g the intersection
+	 * will be considered
+	 *
+	 * \tparam T type of object that the grid store
+	 *
+	 * \param g Grid on which iterate
+	 * \param start starting point
+	 * \param stop end point
+	 *
+	 */
+	template<typename T> grid_key_dx_iterator_sub(const grid_sm<dim,T> & g, const grid_key_dx<dim> & start, const grid_key_dx<dim> & stop)
+	: grid_key_dx_iterator<dim>(g),grid_base(g),gk_start(start), gk_stop(stop)
+	{
+#ifdef DEBUG
+		//! If we are on debug check that the stop grid_key id bigger than the start
+		//! grid_key
+
+		for (unsigned int i = 0 ; i < dim ; i++)
+		{
+			if (gk_start.get(i) > gk_stop.get(i))
+			{
+				warn::pw5();
+			}
+		}
+#endif
+
+		Initialize();
+	}
+
+
+	/*! \brief Constructor require a grid grid<dim,T>
+	 *
+	 * It construct an iterator over an hyper-cube defined by start and stop,
+	 * \warning if start and stop are outside the domain defined by g the intersection
+	 * will be considered
+	 *
+	 * \tparam T type of object that the grid store
+	 *
+	 * \param g info of the grid where we are iterating
+	 * \param m Margin of the domain
+	 *
+	 */
+	template<typename T> grid_key_dx_iterator_sub(const grid_sm<dim,T> & g, const size_t m)
+			: grid_key_dx_iterator<dim>(g),grid_base(g)
+			  {
+		// Initialize the start and stop point
+		for (unsigned int i = 0 ; i < dim ; i++)
+		{
+			gk_start.set(i,m);
+			gk_stop.set(i,g.size(i)-m);
+		}
+
+		//
+		Initialize();
+			  }
+
+	/*! \brief Constructor require a grid grid<dim,T>
+	 *
+	 * It construct an iterator over an hyper-cube defined by start and stop,
+	 *
+	 * \tparam T type of object that the grid store
+	 *
+	 * \param g Grid on which iterate
+	 * \param start starting point
+	 * \param stop end point
+	 *
+	 */
+	template<typename T> grid_key_dx_iterator_sub(const grid_sm<dim,T> & g, const size_t (& start)[dim], const size_t (& stop)[dim])
+			: grid_key_dx_iterator<dim>(g),grid_base(g),gk_start(start), gk_stop(stop)
+			  {
+#ifdef DEBUG
+		//! If we are on debug check that the stop grid_key id bigger than the start
+		//! grid_key
+
+		for (unsigned int i = 0 ; i < dim ; i++)
+		{
+			if (start[i] > stop[i])
+			{
+				warn::pw5();
+			}
+		}
+#endif
+
+		Initialize();
+			  }
+
+	/*! \brief Get the next element
+	 *
+	 * Get the next element
+	 *
+	 * \return the next grid_key
+	 *
+	 */
+
+	grid_key_dx_iterator<dim> & operator++()
+			{
+#ifdef DEBUG
+		if (initialized == false)
+		{std::cerr << "Error: " << __FILE__ << __LINE__ << " using unitialized iterator" << "\n";}
+#endif
+
+		//! increment the first index
+
+		size_t id = this->gk.get(0);
+		this->gk.set_d(0,id+1);
+
+		//! check the overflow of all the index with exception of the last dimensionality
+
+		size_t i = 0;
+		for ( ; i < dim-1 ; i++)
+		{
+			size_t id = this->gk.get(i);
+			if ((long int)id > gk_stop.get(i))
+			{
+				// ! overflow, increment the next index
+
+				this->gk.set_d(i,gk_start.get(i));
+				id = this->gk.get(i+1);
+				this->gk.set_d(i+1,id+1);
+			}
+			else
+			{
+				break;
+			}
+		}
+
+		return *this;
+			}
+
+	/*! \brief Check if there is the next element
+	 *
+	 * Check if there is the next element
+	 *
+	 * \return true if there is the next, false otherwise
+	 *
+	 */
+
+	inline bool isNext()
+	{
+#ifdef DEBUG
+		if (initialized == false)
+		{std::cerr << "Error: " << __FILE__ << __LINE__ << " using unitialized iterator" << "\n";}
+#endif
+
+		if (this->gk.get(dim-1) <= gk_stop.get(dim-1))
+		{
+			//! we did not reach the end of the grid
+
+			return true;
+		}
+
+		//! we reach the end of the grid
+		return false;
+	}
+
+	/*! \brief Return the actual grid key iterator
+	 *
+	 */
+	inline grid_key_dx<dim> get()
+	{
+#ifdef DEBUG
+		if (initialized == false)
+		{std::cerr << "Error: " << __FILE__ << __LINE__ << " using unitialized iterator" << "\n";}
+#endif
+
+		return grid_key_dx_iterator<dim>::get();
+	}
+
+	/*! \brief Reinitialize the iterator
+	 *
+	 * it re-initialize the iterator with the passed grid_key_dx_iterator_sub
+	 * the actual position of the grid_key_dx_iterator_sub is ignored
+	 *
+	 * \param g_s_it grid_key_dx_iterator_sub
+	 *
+	 */
+
+	inline void reinitialize(const grid_key_dx_iterator_sub<dim,warn> & g_s_it)
+	{
+		// Reinitialize the iterator
+
+		grid_key_dx_iterator<dim>::reinitialize(g_s_it);
+		grid_base = g_s_it.grid_base;
+		gk_start = g_s_it.gk_start;
+		gk_stop = g_s_it.gk_stop;
+
+
+#ifdef DEBUG
+		//! If we are on debug check that the stop grid_key id bigger than the start
+		//! grid_key
+
+		for (unsigned int i = 0 ; i < dim ; i++)
+		{
+			if (gk_start.get(i) > gk_stop.get(i))
+			{
+				warn::pw5();
+			}
+		}
+
+		initialized = true;
+#endif
+
+		Initialize();
+	}
+
+	/*! \brief Get the volume spanned by this sub-grid iterator
+	 *
+	 * \return the volume
+	 *
+	 */
+	inline size_t getVolume()
+	{
+		return Box<dim,long int>::getVolumeKey(gk_start.k, gk_stop.k);
+	}
+
+	/*! \brief Reset the iterator (it restart from the beginning)
+	 *
+	 */
+	inline void reset()
+	{
+		//! Initialize to 0 the index
+
+		for (size_t i = 0 ; i < dim ; i++)
+		{this->gk.set_d(i,gk_start.get(i));}
+	}
+
+	/*! \brief Return the grid information related to this grid
+	 *
+	 *
+	 *
+	 */
+	inline const grid_sm<dim,void> & getGridInfo()
+	{
+		return grid_base;
+	}
+};
+
+
+
+#endif /* OPENFPM_DATA_SRC_GRID_ITERATORS_GRID_KEY_DX_ITERATOR_SUB_HPP_ */
diff --git a/src/Grid/iterators/grid_key_dx_iterator_sub_bc.hpp b/src/Grid/iterators/grid_key_dx_iterator_sub_bc.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0a86e533103fdd13f65e942f500a9f4284600936
--- /dev/null
+++ b/src/Grid/iterators/grid_key_dx_iterator_sub_bc.hpp
@@ -0,0 +1,184 @@
+/*
+ * grid_key_dx_iterator_sub_p.hpp
+ *
+ *  Created on: Dec 15, 2015
+ *      Author: i-bird
+ */
+
+#ifndef OPENFPM_DATA_SRC_GRID_ITERATORS_GRID_KEY_DX_ITERATOR_SUB_BC_HPP_
+#define OPENFPM_DATA_SRC_GRID_ITERATORS_GRID_KEY_DX_ITERATOR_SUB_BC_HPP_
+
+
+
+/*! \brief The same as grid_key_dx_iterator_sub_p but with periodic boundary
+ *
+ * In this case the boundaries are periodic
+ *
+ */
+template<unsigned int dim, typename warn=print_warning_on_adjustment<dim>>
+class grid_key_dx_iterator_sub_bc : public grid_key_dx_iterator_sub<dim,warn>
+{
+	// Boundary condition
+	size_t bc[dim];
+
+	// actual iterator box
+	size_t act;
+
+	//! Here we have all the boxes that this iterator produce
+	std::vector<Box<dim,size_t>> boxes;
+
+	/*! \brief Check if the position is valid with the actual boundary conditions
+	 *
+	 * \param key key to check
+	 *
+	 */
+	bool inline check_invalid(const grid_key_dx<dim> & key, const size_t (& bc)[dim])
+	{
+		for (size_t i = 0 ; i < dim ; i++)
+		{
+			if (bc[i] == NON_PERIODIC && key.get(i) != 1)
+			{
+				return true;
+			}
+		}
+
+		return false;
+	}
+
+public:
+
+	/*! \brief Constructor
+	 *
+	 * \param g Grid information
+	 * \param
+	 *
+	 */
+	template<typename T> grid_key_dx_iterator_sub_bc(const grid_sm<dim,T> & g, const grid_key_dx<dim> & start , const grid_key_dx<dim> & stop, const size_t (& bc)[dim])
+	:act(0)
+	{
+		// copy the boundary conditions
+
+		for (size_t i = 0 ; i < dim ; i++)
+			this->bc[i] = bc[i];
+
+		// compile-time array {0,0,0,....}  {2,2,2,...} {1,1,1,...}
+
+		typedef typename generate_array<size_t,dim, Fill_zero>::result NNzero;
+		typedef typename generate_array<size_t,dim, Fill_two>::result NNtwo;
+		typedef typename generate_array<size_t,dim, Fill_three>::result NNthree;
+
+		// Generate the sub-grid iterator
+
+		grid_sm<dim,void> nn(NNthree::data);
+		grid_key_dx_iterator_sub<dim> it(nn,NNzero::data,NNtwo::data);
+
+		// Box base
+		Box<dim,long int> base_b(start,stop);
+
+		// intersect with all the boxes
+		while (it.isNext())
+		{
+			auto key = it.get();
+
+			if (check_invalid(key,bc) == true)
+			{
+				++it;
+				continue;
+			}
+
+			bool intersect;
+
+			// intersection box
+			Box<dim,long int> b_int;
+			Box<dim,long int> b_out;
+
+			for (size_t i = 0 ; i < dim ; i++)
+			{
+				b_int.setLow(i,(key.get(i)-1)*g.getSize()[i]);
+				b_int.setHigh(i,key.get(i)*g.getSize()[i]-1);
+			}
+
+			intersect = base_b.Intersect(b_int,b_out);
+
+			// Bring to 0 and size[i]
+			for (size_t i = 0 ; i < dim ; i++)
+			{
+				if (bc[i] == PERIODIC)
+				{
+					b_out.setLow(i,openfpm::math::positive_modulo(b_out.getLow(i),g.size(i)));
+					b_out.setHigh(i,openfpm::math::positive_modulo(b_out.getHigh(i),g.size(i)));
+				}
+			}
+
+			// if intersect add in the box list
+			if (intersect == true)
+				boxes.push_back(b_out);
+
+			++it;
+		}
+
+		// initialize the first iterator
+		if (boxes.size() > 0)
+			grid_key_dx_iterator_sub<dim,warn>::reinitialize(grid_key_dx_iterator_sub<dim,warn>(g,boxes[0].getKP1(),boxes[0].getKP2()));
+	}
+
+
+	/*! \brief Get the next element
+	 *
+	 * Get the next element
+	 *
+	 * \return the next grid_key
+	 *
+	 */
+
+	grid_key_dx_iterator_sub_bc<dim,warn> & operator++()
+	{
+		grid_key_dx_iterator_sub<dim,warn>::operator++();
+		if (grid_key_dx_iterator_sub<dim,warn>::isNext() == true)
+		{
+			return *this;
+		}
+		else
+		{
+			act++;
+			if (act < boxes.size())
+				grid_key_dx_iterator_sub<dim,warn>::reinitialize(grid_key_dx_iterator_sub<dim,warn>(this->getGridInfo(),boxes[act].getKP1(),boxes[act].getKP2()));
+		}
+
+		return *this;
+	}
+
+	/*! \brief Check if there is the next element
+	 *
+	 * Check if there is the next element
+	 *
+	 * \return true if there is the next, false otherwise
+	 *
+	 */
+
+	inline bool isNext()
+	{
+		return act < boxes.size();
+	}
+
+	/*! \brief Return the actual grid key iterator
+	 *
+	 */
+	inline grid_key_dx<dim> get()
+	{
+		return grid_key_dx_iterator_sub<dim,warn>::get();
+	}
+
+	/*! \brief Reset the iterator (it restart from the beginning)
+	 *
+	 */
+	inline void reset()
+	{
+		act = 0;
+		// initialize the first iterator
+		reinitialize(boxes.get(0).getKP1,boxes.get(0).getKP2);
+	}
+};
+
+
+#endif /* OPENFPM_DATA_SRC_GRID_ITERATORS_GRID_KEY_DX_ITERATOR_SUB_BC_HPP_ */
diff --git a/src/Grid/map_grid.hpp b/src/Grid/map_grid.hpp
index c16aff6fbc55c35b6186d72cbebc748963281585..aa51ea60a790607dfb61d24e3e5092894f3dc091 100755
--- a/src/Grid/map_grid.hpp
+++ b/src/Grid/map_grid.hpp
@@ -37,6 +37,10 @@
 #include "memory/HeapMemory.hpp"
 #include "grid_common.hpp"
 #include "util/se_util.hpp"
+#include "iterators/grid_key_dx_iterator.hpp"
+#include "iterators/grid_key_dx_iterator_sub.hpp"
+#include "iterators/grid_key_dx_iterator_sp.hpp"
+#include "iterators/grid_key_dx_iterator_sub_bc.hpp"
 
 #ifndef CUDA_GPU
 typedef HeapMemory CudaMemory;
@@ -93,7 +97,7 @@ private:
 	//! Is the memory initialized
 	bool is_mem_init = false;
 
-	//! This is an header that store all information related to the grid
+	//! This is a structure that store all information related to the grid and how indexes are linearized
 	grid_sm<dim,T> g1;
 
 	//! Memory layout specification + memory chunk pointer
@@ -953,6 +957,8 @@ public:
 	{
 #ifdef SE_CLASS2
 		return check_whoami(this,8);
+#else
+		return -1;
 #endif
 	}
 };
diff --git a/src/NN/CellList/CellList.hpp b/src/NN/CellList/CellList.hpp
index 1ba2cfc741cd6c16a60dc135d015f3d9aafdfe0e..b5dd32c15a700ffc0d497595495009dc7ec93ed4 100644
--- a/src/NN/CellList/CellList.hpp
+++ b/src/NN/CellList/CellList.hpp
@@ -12,7 +12,7 @@
 #include "CellDecomposer.hpp"
 
 //! Point at witch the cell do a reallocation (it should the the maximum for all the implementations)
-#define CELL_REALLOC 16
+#define CELL_REALLOC 16ul
 
 /* NOTE all the implementations
  *
diff --git a/src/NN/CellList/CellListFast.hpp b/src/NN/CellList/CellListFast.hpp
index e2dfeb67a9c7bbb53b1cab087cea2e1a66280770..66c37366dff99c5ec32288b337bfd130e6ff0a0a 100644
--- a/src/NN/CellList/CellListFast.hpp
+++ b/src/NN/CellList/CellListFast.hpp
@@ -14,25 +14,7 @@
 #include "CellNNIterator.hpp"
 #include "Space/Shape/HyperCube.hpp"
 
-// Compile time array functor needed to generate array at compile-time of type
-// {0,0,0,0,0,.....}
-// {3,3,3,3,3,3,.....}
-
- template<size_t index, size_t N> struct Fill_three {
-    enum { value = 3 };
- };
-
- template<size_t index, size_t N> struct Fill_zero {
-    enum { value = 0 };
- };
-
- template<size_t index, size_t N> struct Fill_two {
-    enum { value = 2 };
- };
-
- template<size_t index, size_t N> struct Fill_one {
-    enum { value = 1 };
- };
+#include "util/common.hpp"
 
 /*! \brief Class for FAST cell list implementation
  *
diff --git a/src/NN/CellList/CellList_test.hpp b/src/NN/CellList/CellList_test.hpp
index 7efa2c2470a4d8da4a3df1df3807d8f1e9df13b5..1734017b0724ed262dc6acfad05546cfd5359a96 100644
--- a/src/NN/CellList/CellList_test.hpp
+++ b/src/NN/CellList/CellList_test.hpp
@@ -37,14 +37,14 @@ template<unsigned int dim, typename T, typename CellS> void Test_cell_s()
 	grid_sm<dim,void> g_info(div);
 
 	// Test force reallocation in case of Cell list fast
-	for (int i = 0 ; i < CELL_REALLOC * 3 ; i++)
+	for (size_t i = 0 ; i < CELL_REALLOC * 3 ; i++)
 	{
 		cl2.add(org,i);
 	}
 
 	// Check the elements
-	BOOST_REQUIRE_EQUAL(cl2.getNelements(cl2.getCell(org)),CELL_REALLOC * 3);
-	for (int i = 0 ; i < CELL_REALLOC * 3 ; i++)
+	BOOST_REQUIRE_EQUAL(cl2.getNelements(cl2.getCell(org)),CELL_REALLOC * 3ul);
+	for (size_t i = 0 ; i < CELL_REALLOC * 3 ; i++)
 	{
 		BOOST_REQUIRE_EQUAL(cl2.get(cl2.getCell(org),i),i);
 	}
@@ -111,7 +111,7 @@ template<unsigned int dim, typename T, typename CellS> void Test_cell_s()
 		size_t cell = cl1.getCell(key);
 		size_t n_ele = cl1.getNelements(cell);
 
-		BOOST_REQUIRE_EQUAL(n_ele,2);
+		BOOST_REQUIRE_EQUAL(n_ele,2ul);
 		BOOST_REQUIRE_EQUAL(cl1.get(cell,1) - cl1.get(cell,0),1);
 
 		++g_it;
@@ -151,7 +151,7 @@ template<unsigned int dim, typename T, typename CellS> void Test_cell_s()
 		auto cell = cl1.getCell(key);
 		size_t n_ele = cl1.getNelements(cell);
 
-		BOOST_REQUIRE_EQUAL(n_ele,1);
+		BOOST_REQUIRE_EQUAL(n_ele,1ul);
 		++g_it;
 	}
 
@@ -184,7 +184,7 @@ template<unsigned int dim, typename T, typename CellS> void Test_cell_s()
 			++NN;
 		}
 
-		BOOST_REQUIRE_EQUAL(total,openfpm::math::pow(3,dim));
+		BOOST_REQUIRE_EQUAL(total,(size_t)openfpm::math::pow(3,dim));
 		++g_it_s;
 	}
 
@@ -214,14 +214,14 @@ BOOST_AUTO_TEST_CASE( CellDecomposer_use )
 	{
 	CellDecomposer_sm<3,double> cd(box,div,0);
 	size_t cell = cd.getCell(p);
-	BOOST_REQUIRE_EQUAL(cell,8*16*16 + 8*16 + 8);
+	BOOST_REQUIRE_EQUAL(cell,(size_t)(8*16*16 + 8*16 + 8));
 	auto key = cd.getCellGrid(p);
 	BOOST_REQUIRE_EQUAL(key.get(0),8);
 	BOOST_REQUIRE_EQUAL(key.get(1),8);
 	BOOST_REQUIRE_EQUAL(key.get(2),8);
 
 	cell = cd.getCell(pp);
-	BOOST_REQUIRE_EQUAL(cell,8*16*16 + 8*16 + 8);
+	BOOST_REQUIRE_EQUAL(cell,(size_t)(8*16*16 + 8*16 + 8));
 	key = cd.getCellGrid(pp);
 	BOOST_REQUIRE_EQUAL(key.get(0),8);
 	BOOST_REQUIRE_EQUAL(key.get(1),8);
@@ -234,14 +234,14 @@ BOOST_AUTO_TEST_CASE( CellDecomposer_use )
 	{
 	CellDecomposer_sm<3,double> cd(box,div,1);
 	size_t cell = cd.getCell(p);
-	BOOST_REQUIRE_EQUAL(cell,9*18*18 + 9*18 + 9);
+	BOOST_REQUIRE_EQUAL(cell,(size_t)(9*18*18 + 9*18 + 9));
 	auto key = cd.getCellGrid(p);
 	BOOST_REQUIRE_EQUAL(key.get(0),9);
 	BOOST_REQUIRE_EQUAL(key.get(1),9);
 	BOOST_REQUIRE_EQUAL(key.get(2),9);
 
 	cell = cd.getCell(pp);
-	BOOST_REQUIRE_EQUAL(cell,9*18*18 + 9*18 + 9);
+	BOOST_REQUIRE_EQUAL(cell,(size_t)(9*18*18 + 9*18 + 9));
 	key = cd.getCellGrid(pp);
 	BOOST_REQUIRE_EQUAL(key.get(0),9);
 	BOOST_REQUIRE_EQUAL(key.get(1),9);
@@ -257,14 +257,14 @@ BOOST_AUTO_TEST_CASE( CellDecomposer_use )
 
 	CellDecomposer_sm< 3,double,shift<3,double> > cd(box,div,sht,1);
 	size_t cell = cd.getCell(p + sht);
-	BOOST_REQUIRE_EQUAL(cell,9*18*18 + 9*18 + 9);
+	BOOST_REQUIRE_EQUAL(cell,(size_t)(9*18*18 + 9*18 + 9));
 	auto key = cd.getCellGrid(p + sht);
 	BOOST_REQUIRE_EQUAL(key.get(0),9);
 	BOOST_REQUIRE_EQUAL(key.get(1),9);
 	BOOST_REQUIRE_EQUAL(key.get(2),9);
 
 	cell = cd.getCell(psht);
-	BOOST_REQUIRE_EQUAL(cell,9*18*18 + 9*18 + 9);
+	BOOST_REQUIRE_EQUAL(cell,(size_t)(9*18*18 + 9*18 + 9));
 	key = cd.getCellGrid(psht);
 	BOOST_REQUIRE_EQUAL(key.get(0),9);
 	BOOST_REQUIRE_EQUAL(key.get(1),9);
@@ -323,7 +323,7 @@ BOOST_AUTO_TEST_CASE( CellDecomposer_get_grid_points )
 	// Get the volume of the box
 	size_t vol = gp.getVolumeKey();
 
-	BOOST_REQUIRE_EQUAL(vol,12);
+	BOOST_REQUIRE_EQUAL(vol,12ul);
 	}
 
 	//////////////////////////
@@ -360,7 +360,7 @@ BOOST_AUTO_TEST_CASE( CellDecomposer_get_grid_points )
 	// Get the volume of the box
 	size_t vol = gp.getVolumeKey();
 
-	BOOST_REQUIRE_EQUAL(vol,4);
+	BOOST_REQUIRE_EQUAL(vol,4ul);
 	}
 }
 
diff --git a/src/Space/Shape/Box.hpp b/src/Space/Shape/Box.hpp
index 26e9693f11babf7d0eb5284646271065a605da25..eca8845fe9282aa0301c58d5e3b4b38c90d36699 100644
--- a/src/Space/Shape/Box.hpp
+++ b/src/Space/Shape/Box.hpp
@@ -235,6 +235,18 @@ public:
 	Box()
 	{}
 
+	/*! \brief Constructor from two points
+	 *
+	 * \param p1 Low point, initialize as a list example {0.0,0.0,0.0}
+	 * \param p2 High point, initialized as a list example {1.0,1.0,1.0}
+	 *
+	 */
+	Box(const Point<dim,T> & p1, const Point<dim,T> & p2)
+	{
+		setP1(p1);
+		setP2(p2);
+	}
+
 	/*! \brief Constructor from initializer list
 	 *
 	 * \param p1 Low point, initialize as a list example {0.0,0.0,0.0}
@@ -313,6 +325,21 @@ public:
 		}
 	}
 
+	/*! \brief constructor from 2 grid_key_dx
+	 *
+	 * \param key1 start point
+	 * \param key2 stop point
+	 *
+	 */
+	inline Box(const grid_key_dx<dim> & key1, const grid_key_dx<dim> & key2)
+	{
+		for (size_t i = 0 ; i < dim ; i++)
+		{
+			setLow(i,key1.get(i));
+			setHigh(i,key2.get(i));
+		}
+	}
+
 	/*! \brief Box constructor from vector::fusion of higher dimension
 	 *
 	 * \param box_data fusion vector from which to construct the vector
@@ -407,12 +434,22 @@ public:
 	inline void set(std::initializer_list<T> p1, std::initializer_list<T> p2)
 	{
 		size_t i = 0;
-	    for(T x : p1)
-	    {setLow(i,x);i++;}
+		for(T x : p1)
+		{
+			setLow(i,x);
+			i++;
+			if (i > dim)
+				break;
+		}
 
-	    i = 0;
-	    for(T x : p2)
-	    {setHigh(i,x);i++;}
+		i = 0;
+		for(T x : p2)
+		{
+			setHigh(i,x);
+			i++;
+			if (i > dim)
+				break;
+		}
 	}
 
 	/*! \brief set the low interval of the box
@@ -598,6 +635,24 @@ public:
 		return *this;
 	}
 
+	/*! \brief Translate the box
+	 *
+	 * \p Point translation vector
+	 *
+	 * \return itself
+	 *
+	 */
+	inline Box<dim,T> & operator+=(const Point<dim,T> & p)
+	{
+		for (size_t i = 0 ; i < dim ; i++)
+		{
+			boost::fusion::at_c<p2>(data)[i] += p.get(i);
+			boost::fusion::at_c<p1>(data)[i] += p.get(i);
+		}
+
+		return *this;
+	}
+
 	/* \brief expand expand the box by a vector
 	 *
 	 * \param vector
diff --git a/src/Space/Shape/HyperCube.hpp b/src/Space/Shape/HyperCube.hpp
index 5e02bb0a551afc2e118d3bafd66049d2337c7b00..8b73e0549ec34916916f36816a3d812c74a08d09 100644
--- a/src/Space/Shape/HyperCube.hpp
+++ b/src/Space/Shape/HyperCube.hpp
@@ -96,6 +96,11 @@ public:
 
 	static std::vector<comb<dim>> getCombinations_R(size_t d)
 	{
+#ifdef SE_CLASS1
+		if (d > dim)
+			std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " " << d << " must be smaller than " << "\n";
+#endif
+
 		// Create an Iterator_g_const
 		// And a vector that store all the combination
 
diff --git a/src/Space/Shape/HyperCube_unit_test.hpp b/src/Space/Shape/HyperCube_unit_test.hpp
index 5cdb31075e49d3768b8f9c8beb814b371bc8a05b..7d55712a714febad856f7aec0d8ef423d922c1b3 100644
--- a/src/Space/Shape/HyperCube_unit_test.hpp
+++ b/src/Space/Shape/HyperCube_unit_test.hpp
@@ -167,8 +167,8 @@ BOOST_AUTO_TEST_CASE( Hyper_cube_use)
 	//! [Get vertex and edge on a line]
 
 	// Check
-	BOOST_REQUIRE_EQUAL(ele[0],2);
-	BOOST_REQUIRE_EQUAL(ele[1],1);
+	BOOST_REQUIRE_EQUAL(ele[0],2ul);
+	BOOST_REQUIRE_EQUAL(ele[1],1ul);
 
 	// Fill the expected vector
 
@@ -199,9 +199,9 @@ BOOST_AUTO_TEST_CASE( Hyper_cube_use)
 	//! [Get vertex edge and surfaces of a square]
 
 	// Check
-	BOOST_REQUIRE_EQUAL(ele[0],4);
-	BOOST_REQUIRE_EQUAL(ele[1],4);
-	BOOST_REQUIRE_EQUAL(ele[2],1);
+	BOOST_REQUIRE_EQUAL(ele[0],4ul);
+	BOOST_REQUIRE_EQUAL(ele[1],4ul);
+	BOOST_REQUIRE_EQUAL(ele[2],1ul);
 
 	// Check combination
 
@@ -237,10 +237,10 @@ BOOST_AUTO_TEST_CASE( Hyper_cube_use)
 	//! [Get vertex edge surfaces and volumes of a cube]
 
 	// 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);
+	BOOST_REQUIRE_EQUAL(ele[0],8ul);
+	BOOST_REQUIRE_EQUAL(ele[1],12ul);
+	BOOST_REQUIRE_EQUAL(ele[2],6ul);
+	BOOST_REQUIRE_EQUAL(ele[3],1ul);
 
 	// Check combination
 
@@ -286,11 +286,11 @@ BOOST_AUTO_TEST_CASE( Hyper_cube_use)
 	check_lin_perm(v1_st);
 
 	// 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);
+	BOOST_REQUIRE_EQUAL(ele[0],16ul);
+	BOOST_REQUIRE_EQUAL(ele[1],32ul);
+	BOOST_REQUIRE_EQUAL(ele[2],24ul);
+	BOOST_REQUIRE_EQUAL(ele[3],8ul);
+	BOOST_REQUIRE_EQUAL(ele[4],1ul);
 
 	// Penteract
 	// Number of vertex
@@ -323,12 +323,12 @@ BOOST_AUTO_TEST_CASE( Hyper_cube_use)
 	check_lin_perm(v5_st);
 
 	// 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);
+	BOOST_REQUIRE_EQUAL(ele[0],32ul);
+	BOOST_REQUIRE_EQUAL(ele[1],80ul);
+	BOOST_REQUIRE_EQUAL(ele[2],80ul);
+	BOOST_REQUIRE_EQUAL(ele[3],40ul);
+	BOOST_REQUIRE_EQUAL(ele[4],10ul);
+	BOOST_REQUIRE_EQUAL(ele[5],1ul);
 
 
 	// Test SubHypercube 2D and 3D
@@ -343,7 +343,7 @@ BOOST_AUTO_TEST_CASE( Hyper_cube_use)
 		//! [Getting the vertex as a sub-hypercube]
 		std::vector<comb<2>> combs = SubHyperCube<2,0>::getCombinations_R(sc2_0[i],0);
 		//! [Getting the vertex as a sub-hypercube]
-		BOOST_REQUIRE_EQUAL(combs.size(),1);
+		BOOST_REQUIRE_EQUAL(combs.size(),1ul);
 		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
 		BOOST_REQUIRE_EQUAL(isValid(combs),true);
 		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc2_0[i]),true);
@@ -360,7 +360,7 @@ BOOST_AUTO_TEST_CASE( Hyper_cube_use)
 		//! [Getting the vertex of the line of the original square]
 		std::vector<comb<2>> combs = SubHyperCube<2,1>::getCombinations_R(sc2_1[i],0);
 		//! [Getting the vertex of the line of the original square]
-		BOOST_REQUIRE_EQUAL(combs.size(),2);
+		BOOST_REQUIRE_EQUAL(combs.size(),2ul);
 		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
 		BOOST_REQUIRE_EQUAL(isValid(combs),true);
 		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc2_1[i]),true);
@@ -369,7 +369,7 @@ BOOST_AUTO_TEST_CASE( Hyper_cube_use)
 		//! [Getting the edge(line) of the line of the original square]
 		combs = SubHyperCube<2,1>::getCombinations_R(sc2_1[i],1);
 		//! [Getting the edge(line) of the line of the original square]
-		BOOST_REQUIRE_EQUAL(combs.size(),1);
+		BOOST_REQUIRE_EQUAL(combs.size(),1ul);
 		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
 		BOOST_REQUIRE_EQUAL(isValid(combs),true);
 		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc2_1[i]),true);
@@ -385,7 +385,7 @@ BOOST_AUTO_TEST_CASE( Hyper_cube_use)
 		//! [Getting the vertices of the square of the original square]
 		std::vector<comb<2>> combs = SubHyperCube<2,2>::getCombinations_R(sc2_2[i],0);
 		//! [Getting the vertices of the square of the original square]
-		BOOST_REQUIRE_EQUAL(combs.size(),4);
+		BOOST_REQUIRE_EQUAL(combs.size(),4ul);
 		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
 		BOOST_REQUIRE_EQUAL(isValid(combs),true);
 		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc2_2[i]),true);
@@ -394,7 +394,7 @@ BOOST_AUTO_TEST_CASE( Hyper_cube_use)
 		//! [Getting the lines of the square of the original square]
 		combs = SubHyperCube<2,2>::getCombinations_R(sc2_2[i],1);
 		//! [Getting the lines of the square of the original square]
-		BOOST_REQUIRE_EQUAL(combs.size(),4);
+		BOOST_REQUIRE_EQUAL(combs.size(),4ul);
 		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
 		BOOST_REQUIRE_EQUAL(isValid(combs),true);
 		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc2_2[i]),true);
@@ -403,7 +403,7 @@ BOOST_AUTO_TEST_CASE( Hyper_cube_use)
 		//! [Getting the square of the square of the original square]
 		combs = SubHyperCube<2,2>::getCombinations_R(sc2_2[i],2);
 		//! [Getting the square of the square of the original square]
-		BOOST_REQUIRE_EQUAL(combs.size(),1);
+		BOOST_REQUIRE_EQUAL(combs.size(),1ul);
 		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
 		BOOST_REQUIRE_EQUAL(isValid(combs),true);
 		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc2_2[i]),true);
@@ -421,7 +421,7 @@ BOOST_AUTO_TEST_CASE( Hyper_cube_use)
 		//! [Getting the vertices of the vertices of the cube]
 		std::vector<comb<3>> combs = SubHyperCube<3,0>::getCombinations_R(sc3_0[i],0);
 		//! [Getting the vertices of the vertices of the cube]
-		BOOST_REQUIRE_EQUAL(combs.size(),1);
+		BOOST_REQUIRE_EQUAL(combs.size(),1ul);
 		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
 		BOOST_REQUIRE_EQUAL(isValid(combs),true);
 //		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc3_0[i]),true);
@@ -437,7 +437,7 @@ BOOST_AUTO_TEST_CASE( Hyper_cube_use)
 		//! [Getting the vertices of the edge of the cube]
 		std::vector<comb<3>> combs = SubHyperCube<3,1>::getCombinations_R(sc3_1[i],0);
 		//! [Getting the vertices of the vertices of the cube]
-		BOOST_REQUIRE_EQUAL(combs.size(),2);
+		BOOST_REQUIRE_EQUAL(combs.size(),2ul);
 		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
 		BOOST_REQUIRE_EQUAL(isValid(combs),true);
 		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc3_1[i]),true);
@@ -446,7 +446,7 @@ BOOST_AUTO_TEST_CASE( Hyper_cube_use)
 		//! [Getting the edges of the edges of the cube]
 		combs = SubHyperCube<3,1>::getCombinations_R(sc3_1[i],1);
 		//! [Getting the edges of the edges of the cube]
-		BOOST_REQUIRE_EQUAL(combs.size(),1);
+		BOOST_REQUIRE_EQUAL(combs.size(),1ul);
 		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
 		BOOST_REQUIRE_EQUAL(isValid(combs),true);
 		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc3_1[i]),true);
@@ -462,7 +462,7 @@ BOOST_AUTO_TEST_CASE( Hyper_cube_use)
 		//! [Getting the vertices of one surface of the cube]
 		std::vector<comb<3>> combs = SubHyperCube<3,2>::getCombinations_R(sc3_2[i],0);
 		//! [Getting the vertices of one surface of the cube]
-		BOOST_REQUIRE_EQUAL(combs.size(),4);
+		BOOST_REQUIRE_EQUAL(combs.size(),4ul);
 		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
 		BOOST_REQUIRE_EQUAL(isValid(combs),true);
 		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc3_2[i]),true);
@@ -471,7 +471,7 @@ BOOST_AUTO_TEST_CASE( Hyper_cube_use)
 		//! [Getting the edges of the surfaces of the cube]
 		combs = SubHyperCube<3,2>::getCombinations_R(sc3_2[i],1);
 		//! [Getting the edges of the surfaces of the cube]
-		BOOST_REQUIRE_EQUAL(combs.size(),4);
+		BOOST_REQUIRE_EQUAL(combs.size(),4ul);
 		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
 		BOOST_REQUIRE_EQUAL(isValid(combs),true);
 		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc3_2[i]),true);
@@ -480,7 +480,7 @@ BOOST_AUTO_TEST_CASE( Hyper_cube_use)
 		//! [Getting the surfaces of the surfaces of the cube]
 		combs = SubHyperCube<3,2>::getCombinations_R(sc3_2[i],2);
 		//! [Getting the surfaces of the surfaces of the cube]
-		BOOST_REQUIRE_EQUAL(combs.size(),1);
+		BOOST_REQUIRE_EQUAL(combs.size(),1ul);
 		BOOST_REQUIRE_EQUAL(isDinstict(combs),true);
 		BOOST_REQUIRE_EQUAL(isValid(combs),true);
 		BOOST_REQUIRE_EQUAL(isSubdecomposition(combs,sc3_2[i]),true);
diff --git a/src/Space/Shape/Point.hpp b/src/Space/Shape/Point.hpp
index 776de3b263817134002e68f8f118c125888c5975..462ab2329c609f3e82ce046a47628739650648a3 100644
--- a/src/Space/Shape/Point.hpp
+++ b/src/Space/Shape/Point.hpp
@@ -82,6 +82,23 @@ template<unsigned int dim ,typename T> class Point
 		return *this;
 	}
 
+	/*! \brief Multiply each components by a constant
+	 *
+	 * \param c constanr
+	 *
+	 */
+	inline Point<dim,T> operator*(T c)
+	{
+		Point<dim,T> result;
+
+		for (size_t i = 0 ; i < dim ; i++)
+		{
+			result.get(i) = get(i) * c;
+		}
+
+		return result;
+	}
+
 	/*! \brief Multiply each components
 	 *
 	 * \param p Point
@@ -104,7 +121,7 @@ template<unsigned int dim ,typename T> class Point
 	 * \param p Point
 	 *
 	 */
-	template<typename aT> inline Point<dim,T> & operator+=(const Point<dim,aT> & p)
+	inline Point<dim,T> & operator+=(const Point<dim,T> & p)
 	{
 		for (size_t i = 0 ; i < dim ; i++)
 		{
@@ -304,16 +321,16 @@ template<unsigned int dim ,typename T> class Point
 			get(i) = static_cast<S>(p.get(i));
 	}
 
-	/*! \brief Constructor from a grid_key_dx<dim>
+	/*! \brief Point constructor
 	 *
-	 * \param key from where to initialize
+	 * \param p encapc Point
 	 *
 	 */
-/*	inline Point(grid_key_dx<dim> key)
+	template <unsigned int d, typename M> inline Point(const encapc<d,Point<dim,T>,M> & p)
 	{
-	    for(size_t i = 0 ; i < dim ; i++)
-	    {get(i) = key.k[i];}
-	}*/
+		for (size_t i = 0 ; i < dim ; i++)
+			get(i) = p.template get<0>()[i];
+	}
 
 	/*! \brief Constructor from a list
 	 *
diff --git a/src/Vector/map_vector.hpp b/src/Vector/map_vector.hpp
index 20e2eca27f5fb5c1beb0c922c512f6fc9fba9cdb..772e9944c74506e120ada12990a45633af28d723 100644
--- a/src/Vector/map_vector.hpp
+++ b/src/Vector/map_vector.hpp
@@ -1066,6 +1066,8 @@ namespace openfpm
 		{
 #ifdef SE_CLASS2
 			return check_whoami(this,8);
+#else
+			return -1;
 #endif
 		}
 	};
diff --git a/src/Vector/map_vector_std.hpp b/src/Vector/map_vector_std.hpp
index 11debea4e851a2d7da3c04be095a7bffd99d6675..05a8f0a9352c287b11a32aec307ca14c7f7cd8ab 100644
--- a/src/Vector/map_vector_std.hpp
+++ b/src/Vector/map_vector_std.hpp
@@ -574,6 +574,8 @@ public:
 	{
 #ifdef SE_CLASS2
 		return check_whoami(this,8);
+#else
+		return -1;
 #endif
 	}
 };
diff --git a/src/Vector/vector_unit_tests.hpp b/src/Vector/vector_unit_tests.hpp
index cf3466d5db001f2ddcdffc10cb4ca1fa8c27f15c..863a7e80090d494cf94ab0a0c29f00f5f0887ba3 100644
--- a/src/Vector/vector_unit_tests.hpp
+++ b/src/Vector/vector_unit_tests.hpp
@@ -135,7 +135,7 @@ BOOST_AUTO_TEST_CASE( vector_std_utility )
 	// Check is zero
 	for (size_t i = 0 ;  i < 16 ; i++)
 	{
-		BOOST_REQUIRE_EQUAL(pb.get(i),0);
+		BOOST_REQUIRE_EQUAL(pb.get(i),0ul);
 	}
 
 }
@@ -260,7 +260,7 @@ BOOST_AUTO_TEST_CASE( object_test_creator )
 	BOOST_REQUIRE_EQUAL(tst , true);
 }
 
-#define V_REM_PUSH 1024
+#define V_REM_PUSH 1024ul
 
 BOOST_AUTO_TEST_CASE(vector_remove )
 {
@@ -291,7 +291,7 @@ BOOST_AUTO_TEST_CASE(vector_remove )
 
 	//! [Create push and multiple remove]
 
-	BOOST_REQUIRE_EQUAL(v1.size(),1020);
+	BOOST_REQUIRE_EQUAL(v1.size(),1020ul);
 	BOOST_REQUIRE_EQUAL(v1.template get<p::x>(0),4);
 
 	{
@@ -304,7 +304,7 @@ BOOST_AUTO_TEST_CASE(vector_remove )
 	v1.remove(rem);
 	}
 
-	BOOST_REQUIRE_EQUAL(v1.size(),1016);
+	BOOST_REQUIRE_EQUAL(v1.size(),1016ul);
 	BOOST_REQUIRE_EQUAL(v1.template get<p::x>(v1.size()-1),1019);
 
 	{
@@ -316,12 +316,12 @@ BOOST_AUTO_TEST_CASE(vector_remove )
 	v1.remove(rem);
 	}
 
-	BOOST_REQUIRE_EQUAL(v1.size(),508);
+	BOOST_REQUIRE_EQUAL(v1.size(),508ul);
 
 	// Check only odd
 	for (size_t i = 0 ; i < v1.size() ; i++)
 	{
-		BOOST_REQUIRE_EQUAL((size_t)v1.template get<p::x>(v1.size()-1) % 2, 1);
+		BOOST_REQUIRE_EQUAL((size_t)v1.template get<p::x>(v1.size()-1) % 2, 1ul);
 	}
 }
 
@@ -342,7 +342,7 @@ BOOST_AUTO_TEST_CASE(vector_clear )
 
 	v1.clear();
 
-	BOOST_REQUIRE_EQUAL(v1.size(),0);
+	BOOST_REQUIRE_EQUAL(v1.size(),0ul);
 
 	for (size_t i = 0 ; i < V_REM_PUSH ; i++)
 	{
diff --git a/src/main.cpp b/src/main.cpp
index 13061f7fde9a449a906e914ceead2e56c8dea1e0..3253d032f782be3a6f6838d615ecd0667ddfb6d2 100755
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -25,3 +25,5 @@
 #include "Space/Shape/HyperCube_unit_test.hpp"
 #include "Graph/graph_unit_tests.hpp"
 #include "Grid/grid_unit_tests.hpp"
+#include "Grid/grid_sm_unit_tests.hpp"
+#include "util/mathutil_unit_test.hpp"
diff --git a/src/util/common.hpp b/src/util/common.hpp
index b174ff52112a423bce91e362e6dff59171fbd12e..1b9927712dda789e1dfb03625c61877a12471deb 100644
--- a/src/util/common.hpp
+++ b/src/util/common.hpp
@@ -16,6 +16,26 @@ namespace std
 }
 
 
+// Compile time array functor needed to generate array at compile-time of type
+// {0,0,0,0,0,.....}
+// {3,3,3,3,3,3,.....}
+
+ template<size_t index, size_t N> struct Fill_three {
+    enum { value = 3 };
+ };
+
+ template<size_t index, size_t N> struct Fill_zero {
+    enum { value = 0 };
+ };
+
+ template<size_t index, size_t N> struct Fill_two {
+    enum { value = 2 };
+ };
+
+ template<size_t index, size_t N> struct Fill_one {
+    enum { value = 1 };
+ };
+
 //! Void structure
 template<typename> struct Void
 {
diff --git a/src/util/mathutil.hpp b/src/util/mathutil.hpp
index 5a03e2572968072aef4e5879bb74491ef73516c3..72953e500339e1ccce32a3a035186b326f56048e 100644
--- a/src/util/mathutil.hpp
+++ b/src/util/mathutil.hpp
@@ -11,11 +11,14 @@ namespace openfpm
 		 *
 		 * calculate the factorial of a number
 		 *
+		 * # Example
+		 * \snippet mathutil_unit_test.hpp factorial usage
+		 *
 		 * \param f number
 		 * \return the factorial
 		 *
 		 */
-		static size_t factorial(size_t f)
+		static inline size_t factorial(size_t f)
 		{
 			size_t fc = 1;
 
@@ -33,9 +36,14 @@ namespace openfpm
 		 *
 		 * n!/(k!(n-k)!)
 		 *
+		 * # Example
+		 * \snippet mathutil_unit_test.hpp Combination usage
+		 *
+		 * \param n
+		 * \param k
+		 *
 		 */
-
-		static size_t C(size_t n, size_t k)
+		static inline size_t C(size_t n, size_t k)
 		{
 			return factorial(n)/(factorial(k)*factorial(n-k));
 		}
@@ -43,12 +51,14 @@ namespace openfpm
 		/*! \brief Round to the nearest bigger power of 2 number
 		 *
 		 * \param n number
-		 *
 		 * \return nearest bigger power of 2 number
 		 *
+		 * # Example
+		 * \snippet mathutil_unit_test.hpp round to big pow
+		 *
+		 *
 		 */
-
-		static size_t round_big_2(size_t n)
+		static inline size_t round_big_2(size_t n)
 		{
 			n--;
 			n |= n >> 1;   // Divide by 2^k for consecutive doublings of k up to 32,
@@ -63,6 +73,9 @@ namespace openfpm
 
 
 		/*! \brief It calculate at compile-time and runtime the power with recursion
+		 *
+		 * # Example
+		 * \snippet mathutil_unit_test.hpp pow
 		 *
 		 * \tparam type of the pow expression
 		 *
@@ -79,15 +92,73 @@ namespace openfpm
 		}
 
 		/* \brief Return the positive modulo of a number
+		 *
+		 * # Example
+		 * \snippet mathutil_unit_test.hpp positive modulo
 		 *
 		 * \param i number
 		 * \param n modulo
 		 *
 		 */
-		inline long int positive_modulo(long int i, long int n)
+		static inline long int positive_modulo(long int i, long int n)
 		{
 		    return (i % n + n) % n;
 		}
+
+		/*! \brief Bound the position to be inside p2 and p1
+		 *
+		 * Given pos = 10.9, p2 = 1.0 and p1 = 0.1 and l = p2 - p1 = 1.0,
+		 * it return 10.9 - (integer)( (10.9 - 0.1) / 1.0) * 1.0
+		 *
+		 * # Example
+		 * \snippet mathutil_unit_test.hpp periodic
+		 *
+		 * \param pos
+		 * \param l
+		 * \param b
+		 *
+		 * \return the bound number
+		 *
+		 */
+		template<typename T> static inline T periodic(const T & pos, const T & p2, const T & p1)
+		{
+			T pos_tmp;
+
+			pos_tmp = pos - (p2 - p1) * (long int)( (pos -p1) / (p2 - p1));
+			pos_tmp += (pos < p1)?(p2 - p1):0;
+
+			return pos_tmp;
+		}
+
+		/*! \brief Bound the position to be inside p2 and p1
+		 *
+		 * It is like periodic but faster
+		 *
+		 * \warning pos should not overshoot the domain more than one time, for example
+		 *          if p2 = 1.1 and p1 = 0.1, pos can be between 2.1 and -0.9
+		 *
+		 * # Example
+		 * \snippet mathutil_unit_test.hpp periodic_l
+		 *
+		 * \param pos
+		 * \param l
+		 * \param b
+		 *
+		 * \return the bound number
+		 *
+		 */
+		template<typename T> static inline T periodic_l(const T & pos, const T & p2, const T & p1)
+		{
+			T pos_tmp = pos;
+
+			if (pos > p2)
+				pos_tmp -= (p2 - p1);
+			else if (pos < p1)
+				pos_tmp += (p2 - p1);
+
+
+			return pos_tmp;
+		}
 	}
 }
 
diff --git a/src/util/mathutil_unit_test.hpp b/src/util/mathutil_unit_test.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..59d67304fb7226bf409550fc059e80170f9818be
--- /dev/null
+++ b/src/util/mathutil_unit_test.hpp
@@ -0,0 +1,116 @@
+/*
+ * mathutil_unit_test.hpp
+ *
+ *  Created on: Dec 23, 2015
+ *      Author: i-bird
+ */
+
+#ifndef OPENFPM_DATA_SRC_UTIL_MATHUTIL_UNIT_TEST_HPP_
+#define OPENFPM_DATA_SRC_UTIL_MATHUTIL_UNIT_TEST_HPP_
+
+#include "mathutil.hpp"
+
+BOOST_AUTO_TEST_SUITE( mathutil_test )
+
+BOOST_AUTO_TEST_CASE( math )
+{
+	//! [factorial usage]
+
+	size_t f = openfpm::math::factorial(0);
+	BOOST_REQUIRE_EQUAL(f,1ul);
+
+	f = openfpm::math::factorial(1);
+	BOOST_REQUIRE_EQUAL(f,1ul);
+
+	f = openfpm::math::factorial(7);
+	BOOST_REQUIRE_EQUAL(f,5040ul);
+
+	//! [factorial usage]
+
+	//! [Combination usage]
+
+	f = openfpm::math::C(7,0);
+	BOOST_REQUIRE_EQUAL(f,1ul);
+
+	f = openfpm::math::C(7,7);
+	BOOST_REQUIRE_EQUAL(f,1ul);
+
+	f = openfpm::math::C(7,2);
+	BOOST_REQUIRE_EQUAL(f,21ul);
+
+	//! [Combination usage]
+
+	//! [round to big pow]
+
+	f = openfpm::math::round_big_2(3);
+	BOOST_REQUIRE_EQUAL(f,4ul);
+
+	f = openfpm::math::round_big_2(7);
+	BOOST_REQUIRE_EQUAL(f,8ul);
+
+	f = openfpm::math::round_big_2(21);
+	BOOST_REQUIRE_EQUAL(f,32ul);
+
+	//! [round to big pow]
+
+	//! [pow]
+
+	f = openfpm::math::pow(3,3);
+	BOOST_REQUIRE_EQUAL(f,27ul);
+
+	f = openfpm::math::pow(4,2);
+	BOOST_REQUIRE_EQUAL(f,16ul);
+
+	f = openfpm::math::pow(2,5);
+	BOOST_REQUIRE_EQUAL(f,32ul);
+
+	//! [pow]
+
+	//! [positive modulo]
+
+	f = openfpm::math::positive_modulo(-1,3);
+	BOOST_REQUIRE_EQUAL(f,2ul);
+
+	f = openfpm::math::positive_modulo(-5,2);
+	BOOST_REQUIRE_EQUAL(f,1ul);
+
+	f = openfpm::math::positive_modulo(-7,5);
+	BOOST_REQUIRE_EQUAL(f,3ul);
+
+	//! [positive modulo]
+
+	//! [periodic]
+
+	float ff = openfpm::math::periodic(10.9,1.1,0.1);
+	BOOST_REQUIRE_CLOSE(ff,0.9,0.0001);
+
+	ff = openfpm::math::periodic(0.0,1.1,0.1);
+	BOOST_REQUIRE_CLOSE(ff,1.0,0.0001);
+
+	ff = openfpm::math::periodic(-7.0,1.1,0.1);
+	BOOST_REQUIRE_CLOSE(ff,1.0,0.0001);
+
+	// in float representation 1.1 is a little bigger than 1.1
+	ff = openfpm::math::periodic(-10.9,1.1,0.1);
+	BOOST_REQUIRE_CLOSE(ff,1.1,0.0001);
+
+	//! [periodic]
+
+	//! [periodic_l]
+
+	ff = openfpm::math::periodic(0.0,1.1,0.1);
+	BOOST_REQUIRE_CLOSE(ff,1.0,0.0001);
+
+	ff = openfpm::math::periodic(1.2,1.1,0.1);
+	BOOST_REQUIRE_CLOSE(ff,0.2,0.0001);
+
+	// in float representation 1.1 is a little bigger than 1.1
+	ff = openfpm::math::periodic(0.9,1.1,0.1);
+	BOOST_REQUIRE_CLOSE(ff,0.9,0.0001);
+
+	//! [periodic_l]
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+#endif /* OPENFPM_DATA_SRC_UTIL_MATHUTIL_UNIT_TEST_HPP_ */
diff --git a/src/util/util_test.hpp b/src/util/util_test.hpp
index dbc71b12f13d00d4304bb854b6a84e5a7e0fbef3..4435a938365b5e3e718354cdbc9a3f547948b199 100644
--- a/src/util/util_test.hpp
+++ b/src/util/util_test.hpp
@@ -323,7 +323,7 @@ BOOST_AUTO_TEST_CASE( generate_array )
 	// ct_test_ce::data is equivalent to constexpr size_t [5] = {5,6,7,8,9}
 
 	const size_t ct_calc = MetaFunc<ct_test_ce::data[0],ct_test_ce::data[1]>::value;
-	BOOST_REQUIRE_EQUAL(ct_calc,11);
+	BOOST_REQUIRE_EQUAL(ct_calc,11ul);
 	//! [constexpr array]
 	}
 
@@ -577,13 +577,13 @@ BOOST_AUTO_TEST_CASE( check_templates_util_function )
 		//! [Usage mul_array_extents]
 
 		size_t mul = array_extents<float>::mul();
-		BOOST_REQUIRE_EQUAL(mul,1);
+		BOOST_REQUIRE_EQUAL(mul,1ul);
 		mul = array_extents<float[3]>::mul();
-		BOOST_REQUIRE_EQUAL(mul,3);
+		BOOST_REQUIRE_EQUAL(mul,3ul);
 		mul = array_extents<float[3][2]>::mul();
-		BOOST_REQUIRE_EQUAL(mul,3*2);
+		BOOST_REQUIRE_EQUAL(mul,3ul*2ul);
 		mul = array_extents<float[3][2][5]>::mul();
-		BOOST_REQUIRE_EQUAL(mul,3*2*5);
+		BOOST_REQUIRE_EQUAL(mul,3ul*2ul*5ul);
 
 		//! [Usage mul_array_extents]
 		}
@@ -601,8 +601,8 @@ BOOST_AUTO_TEST_CASE( check_convert_function )
 
 	Point<2,size_t> p = toPoint<2,size_t>::convert(c);
 
-	BOOST_REQUIRE_EQUAL(p.get(0),1);
-	BOOST_REQUIRE_EQUAL(p.get(1),-1);
+	BOOST_REQUIRE_EQUAL(p.get(0),1ul);
+	BOOST_REQUIRE_EQUAL(p.get(1),(size_t)-1);
 
 	//! [Convert combination to Point]
 	}