From e29e2f097ff9ef0fd0e97316ee5ab0e880eb67e8 Mon Sep 17 00:00:00 2001
From: Pietro Incardona <incardon@mpi-cbg.de>
Date: Fri, 16 Feb 2018 21:50:50 +0100
Subject: [PATCH] Adding missing files
---
src/DMatrix/EMatrix.hpp | 188 ++++++++++++++++++
.../vector_dist_operators_apply_kernel.hpp | 10 +
src/interpolation/interpolation.hpp | 150 ++++++++++----
.../interpolation_unit_tests.hpp | 174 ++++++++++------
4 files changed, 426 insertions(+), 96 deletions(-)
create mode 100644 src/DMatrix/EMatrix.hpp
diff --git a/src/DMatrix/EMatrix.hpp b/src/DMatrix/EMatrix.hpp
new file mode 100644
index 00000000..5ec05b19
--- /dev/null
+++ b/src/DMatrix/EMatrix.hpp
@@ -0,0 +1,188 @@
+/*
+ * EMatrix.hpp
+ *
+ * Created on: Feb 12, 2018
+ * Author: i-bird
+ */
+
+#ifndef OPENFPM_DATA_SRC_SPACE_EMATRIX_HPP_
+#define OPENFPM_DATA_SRC_SPACE_EMATRIX_HPP_
+
+#include "config.h"
+
+#ifdef HAVE_EIGEN
+#include <Eigen/Dense>
+#include "memory/ExtPreAlloc.hpp"
+#include "memory/HeapMemory.hpp"
+#include "Packer_Unpacker/Packer_util.hpp"
+#include "Packer_Unpacker/Packer.hpp"
+#include "Packer_Unpacker/Unpacker.hpp"
+
+/*! \brief it is just a wrapper for Eigen matrix compatible for openfpm serialization
+ *
+ *
+ */
+template<typename _Scalar, int _Rows, int _Cols>
+class EMatrix : public Eigen::Matrix<_Scalar,_Rows,_Cols>
+{
+public:
+
+ /*! \brief It calculate the number of byte required to serialize the object
+ *
+ * \tparam prp list of properties
+ *
+ * \param req reference to the total counter required to pack the information
+ *
+ */
+ template<int ... prp> inline void packRequest(size_t & req) const
+ {
+ // Memory required to serialize the Matrix
+ req += sizeof(_Scalar)*this->rows()*this->cols() + 2*sizeof(_Scalar);
+ }
+
+
+ /*! \brief pack a vector selecting the properties to pack
+ *
+ * \param mem preallocated memory where to pack the vector
+ * \param sts pack-stat info
+ *
+ */
+ template<int ... prp> inline void pack(ExtPreAlloc<HeapMemory> & mem, Pack_stat & sts) const
+ {
+ //Pack the number of rows and colums
+ Packer<size_t, HeapMemory>::pack(mem,this->rows(),sts);
+ Packer<size_t, HeapMemory>::pack(mem,this->cols(),sts);
+
+ mem.allocate(sizeof(_Scalar)*this->rows()*this->cols());
+ void * dst = mem.getPointer();
+ void * src = (void *)this->data();
+
+ memcpy(dst,src,sizeof(_Scalar)*this->rows()*this->cols());
+ sts.incReq();
+ }
+
+ /*! \brief unpack a vector
+ *
+ * \param mem preallocated memory from where to unpack the vector
+ * \param ps unpack-stat info
+ */
+ template<int ... prp> inline void unpack(ExtPreAlloc<HeapMemory> & mem, Unpack_stat & ps)
+ {
+ size_t row;
+ size_t col;
+
+ //Pack the number of rows and colums
+ Unpacker<size_t, HeapMemory>::unpack(mem,row,ps);
+ Unpacker<size_t, HeapMemory>::unpack(mem,col,ps);
+
+ this->resize(row,col);
+
+ void * dst = this->data();
+ void * src = mem.getPointerOffset(ps.getOffset());
+
+ memcpy(dst,src,sizeof(_Scalar)*row*col);
+
+ ps.addOffset(sizeof(_Scalar)*row*col);
+ }
+
+};
+
+// Standard typedef from eigen
+typedef EMatrix< std::complex< double >, 2, 2 > EMatrix2cd;
+typedef EMatrix< std::complex< float >, 2, 2 > EMatrix2cf;
+typedef EMatrix< double, 2, 2 > EMatrix2d;
+typedef EMatrix< float, 2, 2 > EMatrix2f;
+typedef EMatrix< int, 2, 2 > EMatrix2i;
+typedef EMatrix< std::complex< double >, 2, Eigen::Dynamic > EMatrix2Xcd;
+typedef EMatrix< std::complex< float >, 2, Eigen::Dynamic > EMatrix2Xcf;
+typedef EMatrix< double, 2, Eigen::Dynamic > EMatrix2Xd;
+typedef EMatrix< float, 2, Eigen::Dynamic > EMatrix2Xf;
+typedef EMatrix< int, 2, Eigen::Dynamic > EMatrix2Xi;
+typedef EMatrix< std::complex< double >, 3, 3 > EMatrix3cd;
+typedef EMatrix< std::complex< float >, 3, 3 > EMatrix3cf;
+typedef EMatrix< double, 3, 3 > EMatrix3d;
+typedef EMatrix< float, 3, 3 > EMatrix3f;
+typedef EMatrix< int, 3, 3 > EMatrix3i;
+typedef EMatrix< std::complex< double >, 3, Eigen::Dynamic > EMatrix3Xcd;
+typedef EMatrix< std::complex< float >, 3, Eigen::Dynamic > EMatrix3Xcf;
+typedef EMatrix< double, 3, Eigen::Dynamic > EMatrix3Xd;
+typedef EMatrix< float, 3, Eigen::Dynamic > EMatrix3Xf;
+typedef EMatrix< int, 3, Eigen::Dynamic > EMatrix3Xi;
+typedef EMatrix< std::complex< double >, 4, 4 > EMatrix4cd;
+typedef EMatrix< std::complex< float >, 4, 4 > EMatrix4cf;
+typedef EMatrix< double, 4, 4 > EMatrix4d;
+typedef EMatrix< float, 4, 4 > EMatrix4f;
+typedef EMatrix< int, 4, 4 > EMatrix4i;
+typedef EMatrix< std::complex< double >, 4, Eigen::Dynamic > EMatrix4Xcd;
+typedef EMatrix< std::complex< float >, 4, Eigen::Dynamic > EMatrix4Xcf;
+typedef EMatrix< double, 4, Eigen::Dynamic > EMatrix4Xd;
+typedef EMatrix< float, 4, Eigen::Dynamic > EMatrix4Xf;
+typedef EMatrix< int, 4, Eigen::Dynamic > EMatrix4Xi;
+typedef EMatrix< std::complex< double >, Eigen::Dynamic, 2 > EMatrixX2cd;
+typedef EMatrix< std::complex< float >, Eigen::Dynamic, 2 > EMatrixX2cf;
+typedef EMatrix< double, Eigen::Dynamic, 2 > EMatrixX2d;
+typedef EMatrix< float, Eigen::Dynamic, 2 > EMatrixX2f;
+typedef EMatrix< int, Eigen::Dynamic, 2 > EMatrixX2i;
+typedef EMatrix< std::complex< double >, Eigen::Dynamic, 3 > EMatrixX3cd;
+typedef EMatrix< std::complex< float >, Eigen::Dynamic, 3 > EMatrixX3cf;
+typedef EMatrix< double, Eigen::Dynamic, 3 > EMatrixX3d;
+typedef EMatrix< float, Eigen::Dynamic, 3 > EMatrixX3f;
+typedef EMatrix< int, Eigen::Dynamic, 3 > EMatrixX3i;
+typedef EMatrix< std::complex< double >, Eigen::Dynamic, 4 > EMatrixX4cd;
+typedef EMatrix< std::complex< float >, Eigen::Dynamic, 4 > EMatrixX4cf;
+typedef EMatrix< double, Eigen::Dynamic, 4 > EMatrixX4d;
+typedef EMatrix< float, Eigen::Dynamic, 4 > EMatrixX4f;
+typedef EMatrix< int, Eigen::Dynamic, 4 > EMatrixX4i;
+typedef EMatrix< std::complex< double >, Eigen::Dynamic, Eigen::Dynamic > EMatrixXcd;
+typedef EMatrix< std::complex< float >, Eigen::Dynamic, Eigen::Dynamic > EMatrixXcf;
+typedef EMatrix< double, Eigen::Dynamic, Eigen::Dynamic > EMatrixXd;
+typedef EMatrix< float, Eigen::Dynamic, Eigen::Dynamic > EMatrixXf;
+typedef EMatrix< int, Eigen::Dynamic, Eigen::Dynamic > EMatrixXi;
+typedef EMatrix< std::complex< double >, 1, 2 > ERowVector2cd;
+typedef EMatrix< std::complex< float >, 1, 2 > ERowVector2cf;
+typedef EMatrix< double, 1, 2 > ERowVector2d;
+typedef EMatrix< float, 1, 2 > ERowVector2f;
+typedef EMatrix< int, 1, 2 > ERowVector2i;
+typedef EMatrix< std::complex< double >, 1, 3 > ERowVector3cd;
+typedef EMatrix< std::complex< float >, 1, 3 > ERowVector3cf;
+typedef EMatrix< double, 1, 3 > ERowVector3d;
+typedef EMatrix< float, 1, 3 > ERowVector3f;
+typedef EMatrix< int, 1, 3 > ERowVector3i;
+typedef EMatrix< std::complex< double >, 1, 4 > ERowVector4cd;
+typedef EMatrix< std::complex< float >, 1, 4 > ERowVector4cf;
+typedef EMatrix< double, 1, 4 > ERowVector4d;
+typedef EMatrix< float, 1, 4 > ERowVector4f;
+typedef EMatrix< int, 1, 4 > ERowVector4i;
+typedef EMatrix< std::complex< double >, 1, Eigen::Dynamic > ERowVectorXcd;
+typedef EMatrix< std::complex< float >, 1, Eigen::Dynamic > ERowVectorXcf;
+typedef EMatrix< double, 1, Eigen::Dynamic > ERowVectorXd;
+typedef EMatrix< float, 1, Eigen::Dynamic > ERowVectorXf;
+typedef EMatrix< int, 1, Eigen::Dynamic > ERowVectorXi;
+typedef EMatrix< std::complex< double >, 2, 1 > EVector2cd;
+typedef EMatrix< std::complex< float >, 2, 1 > EVector2cf;
+typedef EMatrix< double, 2, 1 > EVector2d;
+typedef EMatrix< float, 2, 1 > EVector2f;
+typedef EMatrix< int, 2, 1 > EVector2i;
+typedef EMatrix< std::complex< double >, 3, 1 > EVector3cd;
+typedef EMatrix< std::complex< float >, 3, 1 > EVector3cf;
+typedef EMatrix< double, 3, 1 > EVector3d;
+typedef EMatrix< float, 3, 1 > EVector3f;
+typedef EMatrix< int, 3, 1 > EVector3i;
+typedef EMatrix< std::complex< double >, 4, 1 > EVector4cd;
+typedef EMatrix< std::complex< float >, 4, 1 > EVector4cf;
+typedef EMatrix< double, 4, 1 > EVector4d;
+typedef EMatrix< float, 4, 1 > EVector4f;
+typedef EMatrix< int, 4, 1 > EVector4i;
+typedef EMatrix< std::complex< double >, Eigen::Dynamic, 1 > EVectorXcd;
+typedef EMatrix< std::complex< float >, Eigen::Dynamic, 1 > EVectorXcf;
+typedef EMatrix< double, Eigen::Dynamic, 1 > EVectorXd;
+typedef EMatrix< float, Eigen::Dynamic, 1 > EVectorXf;
+typedef EMatrix< int, Eigen::Dynamic, 1 > EVectorXi;
+
+#else
+
+// There is not EMatrix if we do not have Eigen
+
+#endif
+
+#endif /* OPENFPM_DATA_SRC_SPACE_EMATRIX_HPP_ */
diff --git a/src/Operators/Vector/vector_dist_operators_apply_kernel.hpp b/src/Operators/Vector/vector_dist_operators_apply_kernel.hpp
index 844360ab..80ed67c6 100644
--- a/src/Operators/Vector/vector_dist_operators_apply_kernel.hpp
+++ b/src/Operators/Vector/vector_dist_operators_apply_kernel.hpp
@@ -8,6 +8,16 @@
#ifndef OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_APPLY_KERNEL_HPP_
#define OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_APPLY_KERNEL_HPP_
+//////// SET of small macro to make to define integral easy
+
+#define DEFINE_INTERACTION_3D(name) struct name \
+{\
+\
+ Point<3,double> value(const Point<3,double> & xp, const Point<3,double> xq)\
+ {
+
+#define END_INTERACTION }\
+ };
/*! \brief is_expression check if a type is simple a type or is just an encapsulation of an expression
*
diff --git a/src/interpolation/interpolation.hpp b/src/interpolation/interpolation.hpp
index 8fd0ca06..f1384dbb 100644
--- a/src/interpolation/interpolation.hpp
+++ b/src/interpolation/interpolation.hpp
@@ -11,6 +11,7 @@
#include "NN/Mem_type/MemFast.hpp"
#include "NN/CellList/CellList.hpp"
#include "Grid/grid_dist_key.hpp"
+#include "Vector/vector_dist_key.hpp"
#define INTERPOLATION_ERROR_OBJECT std::runtime_error("Runtime interpolation error");
@@ -407,8 +408,8 @@ struct inte_calc_impl
* kernel stencil for each local grid (sub-domain)
*
*/
- template<unsigned int prp_g, unsigned int prp_v, unsigned int m2p_or_p2m, unsigned int np_a_int, typename iterator, typename grid>
- static inline void inte_calc(iterator & it,
+ template<unsigned int prp_g, unsigned int prp_v, unsigned int m2p_or_p2m, unsigned int np_a_int, typename grid>
+ static inline void inte_calc(const vect_dist_key_dx & key_p,
vector & vd,
const Box<vector::dims,typename vector::stype> & domain,
int (& ip)[vector::dims][kernel::np],
@@ -422,8 +423,6 @@ struct inte_calc_impl
const CellList<vector::dims,typename vector::stype,Mem_fast<>,shift<vector::dims,typename vector::stype>> & geo_cell,
openfpm::vector<agg_arr<openfpm::math::pow(kernel::np,vector::dims)>> & offsets)
{
- auto key_p = it.get();
-
Point<vector::dims,typename vector::stype> p = vd.getPos(key_p);
// On which sub-domain we interpolate the particle
@@ -541,6 +540,18 @@ class interpolate
//! Type of the calculations
typedef typename vector::stype arr_type;
+ //! the offset for each sub-sub-domain
+ openfpm::vector<agg_arr<openfpm::math::pow(kernel::np,vector::dims)>> offsets;
+
+ //! kernel size
+ size_t sz[vector::dims];
+
+ //! grid spacing
+ typename vector::stype dx[vector::dims];
+
+ //! Simulation domain
+ Box<vector::dims,typename vector::stype> domain;
+
/*! \brief It calculate the interpolation stencil offsets
*
* \param offsets array where to store the linearized offset of the
@@ -624,6 +635,18 @@ public:
++g_sub;
}
}
+
+ for (size_t i = 0 ; i < vector::dims ; i++)
+ {sz[i] = kernel::np;}
+
+ calculate_the_offsets(offsets,sz);
+
+ // calculate spacing
+ for (size_t i = 0 ; i < vector::dims ; i++)
+ {dx[i] = 1.0/gd.spacing(i);}
+
+ // simulation domain
+ domain = vd.getDecomposition().getDomain();
};
/*! \brief Interpolate particles to mesh
@@ -651,23 +674,51 @@ public:
#endif
- Box<vector::dims,typename vector::stype> domain = vd.getDecomposition().getDomain();
+ // point position
+ typename vector::stype xp[vector::dims];
- // grid spacing
- typename vector::stype dx[vector::dims];
+ int ip[vector::dims][kernel::np];
+ typename vector::stype x[vector::dims][kernel::np];
+ typename vector::stype a[vector::dims][kernel::np];
- for (size_t i = 0 ; i < vector::dims ; i++)
- dx[i] = 1.0/gd.spacing(i);
+ typename vector::stype a_int[openfpm::math::pow(kernel::np,vector::dims)];
- size_t sz[vector::dims];
+ auto it = vd.getDomainIterator();
- for (size_t i = 0 ; i < vector::dims ; i++)
- sz[i] = kernel::np;
+ while (it.isNext() == true)
+ {
+ auto key_p = it.get();
- // Precalculate the offset for each sub-sub-domain
- openfpm::vector<agg_arr<openfpm::math::pow(kernel::np,vector::dims)>> offsets;
+ inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_p2m,openfpm::math::pow(kernel::np,vector::dims)>(key_p,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
- calculate_the_offsets(offsets,sz);
+ ++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<unsigned int prp_g, unsigned int prp_v> void m2p(grid & gd, vector & vd)
+ {
+#ifdef SE_CLASS1
+
+ if (!vd.getDecomposition().is_equal_ng(gd.getDecomposition()) )
+ {
+ std::cerr << __FILE__ << ":" << __LINE__ << " Error: the distribution of the vector of particles" <<
+ " and the grid is different. In order to interpolate the two data structure must have the" <<
+ " same decomposition" << std::endl;
+
+ ACTION_ON_ERROR(INTERPOLATION_ERROR_OBJECT)
+ }
+
+#endif
// point position
typename vector::stype xp[vector::dims];
@@ -676,19 +727,24 @@ public:
typename vector::stype x[vector::dims][kernel::np];
typename vector::stype a[vector::dims][kernel::np];
+ // grid_cpu<vector::dims,aggregate<typename vector::stype>> a_int(sz);
+ // a_int.setMemory();
typename vector::stype a_int[openfpm::math::pow(kernel::np,vector::dims)];
auto it = vd.getDomainIterator();
while (it.isNext() == true)
{
- inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_p2m,openfpm::math::pow(kernel::np,vector::dims)>(it,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
+ auto key_p = it.get();
+
+ inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_m2p,openfpm::math::pow(kernel::np,vector::dims)>(key_p,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
++it;
}
}
- /*! \brief Interpolate mesh to particle
+
+ /*! \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
@@ -696,9 +752,10 @@ public:
*
* \param gd grid or mesh
* \param vd particle set
+ * \param p particle
*
*/
- template<unsigned int prp_g, unsigned int prp_v> void m2p(grid & gd, vector & vd)
+ template<unsigned int prp_v, unsigned int prp_g> inline void p2m(vector & vd, grid & gd,const vect_dist_key_dx & p)
{
#ifdef SE_CLASS1
@@ -713,14 +770,6 @@ public:
#endif
- Box<vector::dims,typename vector::stype> domain = vd.getDecomposition().getDomain();
-
- // grid spacing
- typename vector::stype dx[vector::dims];
-
- for (size_t i = 0 ; i < vector::dims ; i++)
- dx[i] = 1.0/gd.spacing(i);
-
// point position
typename vector::stype xp[vector::dims];
@@ -728,28 +777,49 @@ public:
typename vector::stype x[vector::dims][kernel::np];
typename vector::stype a[vector::dims][kernel::np];
- size_t sz[vector::dims];
+ typename vector::stype a_int[openfpm::math::pow(kernel::np,vector::dims)];
- for (size_t i = 0 ; i < vector::dims ; i++)
- sz[i] = kernel::np;
+ inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_p2m,openfpm::math::pow(kernel::np,vector::dims)>(p,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
+ }
- // Precalculate the offset for each sub-sub-domain
- openfpm::vector<agg_arr<openfpm::math::pow(kernel::np,vector::dims)>> offsets;
+ /*! \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
+ * \param p particle
+ *
+ */
+ template<unsigned int prp_g, unsigned int prp_v> inline void m2p(grid & gd, vector & vd, const vect_dist_key_dx & p)
+ {
+#ifdef SE_CLASS1
- calculate_the_offsets(offsets,sz);
+ if (!vd.getDecomposition().is_equal_ng(gd.getDecomposition()) )
+ {
+ std::cerr << __FILE__ << ":" << __LINE__ << " Error: the distribution of the vector of particles" <<
+ " and the grid is different. In order to interpolate the two data structure must have the" <<
+ " same decomposition" << std::endl;
+
+ ACTION_ON_ERROR(INTERPOLATION_ERROR_OBJECT)
+ }
+
+#endif
+
+ // point position
+ typename vector::stype xp[vector::dims];
+
+ int ip[vector::dims][kernel::np];
+ typename vector::stype x[vector::dims][kernel::np];
+ typename vector::stype a[vector::dims][kernel::np];
// grid_cpu<vector::dims,aggregate<typename vector::stype>> a_int(sz);
// a_int.setMemory();
typename vector::stype a_int[openfpm::math::pow(kernel::np,vector::dims)];
- auto it = vd.getDomainIterator();
-
- while (it.isNext() == true)
- {
- inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_m2p,openfpm::math::pow(kernel::np,vector::dims)>(it,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
-
- ++it;
- }
+ inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_m2p,openfpm::math::pow(kernel::np,vector::dims)>(p,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
}
};
diff --git a/src/interpolation/interpolation_unit_tests.hpp b/src/interpolation/interpolation_unit_tests.hpp
index 6a73dc8e..8e2fa56e 100644
--- a/src/interpolation/interpolation_unit_tests.hpp
+++ b/src/interpolation/interpolation_unit_tests.hpp
@@ -82,8 +82,63 @@ template<typename vector, unsigned int mom_p> void momenta_vector(vector & vd,ty
}
}
+template<unsigned int dim, typename T, typename grid, typename vector>
+void interp_test(grid & gd, vector & vd, bool single_particle)
+{
+ // Reset the grid
-BOOST_AUTO_TEST_CASE( interpolation_full_test_2D )
+ auto it2 = gd.getDomainGhostIterator();
+
+ while (it2.isNext())
+ {
+ auto key = it2.get();
+
+ gd.template get<0>(key) = 0.0;
+
+ ++it2;
+ }
+
+ interpolate<vector,grid,mp4_kernel<float>> inte(vd,gd);
+
+ if (single_particle == false)
+ {inte.template p2m<0,0>(vd,gd);}
+ else
+ {
+ auto it = vd.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto p = it.get();
+
+ inte.template p2m<0,0>(vd,gd,p);
+
+ ++it;
+ }
+ }
+
+ T mg[dim];
+ T mv[dim];
+
+ momenta_grid<grid,0>(gd,mg);
+ momenta_vector<vector,0>(vd,mv);
+
+ for (size_t i = 0 ; i < dim ; i++)
+ {BOOST_REQUIRE_CLOSE(mg[i],mv[i],0.001);}
+
+ momenta_grid<grid,1>(gd,mg);
+ momenta_vector<vector,1>(vd,mv);
+
+ for (size_t i = 0 ; i < dim ; i++)
+ {BOOST_REQUIRE_CLOSE(mg[i],mv[i],0.001);}
+
+ momenta_grid<grid,2>(gd,mg);
+ momenta_vector<vector,2>(vd,mv);
+
+ for (size_t i = 0 ; i < dim ; i++)
+ {BOOST_REQUIRE_CLOSE(mg[i],mv[i],0.001);}
+}
+
+BOOST_AUTO_TEST_CASE( interpolation_full_single_test_2D )
{
Box<2,float> domain({0.0,0.0},{1.0,1.0});
size_t sz[2] = {64,64};
@@ -93,8 +148,7 @@ BOOST_AUTO_TEST_CASE( interpolation_full_test_2D )
size_t bc_v[2] = {PERIODIC,PERIODIC};
- {
- vector_dist<2,float,aggregate<float>> vd(4096,domain,bc_v,gv);
+ vector_dist<2,float,aggregate<float>> vd(65536,domain,bc_v,gv);
grid_dist_id<2,float,aggregate<float>> gd(vd.getDecomposition(),sz,gg);
// set one particle on vd
@@ -115,46 +169,52 @@ BOOST_AUTO_TEST_CASE( interpolation_full_test_2D )
vd.map();
- // Reset the grid
+ interp_test<2,float>(gd,vd,true);
+}
- auto it2 = gd.getDomainGhostIterator();
- while (it2.isNext())
- {
- auto key = it2.get();
+BOOST_AUTO_TEST_CASE( interpolation_full_test_2D )
+{
+ Box<2,float> domain({0.0,0.0},{1.0,1.0});
+ size_t sz[2] = {64,64};
- gd.template get<0>(key) = 0.0;
+ Ghost<2,long int> gg(2);
+ Ghost<2,float> gv(0.01);
- ++it2;
- }
+ size_t bc_v[2] = {PERIODIC,PERIODIC};
- interpolate<decltype(vd),decltype(gd),mp4_kernel<float>> inte(vd,gd);
+ {
+ vector_dist<2,float,aggregate<float>> vd(4096,domain,bc_v,gv);
+ grid_dist_id<2,float,aggregate<float>> gd(vd.getDecomposition(),sz,gg);
- inte.p2m<0,0>(vd,gd);
+ // set one particle on vd
- float mg[2];
- float mv[2];
+ auto it = vd.getDomainIterator();
- momenta_grid<decltype(gd),0>(gd,mg);
- momenta_vector<decltype(vd),0>(vd,mv);
+ while (it.isNext())
+ {
+ auto p = it.get();
- BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
- BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
+ vd.getPos(p)[0] = (double)rand()/RAND_MAX;
+ vd.getPos(p)[1] = (double)rand()/RAND_MAX;
- momenta_grid<decltype(gd),1>(gd,mg);
- momenta_vector<decltype(vd),1>(vd,mv);
+ vd.getProp<0>(p) = 5.0;
- BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
- BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
+ ++it;
+ }
- momenta_grid<decltype(gd),2>(gd,mg);
- momenta_vector<decltype(vd),2>(vd,mv);
+ vd.map();
- BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
- BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
+ // Reset the grid
+ interp_test<2,float>(gd,vd,false);
+
+ float mg[2];
+ float mv[2];
auto & v_cl = create_vcluster();
+ interpolate<decltype(vd),decltype(gd),mp4_kernel<float>> inte(vd,gd);
+
// We have to do a ghost get before interpolating m2p
// Before doing mesh to particle particle must be arranged
// into a grid like
@@ -248,8 +308,7 @@ BOOST_AUTO_TEST_CASE( interpolation_full_test_2D )
}
}
-
-BOOST_AUTO_TEST_CASE( interpolation_full_test_3D )
+BOOST_AUTO_TEST_CASE( interpolation_full_single_test_3D )
{
Box<3,double> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
size_t sz[3] = {64,64,64};
@@ -259,7 +318,6 @@ BOOST_AUTO_TEST_CASE( interpolation_full_test_3D )
size_t bc_v[3] = {PERIODIC,PERIODIC,PERIODIC};
- {
vector_dist<3,double,aggregate<double>> vd(65536,domain,bc_v,gv);
grid_dist_id<3,double,aggregate<double>> gd(vd.getDecomposition(),sz,gg);
@@ -283,47 +341,51 @@ BOOST_AUTO_TEST_CASE( interpolation_full_test_3D )
vd.map();
// Reset the grid
+ interp_test<3,double>(gd,vd,true);
+}
- auto it2 = gd.getDomainGhostIterator();
+BOOST_AUTO_TEST_CASE( interpolation_full_test_3D )
+{
+ Box<3,double> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
+ size_t sz[3] = {64,64,64};
- while (it2.isNext())
- {
- auto key = it2.get();
+ Ghost<3,long int> gg(2);
+ Ghost<3,double> gv(0.01);
- gd.template get<0>(key) = 0.0;
+ size_t bc_v[3] = {PERIODIC,PERIODIC,PERIODIC};
- ++it2;
- }
+ {
+ vector_dist<3,double,aggregate<double>> vd(65536,domain,bc_v,gv);
+ grid_dist_id<3,double,aggregate<double>> gd(vd.getDecomposition(),sz,gg);
- interpolate<decltype(vd),decltype(gd),mp4_kernel<double>> inte(vd,gd);
+ // set one particle on vd
- inte.p2m<0,0>(vd,gd);
+ auto it = vd.getDomainIterator();
- double mg[3];
- double mv[3];
+ while (it.isNext())
+ {
+ auto p = it.get();
- momenta_grid<decltype(gd),0>(gd,mg);
- momenta_vector<decltype(vd),0>(vd,mv);
+ vd.getPos(p)[0] = (double)rand()/RAND_MAX;
+ vd.getPos(p)[1] = (double)rand()/RAND_MAX;
+ vd.getPos(p)[2] = (double)rand()/RAND_MAX;
- BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
- BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
- BOOST_REQUIRE_CLOSE(mg[2],mv[2],0.001);
+ vd.getProp<0>(p) = 5.0;
- momenta_grid<decltype(gd),1>(gd,mg);
- momenta_vector<decltype(vd),1>(vd,mv);
+ ++it;
+ }
- BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
- BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
- BOOST_REQUIRE_CLOSE(mg[2],mv[2],0.001);
+ vd.map();
- momenta_grid<decltype(gd),2>(gd,mg);
- momenta_vector<decltype(vd),2>(vd,mv);
+ // Reset the grid
- BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
- BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
- BOOST_REQUIRE_CLOSE(mg[2],mv[2],0.001);
+ // Reset the grid
+ interp_test<3,double>(gd,vd,false);
auto & v_cl = create_vcluster();
+ double mg[3];
+ double mv[3];
+ interpolate<decltype(vd),decltype(gd),mp4_kernel<double>> inte(vd,gd);
// We have to do a ghost get before interpolating m2p
// Before doing mesh to particle particle must be arranged
--
GitLab