From 1214b212de48035277e93f1745f462af58b305aa Mon Sep 17 00:00:00 2001
From: Pietro Incardona <i-bird@localhost.localdomain>
Date: Wed, 13 Jan 2016 20:12:36 -0500
Subject: [PATCH] Added CellList in vector distributed

---
 src/Decomposition/CartDecomposition.hpp |  12 ++
 src/Decomposition/ie_ghost.hpp          |   4 +-
 src/Decomposition/ie_loc_ghost.hpp      |  13 ++
 src/Grid/grid_dist_id.hpp               |   2 +-
 src/Grid/grid_dist_id_unit_test.hpp     |   4 +-
 src/Vector/vector_dist.hpp              |  36 ++++-
 src/Vector/vector_dist_unit_test.hpp    | 169 +++++++++++++++++++++++-
 7 files changed, 221 insertions(+), 19 deletions(-)

diff --git a/src/Decomposition/CartDecomposition.hpp b/src/Decomposition/CartDecomposition.hpp
index 64416216..898e17e1 100644
--- a/src/Decomposition/CartDecomposition.hpp
+++ b/src/Decomposition/CartDecomposition.hpp
@@ -682,6 +682,9 @@ p1[0]<-----+         +----> p2[0]
 		cart.ss_box = ss_box;
 		cart.ghost = g;
 
+		for (size_t i = 0 ; i < dim ; i++)
+			cart.bc[i] = bc[i];
+
 		(static_cast<nn_prcs<dim,T> &>(cart)).create(box_nn_processor, sub_domains);
 		(static_cast<nn_prcs<dim,T> &>(cart)).applyBC(domain,ghost,bc);
 
@@ -720,6 +723,9 @@ p1[0]<-----+         +----> p2[0]
 		cart.bbox = bbox;
 		cart.ss_box = ss_box;
 
+		for (size_t i = 0 ; i < dim ; i++)
+			cart.bc[i] = this->bc[i];
+
 		return cart;
 	}
 
@@ -750,6 +756,9 @@ p1[0]<-----+         +----> p2[0]
 		bbox = cart.bbox;
 		ss_box = cart.ss_box;
 
+		for (size_t i = 0 ; i < dim ; i++)
+			bc[i] = cart.bc[i];
+
 		return *this;
 	}
 
@@ -780,6 +789,9 @@ p1[0]<-----+         +----> p2[0]
 		cart.bbox = bbox;
 		cart.ss_box = ss_box;
 
+		for (size_t i = 0 ; i < dim ; i++)
+			cart.bc[i] = bc[i];
+
 		return *this;
 	}
 
diff --git a/src/Decomposition/ie_ghost.hpp b/src/Decomposition/ie_ghost.hpp
index 333e12e1..b865944a 100644
--- a/src/Decomposition/ie_ghost.hpp
+++ b/src/Decomposition/ie_ghost.hpp
@@ -951,7 +951,7 @@ public:
 		if (getNIGhostBox() != ig.getNIGhostBox())
 			return false;
 
-		for (size_t i = 0 ; i < getNIGhostBox() ; i++)
+		for (size_t i = 0 ; i < proc_int_box.size() ; i++)
 		{
 			if (getProcessorNIGhost(i) != ig.getProcessorNIGhost(i))
 				return false;
@@ -970,7 +970,7 @@ public:
 				return false;
 		}
 
-		for (size_t i = 0 ; i < getNEGhostBox() ; i++)
+		for (size_t i = 0 ; i < proc_int_box.size() ; i++)
 		{
 			if (getProcessorNEGhost(i) != ig.getProcessorNEGhost(i))
 				return false;
diff --git a/src/Decomposition/ie_loc_ghost.hpp b/src/Decomposition/ie_loc_ghost.hpp
index 1b3c53ad..2669fc8b 100644
--- a/src/Decomposition/ie_loc_ghost.hpp
+++ b/src/Decomposition/ie_loc_ghost.hpp
@@ -39,6 +39,9 @@ class ie_loc_ghost
 	 */
 	void create_loc_ghost_ebox(Ghost<dim,T> & ghost, openfpm::vector<SpaceBox<dim,T>> & sub_domains, openfpm::vector<Box_loc_sub<dim,T>> & sub_domains_prc)
 	{
+		comb<dim> zero;
+		zero.zero();
+
 		loc_ghost_box.resize(sub_domains.size());
 
 		// For each sub-domain
@@ -54,6 +57,9 @@ class ie_loc_ghost
 			{
 				size_t rj = sub_domains_prc.get(j).sub;
 
+				if (rj == i && sub_domains_prc.get(j).cmb == zero)
+					continue;
+
 				::Box<dim,T> bi;
 
 				bool intersect = sub_with_ghost.Intersect(sub_domains_prc.get(j).bx,bi);
@@ -88,6 +94,9 @@ class ie_loc_ghost
 	 */
 	void create_loc_ghost_ibox(Ghost<dim,T> & ghost, openfpm::vector<SpaceBox<dim,T>> & sub_domains, openfpm::vector<Box_loc_sub<dim,T>> & sub_domains_prc)
 	{
+		comb<dim> zero;
+		zero.zero();
+
 		loc_ghost_box.resize(sub_domains.size());
 
 		// For each sub-domain
@@ -98,6 +107,10 @@ class ie_loc_ghost
 			{
 				SpaceBox<dim,T> sub_with_ghost = sub_domains_prc.get(j).bx;
 				size_t rj = sub_domains_prc.get(j).sub;
+
+				if (rj == i && sub_domains_prc.get(j).cmb == zero)
+					continue;
+
 				// enlarge the sub-domain with the ghost
 				sub_with_ghost.enlarge(ghost);
 
diff --git a/src/Grid/grid_dist_id.hpp b/src/Grid/grid_dist_id.hpp
index ccd8eec4..0379c9a1 100644
--- a/src/Grid/grid_dist_id.hpp
+++ b/src/Grid/grid_dist_id.hpp
@@ -562,7 +562,7 @@ public:
 	 */
 	size_t size(size_t i) const
 	{
-		return ginfo_v.size();
+		return ginfo_v.size(i);
 	}
 
 	static inline Ghost<dim,float> convert_ghost(const Ghost<dim,long int> & gd,const CellDecomposer_sm<dim,St> & cd_sm)
diff --git a/src/Grid/grid_dist_id_unit_test.hpp b/src/Grid/grid_dist_id_unit_test.hpp
index d4f84fd8..ce59a1d9 100644
--- a/src/Grid/grid_dist_id_unit_test.hpp
+++ b/src/Grid/grid_dist_id_unit_test.hpp
@@ -271,8 +271,8 @@ void Test2D(const Box<2,float> & domain, long int k)
 		grid_key_dx<2> start = dom2.getStart();
 		grid_key_dx<2> stop = dom2.getStop();
 
-		BOOST_REQUIRE_EQUAL(stop.get(0),g_dist.size(0));
-		BOOST_REQUIRE_EQUAL(stop.get(1),g_dist.size(1));
+		BOOST_REQUIRE_EQUAL(stop.get(0),g_dist.size(0)-1);
+		BOOST_REQUIRE_EQUAL(stop.get(1),g_dist.size(1)-1);
 
 		BOOST_REQUIRE_EQUAL(start.get(0),0);
 		BOOST_REQUIRE_EQUAL(start.get(1),0);
diff --git a/src/Vector/vector_dist.hpp b/src/Vector/vector_dist.hpp
index 5cb96f0d..85e0b013 100644
--- a/src/Vector/vector_dist.hpp
+++ b/src/Vector/vector_dist.hpp
@@ -162,18 +162,25 @@ private:
 			size_t p_id = 0;
 
 			// Check if the particle is inside the domain
-//			if (dec.getDomain().isInside(v_pos.get(key)) == true)
+			if (dec.getDomain().isInside(v_pos.get(key)) == true)
 				p_id = dec.processorIDBC(v_pos.get(key));
-//			else
-//				p_id = obp::out(key,v_cl.getProcessUnitID());
+			else
+				p_id = obp::out(key,v_cl.getProcessUnitID());
 
 
 			// Particle to move
 			if (p_id != v_cl.getProcessUnitID())
 			{
-				prc_sz.get(p_id)++;
-				lbl_p.get(p_id).add(key);
-				opart.add(key);
+				if (p_id != -1)
+				{
+					prc_sz.get(p_id)++;
+					lbl_p.get(p_id).add(key);
+					opart.add(key);
+				}
+				else
+				{
+					opart.add(key);
+				}
 			}
 
 			// Add processors and add size
@@ -922,10 +929,25 @@ public:
 	 * \tparam CellL CellList type to construct
 	 *
 	 */
-	template<typename CellL=CellList<dim,St,FAST>> CellL getCellList()
+	template<typename CellL=CellList<dim,St,FAST>> CellL getCellList(St r_cut)
 	{
 		CellL cell_list;
 
+		// calculate the parameters of the cell list
+
+		// get the processor bounding box
+		Box<dim,St> pbox = dec.getProcessorBounds();
+		Box<dim,St> pbox_dom = pbox;
+		pbox_dom -= pbox_dom.getP1();
+
+		size_t div[dim];
+		// Calculate the division array
+
+		for (size_t i = 0 ; i < dim ; i++)
+			div[i] = (pbox.getP2().get(i) - pbox.getP1().get(i))/ r_cut;
+
+		cell_list.Initialize(pbox_dom,div,pbox.getP1());
+
 		// for each particle add the particle to the cell list
 
 		auto it = getDomainIterator();
diff --git a/src/Vector/vector_dist_unit_test.hpp b/src/Vector/vector_dist_unit_test.hpp
index 9d810f51..f0e2b3cb 100644
--- a/src/Vector/vector_dist_unit_test.hpp
+++ b/src/Vector/vector_dist_unit_test.hpp
@@ -508,7 +508,10 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_use_2d )
 
 		// Ghost must populated because we synchronized them
 		if (k > 524288)
+		{
 			BOOST_REQUIRE(nl_cnt != 0);
+			BOOST_REQUIRE(l_cnt > nl_cnt);
+		}
 
 		// Sum all the particles inside the domain
 		v_cl.sum(l_cnt);
@@ -529,7 +532,9 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_use_2d )
 
 		// Ghost must be populated
 		if (k > 524288)
+		{
 			BOOST_REQUIRE(nl_cnt != 0);
+		}
 	}
 }
 
@@ -613,7 +618,10 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_use_3d )
 
 		// Ghost must populated because we synchronized them
 		if (k > 524288)
+		{
 			BOOST_REQUIRE(nl_cnt != 0);
+			BOOST_REQUIRE(l_cnt > nl_cnt);
+		}
 
 		// Sum all the particles inside the domain
 		v_cl.sum(l_cnt);
@@ -632,7 +640,9 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_use_3d )
 
 		// Ghost must be populated
 		if (k > 524288)
+		{
 			BOOST_REQUIRE(nl_cnt != 0);
+		}
 	}
 }
 
@@ -824,7 +834,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_not_periodic_map )
 
 BOOST_AUTO_TEST_CASE( vector_dist_out_of_bound_policy )
 {
-/*	Vcluster & v_cl = *global_v_cluster;
+	Vcluster & v_cl = *global_v_cluster;
 
 	if (v_cl.getProcessingUnits() > 8)
 		return;
@@ -849,14 +859,26 @@ BOOST_AUTO_TEST_CASE( vector_dist_out_of_bound_policy )
 
 	auto it = vd.getIterator();
 
+	size_t cnt = 0;
+
 	while (it.isNext())
 	{
 		auto key = it.get();
 
-		vd.template getPos<s::x>(key)[0] = -0.06;
-		vd.template getPos<s::x>(key)[1] = -0.06;
-		vd.template getPos<s::x>(key)[2] = -0.06;
+		if (cnt < 1)
+		{
+			vd.template getPos<s::x>(key)[0] = -0.06;
+			vd.template getPos<s::x>(key)[1] = -0.06;
+			vd.template getPos<s::x>(key)[2] = -0.06;
+		}
+		else
+		{
+			vd.template getPos<s::x>(key)[0] = 0.06;
+			vd.template getPos<s::x>(key)[1] = 0.06;
+			vd.template getPos<s::x>(key)[2] = 0.06;
+		}
 
+		cnt++;
 		++it;
 	}
 
@@ -864,12 +886,145 @@ BOOST_AUTO_TEST_CASE( vector_dist_out_of_bound_policy )
 
 	// Particles out of the boundary are killed
 
-	size_t cnt = vd.size_local();
+	size_t cnt_l = vd.size_local();
 
-	v_cl.sum(cnt);
+	v_cl.sum(cnt_l);
 	v_cl.execute();
 
-	BOOST_REQUIRE_EQUAL(cnt,0);*/
+	BOOST_REQUIRE_EQUAL(cnt_l,100-v_cl.getProcessingUnits());
+}
+
+BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_interacting_particles )
+{
+
+	typedef Point<3,float> s;
+
+	Vcluster & v_cl = *global_v_cluster;
+
+	if (v_cl.getProcessingUnits() > 8)
+		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);
+
+	size_t nsz[] = {0,32,4};
+	nsz[0] = 4096 * v_cl.getProcessingUnits();
+
+	// 3D test
+	for (size_t i = 0 ; i < 3 ; i++ )
+	{
+		size_t k = nsz[i];
+
+		BOOST_TEST_CHECKPOINT( "Testing 3D random walk vector 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};
+
+		// factor
+		float factor = pow(global_v_cluster->getProcessingUnits()/2.0f,1.0f/3.0f);
+
+		// interaction radius
+		float r_cut = 0.01 / factor;
+
+		// ghost
+		Ghost<3,float> ghost(r_cut);
+
+		// Distributed vector
+		vector_dist<3,float, Point_test<float>, CartDecomposition<3,float> > vd(k,box,bc,ghost);
+
+		auto it = vd.getIterator();
+
+		while (it.isNext())
+		{
+			auto key = it.get();
+
+			vd.template getPos<s::x>(key)[0] = ud(eg);
+			vd.template getPos<s::x>(key)[1] = ud(eg);
+			vd.template getPos<s::x>(key)[2] = ud(eg);
+
+			++it;
+		}
+
+		vd.map();
+
+		// 4 step random walk
+
+		for (size_t j = 0 ; j < 4 ; j++)
+		{
+			auto it = vd.getDomainIterator();
+
+			// Move the particles
+
+			while (it.isNext())
+			{
+				auto key = it.get();
+
+				vd.template getPos<s::x>(key)[0] += 0.02 * ud(eg);
+				vd.template getPos<s::x>(key)[1] += 0.02 * ud(eg);
+				vd.template getPos<s::x>(key)[2] += 0.02 * ud(eg);
+
+				++it;
+			}
+
+			vd.map();
+
+			vd.ghost_get<0>();
+
+			// get the cell list with a cutoff radius
+
+			auto NN = vd.getCellList(0.01 / factor);
+
+			// iterate across the domain particle
+
+			auto it2 = vd.getDomainIterator();
+
+			while (it2.isNext())
+			{
+				auto p = it2.get();
+
+				Point<3,float> xp = vd.getPos<0>(p);
+
+				// reset the force
+
+				vd.getProp<4>(p)[0] = 0;
+				vd.getProp<4>(p)[1] = 0;
+				vd.getProp<4>(p)[2] = 0;
+
+				auto Np = NN.getIterator(NN.getCell(vd.getPos<0>(p)));
+
+				while (Np.isNext())
+				{
+					auto q = Np.get();
+
+					// repulsive
+
+					Point<3,float> xq = vd.getPos<0>(q);
+					Point<3,float> f = (xp - xq);
+
+					float distance = f.norm();
+
+					// Particle should be in the r_cut and
+					// 2 * r_cut range
+
+
+
+					++Np;
+				}
+
+				++it2;
+			}
+
+			// Count the local particles and check that the total number is consistent
+			size_t cnt = total_n_part_lc(vd,bc);
+
+			BOOST_REQUIRE_EQUAL((size_t)k,cnt);
+		}
+	}
 }
 
 BOOST_AUTO_TEST_SUITE_END()
-- 
GitLab