diff --git a/configure.ac b/configure.ac index b6332f602ffe35189a81bc3c6c2d236ec91f4976..fa7c74c37f6ebc5dd986f48d06ac590368b8c5f0 100755 --- a/configure.ac +++ b/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]) diff --git a/src/Draw/DrawParticles.hpp b/src/Draw/DrawParticles.hpp new file mode 100644 index 0000000000000000000000000000000000000000..fad43454f0aaa37819be02f44e44a45627265336 --- /dev/null +++ b/src/Draw/DrawParticles.hpp @@ -0,0 +1,66 @@ +/* + * 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_ */ diff --git a/src/Draw/DrawParticlesGrid.hpp b/src/Draw/DrawParticlesGrid.hpp deleted file mode 100644 index 1767d2eee4686dfbc439425ea406f80d9f94a189..0000000000000000000000000000000000000000 --- a/src/Draw/DrawParticlesGrid.hpp +++ /dev/null @@ -1,75 +0,0 @@ -/* - * 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_ */ diff --git a/src/Draw/DrawParticles_unit_tests.hpp b/src/Draw/DrawParticles_unit_tests.hpp index e08d9c38c31d95cf08a9ab8445033036c40c55b9..9423e637bff6c575817d678d05dae6702e812a7f 100644 --- a/src/Draw/DrawParticles_unit_tests.hpp +++ b/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); } diff --git a/src/Draw/PointIterator.hpp b/src/Draw/PointIterator.hpp index 6a68625cb03ed66899f109503fadf8ec3b0c3e96..c753e7d9342fcc1da11a706d2b06785ee6018b6e 100644 --- a/src/Draw/PointIterator.hpp +++ b/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; + } }; diff --git a/src/Draw/PointIteratorSkin.hpp b/src/Draw/PointIteratorSkin.hpp new file mode 100644 index 0000000000000000000000000000000000000000..d1c9979a6678ac7445eaafe49762b8535f4b9a70 --- /dev/null +++ b/src/Draw/PointIteratorSkin.hpp @@ -0,0 +1,198 @@ +/* + * 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_ */ diff --git a/src/FiniteDifference/FDScheme.hpp b/src/FiniteDifference/FDScheme.hpp index 66930c3ea8ec2cc0e39b4e5690796950d8f02874..93eda7ebf197e4b8f08b53813d45d110927fc67e 100644 --- a/src/FiniteDifference/FDScheme.hpp +++ b/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" diff --git a/src/Makefile.am b/src/Makefile.am index d631e1b8b63ff3e9792375ef401153d2beb7d2e8..d802738e4a88986d8fed6e1501633cc39a4ccaae 100755 --- a/src/Makefile.am +++ b/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 $< diff --git a/src/main.cpp b/src/main.cpp index 30e4ef4b77135a1fd2d083462998cbb36aa5d1d6..bb616a87c5692cafcda541b54e771e2e19318f0e 100755 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,9 +1,9 @@ -#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"