1/*
2 * staggered_util.hpp
3 *
4 * Created on: Aug 19, 2015
5 * Author: i-bird
6 */
7
8#ifndef SRC_GRID_STAGGERED_DIST_GRID_UTIL_HPP_
9#define SRC_GRID_STAGGERED_DIST_GRID_UTIL_HPP_
10
11#include "util/common.hpp"
12#include "VTKWriter/VTKWriter.hpp"
13#include "util/convert.hpp"
14
15
16/*! \brief write a property that has attributes
17 *
18 * \tparam ele object we are writing
19 * \tparam vtk vtk writer
20 * \tparam true in case the basic object has attributes
21 *
22 */
23template<typename ele, typename vtk, bool has_attributes=has_attributes<ele>::value>
24struct vtk_write
25{
26 /*! \brief Add the grid with attributes name
27 *
28 * \param vtk_w VTK writer
29 * \param output where to write
30 * \param i property to write
31 *
32 */
33 vtk_write(vtk vtk_w, const std::string output, const size_t i)
34 {
35 vtk_w.write(output + "_" + ele::attributes::name[i] + ".vtk",ele::attributes::name[i]);
36 }
37};
38
39/*! \brief Add to the vtk writer the key
40 *
41 * \tparam ele object we are writing
42 * \tparam vtk vtk writer
43 * \tparam false in case the basic object has not attributes
44 *
45 */
46template<typename ele, typename vtk>
47struct vtk_write<ele,vtk,false>
48{
49 /*! \brief Add the grid with attributes name
50 *
51 * \param vtk_w VTK writer
52 * \param output where to write
53 * \param i property to write
54 *
55 */
56 vtk_write(vtk vtk_w, const std::string output, const size_t i)
57 {
58 vtk_w.write(output + "_" + std::to_string(i) + ".vtk","attr" + std::to_string(i));
59 }
60};
61
62
63/*! \brief Classes to get the number of components of the properties
64 *
65 */
66template<typename T>
67struct extends
68{
69 /*! \brief Scalar case
70 *
71 * \return 1 component
72 *
73 */
74 static inline size_t mul()
75 {
76 return 1;
77 }
78
79 /*! \brief Dimensionality
80 *
81 * \return 0
82 *
83 */
84 static inline size_t dim()
85 {
86 return 0;
87 }
88};
89
90//! Partial specialization for N=1 1D-Array
91template<typename T,size_t N1>
92struct extends<T[N1]>
93{
94 /*! \brief Vector case return N1 component
95 *
96 * \return N1
97 *
98 */
99 static inline size_t mul()
100 {
101 return N1;
102 }
103
104 /*! Dimensionality 1
105 *
106 * \return 1
107 *
108 */
109 static inline size_t dim()
110 {
111 return 1;
112 }
113};
114
115//! Partial specialization for N=2 2D-Array
116template<typename T,size_t N1,size_t N2>
117struct extends<T[N1][N2]>
118{
119 /*! \brief Matrix case return N1*N2 component
120 *
121 * \return N1*N2
122 *
123 */
124 static inline size_t mul()
125 {
126 return N1 * N2;
127 }
128
129 /*! Dimensionality 2
130 *
131 * \return 2
132 *
133 */
134 static inline size_t dim()
135 {
136 return 2;
137 }
138};
139
140//! Partial specialization for N=3
141template<typename T,size_t N1,size_t N2,size_t N3>
142struct extends<T[N1][N2][N3]>
143{
144 //! number of elements
145 static inline size_t mul()
146 {
147 return N1 * N2 * N3;
148 }
149
150 /*! number of indexes
151 *
152 * \return 3
153 *
154 */
155 static inline size_t dim()
156 {
157 return 3;
158 }
159};
160
161//! Partial specialization for N=4
162template<typename T,size_t N1,size_t N2,size_t N3,size_t N4>
163struct extends<T[N1][N2][N3][N4]>
164{
165 //! number of elements
166 static inline size_t mul()
167 {
168 return N1 * N2 * N3 * N4;
169 }
170
171 /*! number of indexes
172 *
173 * \return 4
174 *
175 */
176 static inline size_t dim()
177 {
178 return 4;
179 }
180};
181
182//! Partial specialization for N=5
183template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5>
184struct extends<T[N1][N2][N3][N4][N5]>
185{
186 //! number of elements
187 static inline size_t mul()
188 {
189 return N1 * N2 * N3 * N4 * N5;
190 }
191
192 /*! number of indexes
193 *
194 * \return 5
195 *
196 */
197 static inline size_t dim()
198 {
199 return 5;
200 }
201};
202
203//! Partial specialization for N=6
204template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6>
205struct extends<T[N1][N2][N3][N4][N5][N6]>
206{
207 //! number of elements
208 static inline size_t mul()
209 {
210 return N1 * N2 * N3 * N4 * N5 * N6;
211 }
212
213 /*! number of indexes
214 *
215 * \return 6
216 *
217 */
218 static inline size_t dim()
219 {
220 return 6;
221 }
222};
223
224//! Partial specialization for N=7
225template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6, size_t N7>
226struct extends<T[N1][N2][N3][N4][N5][N6][N7]>
227{
228 //! number of elements
229 static inline size_t mul()
230 {
231 return N1 * N2 * N3 * N4 * N5 * N6 * N7;
232 }
233
234 /*! number of indexes
235 *
236 * \return 7
237 *
238 */
239 static inline size_t dim()
240 {
241 return 7;
242 }
243};
244
245//! Partial specialization for N=8
246template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6, size_t N7, size_t N8>
247struct extends<T[N1][N2][N3][N4][N5][N6][N7][N8]>
248{
249 /*! number of elements
250 *
251 * \return the number of elements as N1*N2*N3*.........
252 *
253 */
254 static inline size_t mul()
255 {
256 return N1 * N2 * N3 * N4 * N5 * N6 * N7 * N8;
257 }
258
259 /*! number of indexes
260 *
261 * \return 8
262 *
263 */
264 static inline size_t dim()
265 {
266 return 8;
267 }
268};
269
270//! Partial specialization for N=9
271template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6, size_t N7, size_t N8, size_t N9>
272struct extends<T[N1][N2][N3][N4][N5][N6][N7][N8][N9]>
273{
274 /*! number of elements
275 *
276 * \return the number of elements as N1*N2*N3*.........
277 *
278 */
279 static inline size_t mul()
280 {
281 return N1 * N2 * N3 * N4 * N5 * N6 * N7 * N8 * N9;
282 }
283
284 /*! number of indexes
285 *
286 * \return 9
287 *
288 */
289 static inline size_t dim()
290 {
291 return 9;
292 }
293};
294
295//! Partial specialization for N=10
296template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6, size_t N7, size_t N8, size_t N9, size_t N10>
297struct extends<T[N1][N2][N3][N4][N5][N6][N7][N8][N9][N10]>
298{
299 /*! number of elements
300 *
301 * \return the number of elements as N1*N2*N3*.........
302 *
303 */
304 static inline size_t mul()
305 {
306 return N1 * N2 * N3 * N4 * N5 * N6 * N7 * N8 * N9 * N10;
307 }
308
309 /*! number of indexes
310 *
311 * \return 10
312 *
313 */
314 static inline size_t dim()
315 {
316 return 10;
317 }
318};
319
320///////////////////// Copy grid extends
321
322/*! \brief Classes to copy each component into a grid and add to the VTKWriter the grid
323 *
324 * \param T property to write
325 * \param dim dimansionality
326 * \param St type of space
327 * \param VTK VTK writer
328 *
329 */
330template<typename T>
331struct write_stag
332{
333 /*! \brief write the staggered grid
334 *
335 * \tparam p_val property we are going to write
336 * \tparam sg staggered grid type
337 * \tparam v_g vector of grids
338 *
339 * \param st_g staggered grid
340 * \param v_g vector of grids
341 * \param lg local grid of the staggered grid we are writing
342 *
343 */
344 template<unsigned int p_val, typename sg, typename v_g> static inline void write(sg & st_g, v_g & vg,size_t lg)
345 {
346 // Add a grid;
347 vg.add();
348 size_t k = vg.size() - 1;
349
350 // Get the source and destination grid
351 auto & g_src = st_g.get_loc_grid(lg);
352 auto & g_dst = vg.get(k);
353
354 // Set dimensions and memory
355 g_dst.resize(g_src.getGrid().getSize());
356
357 // copy
358
359 auto it = vg.get(k).getIterator();
360
361 while(it.isNext())
362 {
363 g_dst.template get<0>(it.get()) = g_src.template get<p_val>(it.get());
364
365 ++it;
366 }
367 }
368};
369
370/*! \brief for each component add a grid fill it, and add to the VTK writer
371 *
372 * \param T Property to copy
373 * \param N1 number of components
374 *
375 */
376template<typename T,size_t N1>
377struct write_stag<T[N1]>
378{
379 /*! \brief write the staggered grid
380 *
381 * \tparam p_val property we are going to write
382 * \tparam sg staggered grid type
383 * \tparam v_g vector of grids
384 *
385 * \param st_g staggered grid
386 * \param v_g vector of grids
387 * \param lg local grid of the staggered grid we are writing
388 *
389 */
390 template<unsigned int p_val, typename sg, typename v_g> static inline void write(sg & st_g, v_g & vg,size_t lg)
391 {
392 for (size_t i = 0 ; i < N1 ; i++)
393 {
394 // Add a grid;
395 vg.add();
396 size_t k = vg.size() - 1;
397
398 // Get the source and destination grid
399 auto & g_src = st_g.get_loc_grid(lg);
400 auto & g_dst = vg.get(k);
401
402 // Set dimensions and memory
403 g_dst.resize(g_src.getGrid().getSize());
404
405 auto it = vg.get(k).getIterator();
406
407 while(it.isNext())
408 {
409 g_dst.template get<0>(it.get()) = g_src.template get<p_val>(it.get())[i];
410
411 ++it;
412 }
413 }
414 }
415};
416
417//! Partial specialization for N=2 2D-Array
418template<typename T,size_t N1,size_t N2>
419struct write_stag<T[N1][N2]>
420{
421 /*! \brief write the staggered grid
422 *
423 * \tparam p_val property we are going to write
424 * \tparam sg staggered grid type
425 * \tparam v_g vector of grids
426 *
427 * \param st_g staggered grid
428 * \param v_g vector of grids
429 * \param lg local grid of the staggered grid we are writing
430 *
431 */
432 template<unsigned int p_val, typename sg, typename v_g> static inline void write(sg & st_g, v_g & vg,size_t lg)
433 {
434 for (size_t i = 0 ; i < N1 ; i++)
435 {
436 for (size_t j = 0 ; j < N2 ; j++)
437 {
438 // Add a grid;
439 vg.add();
440 size_t k = vg.size() - 1;
441
442 // Set dimensions and memory
443 vg.get(k).resize(st_g.get_loc_grid(lg).getGrid().getSize());
444
445 // copy
446 auto & g_src = st_g.get_loc_grid(lg);
447 auto & g_dst = vg.get(k);
448 auto it = vg.get(k).getIterator();
449
450 while(it.isNext())
451 {
452 g_dst.template get<0>(it.get()) = g_src.template get<p_val>(it.get())[i][j];
453
454 ++it;
455 }
456 }
457 }
458 }
459};
460
461///////////////////// Staggered default positioning ////////////////////////
462
463/*! \brief this class is a functor for "for_each" algorithm
464 *
465 * For each element of the boost::vector the operator() is called.
466 * Is mainly used to produce a default position vector for each
467 * property
468 *
469 * \tparam dim dimensionality
470 * \tparam v boost::fusion::vector of properties
471 * \tparam has_posMask case when v has a position mask
472 *
473 */
474
475template<unsigned int dim, typename v, bool has_pM = has_posMask<v>::value>
476class stag_set_position
477{
478 //! vector containing the position of the properties in the cells (staggered properties are staggered)
479 // within the cell
480 openfpm::vector<comb<dim>> (& pos_prp)[boost::fusion::result_of::size<v>::type::value];
481
482public:
483
484 /*! \brief Constructor
485 *
486 * \param vector of the staggered position (It is going to be filled by this class)
487 *
488 */
489 stag_set_position( openfpm::vector<comb<dim>> (& pos_prp)[boost::fusion::result_of::size<v>::type::value])
490 :pos_prp(pos_prp)
491 {}
492
493 /*! It calculate the staggered position for every property
494 *
495 * \param t property
496 *
497 */
498 template<typename T>
499 void operator()(T& t) const
500 {
501 // This is the type of the object we have to copy
502 typedef typename boost::mpl::at<v,typename boost::mpl::int_<T::value>>::type prop;
503
504 bool val = prop::stag_pos_mask[T::value];
505
506 if (val == false)
507 return;
508
509 // Dimension of the object
510 size_t dim_prp = extends<prop>::dim();
511
512 // It is a scalar
513 if (dim_prp == 0)
514 {
515 comb<dim> c;
516 c.zero();
517
518 // It stay in the center
519 pos_prp[T::value].add(c);
520 }
521 else if (dim_prp == 1)
522 {
523 // It stay on the object of dimension dim-1 (Negative part)
524 for (size_t i = 0 ; i < dim ; i++)
525 {
526 comb<dim> c;
527 c.zero();
528 c.value(i) = -1;
529
530 pos_prp[T::value].add(c);
531 }
532 }
533 else if (dim_prp == 2)
534 {
535 // Create an hypercube object
536 HyperCube<dim> hyp;
537
538 // Diagonal part live in
539 for (size_t i = 0 ; i < dim ; i++)
540 {
541 comb<dim> c1 = pos_prp[T::value-1].get(i);
542 for (size_t j = 0 ; j < dim ; j++)
543 {
544 comb<dim> c2;
545 c2.zero();
546 c2.value(i) = 1;
547
548 comb<dim> c_res = (c1 + c2) & 0x1;
549
550 pos_prp[T::value].add(c_res);
551 }
552 }
553 }
554 else if (dim_prp > 2)
555 {
556 std::cerr << __FILE__ << ":" << __LINE__ << " Tensor of rank bigger than 2 are not supported";
557 }
558 }
559};
560
561///////////////////// Staggered default positioning ////////////////////////
562
563/*! \brief this class is a functor for "for_each" algorithm
564 *
565 * For each element of the boost::vector the operator() is called.
566 * Is mainly used to produce a default position vector for each
567 * property
568 *
569 * \tparam vector of properties
570 *
571 */
572
573template<unsigned int dim, typename v>
574class stag_set_position<dim,v,false>
575{
576private:
577
578 //! vector containing the position of the properties in the cells (staggered properties are staggered)
579 // within the cell
580 openfpm::vector<comb<dim>> (& pos_prp)[boost::fusion::result_of::size<v>::type::value];
581
582
583public:
584
585 /*! \brief Constructor
586 *
587 * \param vector of the staggered position (It is going to be filled by this class)
588 *
589 */
590 stag_set_position( openfpm::vector<comb<dim>> (& pos_prp)[boost::fusion::result_of::size<v>::type::value])
591 :pos_prp(pos_prp)
592 {}
593
594 /*! It calculate the staggered position for every property
595 *
596 * \param t property
597 *
598 */
599 template<typename T>
600 void operator()(T& t) const
601 {
602 // This is the type of the object we have to copy
603 typedef typename boost::mpl::at<v,typename boost::mpl::int_<T::value>>::type prop;
604
605 // Dimension of the object
606 size_t dim_prp = extends<prop>::dim();
607
608 // It is a scalar
609 if (dim_prp == 0)
610 {
611 comb<dim> c;
612 c.zero();
613
614 // It stay in the center
615 pos_prp[T::value].add(c);
616 }
617 else if (dim_prp == 1)
618 {
619 // It stay on the object of dimension dim-1 (Negative part)
620 for (size_t i = 0 ; i < dim ; i++)
621 {
622 comb<dim> c;
623 c.zero();
624 c.getComb()[i] = -1;
625
626 pos_prp[T::value].add(c);
627 }
628 }
629 else if (dim_prp == 2)
630 {
631 // Diagonal part live in
632 for (size_t i = 0 ; i < dim ; i++)
633 {
634 comb<dim> c1 = pos_prp[T::value-1].get(i);
635 for (size_t j = 0 ; j < dim ; j++)
636 {
637 comb<dim> c2;
638 c2.zero();
639 c2.getComb()[j] = 1;
640
641 comb<dim> c_res = (c2 + c1).flip();
642
643 pos_prp[T::value].add(c_res);
644 }
645 }
646 }
647 else if (dim_prp > 2)
648 {
649 std::cerr << __FILE__ << ":" << __LINE__ << " Tensor of rank bigger than 2 are not supported";
650 }
651 }
652};
653
654/*! \brief It create separated grid for each properties to write them into a file
655 *
656 * \tparam dim dimensionality of the grids
657 * \tparam obj type object to print, must be in OpenFPM format
658 *
659 */
660template<unsigned int dim, typename st_grid, typename St>
661class stag_create_and_add_grid
662{
663
664 size_t p_id;
665
666 // staggered grid to write
667 st_grid & st_g;
668
669public:
670
671 /*! \brief Constructor
672 *
673 * \param st_g staggered grid
674 * \param p_id process id
675 *
676 */
677 stag_create_and_add_grid(st_grid & st_g, size_t p_id)
678 :p_id(p_id),st_g(st_g)
679 {}
680
681 template<unsigned int p_val> void out_normal()
682 {
683 // property type
684 typedef typename boost::mpl::at< typename st_grid::value_type::type , typename boost::mpl::int_<p_val> >::type ele;
685
686 // create an openfpm format object from the property type
687 typedef object<typename boost::fusion::vector<ele>> d_object;
688
689 VTKWriter<boost::mpl::pair<grid_cpu<dim, d_object >,St>,VECTOR_GRIDS> vtk_w;
690
691 // Create a vector of grids
692
693 openfpm::vector< grid_cpu<dim, d_object > > vg(st_g.getN_loc_grid());
694
695 // for each domain grid
696 for (size_t i = 0 ; i < vg.size() ; i++)
697 {
698 // Set dimansions and memory
699 vg.get(i).resize(st_g.get_loc_grid(i).getGrid().getSize());
700
701 auto & g_src = st_g.get_loc_grid(i);
702 auto & g_dst = vg.get(i);
703
704 auto it = vg.get(i).getIterator();
705
706 while(it.isNext())
707 {
708 object_si_d< decltype(g_src.get_o(it.get())),decltype(g_dst.get_o(it.get())) ,OBJ_ENCAP,p_val>(g_src.get_o(it.get()),g_dst.get_o(it.get()));
709
710 ++it;
711 }
712
713 Point<dim,St> offset = st_g.getOffset(i);
714 Point<dim,St> spacing = st_g.getSpacing();
715 Box<dim,size_t> dom = st_g.getDomain(i);
716
717 vtk_w.add(g_dst,offset,spacing,dom);
718 }
719
720 vtk_w.write("vtk_grids_st_" + std::to_string(p_id) + "_" + std::to_string(p_val) + ".vtk");
721 }
722
723 template<unsigned int p_val> void out_staggered()
724 {
725 // property type
726 typedef typename boost::mpl::at< typename st_grid::value_type::type , typename boost::mpl::int_<p_val> >::type ele;
727
728 // Eliminate the extends
729 typedef typename std::remove_all_extents<ele>::type r_ele;
730
731 // create an openfpm format object from the property type
732 typedef object<typename boost::fusion::vector<r_ele>> d_object;
733
734 VTKWriter<boost::mpl::pair<grid_cpu<dim, d_object >,St>,VECTOR_ST_GRIDS> vtk_w;
735
736 // Create a vector of grids
737 openfpm::vector< grid_cpu<dim, d_object > > vg;
738 vg.reserve(st_g.getN_loc_grid() * extends<ele>::mul());
739
740 size_t k = 0;
741
742 // for each domain grid
743 for (size_t i = 0 ; i < st_g.getN_loc_grid() ; i++)
744 {
745 write_stag<ele>::template write<p_val, st_grid,openfpm::vector< grid_cpu<dim, d_object > > >(st_g,vg,i);
746
747 // for each component
748 for ( ; k < vg.size() ; k++)
749 {
750 Point<dim,St> offset = st_g.getOffset(i);
751 Point<dim,St> spacing = st_g.getSpacing();
752 Box<dim,size_t> dom = st_g.getDomain(i);
753
754 vtk_w.add(i,vg.get(k),offset,spacing,dom,st_g.c_prp[p_val].get(k));
755 }
756
757 k = vg.size();
758 }
759
760 vtk_write<typename st_grid::value_type,VTKWriter<boost::mpl::pair<grid_cpu<dim, d_object >,St>,VECTOR_ST_GRIDS>> v(vtk_w,"vtk_grids_st_" + std::to_string(p_id),p_val);
761 }
762
763 //! It call the copy function for each property
764 template<typename T>
765 void operator()(T& t)
766 {
767 if (st_g.is_staggered_prop(T::value) == false)
768 out_normal<T::value>();
769 else
770 out_staggered<T::value>();
771 }
772};
773
774/*! \brief Check that the size of the iterators match
775 *
776 * It check the the boxes that the sub iterator defines has same dimensions, for example
777 * if the first sub-iterator, iterate from (1,1) to (5,3) and the second from (2,2) to (6,4)
778 * they match (2,2) to (4,6) they do not match
779 *
780 * \tparam Grid_map type of the map grid
781 * \tparam Grid_dst type of the destination grid
782 *
783 * \param it1 Iterator1
784 * \param it2 Iterator2
785 *
786 * \return true if they match
787 *
788 */
789template<typename Eqs_sys, typename it1_type, typename it2_type> bool checkIterator(const it1_type & it1, const it2_type & it2)
790{
791#ifdef SE_CLASS1
792
793 grid_key_dx<Eqs_sys::dims> it1_k = it1.getStop() - it1.getStart();
794 grid_key_dx<Eqs_sys::dims> it2_k = it2.getStop() - it2.getStart();
795
796 for (size_t i = 0 ; i < Eqs_sys::dims ; i++)
797 {
798 if (it1_k.get(i) != it2_k.get(i))
799 {
800 std::cerr << __FILE__ << ":" << __LINE__ << " error src iterator and destination iterator does not match in size\n";
801 return false;
802 }
803 }
804
805 return true;
806#else
807
808 return true;
809
810#endif
811}
812
813/*! \brief this class is a functor for "for_each" algorithm
814 *
815 * This class is a functor for "for_each" algorithm. For each
816 * element of the boost::vector the operator() is called.
817 * Is mainly used to calculate the interpolation points for each
818 * property in a staggered grid
819 *
820 * \tparam dim Dimensionality
821 * \tparam v_prp_id vector of properties id
822 * \tparam v_prp_type vector with the properties
823 *
824 */
825template<unsigned int dim, unsigned int n_prop, typename v_prp_id, typename v_prp_type>
826struct interp_points
827{
828/*#ifdef SE_CLASS3
829
830 // number of properties we are processing
831 typedef boost::mpl::size<v_prp_id> v_size_old;
832
833 typedef boost::mpl::int_<v_size_old::value-3> v_size;
834
835#else*/
836
837 // number of properties we are processing
838 typedef boost::mpl::size<v_prp_id> v_size;
839
840//#endif
841
842 // interpolation points for each property
843 openfpm::vector<std::vector<comb<dim>>> (& interp_pts)[v_size::value];
844
845 // staggered position for each property
846 const openfpm::vector<comb<dim>> (&stag_pos)[n_prop];
847
848 /*! \brief constructor
849 *
850 * It define the copy parameters.
851 *
852 * \param inter_pts array that for each property contain the interpolation points for each components
853 * \param staggered position for each property and components
854 *
855 */
856 inline interp_points(openfpm::vector<std::vector<comb<dim>>> (& interp_pts)[v_size::value],const openfpm::vector<comb<dim>> (&stag_pos)[n_prop])
857 :interp_pts(interp_pts),stag_pos(stag_pos){};
858
859 //! It call the copy function for each property
860 template<typename T>
861 inline void operator()(T& t)
862 {
863 // This is the type of the object we have to copy
864 typedef typename boost::mpl::at_c<v_prp_type,T::value>::type prp_type;
865 typedef typename boost::mpl::at<v_prp_id,T>::type p_id;
866
867 interp_pts[T::value].resize(stag_pos[p_id::value].size());
868
869 for (size_t i = 0 ; i < stag_pos[p_id::value].size() ; i++)
870 {
871 // Create the interpolation points
872 interp_pts[T::value].get(i) = SubHyperCube<dim,dim - std::rank<prp_type>::value>::getCombinations_R(stag_pos[p_id::value].get(i),0);
873
874 // interp_point are -1,0,1, map the -1 to 0 and 1 to -1
875 for (size_t j = 0 ; j < interp_pts[T::value].get(i).size() ; j++)
876 {
877 for (size_t k = 0 ; k < dim ; k++)
878 interp_pts[T::value].get(i)[j].getComb()[k] = - ((interp_pts[T::value].get(i)[j].getComb()[k] == -1)?0:interp_pts[T::value].get(i)[j].getComb()[k]);
879 }
880 }
881 }
882};
883
884#endif /* SRC_GRID_STAGGERED_DIST_GRID_UTIL_HPP_ */
885