1/*
2 * Average.hpp
3 *
4 * Created on: Nov 18, 2015
5 * Author: Pietro Incardona
6 */
7
8#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_AVERAGE_HPP_
9#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_AVERAGE_HPP_
10
11
12#define CENTRAL 0
13#define CENTRAL_B_ONE_SIDE 1
14#define FORWARD 2
15#define BACKWARD 3
16
17#include "util/mathutil.hpp"
18#include "Vector/map_vector.hpp"
19#include "Grid/comb.hpp"
20#include "FiniteDifference/util/common.hpp"
21#include "util/util_num.hpp"
22
23/*! \brief Average
24 *
25 * \tparam d on which dimension average
26 * \tparam Field which field average
27 * \tparam impl which implementation
28 *
29 */
30template<unsigned int d, typename Field, typename Sys_eqs, unsigned int impl=CENTRAL>
31class Avg
32{
33 /*! \brief Create the row of the Matrix
34 *
35 * \tparam ord
36 *
37 */
38 inline static void value(const 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)
39 {
40 std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
41 }
42
43 /*! \brief Calculate the position where the derivative is calculated
44 *
45 * In case of non staggered case this function just return a null grid_key, in case of staggered,
46 * it calculate how the operator shift the calculation in the cell
47 *
48 * \param pos position where we are calculating the derivative
49 * \param gs Grid info
50 * \param s_pos staggered position of the properties
51 *
52 */
53 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 openfpm::vector<comb<Sys_eqs::dims>> & s_pos = openfpm::vector<comb<Sys_eqs::dims>>())
54 {
55 std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
56 }
57};
58
59/*! \brief Central average scheme on direction i
60 *
61 * \verbatim
62 *
63 * +0.5 +0.5
64 * *---+---*
65 *
66 * \endverbatim
67 *
68 */
69template<unsigned int d, typename arg, typename Sys_eqs>
70class Avg<d,arg,Sys_eqs,CENTRAL>
71{
72 public:
73
74 /*! \brief Calculate which colums of the Matrix are non zero
75 *
76 * stub_or_real it is just for change the argument type when testing, in normal
77 * conditions it is a distributed map
78 *
79 * \param g_map It is the map explained in FDScheme
80 * \param kmap position where the average is calculated
81 * \param gs Grid info
82 * \param cols non-zero colums calculated by the function
83 * \param coeff coefficent (constant in front of the derivative)
84 *
85 * ### Example
86 *
87 * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
88 *
89 */
90 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)
91 {
92 // if the system is staggered the CENTRAL derivative is equivalent to a forward derivative
93 if (is_grid_staggered<Sys_eqs>::value())
94 {
95 Avg<d,arg,Sys_eqs,BACKWARD>::value(g_map,kmap,gs,spacing,cols,coeff);
96 return;
97 }
98
99 long int old_val = kmap.getKeyRef().get(d);
100 kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 1);
101 arg::value(g_map,kmap,gs,spacing,cols,coeff/2);
102 kmap.getKeyRef().set_d(d,old_val);
103
104
105 old_val = kmap.getKeyRef().get(d);
106 kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 1);
107 arg::value(g_map,kmap,gs,spacing,cols,coeff/2);
108 kmap.getKeyRef().set_d(d,old_val);
109 }
110
111
112 /*! \brief Calculate the position where the average is calculated
113 *
114 * It follow the same concept of central derivative
115 *
116 * \param pos position where we are calculating the derivative
117 * \param gs Grid info
118 * \param s_pos staggered position of the properties
119 *
120 */
121 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])
122 {
123 auto arg_pos = arg::position(pos,gs,s_pos);
124 if (is_grid_staggered<Sys_eqs>::value())
125 {
126 if (arg_pos.get(d) == -1)
127 {
128 arg_pos.set_d(d,0);
129 return arg_pos;
130 }
131 else
132 {
133 arg_pos.set_d(d,-1);
134 return arg_pos;
135 }
136 }
137
138 return arg_pos;
139 }
140};
141
142/*! \brief FORWARD average on direction i
143 *
144 * \verbatim
145 *
146 * +0.5 0.5
147 * +------*
148 *
149 * \endverbatim
150 *
151 */
152template<unsigned int d, typename arg, typename Sys_eqs>
153class Avg<d,arg,Sys_eqs,FORWARD>
154{
155 public:
156
157 /*! \brief Calculate which colums of the Matrix are non zero
158 *
159 * stub_or_real it is just for change the argument type when testing, in normal
160 * conditions it is a distributed map
161 *
162 * \param g_map It is the map explained in FDScheme
163 * \param kmap position where the average is calculated
164 * \param gs Grid info
165 * \param cols non-zero colums calculated by the function
166 * \param coeff coefficent (constant in front of the derivative)
167 *
168 * ### Example
169 *
170 * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
171 *
172 */
173 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)
174 {
175
176 long int old_val = kmap.getKeyRef().get(d);
177 kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 1);
178 arg::value(g_map,kmap,gs,spacing,cols,coeff/2);
179 kmap.getKeyRef().set_d(d,old_val);
180
181 // backward
182 arg::value(g_map,kmap,gs,spacing,cols,coeff/2);
183 }
184
185
186 /*! \brief Calculate the position in the cell where the average is calculated
187 *
188 * In case of non staggered case this function just return a null grid_key, in case of staggered,
189 * the FORWARD scheme return the position of the staggered property
190 *
191 * \param pos position where we are calculating the derivative
192 * \param gs Grid info
193 * \param s_pos staggered position of the properties
194 *
195 */
196 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])
197 {
198 return arg::position(pos,gs,s_pos);
199 }
200};
201
202/*! \brief First order BACKWARD derivative on direction i
203 *
204 * \verbatim
205 *
206 * +0.5 0.5
207 * *------+
208 *
209 * \endverbatim
210 *
211 */
212template<unsigned int d, typename arg, typename Sys_eqs>
213class Avg<d,arg,Sys_eqs,BACKWARD>
214{
215 public:
216
217 /*! \brief Calculate which colums of the Matrix are non zero
218 *
219 * stub_or_real it is just for change the argument type when testing, in normal
220 * conditions it is a distributed map
221 *
222 * \param g_map It is the map explained in FDScheme
223 * \param kmap position where the average is calculated
224 * \param gs Grid info
225 * \param cols non-zero colums calculated by the function
226 * \param coeff coefficent (constant in front of the derivative)
227 *
228 * ### Example
229 *
230 * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
231 *
232 */
233 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)
234 {
235 long int old_val = kmap.getKeyRef().get(d);
236 kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 1);
237 arg::value(g_map,kmap,gs,spacing,cols,coeff/2);
238 kmap.getKeyRef().set_d(d,old_val);
239
240 // forward
241 arg::value(g_map,kmap,gs,spacing,cols,coeff/2);
242 }
243
244
245 /*! \brief Calculate the position in the cell where the average is calculated
246 *
247 * In case of non staggered case this function just return a null grid_key, in case of staggered,
248 * the BACKWARD scheme return the position of the staggered property
249 *
250 * \param pos position where we are calculating the derivative
251 * \param gs Grid info
252 * \param s_pos staggered position of the properties
253 *
254 */
255 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])
256 {
257 return arg::position(pos,gs,s_pos);
258 }
259};
260
261#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_AVERAGE_HPP_ */
262