1/*
2 * Kernels.hpp
3 *
4 * Created on: Jan 12, 2016
5 * Author: i-bird
6 */
7
8#ifndef OPENFPM_NUMERICS_SRC_PSE_KERNELS_HPP_
9#define OPENFPM_NUMERICS_SRC_PSE_KERNELS_HPP_
10
11#include <boost/math/constants/constants.hpp>
12
13// Gaussian kernel
14#define KER_GAUSSIAN 1
15
16/*! \brief Implementation of the Laplacian kernels for PSE
17 *
18 * \tparam dim Dimension
19 * \tparam T type
20 * \ord order pf approximation (default 2)
21 * \impl TYPE of kernel
22 *
23 */
24template<unsigned int dim, typename T, unsigned int ord=2, unsigned int impl=KER_GAUSSIAN>
25struct Lap_PSE
26{
27 T epsilon;
28
29 Lap_PSE(T epsilon)
30 :epsilon(epsilon)
31 {}
32
33 /*! \brief From a kernel centered in x, it give the value of the kernel in y
34 *
35 * \param x center of the kernel
36 * \param y where we calculate the kernel
37 *
38 */
39 inline T value(T (&x)[dim], T (&y)[dim])
40 {
41 std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " The laplacian for order:" << ord << " in dimension " << dim << " has not been implemented";
42 return 0.0;
43 }
44
45 /*! \brief From a kernel centered in x, it give the value of the kernel in y
46 *
47 * \param x center of the kernel
48 * \param y where we calculate the kernel
49 *
50 */
51 inline T value(T (&x)[dim], const Point<dim,T> & y)
52 {
53 std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " The laplacian for order:" << ord << " in dimension " << dim << " has not been implemented";
54 return 0.0;
55 }
56
57 /*! \brief From a kernel centered in x, it give the value of the kernel in y
58 *
59 * \param x center of the kernel
60 * \param y where we calculate the kernel
61 *
62 */
63 inline T value(const Point<dim,T> & x, T (&y)[dim])
64 {
65 std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " The laplacian for order:" << ord << " in dimension " << dim << " has not been implemented";
66 return 0.0;
67 }
68
69 /*! \brief From a kernel centered in x, it give the value of the kernel in y
70 *
71 * \param x center of the kernel
72 * \param y where we calculate the kernel
73 *
74 */
75 inline T value(const Point<dim,T> & x, const Point<dim,T> & y)
76 {
77 std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " The laplacian for order:" << ord << " in dimension " << dim << " has not been implemented";
78 return 0.0;
79 }
80};
81
82template<typename T>
83struct Lap_PSE<1,T,2,KER_GAUSSIAN>
84{
85 T epsilon;
86
87 inline Lap_PSE(T epsilon)
88 :epsilon(epsilon)
89 {}
90
91 /*! \brief From a kernel centered in x, it give the value of the kernel in y
92 *
93 * \param x center of the kernel
94 * \param y where we calculate the kernel
95 *
96 */
97 inline T value(T (&x)[1], T (&y)[1])
98 {
99 T d = 0.0;
100 for (size_t i = 0 ; i < 1 ; i++)
101 d += (x[i] - y[i]) * (x[i] - y[i]);
102 d = sqrt(d) / epsilon;
103
104 return T(4.0) / epsilon / boost::math::constants::root_pi<T>() * exp(-d*d);
105 }
106
107 /*! \brief From a kernel centered in x, it give the value of the kernel in y
108 *
109 * \param x center of the kernel
110 * \param y where we calculate the kernel
111 *
112 */
113 inline T value(T (&x)[1], const Point<1,T> & y)
114 {
115 T d = 0.0;
116 for (size_t i = 0 ; i < 1 ; i++)
117 d += (x[i] - y.get(i)) * (x[i] - y.get(i));
118 d = sqrt(d) / epsilon;
119
120 return T(4.0) / epsilon / boost::math::constants::root_pi<T>() * exp(-d*d);
121 }
122
123 /*! \brief From a kernel centered in x, it give the value of the kernel in y
124 *
125 * \param x center of the kernel
126 * \param y where we calculate the kernel
127 *
128 */
129 inline T value(const Point<1,T> & x, T (&y)[1])
130 {
131 T d = 0.0;
132 for (size_t i = 0 ; i < 1 ; i++)
133 d += (x.get(i) - y[i]) * (x.get(i) - y[i]);
134 d = sqrt(d) / epsilon;
135
136 return T(4.0) / epsilon / boost::math::constants::root_pi<T>() * exp(-d*d);
137 }
138
139 /*! \brief From a kernel centered in x, it give the value of the kernel in y
140 *
141 * \param x center of the kernel
142 * \param y where we calculate the kernel
143 *
144 */
145 inline T value(const Point<1,T> & x, const Point<1,T> & y)
146 {
147 T d = 0.0;
148 for (size_t i = 0 ; i < 1 ; i++)
149 d += (x.get(i) - y.get(i)) * (x.get(i) - y.get(i));
150 d = sqrt(d) / epsilon;
151
152 return T(4.0) / epsilon / boost::math::constants::root_pi<T>() * exp(-d*d);
153 }
154};
155
156template<typename T>
157struct Lap_PSE<1,T,4,KER_GAUSSIAN>
158{
159 T epsilon;
160
161 inline Lap_PSE(T epsilon)
162 :epsilon(epsilon)
163 {}
164
165 /*! \brief From a kernel centered in x, it give the value of the kernel in y
166 *
167 * \param x center of the kernel
168 * \param y where we calculate the kernel
169 *
170 */
171 inline T value(T (&x)[1], T (&y)[1])
172 {
173 T d = 0.0;
174 for (size_t i = 0 ; i < 1 ; i++)
175 d += (x[i] - y[i]) * (x[i] - y[i]);
176 d = sqrt(d) / epsilon;
177
178 return T(1.0) / epsilon / boost::math::constants::root_pi<T>() * exp(-d*d) * (-4.0*d*d+10.0);
179 }
180
181 /*! \brief From a kernel centered in x, it give the value of the kernel in y
182 *
183 * \param x center of the kernel
184 * \param y where we calculate the kernel
185 *
186 */
187 inline T value(T (&x)[1], const Point<1,T> & y)
188 {
189 T d = 0.0;
190 for (size_t i = 0 ; i < 1 ; i++)
191 d += (x[i] - y.get(i)) * (x[i] - y.get(i));
192 d = sqrt(d) / epsilon;
193
194 return T(1.0) / epsilon / boost::math::constants::root_pi<T>() * exp(-d*d) * (-4.0*d*d+10.0);
195 }
196
197 /*! \brief From a kernel centered in x, it give the value of the kernel in y
198 *
199 * \param x center of the kernel
200 * \param y where we calculate the kernel
201 *
202 */
203 inline T value(const Point<1,T> & x, T (&y)[1])
204 {
205 T d = 0.0;
206 for (size_t i = 0 ; i < 1 ; i++)
207 d += (x.get(i) - y[i]) * (x.get(i) - y[i]);
208 d = sqrt(d) / epsilon;
209
210 return T(1.0) / epsilon / boost::math::constants::root_pi<T>() * exp(-d*d) * (-4.0*d*d+10.0);
211 }
212
213 /*! \brief From a kernel centered in x, it give the value of the kernel in y
214 *
215 * \param x center of the kernel
216 * \param y where we calculate the kernel
217 *
218 */
219 inline T value(const Point<1,T> & x, const Point<1,T> & y)
220 {
221 T d = 0.0;
222 for (size_t i = 0 ; i < 1 ; i++)
223 d += (x.get(i) - y.get(i)) * (x.get(i) - y.get(i));
224 d = sqrt(d) / epsilon;
225
226 return T(1.0) / epsilon / boost::math::constants::root_pi<T>() * exp(-d*d) * (-4.0*d*d+10.0);
227 }
228};
229
230template<typename T>
231struct Lap_PSE<1,T,6,KER_GAUSSIAN>
232{
233 T epsilon;
234
235 inline Lap_PSE(T epsilon)
236 :epsilon(epsilon)
237 {}
238
239 /*! \brief From a kernel centered in x, it give the value of the kernel in y
240 *
241 * \param x center of the kernel
242 * \param y where we calculate the kernel
243 *
244 */
245 inline T value(T (&x)[1], T (&y)[1])
246 {
247 T d = 0.0;
248 for (size_t i = 0 ; i < 1 ; i++)
249 d += (x[i] - y[i]) * (x[i] - y[i]);
250 d = sqrt(d) / epsilon;
251
252 return T(1.0) / epsilon / boost::math::constants::root_pi<T>() * exp(-d*d) * (2.0*d*d*d*d-14.0*d*d+35.0/2.0);
253 }
254
255 /*! \brief From a kernel centered in x, it give the value of the kernel in y
256 *
257 * \param x center of the kernel
258 * \param y where we calculate the kernel
259 *
260 */
261 inline T value(T (&x)[1], const Point<1,T> & y)
262 {
263 T d = 0.0;
264 for (size_t i = 0 ; i < 1 ; i++)
265 d += (x[i] - y.get(i)) * (x[i] - y.get(i));
266 d = sqrt(d) / epsilon;
267
268 return T(1.0) / epsilon / boost::math::constants::root_pi<T>() * exp(-d*d) * (2.0*d*d*d*d-14.0*d*d+35.0/2.0);
269 }
270
271 /*! \brief From a kernel centered in x, it give the value of the kernel in y
272 *
273 * \param x center of the kernel
274 * \param y where we calculate the kernel
275 *
276 */
277 inline T value(const Point<1,T> & x, T (&y)[1])
278 {
279 T d = 0.0;
280 for (size_t i = 0 ; i < 1 ; i++)
281 d += (x.get(i) - y[i]) * (x.get(i) - y[i]);
282 d = sqrt(d) / epsilon;
283
284 return T(1.0) / epsilon / boost::math::constants::root_pi<T>() * exp(-d*d) * (2.0*d*d*d*d-14.0*d*d+35.0/2.0);
285 }
286
287 /*! \brief From a kernel centered in x, it give the value of the kernel in y
288 *
289 * \param x center of the kernel
290 * \param y where we calculate the kernel
291 *
292 */
293 inline T value(const Point<1,T> & x, const Point<1,T> & y)
294 {
295 T d = 0.0;
296 for (size_t i = 0 ; i < 1 ; i++)
297 d += (x.get(i) - y.get(i)) * (x.get(i) - y.get(i));
298 d = sqrt(d) / epsilon;
299
300 return T(1.0) / epsilon / boost::math::constants::root_pi<T>() * exp(-d*d) * (2.0*d*d*d*d-14.0*d*d+35.0/2.0);
301 }
302};
303
304template<typename T>
305struct Lap_PSE<1,T,8,KER_GAUSSIAN>
306{
307 T epsilon;
308
309 inline Lap_PSE(T epsilon)
310 :epsilon(epsilon)
311 {}
312
313 /*! \brief From a kernel centered in x, it give the value of the kernel in y
314 *
315 * \param x center of the kernel
316 * \param y where we calculate the kernel
317 *
318 */
319 inline T value(T (&x)[1], T (&y)[1])
320 {
321 T d = 0.0;
322 for (size_t i = 0 ; i < 1 ; i++)
323 d += (x[i] - y[i]) * (x[i] - y[i]);
324 d = sqrt(d) / epsilon;
325
326 return T(1.0) / epsilon / boost::math::constants::root_pi<T>() * exp(-d*d) * (-T(2.0)/T(3.0)*d*d*d*d*d*d+9.0*d*d*d*d-63.0/2.0*d*d+105.0/4.0);
327 }
328
329 /*! \brief From a kernel centered in x, it give the value of the kernel in y
330 *
331 * \param x center of the kernel
332 * \param y where we calculate the kernel
333 *
334 */
335 inline T value(T (&x)[1], const Point<1,T> & y)
336 {
337 T d = 0.0;
338 for (size_t i = 0 ; i < 1 ; i++)
339 d += (x[i] - y.get(i)) * (x[i] - y.get(i));
340 d = sqrt(d) / epsilon;
341
342 return T(1.0) / epsilon / boost::math::constants::root_pi<T>() * exp(-d*d) * (-T(2.0)/T(3.0)*d*d*d*d*d*d+9.0*d*d*d*d-63.0/2.0*d*d+105.0/4.0);
343 }
344
345 /*! \brief From a kernel centered in x, it give the value of the kernel in y
346 *
347 * \param x center of the kernel
348 * \param y where we calculate the kernel
349 *
350 */
351 inline T value(const Point<1,T> & x, T (&y)[1])
352 {
353 T d = 0.0;
354 for (size_t i = 0 ; i < 1 ; i++)
355 d += (x.get(i) - y[i]) * (x.get(i) - y[i]);
356 d = sqrt(d) / epsilon;
357
358 return T(1.0) / epsilon / boost::math::constants::root_pi<T>() * exp(-d*d) * (-T(2.0)/T(3.0)*d*d*d*d*d*d+9.0*d*d*d*d-63.0/2.0*d*d+105.0/4.0);
359 }
360
361 /*! \brief From a kernel centered in x, it give the value of the kernel in y
362 *
363 * \param x center of the kernel
364 * \param y where we calculate the kernel
365 *
366 */
367 inline T value(const Point<1,T> & x, const Point<1,T> & y)
368 {
369 T d = 0.0;
370 for (size_t i = 0 ; i < 1 ; i++)
371 d += (x.get(i) - y.get(i)) * (x.get(i) - y.get(i));
372 d = sqrt(d) / epsilon;
373
374 return T(1.0) / epsilon / boost::math::constants::root_pi<T>() * exp(-d*d) * (-T(2.0)/T(3.0)*d*d*d*d*d*d+9.0*d*d*d*d-63.0/2.0*d*d+105.0/4.0);
375 }
376};
377
378#endif /* OPENFPM_NUMERICS_SRC_PSE_KERNELS_HPP_ */
379