1/*
2 * mp4_interpolation.hpp
3 *
4 * Created on: May 4, 2017
5 * Author: i-bird
6 */
7
8#ifndef OPENFPM_NUMERICS_SRC_INTERPOLATION_INTERPOLATION_HPP_
9#define OPENFPM_NUMERICS_SRC_INTERPOLATION_INTERPOLATION_HPP_
10
11#include "NN/Mem_type/MemFast.hpp"
12#include "NN/CellList/CellList.hpp"
13#include "Grid/grid_dist_key.hpp"
14#include "Vector/vector_dist_key.hpp"
15
16#define INTERPOLATION_ERROR_OBJECT std::runtime_error("Runtime interpolation error");
17
18constexpr int inte_m2p = 0;
19constexpr int inte_p2m = 1;
20
21/*! \brief It store the offsets of the interpolation points
22 *
23 * \tparam n_ele number of interpolation points
24 * \tparam T type in general an
25 *
26 */
27template<unsigned int n_ele>
28struct agg_arr
29{
30 //! offset of the interpolation points
31 size_t ele[n_ele];
32};
33
34/*! \brief multiply the src by coeff for several types T
35 *
36 * \tparam type T
37 *
38 */
39template<typename T>
40struct mul_inte
41{
42 /*! \brief multiply the src by coeff for several types T
43 *
44 * \param result the result of the multiplication
45 * \param coeff coefficent to use for of the multiplication
46 * \param src source
47 *
48 */
49 inline static void value(T & result, const T & coeff, const T & src)
50 {
51 result += coeff * src;
52 }
53};
54
55/*! \brief multiply the src by coeff for several types T
56 *
57 * \tparam type T
58 *
59 */
60template<typename T, unsigned int N1>
61struct mul_inte<T[N1]>
62{
63 /*! \brief multiply the src by coeff for several types T
64 *
65 * \param result the result of the multiplication
66 * \param coeff coefficent to use for of the multiplication
67 * \param src source
68 *
69 */
70 inline static void value(T (& result)[N1], const T & coeff, const T (& src)[N1])
71 {
72 for (size_t i = 0 ; i < N1 ; i++)
73 result[i] += coeff * src[i];
74 }
75};
76
77/*! \brief multiply the src by coeff for several types T
78 *
79 * \tparam type T
80 *
81 */
82template<typename T, unsigned int N1, unsigned int N2>
83struct mul_inte<T[N1][N2]>
84{
85 /*! \brief multiply the src by coeff for several types T
86 *
87 * \param result the result of the multiplication
88 * \param coeff coefficent to use for of the multiplication
89 * \param src source
90 *
91 */
92 inline static void value(T (& result)[N1][N2], const T & coeff, const T (& src)[N1][N2])
93 {
94 for (size_t i = 0 ; i < N1 ; i++)
95 for (size_t j = 0 ; j < N2 ; j++)
96 result[i][j] += coeff * src[i][j];
97 }
98};
99
100/*! \brief Class that select the operation to do differently if we are doing Mesh to particle (m2p) or particle to mesh (p2m)
101 *
102 * \tparam prp_g property of the grid to interpolate
103 * \tparam prp_v property of the vector to interpolate
104 * \param M2P or P2M
105 *
106 */
107template<unsigned int np, unsigned int prp_g, unsigned int prp_v,unsigned int m2p_or_p2m>
108struct inte_template
109{
110 /*! \brief Evaluate the interpolation
111 *
112 * \tparam np number of kernel points in one direction
113 * \tparam grid type of grid
114 * \tparam vector type of vector
115 * \tparam iterator type of the iterator
116 *
117 * \param gd grid for interpolation
118 * \param vd vector for interpolation
119 * \param k_dist grid key grid point for interpolation
120 * \param key_p particle for interpolation
121 * \param a_int interpolation coefficients pre-calculated
122 * \param key indicate which pre-calculated coefficient we have to use
123 *
124 */
125 template<unsigned int np_a_int, typename grid, typename vector, typename iterator> inline static void value(grid & gd,
126 vector & vd,
127 const grid_dist_lin_dx & k_dist,
128 iterator & key_p,
129 typename vector::stype (& a_int)[np_a_int],
130 const size_t & key)
131 {
132 mul_inte<typename std::remove_reference<decltype(gd.template get<prp_g>(k_dist))>::type>::value(gd.template get<prp_g>(k_dist),a_int[key],vd.template getProp<prp_v>(key_p));
133 }
134};
135
136/*! \brief Class that select the operation to do differently if we are doing Mesh to particle (m2p) or particle to mesh (p2m)
137 *
138 * \tparam prp_g property of the grid to interpolate
139 * \tparam prp_v property of the vector to interpolate
140 * \param M2P or P2M
141 *
142 */
143template<unsigned int np, unsigned int prp_g, unsigned int prp_v>
144struct inte_template<np,prp_g,prp_v,inte_m2p>
145{
146 /*! \brief Evaluate the interpolation
147 *
148 * \tparam grid type of grid
149 * \tparam vector type of vector
150 * \tparam iterator type of the iterator
151 * \tparam key_type type of the key
152 *
153 * \param gd grid for interpolation
154 * \param vd vector for interpolation
155 * \param k_dist grid key grid point for interpolation
156 * \param key_p particle for interpolation
157 * \param a_int interpolation coefficients pre-calculated
158 * \param key indicate which pre-calculated coefficient we have to use
159 *
160 */
161 template<unsigned int np_a_int, typename grid, typename vector, typename iterator> inline static void value(grid & gd,
162 vector & vd,
163 const grid_dist_lin_dx & k_dist,
164 iterator & key_p,
165 typename vector::stype (& a_int)[np_a_int],
166 const size_t & key)
167 {
168 mul_inte<typename std::remove_reference<decltype(gd.template get<prp_g>(k_dist))>::type>::value(vd.template getProp<prp_v>(key_p),a_int[key],gd.template get<prp_g>(k_dist));
169 }
170};
171
172/*! \brief Calculate aint
173 *
174 * This class store
175 *
176 */
177template<unsigned int dims, typename vector, unsigned int np>
178struct calculate_aint
179{
180 /*! \brief Calculate the coefficients of the interpolation a_int for one particles
181 * having the 1D kernel values
182 *
183 * \param sz size of stencil for the interpolation
184 * \param a_int coefficients on the stencil points
185 * \param a coefficients of the kernel for each direction, consider
186 * that for 3 dimensions the kernel is the multiplication
187 * the 1D kernel on each direction. The "a" array store the
188 * calculated coefficient of the 1D kernel on each direction
189 *
190 */
191 static inline void value(size_t (& sz)[vector::dims],
192 typename vector::stype a_int[openfpm::math::pow(np,vector::dims)],
193 typename vector::stype (& a)[vector::dims][np])
194 {
195 grid_sm<vector::dims,void> gs(sz);
196 grid_key_dx_iterator<vector::dims> kit(gs);
197
198 size_t s = 0;
199 while (kit.isNext())
200 {
201 auto key = kit.get();
202
203 a_int[s] = 1;
204
205 for (size_t i = 0 ; i < vector::dims ; i++)
206 a_int[s] *= a[i][key.get(i)];
207
208 s++;
209 ++kit;
210 }
211 }
212};
213
214/*! \brief Calculate aint 2D
215 *
216 *
217 */
218template<typename vector, unsigned int np>
219struct calculate_aint<2,vector,np>
220{
221
222 /*! \brief Calculate the coefficients of the interpolation a_int for one particles
223 * having the 1D kernel values
224 *
225 * \param sz size of stencil for the interpolation
226 * \param a_int coefficients on the stencil points
227 * \param a coefficients of the kernel for each direction, consider
228 * that for 3 dimensions the kernel is the multiplication
229 * the 1D kernel on each direction. The "a" array store the
230 * calculated coefficient of the 1D kernel on each direction
231 *
232 */
233 static inline void value(size_t (& sz)[vector::dims],
234 typename vector::stype a_int[openfpm::math::pow(np,vector::dims)],
235 typename vector::stype (& a)[vector::dims][np])
236 {
237 size_t s = 0;
238 for (size_t i = 0 ; i < np ; i++)
239 {
240 for (size_t j = 0 ; j < np ; j++)
241 {
242 a_int[s] = a[0][j]*a[1][i];
243
244 s++;
245 }
246 }
247 }
248};
249
250/*! \brief Calculate aint
251 *
252 *
253 */
254template<typename vector, unsigned int np>
255struct calculate_aint<3,vector,np>
256{
257 /*! \brief Calculate the coefficients of the interpolation a_int for one particles
258 * having the 1D kernel values
259 *
260 * \param sz size of stencil for the interpolation
261 * \param a_int coefficients on the stencil points
262 * \param a coefficients of the kernel for each direction, consider
263 * that for 3 dimensions the kernel is the multiplication
264 * the 1D kernel on each direction. The "a" array store the
265 * calculated coefficient of the 1D kernel on each direction
266 *
267 */
268 static inline void value(size_t (& sz)[vector::dims],
269 typename vector::stype a_int[openfpm::math::pow(np,vector::dims)],
270 typename vector::stype (& a)[vector::dims][np])
271 {
272 size_t s = 0;
273 for (size_t i = 0 ; i < np ; i++)
274 {
275 for (size_t j = 0 ; j < np ; j++)
276 {
277 for (size_t k = 0 ; k < np ; k++)
278 {
279 a_int[s] = a[0][k]*a[1][j]*a[2][i];
280
281 s++;
282 }
283 }
284 }
285 }
286};
287
288/*! \brief return the sub-domain where this particle must be interpolated
289 *
290 * \param p particle position
291 *
292 * \return the sub-domain id
293 *
294 */
295template<typename vector, typename grid>
296inline size_t getSub(Point<vector::dims,typename vector::stype> & p,
297 const CellList<vector::dims,typename vector::stype,Mem_fast<>,shift<vector::dims,typename vector::stype>> & geo_cell,
298 grid & gd)
299{
300 size_t cell = geo_cell.getCell(p);
301
302 for (size_t i = 0 ; i < geo_cell.getNelements(cell) ; i++)
303 {
304 size_t ns = geo_cell.get(cell,i);
305
306 if (gd.getDecomposition().getSubDomain(ns).isInside(p))
307 return ns;
308 }
309
310 // If we end up here it mean that the particle decomposition and the grid decomposition are at machine precision level
311 // different. To recover we shift the particle of a machine precision correction inside the domain.
312
313 typename vector::stype dx[vector::dims];
314 typename vector::stype best_dx[vector::dims];
315 typename vector::stype best_tot_dx;
316 long int best_sub;
317
318 Box<vector::dims,typename vector::stype> sub_dom = gd.getDecomposition().getSubDomain(0);
319
320 for (size_t j = 0 ; j < vector::dims ; j++)
321 {
322 if (sub_dom.getLow(j) > p.get(j))
323 {dx[j] = 2*(sub_dom.getLow(j) - p.get(j));}
324 else if (sub_dom.getHigh(j) < p.get(j))
325 {dx[j] = 2*(sub_dom.getHigh(j) - p.get(j));}
326 else {dx[j] = 0;}
327 }
328
329 typename vector::stype tot_dx = 0.0;
330
331 for (size_t j = 0 ; j < vector::dims ; j++)
332 {tot_dx += fabs(dx[j]);}
333
334 best_tot_dx = tot_dx;
335 for (size_t j = 0 ; j < vector::dims ; j++)
336 {best_dx[j] = dx[j];}
337 best_sub = 0;
338
339 for (size_t i = 1 ; i < gd.getDecomposition().getNSubDomain() ; i++)
340 {
341 Box<vector::dims,typename vector::stype> sub_dom = gd.getDecomposition().getSubDomain(i);
342
343 for (size_t j = 0 ; j < vector::dims ; j++)
344 {
345 if (sub_dom.getLow(j) > p.get(j))
346 {dx[j] = 2*(sub_dom.getLow(j) - p.get(j));}
347 else if (sub_dom.getHigh(j) < p.get(j))
348 {dx[j] = 2*(sub_dom.getHigh(j) - p.get(j));}
349 else {dx[j] = 0;}
350 }
351
352 typename vector::stype tot_dx = 0.0;
353
354 for (size_t j = 0 ; j < vector::dims ; j++)
355 {tot_dx += fabs(dx[j]);}
356
357 if (tot_dx < best_tot_dx)
358 {
359 best_tot_dx = tot_dx;
360 for (size_t j = 0 ; j < vector::dims ; j++)
361 {best_dx[j] = dx[j];}
362 best_sub = i;
363 }
364 }
365
366 // shift xp
367
368 for (size_t j = 0 ; j < vector::dims ; j++)
369 {p.get(j) += best_dx[j];}
370
371 return best_sub;
372}
373
374/*! \brief calculate the interpolation for one point
375 *
376 * \tparam vector of particles
377 * \tparam kernel type
378 *
379 */
380template<typename vector,typename kernel>
381struct inte_calc_impl
382{
383 //! Type of the calculations
384 typedef typename vector::stype arr_type;
385
386 /*! \brief M2P or P2M for one particle
387 *
388 * \tparam prp_g property to interpolate from(M2P) or to(P2M) for grid
389 * \tparam prp_v property to interpolate to(M2P) or from(P2M) for vector
390 * \tparam m2p_or_p2m mesh to particle or mesh to particle interpolation
391 *
392 * \param it iterator used to retrieve the particle p for interpolation
393 * \param vd vector of particles
394 * \param domain simulation domain
395 * \param ip index of the grid on each direction (1D) used for interpolation
396 * \param gd interpolation grid
397 * \param dx spacing on each direction
398 * \param xp Point that store the position of xp
399 * \param a_int coefficients on the stencil points
400 * \param a coefficients of the kernel for each direction, consider
401 * that for 3 dimensions the kernel is the multiplication
402 * the 1D kernel on each direction. The "a" array store the
403 * calculated coefficient of the 1D kernel on each direction
404 * \param x position of
405 * \param sz grid size
406 * \param geo_cell cell list to convert particle position into sub-domain id
407 * \param offsets array where are stored the linearized offset of the
408 * kernel stencil for each local grid (sub-domain)
409 *
410 */
411 template<unsigned int prp_g, unsigned int prp_v, unsigned int m2p_or_p2m, unsigned int np_a_int, typename grid>
412 static inline void inte_calc(const vect_dist_key_dx & key_p,
413 vector & vd,
414 const Box<vector::dims,typename vector::stype> & domain,
415 int (& ip)[vector::dims][kernel::np],
416 grid & gd,
417 const typename vector::stype (& dx)[vector::dims],
418 typename vector::stype (& xp)[vector::dims],
419 typename vector::stype (& a_int)[np_a_int],
420 typename vector::stype (& a)[vector::dims][kernel::np],
421 typename vector::stype (& x)[vector::dims][kernel::np],
422 size_t (& sz)[vector::dims],
423 const CellList<vector::dims,typename vector::stype,Mem_fast<>,shift<vector::dims,typename vector::stype>> & geo_cell,
424 openfpm::vector<agg_arr<openfpm::math::pow(kernel::np,vector::dims)>> & offsets)
425 {
426 Point<vector::dims,typename vector::stype> p = vd.getPos(key_p);
427
428 // On which sub-domain we interpolate the particle
429 size_t sub = getSub<vector>(p,geo_cell,gd);
430
431 typename vector::stype x0[vector::dims];
432
433 // calculate the position of the particle in the grid
434 // coordinated
435 for (size_t i = 0 ; i < vector::dims ; i++)
436 {x0[i] = (p.get(i)-domain.getLow(i))*dx[i];}
437
438 // convert into integer
439 for (size_t i = 0 ; i < vector::dims ; i++)
440 {ip[i][0] = (int)x0[i];}
441
442 // convert the global grid position into local grid position
443 grid_key_dx<vector::dims> base;
444
445 for (size_t i = 0 ; i < vector::dims ; i++)
446 {base.set_d(i,ip[i][0] - gd.getLocalGridsInfo().get(sub).origin.get(i) - (long int)kernel::np/2 + 1);}
447
448 // convenient grid of points
449
450 for (size_t j = 0 ; j < kernel::np-1 ; j++)
451 {
452 for (size_t i = 0 ; i < vector::dims ; i++)
453 {ip[i][j+1] = (int)ip[i][j]+1;}
454 }
455
456 for (size_t i = 0 ; i < vector::dims ; i++)
457 xp[i] = x0[i] - ip[i][0];
458
459 for (long int j = 0 ; j < kernel::np ; j++)
460 {
461 for (size_t i = 0 ; i < vector::dims ; i++)
462 {x[i][j] = - xp[i] + typename vector::stype((long int)j - (long int)kernel::np/2 + 1);}
463 }
464
465 for (size_t j = 0 ; j < kernel::np ; j++)
466 {
467 for (size_t i = 0 ; i < vector::dims ; i++)
468 {a[i][j] = kernel::value(x[i][j],j);}
469 }
470
471 calculate_aint<vector::dims,vector,kernel::np>::value(sz,a_int,a);
472
473 grid_dist_lin_dx k_dist_lin;
474 k_dist_lin.setSub(sub);
475
476 size_t k = 0;
477 auto & gs_info = gd.get_loc_grid(sub).getGrid();
478
479 size_t lin_base = gs_info.LinId(base);
480
481 for (size_t i = 0 ; i < openfpm::math::pow(kernel::np,vector::dims) ; i++)
482 {
483 size_t lin = offsets.get(sub).ele[k] + lin_base;
484 k_dist_lin.getKeyRef() = lin;
485
486 inte_template<kernel::np,prp_g,prp_v,m2p_or_p2m>::value(gd,vd,k_dist_lin,key_p,a_int,k);
487
488 k++;
489 }
490 }
491};
492
493/*! \brief Main class for interpolation Particle to mest p2m and Mesh to particle m2p
494 *
495 * This function is the main class to interpolate from particle to mesh and mesh to particle
496 *
497 * \tparam vector type of vector for interpolation
498 * \tparam grid type of grid for interpolation
499 * \tparam interpolation kernel
500 *
501 */
502template<typename vector,typename grid, typename kernel>
503class interpolate
504{
505 //! Cell list used to convert particles position to sub-domain
506 CellList<vector::dims,typename vector::stype,Mem_fast<>,shift<vector::dims,typename vector::stype>> geo_cell;
507
508 /*! Structure to order boxes by volume
509 *
510 *
511 *
512 */
513 struct Box_vol
514 {
515 //! Box
516 Box<vector::dims,size_t> bv;
517
518 //! Calculated volume
519 size_t vol;
520
521 /*! \brief operator to reorder
522 *
523 * \param bv box to compare with
524 *
525 * \return true if bv has volume bigger volume
526 *
527 */
528 bool operator<(const Box_vol & bv)
529 {
530 return vol < bv.vol;
531 }
532 };
533
534 //! particles
535 vector & vd;
536
537 //! grid
538 grid & gd;
539
540 //! Type of the calculations
541 typedef typename vector::stype arr_type;
542
543 //! the offset for each sub-sub-domain
544 openfpm::vector<agg_arr<openfpm::math::pow(kernel::np,vector::dims)>> offsets;
545
546 //! kernel size
547 size_t sz[vector::dims];
548
549 //! grid spacing
550 typename vector::stype dx[vector::dims];
551
552 //! Simulation domain
553 Box<vector::dims,typename vector::stype> domain;
554
555 /*! \brief It calculate the interpolation stencil offsets
556 *
557 * \param offsets array where to store the linearized offset of the
558 * kernel stencil for each local-grid (sub-domain)
559 * \param sz kernel stencil points in each direction
560 *
561 */
562 void calculate_the_offsets(openfpm::vector<agg_arr<openfpm::math::pow(kernel::np,vector::dims)>> & offsets, size_t (& sz)[vector::dims])
563 {
564 offsets.resize(gd.getN_loc_grid());
565
566 for (size_t i = 0 ; i < offsets.size() ; i++)
567 {
568 auto & loc_grid = gd.get_loc_grid(i);
569 const grid_sm<vector::dims,void> & gs = loc_grid.getGrid();
570
571 grid_sm<vector::dims,void> gs2(sz);
572 grid_key_dx_iterator<vector::dims> kit2(gs2);
573
574 size_t k = 0;
575
576 while (kit2.isNext())
577 {
578 auto key = kit2.get();
579
580 long int lin = gs.LinId(key);
581
582 offsets.get(i).ele[k] = lin;
583
584 ++k;
585 ++kit2;
586 }
587 }
588 }
589
590public:
591
592 /*! \brief construct an interpolation object between a grid and a vector
593 *
594 * When possible and easy to do we suggest to retain the object
595 *
596 * \param vd interpolation vector
597 * \param gd interpolation grid
598 *
599 */
600 interpolate(vector & vd, grid & gd)
601 :vd(vd),gd(gd)
602 {
603 // get the processor bounding box in grid units
604 Box<vector::dims,typename vector::stype> bb = gd.getDecomposition().getProcessorBounds();
605 Box<vector::dims,typename vector::stype> bunit = gd.getDecomposition().getCellDecomposer().getCellBox();
606
607 size_t div[vector::dims];
608
609 for (size_t i = 0 ; i < vector::dims ; i++)
610 div[i] = (bb.getHigh(i) - bb.getLow(i)) / bunit.getHigh(i);
611
612 geo_cell.Initialize(bb,div);
613
614 // Now draw the domain into the cell list
615
616 auto & dec = gd.getDecomposition();
617
618 for (size_t i = 0 ; i < dec.getNSubDomain() ; i++)
619 {
620 const Box<vector::dims,typename vector::stype> & bx = dec.getSubDomain(i);
621
622 // get the cells this box span
623 const grid_key_dx<vector::dims> p1 = geo_cell.getCellGrid(bx.getP1());
624 const grid_key_dx<vector::dims> p2 = geo_cell.getCellGrid(bx.getP2());
625
626 // Get the grid and the sub-iterator
627 auto & gi = geo_cell.getGrid();
628 grid_key_dx_iterator_sub<vector::dims> g_sub(gi,p1,p2);
629
630 // add the box-id to the cell list
631 while (g_sub.isNext())
632 {
633 auto key = g_sub.get();
634 geo_cell.addCell(gi.LinId(key),i);
635 ++g_sub;
636 }
637 }
638
639 for (size_t i = 0 ; i < vector::dims ; i++)
640 {sz[i] = kernel::np;}
641
642 calculate_the_offsets(offsets,sz);
643
644 // calculate spacing
645 for (size_t i = 0 ; i < vector::dims ; i++)
646 {dx[i] = 1.0/gd.spacing(i);}
647
648 // simulation domain
649 domain = vd.getDecomposition().getDomain();
650 };
651
652 /*! \brief Interpolate particles to mesh
653 *
654 * Most of the time the particle set and the mesh are the same
655 * as the one used in the constructor. They also can be different
656 * as soon as they have the same decomposition
657 *
658 * \param gd grid or mesh
659 * \param vd particle set
660 *
661 */
662 template<unsigned int prp_v, unsigned int prp_g> void p2m(vector & vd, grid & gd)
663 {
664#ifdef SE_CLASS1
665
666 if (!vd.getDecomposition().is_equal_ng(gd.getDecomposition()) )
667 {
668 std::cerr << __FILE__ << ":" << __LINE__ << " Error: the distribution of the vector of particles" <<
669 " and the grid is different. In order to interpolate the two data structure must have the" <<
670 " same decomposition" << std::endl;
671
672 ACTION_ON_ERROR(INTERPOLATION_ERROR_OBJECT)
673 }
674
675#endif
676
677 // point position
678 typename vector::stype xp[vector::dims];
679
680 int ip[vector::dims][kernel::np];
681 typename vector::stype x[vector::dims][kernel::np];
682 typename vector::stype a[vector::dims][kernel::np];
683
684 typename vector::stype a_int[openfpm::math::pow(kernel::np,vector::dims)];
685
686 auto it = vd.getDomainIterator();
687
688 while (it.isNext() == true)
689 {
690 auto key_p = it.get();
691
692 inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_p2m,openfpm::math::pow(kernel::np,vector::dims)>(key_p,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
693
694 ++it;
695 }
696 }
697
698 /*! \brief Interpolate mesh to particle
699 *
700 * Most of the time the particle set and the mesh are the same
701 * as the one used in the constructor. They also can be different
702 * as soon as they have the same decomposition
703 *
704 * \param gd grid or mesh
705 * \param vd particle set
706 *
707 */
708 template<unsigned int prp_g, unsigned int prp_v> void m2p(grid & gd, vector & vd)
709 {
710#ifdef SE_CLASS1
711
712 if (!vd.getDecomposition().is_equal_ng(gd.getDecomposition()) )
713 {
714 std::cerr << __FILE__ << ":" << __LINE__ << " Error: the distribution of the vector of particles" <<
715 " and the grid is different. In order to interpolate the two data structure must have the" <<
716 " same decomposition" << std::endl;
717
718 ACTION_ON_ERROR(INTERPOLATION_ERROR_OBJECT)
719 }
720
721#endif
722
723 // point position
724 typename vector::stype xp[vector::dims];
725
726 int ip[vector::dims][kernel::np];
727 typename vector::stype x[vector::dims][kernel::np];
728 typename vector::stype a[vector::dims][kernel::np];
729
730 // grid_cpu<vector::dims,aggregate<typename vector::stype>> a_int(sz);
731 // a_int.setMemory();
732 typename vector::stype a_int[openfpm::math::pow(kernel::np,vector::dims)];
733
734 auto it = vd.getDomainIterator();
735
736 while (it.isNext() == true)
737 {
738 auto key_p = it.get();
739
740 inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_m2p,openfpm::math::pow(kernel::np,vector::dims)>(key_p,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
741
742 ++it;
743 }
744 }
745
746
747 /*! \brief Interpolate particles to mesh
748 *
749 * Most of the time the particle set and the mesh are the same
750 * as the one used in the constructor. They also can be different
751 * as soon as they have the same decomposition
752 *
753 * \param gd grid or mesh
754 * \param vd particle set
755 * \param p particle
756 *
757 */
758 template<unsigned int prp_v, unsigned int prp_g> inline void p2m(vector & vd, grid & gd,const vect_dist_key_dx & p)
759 {
760#ifdef SE_CLASS1
761
762 if (!vd.getDecomposition().is_equal_ng(gd.getDecomposition()) )
763 {
764 std::cerr << __FILE__ << ":" << __LINE__ << " Error: the distribution of the vector of particles" <<
765 " and the grid is different. In order to interpolate the two data structure must have the" <<
766 " same decomposition" << std::endl;
767
768 ACTION_ON_ERROR(INTERPOLATION_ERROR_OBJECT)
769 }
770
771#endif
772
773 // point position
774 typename vector::stype xp[vector::dims];
775
776 int ip[vector::dims][kernel::np];
777 typename vector::stype x[vector::dims][kernel::np];
778 typename vector::stype a[vector::dims][kernel::np];
779
780 typename vector::stype a_int[openfpm::math::pow(kernel::np,vector::dims)];
781
782 /* coverty[uninit_use_in_call] */
783 inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_p2m,openfpm::math::pow(kernel::np,vector::dims)>(p,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
784 }
785
786 /*! \brief Interpolate mesh to particle
787 *
788 * Most of the time the particle set and the mesh are the same
789 * as the one used in the constructor. They also can be different
790 * as soon as they have the same decomposition
791 *
792 * \param gd grid or mesh
793 * \param vd particle set
794 * \param p particle
795 *
796 */
797 template<unsigned int prp_g, unsigned int prp_v> inline void m2p(grid & gd, vector & vd, const vect_dist_key_dx & p)
798 {
799#ifdef SE_CLASS1
800
801 if (!vd.getDecomposition().is_equal_ng(gd.getDecomposition()) )
802 {
803 std::cerr << __FILE__ << ":" << __LINE__ << " Error: the distribution of the vector of particles" <<
804 " and the grid is different. In order to interpolate the two data structure must have the" <<
805 " same decomposition" << std::endl;
806
807 ACTION_ON_ERROR(INTERPOLATION_ERROR_OBJECT)
808 }
809
810#endif
811
812 // point position
813 typename vector::stype xp[vector::dims];
814
815 int ip[vector::dims][kernel::np];
816 typename vector::stype x[vector::dims][kernel::np];
817 typename vector::stype a[vector::dims][kernel::np];
818
819 // grid_cpu<vector::dims,aggregate<typename vector::stype>> a_int(sz);
820 // a_int.setMemory();
821 typename vector::stype a_int[openfpm::math::pow(kernel::np,vector::dims)];
822
823 inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_m2p,openfpm::math::pow(kernel::np,vector::dims)>(p,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
824 }
825
826};
827
828#endif /* OPENFPM_NUMERICS_SRC_INTERPOLATION_INTERPOLATION_HPP_ */
829