Fixing draw particles component

configure.ac

@@ -4,7 +4,7 @@
 ## Take all the options with the exception of --enable-install-req
 
 AC_PREREQ(2.59)
-AC_INIT(FULL-PACKAGE-NAME, VERSION, BUG-REPORT-ADDRESS)
+AC_INIT(OpenFPM_numerics, 0.8.0, BUG-REPORT-ADDRESS)
 AC_CANONICAL_SYSTEM
 AC_CONFIG_SRCDIR([src/main.cpp])
 

src/Draw/DrawParticles.hpp

+/*
+ * DrawParticles.hpp
+ *
+ *  Created on: Jan 5, 2017
+ *      Author: i-bird
+ */
+
+#ifndef OPENFPM_NUMERICS_SRC_DRAW_DRAWPARTICLES_HPP_
+#define OPENFPM_NUMERICS_SRC_DRAW_DRAWPARTICLES_HPP_
+
+#include "PointIterator.hpp"
+#include "PointIteratorSkin.hpp"
+#include "Vector/vector_dist.hpp"
+
+class DrawParticles
+{
+public:
+
+	template<unsigned int dim, typename T, typename aggr, typename Decomposition> static PointIteratorSkin<dim,T,Decomposition>
+	DrawSkin(vector_dist<dim,T,aggr,Decomposition> & vd,
+			 size_t (& sz)[dim],
+			 Box<dim,T> & domain,
+			 Box<dim,T> & sub_A,
+			 Box<dim,T> & sub_B)
+	{
+		size_t bc[dim];
+
+		for (size_t i = 0 ; i < dim ; i++)
+			bc[i] = NON_PERIODIC;
+
+		return PointIteratorSkin<dim,T,Decomposition>(vd.getDecomposition(),sz,domain,sub_A, sub_B, bc);
+	}
+
+	template<unsigned int dim, typename T, typename aggr, typename Decomposition> static PointIteratorSkin<dim,T,Decomposition>
+	DrawSkin(vector_dist<dim,T,aggr,Decomposition> & vd,
+			 size_t (& sz)[dim],
+			 Box<dim,T> & domain,
+			 openfpm::vector<Box<dim,T>> & sub_A,
+			 Box<dim,T> & sub_B)
+	{
+		size_t bc[dim];
+
+		for (size_t i = 0 ; i < dim ; i++)
+			bc[i] = NON_PERIODIC;
+
+		PointIteratorSkin<dim,T,Decomposition> it(vd.getDecomposition(),sz,domain,sub_A.get(0), sub_B, bc);
+
+		for (size_t i = 1 ; i < sub_A.size() ; i++)
+			it.addBoxA(Box<dim,T>(sub_A.get(i)));
+
+		return it;
+	}
+
+	template<unsigned int dim, typename T, typename aggr, typename Decomposition> static PointIterator<dim,T,Decomposition>
+	DrawBox(vector_dist<dim,T,aggr,Decomposition> & vd,
+			 size_t (& sz)[dim],
+			 Box<dim,T> & domain,
+			 Box<dim,T> & sub)
+	{
+		return PointIterator<dim,T,Decomposition>(vd.getDecomposition(),sz,domain,sub);
+	}
+
+};
+
+
+#endif /* OPENFPM_NUMERICS_SRC_DRAW_DRAWPARTICLES_HPP_ */

src/Draw/DrawParticlesGrid.hpp

-/*
- * DrawParticlesGrid.hpp
- *
- *  Created on: Jan 3, 2017
- *      Author: i-bird
- */
-
-#ifndef OPENFPM_NUMERICS_SRC_DRAW_DRAWPARTICLESGRID_HPP_
-#define OPENFPM_NUMERICS_SRC_DRAW_DRAWPARTICLESGRID_HPP_
-
-/*! \brief this class draw particles on subset of grid-like position
- *
- \verbatim
-                  Point B
-                    |                 (0.6,0.6)
-  +   +   +   +   + |  +   +   +   +
-    +---------------+
-  + | +   +   +   + | +   +   +   +
-    |               |
-  + | +   +   +   + | +   +   +   +
-    |               |
-  + | +   +   +   + | +   +   +   +
-    |               |
-  + | +   +   +   + | +   +   +   +
-    |               |
-  + | +   +   +   + | +   +   +   +
-    |               |
-  + | +   +   +   + | +   +   +   +
-    |               |
-  + | +   +   +   + | +   +   +   +
-    +---------------+
-  +   +   +   +   + ^ +   +   +   +
-(-1.0,-1.0)         |
-                    |
-                   (Point A)
-
-
- \endverbatim
- *
- *
- *  Suppose we have a grid 9x9 from (-1.0,-1.0 to 0.6,0.6)
- *
- *  Defined a Box (-0.9,-0.9 to -0.1,0.5)
- *
- *  This class will return the points
- *
- *  (-0.8 -0.8)
- *  (-0.6 -0.8)
- *  (-0.4 -0.8)
- *  ...
- *  (-0.1,-0.8) Point A
- *  ...
- *  (-0.1,0.5) Point B
- *
- */
-template<unsigned int dim, typename T>
-class DrawParticlesGrid
-{
-	/*! \brief Draw Particles
-	 *
-	 * \param sp grid spacing
-	 *
-	 */
-	DrawParticlesGrid(T (& sp)[dim], Box<dim,T> & domain)
-	{
-	}
-
-	PointIterator getPointIterator()
-	{
-
-	}
-};
-
-
-#endif /* OPENFPM_NUMERICS_SRC_DRAW_DRAWPARTICLESGRID_HPP_ */

src/Draw/DrawParticles_unit_tests.hpp

@@ -9,15 +9,20 @@
 #define OPENFPM_NUMERICS_SRC_DRAW_DRAWPARTICLES_UNIT_TESTS_HPP_
 
 #include "PointIterator.hpp"
+#include "PointIteratorSkin.hpp"
+#include "DrawParticles.hpp"
 
 BOOST_AUTO_TEST_SUITE( draw_particles )
 
 BOOST_AUTO_TEST_CASE(point_iterator)
 {
-	size_t sz[] = {30,17,28};
+	size_t sz[] = {23,27,20};
 
-	Box<3,double> domain({-1.2,0.5,-0.6},{1.0,3.1,1.312});
-	Box<3,double> sub_domain({-0.2,0.8,0.2},{1.0,1.1,1.0});
+	Box<3,double> domain({-1.2,0.5,-0.6},{1.0,3.1,1.3});
+	Box<3,double> sub_domain({-0.15,0.75,0.15},{1.05,1.15,1.05});
+	Box<3,double> sub_real({-0.1,0.80,0.2},{1.0,1.1,1.0});
+
+	size_t sz_sub[] = {12,4,9};
 
 	// Boundary conditions
 	size_t bc[3]={NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
@@ -27,7 +32,7 @@ BOOST_AUTO_TEST_CASE(point_iterator)
 
 	vector_dist<3,double,aggregate<double>> vd(0,domain,bc,ghost);
 
-	PointIterator<3,double,CartDecomposition<3,double>> p(vd.getDecomposition(),sz,domain,sub_domain);
+	auto p = DrawParticles::DrawBox(vd,sz,domain,sub_domain);
 
 	size_t cnt = 0;
 
@@ -50,8 +55,67 @@ BOOST_AUTO_TEST_CASE(point_iterator)
 	Vcluster & v_cl = create_vcluster();
 
 	v_cl.sum(cnt);
+	v_cl.execute();
+
+	BOOST_REQUIRE_EQUAL(cnt,sz_sub[0]*sz_sub[1]*sz_sub[2]);
+	BOOST_REQUIRE_EQUAL(good,true);
+
+	Box<3,double> marg = p.getBoxMargins();
+
+	BOOST_REQUIRE_CLOSE(marg.getLow(0), sub_real.getLow(0) ,0.0001);
+	BOOST_REQUIRE_CLOSE(marg.getLow(1), sub_real.getLow(1) ,0.0001);
+	BOOST_REQUIRE_CLOSE(marg.getLow(2), sub_real.getLow(2) ,0.0001);
+
+	BOOST_REQUIRE_CLOSE(marg.getHigh(0), sub_real.getHigh(0) ,0.0001);
+	BOOST_REQUIRE_CLOSE(marg.getHigh(1), sub_real.getHigh(1) ,0.0001);
+	BOOST_REQUIRE_CLOSE(marg.getHigh(2), sub_real.getHigh(2) ,0.0001);
+}
+
+BOOST_AUTO_TEST_CASE(point_iterator_skin)
+{
+	size_t sz[] = {23,27,20};
+
+	Box<3,double> domain({-1.2,0.5,-0.6},{1.0,3.1,1.3});
+	Box<3,double> sub_domainA({-0.15,0.75,0.15},{1.05,1.15,1.05});
+	Box<3,double> sub_domainB({-0.25,0.65,0.05},{0.95,1.05,1.05});
+
+	// Boundary conditions
+	size_t bc[3]={NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
+
+	// ghost, big enough to contain the interaction radius
+	Ghost<3,float> ghost(0.01);
+
+	vector_dist<3,double,aggregate<double>> vd(0,domain,bc,ghost);
+
+	auto p = DrawParticles::DrawSkin(vd,sz,domain,sub_domainA, sub_domainB);
+
+	size_t cnt = 0;
+
+	bool good = true;
+	while (p.isNext())
+	{
+		vd.add();
+
+		vd.getLastPos()[0] = p.get().get(0);
+		vd.getLastPos()[1] = p.get().get(1);
+		vd.getLastPos()[2] = p.get().get(2);
+
+		if (!((p.get().get(0) > -0.21 && p.get().get(0) < -0.19) || (p.get().get(1) > 0.69 && p.get().get(1) < 0.71) || (p.get().get(2) > 0.09 && p.get().get(2) < 0.11)) )
+			good &= false;
+
+		if (sub_domainA.isInside(p.get()))
+			good &= false;
+
+		cnt++;
+
+		++p;
+	}
+
+	Vcluster & v_cl = create_vcluster();
+
+	v_cl.sum(cnt);
+	v_cl.execute();
 
-	BOOST_REQUIRE_EQUAL(cnt,sz[0]*sz[1]*sz[2]);
 	BOOST_REQUIRE_EQUAL(good,true);
 }
 

src/Draw/PointIterator.hpp

@@ -63,15 +63,17 @@ class PointIterator: protected grid_dist_id_iterator_dec<Decomposition>
 	//! sub_domain
 	Box<dim,T> sub_domain;
 
+	//! domain
+	Box<dim,T> domain;
+
 	//! Spacing
 	T sp[dim];
 
 	static grid_key_dx<dim> getStart(size_t (& gs)[dim], Box<dim,T> & dom, Box<dim,T> & sub_dom, T (& sp)[dim])
 	{
 		for (size_t i = 0 ; i < dim ; i++)
-			sp[i] = (dom.getHigh(i) - dom.getLow(i)) / gs[i];
+			sp[i] = (dom.getHigh(i) - dom.getLow(i)) / (gs[i] -1);
 
-		Point<dim,T> pstart = sub_dom.getP1();
 		grid_key_dx<dim> pkey;
 
 		for (size_t i = 0 ; i < dim ; i++)
@@ -83,17 +85,24 @@ class PointIterator: protected grid_dist_id_iterator_dec<Decomposition>
 	static grid_key_dx<dim> getStop(size_t (& gs)[dim], Box<dim,T> & dom, Box<dim,T> & sub_dom, T (& sp)[dim])
 	{
 		for (size_t i = 0 ; i < dim ; i++)
-			sp[i] = (dom.getHigh(i) - dom.getLow(i)) / gs[i];
+			sp[i] = (dom.getHigh(i) - dom.getLow(i)) / (gs[i] - 1);
 
-		Point<dim,T> pstop = sub_dom.getP2();
 		grid_key_dx<dim> pkey;
 
 		for (size_t i = 0 ; i < dim ; i++)
-			pkey.set_d(i,std::ceil( (sub_dom.getLow(i) - dom.getLow(i)) / sp[i]));
+			pkey.set_d(i,std::floor( (sub_dom.getHigh(i) - dom.getLow(i)) / sp[i]));
 
 		return pkey;
 	}
 
+	void calculateAp()
+	{
+		grid_key_dx<dim> key = grid_dist_id_iterator_dec<Decomposition>::get();
+
+		for (size_t i = 0 ; i < dim ; i++)
+			ap.get(i) = key.get(i) * sp[i] + domain.getLow(i);
+	}
+
 public:
 
 	/*! \brief Draw Particles
@@ -102,11 +111,9 @@ public:
 	 *
 	 */
 	PointIterator( Decomposition & dec, size_t (& sz)[dim], Box<dim,T> & domain, Box<dim,T> & sub_domain)
-	:grid_dist_id_iterator_dec<Decomposition>(dec, sz, getStart(sz,domain,sub_domain,sp), getStop(sz,domain,sub_domain,sp)),sub_domain(sub_domain)
+	:grid_dist_id_iterator_dec<Decomposition>(dec, sz, getStart(sz,domain,sub_domain,sp), getStop(sz,domain,sub_domain,sp)),sub_domain(sub_domain),domain(domain)
 	{
-		// Calculate the spacing
-		for (size_t i = 0 ; i < dim ; i++)
-			sp[i] = domain.getLow(i) / sz[i];
+		calculateAp();
 	}
 
 	/*! \Return the actual point
@@ -126,10 +133,7 @@ public:
 	PointIterator & operator++()
 	{
 		grid_dist_id_iterator_dec<Decomposition>::operator++();
-		grid_key_dx<dim> key = grid_dist_id_iterator_dec<Decomposition>::get();
-
-		for (size_t i = 0 ; i < dim ; i++)
-			ap.get(i) = key.get(i) * sp[i] + sub_domain.getLow(i);
+		calculateAp();
 
 		return *this;
 	}
@@ -138,6 +142,31 @@ public:
 	{
 		return grid_dist_id_iterator_dec<Decomposition>::isNext();
 	}
+
+	/*! \brief Return the real Margin of the box
+	 *
+	 * For example consider to have a domain (0.0,0.0) to (1.0,1.0)
+	 * 11 points in each direction (spacing 0.1). if we have a sub-domain
+	 * (0.15,0.15) to (0.55,0.55) getBoxMargins return the box (0.2,0.2) to
+	 * (0.5,0.5). the smallest box that enclose the points that the point
+	 * iterator is going to give
+	 *
+	 */
+	Box<dim,T> getBoxMargins()
+	{
+		Box<dim,T> box;
+
+		grid_key_dx<dim> start = grid_dist_id_iterator_dec<Decomposition>::getStart();
+		grid_key_dx<dim> stop = grid_dist_id_iterator_dec<Decomposition>::getStop();
+
+		for (size_t i = 0 ; i < dim ; i++)
+		{
+			box.setLow(i, start.get(i)*sp[i] + domain.getLow(i));
+			box.setHigh(i, stop.get(i)*sp[i] + domain.getLow(i));
+		}
+
+		return box;
+	}
 };
 
 

src/Draw/PointIteratorSkin.hpp

+/*
+ * PointIteratorSkin.hpp
+ *
+ *  Created on: Jan 4, 2017
+ *      Author: i-bird
+ */
+
+#ifndef OPENFPM_NUMERICS_SRC_DRAW_POINTITERATORSKIN_HPP_
+#define OPENFPM_NUMERICS_SRC_DRAW_POINTITERATORSKIN_HPP_
+
+#include "Grid/Iterators/grid_dist_id_iterator_dec_skin.hpp"
+
+#define RETURN_A 1
+#define RETURN_B 2
+
+/*! \brief this class draw particles on subset of grid-like position
+ *
+ \verbatim
+             (Box B)
+                |                    (0.6,0.6)
+  +   +   +   + | +   +   +   +   +
+    +---+-------+---+
+  + | * | +   + | * | +   +   +   +
+    |   |       |
+  + | * | +   + | * | +   +   +   +
+    |   |       |   |
+  + | * | +   + | * | +   +   +   +
+    |   |       |   |
+  + | * | +   + | * | +   +   +   +
+    |   |       |   |
+  + | * | +   + | * | +   +   +   +
+    |   |       |   |
+  + | * | +   + | * | +   +   +   +
+    |   +-------+   |
+  + | *   *   *   * | +   +   +   +
+    +---------------+
+  +   +   +   +   + | +   +   +   +
+(-1.0,-1.0)         |
+                    |
+                 (Box A)
+
+
+ \endverbatim
+ *
+ *
+ *  Suppose we have a grid 9x9 from (-1.0,-1.0 to 0.6,0.6)
+ *
+ *  Defined a Box A (-0.9,-0.9 to -0.1,0.5) and a box B (-0.7, -0.7 to -0.3,0.5)
+ *
+ *  This class will return the points indicated with *
+ *
+ */
+template<unsigned int dim, typename T, typename Decomposition>
+class PointIteratorSkin: protected grid_dist_id_iterator_dec_skin<Decomposition>
+{
+	//! Actual point
+	Point<dim,T> ap;
+
+	//! sub_domain (required to filter out points)
+	openfpm::vector<Box<dim,T>> sub_domainA;
+
+	//! domain
+	Box<dim,T> domain;
+
+	//! Spacing
+	T sp[dim];
+
+	static Box<dim,long int> getAB(size_t (& gs)[dim], const Box<dim,T> & dom, const Box<dim,T> & sub_domA , const Box<dim,T> & sub_domB, T (& sp)[dim], size_t AB)
+	{
+		for (size_t i = 0 ; i < dim ; i++)
+			sp[i] = (dom.getHigh(i) - dom.getLow(i)) / (gs[i] -1);
+
+		Box<dim,long int> box;
+
+		for (size_t i = 0 ; i < dim ; i++)
+		{
+			size_t Ast = std::ceil( (sub_domA.getLow(i) - dom.getLow(i)) / sp[i] ) - 1;
+			size_t Asp = std::floor( (sub_domA.getHigh(i) - dom.getLow(i)) / sp[i] ) + 1;
+
+			size_t Bst = std::ceil( (sub_domB.getLow(i) - dom.getLow(i)) / sp[i] );
+			size_t Bsp = std::floor( (sub_domB.getHigh(i) - dom.getLow(i)) / sp[i] );
+
+			// grid_dist_id_iterator_dec_skin only work if A is contained into B
+			Ast = (Ast < Bst)?Bst:Ast;
+			Asp = (Asp > Bsp)?Bsp:Asp;
+
+			if (AB == RETURN_A)
+			{
+				box.setLow(i,Ast);
+				box.setHigh(i,Asp);
+			}
+			else
+			{
+				box.setLow(i,Bst);
+				box.setHigh(i,Bsp);
+			}
+		}
+
+		return box;
+	}
+
+	void calculateAp()
+	{
+		grid_key_dx<dim> key = grid_dist_id_iterator_dec_skin<Decomposition>::get();
+
+		for (size_t i = 0 ; i < dim ; i++)
+			ap.get(i) = key.get(i) * sp[i] + domain.getLow(i);
+	}
+
+
+	/*! it check that the actual point is not inside B
+	 *
+	 * \return true if the point is not inside B
+	 *
+	 */
+	bool isValidPoint()
+	{
+		bool valid = true;
+
+		for (size_t i = 0 ; i < sub_domainA.size() ; i++)
+		{
+			if (Box<dim,T>(sub_domainA.get(i)).isInside(ap) == true)
+				valid = false;
+		}
+
+		return valid;
+	}
+
+public:
+
+	/*! \brief Draw Particles
+	 *
+	 * \param sp grid spacing
+	 *
+	 */
+	PointIteratorSkin( Decomposition & dec, size_t (& sz)[dim], Box<dim,T> & domain, const Box<dim,T> & sub_A, const Box<dim,T> & sub_B, size_t (& bc)[dim])
+	:grid_dist_id_iterator_dec_skin<Decomposition>(dec, sz, getAB(sz,domain,sub_A,sub_B,sp,RETURN_A), getAB(sz,domain,sub_A,sub_B,sp,RETURN_B), bc),domain(domain)
+	{
+		sub_domainA.add(sub_A);
+		calculateAp();
+	}
+
+	/*! \Return the actual point
+	 *
+	 *
+	 */
+	Point<dim,T> & get()
+	{
+		return ap;
+	}
+
+	/*! \brief Next point
+	 *
+	 * \return itself
+	 *
+	 */
+	PointIteratorSkin & operator++()
+	{
+		grid_dist_id_iterator_dec_skin<Decomposition>::operator++();
+		calculateAp();
+
+		while (grid_dist_id_iterator_dec_skin<Decomposition>::isNext() && isValidPoint() == false)
+		{
+			grid_dist_id_iterator_dec_skin<Decomposition>::operator++();
+			calculateAp();
+		}
+
+		return *this;
+	}
+
+	bool isNext()
+	{
+		return grid_dist_id_iterator_dec_skin<Decomposition>::isNext();
+	}
+
+	PointIteratorSkin<dim,T,Decomposition> & operator=(PointIteratorSkin<dim,T,Decomposition> & p)
+	{
+		grid_dist_id_iterator_dec_skin<Decomposition>::operator=(p);
+
+		ap = p.ap;
+		sub_domainA = p.sub_domainA;
+		domain = p.domain;
+
+		for (size_t i = 0 ; i < dim; i++)
+			sp[i] = p.sp[i];
+
+		return *this;
+	}
+
+	void addBoxA(const Box<dim,double> & BoxA)
+	{
+		sub_domainA.add(BoxA);
+	}
+};
+
+
+
+#endif /* OPENFPM_NUMERICS_SRC_DRAW_POINTITERATORSKIN_HPP_ */

src/FiniteDifference/FDScheme.hpp

@@ -10,7 +10,7 @@
 
 #include "../Matrix/SparseMatrix.hpp"
 #include "Grid/grid_dist_id.hpp"
-#include "Grid/grid_dist_id_iterator_sub.hpp"
+#include "Grid/Iterators/grid_dist_id_iterator_sub.hpp"
 #include "eq.hpp"
 #include "data_type/scalar.hpp"
 #include "NN/CellList/CellDecomposer.hpp"

src/Makefile.am

@@ -11,7 +11,8 @@ Vector/Vector_eigen.hpp Vector/Vector_petsc.hpp Vector/Vector_util.hpp Vector/Ve
 Solvers/umfpack_solver.hpp Solvers/petsc_solver.hpp \
 util/petsc_util.hpp util/linalgebra_lib.hpp util/util_num.hpp util/grid_dist_testing.hpp  \
 FiniteDifference/Average.hpp FiniteDifference/Derivative.hpp FiniteDifference/eq.hpp FiniteDifference/FDScheme.hpp FiniteDifference/Laplacian.hpp FiniteDifference/mul.hpp FiniteDifference/sum.hpp FiniteDifference/util/common.hpp \
-PSE/Kernels.hpp PSE/Kernels_test_util.hpp Operators/Vector/vector_dist_operators_extensions.hpp Operators/Vector/vector_dist_operators.hpp Operators/Vector/vector_dist_operators_apply_kernel.hpp Operators/Vector/vector_dist_operators_functions.hpp Operators/Vector/vector_dist_operator_assign.hpp
+PSE/Kernels.hpp PSE/Kernels_test_util.hpp Operators/Vector/vector_dist_operators_extensions.hpp Operators/Vector/vector_dist_operators.hpp Operators/Vector/vector_dist_operators_apply_kernel.hpp Operators/Vector/vector_dist_operators_functions.hpp Operators/Vector/vector_dist_operator_assign.hpp \
+Draw/DrawParticles.hpp Draw/PointIterator.hpp Draw/PointIteratorSkin.hpp
 
 .cu.o :
 	$(NVCC) $(NVCCFLAGS) -o $@ -c $<

src/main.cpp

-#include "config.h"
 #define BOOST_TEST_MODULE "C++ test module for OpenFPM_numerics project"
 #include <boost/test/included/unit_test.hpp>
 
 #include "unit_test_init_cleanup.hpp"
 
+#include "config.h"
 #include "Vector/Vector_unit_tests.hpp"
 #include "Matrix/SparseMatrix_unit_tests.hpp"
 #include "FiniteDifference/FDScheme_unit_tests.hpp"