Vector_eigen.hpp 6.74 KB
Newer Older
1 2 3 4 5 6 7 8 9 10
/*
 * Vector_eigen.hpp
 *
 *  Created on: Nov 27, 2015
 *      Author: i-bird
 */

#ifndef OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_EIGEN_HPP_
#define OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_EIGEN_HPP_

11 12 13 14 15 16
#include <type_traits>
#include "util/mul_array_extents.hpp"
#include <fstream>
#include "Vector_eigen_util.hpp"
#include "Grid/staggered_dist_grid_util.hpp"
#include <boost/mpl/vector_c.hpp>
17 18


19 20
template<typename T>
class Vector<T,Eigen::Matrix<T, Eigen::Dynamic, 1>>
21
{
22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
	Eigen::Matrix<T, Eigen::Dynamic, 1> v;

	/*! \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)
40
	{
41 42 43 44 45 46
#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++)
47
		{
48 49 50 51 52
			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;
			}
53 54
		}

55 56
		return true;
#else
57

58
		return true;
59

60 61
#endif
	}
62 63


64
	/*! \brief Copy in a staggered grid
65 66 67
	 *
	 *
	 */
68 69 70 71
	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)
	{
		// get the map
		const auto & g_map = sc.getMap();
72

73 74
		// get the grid padding
		const Padding<scheme::Sys_eqs_typ::dims> & pd = sc.getPadding();
75

76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111
		// 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);

		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
112 113 114
	 *
	 *
	 */
115
	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)
116
	{
117 118 119 120 121 122 123 124
		// 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)];
125

Pietro Incardona's avatar
Pietro Incardona committed
126 127 128
		// interpolation points for each property
		openfpm::vector<std::vector<comb<scheme::Sys_eqs_typ::dims>>> interp_pos[sizeof...(pos)];

129 130 131 132 133 134
		// 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);

Pietro Incardona's avatar
Pietro Incardona committed
135 136 137
		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);

138 139 140 141 142 143 144 145 146 147 148 149 150 151
		// 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);

		while (g_map_it.isNext() == true)
152
		{
153 154
			typedef typename to_boost_vmpl<pos...>::type vid;
			typedef boost::mpl::size<vid> v_size;
155

156
			auto key_src = g_map_it.get();
157

158 159
			// destination point
			auto key_dst = g_dst_it.get();
160

161 162
			// Transform this id into an id for the Eigen vector

Pietro Incardona's avatar
Pietro Incardona committed
163
			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);
164 165 166 167 168 169 170

			boost::mpl::for_each_ref<boost::mpl::range_c<int,0,v_size::value>>(cp);

			++g_map_it;
			++g_dst_it;
		}
	}
171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219

public:

	/*! \brief Resize the Vector
	 *
	 * \param row numbers of row
	 *
	 */
	void resize(size_t row)
	{
		v.resize(row);
	}

	/*! \brief Return a reference to the vector element
	 *
	 * \param i element
	 *
	 * \return reference to the element vector
	 *
	 */
	T & get(size_t i)
	{
		return v(i);
	}

	/*! \brief Get the Eigen Vector object
	 *
	 * \return the Eigen Vector
	 *
	 */
	const Eigen::Matrix<T, Eigen::Dynamic, 1> & getVec() const
	{
		return v;
	}

	/*! \brief Get the Eigen Vector object
	 *
	 * \return the Eigen Vector
	 *
	 */
	Eigen::Matrix<T, Eigen::Dynamic, 1> & getVec()
	{
		return v;
	}

	/*! \brief Copy the vector into the grid
	 *
	 *
	 */
220
	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)
221
	{
222
		if (is_grid_staggered<typename scheme::Sys_eqs_typ>::value())
223
		{
224 225 226 227 228 229
			if (g_dst.is_staggered() == true)
				copy_staggered_to_staggered<scheme,Grid_dst,pos...>(sc,start,stop,g_dst);
			else
				copy_staggered_to_normal<scheme,Grid_dst,pos...>(sc,start,stop,g_dst);
		}
	}
230

231 232 233 234 235 236 237 238 239
	/*! \brief Load from file
	 *
	 *
	 *
	 */
	void fromFile(std::string file)
	{
		std::ifstream inputf;
		inputf.open(file);
240

241 242 243 244
		for (size_t i = 0 ; i < v.size() ; i++)
			inputf >> v(i);

		inputf.close();
245 246 247 248 249 250

	}
};


#endif /* OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_EIGEN_HPP_ */