diff --git a/src/Grid/grid_dist_id.hpp b/src/Grid/grid_dist_id.hpp
index c16ae0f8fdf5843aabdeee22949c42995a31b2a1..208b654f9eb2edca1d9607b96ff49db1be20770d 100644
--- a/src/Grid/grid_dist_id.hpp
+++ b/src/Grid/grid_dist_id.hpp
@@ -200,7 +200,6 @@ class grid_dist_id
{
// Get the internal ghost boxes and transform into grid units
::Box<dim,St> ib_dom = dec.getProcessorIGhostBox(i,j);
- ib_dom -= cd_sm.getOrig();
::Box<dim,long int> ib = cd_sm.convertDomainSpaceIntoGridUnits(ib_dom,dec.periodicity());
// Check if ib is valid if not it mean that the internal ghost does not contain information so skip it
@@ -321,7 +320,6 @@ class grid_dist_id
{
// Get the internal ghost boxes and transform into grid units
::Box<dim,St> ib_dom = dec.getLocalIGhostBox(i,j);
- ib_dom -= cd_sm.getOrig();
::Box<dim,long int> ib = cd_sm.convertDomainSpaceIntoGridUnits(ib_dom,dec.periodicity());
// Check if ib is valid if not it mean that the internal ghost does not contain information so skip it
diff --git a/src/Grid/grid_dist_id_unit_test.cpp b/src/Grid/grid_dist_id_unit_test.cpp
index 22cd749ccc7e22b7caf1c72ab2cc6554928f95c3..ff365154648a7188f1dded515cc5c3d16263992c 100644
--- a/src/Grid/grid_dist_id_unit_test.cpp
+++ b/src/Grid/grid_dist_id_unit_test.cpp
@@ -78,7 +78,7 @@ BOOST_AUTO_TEST_CASE( grid_dist_id_domain_grid_unit_converter3D_test)
{
// Get the local hyper-cube
SpaceBox<3,float> sub = dec.getSubDomain(i);
- sub -= domain.getP1();
+// sub -= domain.getP1();
Box<3,size_t> g_box = g_dist.getCellDecomposer().convertDomainSpaceIntoGridUnits(sub,bc);
diff --git a/src/Grid/grid_dist_util.hpp b/src/Grid/grid_dist_util.hpp
index bba37eb0335b80bd7e99aeb03f9be68a185861d4..d16bf8bb1d073e85edd34b1a9d1778f61b823ea6 100644
--- a/src/Grid/grid_dist_util.hpp
+++ b/src/Grid/grid_dist_util.hpp
@@ -78,9 +78,7 @@ template<int dim, typename Decomposition> inline void create_gdb_ext(openfpm::ve
// Get the local sub-domain (Grid conversion must be done with the domain P1 equivalent to 0.0)
// consider that the sub-domain with point P1 equivalent to the domain P1 is a (0,0,0) in grid unit
SpaceBox<Decomposition::dims, typename Decomposition::stype> sp = dec.getSubDomain(i);
- sp -= cd_sm.getOrig();
SpaceBox<Decomposition::dims, typename Decomposition::stype> sp_g = dec.getSubDomainWithGhost(i);
- sp_g -= cd_sm.getOrig();
// Convert from SpaceBox<dim,St> to SpaceBox<dim,long int>
SpaceBox<Decomposition::dims,long int> sp_t = cd_sm.convertDomainSpaceIntoGridUnits(sp,dec.periodicity());
diff --git a/src/Vector/vector_dist.hpp b/src/Vector/vector_dist.hpp
index d2efb1f119d52055d5d0f880182f89bee7afb52d..305e2d9024a79d95bed78057c6ac5228a5a80331 100644
--- a/src/Vector/vector_dist.hpp
+++ b/src/Vector/vector_dist.hpp
@@ -347,6 +347,48 @@ public:
return v_prp.template get<id>(g_m - 1);
}
+ /*! \brief Construct a cell list symmetric based on a cut of radius
+ *
+ * \tparam CellL CellList type to construct
+ *
+ * \param r_cut interation radius, or size of each cell
+ *
+ * \return the Cell list
+ *
+ */
+ template<typename CellL = CellList<dim, St, FAST, shift<dim, St> > > CellL getCellListSym(St r_cut)
+ {
+ // Cell list
+ CellL cell_list;
+ size_t div[dim];
+
+ // Calculate the Cell list division for this CellList
+ CellDecomposer<dim,St,shift<dim,St>> cd_sm;
+
+ for (size_t i = 0 ; i < dim ; i++)
+ div[i] = (domain.getHigh(i) - domain.getLow(i)) / r_cut;
+
+ size_t pad = 0;
+ Ghost g = getDecomposition().getGhost();
+ g.magnify(1.013);
+
+ // Calculate the maximum padding
+ for (size_t i = 0 ; i < dim ; i++)
+ {
+ size_t tmp = std::ceil(fabs(g.getLow(i)) / r_cut);
+ pad = (pad > tmp)?pad:tmp;
+ }
+
+ cd_sm.Initialize(domain,div,pad);
+
+ // get the processor bounding box
+ Box<dim, St> pbox = getDecomposition().getProcessorBounds();
+
+ cell_list.Initialize(cd_sm, pbox);
+
+ updateCellList(cell_list);
+ }
+
/*! \brief Construct a cell list starting from the stored particles
*
* \tparam CellL CellList type to construct
diff --git a/src/Vector/vector_dist_unit_test.hpp b/src/Vector/vector_dist_unit_test.hpp
index 88a7275a268faa0eebc47208390061bb33de8a76..ad982528eed45dd6924925ded3691977ff346ea7 100644
--- a/src/Vector/vector_dist_unit_test.hpp
+++ b/src/Vector/vector_dist_unit_test.hpp
@@ -1817,6 +1817,140 @@ BOOST_AUTO_TEST_CASE( vector_of_vector_dist )
BOOST_REQUIRE_EQUAL(cnt,4*4096ul);
}
+BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
+{
+ Vcluster & v_cl = create_vcluster();
+
+ if (v_cl.getProcessingUnits() > 24)
+ return;
+
+ // set the seed
+ // create the random generator engine
+ std::srand(v_cl.getProcessUnitID());
+ std::default_random_engine eg;
+ std::uniform_real_distribution<float> ud(0.0f, 1.0f);
+
+ long int k = 524288 * v_cl.getProcessingUnits();
+
+ long int big_step = k / 4;
+ big_step = (big_step == 0)?1:big_step;
+
+ print_test("Testing 3D periodic vector symmetric cell-list k=",k);
+ BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric cell-list k=" << k );
+
+ Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
+
+ // Boundary conditions
+ size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+
+ float r_cut = 0.1;
+
+ // ghost
+ Ghost<3,float> ghost(r_cut);
+
+ typedef aggregate<size_t,size_t> part_prop;
+
+ // Distributed vector
+ vector_dist<3,float, part_prop > vd(k,box,bc,ghost);
+
+ auto it = vd.getIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ vd.getPos(key)[0] = ud(eg);
+ vd.getPos(key)[1] = ud(eg);
+ vd.getPos(key)[2] = ud(eg);
+
+ // Fill some properties randomly
+
+ vd.getProp<0>(key) = 0.0;
+
+ ++it;
+ }
+
+ vd.map();
+
+ // sync the ghost
+ vd.ghost_get<0>();
+
+ auto NN = vd.getCellList(0.1);
+
+ auto p_it = vd.getDomainIterator();
+
+ while (p_it.isNext())
+ {
+ auto p = p_it.get();
+
+ Point<3,float> xp = vd.getPos(p);
+
+ auto Np = NN.getIterator(NN.getCell(vd.getPos(p)));
+
+ while (Np.isNext())
+ {
+ auto q = Np.get();
+
+ // repulsive
+
+ Point<3,float> xq = vd.getPos(q);
+ Point<3,float> f = (xp - xq);
+
+ float distance = f.norm();
+
+ // Particle should be inside 2 * r_cut range
+
+ if (distance < r_cut )
+ vd.getProp<0>(p)++;
+
+ ++Np;
+ }
+
+ ++p_it;
+ }
+
+ // We now try symmetric Cell-list
+
+ auto NN2 = vd.getCellListSym(0.1);
+
+ auto p_it2 = vd.getDomainIterator();
+
+ while (p_it2.isNext())
+ {
+ auto p = p_it2.get();
+
+ Point<3,float> xp = vd.getPos(p);
+
+ auto Np = NN2.getIterator(NN2.getCell(vd.getPos(p)));
+
+ while (Np.isNext())
+ {
+ auto q = Np.get();
+
+ // repulsive
+
+ Point<3,float> xq = vd.getPos(q);
+ Point<3,float> f = (xp - xq);
+
+ float distance = f.norm();
+
+ // Particle should be inside r_cut range
+
+ if (distance < r_cut )
+ {
+ vd.getProp<1>(p)++;
+ vd.getProp<1>(q)++;
+ }
+
+ ++Np;
+ }
+
+ ++p_it;
+ }
+
+ vd.ghost_put<add,1>();
+}
+
#include "vector_dist_cell_list_tests.hpp"
#include "vector_dist_NN_tests.hpp"
#include "vector_dist_complex_prp_unit_test.hpp"