diff --git a/images/Makefile.am b/images/Makefile.am
index 62b67035a449140148e57b9baf38032dfe91a1b9..4634f73c345017c981ad1820beae52fbfad84140 100644
--- a/images/Makefile.am
+++ b/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
diff --git a/images/vector.cpp b/images/vector.cpp
index 70395b09008a7156d19bc03aceaca8b95f9992b8..a1dddd3100fb48e0edededc1d821c14b8c029a20 100644
--- a/images/vector.cpp
+++ b/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
diff --git a/images/vector_ghost.py b/images/vector_ghost.py
new file mode 100644
index 0000000000000000000000000000000000000000..bc1c2d78e5da027c510b2707099706b5cc3bbbe8
--- /dev/null
+++ b/images/vector_ghost.py
@@ -0,0 +1,216 @@
+#### 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(...).
diff --git a/openfpm_numerics b/openfpm_numerics
index 834b07c439b8cb308f535fcf43bddf9cc9ac7823..6c6b3e43a260d6196e8b4aecacf7058068378533 160000
--- a/openfpm_numerics
+++ b/openfpm_numerics
@@ -1 +1 @@
-Subproject commit 834b07c439b8cb308f535fcf43bddf9cc9ac7823
+Subproject commit 6c6b3e43a260d6196e8b4aecacf7058068378533
diff --git a/src/Decomposition/CartDecomposition.hpp b/src/Decomposition/CartDecomposition.hpp
index 8454d6d349e0206cf69b47080269956b965374a2..bf2b8fbb06b6389f184fd83c87bae054c0c4825e 100755
--- a/src/Decomposition/CartDecomposition.hpp
+++ b/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)
diff --git a/src/Grid/grid_dist_id.hpp b/src/Grid/grid_dist_id.hpp
index 9902cb07a74e2d92af4097a32f5cfd1030f99a92..18c5d5a2afc3a9a1eb4f7678b3a0c6c96e04da60 100644
--- a/src/Grid/grid_dist_id.hpp
+++ b/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);
diff --git a/src/Grid/grid_dist_id_iterator.hpp b/src/Grid/grid_dist_id_iterator.hpp
index 7c32710e6417b25cf0c46335f2e20616705ec02d..571b6d6acce46c48aeee78c6855821668447b281 100644
--- a/src/Grid/grid_dist_id_iterator.hpp
+++ b/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
diff --git a/src/Grid/grid_dist_id_unit_test_ext_dom.hpp b/src/Grid/grid_dist_id_unit_test_ext_dom.hpp
index ea91ca57e149996059c1c0a89ce2371314bb5b6e..572dde5d0f052e7403f8ef378ff92bdd40bb6778 100644
--- a/src/Grid/grid_dist_id_unit_test_ext_dom.hpp
+++ b/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
diff --git a/src/Grid/grid_dist_key.hpp b/src/Grid/grid_dist_key.hpp
index 7d69d5b20551fcc6c8a6c309dd87cf15c0d615f1..f3de75096e9bdd3d319ff7f614647af4896eef12 100644
--- a/src/Grid/grid_dist_key.hpp
+++ b/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++)
diff --git a/src/Grid/staggered_dist_grid.hpp b/src/Grid/staggered_dist_grid.hpp
index 75809182460dbe54de6068b3885d680ee136f7d1..e82687f7b8b97c700fce5a9b77faa68839bb1d3c 100644
--- a/src/Grid/staggered_dist_grid.hpp
+++ b/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);
}
+ /*! \brief Copy the staggered grid into a normal one
+ *
+ *
+ * \tparam Grid_dst type of the destination Grid
+ * \tparam pos destination grid properties to fill
+ *
+ */
+ template<typename Grid_dst ,unsigned int ... pos> bool to_normal(Grid_dst & g_dst, const Padding<dim> & pd, const long int (& start)[dim], const long int (& stop)[dim])
+ {
+ // interpolation points for each property
+ openfpm::vector<std::vector<comb<dim>>> interp_pos[sizeof...(pos)];
+
+ typedef boost::mpl::vector_c<unsigned int,pos ... > v_pos_type;
+
+ interp_points<dim,v_pos_type,typename Grid_dst::value_type::type> itp(interp_pos,c_prp);
+ boost::mpl::for_each_ref<v_pos_type>(itp);
+
+ // shift the start and stop by the padding
+ grid_key_dx<dim> start_k = grid_key_dx<dim>(start);
+ grid_key_dx<dim> stop_k = grid_key_dx<dim>(stop);
+
+ // sub-grid iterator over the grid map
+ auto g_map_it = this->getSubDomainIterator(start_k,stop_k);
+
+ // Iterator over the destination grid
+ auto g_dst_it = g_dst.getDomainIterator();
+
+ // Check that the 2 iterator has the same size
+ checkIterator<St,decltype(g_map_it),decltype(g_dst_it)>(g_map_it,g_dst_it);
+
+ while (g_map_it.isNext() == true)
+ {
+ typedef typename to_boost_vmpl<pos...>::type vid;
+ typedef boost::mpl::size<vid> v_size;
+
+ auto key_src = g_map_it.get();
+
+ // destination point
+ auto key_dst = g_dst_it.get();
+
+ // Transform this id into an id for the Eigen vector
+
+ interp_ele<vid,Grid_dst,typename std::remove_reference<decltype(*this)>::type,sizeof...(pos)> cp(key_dst,g_dst,*this,key_src,interp_pos);
+
+ // For each property in the target grid
+ boost::mpl::for_each_ref<boost::mpl::range_c<int,0,v_size::value>>(cp);
+
+ ++g_map_it;
+ ++g_dst_it;
+ }
+
+ return true;
+ }
+
/*! \brief Get the staggered positions
*
* \return The vector of the staggered positions
@@ -131,6 +206,16 @@ public:
return c_prp[prp].size() != 0;
}
+ /*! \brief Return if the grid is staggered
+ *
+ * \return true
+ *
+ */
+ bool is_staggered()
+ {
+ return true;
+ }
+
friend class stag_create_and_add_grid<dim,staggered_grid_dist<dim,St,T,Decomposition,Memory,device_grid>,St>;
};
diff --git a/src/Grid/staggered_dist_grid_util.hpp b/src/Grid/staggered_dist_grid_util.hpp
index 5a269985abadcdf4efd169f8d55e733100ad1122..c7d20ec6a1902d11264aa8ceeb7c67119ff368f8 100644
--- a/src/Grid/staggered_dist_grid_util.hpp
+++ b/src/Grid/staggered_dist_grid_util.hpp
@@ -658,4 +658,102 @@ public:
}
};
+/*! \brief Check that the size of the iterators match
+ *
+ * It check the the boxes that the sub iterator defines has same dimensions, for example
+ * if the first sub-iterator, iterate from (1,1) to (5,3) and the second from (2,2) to (6,4)
+ * they match (2,2) to (4,6) they do not match
+ *
+ * \tparam Grid_map type of the map grid
+ * \tparam Grid_dst type of the destination grid
+ *
+ * \param it1 Iterator1
+ * \param it2 Iterator2
+ *
+ * \return true if they match
+ *
+ */
+template<typename Eqs_sys, typename it1_type, typename it2_type> bool checkIterator(const it1_type & it1, const it2_type & it2)
+{
+#ifdef SE_CLASS1
+
+ grid_key_dx<Eqs_sys::dims> it1_k = it1.getStop() - it1.getStart();
+ grid_key_dx<Eqs_sys::dims> it2_k = it2.getStop() - it2.getStart();
+
+ for (size_t i = 0 ; i < Eqs_sys::dims ; i++)
+ {
+ if (it1_k.get(i) != it2_k.get(i))
+ {
+ std::cerr << __FILE__ << ":" << __LINE__ << " error src iterator and destination iterator does not match in size\n";
+ return false;
+ }
+ }
+
+ return true;
+#else
+
+ return true;
+
+#endif
+}
+
+/*! \brief this class is a functor for "for_each" algorithm
+ *
+ * This class is a functor for "for_each" algorithm. For each
+ * element of the boost::vector the operator() is called.
+ * Is mainly used to calculate the interpolation points for each
+ * property in a staggered grid
+ *
+ * \tparam dim Dimensionality
+ * \tparam v_prp_id vector of properties id
+ * \tparam v_prp_type vector with the properties
+ *
+ */
+template<unsigned int dim, typename v_prp_id, typename v_prp_type>
+struct interp_points
+{
+ // number of properties we are processing
+ typedef boost::mpl::size<v_prp_id> v_size;
+
+ // interpolation points for each property
+ openfpm::vector<std::vector<comb<dim>>> (& interp_pts)[v_size::value];
+
+ // staggered position for each property
+ const openfpm::vector<comb<dim>> (&stag_pos)[v_size::value];
+
+ /*! \brief constructor
+ *
+ * It define the copy parameters.
+ *
+ * \param inter_pts array that for each property contain the interpolation points for each components
+ * \param staggered position for each property and components
+ *
+ */
+ inline interp_points(openfpm::vector<std::vector<comb<dim>>> (& interp_pts)[v_size::value],const openfpm::vector<comb<dim>> (&stag_pos)[v_size::value])
+ :interp_pts(interp_pts),stag_pos(stag_pos){};
+
+ //! It call the copy function for each property
+ template<typename T>
+ inline void operator()(T& t)
+ {
+ // This is the type of the object we have to copy
+ typedef typename boost::mpl::at_c<v_prp_type,T::value>::type prp_type;
+
+ interp_pts[T::value].resize(stag_pos[T::value].size());
+
+ for (size_t i = 0 ; i < stag_pos[T::value].size() ; i++)
+ {
+ // Create the interpolation points
+ interp_pts[T::value].get(i) = SubHyperCube<dim,dim - std::rank<prp_type>::value>::getCombinations_R(stag_pos[T::value].get(i),0);
+
+ // interp_point are -1,0,1, map the -1 to 0 and 1 to -1
+ for (size_t j = 0 ; j < interp_pts[T::value].get(i).size() ; j++)
+ {
+ for (size_t k = 0 ; k < dim ; k++)
+ interp_pts[T::value].get(i)[j].getComb()[k] = - ((interp_pts[T::value].get(i)[j].getComb()[k] == -1)?0:interp_pts[T::value].get(i)[j].getComb()[k]);
+ }
+ }
+ }
+};
+
#endif /* SRC_GRID_STAGGERED_DIST_GRID_UTIL_HPP_ */
diff --git a/src/Vector/vector_dist.hpp b/src/Vector/vector_dist.hpp
index 44540b8e9e3ee878b767819a64166fc492d04bfe..ab06e71d64f6d68c860d26e673e31290d5f9957c 100644
--- a/src/Vector/vector_dist.hpp
+++ b/src/Vector/vector_dist.hpp
@@ -262,8 +262,8 @@ private:
| v | (0,0) --> 4 | v |
| 3 | | 5 |
| | | |
- B | | | A |
- * | | | * |
+ B | | | A |
+ * | | | * |
| | | |
| | | |
| | | |
diff --git a/src/gargabe.hpp b/src/gargabe.hpp
index 20e091ad5b9e72b1b994f9a9e3143e9fc5eeb4f2..eb05c7e75b013c0680563c01326c4ba079180d9d 100644
--- a/src/gargabe.hpp
+++ b/src/gargabe.hpp
@@ -621,4 +621,242 @@ public:
virtual void destroy(size_t n) = 0;
};*/
+
+
+
+
+
+
+
+
+/*! \brief Impose an operator
+ *
+ * This function impose an operator on a particular grid region to produce the system
+ *
+ * Ax = b
+ *
+ * ## Stokes equation, lid driven cavity with one splipping wall
+ *
+ * \param op Operator to impose (A term)
+ * \param num right hand side of the term (b term)
+ * \param id Equation id in the system that we are imposing
+ * \param it_d iterator that define where you want to impose
+ *
+ */
+template<typename T> void impose(const T & op , typename Sys_eqs::stype num ,long int id ,grid_dist_iterator_sub<Sys_eqs::dims,typename g_map_type::d_grid> it_d, bool skip_first = false)
+{
+ //////////////////////// DEBUG /////////////////
+
+ SparseMatrix<double,int> Al;
+ Al.load("debug_matrix_single_processor");
+
+ // Construct the map 3 processors 1 processors
+
+ std::unordered_map<size_t,size_t> map_row;
+
+ auto it2 = g_map.getDomainGhostIterator();
+ auto ginfo = g_map.getGridInfoVoid();
+
+ while (it2.isNext())
+ {
+ auto key = it2.get();
+ auto key_g = g_map.getGKey(key);
+ key_g += pd.getKP1();
+
+ // To linearize must be positive
+ bool is_negative = false;
+ for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
+ {
+ if (key_g.get(i) < 0)
+ is_negative = true;
+ }
+
+ if (is_negative == true)
+ {
+ ++it2;
+ continue;
+ }
+
+ // Carefull g map is extended, so the original (0,0) is shifted in g_map by
+
+ if (g_map.template get<0>(key) == 7)
+ {
+ int debug = 0;
+ debug++;
+ }
+
+ map_row[g_map.template get<0>(key)] = ginfo.LinId(key_g);
+
+ ++it2;
+ }
+
+ ////////////////////////////////////////////////
+
+ Vcluster & v_cl = *global_v_cluster;
+
+ openfpm::vector<triplet> & trpl = A.getMatrixTriplets();
+
+ auto it = it_d;
+ grid_sm<Sys_eqs::dims,void> gs = g_map.getGridInfoVoid();
+
+ std::unordered_map<long int,float> cols;
+
+ // resize b if needed
+ b.resize(Sys_eqs::nvar * g_map.size());
+
+ bool is_first = skip_first;
+
+ // iterate all the grid points
+ while (it.isNext())
+ {
+ if (is_first == true && v_cl.getProcessUnitID() == 0)
+ {
+ ++it;
+ is_first = false;
+ continue;
+ }
+ else
+ is_first = false;
+
+ // get the position
+ auto key = it.get();
+
+ // Calculate the non-zero colums
+ T::value(g_map,key,gs,spacing,cols,1.0);
+
+ //////////// DEBUG //////////////////
+
+ auto g_calc_pos = g_map.getGKey(key);
+ g_calc_pos += pd.getKP1();
+
+ /////////////////////////////////////
+
+ // create the triplet
+
+ for ( auto it = cols.begin(); it != cols.end(); ++it )
+ {
+ trpl.add();
+ trpl.last().row() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
+ trpl.last().col() = it->first;
+ trpl.last().value() = it->second;
+
+ ///////////// DEBUG ///////////////////////
+
+ auto ginfo = g_map.getGridInfoVoid();
+
+ size_t r = (trpl.last().row() / Sys_eqs::nvar);
+ size_t r_rest = (trpl.last().row() % Sys_eqs::nvar);
+ size_t c = (trpl.last().col() / Sys_eqs::nvar);
+ size_t c_rest = (trpl.last().col() % Sys_eqs::nvar);
+ double val = trpl.last().value();
+
+ // Transform
+
+ size_t rf = map_row[r] * 3 + r_rest;
+ size_t cf = map_row[c] * 3 + c_rest;
+
+ auto position_row = ginfo.InvLinId(rf / 3);
+ auto position_col = ginfo.InvLinId(cf / 3);
+
+ double valf = Al.getValue(rf,cf);
+
+ if (val != valf)
+ {
+ int debug = 0;
+ debug++;
+ }
+
+ ///////////////////////////////////////////
+
+// std::cout << "(" << trpl.last().row() << "," << trpl.last().col() << "," << trpl.last().value() << ")" << "\n";
+ }
+
+ b(g_map.template get<0>(key)*Sys_eqs::nvar + id) = num;
+
+ cols.clear();
+// std::cout << "\n";
+
+ // if SE_CLASS1 is defined check the position
+#ifdef SE_CLASS1
+// T::position(key,gs,s_pos);
+#endif
+
+ ++row;
+ ++row_b;
+ ++it;
+ }
+}
+
+typename Sys_eqs::SparseMatrix_type A;
+
+/*! \brief produce the Matrix
+ *
+ * \return the Sparse matrix produced
+ *
+ */
+typename Sys_eqs::SparseMatrix_type & getA()
+{
+#ifdef SE_CLASS1
+ consistency();
+#endif
+ A.resize(g_map.size()*Sys_eqs::nvar,g_map.size()*Sys_eqs::nvar);
+
+ ///////////////// DEBUG SAVE //////////////////
+
+// A.save("debug_matrix_single_processor");
+
+ ////////////////////////////////////////////////
+
+ return A;
+
+}
+
+
+typename Sys_eqs::SparseMatrix_type A;
+
+/*! \brief produce the Matrix
+ *
+ * \return the Sparse matrix produced
+ *
+ */
+typename Sys_eqs::SparseMatrix_type & getA()
+{
+#ifdef SE_CLASS1
+ consistency();
+#endif
+ A.resize(g_map.size()*Sys_eqs::nvar,g_map.size()*Sys_eqs::nvar);
+
+ ///////////////// DEBUG SAVE //////////////////
+
+// A.save("debug_matrix_single_processor");
+
+ ////////////////////////////////////////////////
+
+ return A;
+
+}
+
+
+/*! \brief produce the B vector
+ *
+ * \return the vector produced
+ *
+ */
+typename Sys_eqs::Vector_type & getB()
+{
+#ifdef SE_CLASS1
+ consistency();
+#endif
+
+ // size of the matrix
+// B.resize(g_map.size()*Sys_eqs::nvar);
+
+ // copy the vector
+// for (size_t i = 0; i < row_b; i++)
+// B.insert(i,b.get(i));
+
+ return b;
+}
+};
+
#endif /* GARGABE_HPP_ */