Adding ghost particles images

images/Makefile.am

@@ -48,6 +48,7 @@ images : cart_dec metis_dec dom_box
 	pvbatch vector_scal_vect.py
 	pvbatch vector_particles.py
 	pvbatch particles_maps.py
+	pvbatch vector_ghost.py
 	dot -Tsvg openfpm.dot -o generated/openfpm.svg
 	avconv -i generated/particles_mooving.ogv -f mp4 generated/particles_mooving.mp4
 	avconv -i generated/particles_mooving_prc.ogv -f mp4 generated/particles_mooving_prc.mp4

images/vector.cpp

@@ -61,6 +61,10 @@ int main(int argc, char* argv[])
 
 	vd.write("Vector/vector_after_map");
 	
+	vd.ghost_get<0>();
+
+	vd.write("Vector/vector_ghost_fill");
+
 	vd.getDecomposition().write("Vector/vect_decomposition");
 
 	// move the particles

images/vector_ghost.py

+#### import the simple module from the paraview
+from paraview.simple import *
+#### disable automatic camera reset on 'Show'
+paraview.simple._DisableFirstRenderCameraReset()
+
+# create a new 'CSV Reader'
+vector_ghost_fill = CSVReader(FileName=['Vector/vector_ghost_fill0.csv'])
+
+# get animation scene
+animationScene1 = GetAnimationScene()
+
+# update animation scene based on data timesteps
+animationScene1.UpdateAnimationUsingDataTimeSteps()
+
+# get active view
+renderView1 = GetActiveViewOrCreate('RenderView')
+# uncomment following to set a specific view size
+renderView1.ViewSize = [400, 400]
+
+# get layout
+viewLayout1 = GetLayout()
+
+# Create a new 'SpreadSheet View'
+spreadSheetView1 = CreateView('SpreadSheetView')
+spreadSheetView1.BlockSize = 1024L
+# uncomment following to set a specific view size
+# spreadSheetView1.ViewSize = [400, 400]
+
+# place view in the layout
+viewLayout1.AssignView(2, spreadSheetView1)
+
+# show data in view
+vector_ghost_fillDisplay = Show(vector_ghost_fill, spreadSheetView1)
+
+# destroy spreadSheetView1
+Delete(spreadSheetView1)
+del spreadSheetView1
+
+# close an empty frame
+viewLayout1.Collapse(2)
+
+# set active view
+SetActiveView(renderView1)
+
+# create a new 'Table To Points'
+tableToPoints1 = TableToPoints(Input=vector_ghost_fill)
+
+# Properties modified on tableToPoints1
+tableToPoints1.XColumn = 'x[0]'
+tableToPoints1.YColumn = 'x[1]'
+tableToPoints1.a2DPoints = 1
+
+# show data in view
+#tableToPoints1Display = Show(tableToPoints1, renderView1)
+
+# reset view to fit data
+renderView1.ResetCamera()
+
+#changing interaction mode based on data extents
+renderView1.InteractionMode = '2D'
+renderView1.CameraPosition = [0.499922575, 0.18749537000000002, 10000.0]
+renderView1.CameraFocalPoint = [0.499922575, 0.18749537000000002, 0.0]
+
+# create a new 'Legacy VTK Reader'
+vect_decompositionexternal_ghost_0vtk = LegacyVTKReader(FileNames=['Vector/vect_decompositionexternal_ghost_0.vtk'])
+
+# destroy vect_decompositionexternal_ghost_0vtk
+Delete(vect_decompositionexternal_ghost_0vtk)
+del vect_decompositionexternal_ghost_0vtk
+
+# create a new 'CSV Reader'
+vector_after_map0csv = CSVReader(FileName=['Vector/vector_after_map0.csv'])
+
+# Create a new 'SpreadSheet View'
+spreadSheetView1 = CreateView('SpreadSheetView')
+spreadSheetView1.BlockSize = 1024L
+# uncomment following to set a specific view size
+# spreadSheetView1.ViewSize = [400, 400]
+
+# place view in the layout
+viewLayout1.AssignView(2, spreadSheetView1)
+
+# show data in view
+vector_after_map0csvDisplay = Show(vector_after_map0csv, spreadSheetView1)
+
+# destroy spreadSheetView1
+Delete(spreadSheetView1)
+del spreadSheetView1
+
+# close an empty frame
+viewLayout1.Collapse(2)
+
+# set active view
+SetActiveView(renderView1)
+
+# create a new 'Table To Points'
+tableToPoints2 = TableToPoints(Input=vector_after_map0csv)
+
+# Properties modified on tableToPoints2
+tableToPoints2.XColumn = 'x[0]'
+tableToPoints2.YColumn = 'x[1]'
+tableToPoints2.a2DPoints = 1
+
+# show data in view
+tableToPoints2Display = Show(tableToPoints2, renderView1)
+
+tableToPoints2Display.PointSize = 4.0
+
+# create a new 'Legacy VTK Reader'
+vect_decompositionsubdomains_0vtk = LegacyVTKReader(FileNames=['Vector/vect_decompositionsubdomains_0.vtk'])
+
+# show data in view
+vect_decompositionsubdomains_0vtkDisplay = Show(vect_decompositionsubdomains_0vtk, renderView1)
+
+# show color bar/color legend
+vect_decompositionsubdomains_0vtkDisplay.SetScalarBarVisibility(renderView1, True)
+
+# get color transfer function/color map for 'data'
+dataLUT = GetColorTransferFunction('data')
+
+# get opacity transfer function/opacity map for 'data'
+dataPWF = GetOpacityTransferFunction('data')
+
+# hide data in view
+Hide(tableToPoints1, renderView1)
+
+# change representation type
+vect_decompositionsubdomains_0vtkDisplay.SetRepresentationType('Points')
+
+# change representation type
+vect_decompositionsubdomains_0vtkDisplay.SetRepresentationType('Wireframe')
+
+# turn off scalar coloring
+ColorBy(vect_decompositionsubdomains_0vtkDisplay, None)
+
+# set active source
+SetActiveSource(tableToPoints2)
+
+# change representation type
+tableToPoints2Display.SetRepresentationType('Points')
+
+# current camera placement for renderView1
+renderView1.InteractionMode = '2D'
+renderView1.CameraPosition = [0.45889311126344623, 0.2866467748466169, 10000.0]
+renderView1.CameraFocalPoint = [0.45889311126344623, 0.2866467748466169, 0.0]
+renderView1.CameraParallelScale = 0.1742200032569875
+
+# create a new 'Legacy VTK Reader'
+vect_decompositionexternal_ghost_0vtk = LegacyVTKReader(FileNames=['Vector/vect_decompositionexternal_ghost_0.vtk'])
+
+# show data in view
+vect_decompositionexternal_ghost_0vtkDisplay = Show(vect_decompositionexternal_ghost_0vtk, renderView1)
+
+# show color bar/color legend
+vect_decompositionexternal_ghost_0vtkDisplay.SetScalarBarVisibility(renderView1, False)
+
+# hide data in view
+Hide(vect_decompositionexternal_ghost_0vtk, renderView1)
+
+# current camera placement for renderView1
+renderView1.InteractionMode = '2D'
+renderView1.CameraPosition = [0.4322940474199577, 0.3362326839869475, 10000.0]
+renderView1.CameraFocalPoint = [0.4322940474199577, 0.3362326839869475, 0.0]
+renderView1.CameraParallelScale = 0.08127491729954842
+
+# save screenshot
+WriteImage('generated/vector_one_p_zoom.jpg')
+
+# set active source
+SetActiveSource(vect_decompositionexternal_ghost_0vtk)
+
+# show data in view
+vect_decompositionexternal_ghost_0vtkDisplay = Show(vect_decompositionexternal_ghost_0vtk, renderView1)
+
+# show color bar/color legend
+vect_decompositionexternal_ghost_0vtkDisplay.SetScalarBarVisibility(renderView1, True)
+
+# current camera placement for renderView1
+renderView1.InteractionMode = '2D'
+renderView1.CameraPosition = [0.4322940474199577, 0.3362326839869475, 10000.0]
+renderView1.CameraFocalPoint = [0.4322940474199577, 0.3362326839869475, 0.0]
+renderView1.CameraParallelScale = 0.08127491729954842
+
+# save screenshot
+WriteImage('generated/vector_one_p_zoom_ghost.jpg')
+
+# set active source
+SetActiveSource(tableToPoints1)
+
+# show data in view
+tableToPoints1Display = Show(tableToPoints1, renderView1)
+
+# change representation type
+#tableToPoints1Display.SetRepresentationType('Points')
+
+# Properties modified on tableToPoints1Display
+tableToPoints1Display.PointSize = 4.0
+
+# change solid color
+tableToPoints1Display.AmbientColor = [1.0, 0.0, 0.0]
+
+# current camera placement for renderView1
+renderView1.InteractionMode = '2D'
+renderView1.CameraPosition = [0.4322940474199577, 0.3362326839869475, 10000.0]
+renderView1.CameraFocalPoint = [0.4322940474199577, 0.3362326839869475, 0.0]
+renderView1.CameraParallelScale = 0.08127491729954842
+
+#### saving camera placements for all active views
+
+# save screenshot
+WriteImage('generated/vector_one_p_zoom_ghost_part.jpg')
+
+
+#### uncomment the following to render all views
+# RenderAllViews()
+# alternatively, if you want to write images, you can use SaveScreenshot(...). 

openfpm_numerics @ 6c6b3e43

-Subproject commit 834b07c439b8cb308f535fcf43bddf9cc9ac7823
+Subproject commit 6c6b3e43a260d6196e8b4aecacf7058068378533

src/Decomposition/CartDecomposition.hpp

@@ -938,8 +938,10 @@ public:
 
 	/*! \brief Set the parameter of the decomposition
 	 *
-	 * \param div_ storing into how many domain to decompose on each dimension
+	 * \param div_ storing into how many sub-sub-domains to decompose on each dimension
 	 * \param domain_ domain to decompose
+	 * \param bc_ boundary conditions
+	 * \param ghost Ghost size
 	 *
 	 */
 	void setParameters(const size_t (& div_)[dim], ::Box<dim,T> domain_, const size_t (& bc)[dim] ,const Ghost<dim,T> & ghost)

src/Grid/grid_dist_id.hpp

@@ -545,6 +545,37 @@ public:
 	// value_type
 	typedef T value_type;
 
+	// Type of space
+	typedef St stype;
+
+	// Type of Memory
+	typedef Memory memory_type;
+
+	// Type of device grid
+	typedef device_grid device_grid_type;
+
+	// Number of dimensions
+	static const unsigned int dims = dim;
+
+	/*! \brief Get the domain where the grid is defined
+	 *
+	 *
+	 */
+	inline const Box<dim,St> getDomain() const
+	{
+		return domain;
+	}
+
+    /*! \brief Get the spacing of the grid in direction i
+     *
+     * \return the spacing
+     *
+     */
+    inline St spacing(size_t i) const
+    {
+    	return cd_sm.getCellBox().getHigh(i);
+    }
+
 	/*! \brief Return the total number of points in the grid
 	 *
 	 * \return number of points
@@ -593,11 +624,12 @@ public:
 	 * something similar, but because of rounding-off error it can happen that it is not perfectly overlapping
 	 *
 	 * \param g previous grid
+	 * \param Ghost part in grid units
 	 * \param ext extension of the grid (must be positive on every direction)
 	 *
 	 */
-	grid_dist_id(const grid_dist_id<dim,St,T,typename Decomposition::base_type,Memory,device_grid> & g, Box<dim,size_t> ext)
-	:ghost(g.getDecomposition().getGhost()),dec(*global_v_cluster),v_cl(*global_v_cluster)
+	template<typename H> grid_dist_id(const grid_dist_id<dim,St,H,typename Decomposition::base_type,Memory,grid_cpu<dim,H>> & g, const Ghost<dim,long int> & gh, Box<dim,size_t> ext)
+	:dec(*global_v_cluster),v_cl(*global_v_cluster)
 	{
 #ifdef SE_CLASS2
 		check_new(this,8,GRID_DIST_EVENT,4);
@@ -605,26 +637,42 @@ public:
 
 		this->dec.incRef();
 
-		InitializeCellDecomposer(g.cd_sm,ext);
+		size_t ext_dim[dim];
+		for (size_t i = 0 ; i < dim ; i++) {ext_dim[i] = g.getGridInfoVoid().size(i) + ext.getKP1().get(i) + ext.getKP2().get(i);}
+
+		// Set the grid info of the extended grid
+		ginfo.setDimensions(ext_dim);
+		ginfo_v.setDimensions(ext_dim);
+
+		InitializeCellDecomposer(g.getCellDecomposer(),ext);
+
+		ghost = convert_ghost(gh,cd_sm);
 
 		// Extend the grid by the extension part and calculate the domain
 
 		for (size_t i = 0 ; i < dim ; i++)
 		{
-			g_sz[i] = g.g_sz[i] + ext.getLow(i) + ext.getHigh(i);
+			g_sz[i] = g.size(i) + ext.getLow(i) + ext.getHigh(i);
 
-			this->domain.setLow(i,g.domain.getLow(i) - ext.getLow(i) * g.spacing(i) - g.spacing(i) / 2.0);
-			this->domain.setHigh(i,g.domain.getHigh(i) + ext.getHigh(i) * g.spacing(i) + g.spacing(i) / 2.0);
+			this->domain.setLow(i,g.getDomain().getLow(i) - ext.getLow(i) * g.spacing(i) - g.spacing(i) / 2.0);
+			this->domain.setHigh(i,g.getDomain().getHigh(i) + ext.getHigh(i) * g.spacing(i) + g.spacing(i) / 2.0);
 		}
 
-		dec.setParameters(g.getDecomposition(),g.getDecomposition().getGhost(),this->domain);
+		dec.setParameters(g.getDecomposition(),ghost,this->domain);
 
-		InitializeStructures(g_sz);
+		InitializeStructures(g.getGridInfoVoid().getSize());
 	}
 
-    //! constructor
+    /*! It constructs a grid of a specified size, defined on a specified Box space, forcing to follow a specified decomposition and with a specified ghost size
+     *
+     * \param dec Decomposition
+     * \param g_sz grid size on each dimension
+     * \param domain Box that contain the grid
+     * \param ghost Ghost part
+     *
+     */
     grid_dist_id(const Decomposition & dec, const size_t (& g_sz)[dim], const Box<dim,St> & domain, const Ghost<dim,St> & ghost)
-    :domain(domain),ghost(ghost),dec(dec),v_cl(*global_v_cluster)
+    :domain(domain),ghost(ghost),dec(dec),v_cl(*global_v_cluster),ginfo(g_sz),ginfo_v(g_sz)
 	{
 		// Increment the reference counter of the decomposition
 		this->dec.incRef();
@@ -633,9 +681,16 @@ public:
 		InitializeStructures(g_sz);
 	}
 
-    //! constructor
+    /*! It constructs a grid of a specified size, defined on a specified Box space, forcing to follow a specified decomposition and with a specified ghost size
+     *
+     * \param dec Decomposition
+     * \param g_sz grid size on each dimension
+     * \param domain Box that contain the grid
+     * \param ghost Ghost part
+     *
+     */
     grid_dist_id(Decomposition && dec, const size_t (& g_sz)[dim], const Box<dim,St> & domain, const Ghost<dim,St> & ghost)
-    :domain(domain),ghost(ghost),dec(dec),v_cl(*global_v_cluster)
+    :domain(domain),ghost(ghost),dec(dec),ginfo(g_sz),ginfo_v(g_sz),v_cl(*global_v_cluster)
 	{
 #ifdef SE_CLASS2
 		check_new(this,8,GRID_DIST_EVENT,4);
@@ -645,23 +700,16 @@ public:
 		InitializeStructures(g_sz);
 	}
 
-    /*! \brief Get the spacing of the grid in direction i
+    /*! It constructs a grid of a specified size, defined on a specified Box space, forcing to follow a specified decomposition, and having a specified ghost size
      *
-     * \return the spacing
+     * \param dec Decomposition
+     * \param g_sz grid size on each dimension
+     * \param domain Box that contain the grid
+     * \param ghost Ghost part (given in grid units)
+     *
+     * \warning In very rare case the ghost part can be one point bigger than the one specified
      *
      */
-    inline St spacing(size_t i) const
-    {
-    	return cd_sm.getCellBox().getHigh(i);
-    }
-
-	/*! \brief Constrcuctor
-	 *
-	 * \param g_sz array with the grid size on each dimension
-	 * \param domain domain where this grid live
-	 * \param g Ghost given in grid units
-	 *
-	 */
 	grid_dist_id(const Decomposition & dec, const size_t (& g_sz)[dim],const Box<dim,St> & domain, const Ghost<dim,long int> & g)
 	:domain(domain),dec(*global_v_cluster),v_cl(*global_v_cluster),ginfo(g_sz),ginfo_v(g_sz)
 	{
@@ -678,13 +726,16 @@ public:
 		InitializeStructures(g_sz);
 	}
 
-	/*! \brief Constrcuctor
-	 *
-	 * \param g_sz array with the grid size on each dimension
-	 * \param domain domain where this grid live
-	 * \param g Ghost given in grid units
-	 *
-	 */
+    /*! It construct a grid of a specified size, defined on a specified Box space, forcing to follow a specified decomposition, and having a specified ghost size
+     *
+     * \param dec Decomposition
+     * \param g_sz grid size on each dimension
+     * \param domain Box that contain the grid
+     * \param ghost Ghost part (given in grid units)
+     *
+     * \warning In very rare case the ghost part can be one point bigger than the one specified
+     *
+     */
 	grid_dist_id(Decomposition && dec, const size_t (& g_sz)[dim],const Box<dim,St> & domain, const Ghost<dim,long int> & g)
 	:domain(domain),dec(dec),v_cl(*global_v_cluster),ginfo(g_sz),ginfo_v(g_sz)
 	{
@@ -699,13 +750,15 @@ public:
 		InitializeStructures(g_sz);
 	}
 
-	/*! \brief Constrcuctor
-	 *
-	 * \param g_sz array with the grid size on each dimension
-	 * \param domain domain where this grid live
-	 * \param g Ghost
-	 *
-	 */
+    /*! It construct a grid of a specified size, defined on a specified Box space, and having a specified ghost size
+     *
+     * \param g_sz grid size on each dimension
+     * \param domain Box that contain the grid
+     * \param ghost Ghost part (given in grid units)
+     *
+     * \warning In very rare case the ghost part can be one point bigger than the one specified
+     *
+     */
 	grid_dist_id(const size_t (& g_sz)[dim],const Box<dim,St> & domain, const Ghost<dim,St> & g)
 	:domain(domain),ghost(g),dec(*global_v_cluster),v_cl(*global_v_cluster),ginfo(g_sz),ginfo_v(g_sz)
 	{
@@ -720,13 +773,16 @@ public:
 		InitializeStructures(g_sz);
 	}
 
-	/*! \brief Constrcuctor
-	 *
-	 * \param g_sz array with the grid size on each dimension
-	 * \param domain domain where this grid live
-	 * \param g Ghost given in grid units
-	 *
-	 */
+    /*! It construct a grid of a specified size, defined on a specified Box space,  and having a specified ghost size
+     *
+     * \param dec Decomposition
+     * \param g_sz grid size on each dimension
+     * \param domain Box that contain the grid
+     * \param ghost Ghost part of the domain (given in grid units)
+     *
+     * \warning In very rare case the ghost part can be one point bigger than the one specified
+     *
+     */
 	grid_dist_id(const size_t (& g_sz)[dim],const Box<dim,St> & domain, const Ghost<dim,long int> & g)
 	:domain(domain),dec(*global_v_cluster),v_cl(*global_v_cluster),ginfo(g_sz),ginfo_v(g_sz)
 	{
@@ -748,6 +804,8 @@ public:
 	 * \param domain domain where this grid live
 	 * \param g Ghost given in grid units
 	 *
+	 * \warning In very rare case the ghost part can be one point bigger than the one specified
+	 *
 	 */
 	grid_dist_id(const Decomposition & dec, const std::vector<size_t> & g_sz,const Box<dim,St> & domain, const Ghost<dim,long int> & g)
 	:grid_dist_id(dec,*static_cast<const size_t(*) [dim]>(static_cast<const void*>(&g_sz[0])),domain,g)
@@ -761,6 +819,8 @@ public:
 	 * \param domain domain where this grid live
 	 * \param g Ghost given in grid units
 	 *
+	 * \warning In very rare case the ghost part can be one point bigger than the one specified
+	 *
 	 */
 	grid_dist_id(Decomposition && dec,const std::vector<size_t> & g_sz,const Box<dim,St> & domain, const Ghost<dim,long int> & g)
 	:grid_dist_id(dec, *static_cast<const size_t(*) [dim]>(static_cast<const void*>(&g_sz[0])) , domain, g)
@@ -824,7 +884,7 @@ public:
 	 * \return the cell decomposer
 	 *
 	 */
-	const CellDecomposer_sm<dim,St,shift<dim,St>> & getCellDecomposer()
+	const CellDecomposer_sm<dim,St,shift<dim,St>> & getCellDecomposer() const
 	{
 #ifdef SE_CLASS2
 		check_valid(this,8);
@@ -856,7 +916,7 @@ public:
 	 * \return The size of the local domain
 	 *
 	 */
-	size_t getLocalDomainSize()
+	size_t getLocalDomainSize() const
 	{
 #ifdef SE_CLASS2
 		check_valid(this,8);
@@ -889,7 +949,7 @@ public:
 	 * \return the iterator
 	 *
 	 */
-	grid_dist_iterator<dim,device_grid,FREE> getDomainIterator()
+	grid_dist_iterator<dim,device_grid,FREE> getDomainIterator() const
 	{
 #ifdef SE_CLASS2
 		check_valid(this,8);

src/Grid/grid_dist_id_iterator.hpp

@@ -88,7 +88,7 @@ class grid_dist_iterator<dim,device_grid,FREE>
 	const openfpm::vector<device_grid> & gList;
 
 	//! Extension of each grid: domain and ghost + domain
-	openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext;
+	const openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext;
 
 	//! Actual iterator
 	grid_key_dx_iterator_sub<dim> a_it;
@@ -118,7 +118,7 @@ class grid_dist_iterator<dim,device_grid,FREE>
 	 * \param gk std::vector of the local grid
 	 *
 	 */
-	grid_dist_iterator(const openfpm::vector<device_grid> & gk, openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext, grid_key_dx<dim> stop)
+	grid_dist_iterator(const openfpm::vector<device_grid> & gk, const openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext, grid_key_dx<dim> stop)
 	:g_c(0),gList(gk),gdb_ext(gdb_ext),stop(stop)
 	{
 		// Initialize the current iterator

src/Grid/grid_dist_id_unit_test_ext_dom.hpp

@@ -50,12 +50,15 @@ void Test3D_extended_grid(const Box<3,float> & domain, long int k)
 		// Extend the grid by 2 points
 		Box<3,size_t> ext({2,2,2},{2,2,2});
 
+		// Ghost size of 1 grid point
+		Ghost<3,long int> gp(1);
+
 		// another grid perfectly overlapping the previous, extended by 2 points
-		grid_dist_id<3, float, aggregate<size_t[3],size_t>, CartDecomposition_ext<3,float>> g_dist2(g_dist1,ext);
+		grid_dist_id<3, float, aggregate<size_t[3],size_t>, CartDecomposition_ext<3,float>> g_dist2(g_dist1,gp,ext);
 
-		// Given an iterator on grid 1
+		// Given an iterator of the grid 1
 		auto dom_g1 = g_dist1.getDomainIterator();
-		// And a sub-iterator on grid 2 overlapping grid 1
+		// And a sub-iterator of grid 2 overlapping grid 1
 		auto dom_g2 = g_dist2.getSubDomainIterator({0,0,0},{k-1,k-1,k-1});
 
 		// the 2 iterator must match

src/Grid/grid_dist_key.hpp

@@ -77,7 +77,7 @@ public:
 	 * \return new key
 	 *
 	 */
-	inline grid_dist_key_dx<dim> move(size_t i,size_t s)
+	inline grid_dist_key_dx<dim> move(size_t i,size_t s) const
 	{
 		grid_key_dx<dim> key = getKey();
 		key.set_d(i,key.get(i) + s);
@@ -91,7 +91,7 @@ public:
 	 * \return new key
 	 *
 	 */
-	inline grid_dist_key_dx<dim> move(const comb<dim> & c)
+	inline grid_dist_key_dx<dim> move(const comb<dim> & c) const
 	{
 		grid_key_dx<dim> key = getKey();
 		for (size_t i = 0 ; i < dim ; i++)

src/Grid/staggered_dist_grid.hpp

@@ -11,7 +11,8 @@
 #include "Grid/grid_dist_id.hpp"
 #include "staggered_dist_grid_util.hpp"
 #include "VTKWriter/VTKWriter.hpp"
-
+#include "staggered_dist_grid_copy.hpp"
+#include <boost/mpl/vector_c.hpp>
 
 /*! \brief Implementation of the staggered grid
  *
@@ -67,6 +68,26 @@ public:
 
 	typedef T value_type;
 
+	// Number of dimensions
+	static const unsigned int dims = dim;
+
+	/*! \brief This constructor is special, it construct an expanded grid that perfectly overlap with the previous
+	 *
+	 * The key-word here is "perfectly overlap". Using the default constructor you could create
+	 * something similar, but because of rounding-off error it can happen that it is not perfectly overlapping
+	 *
+	 * \param g previous grid
+	 * \param Ghost part in grid units
+	 * \param ext extension of the grid (must be positive on every direction)
+	 *
+	 */
+	template<typename H> staggered_grid_dist(const grid_dist_id<dim,St,H,typename Decomposition::base_type,Memory,grid_cpu<dim,H>> & g,
+			                          const Ghost<dim,long int> & gh,
+									  Box<dim,size_t> ext)
+	:grid_dist_id<dim,St,T,Decomposition,Memory,device_grid>(g,gh,ext)
+	{
+	}
+
 	staggered_grid_dist(const size_t (& g_sz)[dim], const Box<dim,St> & domain, const Ghost<dim,St> & ghost)
 	:grid_dist_id<dim,St,T,Decomposition,Memory,device_grid>(g_sz,domain,ghost)
 	{}
@@ -97,6 +118,60 @@ public:
 		boost::mpl::for_each_ref< boost::mpl::range_c<int,0,T::max_prop> >(ssp);
 	}