diff --git a/configure.ac b/configure.ac
index bfa1a6701ce24e9a1e83db7a6f361578b0b64dba..3e1da5c5f3efef316e7bc0476c870cfa054506bc 100755
--- a/configure.ac
+++ b/configure.ac
@@ -118,7 +118,7 @@ fi
####### include OpenFPM_devices include path
-INCLUDES_PATH+="-I. -Isrc/config/ -I../openfpm_io/src -I../openfpm_data/src -I../openfpm_devices/src -I../openfpm_vcluster/src/"
+INCLUDES_PATH+="-I. -Iconfig/ -I../openfpm_io/src -I../openfpm_data/src -I../openfpm_devices/src -I../openfpm_vcluster/src/"
###### Check for se-class1
diff --git a/install b/install
index ecc631cb625ce8f1eb9f4722198b06438d3d9ffb..5ea13c9b4de039498ea5cfbda54ed3b1e8dade5e 100755
--- a/install
+++ b/install
@@ -130,23 +130,6 @@ EIGEN_installed=0
blas_options=""
conf_err=1
-## MPI
-
-if [ x"$MPI_installation_required" == x"yes" ]; then
- ./script/install_MPI.sh $i_dir $compiler_opt
- export PATH="$i_dir/MPI/bin:$PATH"
- configure_options="$configure_options CXX=mpic++ "
- MPI_installed=1
-else
- command -v mpic++ >/dev/null 2>&1
- if [ $? -eq 0 ]; then
- configure_options="$configure_options CXX=mpic++ "
- fi
-fi
-
-
-echo "./configure $options $configure_options"
-
if [ $install_req -eq 0 ]; then
./configure $options $configure_options "$blas_options"
else
diff --git a/openfpm_data b/openfpm_data
index 24cff9e625585d446999e55c9a5bcf173333aea8..fb93d17019fdfb1f7c65c463143f5100d506d6d2 160000
--- a/openfpm_data
+++ b/openfpm_data
@@ -1 +1 @@
-Subproject commit 24cff9e625585d446999e55c9a5bcf173333aea8
+Subproject commit fb93d17019fdfb1f7c65c463143f5100d506d6d2
diff --git a/openfpm_io b/openfpm_io
index 9154ba930b68bd58a065f1809b037e0a2320c2c8..c54eca6ee2e0bbfefda778df50b369ce9a90323f 160000
--- a/openfpm_io
+++ b/openfpm_io
@@ -1 +1 @@
-Subproject commit 9154ba930b68bd58a065f1809b037e0a2320c2c8
+Subproject commit c54eca6ee2e0bbfefda778df50b369ce9a90323f
diff --git a/openfpm_numerics b/openfpm_numerics
index dfb8b8f7b56e455c3e1db0ff39e5d513055eb472..7be6531420de26986cd584e3bec0f5ee51afdef1 160000
--- a/openfpm_numerics
+++ b/openfpm_numerics
@@ -1 +1 @@
-Subproject commit dfb8b8f7b56e455c3e1db0ff39e5d513055eb472
+Subproject commit 7be6531420de26986cd584e3bec0f5ee51afdef1
diff --git a/src/Decomposition/CartDecomposition.hpp b/src/Decomposition/CartDecomposition.hpp
index 9981a969923af8477797b8a7f76243757172a805..dd7736e6612cf340e63636c8696f824b6605fc01 100644
--- a/src/Decomposition/CartDecomposition.hpp
+++ b/src/Decomposition/CartDecomposition.hpp
@@ -41,10 +41,10 @@
* \tparam Domain Structure that contain the information of your physical domain
*
* Given an N-dimensional space, this class decompose the space into a Cartesian grid of small
- * sub-sub-domain. At each sub-sub-domain is assigned an id that identify which processor is
- * going to take care of that part of space (in general the space assigned to a processor is
- * simply connected), a second step merge several sub-sub-domain with same id into bigger region
- * sub-domain with the id. Each sub-domain has an extended space called ghost part
+ * sub-sub-domain. To each sub-sub-domain is assigned an id that identify at which processor is
+ * assigned (in general the union of all the sub-sub-domain assigned to a processor is
+ * simply connected space), a second step merge several sub-sub-domain with same id into bigger region
+ * sub-domain. Each sub-domain has an extended space called ghost part
*
* Assuming that VCluster.getProcessUnitID(), equivalent to the MPI processor rank, return the processor local
* processor id, we define
@@ -468,18 +468,18 @@ public:
+----------------------------------------------------+
| |
| Processor 8 |
-| Sub-domain 0 +-----------------------------------+
+| Sub+domain 0 +-----------------------------------+
| | |
| | |
++--------------+---+---------------------------+----+ Processor 9 |
| | | B8_0 | | Subdomain 0 |
| +------------------------------------+ |
| | | | | |
- | | | XXXXXXXXXXXXX XX |B9_0| |
- | | B | X Processor 10 X | | |
- | Processor 5 | 5 | X Sub-domain 0 X | | |
- | Subdomain 0 | _ | X X +----------------------------------------+
- | | 0 | XXXXXXXXXXXXXXXX | | |
+ | | | |B9_0| |
+ | | B | Local processor | | |
+ | Processor 5 | 5 | Subdomain 0 | | |
+ | Subdomain 0 | _ | +----------------------------------------+
+ | | 0 | | | |
| | | | | |
| | | | | Processor 9 |
| | | |B9_1| Subdomain 1 |
@@ -490,37 +490,36 @@ public:
| |
+-----------------------------------+
+
\endverbatim
and also
G8_0 G9_0 G9_1 G5_0 (External ghost boxes)
- \verbatim
-
-+----------------------------------------------------+
-| |
-| Processor 8 |
-| Sub-domain 0 +-----------------------------------+
-| +---------------------------------------------+ |
-| | G8_0 | | |
-++--------------+------------------------------------+ | Processor 9 |
- | | | | | Subdomain 0 |
- | | | |G9_0| |
- | | | | | |
- | | | XXXXXXXXXXXXX XX | | |
- | | | X Processor 10 X | | |
- | Processor|5 | X Sub-domain 0 X | | |
- | Subdomain|0 | X X +-----------------------------------+
- | | | XXXXXXXXXXXXXXXX | | |
- | | G | | | |
- | | 5 | | | Processor 9 |
- | | | | | | Subdomain 1 |
- | | 0 | |G9_1| |
- | | | | | |
- | | | | | |
- +--------------+------------------------------------+ | |
- | | | |
- +----------------------------------------+----+------------------------------+
+ +----------------------------------------------------+
+ | Processor 8 |
+ | Subdomain 0 +-----------------------------------+
+ | | |
+ | +---------------------------------------------+ |
+ | | G8_0 | | |
++-----+---------------+------------------------------------+ | Processor 9 |
+| | | | | Subdomain 0 |
+| | | |G9_0| |
+| | | | | |
+| | | | | |
+| | | Local processor | | |
+| Processor 5 | | Sub+domain 0 | | |
+| Subdomain 0 | | +-----------------------------------+
+| | | | | |
+| | G | | | |
+| | 5 | | | Processor 9 |
+| | | | | | Subdomain 1 |
+| | 0 | |G9_1| |
+| | | | | |
+| | | | | |
++---------------------+------------------------------------+ | |
+ | | | |
+ +----------------------------------------+----+------------------------------+
\endverbatim
diff --git a/src/Grid/grid_dist_id.hpp b/src/Grid/grid_dist_id.hpp
index bd1e22a4b767a872ca5132888e92e4cd74cacf40..f7f436fe1c3dc788e40dd89fe042921a5511dac4 100644
--- a/src/Grid/grid_dist_id.hpp
+++ b/src/Grid/grid_dist_id.hpp
@@ -417,7 +417,33 @@ class grid_dist_id
{
}
- /*! \brief Initialize the Cell decomposer of the grid
+ void write_ie_boxes(std::string output)
+ {
+ // Write internal ghost box
+ VTKWriter<openfpm::vector<::Box<dim,size_t>>,VECTOR_BOX> vtk_box1;
+
+ openfpm::vector< openfpm::vector< ::Box<dim,size_t> > > boxes;
+
+ //! Carefully we have to ensure that boxes does not reallocate inside the for loop
+ boxes.reserve(ig_box.size());
+
+ //! Write internal ghost in grid units (Color encoded)
+ for (size_t p = 0 ; p < ig_box.size() ; p++)
+ {
+ boxes.add();
+
+ // Create a vector of boxes
+ for (size_t j = 0 ; j < ig_box.get(p).bid.size() ; j++)
+ {
+ boxes.last().add(ig_box.get(p).bid.get(j).box);
+ }
+
+ vtk_box1.add(boxes.last());
+ }
+ vtk_box1.write(output + std::string("internal_ghost_") + std::to_string(v_cl.getProcessUnitID()) + std::string(".vtk"));
+ }
+
+ /*! \brief Initialize the Cell decomposer of the grid
*
*
*/
@@ -476,6 +502,28 @@ class grid_dist_id
Create();
}
+protected:
+
+ /*! \brief Get the point where it start the origin of the grid in the sub-domain i
+ *
+ * \return the point
+ *
+ */
+ Point<dim,St> getOffset(size_t i)
+ {
+ return Point<dim,St>(gdb_ext.get(i).origin) * cd_sm.getCellBox().getP2();
+ }
+
+ /*! \brief Given a local sub-domain i with a local grid Domain + ghost return the part of the local grid that is domain
+ *
+ * \return the Box defining the domain in the local grid
+ *
+ */
+ Box<dim,size_t> getDomain(size_t i)
+ {
+ return gdb_ext.get(i).Dbox;
+ }
+
public:
// Which kind of grid the structure store
@@ -1143,6 +1191,16 @@ public:
}
}
+ /*! \brief Get the spacing on each dimension
+ *
+ * \param get the spacing
+ *
+ */
+ Point<dim,St> getSpacing()
+ {
+ return cd_sm.getCellBox().getP2();
+ }
+
/*! \brief Convert a g_dist_key_dx into a global key
*
* \see grid_dist_key_dx
@@ -1178,37 +1236,36 @@ public:
VTKWriter<boost::mpl::pair<device_grid,float>,VECTOR_GRIDS> vtk_g;
for (size_t i = 0 ; i < loc_grid.size() ; i++)
{
- Point<dim,St> offset = Point<dim,St>(gdb_ext.get(i).origin) * cd_sm.getCellBox().getP2();
+ Point<dim,St> offset = getOffset(i);
vtk_g.add(loc_grid.get(i),offset,cd_sm.getCellBox().getP2(),gdb_ext.get(i).Dbox);
}
vtk_g.write(output + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk");
- // Write internal ghost box
- VTKWriter<openfpm::vector<::Box<dim,size_t>>,VECTOR_BOX> vtk_box1;
-
- openfpm::vector< openfpm::vector< ::Box<dim,size_t> > > boxes;
-
- //! Carefully we have to ensure that boxes does not reallocate inside the for loop
- boxes.reserve(ig_box.size());
-
- //! Write internal ghost in grid units (Color encoded)
- for (size_t p = 0 ; p < ig_box.size() ; p++)
- {
- boxes.add();
+ write_ie_boxes(output);
- // Create a vector of boxes
- for (size_t j = 0 ; j < ig_box.get(p).bid.size() ; j++)
- {
- boxes.last().add(ig_box.get(p).bid.get(j).box);
- }
-
- vtk_box1.add(boxes.last());
- }
- vtk_box1.write(output + std::string("internal_ghost_") + std::to_string(v_cl.getProcessUnitID()) + std::string(".vtk"));
+ return true;
+ }
- vtk_g.write("vtk_grids.vtk");
+ /*! \brief Get the i sub-domain grid
+ *
+ * \param i sub-domain
+ *
+ * \return local grid
+ *
+ */
+ device_grid & get_loc_grid(size_t i)
+ {
+ return loc_grid.get(i);
+ }
- return true;
+ /*! \brief Return the number of local grid
+ *
+ * \return the number of local grid
+ *
+ */
+ size_t getN_loc_grid()
+ {
+ return loc_grid.size();
}
};
diff --git a/src/Grid/grid_dist_id_iterator.hpp b/src/Grid/grid_dist_id_iterator.hpp
index 56ba2ea60a56521d979f984a0f267683307c966d..e63cb25e77ebe39a0f6263099b743e7a2094a37d 100644
--- a/src/Grid/grid_dist_id_iterator.hpp
+++ b/src/Grid/grid_dist_id_iterator.hpp
@@ -85,7 +85,7 @@ class grid_dist_iterator<dim,device_grid,FREE>
size_t g_c;
//! List of the grids we are going to iterate
- Vcluster_object_array<device_grid> & gList;
+ const Vcluster_object_array<device_grid> & gList;
//! Extension of each grid: domain and ghost + domain
openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext;
@@ -243,7 +243,7 @@ class grid_dist_iterator<dim,device_grid,FIXED>
size_t g_c;
//! List of the grids we are going to iterate
- Vcluster_object_array<device_grid> & gList;
+ const Vcluster_object_array<device_grid> & gList;
//! Extension of each grid: domain and ghost + domain
const openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext;
diff --git a/src/Grid/grid_dist_key.hpp b/src/Grid/grid_dist_key.hpp
index 832a7d8338745cff2d7b6df65cee8019bd060e66..7d69d5b20551fcc6c8a6c309dd87cf15c0d615f1 100644
--- a/src/Grid/grid_dist_key.hpp
+++ b/src/Grid/grid_dist_key.hpp
@@ -40,6 +40,7 @@ public:
return key;
}
+
/*! \brief Get the reference key
*
* \return the local key
@@ -83,10 +84,39 @@ public:
return grid_dist_key_dx<dim>(getSub(),key);
}
+ /*! \brief Create a new key moving the old one
+ *
+ * \param c where to move for each component
+ *
+ * \return new key
+ *
+ */
+ inline grid_dist_key_dx<dim> move(const comb<dim> & c)
+ {
+ grid_key_dx<dim> key = getKey();
+ for (size_t i = 0 ; i < dim ; i++)
+ key.set_d(i,key.get(i) + c[i]);
+ return grid_dist_key_dx<dim>(getSub(),key);
+ }
+
inline grid_dist_key_dx(int g_c, const grid_key_dx<dim> & key)
:g_c(g_c),key(key)
{
}
+
+ std::string to_string()
+ {
+ std::stringstream str;
+
+ str << "sub_domain=" << g_c << " ";
+
+ for (size_t i = 0 ; i < dim ; i++)
+ str << "x[" << i << "]=" << key.get(i) << " ";
+
+ str << "\n";
+
+ return str.str();
+ }
};
#endif
diff --git a/src/Grid/staggered_dist_grid.hpp b/src/Grid/staggered_dist_grid.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a1bedc4127b79d914c3f9514ea045c97639928a9
--- /dev/null
+++ b/src/Grid/staggered_dist_grid.hpp
@@ -0,0 +1,127 @@
+/*
+ * staggered_grid.hpp
+ *
+ * Created on: Aug 19, 2015
+ * Author: i-bird
+ */
+
+#ifndef SRC_GRID_STAGGERED_DIST_GRID_HPP_
+#define SRC_GRID_STAGGERED_DIST_GRID_HPP_
+
+#include "Grid/grid_dist_id.hpp"
+#include "staggered_dist_grid_util.hpp"
+#include "VTKWriter.hpp"
+
+
+/*! \brief Implementation of the staggered grid
+ *
+ * \param dim Dimensionality of the staggered grid
+ * \param ele elements object on each dimensional objects, must be a stag_elements
+ *
+ * \verbatim
+
+ +--#--+--#--+--#--+--#--+--#--+--#--+
+ | | | | | | |
+ # * # * # * # * # * # * #
+ | | | | | | |
+ +--#--+--#--+--#--+--#--+--#--+--#--+
+ | | | | | | |
+ # * # * # * # * # * # * #
+ | | | | | | |
+ +--#--+--#--+--#--+--#--+--#--+--#--+
+ | | | | | | |
+ # * # * # * # * # * # * #
+ | | | | | | |
+ +--#--+--#--+--#--+--#--+--#--+--#--+
+ | | | | | | |
+ # * # * # * # * # * # * #
+ | | | | | | |
+ +--#--+--#--+--#--+--#--+--#--+--#--+
+ | | | | | | |
+ # * # * # * # * # * # * #
+ | | | | | | |
+ +--#--+--#--+--#--+--#--+--#--+--#--+
+
+\endverbatim
+
+ In the case of a 2D staggered grid we have 3 (in general dim+1 ) elements
+
+ + = vertex
+ # = edge
+ * = volume
+
+ ele = stag_ele<scalar<float>,Point_test<float>,scalar<float>>
+
+ It place a scalar on (*) an object Point_test<float> on (#) and an object scalar<float> on (+)
+
+ *
+ *
+ *
+ */
+template<unsigned int dim, typename St, typename T, typename Decomposition,typename Memory=HeapMemory , typename device_grid=grid_cpu<dim,T>>
+class staggered_grid_dist : public grid_dist_id<dim,St,T,Decomposition,Memory,device_grid>
+{
+
+public:
+
+ typedef T value_type;
+
+ 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)
+ {}
+
+ openfpm::vector<comb<dim>> c_prp[T::max_prop];
+
+ /*! \brief Set the staggered positions
+ *
+ *
+ */
+ template<unsigned int p> void setStagPosition(openfpm::vector<comb<dim>> & cmb)
+ {
+#ifdef SE_CLASS1
+ if (extends< typename boost::mpl::at<typename T::type,boost::mpl::int_<p> >::type >::mul() != cmb.size())
+ std::cerr << __FILE__ << ":" << __LINE__ << " error properties has " << extends< typename boost::mpl::at<typename T::type,boost::mpl::int_<p> >::type >::mul() << " components, but " << cmb.size() << "has been defined \n";
+#endif
+ c_prp.get(p) = cmb;
+ }
+
+ /*! \brief It set all the properties defined to be staggered on the default location
+ *
+ */
+ void setDefaultStagPosition()
+ {
+ // for each properties
+
+ stag_set_position<dim,typename T::type> ssp(c_prp);
+
+ boost::mpl::for_each_ref< boost::mpl::range_c<int,0,T::max_prop> >(ssp);
+ }
+
+ /*! \brief Write a vtk file with the information of the staggered grid
+ *
+ * \param str vtk output file
+ *
+ */
+ void write(std::string str)
+ {
+ stag_create_and_add_grid<dim,staggered_grid_dist<dim,St,T,Decomposition,Memory,device_grid>,St> sgw(*this, this->getVC().getProcessUnitID());
+
+ boost::mpl::for_each_ref< boost::mpl::range_c<int,0,T::max_prop> >(sgw);
+ }
+
+ /*! \brief Return if the properties is a staggered property or not
+ *
+ * \param prp property to check
+ *
+ * \return true if the property is staggered
+ *
+ */
+ bool is_staggered_prop(size_t prp)
+ {
+ return c_prp[prp].size() != 0;
+ }
+
+ friend class stag_create_and_add_grid<dim,staggered_grid_dist<dim,St,T,Decomposition,Memory,device_grid>,St>;
+};
+
+#endif /* SRC_GRID_STAGGERED_DIST_GRID_HPP_ */
diff --git a/src/Grid/staggered_dist_grid_util.hpp b/src/Grid/staggered_dist_grid_util.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..fc9cf783eb8588fa13eb3578907924be5bc3ac77
--- /dev/null
+++ b/src/Grid/staggered_dist_grid_util.hpp
@@ -0,0 +1,661 @@
+/*
+ * staggered_util.hpp
+ *
+ * Created on: Aug 19, 2015
+ * Author: i-bird
+ */
+
+#ifndef SRC_GRID_STAGGERED_DIST_GRID_UTIL_HPP_
+#define SRC_GRID_STAGGERED_DIST_GRID_UTIL_HPP_
+
+#include "util/common.hpp"
+#include "VTKWriter.hpp"
+#include "util/convert.hpp"
+
+
+/*! \brief write a property that has attributes
+ *
+ * \tparam ele object we are writing
+ * \tparam vtk vtk writer
+ * \tparam true in case the basic object has attributes
+ *
+ */
+template<typename ele, typename vtk, bool has_attributes=has_attributes<ele>::value>
+struct vtk_write
+{
+ /*! \brief Add the grid with attributes name
+ *
+ * \param vtk_w VTK writer
+ * \param output where to write
+ * \param i property to write
+ *
+ */
+ vtk_write(vtk vtk_w, const std::string output, const size_t i)
+ {
+ vtk_w.write(output + "_" + ele::attributes::name[i] + ".vtk",ele::attributes::name[i]);
+ }
+};
+
+/*! \brief Add to the vtk writer the key
+ *
+ * \tparam ele object we are writing
+ * \tparam vtk vtk writer
+ * \tparam false in case the basic object has not attributes
+ *
+ */
+template<typename ele, typename vtk>
+struct vtk_write<ele,vtk,false>
+{
+ /*! \brief Add the grid with attributes name
+ *
+ * \param vtk_w VTK writer
+ * \param output where to write
+ * \param i property to write
+ *
+ */
+ vtk_write(vtk vtk_w, const std::string output, const size_t i)
+ {
+ vtk_w.write(output + "_" + std::to_string(i) + ".vtk","attr" + std::to_string(i));
+ }
+};
+
+
+/*! \brief Classes to get the number of components of the properties
+ *
+ */
+template<typename T>
+struct extends
+{
+ static inline size_t mul()
+ {
+ return 1;
+ }
+
+ static inline size_t dim()
+ {
+ return 0;
+ }
+};
+
+//! Partial specialization for N=1 1D-Array
+template<typename T,size_t N1>
+struct extends<T[N1]>
+{
+ static inline size_t mul()
+ {
+ return N1;
+ }
+
+ static inline size_t dim()
+ {
+ return 1;
+ }
+};
+
+//! Partial specialization for N=2 2D-Array
+template<typename T,size_t N1,size_t N2>
+struct extends<T[N1][N2]>
+{
+ static inline size_t mul()
+ {
+ return N1 * N2;
+ }
+
+ static inline size_t dim()
+ {
+ return 2;
+ }
+};
+
+//! Partial specialization for N=3
+template<typename T,size_t N1,size_t N2,size_t N3>
+struct extends<T[N1][N2][N3]>
+{
+ static inline size_t mul()
+ {
+ return N1 * N2 * N3;
+ }
+
+ static inline size_t dim()
+ {
+ return 3;
+ }
+};
+
+//! Partial specialization for N=4
+template<typename T,size_t N1,size_t N2,size_t N3,size_t N4>
+struct extends<T[N1][N2][N3][N4]>
+{
+ static inline size_t mul()
+ {
+ return N1 * N2 * N3 * N4;
+ }
+
+ static inline size_t dim()
+ {
+ return 4;
+ }
+};
+
+//! Partial specialization for N=5
+template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5>
+struct extends<T[N1][N2][N3][N4][N5]>
+{
+ static inline size_t mul()
+ {
+ return N1 * N2 * N3 * N4 * N5;
+ }
+
+ static inline size_t dim()
+ {
+ return 5;
+ }
+};
+
+//! Partial specialization for N=6
+template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6>
+struct extends<T[N1][N2][N3][N4][N5][N6]>
+{
+ static inline size_t mul()
+ {
+ return N1 * N2 * N3 * N4 * N5 * N6;
+ }
+
+ static inline size_t dim()
+ {
+ return 6;
+ }
+};
+
+//! Partial specialization for N=7
+template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6, size_t N7>
+struct extends<T[N1][N2][N3][N4][N5][N6][N7]>
+{
+ static inline size_t mul()
+ {
+ return N1 * N2 * N3 * N4 * N5 * N6 * N7;
+ }
+
+ static inline size_t dim()
+ {
+ return 7;
+ }
+};
+
+//! Partial specialization for N=8
+template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6, size_t N7, size_t N8>
+struct extends<T[N1][N2][N3][N4][N5][N6][N7][N8]>
+{
+ static inline size_t mul()
+ {
+ return N1 * N2 * N3 * N4 * N5 * N6 * N7 * N8;
+ }
+
+ static inline size_t dim()
+ {
+ return 8;
+ }
+};
+
+//! Partial specialization for N=9
+template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6, size_t N7, size_t N8, size_t N9>
+struct extends<T[N1][N2][N3][N4][N5][N6][N7][N8][N9]>
+{
+ static inline size_t mul()
+ {
+ return N1 * N2 * N3 * N4 * N5 * N6 * N7 * N8 * N9;
+ }
+
+ static inline size_t dim()
+ {
+ return 9;
+ }
+};
+
+//! Partial specialization for N=10
+template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6, size_t N7, size_t N8, size_t N9, size_t N10>
+struct extends<T[N1][N2][N3][N4][N5][N6][N7][N8][N9][N10]>
+{
+ static inline size_t mul()
+ {
+ return N1 * N2 * N3 * N4 * N5 * N6 * N7 * N8 * N9 * N10;
+ }
+
+ static inline size_t dim()
+ {
+ return 10;
+ }
+};
+
+///////////////////// Copy grid extends
+
+/*! \brief Classes to copy each component into a grid and add to the VTKWriter the grid
+ *
+ * \param T property to write
+ * \param dim dimansionality
+ * \param St type of space
+ * \param VTKW VTK writer
+ * \param
+ *
+ */
+template<typename T>
+struct write_stag
+{
+ /*! \brieg write the staggered grid
+ *
+ * \tparam p_val property we are going to write
+ * \tparam sg staggered grid type
+ * \tparam v_g vector of grids
+ *
+ * \param st_g staggered grid
+ * \param v_g vector of grids
+ * \param i local grid of the staggered grid we are writing
+ *
+ */
+ template<unsigned int p_val, typename sg, typename v_g> static inline void write(sg & st_g, v_g & vg,size_t lg)
+ {
+ // Add a grid;
+ vg.add();
+ size_t k = vg.size() - 1;
+
+ // Get the source and destination grid
+ auto & g_src = st_g.get_loc_grid(lg);
+ auto & g_dst = vg.get(k);
+
+ // Set dimensions and memory
+ g_dst.resize(g_src.getGrid().getSize());
+
+ // copy
+
+ auto it = vg.get(k).getIterator();
+
+ while(it.isNext())
+ {
+ g_dst.template get<0>(it.get()) = g_src.template get<p_val>(it.get());
+
+ ++it;
+ }
+ }
+};
+
+/*! \brief for each component add a grid fill it, and add to the VTK writer
+ *
+ * \param T Property to copy
+ * \param N1 number of components
+ *
+ */
+template<typename T,size_t N1>
+struct write_stag<T[N1]>
+{
+ /*! \brieg write the staggered grid
+ *
+ * \tparam p_val property we are going to write
+ * \tparam sg staggered grid type
+ * \tparam v_g vector of grids
+ *
+ * \param st_g staggered grid
+ * \param v_g vector of grids
+ * \param i local grid of the staggered grid we are writing
+ *
+ */
+ template<unsigned int p_val, typename sg, typename v_g> static inline void write(sg & st_g, v_g & vg,size_t lg)
+ {
+ for (size_t i = 0 ; i < N1 ; i++)
+ {
+ // Add a grid;
+ vg.add();
+ size_t k = vg.size() - 1;
+
+ // Get the source and destination grid
+ auto & g_src = st_g.get_loc_grid(lg);
+ auto & g_dst = vg.get(k);
+
+ // Set dimensions and memory
+ g_dst.resize(g_src.getGrid().getSize());
+
+ auto it = vg.get(k).getIterator();
+
+ while(it.isNext())
+ {
+ g_dst.template get<0>(it.get()) = g_src.template get<p_val>(it.get())[i];
+
+ ++it;
+ }
+ }
+ }
+};
+
+//! Partial specialization for N=2 2D-Array
+template<typename T,size_t N1,size_t N2>
+struct write_stag<T[N1][N2]>
+{
+ /*! \brieg write the staggered grid
+ *
+ * \tparam p_val property we are going to write
+ * \tparam sg staggered grid type
+ * \tparam v_g vector of grids
+ *
+ * \param st_g staggered grid
+ * \param v_g vector of grids
+ * \param i local grid of the staggered grid we are writing
+ *
+ */
+ template<unsigned int p_val, typename sg, typename v_g> static inline void write(sg & st_g, v_g & vg,size_t lg)
+ {
+ for (size_t i = 0 ; i < N1 ; i++)
+ {
+ for (size_t j = 0 ; j < N2 ; j++)
+ {
+ // Add a grid;
+ vg.add();
+ size_t k = vg.size() - 1;
+
+ // Set dimensions and memory
+ vg.get(k).resize(st_g.get_loc_grid(lg).getGrid().getSize());
+
+ // copy
+ auto & g_src = st_g.get_loc_grid(lg);
+ auto & g_dst = vg.get(k);
+ auto it = vg.get(k).getIterator();
+
+ while(it.isNext())
+ {
+ g_dst.template get<0>(it.get()) = g_src.template get<p_val>(it.get())[i][j];
+
+ ++it;
+ }
+ }
+ }
+ }
+};
+
+///////////////////// Staggered default positioning ////////////////////////
+
+/*! \brief this class is a functor for "for_each" algorithm
+ *
+ * For each element of the boost::vector the operator() is called.
+ * Is mainly used to produce a default position vector for each
+ * property
+ *
+ * \tparam dim dimensionality
+ * \tparam v boost::fusion::vector of properties
+ * \tparam has_posMask case when v has a position mask
+ *
+ */
+
+template<unsigned int dim, typename v, bool has_pM = has_posMask<v>::value>
+class stag_set_position
+{
+ openfpm::vector<comb<dim>> (& pos_prp)[boost::fusion::result_of::size<v>::type::value];
+
+public:
+
+ stag_set_position( openfpm::vector<comb<dim>> (& pos_prp)[boost::fusion::result_of::size<v>::type::value])
+ :pos_prp(pos_prp)
+ {}
+
+ //! It call the copy function for each property
+ template<typename T>
+ void operator()(T& t) const
+ {
+ // This is the type of the object we have to copy
+ typedef typename boost::mpl::at<v,typename boost::mpl::int_<T::value>>::type prop;
+
+ bool val = prop::stag_pos_mask[T::value];
+
+ if (val == false)
+ return;
+
+ // Dimension of the object
+ size_t dim_prp = extends<prop>::dim();
+
+ // It is a scalar
+ if (dim_prp == 0)
+ {
+ comb<dim> c;
+ c.zero();
+
+ // It stay in the center
+ pos_prp[T::value].add(c);
+ }
+ else if (dim_prp == 1)
+ {
+ // It stay on the object of dimension dim-1 (Negative part)
+ for (size_t i = 0 ; i < dim ; i++)
+ {
+ comb<dim> c;
+ c.zero();
+ c.value(i) = -1;
+
+ pos_prp[T::value].add(c);
+ }
+ }
+ else if (dim_prp == 2)
+ {
+ // Create an hypercube object
+ HyperCube<dim> hyp;
+
+ // Diagonal part live in
+ for (size_t i = 0 ; i < dim ; i++)
+ {
+ comb<dim> c1 = pos_prp[T::value-1].get(i);
+ for (size_t j = 0 ; j < dim ; j++)
+ {
+ comb<dim> c2;
+ c2.zero();
+ c2.value(i) = 1;
+
+ comb<dim> c_res = (c1 + c2) & 0x1;
+
+ pos_prp[T::value].add(c_res);
+ }
+ }
+ }
+ else if (dim_prp > 2)
+ {
+ std::cerr << __FILE__ << ":" << __LINE__ << " Tensor of rank bigger than 2 are not supported";
+ }
+ }
+};
+
+///////////////////// Staggered default positioning ////////////////////////
+
+/*! \brief this class is a functor for "for_each" algorithm
+ *
+ * For each element of the boost::vector the operator() is called.
+ * Is mainly used to produce a default position vector for each
+ * property
+ *
+ * \tparam vector of properties
+ *
+ */
+
+template<unsigned int dim, typename v>
+class stag_set_position<dim,v,false>
+{
+private:
+ openfpm::vector<comb<dim>> (& pos_prp)[boost::fusion::result_of::size<v>::type::value];
+
+
+public:
+ stag_set_position( openfpm::vector<comb<dim>> (& pos_prp)[boost::fusion::result_of::size<v>::type::value])
+ :pos_prp(pos_prp)
+ {}
+
+ //! It call the copy function for each property
+ template<typename T>
+ void operator()(T& t) const
+ {
+ // This is the type of the object we have to copy
+ typedef typename boost::mpl::at<v,typename boost::mpl::int_<T::value>>::type prop;
+
+ // Dimension of the object
+ size_t dim_prp = extends<prop>::dim();
+
+ // It is a scalar
+ if (dim_prp == 0)
+ {
+ comb<dim> c;
+ c.zero();
+
+ // It stay in the center
+ pos_prp[T::value].add(c);
+ }
+ else if (dim_prp == 1)
+ {
+ // It stay on the object of dimension dim-1 (Negative part)
+ for (size_t i = 0 ; i < dim ; i++)
+ {
+ comb<dim> c;
+ c.zero();
+ c.getComb()[i] = -1;
+
+ pos_prp[T::value].add(c);
+ }
+ }
+ else if (dim_prp == 2)
+ {
+ // Diagonal part live in
+ for (size_t i = 0 ; i < dim ; i++)
+ {
+ comb<dim> c1 = pos_prp[T::value-1].get(i);
+ for (size_t j = 0 ; j < dim ; j++)
+ {
+ comb<dim> c2;
+ c2.zero();
+ c2.getComb()[j] = 1;
+
+ comb<dim> c_res = (c2 + c1).flip();
+
+ pos_prp[T::value].add(c_res);
+ }
+ }
+ }
+ else if (dim_prp > 2)
+ {
+ std::cerr << __FILE__ << ":" << __LINE__ << " Tensor of rank bigger than 2 are not supported";
+ }
+ }
+};
+
+/*! \brief It create separated grid for each properties to write them into a file
+ *
+ * \tparam dim dimensionality of the grids
+ * \tparam obj type object to print, must be in OpenFPM format
+ *
+ */
+template<unsigned int dim, typename st_grid, typename St>
+class stag_create_and_add_grid
+{
+
+ size_t p_id;
+
+ // staggered grid to write
+ st_grid & st_g;
+
+public:
+
+ /*! \brief Constructor
+ *
+ * \param st_g staggered grid
+ * \param p_id process id
+ *
+ */
+ stag_create_and_add_grid(st_grid & st_g, size_t p_id)
+ :p_id(p_id),st_g(st_g)
+ {}
+
+ template<unsigned int p_val> void out_normal()
+ {
+ // property type
+ typedef typename boost::mpl::at< typename st_grid::value_type::type , typename boost::mpl::int_<p_val> >::type ele;
+
+ // create an openfpm format object from the property type
+ typedef object<typename boost::fusion::vector<ele>> d_object;
+
+ VTKWriter<boost::mpl::pair<grid_cpu<dim, d_object >,St>,VECTOR_GRIDS> vtk_w;
+
+ // Create a vector of grids
+
+ openfpm::vector< grid_cpu<dim, d_object > > vg(st_g.getN_loc_grid());
+
+ // for each domain grid
+ for (size_t i = 0 ; i < vg.size() ; i++)
+ {
+ // Set dimansions and memory
+ vg.get(i).resize(st_g.get_loc_grid(i).getGrid().getSize());
+
+ auto & g_src = st_g.get_loc_grid(i);
+ auto & g_dst = vg.get(i);
+
+ auto it = vg.get(i).getIterator();
+
+ while(it.isNext())
+ {
+ object_si_d< decltype(g_src.get_o(it.get())),decltype(g_dst.get_o(it.get())) ,OBJ_ENCAP,p_val>(g_src.get_o(it.get()),g_dst.get_o(it.get()));
+
+ ++it;
+ }
+
+ Point<dim,St> offset = st_g.getOffset(i);
+ Point<dim,St> spacing = st_g.getSpacing();
+ Box<dim,size_t> dom = st_g.getDomain(i);
+
+ vtk_w.add(g_dst,offset,spacing,dom);
+ }
+
+ vtk_w.write("vtk_grids_st_" + std::to_string(p_id) + "_" + std::to_string(p_val) + ".vtk");
+ }
+
+ template<unsigned int p_val> void out_staggered()
+ {
+ // property type
+ typedef typename boost::mpl::at< typename st_grid::value_type::type , typename boost::mpl::int_<p_val> >::type ele;
+
+ // Eliminate the extends
+ typedef typename std::remove_all_extents<ele>::type r_ele;
+
+ // create an openfpm format object from the property type
+ typedef object<typename boost::fusion::vector<r_ele>> d_object;
+
+ VTKWriter<boost::mpl::pair<grid_cpu<dim, d_object >,St>,VECTOR_ST_GRIDS> vtk_w;
+
+ // Create a vector of grids
+ openfpm::vector< grid_cpu<dim, d_object > > vg;
+ vg.reserve(st_g.getN_loc_grid() * extends<ele>::mul());
+
+ size_t k = 0;
+
+ // for each domain grid
+ for (size_t i = 0 ; i < st_g.getN_loc_grid() ; i++)
+ {
+ write_stag<ele>::template write<p_val, st_grid,openfpm::vector< grid_cpu<dim, d_object > > >(st_g,vg,i);
+
+ // for each component
+ for ( ; k < vg.size() ; k++)
+ {
+ Point<dim,St> offset = st_g.getOffset(i);
+ Point<dim,St> spacing = st_g.getSpacing();
+ Box<dim,size_t> dom = st_g.getDomain(i);
+
+ vtk_w.add(i,vg.get(k),offset,spacing,dom,st_g.c_prp[p_val].get(k));
+ }
+
+ k = vg.size();
+ }
+
+ vtk_write<typename st_grid::value_type,VTKWriter<boost::mpl::pair<grid_cpu<dim, d_object >,St>,VECTOR_ST_GRIDS>> v(vtk_w,"vtk_grids_st_" + std::to_string(p_id),p_val);
+ }
+
+ //! It call the copy function for each property
+ template<typename T>
+ void operator()(T& t)
+ {
+ if (st_g.is_staggered_prop(T::value) == false)
+ out_normal<T::value>();
+ else
+ out_staggered<T::value>();
+ }
+};
+
+#endif /* SRC_GRID_STAGGERED_DIST_GRID_UTIL_HPP_ */
diff --git a/src/Grid/staggered_grid_dist_unit_test.hpp b/src/Grid/staggered_grid_dist_unit_test.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d7fcaeb035f28d13d9790e50bbdab41a98b1a455
--- /dev/null
+++ b/src/Grid/staggered_grid_dist_unit_test.hpp
@@ -0,0 +1,71 @@
+/*
+ * staggered_grid_unit_test.hpp
+ *
+ * Created on: Aug 20, 2015
+ * Author: i-bird
+ */
+
+#ifndef SRC_GRID_STAGGERED_GRID_DIST_UNIT_TEST_HPP_
+#define SRC_GRID_STAGGERED_GRID_DIST_UNIT_TEST_HPP_
+
+#include "staggered_dist_grid.hpp"
+#include "Point_test.hpp"
+
+BOOST_AUTO_TEST_SUITE( staggered_grid_dist_id_test )
+
+
+BOOST_AUTO_TEST_CASE( staggered_grid_dist_unit_test)
+{
+ typedef Point2D_test<float> p;
+
+ Vcluster & v_cl = *global_v_cluster;
+
+ // Domain
+ Box<2,float> domain({0.0,0.0},{1.0,1.0});
+
+ size_t k = 1024;
+
+/* for (size_t k = 1024 ; k >= 2 ; k--)
+ {*/
+ BOOST_TEST_CHECKPOINT( "Testing grid k=" << k );
+
+ // grid size
+ size_t sz[2];
+ sz[0] = k;
+ sz[1] = k;
+
+ // Ghost
+ Ghost<2,float> g(0.01);
+
+ staggered_grid_dist<2,float,Point2D_test<float>,CartDecomposition<2,float>> sg(sz,domain,g);
+
+ sg.setDefaultStagPosition();
+
+ auto it = sg.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ grid_key_dx<2> keyg = sg.getGKey(key);
+
+ sg.template get<p::s>(key) = 1;
+
+ sg.template get<p::v>(key)[0] = 0;
+ sg.template get<p::v>(key)[1] = 1;
+
+ sg.template get<p::t>(key)[0][0] = 0;
+ sg.template get<p::t>(key)[0][1] = 1;
+ sg.template get<p::t>(key)[1][0] = 2;
+ sg.template get<p::t>(key)[1][1] = 3;
+
+ ++it;
+ }
+
+ sg.write("stag_test.vtk");
+/* }*/
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+#endif /* SRC_GRID_STAGGERED_GRID_DIST_UNIT_TEST_HPP_ */
diff --git a/src/main.cpp b/src/main.cpp
index 529be072e4382f17919893b0069638912b2161e9..4e949b15552d7c5fce5a38bd1d323df62ed4fc1c 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -18,6 +18,7 @@
#include "util.hpp"
#include "unit_test_init_cleanup.hpp"
+#include "Grid/staggered_grid_dist_unit_test.hpp"
#include "Decomposition/CartDecomposition_unit_test.hpp"
#include "Decomposition/ORB_unit_test.hpp"
#include "Graph/CartesianGraphFactory_unit_test.hpp"