diff --git a/src/Equations/eq_unit_tests.hpp b/src/Equations/eq_unit_tests.hpp
deleted file mode 100644
index ba30b638addcb3ffeac9642512a8a2c5d89b4c63..0000000000000000000000000000000000000000
--- a/src/Equations/eq_unit_tests.hpp
+++ /dev/null
@@ -1,86 +0,0 @@
-/*
- * eq_unit_tests.hpp
- *
- * Created on: Sep 18, 2015
- * Author: i-bird
- */
-
-#ifndef OPENFPM_NUMERICS_SRC_EQUATIONS_EQ_UNIT_TESTS_HPP_
-#define OPENFPM_NUMERICS_SRC_EQUATIONS_EQ_UNIT_TESTS_HPP_
-
-#include "FiniteDifference/FDScheme.hpp"
-
-BOOST_AUTO_TEST_SUITE( eq_test )
-
-struct Sys_p
-{
- static const unsigned int dims = 3;
- typedef float stype;
-};
-
-// define coordinates id
-
-constexpr unsigned int x = 0;
-constexpr unsigned int y = 1;
-
-// define constants
-
-constexpr unsigned int eta = 0;
-constexpr unsigned int zeta = 1;
-constexpr unsigned int gamma = 2;
-constexpr unsigned int nu = 3;
-constexpr unsigned int lambda = 4;
-constexpr unsigned int dmu = 5;
-constexpr unsigned int sigma[2][2] = {{6,7},{8,9}};
-constexpr unsigned int p[] = {10,11};
-constexpr unsigned int del_nu = 12;
-constexpr unsigned int g_ext = 13;
-
-// define variables fields
-
-// velocity
-constexpr unsigned int V[2] = {0,1};
-
-// pressure
-constexpr unsigned int P = 2;
-
-// model Eq1 active Gel
-
-//// First level expression
-
-
-
-/*
- *
- * p_x^2
- *
- */
-//typedef mul<p[x],p[x],Sys_p> px2;
-
-
-/*
- *
- * p_y^2
- *
- */
-//typedef mul<p[y],p[y],Sys_p> py2;
-
-
-/*
- * p_x * p_y
- *
- */
-//typedef mul<p[x],p[y],Sys_p> px_py;
-
-
-
-BOOST_AUTO_TEST_CASE( eq_test_use)
-{
- // create an equation object
-
-
-}
-
-BOOST_AUTO_TEST_SUITE_END()
-
-#endif /* OPENFPM_NUMERICS_SRC_EQUATIONS_EQ_UNIT_TESTS_HPP_ */
diff --git a/src/FiniteDifference/Average.hpp b/src/FiniteDifference/Average.hpp
index 2d578bd592c8540c75cdfb21dfa202d5d885feba..8e92b279e85a5f1071675bb1bae0e0f7da66667a 100644
--- a/src/FiniteDifference/Average.hpp
+++ b/src/FiniteDifference/Average.hpp
@@ -87,7 +87,7 @@ class Avg<d,arg,Sys_eqs,CENTRAL>
* \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
*
*/
- inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
+ inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
{
// if the system is staggered the CENTRAL derivative is equivalent to a forward derivative
if (is_grid_staggered<Sys_eqs>::value())
@@ -170,7 +170,7 @@ class Avg<d,arg,Sys_eqs,FORWARD>
* \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
*
*/
- inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
+ inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
{
long int old_val = kmap.getKeyRef().get(d);
@@ -230,7 +230,7 @@ class Avg<d,arg,Sys_eqs,BACKWARD>
* \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
*
*/
- inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims], std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
+ inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims], std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
{
long int old_val = kmap.getKeyRef().get(d);
kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 1);
diff --git a/src/FiniteDifference/Derivative.hpp b/src/FiniteDifference/Derivative.hpp
index dfd295ba165b26e30306cd8a19c5dafae074a667..b356276fddb2333f154ed75ff53ae70b14584e4a 100644
--- a/src/FiniteDifference/Derivative.hpp
+++ b/src/FiniteDifference/Derivative.hpp
@@ -90,7 +90,7 @@ class D<d,arg,Sys_eqs,CENTRAL>
* \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
*
*/
- inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
+ inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
{
// if the system is staggered the CENTRAL derivative is equivalent to a forward derivative
if (is_grid_staggered<Sys_eqs>::value())
@@ -205,7 +205,7 @@ public:
* \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
*
*/
- static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
+ static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
{
#ifdef SE_CLASS1
if (Sys_eqs::boundary[d] == PERIODIC)
@@ -318,7 +318,7 @@ class D<d,arg,Sys_eqs,FORWARD>
* \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
*
*/
- inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] ,std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
+ inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] ,std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
{
long int old_val = kmap.getKeyRef().get(d);
@@ -374,7 +374,7 @@ class D<d,arg,Sys_eqs,BACKWARD>
* \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
*
*/
- inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims], std::unordered_map<long int,typename Sys_eqs::stype > & cols , typename Sys_eqs::stype coeff)
+ inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims], std::unordered_map<long int,typename Sys_eqs::stype > & cols , typename Sys_eqs::stype coeff)
{
long int old_val = kmap.getKeyRef().get(d);
diff --git a/src/FiniteDifference/FDScheme.hpp b/src/FiniteDifference/FDScheme.hpp
index 64d9ad79fd50fabf7e79210292f92fc202ccbfc1..5e5119ec03a15edfcc7520c41366f4618eb8eea5 100644
--- a/src/FiniteDifference/FDScheme.hpp
+++ b/src/FiniteDifference/FDScheme.hpp
@@ -124,7 +124,7 @@ class FDScheme
public:
// Distributed grid map
- typedef grid_dist_id<Sys_eqs::dims,typename Sys_eqs::stype,scalar<size_t>,typename Sys_eqs::b_grid::decomposition> g_map_type;
+ typedef grid_dist_id<Sys_eqs::dims,typename Sys_eqs::stype,scalar<size_t>,typename Sys_eqs::b_grid::decomposition::extended_type> g_map_type;
typedef Sys_eqs Sys_eqs_typ;
@@ -135,7 +135,8 @@ private:
typedef typename Sys_eqs::SparseMatrix_type::triplet_type triplet;
- openfpm::vector<typename Sys_eqs::stype> b;
+ // Vector b
+ typename Sys_eqs::Vector_type b;
// Domain Grid informations
const grid_sm<Sys_eqs::dims,void> & gs;
@@ -325,13 +326,14 @@ public:
* \param pd Padding, how many points out of boundary are present
* \param domain extension of the domain
* \param gs grid infos where Finite differences work
+ * \param stencil maximum extension of the stencil on each directions
* \param dec Decomposition of the domain
*
*/
- FDScheme(Padding<Sys_eqs::dims> & pd, const Box<Sys_eqs::dims,typename Sys_eqs::stype> & domain, const grid_sm<Sys_eqs::dims,void> & gs, const typename Sys_eqs::b_grid::decomposition & dec)
- :pd(pd),gs(gs),g_map(dec,padding(gs.getSize(),pd),domain,Ghost<Sys_eqs::dims,typename Sys_eqs::stype>(1)),row(0),row_b(0)
+ FDScheme(Padding<Sys_eqs::dims> & pd, const Ghost<Sys_eqs::dims,long int> & stencil, const Box<Sys_eqs::dims,typename Sys_eqs::stype> & domain, const grid_sm<Sys_eqs::dims,void> & gs, const typename Sys_eqs::b_grid & b_g)
+ :pd(pd),gs(gs),g_map(b_g,stencil,pd),row(0),row_b(0)
{
- Vcluster & v_cl = dec.getVC();
+ Vcluster & v_cl = b_g.getDecomposition().getVC();
// Calculate the size of the local domain
size_t sz = g_map.getLocalDomainSize();
@@ -364,10 +366,7 @@ public:
}
// sync the ghost
-
- g_map.write("g_map_before_ghost.vtk");
g_map.template ghost_get<0>();
- g_map.write("g_map_after_ghost.vtk");
// Create a CellDecomposer and calculate the spacing
@@ -399,8 +398,8 @@ public:
template<typename T> void impose(const T & op , typename Sys_eqs::stype num ,long int id ,const long int (& start)[Sys_eqs::dims], const long int (& stop)[Sys_eqs::dims], bool skip_first = false)
{
// add padding to start and stop
- grid_key_dx<Sys_eqs::dims> start_k = grid_key_dx<Sys_eqs::dims>(start) + pd.getKP1();
- grid_key_dx<Sys_eqs::dims> stop_k = grid_key_dx<Sys_eqs::dims>(stop) + pd.getKP1();
+ grid_key_dx<Sys_eqs::dims> start_k = grid_key_dx<Sys_eqs::dims>(start);
+ grid_key_dx<Sys_eqs::dims> stop_k = grid_key_dx<Sys_eqs::dims>(stop);
auto it = g_map.getSubDomainIterator(start_k,stop_k);
@@ -423,6 +422,8 @@ public:
*/
template<typename T> void impose(const T & op , typename Sys_eqs::stype num ,long int id ,grid_dist_iterator_sub<Sys_eqs::dims,typename g_map_type::d_grid> it_d, bool skip_first = false)
{
+ Vcluster & v_cl = *global_v_cluster;
+
openfpm::vector<triplet> & trpl = A.getMatrixTriplets();
auto it = it_d;
@@ -431,19 +432,21 @@ public:
std::unordered_map<long int,float> cols;
// resize b if needed
- b.resize(Sys_eqs::nvar * gs.size());
+ b.resize(Sys_eqs::nvar * g_map.size());
bool is_first = skip_first;
// iterate all the grid points
while (it.isNext())
{
- if (is_first == true)
+ if (is_first == true && v_cl.getProcessUnitID() == 0)
{
++it;
is_first = false;
continue;
}
+ else
+ is_first = false;
// get the position
auto key = it.get();
@@ -463,7 +466,7 @@ public:
// std::cout << "(" << trpl.last().row() << "," << trpl.last().col() << "," << trpl.last().value() << ")" << "\n";
}
- b.get(g_map.template get<0>(key)*Sys_eqs::nvar + id) = num;
+ b(g_map.template get<0>(key)*Sys_eqs::nvar + id) = num;
cols.clear();
// std::cout << "\n";
@@ -491,14 +494,12 @@ public:
#ifdef SE_CLASS1
consistency();
#endif
- A.resize(row,row);
+ A.resize(g_map.size()*Sys_eqs::nvar,g_map.size()*Sys_eqs::nvar);
return A;
}
- typename Sys_eqs::Vector_type B;
-
/*! \brief produce the B vector
*
* \return the vector produced
@@ -510,13 +511,7 @@ public:
consistency();
#endif
- B.resize(row_b);
-
- // copy the vector
- for (size_t i = 0; i < row_b; i++)
- B.insert(i,b.get(i));
-
- return B;
+ return b;
}
};
diff --git a/src/FiniteDifference/FDScheme_unit_tests.hpp b/src/FiniteDifference/FDScheme_unit_tests.hpp
index 9079bdc4ca0a30098cd7c271e17659417bfab4e3..a9ca5fdb1d51384d37d88285c3d75e3379e9a11a 100644
--- a/src/FiniteDifference/FDScheme_unit_tests.hpp
+++ b/src/FiniteDifference/FDScheme_unit_tests.hpp
@@ -14,6 +14,7 @@
#include "util/grid_dist_testing.hpp"
#include "FiniteDifference/Average.hpp"
#include "FiniteDifference/sum.hpp"
+#include "eq.hpp"
constexpr unsigned int x = 0;
constexpr unsigned int y = 1;
diff --git a/src/FiniteDifference/Laplacian.hpp b/src/FiniteDifference/Laplacian.hpp
index d40d91f048bc6545efd79947b63c449ab18e8922..b14f90555403e69482bf18e9e95289ba1077b40a 100644
--- a/src/FiniteDifference/Laplacian.hpp
+++ b/src/FiniteDifference/Laplacian.hpp
@@ -33,7 +33,7 @@ class Lap
* \snippet FDScheme_unit_tests.hpp Laplacian usage
*
*/
- inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
+ inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
{
std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
}
@@ -88,7 +88,7 @@ public:
* \snippet FDScheme_unit_tests.hpp Laplacian usage
*
*/
- inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
+ inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
{
// for each dimension
for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
diff --git a/src/FiniteDifference/eq.hpp b/src/FiniteDifference/eq.hpp
index f4af2d96074dca9ffa1628647a55fd97e7ffeb47..4092517538954e4ce073c2fb301b82fe6bc0a55a 100644
--- a/src/FiniteDifference/eq.hpp
+++ b/src/FiniteDifference/eq.hpp
@@ -81,7 +81,7 @@ struct pos_val
template<unsigned int f, typename Sys_eqs>
class Field
{
- typedef typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type map_grid;
+ typedef typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type::extended_type>::type map_grid;
public:
diff --git a/src/FiniteDifference/eq_unit_test.hpp b/src/FiniteDifference/eq_unit_test.hpp
index 6ebf2ab87d16538e09dc3c9d1bbed524c05dacb5..8aff09c951f0a2856a4e017ee88c3993959d02ce 100644
--- a/src/FiniteDifference/eq_unit_test.hpp
+++ b/src/FiniteDifference/eq_unit_test.hpp
@@ -18,6 +18,7 @@
#include "Vector/Vector.hpp"
#include "Solvers/umfpack_solver.hpp"
#include "data_type/aggregate.hpp"
+#include "FiniteDifference/FDScheme.hpp"
BOOST_AUTO_TEST_SUITE( eq_test_suite )
@@ -139,16 +140,18 @@ typedef Avg<y,v_x,lid_nn,FORWARD> avg_vx_f;
BOOST_AUTO_TEST_CASE(lid_driven_cavity)
{
+ Vcluster & v_cl = *global_v_cluster;
+
//! [lid-driven cavity 2D]
// Domain, a rectangle
- Box<2,float> domain({0.0,0.0},{1.0,1.0});
+ Box<2,float> domain({0.0,0.0},{3.0,1.0});
// Ghost (Not important in this case but required)
Ghost<2,float> g(0.01);
// Grid points on x=256 and y=64
- long int sz[] = {8,8};
+ long int sz[] = {256,64};
size_t szu[2];
szu[0] = (size_t)sz[0];
szu[1] = (size_t)sz[1];
@@ -165,8 +168,11 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
// Distributed grid that store the solution
grid_dist_id<2,float,aggregate<float[2],float>,CartDecomposition<2,float>> g_dist(szu,domain,g);
+ // Ghost stencil
+ Ghost<2,long int> stencil_max(1);
+
// Finite difference scheme
- FDScheme<lid_nn> fd(pd,domain,g_dist.getGridInfo(),g_dist.getDecomposition());
+ FDScheme<lid_nn> fd(pd, stencil_max, domain, g_dist.getGridInfo(), g_dist);
// Here we impose the equation, we start from the incompressibility Eq imposed in the bulk with the
// exception of the first point {0,0} and than we set P = 0 in {0,0}, why we are doing this is again
@@ -218,12 +224,35 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
//! [lid-driven cavity 2D]
- g_dist.write("lid_driven_cavity");
-
-
- // Check that match
- bool test = compare("lid_driven_cavity_grid_0_test.vtk","lid_driven_cavity_grid_0.vtk");
- BOOST_REQUIRE_EQUAL(test,true);
+ g_dist.write("lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()));
+
+ if (v_cl.getProcessUnitID() == 0)
+ {
+ if (v_cl.getProcessingUnits() == 1)
+ {
+ // Check that match
+ bool test = compare("lid_driven_cavity_p1_grid_0_test.vtk","lid_driven_cavity_grid_0.vtk");
+ BOOST_REQUIRE_EQUAL(test,true);
+ }
+ else if (v_cl.getProcessingUnits() == 2)
+ {
+ // Check that match
+ bool test = compare("lid_driven_cavity_p2_grid_0_test.vtk","lid_driven_cavity_p2_grid_0.vtk");
+ BOOST_REQUIRE_EQUAL(test,true);
+ test = compare("lid_driven_cavity_p2_grid_1_test.vtk","lid_driven_cavity_p2_grid_1.vtk");
+ BOOST_REQUIRE_EQUAL(test,true);
+ }
+ else if (v_cl.getProcessingUnits() == 3)
+ {
+ // Check that match
+ bool test = compare("lid_driven_cavity_p3_grid_0_test.vtk","lid_driven_cavity_p3_grid_0.vtk");
+ BOOST_REQUIRE_EQUAL(test,true);
+ test = compare("lid_driven_cavity_p3_grid_1_test.vtk","lid_driven_cavity_p3_grid_1.vtk");
+ BOOST_REQUIRE_EQUAL(test,true);
+ test = compare("lid_driven_cavity_p3_grid_2_test.vtk","lid_driven_cavity_p3_grid_2.vtk");
+ BOOST_REQUIRE_EQUAL(test,true);
+ }
+ }
}
BOOST_AUTO_TEST_SUITE_END()
diff --git a/src/FiniteDifference/eq_unit_test_3d.hpp b/src/FiniteDifference/eq_unit_test_3d.hpp
index 57d8a6eaa0e8815a505ff808c069c3389f5077e1..2bae893de9f363195cef9a4032b8bd97adc2e7e7 100644
--- a/src/FiniteDifference/eq_unit_test_3d.hpp
+++ b/src/FiniteDifference/eq_unit_test_3d.hpp
@@ -130,6 +130,8 @@ typedef Avg<x,v_z,lid_nn_3d,FORWARD> avg_x_vz_f;
BOOST_AUTO_TEST_CASE(lid_driven_cavity)
{
+ Vcluster & v_cl = *global_v_cluster;
+
// Domain
Box<3,float> domain({0.0,0.0,0.0},{3.0,1.0,1.0});
@@ -154,8 +156,11 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
// Initialize openfpm
grid_dist_id<3,float,aggregate<float[3],float>,CartDecomposition<3,float>> g_dist(szu,domain,g);
+ // Ghost stencil
+ Ghost<3,long int> stencil_max(1);
+
// Distributed grid
- FDScheme<lid_nn_3d> fd(pd,domain,g_dist.getGridInfo(),g_dist.getDecomposition());
+ FDScheme<lid_nn_3d> fd(pd,stencil_max,domain,g_dist.getGridInfo(),g_dist);
// start and end of the bulk
@@ -228,11 +233,35 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
// Bring the solution to grid
x.copy<FDScheme<lid_nn_3d>,decltype(g_dist),velocity,pressure>(fd,{0,0},{sz[0]-1,sz[1]-1,sz[2]-1},g_dist);
- g_dist.write("lid_driven_cavity_3d");
-
- // Check that match
- bool test = compare("lid_driven_cavity_3d_grid_0.vtk","lid_driven_cavity_3d_grid_0_test.vtk");
- BOOST_REQUIRE_EQUAL(test,true);
+ g_dist.write("lid_driven_cavity_3d_p" + std::to_string(v_cl.getProcessingUnits()));
+
+ if (v_cl.getProcessUnitID() == 0)
+ {
+ if (v_cl.getProcessingUnits() == 1)
+ {
+ // Check that match
+ bool test = compare("lid_driven_cavity_3d_p1_grid_0_test.vtk","lid_driven_cavity_3d_p1_grid_0.vtk");
+ BOOST_REQUIRE_EQUAL(test,true);
+ }
+ else if (v_cl.getProcessingUnits() == 2)
+ {
+ // Check that match
+ bool test = compare("lid_driven_cavity_p2_grid_0_test.vtk","lid_driven_cavity_p2_grid_0.vtk");
+ BOOST_REQUIRE_EQUAL(test,true);
+ test = compare("lid_driven_cavity_p2_grid_1_test.vtk","lid_driven_cavity_p2_grid_1.vtk");
+ BOOST_REQUIRE_EQUAL(test,true);
+ }
+ else if (v_cl.getProcessingUnits() == 3)
+ {
+ // Check that match
+ bool test = compare("lid_driven_cavity_p3_grid_0_test.vtk","lid_driven_cavity_p3_grid_0.vtk");
+ BOOST_REQUIRE_EQUAL(test,true);
+ test = compare("lid_driven_cavity_p3_grid_1_test.vtk","lid_driven_cavity_p3_grid_1.vtk");
+ BOOST_REQUIRE_EQUAL(test,true);
+ test = compare("lid_driven_cavity_p3_grid_2_test.vtk","lid_driven_cavity_p3_grid_2.vtk");
+ BOOST_REQUIRE_EQUAL(test,true);
+ }
+ }
}
BOOST_AUTO_TEST_SUITE_END()
diff --git a/src/FiniteDifference/mul.hpp b/src/FiniteDifference/mul.hpp
index 15ae6056eab90ef92d46ed0270e11808aefdba8e..eab2e07e644fa1deb607edcd06ea4bddb010d2e1 100644
--- a/src/FiniteDifference/mul.hpp
+++ b/src/FiniteDifference/mul.hpp
@@ -51,7 +51,7 @@ struct const_mul_functor_value
const grid_sm<last::dims,void> & gs;
// grid mapping
- const grid_dist_id<last::dims,typename last::stype,scalar<size_t>,typename last::b_grid::decomposition> & g_map;
+ const grid_dist_id<last::dims,typename last::stype,scalar<size_t>,typename last::b_grid::decomposition::extended_type> & g_map;
// grid position
grid_dist_key_dx<last::dims> & kmap;
@@ -65,7 +65,7 @@ struct const_mul_functor_value
/*! \brief constructor
*
*/
- const_mul_functor_value(const grid_dist_id<last::dims,typename last::stype,scalar<size_t>,typename last::b_grid::decomposition> & g_map, grid_dist_key_dx<last::dims> & kmap, const grid_sm<last::dims,void> & gs, typename last::stype (& spacing)[last::dims], std::unordered_map<long int,typename last::stype> & cols, typename last::stype coeff)
+ const_mul_functor_value(const grid_dist_id<last::dims,typename last::stype,scalar<size_t>,typename last::b_grid::decomposition::extended_type> & g_map, grid_dist_key_dx<last::dims> & kmap, const grid_sm<last::dims,void> & gs, typename last::stype (& spacing)[last::dims], std::unordered_map<long int,typename last::stype> & cols, typename last::stype coeff)
:cols(cols),gs(gs),g_map(g_map),kmap(kmap),coeff(coeff),spacing(spacing)
{};
@@ -113,7 +113,7 @@ struct mul
* \param coeff coefficent (constant in front of the derivative)
*
*/
- inline static void value(const grid_dist_id<Sys_eqs::dims,typename Sys_eqs::stype,scalar<size_t>,typename Sys_eqs::b_grid::decomposition> & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap, const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
+ inline static void value(const grid_dist_id<Sys_eqs::dims,typename Sys_eqs::stype,scalar<size_t>,typename Sys_eqs::b_grid::decomposition::extended_type> & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap, const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
{
const_mul_functor_value<v_expr> mfv(g_map,kmap,gs,spacing,cols,coeff);
diff --git a/src/FiniteDifference/sum.hpp b/src/FiniteDifference/sum.hpp
index 5e485e88f24f99c9b6c9025a896019ad1be1102c..aa80f6b9144d5e150c47e7ee1e757ba189a913c8 100644
--- a/src/FiniteDifference/sum.hpp
+++ b/src/FiniteDifference/sum.hpp
@@ -29,7 +29,7 @@ struct sum_functor_value
const grid_sm<last::dims,void> & gs;
// grid mapping
- const typename stub_or_real<last,last::dims,typename last::stype,typename last::b_grid::decomposition>::type & g_map;
+ const typename stub_or_real<last,last::dims,typename last::stype,typename last::b_grid::decomposition::extended_type>::type & g_map;
// grid position
grid_dist_key_dx<last::dims> & kmap;
@@ -47,7 +47,7 @@ struct sum_functor_value
/*! \brief constructor
*
*/
- sum_functor_value(const typename stub_or_real<last,last::dims,typename last::stype,typename last::b_grid::decomposition>::type & g_map, grid_dist_key_dx<last::dims> & kmap, const grid_sm<last::dims,void> & gs, typename last::stype (& spacing)[last::dims], std::unordered_map<long int,typename last::stype> & cols, typename last::stype coeff)
+ sum_functor_value(const typename stub_or_real<last,last::dims,typename last::stype,typename last::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<last::dims> & kmap, const grid_sm<last::dims,void> & gs, typename last::stype (& spacing)[last::dims], std::unordered_map<long int,typename last::stype> & cols, typename last::stype coeff)
:cols(cols),gs(gs),g_map(g_map),kmap(kmap),key(key),coeff(coeff),spacing(spacing)
{};
@@ -93,7 +93,7 @@ struct sum
* \snippet FDScheme_unit_tests.hpp sum example
*
*/
- inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap, const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
+ inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap, const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
{
// Sum functor
sum_functor_value<v_expr> sm(g_map,kmap,gs,spacing,cols,coeff);
@@ -129,7 +129,7 @@ struct minus
* \snipper FDScheme_unit_tests.hpp minus example
*
*/
- inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap, const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
+ inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap, const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
{
arg::value(g_map,kmap,gs,spacing,cols,-coeff);
}
diff --git a/src/Matrix/SparseMatrix_Eigen.hpp b/src/Matrix/SparseMatrix_Eigen.hpp
index 85bb274cb3256129531908a55b43b3877a77c3a8..87f18d6bbd0d0a450097729fd91031ad3b9d8474 100644
--- a/src/Matrix/SparseMatrix_Eigen.hpp
+++ b/src/Matrix/SparseMatrix_Eigen.hpp
@@ -214,6 +214,103 @@ public:
{
return mat.coeffRef(i,j);
}
+
+ /*! \brief Save this object into file
+ *
+ * \param file filename
+ *
+ * \return true if succed
+ *
+ */
+ bool save(const std::string & file) const
+ {
+ std::vector<size_t> pap_prp;
+
+ Packer<decltype(trpl),HeapMemory>::packRequest(trpl,pap_prp);
+
+ // Calculate how much preallocated memory we need to pack all the objects
+ size_t req = ExtPreAlloc<HeapMemory>::calculateMem(pap_prp);
+
+ // allocate the memory
+ HeapMemory pmem;
+ pmem.allocate(req);
+ ExtPreAlloc<HeapMemory> mem(pap_prp,pmem);
+
+ //Packing
+
+ Pack_stat sts;
+ Packer<decltype(trpl),HeapMemory>::pack(mem,trpl,sts);
+
+ // Save into a binary file
+ std::ofstream dump (file, std::ios::out | std::ios::binary);
+ if (dump.is_open() == false)
+ return false;
+ dump.write ((const char *)pmem.getPointer(), pmem.size());
+
+ return true;
+ }
+
+ /*! \brief Load this object from file
+ *
+ * \param file filename
+ *
+ * \return true if succed
+ *
+ */
+ bool load(const std::string & file)
+ {
+ std::ifstream fs (file, std::ios::in | std::ios::binary | std::ios::ate );
+ if (fs.is_open() == false)
+ return false;
+
+ // take the size of the file
+ size_t sz = fs.tellg();
+
+ fs.close();
+
+ // reopen the file without ios::ate to read
+ std::ifstream input (file, std::ios::in | std::ios::binary );
+ if (input.is_open() == false)
+ return false;
+
+ // Create the HeapMemory and the ExtPreAlloc memory
+ std::vector<size_t> pap_prp;
+ pap_prp.push_back(sz);
+ HeapMemory pmem;
+ ExtPreAlloc<HeapMemory> mem(pap_prp,pmem);
+
+ // read
+ input.read((char *)pmem.getPointer(), sz);
+
+ //close the file
+ input.close();
+
+ //Unpacking
+ Unpack_stat ps;
+
+ Unpacker<decltype(trpl),HeapMemory>::unpack(mem,trpl,ps);
+
+ return true;
+ }
+
+ /*! \brief Get the value from triplet
+ *
+ * \warning It is extremly slow because it do a full search across the triplets elements
+ *
+ * \param r row
+ * \param c colum
+ *
+ */
+ T getValue(size_t r, size_t c)
+ {
+ for (size_t i = 0 ; i < trpl.size() ; i++)
+ {
+ if (r == (size_t)trpl.get(i).row() && c == (size_t)trpl.get(i).col())
+ return trpl.get(i).value();
+ }
+
+ return 0;
+ }
};
diff --git a/src/Vector/Vector_eigen.hpp b/src/Vector/Vector_eigen.hpp
index 4ef615b6c7eb9c825515ef6ea3523038d0f11078..57e1c690bd267d10d5706f0975815e808096ef38 100644
--- a/src/Vector/Vector_eigen.hpp
+++ b/src/Vector/Vector_eigen.hpp
@@ -12,7 +12,7 @@
#include "util/mul_array_extents.hpp"
#include <fstream>
#include "Vector_eigen_util.hpp"
-#include "Grid/staggered_dist_grid_util.hpp"
+#include "Grid/staggered_dist_grid.hpp"
#include "Space/Ghost.hpp"
#include "FiniteDifference/util/common.hpp"
#include <boost/mpl/vector_c.hpp>
@@ -88,146 +88,45 @@ class Vector<T,Eigen::Matrix<T, Eigen::Dynamic, 1>>
//size of each chunk
mutable openfpm::vector<size_t> sz;
- /*! \brief Check that the size of the iterators match
- *
- * It check the the boxes that the sub iterator defines has same dimensions, for example
- * if the first sub-iterator, iterate from (1,1) to (5,3) and the second from (2,2) to (6,4)
- * they match (2,2) to (4,6) they do not match
- *
- * \tparam Grid_map type of the map grid
- * \tparam Grid_dst type of the destination grid
- *
- * \param it1 Iterator1
- * \param it2 Iterator2
- *
- * \return true if they match
- *
- */
- template<typename Eqs_sys, typename it1_type, typename it2_type> bool checkIterator(const it1_type & it1, const it2_type & it2)
- {
-#ifdef SE_CLASS1
-
- grid_key_dx<Eqs_sys::dims> it1_k = it1.getStop() - it1.getStart();
- grid_key_dx<Eqs_sys::dims> it2_k = it2.getStop() - it2.getStart();
-
- for (size_t i = 0 ; i < Eqs_sys::dims ; i++)
- {
- if (it1_k.get(i) != it2_k.get(i))
- {
- std::cerr << __FILE__ << ":" << __LINE__ << " error src iterator and destination iterator does not match in size\n";
- return false;
- }
- }
-
- return true;
-#else
-
- return true;
-
-#endif
- }
-
-
/*! \brief Copy in a staggered grid
*
*
*/
- template<typename scheme, typename Grid_dst ,unsigned int ... pos> void copy_staggered_to_staggered(scheme & sc,const long int (& start)[scheme::Sys_eqs_typ::dims], const long int (& stop)[scheme::Sys_eqs_typ::dims], Grid_dst & g_dst)
+ template<typename scheme, typename Grid_dst ,unsigned int ... pos> void copy_staggered_to_staggered(scheme & sc, Grid_dst & g_dst)
{
// get the map
const auto & g_map = sc.getMap();
- // get the grid padding
- const Padding<scheme::Sys_eqs_typ::dims> & pd = sc.getPadding();
-
- // shift the start and stop by the padding
- grid_key_dx<scheme::Sys_eqs_typ::dims> start_k = grid_key_dx<scheme::Sys_eqs_typ::dims>(start) + pd.getKP1();
- grid_key_dx<scheme::Sys_eqs_typ::dims> stop_k = grid_key_dx<scheme::Sys_eqs_typ::dims>(stop) + pd.getKP1();
-
- // sub-grid iterator over the grid map
- auto g_map_it = g_map.getSubDomainIterator(start_k,stop_k);
-
- // Iterator over the destination grid
- auto g_dst_it = g_dst.getDomainIterator();
-
- // Check that the 2 iterator has the same size
- checkIterator<typename scheme::Sys_eqs_typ,decltype(g_map_it),decltype(g_dst_it)>(g_map_it,g_dst_it);
+ // check that g_dst is staggered
+ if (is_grid_staggered<Grid_dst>::value() == false)
+ std::cerr << __FILE__ << ":" << __LINE__ << " The destination grid must be staggered " << std::endl;
- while (g_map_it.isNext() == true)
- {
- typedef typename to_boost_vmpl<pos...>::type vid;
- typedef boost::mpl::size<vid> v_size;
-
- auto key_src = g_map_it.get();
- size_t lin_id = g_map.template get<0>(key_src);
-
- // destination point
- auto key_dst = g_dst_it.get();
-
- // Transform this id into an id for the Eigen vector
-
- copy_ele<typename scheme::Sys_eqs_typ,Grid_dst,Eigen::Matrix<T, Eigen::Dynamic, 1>> cp(key_dst,g_dst,v,lin_id,g_map.size());
-
- boost::mpl::for_each_ref<boost::mpl::range_c<int,0,v_size::value>>(cp);
-
- ++g_map_it;
- ++g_dst_it;
- }
- }
-
- /*! \brief Copy in a normal grid
- *
- *
- */
- template<typename scheme, typename Grid_dst ,unsigned int ... pos> void copy_staggered_to_normal(scheme & sc,const long int (& start)[scheme::Sys_eqs_typ::dims], const long int (& stop)[scheme::Sys_eqs_typ::dims], Grid_dst & g_dst)
- {
- // get the map
- const auto & g_map = sc.getMap();
-
- // get the grid padding
- const Padding<scheme::Sys_eqs_typ::dims> & pd = sc.getPadding();
-
- // set the staggered position of the property
- openfpm::vector<comb<scheme::Sys_eqs_typ::dims>> stag_pos[sizeof...(pos)];
-
- // interpolation points for each property
- openfpm::vector<std::vector<comb<scheme::Sys_eqs_typ::dims>>> interp_pos[sizeof...(pos)];
-
- // fill the staggered position vector for each property
- stag_set_position<scheme::Sys_eqs_typ::dims,typename Grid_dst::value_type::type> ssp(stag_pos);
-
- typedef boost::mpl::vector_c<unsigned int,pos ... > v_pos_type;
- boost::mpl::for_each_ref<v_pos_type>(ssp);
-
- interp_points<scheme::Sys_eqs_typ::dims,v_pos_type,typename Grid_dst::value_type::type> itp(interp_pos,stag_pos);
- boost::mpl::for_each_ref<v_pos_type>(itp);
+#ifdef SE_CLASS1
- // shift the start and stop by the padding
- grid_key_dx<scheme::Sys_eqs_typ::dims> start_k = grid_key_dx<scheme::Sys_eqs_typ::dims>(start) + pd.getKP1();
- grid_key_dx<scheme::Sys_eqs_typ::dims> stop_k = grid_key_dx<scheme::Sys_eqs_typ::dims>(stop) + pd.getKP1();
+ if (g_map.getLocalDomainSize() != g_dst.getLocalDomainSize())
+ std::cerr << __FILE__ << ":" << __LINE__ << " The staggered and destination grid in size does not match " << std::endl;
+#endif
// sub-grid iterator over the grid map
- auto g_map_it = g_map.getSubDomainIterator(start_k,stop_k);
+ auto g_map_it = g_map.getDomainIterator();
// Iterator over the destination grid
auto g_dst_it = g_dst.getDomainIterator();
- // Check that the 2 iterator has the same size
- checkIterator<typename scheme::Sys_eqs_typ,decltype(g_map_it),decltype(g_dst_it)>(g_map_it,g_dst_it);
-
while (g_map_it.isNext() == true)
{
typedef typename to_boost_vmpl<pos...>::type vid;
typedef boost::mpl::size<vid> v_size;
auto key_src = g_map_it.get();
+ size_t lin_id = g_map.template get<0>(key_src);
// destination point
auto key_dst = g_dst_it.get();
// Transform this id into an id for the Eigen vector
- interp_ele<scheme,Grid_dst,Eigen::Matrix<T, Eigen::Dynamic, 1>,T,sizeof...(pos)> cp(key_dst,g_dst,v,key_src,g_map,interp_pos);
+ copy_ele<typename scheme::Sys_eqs_typ,Grid_dst,typename std::remove_reference<decltype(*this)>::type> cp(key_dst,g_dst,*this,lin_id,g_map.size());
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,v_size::value>>(cp);
@@ -333,8 +232,7 @@ public:
/*! \brief Return a reference to the vector element
*
* \param i element
- *
- * \return reference to the element vector
+ * \param val value
*
*/
void insert(size_t i, T val)
@@ -348,6 +246,42 @@ public:
row_val.last().value() = val;
}
+ /*! \brief Return a reference to the vector element
+ *
+ * \param i element
+ *
+ * \return reference to the element vector
+ *
+ */
+ inline T & insert(size_t i)
+ {
+ row_val.add();
+
+ // Map
+ map[i] = row_val.size()-1;
+
+ row_val.last().row() = i;
+ return row_val.last().value();
+ }
+
+ /*! \brief Return a reference to the vector element
+ *
+ * \param i element
+ *
+ * \return reference to the element vector
+ *
+ */
+ inline const T & insert(size_t i) const
+ {
+ row_val.add();
+
+ // Map
+ map[i] = row_val.size()-1;
+
+ row_val.last().row() = i;
+ return row_val.last().value();
+ }
+
/*! \brief Return a reference to the vector element
*
* \warning The element must exist
@@ -357,24 +291,37 @@ public:
* \return reference to the element vector
*
*/
- T & operator()(size_t i)
+ const T & operator()(size_t i) const
{
// Search if exist
std::unordered_map<size_t,size_t>::iterator it = map.find(i);
- if ( it == map.end() )
- {
- // Does not exist
+ if ( it != map.end() )
+ return row_val.get(it->second).value();
- std::cerr << __FILE__ << ":" << __LINE__ << " value does not exist " << std::endl;
+ return insert(i);
+ }
- return invalid;
- }
- else
+ /*! \brief Return a reference to the vector element
+ *
+ * \warning The element must exist
+ *
+ * \param i element
+ *
+ * \return reference to the element vector
+ *
+ */
+ T & operator()(size_t i)
+ {
+ // Search if exist
+
+ std::unordered_map<size_t,size_t>::iterator it = map.find(i);
+
+ if ( it != map.end() )
return row_val.get(it->second).value();
- return invalid;
+ return insert(i);
}
/*! \brief Get the Eigen Vector object
@@ -405,6 +352,17 @@ public:
/*! \brief Copy the vector into the grid
*
+ * ## Copy from the vector to the destination grid
+ * \snippet eq_unit_tests.hpp
+ *
+ * \tparam scheme Discretization scheme
+ * \tparam Grid_dst type of the target grid
+ * \tparam pos target properties
+ *
+ * \param scheme Discretization scheme
+ * \param start point
+ * \param stop point
+ * \param g_dst Destination grid
*
*/
template<typename scheme, typename Grid_dst ,unsigned int ... pos> void copy(scheme & sc,const long int (& start)[scheme::Sys_eqs_typ::dims], const long int (& stop)[scheme::Sys_eqs_typ::dims], Grid_dst & g_dst)
@@ -412,9 +370,19 @@ public:
if (is_grid_staggered<typename scheme::Sys_eqs_typ>::value())
{
if (g_dst.is_staggered() == true)
- copy_staggered_to_staggered<scheme,Grid_dst,pos...>(sc,start,stop,g_dst);
+ copy_staggered_to_staggered<scheme,Grid_dst,pos...>(sc,g_dst);
else
- copy_staggered_to_normal<scheme,Grid_dst,pos...>(sc,start,stop,g_dst);
+ {
+ // Create a temporal staggered grid and copy the data there
+ auto & g_map = sc.getMap();
+ staggered_grid_dist<Grid_dst::dims,typename Grid_dst::stype,typename Grid_dst::value_type,typename Grid_dst::decomposition::extended_type, typename Grid_dst::memory_type, typename Grid_dst::device_grid_type> stg(g_dst,g_map.getDecomposition().getGhost(),sc.getPadding());
+ stg.setDefaultStagPosition();
+ copy_staggered_to_staggered<scheme,decltype(stg),pos...>(sc,stg);
+
+ // sync the ghost and interpolate to the normal grid
+ stg.template ghost_get<pos...>();
+ stg.template to_normal<Grid_dst,pos...>(g_dst,sc.getPadding(),start,stop);
+ }
}
}
diff --git a/src/Vector/Vector_eigen_util.hpp b/src/Vector/Vector_eigen_util.hpp
index 221a8600be5e5e4bb819c62ff85b9c2bb294b67e..aebff63faa823cf7cbbafa30068f28b1c2b08218 100644
--- a/src/Vector/Vector_eigen_util.hpp
+++ b/src/Vector/Vector_eigen_util.hpp
@@ -49,145 +49,23 @@ struct copy_ele_sca_array<copy_type,T,Ev,Eqs_sys,1>
}
};
-/*! \brief Add scalar elements
- *
- * \tparam copy_type Type that should be copied
- * \tparam T property id to copy
- * \tparam Ev Type of source the Vector
- * \tparam sa dimensionality of the array 0 is a scalar
- *
- */
-template<typename copy_type, typename T, typename Ev, typename scheme, int sa>
-struct interp_ele_sca_array
-{
- template<typename Grid> inline static void interp(Grid & grid_dst, const grid_dist_key_dx<scheme::Sys_eqs_typ::dims> & key, const Ev & x, grid_dist_key_dx<scheme::Sys_eqs_typ::dims> key_src, const openfpm::vector<std::vector<comb<scheme::Sys_eqs_typ::dims>>> & interp_pos, const typename scheme::g_map_type & g_map,size_t base_id)
- {
- copy_type division = 0.0;
-
- for (size_t i = 0 ; i < interp_pos.get(0).size() ; i++)
- {
- auto key_m = key_src.move(interp_pos.get(0)[i]);
- size_t lin_id = g_map.template get<0>(key_m);
-
- grid_dst.template get<T::value>(key) += x(lin_id * scheme::Sys_eqs_typ::nvar + base_id);
-
- division += 1.0;
- }
- grid_dst.template get<T::value>(key) /= division;
- }
-};
-
-/*! \brief Add 1D array elements
- *
- * spacialization in case of 1D array
- *
- * \tparam copy_type Type that should be copied
- * \tparam T property id to copy
- * \tparam Ev Type of source the Vector
- *
- */
-template<typename copy_type, typename T, typename Ev, typename scheme>
-struct interp_ele_sca_array<copy_type,T,Ev,scheme,1>
-{
- template<typename Grid> inline static void interp(Grid & grid_dst, const grid_dist_key_dx<scheme::Sys_eqs_typ::dims> & key, const Ev & x, grid_dist_key_dx<scheme::Sys_eqs_typ::dims> key_src, const openfpm::vector<std::vector<comb<scheme::Sys_eqs_typ::dims>>> & interp_pos, const typename scheme::g_map_type & g_map,size_t base_id)
- {
- typename std::remove_all_extents<copy_type>::type division;
-
- for (size_t j = 0 ; j < std::extent<copy_type>::value ; j++)
- {
- division = 0.0;
- for (size_t i = 0 ; i < interp_pos.get(j).size() ; i++)
- {
- auto key_m = key_src.move(interp_pos.get(j)[i]);
- size_t lin_id = g_map.template get<0>(key_m);
-
- grid_dst.template get<T::value>(key)[j] += x(lin_id * scheme::Sys_eqs_typ::nvar + base_id);
-
- division += 1.0;
- }
- grid_dst.template get<T::value>(key)[j] /= division;
- base_id++;
- }
- }
-};
-
/*! \brief this class is a functor for "for_each" algorithm
*
* This class is a functor for "for_each" algorithm. For each
* element of the boost::vector the operator() is called.
- * Is mainly used to calculate the interpolation points for each
- * property in a staggered grid
+ * Is mainly used to copy from the Vector to the grid target
+ * \note properties can be scalars or arrays of C++ primitives
*
- * \tparam dim Dimensionality
- * \tparam v_prp_id vector of properties id
- * \tparam v_prp_type vector with the properties
- *
- */
-template<unsigned int dim, typename v_prp_id, typename v_prp_type>
-struct interp_points
-{
- // number of properties we are processing
- typedef boost::mpl::size<v_prp_id> v_size;
-
- // interpolation points for each property
- openfpm::vector<std::vector<comb<dim>>> (& interp_pts)[v_size::value];
-
- // staggered position for each property
- const openfpm::vector<comb<dim>> (&stag_pos)[v_size::value];
-
- /*! \brief constructor
- *
- * It define the copy parameters.
- *
- * \param inter_pts array that for each property contain the interpolation points for each components
- * \param staggered position for each property and components
- *
- */
- inline interp_points(openfpm::vector<std::vector<comb<dim>>> (& interp_pts)[v_size::value],const openfpm::vector<comb<dim>> (&stag_pos)[v_size::value])
- :interp_pts(interp_pts),stag_pos(stag_pos){};
-
- //! It call the copy function for each property
- template<typename T>
- inline void operator()(T& t)
- {
- // This is the type of the object we have to copy
- typedef typename boost::mpl::at_c<v_prp_type,T::value>::type prp_type;
-
- interp_pts[T::value].resize(stag_pos[T::value].size());
-
- for (size_t i = 0 ; i < stag_pos[T::value].size() ; i++)
- {
- // Create the interpolation points
- interp_pts[T::value].get(i) = SubHyperCube<dim,dim - std::rank<prp_type>::value>::getCombinations_R(stag_pos[T::value].get(i),0);
-
- // interp_point are -1,0,1, map the -1 to 0 and 1 to -1
- for (size_t j = 0 ; j < interp_pts[T::value].get(i).size() ; j++)
- {
- for (size_t k = 0 ; k < dim ; k++)
- interp_pts[T::value].get(i)[j].getComb()[k] = - ((interp_pts[T::value].get(i)[j].getComb()[k] == -1)?0:interp_pts[T::value].get(i)[j].getComb()[k]);
- }
- }
- }
-};
-
-/*! \brief this class is a functor for "for_each" algorithm
- *
- * This class is a functor for "for_each" algorithm. For each
- * element of the boost::vector the operator() is called.
- * Is mainly used to copy from the Vector_eigen to the grid target
- * in a generic way for a generic object T
- * with variables number of properties scalar or array of C++ primitives
- *
- * \tparam dim Dimensionality
- * \tparam S type of grid
- * \tparam Ev Eigen vector type
+ * \tparam Eqs_sys System of equation information
+ * \tparam S type of destination grid
+ * \tparam Ev type of vector
*
*/
template<typename Eqs_sys, typename S, typename Ev>
struct copy_ele
{
- //! destination element inside the grid
+ //! destination grid element
const grid_dist_key_dx<Eqs_sys::dims> key;
//! destination grid
@@ -196,10 +74,9 @@ struct copy_ele
//! source element inside the Eigen vector
size_t lin_id;
- //! property id inside the Eigen vector, each property processed increase such id
+ //! counter
size_t prp_id;
- //! number of the elements in the Eigen vector divided the number of variables
//! It is basically the number of grid points the Eigen vector has inside
size_t gs_size;
@@ -210,8 +87,8 @@ struct copy_ele
*
* It define the copy parameters.
*
- * \param key which element we are modifying
- * \param grid_dst grid we are updating
+ * \param key destination position
+ * \param grid_dst grid destination
* \param v Source Eigen vector
*
*/
@@ -243,82 +120,5 @@ struct copy_ele
}
};
-/*! \brief this class is a functor for "for_each" algorithm
- *
- * This class is a functor for "for_each" algorithm. For each
- * element of the boost::vector the operator() is called.
- * Is mainly used to interpolate from the staggered Vector_eigen to the grid target
- * in a generic way for a generic object T
- * with variables number of properties scalar or array of C++ primitives
- *
- * \tparam scheme Discretization scheme
- * \tparam S type of grid
- * \tparam Ev Eigen vector type
- *
- */
-template<typename scheme, typename S, typename Ev, typename base_type, unsigned int nst_pos>
-struct interp_ele
-{
- //! destination element inside the grid
- const grid_dist_key_dx<scheme::Sys_eqs_typ::dims> key;
-
- //! destination grid
- S & grid_dst;
-
- //! source element inside the Eigen vector
- grid_dist_key_dx<scheme::Sys_eqs_typ::dims> key_src;
-
- //! g_map
- const typename scheme::g_map_type & g_map;
-
- //! For each properties [] for each components (openfpm::vector) interpolants points positions (std::vector<comb>)
- openfpm::vector<std::vector<comb<scheme::Sys_eqs_typ::dims>>> (&interp_pos)[nst_pos];
-
- //! property id inside the Eigen vector, each property processed increase such id
- size_t prp_id;
-
- //! source Eigen vector
- const Ev & x;
-
- //! Value type of the vector
- base_type division;
-
- /*! \brief constructor
- *
- * It define the interpolation parameters.
- *
- * \param key which element we are modifying
- * \param grid_dst grid we are updating
- * \param v Source Eigen vector
- *
- */
- inline interp_ele(const grid_dist_key_dx<scheme::Sys_eqs_typ::dims> & key, S & grid_dst, const Ev & x, const grid_dist_key_dx<scheme::Sys_eqs_typ::dims> & key_src, const typename scheme::g_map_type & g_map, openfpm::vector<std::vector<comb<scheme::Sys_eqs_typ::dims>>> (&interp_pos)[nst_pos])
- :key(key),grid_dst(grid_dst),key_src(key_src),g_map(g_map),interp_pos(interp_pos),prp_id(0),x(x){};
-
-
-#ifdef SE_CLASS1
- /*! \brief Constructor
- *
- * Calling this constructor produce an error. This class store the reference of the object,
- * this mean that the object passed must not be a temporal object
- *
- */
- inline interp_ele(grid_dist_key_dx<scheme::Sys_eqs_typ::dims> & key, S & grid_dst, Ev && x)
- :key(key),grid_dst(grid_dst),x(x),interp_pos(interp_pos)
- {std::cerr << "Error: " <<__FILE__ << ":" << __LINE__ << " Passing a temporal object";};
-#endif
-
- //! It call the copy function for each property
- template<typename T>
- inline void operator()(T& t)
- {
- // This is the type of the object we have to copy
- typedef typename boost::mpl::at_c<typename S::value_type::type,T::value>::type copy_type;
-
- interp_ele_sca_array<copy_type,T,Ev,scheme,std::rank<copy_type>::value>::interp(grid_dst,key,x,key_src,interp_pos[T::value],g_map,prp_id);
-
- prp_id += array_extents<copy_type>::mul();
- }
-};
#endif /* OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_EIGEN_UTIL_HPP_ */
diff --git a/src/main.cpp b/src/main.cpp
index 8a63fcf886529ad092a5e275a6a60affdd764194..5273f673c29eed21fd51e5c101c16daddff2e759 100755
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -6,7 +6,6 @@
#include "Vector/Vector_unit_tests.hpp"
#include "Matrix/SparseMatrix_unit_tests.hpp"
-#include "Equations/eq_unit_tests.hpp"
#include "FiniteDifference/FDScheme_unit_tests.hpp"
#include "FiniteDifference/util/common_test.hpp"
#include "FiniteDifference/eq_unit_test.hpp"
diff --git a/src/util/grid_dist_testing.hpp b/src/util/grid_dist_testing.hpp
index 39a87e67fd45dcc1991e1842a341653cd0415a31..f49359d6114876ad9206466721ddf4216f59da33 100644
--- a/src/util/grid_dist_testing.hpp
+++ b/src/util/grid_dist_testing.hpp
@@ -8,6 +8,8 @@
#ifndef OPENFPM_NUMERICS_SRC_UTIL_GRID_DIST_TESTING_HPP_
#define OPENFPM_NUMERICS_SRC_UTIL_GRID_DIST_TESTING_HPP_
+#include "data_type/scalar.hpp"
+
template<unsigned int dim>
class grid_dist_testing
{