Laplacian.hpp 6.69 KB
Newer Older
incardon's avatar
incardon committed
1 2 3 4 5 6 7 8 9 10
/*
 * Laplacian.hpp
 *
 *  Created on: Oct 5, 2015
 *      Author: i-bird
 */

#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_LAPLACIAN_HPP_
#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_LAPLACIAN_HPP_

11 12 13
#include "FD_util_include.hpp"
#include "util/util_num.hpp"
#include "FiniteDifference/eq.hpp"
incardon's avatar
incardon committed
14 15 16 17 18 19 20 21 22 23

/*! \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
{
24
	/*! \brief Calculate which colums of the Matrix are non zero
incardon's avatar
incardon committed
25
	 *
26 27 28 29 30 31 32 33 34 35 36
	 * stub_or_real it is just for change the argument type when testing, in normal
	 * conditions it is a distributed map
	 *
	 * \param pos position where the laplacian is calculated
	 * \param gs Grid info
	 * \param cols non-zero colums calculated by the function
	 * \param coeff coefficent (constant in front of the derivative)
	 *
	 * ### Example
	 *
	 * \snippet FDScheme_unit_tests.hpp Laplacian usage
incardon's avatar
incardon committed
37 38
	 *
	 */
39
	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)
incardon's avatar
incardon committed
40 41 42 43
	{
		std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
	}

44
	/*! \brief Calculate the position where the derivative is calculated
incardon's avatar
incardon committed
45
	 *
46 47 48
	 * \param position where we are calculating the derivative
	 * \param gs Grid info
	 * \param s_pos staggered position of the properties
incardon's avatar
incardon committed
49 50 51 52 53 54 55 56
	 *
	 */
	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";
	}
};

57 58 59 60 61 62 63 64 65 66 67 68 69
/*! \brief Laplacian second order approximation CENTRAL Scheme
 *
 * \verbatim

          1.0
           *
           | -4.0
   1.0 *---+---* 1.0
           |
           *
          1.0

 * \endverbatim
incardon's avatar
incardon committed
70 71 72 73 74 75 76 77
 *
 *
 */
template<typename arg, typename Sys_eqs>
class Lap<arg,Sys_eqs,CENTRAL>
{
public:

78 79 80 81 82

	/*! \brief Calculate which colums of the Matrix are non zero
	 *
	 * stub_or_real it is just for change the argument type when testing, in normal
	 * conditions it is a distributed map
incardon's avatar
incardon committed
83
	 *
84 85 86 87 88 89 90 91
	 * \param pos position where the laplacian is calculated
	 * \param gs Grid info
	 * \param cols non-zero colums calculated by the function
	 * \param coeff coefficent (constant in front of the derivative)
	 *
	 * ### Example
	 *
	 * \snippet FDScheme_unit_tests.hpp Laplacian usage
incardon's avatar
incardon committed
92 93
	 *
	 */
94
	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)
incardon's avatar
incardon committed
95 96 97 98
	{
		// for each dimension
		for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
		{
99 100
			long int old_val = kmap.getKeyRef().get(i);
			kmap.getKeyRef().set_d(i, kmap.getKeyRef().get(i) + 1);
101
			arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[i]/spacing[i]);
102
			kmap.getKeyRef().set_d(i,old_val);
incardon's avatar
incardon committed
103

104 105
			old_val = kmap.getKeyRef().get(i);
			kmap.getKeyRef().set_d(i, kmap.getKeyRef().get(i) - 1);
106
			arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[i]/spacing[i]);
107
			kmap.getKeyRef().set_d(i,old_val);
incardon's avatar
incardon committed
108

109 110
			arg::value(g_map,kmap,gs,spacing,cols, - 2.0 * coeff/spacing[i]/spacing[i]);
		}
incardon's avatar
incardon committed
111 112 113 114 115
	}


	/*! \brief Calculate the position where the derivative is calculated
	 *
116 117 118 119 120 121
	 * In case of non staggered case this function just return a null grid_key, in case of staggered,
	 * the CENTRAL Laplacian scheme return the position of the staggered property
	 *
	 * \param position where we are calculating the derivative
	 * \param gs Grid info
	 * \param s_pos staggered position of the properties
incardon's avatar
incardon committed
122 123
	 *
	 */
124
	inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, const comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
incardon's avatar
incardon committed
125
	{
126
		return arg::position(pos,gs,s_pos);
incardon's avatar
incardon committed
127 128 129 130
	}
};


incardon's avatar
incardon committed
131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 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
/*! \brief Laplacian second order approximation CENTRAL Scheme (with central derivative in the single)
 *
 * \verbatim

          1.0
           *
           | -4.0
   1.0 *---+---* 1.0
           |
           *
          1.0

 * \endverbatim
 *
 *
 */
template<typename arg, typename Sys_eqs>
class Lap<arg,Sys_eqs,CENTRAL_SYM>
{
public:


	/*! \brief Calculate which colums of the Matrix are non zero
	 *
	 * stub_or_real it is just for change the argument type when testing, in normal
	 * conditions it is a distributed map
	 *
	 * \param pos position where the laplacian is calculated
	 * \param gs Grid info
	 * \param cols non-zero colums calculated by the function
	 * \param coeff coefficent (constant in front of the derivative)
	 *
	 * ### Example
	 *
	 * \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::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++)
		{
			long int old_val = kmap.getKeyRef().get(i);
			kmap.getKeyRef().set_d(i, kmap.getKeyRef().get(i) + 2);
			arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[i]/spacing[i]/4.0);
			kmap.getKeyRef().set_d(i,old_val);

			old_val = kmap.getKeyRef().get(i);
			kmap.getKeyRef().set_d(i, kmap.getKeyRef().get(i) - 2);
			arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[i]/spacing[i]/4.0);
			kmap.getKeyRef().set_d(i,old_val);

			arg::value(g_map,kmap,gs,spacing,cols, - 2.0 * coeff/spacing[i]/spacing[i]/4.0);
		}
	}


	/*! \brief Calculate the position where the derivative is calculated
	 *
	 * In case of non staggered case this function just return a null grid_key, in case of staggered,
	 * the CENTRAL Laplacian scheme return the position of the staggered property
	 *
	 * \param position where we are calculating the derivative
	 * \param gs Grid info
	 * \param s_pos staggered position of the properties
	 *
	 */
	inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, const comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
	{
		return arg::position(pos,gs,s_pos);
	}
};

incardon's avatar
incardon committed
204
#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_LAPLACIAN_HPP_ */