From b8b1616b21a1d2f1daa93d887071504fbedb3494 Mon Sep 17 00:00:00 2001 From: Pietro Incardona Date: Sun, 2 Jul 2017 15:48:20 +0200 Subject: [PATCH] Interpolation new version --- src/interpolation/interpolation.hpp | 454 ++++++++++++++++++++++------ 1 file changed, 358 insertions(+), 96 deletions(-) diff --git a/src/interpolation/interpolation.hpp b/src/interpolation/interpolation.hpp index da5ecef..d96a957 100644 --- a/src/interpolation/interpolation.hpp +++ b/src/interpolation/interpolation.hpp @@ -16,18 +16,48 @@ constexpr int inte_m2p = 0; constexpr int inte_p2m = 1; +template +struct agg_arr +{ + T ele[n_ele]; +}; + +/*! \brief multiply the src by coeff for several types T + * + * \tparam type T + * + */ template struct mul_inte { + /*! \brief multiply the src by coeff for several types T + * + * \param result the result of the multiplication + * \param coeff coefficent to use for of the multiplication + * \param src source + * + */ inline static void value(T & result, const T & coeff, const T & src) { result += coeff * src; } }; +/*! \brief multiply the src by coeff for several types T + * + * \tparam type T + * + */ template struct mul_inte { + /*! \brief multiply the src by coeff for several types T + * + * \param result the result of the multiplication + * \param coeff coefficent to use for of the multiplication + * \param src source + * + */ inline static void value(T (& result)[N1], const T & coeff, const T (& src)[N1]) { for (size_t i = 0 ; i < N1 ; i++) @@ -35,9 +65,21 @@ struct mul_inte } }; +/*! \brief multiply the src by coeff for several types T + * + * \tparam type T + * + */ template struct mul_inte { + /*! \brief multiply the src by coeff for several types T + * + * \param result the result of the multiplication + * \param coeff coefficent to use for of the multiplication + * \param src source + * + */ inline static void value(T (& result)[N1][N2], const T & coeff, const T (& src)[N1][N2]) { for (size_t i = 0 ; i < N1 ; i++) @@ -46,34 +88,158 @@ struct mul_inte } }; -template +/*! \brief Class that select the operation to do differently if we are doing Mesh to particle (m2p) or particle to mesh (p2m) + * + * \tparam prp_g property of the grid to interpolate + * \tparam prp_v property of the vector to interpolate + * \param M2P or P2M + * + */ +template struct inte_template { - template inline static void value(grid & gd, + /*! \brief Evaluate the interpolation + * + * \tparam np number of kernel points in one direction + * \tparam grid type of grid + * \tparam vector type of vector + * \tparam iterator type of the iterator + * + * \param gd grid for interpolation + * \param vd vector for interpolation + * \param k_dist grid key grid point for interpolation + * \param key_p particle for interpolation + * \param a_int interpolation coefficients pre-calculated + * \param key indicate which pre-calculated coefficient we have to use + * + */ + template inline static void value(grid & gd, vector & vd, - const grid_dist_key_dx & k_dist, + const grid_dist_lin_dx & k_dist, iterator & key_p, - const grid_cpu> & a_int, - const key_type & key) + typename vector::stype (& a_int)[openfpm::math::pow(vector::dims,np)], + const size_t & key) { - mul_inte(k_dist))>::type>::value(gd.template get(k_dist),a_int.template get<0>(key),vd.template getProp(key_p)); + mul_inte(k_dist))>::type>::value(gd.template get(k_dist),a_int[key],vd.template getProp(key_p)); } }; -template -struct inte_template +/*! \brief Class that select the operation to do differently if we are doing Mesh to particle (m2p) or particle to mesh (p2m) + * + * \tparam prp_g property of the grid to interpolate + * \tparam prp_v property of the vector to interpolate + * \param M2P or P2M + * + */ +template +struct inte_template { - template inline static void value(grid & gd, + /*! \brief Evaluate the interpolation + * + * \tparam grid type of grid + * \tparam vector type of vector + * \tparam iterator type of the iterator + * \tparam key_type type of the key + * + * \param gd grid for interpolation + * \param vd vector for interpolation + * \param k_dist grid key grid point for interpolation + * \param key_p particle for interpolation + * \param a_int interpolation coefficients pre-calculated + * \param key indicate which pre-calculated coefficient we have to use + * + */ + template inline static void value(grid & gd, vector & vd, - const grid_dist_key_dx & k_dist, + const grid_dist_lin_dx & k_dist, iterator & key_p, - const grid_cpu> & a_int, - const key_type & key) + typename vector::stype (& a_int)[openfpm::math::pow(vector::dims,np)], + const size_t & key) + { + mul_inte(k_dist))>::type>::value(vd.template getProp(key_p),a_int[key],gd.template get(k_dist)); + } +}; + +/*! \brief Calculate aint + * + * + */ +template +struct calculate_aint +{ + static inline void value(size_t (& sz)[vector::dims], + typename vector::stype a_int[openfpm::math::pow(vector::dims,np)], + typename vector::stype (& a)[vector::dims][np]) + { + grid_sm gs(sz); + grid_key_dx_iterator kit(gs); + + size_t s = 0; + while (kit.isNext()) + { + auto key = kit.get(); + + a_int[s] = 1; + + for (size_t i = 0 ; i < vector::dims ; i++) + a_int[s] *= a[i][key.get(i)]; + + s++; + ++kit; + } + } +}; + +/*! \brief Calculate aint 2D + * + * + */ +template +struct calculate_aint<2,vector,np> +{ + static inline void value(size_t (& sz)[vector::dims], + typename vector::stype a_int[openfpm::math::pow(vector::dims,np)], + typename vector::stype (& a)[vector::dims][np]) { - mul_inte(k_dist))>::type>::value(vd.template getProp(key_p),a_int.template get<0>(key),gd.template get(k_dist)); + size_t s = 0; + for (size_t i = 0 ; i < np ; i++) + { + for (size_t j = 0 ; j < np ; j++) + { + a_int[s] = a[0][j]*a[1][i]; + + s++; + } + } } }; +/*! \brief Calculate aint + * + * + */ +template +struct calculate_aint<3,vector,np> +{ + static inline void value(size_t (& sz)[vector::dims], + typename vector::stype a_int[openfpm::math::pow(vector::dims,np)], + typename vector::stype (& a)[vector::dims][np]) + { + size_t s = 0; + for (size_t i = 0 ; i < np ; i++) + { + for (size_t j = 0 ; j < np ; j++) + { + for (size_t k = 0 ; k < np ; k++) + { + a_int[s] = a[0][k]*a[1][j]*a[2][i]; + + s++; + } + } + } + } +}; /*! \brief return the sub-domain where this particle must be interpolated * @@ -105,113 +271,133 @@ inline size_t getSub(Point & p, return (size_t)-1; } - -template - inline void inte_calc(iterator & it, - vector & vd, - const Box & domain, - int (& ip)[vector::dims][kernel::np], - grid & gd, - const typename vector::stype (& dx)[vector::dims], - typename vector::stype (& xp)[vector::dims], - const grid_inte & a_int, - typename vector::stype (& a)[vector::dims][kernel::np], - typename vector::stype (& x)[vector::dims][kernel::np], - size_t (& sz)[vector::dims], - const CellList,shift> & geo_cell) +template +struct inte_calc_impl { - auto key_p = it.get(); - - Point p = vd.getPos(key_p); - - // On which sub-domain we interpolate the particle - size_t sub = getSub(p,geo_cell,gd); + //! Type of the calculations + typedef typename vector::stype arr_type; + + template + static inline void inte_calc(iterator & it, + vector & vd, + const Box & domain, + int (& ip)[vector::dims][kernel::np], + grid & gd, + const typename vector::stype (& dx)[vector::dims], + typename vector::stype (& xp)[vector::dims], + typename vector::stype (& a_int)[openfpm::math::pow(vector::dims,kernel::np)], + typename vector::stype (& a)[vector::dims][kernel::np], + typename vector::stype (& x)[vector::dims][kernel::np], + size_t (& sz)[vector::dims], + const CellList,shift> & geo_cell, + openfpm::vector> & offsets) + { + auto key_p = it.get(); - typename vector::stype x0[vector::dims]; + Point p = vd.getPos(key_p); - // calculate the position of the particle in the grid - // coordinated - for (size_t i = 0 ; i < vector::dims ; i++) - x0[i] = (p.get(i)-domain.getLow(i))*dx[i]; + // On which sub-domain we interpolate the particle + size_t sub = getSub(p,geo_cell,gd); - // convert into integer - for (size_t i = 0 ; i < vector::dims ; i++) - ip[i][0] = (int)x0[i]; + typename vector::stype x0[vector::dims]; - // convert the global grid position into local grid position - grid_key_dx base; + // calculate the position of the particle in the grid + // coordinated + for (size_t i = 0 ; i < vector::dims ; i++) + x0[i] = (p.get(i)-domain.getLow(i))*dx[i]; - for (size_t i = 0 ; i < vector::dims ; i++) - base.set_d(i,ip[i][0] - gd.getLocalGridsInfo().get(sub).origin.get(i) - (long int)kernel::np/2 + 1); + // convert into integer + for (size_t i = 0 ; i < vector::dims ; i++) + ip[i][0] = (int)x0[i]; - // convenient grid of points + // convert the global grid position into local grid position + grid_key_dx base; - for (size_t j = 0 ; j < kernel::np-1 ; j++) - { for (size_t i = 0 ; i < vector::dims ; i++) - ip[i][j+1] = (int)ip[i][j]+1; - } + base.set_d(i,ip[i][0] - gd.getLocalGridsInfo().get(sub).origin.get(i) - (long int)kernel::np/2 + 1); - for (size_t i = 0 ; i < vector::dims ; i++) - xp[i] = x0[i] - ip[i][0]; + // convenient grid of points - for (long int j = 0 ; j < kernel::np ; j++) - { - for (size_t i = 0 ; i < vector::dims ; i++) - x[i][j] = - xp[i] + typename vector::stype((long int)j - (long int)kernel::np/2 + 1); - } + for (size_t j = 0 ; j < kernel::np-1 ; j++) + { + for (size_t i = 0 ; i < vector::dims ; i++) + ip[i][j+1] = (int)ip[i][j]+1; + } - for (size_t j = 0 ; j < kernel::np ; j++) - { for (size_t i = 0 ; i < vector::dims ; i++) - a[i][j] = kernel::value(x[i][j],j); - } + xp[i] = x0[i] - ip[i][0]; - grid_sm gs(sz); - grid_key_dx_iterator kit(gs); - - while (kit.isNext()) - { - auto key = kit.get(); + for (long int j = 0 ; j < kernel::np ; j++) + { + for (size_t i = 0 ; i < vector::dims ; i++) + x[i][j] = - xp[i] + typename vector::stype((long int)j - (long int)kernel::np/2 + 1); + } - a_int.template get<0>(key) = 1; + for (size_t j = 0 ; j < kernel::np ; j++) + { + for (size_t i = 0 ; i < vector::dims ; i++) + a[i][j] = kernel::value(x[i][j],j); + } - for (size_t i = 0 ; i < vector::dims ; i++) - a_int.template get<0>(key) *= a[i][key.get(i)]; + calculate_aint::value(sz,a_int,a); - ++kit; - } + grid_dist_lin_dx k_dist_lin; + k_dist_lin.setSub(sub); - grid_key_dx_iterator kit2(gs); - grid_dist_key_dx k_dist; - k_dist.setSub(sub); + size_t k = 0; + auto & gs_info = gd.get_loc_grid(sub).getGrid(); - while (kit2.isNext()) - { - auto key = kit2.get(); + size_t lin_base = gs_info.LinId(base); - for (size_t i = 0 ; i < vector::dims ; i++) - k_dist.getKeyRef().set_d(i,key.get(i) + base.get(i)); + #pragma omp simd + for (size_t i = 0 ; i < openfpm::math::pow(vector::dims,kernel::np) ; i++) + { + size_t lin = offsets.get(sub).ele[k] + lin_base; + k_dist_lin.getKeyRef() = lin; - inte_template::value(gd,vd,k_dist,key_p,a_int,key); + inte_template::value(gd,vd,k_dist_lin,key_p,a_int,k); - ++kit2; + k++; + } } -} - +}; +/*! \brief Main class for interpolation Particle to mest p2m and Mesh to particle m2p + * + * This function is the main class to interpolate from particle to mesh and mesh to particle + * + * \tparam vector type of vector for interpolation + * \tparam grid type of grid for interpolation + * \tparam interpolation kernel + * + */ template class interpolate { + //! Cell list used to convert particles position to sub-domain CellList,shift> geo_cell; + /*! Structure to order boxes by volume + * + * + * + */ struct Box_vol { + //! Box Box bv; + //! Calculated volume size_t vol; - void operator<(const Box_vol & bv) + /*! \brief operator to reorder + * + * \param bv box to compare with + * + * \return true if bv has volume bigger volume + * + */ + bool operator<(const Box_vol & bv) { return vol < bv.vol; } @@ -223,9 +409,53 @@ class interpolate //! grid grid & gd; + //! Type of the calculations + typedef typename vector::stype arr_type; + + /*! \brief It calculate the interpolation stencil offsets + * + * \param array where to store offsets + * \param sz kernel stencil points in each direction + * + */ + void calculate_the_offsets(openfpm::vector> & offsets, size_t (& sz)[vector::dims]) + { + offsets.resize(gd.getN_loc_grid()); + + for (size_t i = 0 ; i < offsets.size() ; i++) + { + auto & loc_grid = gd.get_loc_grid(i); + const grid_sm & gs = loc_grid.getGrid(); + + grid_sm gs2(sz); + grid_key_dx_iterator kit2(gs2); + + size_t k = 0; + + while (kit2.isNext()) + { + auto key = kit2.get(); + + long int lin = gs.LinId(key); + + offsets.get(i).ele[k] = lin; + + ++k; + ++kit2; + } + } + } public: + /*! \brief construct an interpolation object between a grid and a vector + * + * When possible and easy to do we suggest to retain the object + * + * \param vd interpolation vector + * \param gd interpolation grid + * + */ interpolate(vector & vd, grid & gd) :vd(vd),gd(gd) { @@ -266,6 +496,16 @@ public: } }; + /*! \brief Interpolate particles to mesh + * + * Most of the time the particle set and the mesh are the same + * as the one used in the constructor. They also can be different + * as soon as they have the same decomposition + * + * \param gd grid or mesh + * \param vd particle set + * + */ template void p2m(vector & vd, grid & gd) { #ifdef SE_CLASS1 @@ -289,6 +529,16 @@ public: for (size_t i = 0 ; i < vector::dims ; i++) dx[i] = 1.0/gd.spacing(i); + size_t sz[vector::dims]; + + for (size_t i = 0 ; i < vector::dims ; i++) + sz[i] = kernel::np; + + // Precalculate the offset for each sub-sub-domain + openfpm::vector> offsets; + + calculate_the_offsets(offsets,sz); + // point position typename vector::stype xp[vector::dims]; @@ -296,24 +546,30 @@ public: typename vector::stype x[vector::dims][kernel::np]; typename vector::stype a[vector::dims][kernel::np]; - size_t sz[vector::dims]; - - for (size_t i = 0 ; i < vector::dims ; i++) - sz[i] = kernel::np; - - grid_cpu> a_int(sz); - a_int.setMemory(); +// grid_cpu> a_int(sz); +// a_int.setMemory(); + typename vector::stype a_int[openfpm::math::pow(vector::dims,kernel::np)]; auto it = vd.getDomainIterator(); while (it.isNext() == true) { - inte_calc(it,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell); + inte_calc_impl::template inte_calc(it,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets); ++it; } } + /*! \brief Interpolate mesh to particle + * + * Most of the time the particle set and the mesh are the same + * as the one used in the constructor. They also can be different + * as soon as they have the same decomposition + * + * \param gd grid or mesh + * \param vd particle set + * + */ template void m2p(grid & gd, vector & vd) { #ifdef SE_CLASS1 @@ -349,14 +605,20 @@ public: for (size_t i = 0 ; i < vector::dims ; i++) sz[i] = kernel::np; - grid_cpu> a_int(sz); - a_int.setMemory(); + // Precalculate the offset for each sub-sub-domain + openfpm::vector> offsets; + + calculate_the_offsets(offsets,sz); + + // grid_cpu> a_int(sz); + // a_int.setMemory(); + typename vector::stype a_int[openfpm::math::pow(vector::dims,kernel::np)]; auto it = vd.getDomainIterator(); while (it.isNext() == true) { - inte_calc(it,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell); + inte_calc_impl::template inte_calc(it,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets); ++it; } -- GitLab