Adding missing files

src/FiniteDifference/Derivative.hpp

+/*
+ * Derivative.hpp
+ *
+ *  Created on: Oct 5, 2015
+ *      Author: i-bird
+ */
+
+#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_DERIVATIVE_HPP_
+#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_DERIVATIVE_HPP_
+
+#define CENTRAL 0
+#define CENTRAL_B_ONE_SIDE 1
+
+#include "util/mathutil.hpp"
+#include "Vector/map_vector.hpp"
+#include "Grid/comb.hpp"
+#include "FiniteDifference/util/common.hpp"
+
+/*! \brief Derivative second order on h (spacing)
+ *
+ * \tparam d on which dimension derive
+ * \tparam Field which field derive
+ * \tparam impl which implementation
+ *
+ */
+template<unsigned int d, typename Field, typename Sys_eqs, unsigned int impl=CENTRAL>
+class D
+{
+	/*! \brief Create the row of the Matrix
+	 *
+	 * \tparam ord
+	 *
+	 */
+	inline static std::unordered_map<long int,typename Sys_eqs::stype> value(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs)
+	{
+		std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
+	}
+
+	/*! \brief Calculate the position where the derivative is calculated
+	 *
+	 * In case on non staggered case this function just return pos, in case of staggered,
+	 *  it calculate where the operator is calculated on a staggered grid
+	 *
+	 */
+	inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs)
+	{
+		std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
+	}
+};
+
+/*! \brief Derivative on direction i
+ *
+ *
+ */
+template<unsigned int d, typename arg, typename Sys_eqs>
+class D<d,arg,Sys_eqs,CENTRAL>
+{
+	public:
+
+	/*! \brief fill the row
+	 *
+	 *
+	 */
+	inline static void value(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
+	{
+		// forward
+		if (Sys_eqs::boundary[d] == PERIODIC )
+		{
+			long int old_val = pos.get(d);
+			pos.set_d(d, openfpm::math::positive_modulo(pos.get(d) + 1, gs.size(d)));
+			arg::value(pos,gs,cols,coeff);
+			pos.set_d(d,old_val);
+		}
+		else
+		{
+			long int old_val = pos.get(d);
+			pos.set_d(d, pos.get(d) + 1);
+			arg::value(pos,gs,cols,coeff);
+			pos.set_d(d,old_val);
+		}
+
+		// backward
+		if (Sys_eqs::boundary[d] == PERIODIC )
+		{
+			long int old_val = pos.get(d);
+			pos.set_d(d, openfpm::math::positive_modulo(pos.get(d) - 1, gs.size(d)));
+			arg::value(pos,gs,cols,-coeff);
+			pos.set_d(d,old_val);
+		}
+		else
+		{
+			long int old_val = pos.get(d);
+			pos.set_d(d, pos.get(d) - 1);
+			arg::value(pos,gs,cols,-coeff);
+			pos.set_d(d,old_val);
+		}
+	}
+
+
+	/*! \brief Calculate the position where the derivative is calculated
+	 *
+	 * In case on non staggered case this function just return pos, in case of staggered,
+	 *  it calculate where the operator is calculated on a staggered grid
+	 *
+	 *  \param pos from the position
+	 *  \param fld Field we are deriving, if not provided the function just return pos
+	 *  \param s_pos position of the properties in the staggered grid
+	 *
+	 */
+	inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, long int fld = -1, const openfpm::vector<comb<Sys_eqs::dims>> & s_pos = openfpm::vector<comb<Sys_eqs::dims>>())
+	{
+		if (is_grid_staggered<Sys_eqs>::value())
+		{
+			if (fld == -1)
+				return pos;
+
+			if (s_pos[fld][d] == 1)
+			{
+				grid_key_dx<Sys_eqs::dims> ret = pos;
+				ret.set_d(d,0);
+				return pos;
+			}
+			else
+			{
+				grid_key_dx<Sys_eqs::dims> ret = pos;
+				ret.set_d(d,1);
+				return pos;
+			}
+		}
+
+		return pos;
+	}
+};
+
+
+/*! \brief Derivative on direction i
+ *
+ *
+ */
+template<unsigned int d, typename arg, typename Sys_eqs>
+class D<d,arg,Sys_eqs,CENTRAL_B_ONE_SIDE>
+{
+public:
+
+	/*! \brief fill the row
+	 *
+	 *
+	 */
+	static void value(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, std::unordered_map<long int,typename Sys_eqs::stype> & cols, typename Sys_eqs::stype coeff)
+	{
+#ifdef SE_CLASS1
+		if (Sys_eqs::boundary[d] == PERIODIC)
+			std::cerr << __FILE__ << ":" << __LINE__ << " error, it make no sense use one sided derivation with periodic boundary\n";
+#endif
+
+		if (pos.get(d) == (long int)gs.size(d)-1 )
+		{
+			arg::value(pos,gs,cols,1.5*coeff);
+
+			long int old_val = pos.get(d);
+			pos.set_d(d, pos.get(d) - 1);
+			arg::value(pos,gs,cols,-2.0*coeff);
+			pos.set_d(d,old_val);
+
+			old_val = pos.get(d);
+			pos.set_d(d, pos.get(d) - 2);
+			arg::value(pos,gs,cols,0.5*coeff);
+			pos.set_d(d,old_val);
+		}
+		else if (pos.get(d) == 0)
+		{
+			arg::value(pos,gs,cols,-1.5*coeff);
+
+			long int old_val = pos.get(d);
+			pos.set_d(d, pos.get(d) + 1);
+			arg::value(pos,gs,cols,2.0*coeff);
+			pos.set_d(d,old_val);
+
+			old_val = pos.get(d);
+			pos.set_d(d, pos.get(d) + 2);
+			arg::value(pos,gs,cols,-0.5*coeff);
+			pos.set_d(d,old_val);
+		}
+		else
+		{
+			long int old_val = pos.get(d);
+			pos.set_d(d, pos.get(d) + 1);
+			arg::value(pos,gs,cols,coeff);
+			pos.set_d(d,old_val);
+
+			old_val = pos.get(d);
+			pos.set_d(d, pos.get(d) - 1);
+			arg::value(pos,gs,cols,-coeff);
+			pos.set_d(d,old_val);
+		}
+	}
+
+	/*! \brief Calculate the position where the derivative is calculated
+	 *
+	 * In case on non staggered case this function just return pos, in case of staggered,
+	 *  it calculate where the operator is calculated on a staggered grid
+	 *
+	 */
+	inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, long int fld = 0, const openfpm::vector<comb<Sys_eqs::dims>> & s_pos = openfpm::vector<comb<Sys_eqs::dims>>())
+	{
+		if (is_grid_staggered<Sys_eqs>::type::value)
+		{
+			if (fld == -1)
+				return pos;
+
+			if (s_pos[fld][d] == 1)
+			{
+				grid_key_dx<Sys_eqs::dims> ret = pos;
+				ret.set_d(d,0);
+				return pos;
+			}
+			else
+			{
+				grid_key_dx<Sys_eqs::dims> ret = pos;
+				ret.set_d(d,1);
+				return pos;
+			}
+		}
+
+		return pos;
+	}
+};
+
+
+
+#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_DERIVATIVE_HPP_ */

src/FiniteDifference/FDScheme.hpp

+/*
+ * FiniteDifferences.hpp
+ *
+ *  Created on: Sep 17, 2015
+ *      Author: i-bird
+ */
+
+#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_HPP_
+#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_HPP_
+
+#include "Grid/grid_dist_id.hpp"
+#include "Matrix/Matrix.hpp"
+#include "eq.hpp"
+
+/*! \brief Finite Differences
+ *
+ * \tparam dim Dimensionality of the finite differences scheme
+ *
+ */
+
+template<typename Sys_eqs>
+class FDScheme
+{
+	openfpm::vector<triplet<typename Sys_eqs::stype>> trpl;
+
+	/*! \brief Impose an operator
+	 *
+	 *
+	 *
+	 */
+/*	template<typename T> void impose(T & op, grid_key_dx_dist_iterator_sub & it)
+	{
+		size_t lin = 0;
+
+		std::unordered_map<long int,float> cols;
+
+		// iterate all the grid points
+		while (it.isNext())
+		{
+			// get the position
+			auto key = it.get();
+
+			// convert into global coordinate the position
+			auto keyg = g.getGKey(key);
+
+			// Convert local key into global key
+			T::value(op,cols);
+
+			// create the triplet
+
+			for ( auto it = cols.begin(); it != cols.end(); ++it )
+			{
+				trpl.add();
+				trpl.last().i = lin;
+				trpl.last().j = it->first;
+				trpl.last().value = it->second;
+			}
+
+			std::cout << " " << it->first << ":" << it->second;
+
+			cols.clear();
+
+			++loc_cnt;
+			++it;
+		}
+	}*/
+
+	/*! \brief produce the Matrix
+	 *
+	 *  \tparam Syst_eq System of equations, or a single equation
+	 *
+	 */
+	template<typename Grid> SparseMatrix<typename Sys_eqs::stype> getMatrix(const Grid & g, const ConstField(& fld)[Sys_eqs::num_cfields::value] )
+	{
+		// iterate across the full domain
+
+		auto it = g.getDomainIterator();
+		auto g_sm = g.getGrid();
+
+		Sys_eqs eqs();
+
+		// iterate all the grid points
+		while (it.isNext())
+		{
+			// get the position
+			auto key = it.get();
+
+			// convert into global coordinate the position
+			auto keyg = g.getGKey(key);
+
+			// Convert local key into global key
+			size_t row = g_sm.LidId(keyg);
+
+			eqs.value(keyg,trpl);
+
+			++it;
+		}
+	}
+};
+
+#define EQS_FIELDS 0
+#define EQS_SPACE 1
+
+
+#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_HPP_ */

src/FiniteDifference/Laplacian.hpp

+/*
+ * Laplacian.hpp
+ *
+ *  Created on: Oct 5, 2015
+ *      Author: i-bird
+ */
+
+#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_LAPLACIAN_HPP_
+#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_LAPLACIAN_HPP_
+
+
+/*! \brief Laplacian second order on h (spacing)
+ *
+ * \tparam Field which field derive
+ * \tparam impl which implementation
+ *
+ */
+template<typename arg, typename Sys_eqs, unsigned int impl=CENTRAL>
+class Lap
+{
+	/*! \brief Create the row of the Matrix
+	 *
+	 * \tparam ord
+	 *
+	 */
+	inline static std::unordered_map<long int,typename Sys_eqs::stype> value(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs)
+	{
+		std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
+	}
+
+	/*! \brief Calculate the position where the laplacian is calculated
+	 *
+	 * In case on non staggered case this function just return pos, in case of staggered,
+	 *  it calculate where the operator is calculated on a staggered grid
+	 *
+	 */
+	inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs)
+	{
+		std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
+	}
+};
+
+/*! \brief Derivative on direction i
+ *
+ *
+ */
+template<typename arg, typename Sys_eqs>
+class Lap<arg,Sys_eqs,CENTRAL>
+{
+public:
+
+	/*! \brief fill the row
+	 *
+	 *
+	 */
+	inline static void value(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, 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++)
+		{
+			// forward
+			if (Sys_eqs::boundary[i] == PERIODIC )
+			{
+				long int old_val = pos.get(i);
+				pos.set_d(i, openfpm::math::positive_modulo(pos.get(i) + 1, gs.size(i)));
+				arg::value(pos,gs,cols,coeff);
+				pos.set_d(i,old_val);
+			}
+			else
+			{
+				long int old_val = pos.get(i);
+				pos.set_d(i, pos.get(i) + 1);
+				arg::value(pos,gs,cols,coeff);
+				pos.set_d(i,old_val);
+			}
+
+			// backward
+			if (Sys_eqs::boundary[i] == PERIODIC )
+			{
+				long int old_val = pos.get(i);
+				pos.set_d(i, openfpm::math::positive_modulo(pos.get(i) - 1, gs.size(i)));
+				arg::value(pos,gs,cols,coeff);
+				pos.set_d(i,old_val);
+			}
+			else
+			{
+				long int old_val = pos.get(i);
+				pos.set_d(i, pos.get(i) - 1);
+				arg::value(pos,gs,cols,coeff);
+				pos.set_d(i,old_val);
+			}
+		}
+
+		arg::value(pos,gs,cols, - 2.0 * Sys_eqs::dims * coeff);
+	}
+
+
+	/*! \brief Calculate the position where the derivative is calculated
+	 *
+	 * In case of non staggered case this function just return pos, in case of staggered,
+	 * it calculate where the operator is calculated on a staggered grid
+	 *
+	 */
+	inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, const openfpm::vector<comb<Sys_eqs::dims>> & s_pos, long int & fld)
+	{
+		return pos;
+	}
+};
+
+
+#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_LAPLACIAN_HPP_ */

src/FiniteDifference/System.hpp

+/*
+ * System.hpp
+ *
+ *  Created on: Oct 5, 2015
+ *      Author: i-bird
+ */
+
+#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_SYSTEM_HPP_
+#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_SYSTEM_HPP_
+
+
+/*! \brief System of equations
+ *
+ * This class model a system of equations
+ *
+ * \tparam dim Dimensionality of the system
+ * \tparam nvf number of variable fields
+ * \tparam ncf number of constants fields
+ * \tparam eqs sets of equations
+ *
+ */
+template<unsigned int dim, unsigned int nvf, unsigned int ncf, typename ... eqs>
+class System
+{
+	// Define the number of constant fields
+	typedef num_cfields boost::mpl::int_<nf>;
+
+	// Define the number of variable fields
+	typedef num_vfields boost::mpl::int_<nf>;
+
+	// set of equations as boost::mpl::vector
+	typedef eqs_v make_vactor<eqs>;
+
+	/*! \brief Create the row of the Matrix
+	 *
+	 * \tparam ord
+	 *
+	 */
+	template<unsigned int ord=EQS_FIELDS> void value(const grid_key_dx_dist<dim> & pos)
+	{
+		if (EQS_FIELDS)
+			value_f(pos);
+		else
+			value_s(pos);
+	}
+
+	/*! \brief fill the row
+	 *
+	 *
+	 */
+	template<unsigned int eq_id> void value_s(grid_key_dx_dist<dim> & it)
+	{
+		boost::mpl::at<eqs_v,boost::mpl::int_<eq_id>>::type eq;
+
+		eq.value_s(it);
+	}
+
+	/*! \brief fill the row
+	 *
+	 *
+	 */
+	template<unsigned int eq_id> void value_f(grid_key_dx_dist<dim> & it)
+	{
+		boost::mpl::at<eqs_v,boost::mpl::int_<eq_id>>::type eq;
+
+		eq.value_f(it);
+	}
+};
+
+
+#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_SYSTEM_HPP_ */

src/FiniteDifference/eq.hpp

+/*
+ * eq.hpp
+ *
+ *  Created on: Oct 5, 2015
+ *      Author: i-bird
+ */
+
+#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_EQ_HPP_
+#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_EQ_HPP_
+
+#define EQS_FIELD 0
+#define EQS_POS 1
+
+#define STAGGERED_GRID 1
+#define NORMAL_GRID 0
+
+#define PERIODIC true
+#define NON_PERIODIC false
+
+/*! \brief Equation
+ *
+ * It model an equation like expr1 = expr2
+ *
+ * \tparam expr1
+ * \tparam expr2
+ *
+ */
+template<typename expr1,typename expr2,typename Sys_eqs>
+class Eq
+{
+	/*! \brief Create the row of the Matrix
+	 *
+	 * \tparam ord
+	 *
+	 */
+	template<unsigned int ord=EQS_FIELD> static void value(const grid_key_dx<Sys_eqs::dims> & pos)
+	{
+		if (EQS_FIELD)
+			value_f(pos);
+		else
+			value_s(pos);
+	}
+
+	/*! \brief fill the row
+	 *
+	 *
+	 */
+	static openfpm::vector<cval<typename Sys_eqs::stype>> value_s(grid_key_dx<Sys_eqs::dims> & it)
+	{
+		return expr1::value_s(it) - expr2::value_s(it);
+	}
+
+	/*! \brief fill the row
+	 *
+	 *
+	 */
+	static void value_f(grid_key_dx<Sys_eqs::dims> & it)
+	{
+		return expr1::value_s(it) - expr2::value_s(it);
+	}
+};
+
+
+
+/*! \brief It model an expression expr1 - expr2
+ *
+ * \tparam expr1
+ * \tparam expr2
+ *
+ */
+template<typename expr1,typename expr2, typename Sys_eqs>
+class minus
+{
+	/*! \brief Create the row of the Matrix
+	 *
+	 * \tparam ord
+	 *
+	 */
+	template<unsigned int ord=EQS_FIELD> static void value(const grid_key_dx<Sys_eqs::dims> & pos)
+	{
+		if (EQS_FIELD)
+			value_f(pos);
+		else
+			value_s(pos);
+	}
+
+	/*! \brief fill the row
+	 *
+	 *
+	 */
+	static openfpm::vector<triplet<typename Sys_eqs::stype>> value_s(grid_key_dx<Sys_eqs::dims> & it)
+	{
+		return expr1::value_s(it) - expr2::value_s(it);
+	}
+
+	/*! \brief fill the row
+	 *
+	 *
+	 */
+	static void value_f(grid_key_dx<Sys_eqs::dims> & it)
+	{
+		return expr1::value_s(it) - expr2::value_s(it);
+	}
+};
+
+/*! \brief It model an expression expr1 * expr2
+ *
+ * \warning expr1 MUST be a constant expression
+ *
+ * \tparam expr1
+ * \tparam expr2
+ *
+ */
+template<typename expr1,typename expr2, typename Sys_eqs>
+class mul