Small fixes with the new structure

src/FiniteDifference/FDScheme.hpp

@@ -365,7 +365,9 @@ 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
 
@@ -512,7 +514,7 @@ public:
 
 		// copy the vector
 		for (size_t i = 0; i < row_b; i++)
-			B.get(i) = b.get(i);
+			B.insert(i,b.get(i));
 
 		return B;
 	}

src/FiniteDifference/eq.hpp

@@ -11,9 +11,6 @@
 #define EQS_FIELD 0
 #define EQS_POS 1
 
-#define STAGGERED_GRID 1
-#define NORMAL_GRID 0
-
 //#define PERIODIC true
 //#define NON_PERIODIC false
 

src/FiniteDifference/eq_unit_test.hpp

@@ -142,13 +142,13 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
 	//! [lid-driven cavity 2D]
 
 	// Domain, a rectangle
-	Box<2,float> domain({0.0,0.0},{3.0,1.0});
+	Box<2,float> domain({0.0,0.0},{1.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[] = {256,64};
+	long int sz[] = {8,8};
 	size_t szu[2];
 	szu[0] = (size_t)sz[0];
 	szu[1] = (size_t)sz[1];

src/FiniteDifference/util/common.hpp

@@ -8,6 +8,8 @@
 #ifndef OPENFPM_NUMERICS_SRC_UTIL_COMMON_HPP_
 #define OPENFPM_NUMERICS_SRC_UTIL_COMMON_HPP_
 
+#define STAGGERED_GRID 1
+#define NORMAL_GRID 0
 
 template<typename T, typename Sfinae = void>
 struct has_grid_type: std::false_type {};

src/Matrix/SparseMatrix_Eigen.hpp

@@ -21,30 +21,49 @@
  *
  */
 template<typename T>
-struct triplet<T,EIGEN_TRIPLET>
+class triplet<T,EIGEN_TRIPLET>
 {
-	Eigen::Triplet<T,long int> triplet;
+	Eigen::Triplet<T,long int> triplet_ei;
+
+public:
 
 	long int & row()
 	{
-		size_t ptr = (size_t)&triplet.row();
+		size_t ptr = (size_t)&triplet_ei.row();
 
 		return (long int &)*(long int *)ptr;
 	}
 
 	long int & col()
 	{
-		size_t ptr = (size_t)&triplet.col();
+		size_t ptr = (size_t)&triplet_ei.col();
 
 		return (long int &)*(long int *)ptr;
 	}
 
 	T & value()
 	{
-		size_t ptr = (size_t)&triplet.value();
+		size_t ptr = (size_t)&triplet_ei.value();
 
 		return (T &)*(T *)ptr;
 	}
+
+	/*! \brief Constructor from row, colum and value
+	 *
+	 * \param i row
+	 * \param j colum
+	 * \param val value
+	 *
+	 */
+	triplet(long int i, long int j, T val)
+	{
+		row() = i;
+		col() = j;
+		value() = val;
+	}
+
+	// Default constructor
+	triplet()	{};
 };
 
 /* ! \brief Sparse Matrix implementation, that map over Eigen
@@ -73,8 +92,6 @@ private:
 
 	/*! \brief Assemble the matrix
 	 *
-	 * \param trpl Matrix in form of triplets
-	 * \param mat Matrix to assemble
 	 *
 	 */
 	void assemble()
@@ -187,6 +204,16 @@ public:
 	{
 		mat.resize(row,col);
 	}
+
+	/*! \brief Get the row i and the colum j of the Matrix
+	 *
+	 * \return the value of the matrix at row i colum j
+	 *
+	 */
+	T operator()(id_t i, id_t j)
+	{
+		return mat.coeffRef(i,j);
+	}
 };
 
 

src/Solvers/umfpack_solver.hpp

@@ -46,16 +46,27 @@ public:
 	 */
 	template<typename impl> static Vector<double> solve(SparseMatrix<double,int,impl> & A, const Vector<double> & b, size_t opt = UMFPACK_NONE)
 	{
-		Vcluster & v_cl = *global_v_cluster;
+		Vcluster & vcl = *global_v_cluster;
 
 		Vector<double> x;
+
 		// only master processor solve
+		Eigen::UmfPackLU<Eigen::SparseMatrix<double,0,int> > solver;
 
-		if (v_cl.getProcessUnitID() == 0)
-		{
-			Eigen::UmfPackLU<Eigen::SparseMatrix<double,0,int> > solver;
+		// Collect the matrix on master
+		auto mat_ei = A.getMat();
+
+		Eigen::Matrix<double, Eigen::Dynamic, 1> x_ei;
+
+		// Collect the vector on master
+		auto b_ei = b.getVec();
 
-			solver.compute(A.getMat());
+		// Copy b into x, this also copy the information on how to scatter back the information on x
+		x = b;
+
+		if (vcl.getProcessUnitID() == 0)
+		{
+			solver.compute(mat_ei);
 
 			if(solver.info()!=Eigen::Success)
 			{
@@ -64,14 +75,14 @@ public:
 				return x;
 			}
 
-			x.getVec() = solver.solve(b.getVec());
+			x_ei = solver.solve(b_ei);
 
 			if (opt & SOLVER_PRINT_RESIDUAL_NORM_INFINITY)
 			{
-				Vector<double> res;
-				res.getVec() = A.getMat() * x.getVec() - b.getVec();
+				Eigen::Matrix<double, Eigen::Dynamic, 1> res;
+				res = mat_ei * x_ei - b_ei;
 
-				std::cout << "Infinity norm: " << res.getVec().lpNorm<Eigen::Infinity>() << "\n";
+				std::cout << "Infinity norm: " << res.lpNorm<Eigen::Infinity>() << "\n";
 			}
 
 			if (opt & SOLVER_PRINT_DETERMINANT)
@@ -79,9 +90,12 @@ public:
 				std::cout << " Determinant: " << solver.determinant() << "\n";
 			}
 
-			// Vector is only on master, scatter back the information
-			x.scatter();
+			x = x_ei;
 		}
+
+		// Vector is only on master, scatter back the information
+		x.scatter();
+
 		return x;
 	}
 };

src/Vector/Vector.hpp

@@ -8,6 +8,19 @@
 #ifndef OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_HPP_
 #define OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_HPP_
 
+#include <Eigen/Sparse>
+
+/*! \brief It store one row value of a vector
+ *
+ * Given a row, store a value
+ *
+ *
+ */
+template<typename T, unsigned int impl>
+class rval
+{
+};
+
 
 /* ! \brief Sparse Matrix implementation
  *

src/Vector/Vector_eigen.hpp

@@ -13,13 +13,80 @@
 #include <fstream>
 #include "Vector_eigen_util.hpp"
 #include "Grid/staggered_dist_grid_util.hpp"
+#include "Space/Ghost.hpp"
+#include "FiniteDifference/util/common.hpp"
 #include <boost/mpl/vector_c.hpp>
+#include <unordered_map>
 
+#define EIGEN_RVAL 1
+
+/*! \brief It store one row value of a vector
+ *
+ * Given a row, store a value
+ *
+ *
+ */
+template<typename T>
+class rval<T,EIGEN_RVAL>
+{
+	// row
+	long int r;
+
+	// value
+	T val;
+
+public:
+
+	// Get the row
+	long int & row()
+	{
+		return r;
+	}
+
+	// Get the value
+	T & value()
+	{
+		return val;
+	}
+
+	/*! \brief Default constructor
+	 *
+	 */
+	rval()	{}
+
+	/*! \brief Constructor from row, colum and value
+	 *
+	 * \param i row
+	 * \param val value
+	 *
+	 */
+	rval(long int i, T val)
+	{
+		row() = i;
+		value() = val;
+	}
+};
 
 template<typename T>
 class Vector<T,Eigen::Matrix<T, Eigen::Dynamic, 1>>
 {
-	Eigen::Matrix<T, Eigen::Dynamic, 1> v;
+	mutable Eigen::Matrix<T, Eigen::Dynamic, 1> v;
+
+	// row value vector
+	mutable openfpm::vector<rval<T,EIGEN_RVAL>> row_val;
+	mutable openfpm::vector<rval<T,EIGEN_RVAL>> row_val_recv;
+
+	// global to local map
+	mutable std::unordered_map<size_t,size_t> map;
+
+	// invalid
+	T invalid;
+
+	// Processors from where we gather
+	mutable openfpm::vector<size_t> prc;
+
+	//size of each chunk
+	mutable openfpm::vector<size_t> sz;
 
 	/*! \brief Check that the size of the iterators match
 	 *
@@ -169,8 +236,90 @@ class Vector<T,Eigen::Matrix<T, Eigen::Dynamic, 1>>
 		}
 	}
 
+
+	/*! \brief Here we collect the full matrix on master
+	 *
+	 */
+	void collect() const
+	{
+		Vcluster & vcl = *global_v_cluster;
+
+		row_val_recv.clear();
+
+		// here we collect all the triplet in one array on the root node
+		vcl.SGather(row_val,row_val_recv,prc,sz,0);
+
+		if (vcl.getProcessUnitID() != 0)
+			row_val.resize(0);
+		else
+			row_val.swap(row_val_recv);
+
+		build_map();
+	}
+
+	/*! \brief Set the Eigen internal vector
+	 *
+	 *
+	 */
+	void setEigen() const
+	{
+		// set the vector
+
+		for (size_t i = 0 ; i < row_val.size() ; i++)
+			v[row_val.get(i).row()] = row_val.get(i).value();
+	}
+
+	/*! \brief Build the map
+	 *
+	 *
+	 */
+	void build_map() const
+	{
+		map.clear();
+
+		for (size_t i = 0 ; i < row_val.size() ; i++)
+			map[row_val.get(i).row()] = i;
+	}
+
 public:
 
+	/*! \brief Copy the vector
+	 *
+	 * \param v vector to copy
+	 *
+	 */
+	Vector(const Vector<T> & v)
+	{
+		this->operator=(v);
+	}
+
+	/*! \brief Copy the vector
+	 *
+	 * \param v vector to copy
+	 *
+	 */
+	Vector(const Vector<T> && v)
+	{
+		this->operator=(v);
+	}
+
+	/*! \brief Create a vector with n elements
+	 *
+	 * \param n number of elements in the vector
+	 *
+	 */
+	Vector(size_t n)
+	{
+		resize(n);
+	}
+
+	/*! \brief Create a vector with 0 elements
+	 *
+	 */
+	Vector()
+	{
+	}
+
 	/*! \brief Resize the Vector
 	 *
 	 * \param row numbers of row
@@ -188,9 +337,44 @@ public:
 	 * \return reference to the element vector
 	 *
 	 */
-	T & get(size_t i)
+	void insert(size_t i, T val)
 	{
-		return v(i);
+		row_val.add();
+
+		// Map
+		map[i] = row_val.size()-1;
+
+		row_val.last().row() = i;
+		row_val.last().value() = val;
+	}
+
+	/*! \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() )
+		{
+			// Does not exist
+
+			std::cerr << __FILE__ << ":" << __LINE__ << " value does not exist " << std::endl;
+
+			return invalid;
+		}
+		else
+			return row_val.get(it->second).value();
+
+		return invalid;
 	}
 
 	/*! \brief Get the Eigen Vector object
@@ -200,6 +384,9 @@ public:
 	 */
 	const Eigen::Matrix<T, Eigen::Dynamic, 1> & getVec() const
 	{
+		collect();
+		setEigen();
+
 		return v;
 	}
 
@@ -210,6 +397,9 @@ public:
 	 */
 	Eigen::Matrix<T, Eigen::Dynamic, 1> & getVec()
 	{
+		collect();
+		setEigen();
+
 		return v;
 	}
 
@@ -228,30 +418,28 @@ public:
 		}
 	}
 
-	/*! \brief Collect the vector into one processor
-	 *
-	 * Eigen does not have a real parallel vector, so in order to work we have to collect
-	 * the vector in one processor
-	 *
-	 * This function collect the information from all the other processor into one
-	 *
-	 */
-	void collect()
-	{
-
-	}
-
 	/*! \brief Scatter the vector information to the other processors
 	 *
 	 * Eigen does not have a real parallel vector, so in order to work we have to scatter
 	 * the vector from one processor to the other
 	 *
-	 * This function scatter the information to all the other processors
 	 *
 	 */
 	void scatter()
 	{
+		row_val_recv.clear();
+		Vcluster & vcl = *global_v_cluster;
+
+		vcl.SScatter(row_val,row_val_recv,prc,sz,0);
 
+		// if we do not receive anything a previous collect has not been performed
+		// and so nothing is scattered
+		if (row_val_recv.size() != 0)
+		{
+			row_val.clear();
+			row_val.add(row_val_recv);
+			build_map();
+		}
 	}
 
 	/*! \brief Load from file
@@ -270,6 +458,51 @@ public:
 		inputf.close();
 
 	}
+
+	/*! \brief Copy the vector
+	 *
+	 * \param v vector to copy
+	 *
+	 */
+	Vector<T> & operator=(const Vector<T> & v)
+	{
+		prc = v.prc;
+		sz = v.sz;
+		map = v.map;
+		row_val = v.row_val;
+
+		return *this;
+	}
+
+	/*! \brief Copy the vector
+	 *
+	 * \param v vector to copy
+	 *
+	 */
+	Vector<T> & operator=(const Vector<T> && v)
+	{
+		prc = v.prc;
+		sz = v.sz;
+		map = v.map;
+		row_val = v.row_val;
+
+		return *this;
+	}
+
+	/*! \brief Copy the vector (it is used for special purpose)
+	 *
+	 * \warning v MUST contain at least all the elements of the vector
+	 *
+	 * \param v base eigen vector to copy
+	 *
+	 */
+	Vector<T> & operator=(Eigen::Matrix<T, Eigen::Dynamic, 1> & v)
+	{
+		for (size_t i = 0 ; i < row_val.size() ; i++)
+			row_val.get(i).value() = v(row_val.get(i).row());
+
+		return *this;
+	}
 };
 
 

src/Vector/Vector_eigen_util.hpp

@@ -8,6 +8,8 @@
 #ifndef OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_EIGEN_UTIL_HPP_
 #define OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_EIGEN_UTIL_HPP_
 
+#include "Grid/grid_dist_key.hpp"
+#include "Space/Shape/HyperCube.hpp"
 
 /*!	\brief Copy scalar elements
  *

src/main.cpp

@@ -4,6 +4,8 @@
 
 #include "unit_test_init_cleanup.hpp"
 
+#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"