Actual numerics branch

configure.ac

@@ -7,7 +7,7 @@ AC_PREREQ(2.59)
 AC_INIT(FULL-PACKAGE-NAME, VERSION, BUG-REPORT-ADDRESS)
 AC_CANONICAL_SYSTEM
 AC_CONFIG_SRCDIR([src/main.cpp])
-AC_CONFIG_SUBDIRS([openfpm_data openfpm_devices openfpm_vcluster openfpm_io])
+AC_CONFIG_SUBDIRS([openfpm_data openfpm_devices openfpm_vcluster openfpm_io openfpm_numerics])
 
 #### Adding --with-pdata option and openfpm_pdata to prefix folder
 

install

@@ -210,8 +210,8 @@ fi
 
 ### Create example.mk
 install_base=$(cat install_dir)
-echo "INCLUDE_PATH=-I. -I$install_base/openfpm_pdata/include/config -I$install_base/openfpm_pdata/include -I$install_base/openfpm_data/include -I$install_base/openfpm_vcluster/include -I$install_base/openfpm_io/include -I$install_base/openfpm_devices/include -I$i_dir/METIS/include -I$i_dir/BOOST/include" > example.mk
-echo "LIBS_PATH=-L$install_base/openfpm_devices/lib -L$install_base/openfpm_vcluster/lib -L$i_dir/METIS/lib -L$i_dir/BOOST/lib " >> example.mk
+echo "INCLUDE_PATH=-I. -I$install_base/openfpm_numerics/include -I$install_base/openfpm_pdata/include/config -I$install_base/openfpm_pdata/include -I$install_base/openfpm_data/include -I$install_base/openfpm_vcluster/include -I$install_base/openfpm_io/include -I$install_base/openfpm_devices/include -I$i_dir/METIS/include -I$i_dir/BOOST/include" > example.mk
+echo "LIBS_PATH= -L$install_base/openfpm_devices/lib -L$install_base/openfpm_vcluster/lib -L$i_dir/METIS/lib -L$i_dir/BOOST/lib " >> example.mk
 echo "LIBS=-lvcluster -lofpmmemory -lmetis -lboost_iostreams" >> example.mk
 echo "LIBS_SE2=-lvcluster -lofpmmemory_se2 -lmetis -lboost_iostreams" >> example.mk
 cp example.mk src/example.mk

src/Decomposition/CartDecomposition.hpp

@@ -505,7 +505,7 @@ public:
 	 * \param p Point to apply the boundary condition
 	 *
 	 */
-	void applyPointBC(Point<dim,T> & pt) const
+	bool applyPointBC(Point<dim,T> & pt) const
 	{
 		for (size_t i = 0 ; i < dim ; i++)
 		{
@@ -802,9 +802,9 @@ p1[0]<-----+         +----> p2[0]
 	 * \return processorID
 	 *
 	 */
-	template<typename Mem> size_t inline processorID(encapc<1, Point<dim,T>, Mem> p)
+	template<typename Mem, typename ofb> size_t inline processorID(encapc<1, Point<dim,T>, Mem> p)
 	{
-		return fine_s.get(cd.getCell(p));
+		return fine_s.get(cd.template getCell<ofb>(p));
 	}
 
 	/*! \brief Given a point return in which processor the particle should go
@@ -849,7 +849,7 @@ p1[0]<-----+         +----> p2[0]
 	 * \return processorID
 	 *
 	 */
-	size_t inline processorIDBC(const Point<dim,T> &p) const
+	template<typename ofb> size_t inline processorIDBC(const Point<dim,T> &p) const
 	{
 		Point<dim,T> pt = p;
 		applyPointBC(pt);
@@ -864,7 +864,7 @@ p1[0]<-----+         +----> p2[0]
 	 * \return processorID
 	 *
 	 */
-	size_t inline processorIDBC(const T (&p)[dim]) const
+	template<typename ofb> size_t inline processorIDBC(const T (&p)[dim]) const
 	{
 		Point<dim,T> pt = p;
 		applyPointBC(pt);

src/Decomposition/CartDecomposition_unit_test.hpp

@@ -100,8 +100,6 @@ BOOST_AUTO_TEST_CASE( CartDecomposition_non_periodic_test)
 	ret = dec3.check_consistency();
 	BOOST_REQUIRE_EQUAL(ret,true);
 
-	dec3.write("output_dec3_1.vtk");
-
 	// Check that dec3 is equal to dec2 with the exception of the ghost part
 	ret = dec3.is_equal_ng(dec2);
 	BOOST_REQUIRE_EQUAL(ret,true);

src/Decomposition/ie_ghost.hpp

@@ -393,6 +393,51 @@ protected:
 
 public:
 
+	//! Default constructor
+	ie_ghost() {};
+
+	//! Copy constructor
+	ie_ghost(const ie_ghost<dim,T> & ie)
+	{
+		this->operator =(ie);
+	}
+
+	//! Copy constructor
+	ie_ghost(ie_ghost<dim,T> && ie)
+	{
+		this->operator=(ie);
+	}
+
+	//! Copy operator
+	inline ie_ghost<dim,T> & operator=(ie_ghost<dim,T> && ie)
+	{
+		box_nn_processor_int.swap(ie.box_nn_processor_int);
+		proc_int_box.swap(ie.proc_int_box);
+		vb_ext.swap(ie.vb_ext);
+		vb_int.swap(ie.vb_int);
+		geo_cell.swap(ie.geo_cell);
+		shifts.swap(ie.shifts);
+		ids_p.swap(ie.ids_p);
+		ids.swap(ie.ids);
+
+		return *this;
+	}
+
+	//! Copy operator
+	inline ie_ghost<dim,T> & operator=(const ie_ghost<dim,T> & ie)
+	{
+		box_nn_processor_int = ie.box_nn_processor_int;
+		proc_int_box = ie.proc_int_box;
+		vb_ext = ie.vb_ext;
+		vb_int = ie.vb_int;
+		geo_cell = geo_cell;
+		shifts = ie.shifts;
+		ids_p = ie.ids_p;
+		ids = ie.ids;
+
+		return *this;
+	}
+
 	/*! It return the shift vector
 	 *
 	 * Consider a domain with some ghost, at the border of the domain the

src/Vector/vector_dist.hpp

@@ -21,6 +21,7 @@
 #include "memory/ExtPreAlloc.hpp"
 #include "CSVWriter.hpp"
 #include "Decomposition/common.hpp"
+#include "Vector/vector_dist_ofb.hpp"
 
 #define V_SUB_UNIT_FACTOR 64
 
@@ -138,7 +139,7 @@ private:
 	 * \param opart id of the particles to send
 	 *
 	 */
-	void labelParticleProcessor(openfpm::vector<openfpm::vector<size_t>> & lbl_p, openfpm::vector<size_t> & prc_sz, openfpm::vector<size_t> & opart)
+	template<typename obp> void labelParticleProcessor(openfpm::vector<openfpm::vector<size_t>> & lbl_p, openfpm::vector<size_t> & prc_sz, openfpm::vector<size_t> & opart)
 	{
 		// reset lbl_p
 		lbl_p.resize(v_cl.getProcessingUnits());
@@ -158,7 +159,14 @@ private:
 			// Apply the boundary conditions
 			dec.applyPointBC(v_pos.get(key));
 
-			size_t p_id = dec.processorIDBC(v_pos.get(key));
+			size_t p_id = 0;
+
+			// Check if the particle is inside the domain
+//			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());
+
 
 			// Particle to move
 			if (p_id != v_cl.getProcessUnitID())
@@ -674,7 +682,7 @@ public:
 	 */
 	size_t size_local()
 	{
-		return v_pos.get(0).size();
+		return g_m;
 	}
 
 	/*! \brief Get the position of an element
@@ -707,13 +715,15 @@ public:
 	}
 
 	/*! \brief It move all the particles that does not belong to the local processor to the respective processor
+	 *
+	 * \tparam out of bound policy it specify what to do when the particles are detected out of bound
 	 *
 	 * In general this function is called after moving the particles to move the
 	 * elements out the local processor. Or just after initialization if each processor
 	 * contain non local particles
 	 *
 	 */
-	void map()
+	template<typename obp=KillParticle> void map()
 	{
 		// outgoing particles-id
 		openfpm::vector<size_t> out_part;
@@ -729,7 +739,7 @@ public:
 		v_prp.resize(g_m);
 
 		// Contain the processor id of each particle (basically where they have to go)
-		labelParticleProcessor(opart,prc_sz,out_part);
+		labelParticleProcessor<obp>(opart,prc_sz,out_part);
 
 		// Calculate the sending buffer size for each processor, put this information in
 		// a contiguous buffer
@@ -907,6 +917,30 @@ public:
 		return vd->recv_mem_gm.get(i).getPointer();
 	}
 
+	/*! \brief Construct a cell list starting from the stored particles
+	 *
+	 * \tparam CellL CellList type to construct
+	 *
+	 */
+	template<typename CellL=CellList<dim,St,FAST>> CellL getCellList()
+	{
+		CellL cell_list;
+
+		// for each particle add the particle to the cell list
+
+		auto it = getDomainIterator();
+
+		while (it.isNext())
+		{
+			auto key = it.get();
+
+			cell_list.add(this->template getPos<0>(key),key.getKey());
+
+			++it;
+		}
+
+		return cell_list;
+	}
 
 	/*! \brief Get the iterator across the position of the particles
 	 *

src/Vector/vector_dist_unit_test.hpp

@@ -822,6 +822,56 @@ 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;
+
+	if (v_cl.getProcessingUnits() > 8)
+		return;
+
+	typedef Point<3,float> s;
+
+	Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
+
+	// Boundary conditions
+	size_t bc[3]={NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
+
+	// factor
+	float factor = pow(global_v_cluster->getProcessingUnits()/2.0f,1.0f/3.0f);
+
+	// ghost
+	Ghost<3,float> ghost(0.05 / factor);
+
+	// Distributed vector
+	vector_dist<3,float, Point_test<float>, CartDecomposition<3,float> > vd(100,box,bc,ghost);
+
+	// put particles at out of the boundary, they must be detected and and killed
+
+	auto it = vd.getIterator();
+
+	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;
+
+		++it;
+	}
+
+	vd.map();
+
+	// Particles out of the boundary are killed
+
+	size_t cnt = vd.size_local();
+
+	v_cl.sum(cnt);
+	v_cl.execute();
+
+	BOOST_REQUIRE_EQUAL(cnt,0);*/
+}
+
 BOOST_AUTO_TEST_SUITE_END()
 
 #endif /* VECTOR_DIST_UNIT_TEST_HPP_ */