1
2#ifndef BOX_HPP_
3#define BOX_HPP_
4
5
6#include "Space/Shape/Sphere.hpp"
7#include <boost/fusion/sequence/intrinsic/at_c.hpp>
8#include "Grid/grid_key.hpp"
9#include "memory_ly/Encap.hpp"
10#include <sstream>
11
12#define PERIODIC 1
13#define NON_PERIODIC 0
14
15/*! \brief It define if we want the upper base or the down base (Lower or upper)
16 * extreme of the interval
17 *
18 */
19
20enum Base
21{
22 UP,
23 DOWN
24};
25
26/*! \brief This class represent an N-dimensional box
27 *
28 * The box is defined by two points p2 and p1
29 *
30 \verbatim
31
32 +---------+ p2
33 | |
34 | |
35 | |
36 | |
37 | |
38 p1 +---------+
39
40 \endverbatim
41 *
42 * \tparam dim dimensionality of the Box
43 * \tparam T type of space ... double float int size_t
44 *
45 * ### Expand the box with some spacing
46 * \snippet Box_unit_tests.hpp expand the box with some spacing
47 * ### Create an enclosing box
48 * \snippet Box_unit_tests.hpp create an enclosing box
49 * ### Create the smallest boxes between several boxes
50 * \snippet Box_unit_tests.hpp Create the smallest boxes between several boxes
51 * ### Enlarge the box
52 * \snippet Box_unit_tests.hpp Enlarge the box
53 * ### Enlarge the box with fixed P1
54 * \snippet Box_unit_tests.hpp Enlarge the box with fixed P1
55 *
56 * \see SpaceBox
57 *
58 */
59template<unsigned int dim , typename T>
60class Box
61{
62public:
63
64 //! boost fusion that store the point
65 typedef boost::fusion::vector<T[dim],T[dim]> type;
66 //! type of the box
67 typedef T btype;
68
69 //! Indicate that this is a box
70 typedef int yes_is_box;
71
72 //! It store the two point bounding the box
73 type data;
74
75 //! Low point
76 static const unsigned int p1 = 0;
77 //! High point
78 static const unsigned int p2 = 1;
79 //! Maximum number of properties
80 static const unsigned int max_prop = 2;
81
82 //! dimensionality of the box
83 static const unsigned int dims = dim;
84
85 /*! \brief Intersect
86 *
87 * Intersect two boxes and return the result boxes, if the boxes does not intersect, return false
88 *
89 * \param b box to intersect with
90 * \param b_out box result of the intersection
91 *
92 * \return true if they intersect
93 *
94 */
95 __device__ __host__ bool Intersect(const Box<dim,T> & b, Box<dim,T> & b_out) const
96 {
97 // check if p1 of b is smaller than
98
99 for (size_t i = 0 ; i < dim ; i++)
100 {
101 if (getLow(i) <= b.getLow(i))
102 b_out.setLow(i,b.getLow(i));
103 else if (getLow(i) <= b.getHigh(i))
104 b_out.setLow(i,getLow(i));
105 else
106 return false;
107
108 if (getHigh(i) >= b.getHigh(i))
109 b_out.setHigh(i,b.getHigh(i));
110 else if (getHigh(i) >= b.getLow(i))
111 b_out.setHigh(i,getHigh(i));
112 else
113 return false;
114 }
115 return true;
116 }
117
118
119
120 /*! \brief Intersect
121 *
122 * Intersect two boxes and return the result boxes, if the boxes does not intersect, return false
123 *
124 * \param e_b encapsulator box to intersect with
125 * \param b_out box result of the intersection
126 *
127 * \return true if they intersect
128 *
129 */
130 template<typename Mem>
131 __device__ __host__
132 bool Intersect(const encapc<1,Box<dim,T>,Mem> & e_b, Box<dim,T> & b_out) const
133 {
134 return Intersect(Box<dim,T>(e_b),b_out);
135 }
136
137 /*!
138 *
139 * \brief Check if the sphere intersect the box
140 *
141 * \verbatim
142
143 p1 _____
144 | |
145 p. | |
146 |____|
147 p2
148
149 \endverbatim
150 * Given a point p and we search for the distance of the nearest point of the box
151 * from p. In this case the distance is p1.x-p0.x the Y dimension is alligned with p
152 * In general for each dimension we check if the point is in the interval if it is
153 * we do not accumulate, otherwise we accumulate the smallest between (p1-p0) (p2-p0).
154 *
155 * if the distance of the nearest point is smaller than the radius we have an intersection
156 *
157 *
158 *
159 * \tparam distance functor
160 * \param sphere to check the intersection
161 * \return true if intersect false otherwise
162 *
163 */
164 template <typename distance> bool Intersect(Sphere<dim,T> & sphere)
165 {
166 // distance functor
167 distance dist;
168
169 // Get the nearest point of the box from the center of the sphere
170 typename distance::ResultType distance_r = 0;
171
172 for (size_t i = 0 ; i < dim ; i++)
173 {
174
175 // if the center of the sphere on dimension i is not in the i box interval
176 // do not accumulate, otherwise accumulate from the nearest point on that
177 // dimension
178 if (boost::fusion::at_c<p1>(data)[i] < sphere.center(i))
179 {
180 // accumulate the distance from p1
181 distance_r += dist.accum_dist(sphere.center(i),boost::fusion::at_c<p1>(data)[i],i);
182 }
183 else if ( boost::fusion::at_c<p2>(data)[i] <= sphere.center(i))
184 {
185 // accumulate the distance from p2
186 distance_r += dist.accum_dist(sphere.center(i),boost::fusion::at_c<p2>(data)[i],i);
187 }
188 }
189
190 // return if there is intersection
191 return distance_r < sphere.radius();
192 }
193
194 /* \brief return the minimum distance between two boxes
195 *
196 * Consider two box in space like in the figure below. Given two points P and Q in the two boxes 1 and 2.
197 * The minimum distance, is the distance that satify min( dist(P,Q) ). There P and Q is the standard euclidean distance
198 *
199 * \verbatim b box 2
200 *
201 * \verbatim
202
203
204
205 _____ p2
206 | |
207 | | Box 2
208 |p1__|
209 /
210 /
211 / <----------- min distance
212 /
213 _____/
214 | p2|
215 | | Box 1
216 |____|
217 p1
218
219 \endverbatim
220 *
221 *
222 */
223 T min_distance(Box<dim,T> & b) const
224 {
225 Point<dim,T> dist;
226
227 for (unsigned int i = 0 ; i < dim ; i++)
228 {
229 if (getHigh(i) <= b.getLow(i))
230 {dist.get(i) = b.getLow(i) - getHigh(i);}
231 else if (b.getHigh(i) <= getLow(i))
232 {dist.get(i) = getLow(i) - b.getHigh(i);}
233 else
234 {
235 T d1 = fabs(getHigh(i) - b.getLow(i) );
236 T d2 = fabs(getLow(i) - b.getLow(i));
237
238 if (d1 <= d2)
239 {
240 T d3 = fabs(getHigh(i) - b.getHigh(i));
241 dist.get(i) = (d3 < d1)?d3:d1;
242 }
243 else
244 {
245 T d3 = fabs(getHigh(i) - b.getHigh(i));
246 dist.get(i) = (d3 < d1)?d3:d2;
247 }
248 }
249 }
250
251 return dist.norm();
252 }
253
254 /*! \brief Get the coordinate of the bounding point
255 *
256 * \tparam b integer define (lower or upper interval)
257 * \param i i-coordinate of the point
258 * \return the i-coordinate of the bounding point
259 *
260 */
261
262 template<unsigned int b> T getBase(const unsigned int i) const
263 {
264 return boost::fusion::at_c<b>(data)[i];
265 }
266
267 /*! \brief Get the the box of dimensionality dim-1 (it eliminate the last dimension)
268 *
269 * \return Return the sub-box of dimension dim-1
270 *
271 */
272
273 Box<dim-1,T> getSubBox()
274 {
275 return Box<dim-1,T>(data);
276 }
277
278 /*! \brief Operator= between boxes
279 *
280 * Operator= between boxes of the same size
281 *
282 * \param box is the box that store the interval
283 * \return itself
284 *
285 */
286
287 __device__ __host__ Box<dim,T> & operator=(const Box<dim,T> & box)
288 {
289 for(size_t i = 0 ; i < dim ; i++)
290 {setLow(i,box.getLow(i));}
291
292 for(size_t i = 0 ; i < dim ; i++)
293 {setHigh(i,box.getHigh(i));}
294
295 // return itself
296 return *this;
297 }
298
299 public:
300
301 //! default constructor
302 __device__ __host__ Box()
303 {}
304
305 /*! \brief Constructor from two points
306 *
307 * \param p1 Low point, initialize as a list example {0.0,0.0,0.0}
308 * \param p2 High point, initialized as a list example {1.0,1.0,1.0}
309 *
310 */
311 Box(const Point<dim,T> & p1, const Point<dim,T> & p2)
312 {
313 setP1(p1);
314 setP2(p2);
315 }
316
317 /*! \brief Constructor from initializer list
318 *
319 * \param p1 Low point, initialize as a list example {0.0,0.0,0.0}
320 * \param p2 High point, initialized as a list example {1.0,1.0,1.0}
321 *
322 */
323 Box(std::initializer_list<T> p1, std::initializer_list<T> p2)
324 {
325 set(p1,p2);
326 }
327
328 /*! \brief Box constructor from a box
329 *
330 * \param high array indicating the coordinates of the low point
331 * \param low array indicating the coordinates of the high point
332 *
333 */
334 inline Box(T * low, T * high)
335 {
336 // copy all the data
337
338 for (int i = 0 ; i < dim ; i++)
339 {
340 // p1 is low p2 is high
341
342 boost::fusion::at_c<Box::p1>(data)[i] = low[i];
343 boost::fusion::at_c<Box::p2>(data)[i] = high[i];
344 }
345 }
346
347 /*! \brief Box constructor from a box
348 *
349 * \param box from which to construct
350 *
351 */
352 __device__ __host__ inline Box(const Box<dim,T> & box)
353 {
354 // we copy the data
355
356 for (size_t i = 0 ; i < dim ; i++)
357 {
358 boost::fusion::at_c<p1>(data)[i] = boost::fusion::at_c<p1>(box.data)[i];
359 boost::fusion::at_c<p2>(data)[i] = boost::fusion::at_c<p2>(box.data)[i];
360 }
361 }
362
363 /*! \brief Box constructor from vector::fusion
364 *
365 * \param box_data from which to construct
366 *
367 */
368 explicit inline Box(type box_data)
369 {
370 // we copy the data
371
372 for (size_t i = 0 ; i < dim ; i++)
373 {
374 boost::fusion::at_c<p1>(data)[i] = boost::fusion::at_c<p1>(box_data)[i];
375 boost::fusion::at_c<p2>(data)[i] = boost::fusion::at_c<p2>(box_data)[i];
376 }
377 }
378
379 /*! \brief Box constructor from an array reference
380 *
381 * \param box_data array from which to construct the box
382 *
383 */
384 inline Box(T (& box_data)[dim])
385 {
386 // we copy the data
387
388 for (size_t i = 0 ; i < dim ; i++)
389 {
390 boost::fusion::at_c<p1>(data)[i] = 0;
391 boost::fusion::at_c<p2>(data)[i] = box_data[i];
392 }
393 }
394
395 /*! \brief constructor from 2 grid_key_dx
396 *
397 * \param key1 start point
398 * \param key2 stop point
399 *
400 */
401 inline Box(const grid_key_dx<dim> & key1, const grid_key_dx<dim> & key2)
402 {
403 for (size_t i = 0 ; i < dim ; i++)
404 {
405 setLow(i,key1.get(i));
406 setHigh(i,key2.get(i));
407 }
408 }
409
410 /*! \brief Box constructor from vector::fusion of higher dimension
411 *
412 * \param box_data fusion vector from which to construct the vector
413 *
414 */
415 template<unsigned int dimS>
416 __device__ __host__ inline Box(boost::fusion::vector<T[dimS],T[dimS]> & box_data)
417 {
418 // we copy the data
419
420 for (size_t i = 0 ; i < dim ; i++)
421 {
422 boost::fusion::at_c<p1>(data)[i] = boost::fusion::at_c<p1>(box_data)[i];
423 boost::fusion::at_c<p2>(data)[i] = boost::fusion::at_c<p2>(box_data)[i];
424 }
425 }
426
427 /*! \brief Box constructor from encapsulated box
428 *
429 * \param b box from which to construct the vector (encapsulated)
430 *
431 */
432 template<typename Mem>
433 __device__ __host__
434 inline Box(const encapc<1,Box<dim,T>,Mem> & b)
435 {
436 // we copy the data
437
438 for (size_t i = 0 ; i < dim ; i++)
439 {
440 boost::fusion::at_c<p1>(data)[i] = b.template get<p1>()[i];
441 boost::fusion::at_c<p2>(data)[i] = b.template get<p2>()[i];
442 }
443 }
444
445 /*! \brief constructor from a Box of different type
446 *
447 * \param b box
448 *
449 */
450 template <typename S>
451 __device__ __host__
452 inline Box(const Box<dim,S> & b)
453 {
454 for (size_t d = 0 ; d < dim ; d++)
455 {this->setLow(d,b.getLow(d));}
456
457 for (size_t d = 0 ; d < dim ; d++)
458 {this->setHigh(d,b.getHigh(d));}
459 }
460
461 /*! \brief Divide component wise each box points with a point
462 *
463 * \param p point
464 *
465 * \return itself
466 *
467 */
468 inline Box<dim,T> & operator/=(const Point<dim,T> & p)
469 {
470 for (size_t i = 0 ; i < dim ; i++)
471 {
472 setLow(i, getLow(i)/p.get(i));
473 setHigh(i, getHigh(i)/p.get(i));
474 }
475 return *this;
476 }
477
478 /*! \brief Multiply component wise each box points with a point
479 *
480 * \param p point
481 *
482 * \return the result box
483 *
484 */
485 inline Box<dim,T> operator*(const Point<dim,T> & p)
486 {
487 Box<dim,T> ret;
488
489 for (size_t i = 0 ; i < dim ; i++)
490 {
491 ret.setLow(i, getLow(i)*p.get(i));
492 ret.setHigh(i, getHigh(i)*p.get(i));
493 }
494 return ret;
495 }
496
497 /*! \brief Constructor from initializer list
498 *
499 * Constructor from initializer list
500 *
501 * \param p1 Low point, initialize as a list example {0.0,0.0,0.0}
502 * \param p2 High point, initialized as a list example {1.0,1.0,1.0}
503 *
504 */
505
506 inline void set(std::initializer_list<T> p1, std::initializer_list<T> p2)
507 {
508 size_t i = 0;
509 for(T x : p1)
510 {
511 setLow(i,x);
512 i++;
513 if (i >= dim)
514 break;
515 }
516
517 i = 0;
518 for(T x : p2)
519 {
520 setHigh(i,x);
521 i++;
522 if (i >= dim)
523 break;
524 }
525 }
526
527 /*! \brief set the low interval of the box
528 *
529 * \param i dimension
530 * \param val value to set
531 *
532 */
533 __device__ __host__ inline void setLow(int i, T val)
534 {
535 boost::fusion::at_c<p1>(data)[i] = val;
536 }
537
538 /*! \brief set the high interval of the box
539 *
540 * \param i dimension
541 * \param val value to set
542 *
543 */
544 __device__ __host__ inline void setHigh(int i, T val)
545 {
546 boost::fusion::at_c<p2>(data)[i] = val;
547 }
548
549 /*! \brief get the i-coordinate of the low bound interval of the box
550 *
551 * \param i dimension
552 *
553 * \return i-coordinate
554 *
555 */
556 __device__ __host__ inline T getLow(int i) const
557 {
558 return boost::fusion::at_c<p1>(data)[i];
559 }
560
561 /*! \brief get the high interval of the box
562 *
563 * \param i dimension
564 * \return i coordinate of the high interval
565 *
566 */
567 __device__ __host__ inline T getHigh(int i) const
568 {
569 return boost::fusion::at_c<p2>(data)[i];
570 }
571
572 /*! \brief Set the point P1 of the box
573 *
574 * \param p1 point
575 *
576 */
577 inline void setP1(const grid_key_dx<dim> & p1)
578 {
579 for (size_t i = 0 ; i < dim ; i++)
580 setLow(i,p1.get(i));
581 }
582
583 /*! \brief Set the point P2 of the box
584 *
585 * \param p2 point
586 *
587 */
588 inline void setP2(const grid_key_dx<dim> & p2)
589 {
590 for (size_t i = 0 ; i < dim ; i++)
591 setHigh(i,p2.get(i));
592 }
593
594 /*! \brief Set the point P1 of the box
595 *
596 * \param p1 point
597 *
598 */
599 inline void setP1(const Point<dim,T> & p1)
600 {
601 for (size_t i = 0 ; i < dim ; i++)
602 setLow(i,p1.get(i));
603 }
604
605 /*! \brief Set the point P2 of the box
606 *
607 * \param p2 point
608 *
609 */
610 inline void setP2(const Point<dim,T> & p2)
611 {
612 for (size_t i = 0 ; i < dim ; i++)
613 setHigh(i,p2.get(i));
614 }
615
616 /*! \brief Get the box enclosing this Box
617 *
618 * basically return itself
619 *
620 * \return itself
621 *
622 */
623 Box<dim,T> & getBox()
624 {
625 return *this;
626 }
627
628 /*! \brief Get the box enclosing this Box
629 *
630 * basically return itself
631 *
632 * \return itself
633 *
634 */
635 const Box<dim,T> & getBox() const
636 {
637 return *this;
638 }
639
640 /*! \brief Get the internal boost::fusion::vector that store the data
641 *
642 * \return the internal boost::fusion::vector that store the data
643 *
644 */
645
646 type & getVector()
647 {
648 return data;
649 }
650
651 /*! \brief Get the point p1 as grid_key_dx
652 *
653 * \return the key
654 *
655 */
656 grid_key_dx<dim> getKP1() const
657 {
658 // grid key to return
659 grid_key_dx<dim> ret(boost::fusion::at_c<p1>(data));
660
661 return ret;
662 }
663
664 /*! \brief Get the point p12 as grid_key_dx
665 *
666 * \return the key
667 *
668 */
669 grid_key_dx<dim> getKP2() const
670 {
671 // grid key to return
672 grid_key_dx<dim> ret(boost::fusion::at_c<p2>(data));
673
674 return ret;
675 }
676
677 /*! \brief Get the point p1 as grid_key_dx
678 *
679 * \return the key
680 *
681 */
682 grid_key_dx<dim,int> getKP1int() const
683 {
684 // grid key to return
685 grid_key_dx<dim,int> ret(boost::fusion::at_c<p1>(data));
686
687 return ret;
688 }
689
690 /*! \brief Get the point p12 as grid_key_dx
691 *
692 * \return the key
693 *
694 */
695 grid_key_dx<dim,int> getKP2int() const
696 {
697 // grid key to return
698 grid_key_dx<dim,int> ret(boost::fusion::at_c<p2>(data));
699
700 return ret;
701 }
702
703 /*! \brief Get the point p1
704 *
705 * \return the point p1
706 *
707 */
708 inline Point<dim,T> getP1() const
709 {
710 // grid key to return
711 Point<dim,T> ret(boost::fusion::at_c<p1>(data));
712
713 return ret;
714 }
715
716 /*! \brief Get the point p2
717 *
718 * \return the point p2
719 *
720 */
721
722 inline Point<dim,T> getP2() const
723 {
724 // grid key to return
725 Point<dim,T> ret(boost::fusion::at_c<p2>(data));
726
727 return ret;
728 }
729
730 /*! \brief Translate the box
731 *
732 * \param p Point translation vector
733 *
734 * \return itself
735 *
736 */
737 inline Box<dim,T> & operator-=(const Point<dim,T> & p)
738 {
739 for (size_t i = 0 ; i < dim ; i++)
740 {
741 boost::fusion::at_c<p2>(data)[i] -= p.get(i);
742 boost::fusion::at_c<p1>(data)[i] -= p.get(i);
743 }
744
745 return *this;
746 }
747
748 /*! \brief Translate the box
749 *
750 * \param p Point translation vector
751 *
752 * \return itself
753 *
754 */
755 inline Box<dim,T> & operator+=(const Point<dim,T> & p)
756 {
757 for (size_t i = 0 ; i < dim ; i++)
758 {
759 boost::fusion::at_c<p2>(data)[i] += p.get(i);
760 boost::fusion::at_c<p1>(data)[i] += p.get(i);
761 }
762
763 return *this;
764 }
765
766 /*! \brief Translate the box
767 *
768 * \param p Point translation vector
769 *
770 * \return the translated box
771 *
772 */
773 inline Box<dim,T> operator+(const Point<dim,T> & p)
774 {
775 Box<dim,T> b;
776
777 for (size_t i = 0 ; i < dim ; i++)
778 {
779 b.setHigh(i,boost::fusion::at_c<p2>(data)[i] + p.get(i));
780 b.setLow(i,boost::fusion::at_c<p1>(data)[i] + p.get(i));
781 }
782
783 return b;
784 }
785
786 /*! \brief expand the box by a vector
787 *
788 * only P2 is expanded
789 *
790 * \param exp expand vector
791 *
792 */
793 inline void expand(T (& exp)[dim])
794 {
795 for (size_t i = 0 ; i < dim ; i++)
796 {
797 boost::fusion::at_c<p2>(data)[i] = boost::fusion::at_c<p2>(data)[i] + exp[i];
798 }
799 }
800
801 /*! \brief Enlarge the box with ghost margin
802 *
803 * \param gh spacing of the margin to enlarge
804 *
805 *
806 *\verbatim
807 ^ gh.p2[1]
808 |
809 |
810 +----+----+
811 | |
812 | |
813 gh.p1[0]<-----+ +----> gh.p2[0]
814 | |
815 | |
816 +----+----+
817 |
818 v gh.p1[1]
819
820 \endverbatim
821 *
822 */
823 void enlarge(const Box<dim,T> & gh)
824 {
825 typedef ::Box<dim,T> g;
826
827 for (size_t j = 0 ; j < dim ; j++)
828 {
829 this->setLow(j,this->template getBase<g::p1>(j) + gh.template getBase<g::p1>(j));
830 this->setHigh(j,this->template getBase<g::p2>(j) + gh.template getBase<g::p2>(j));
831 }
832 }
833
834 /*! \brief Enlarge the box with ghost margin keeping fix the point P1
835 *
836 * \param gh spacing of the margin to enlarge
837 *
838 *
839 *\verbatim
840 ^ gh.p2[1]
841 |
842 |
843 +----+----+
844 | |
845 | |
846 gh.p1[0]<-----+ +----> gh.p2[0]
847 | |
848 | |
849 +----+----+
850 |
851 v gh.p1[1]
852
853 \endverbatim
854 *
855 */
856 template<typename S> inline void enlarge_fix_P1(Box<dim,S> & gh)
857 {
858 typedef ::Box<dim,T> g;
859
860 for (size_t j = 0 ; j < dim ; j++)
861 {
862 this->setHigh(j,this->template getBase<g::p2>(j) + gh.template getBase<g::p2>(j) - gh.template getBase<g::p1>(j));
863 }
864 }
865
866 /*! \brief Invalidate the box
867 *
868 * Bring the state of this box in a way that isValid return false
869 *
870 */
871 void invalidate()
872 {
873 for (size_t j = 0 ; j < dim ; j++)
874 {
875 this->setLow(j,this->getHigh(j)+1);
876 }
877 }
878
879
880 /*! \brief Magnify the box
881 *
882 * For example 1.001 enlarge the box of 0.1% on each direction
883 *
884 * \warning P1 is mooved if not zero
885 *
886 * \param mg Magnification factor
887 *
888 */
889 void magnify(T mg)
890 {
891 typedef ::Box<dim,T> g;
892
893 for (size_t j = 0 ; j < dim ; j++)
894 {
895 this->setLow(j,mg * this->template getBase<g::p1>(j));
896 this->setHigh(j,mg * this->template getBase<g::p2>(j));
897 }
898 }
899
900 /*! \brief Magnify the box by a factor keeping fix the point P1
901 *
902 * For example 1.001 enlarge the box of 0.1% on each direction
903 *
904 * \param mg Magnification factor
905 *
906 */
907 inline void magnify_fix_P1(T mg)
908 {
909 typedef ::Box<dim,T> g;
910
911 for (size_t j = 0 ; j < dim ; j++)
912 {
913 this->setHigh(j,this->template getBase<g::p1>(j) + mg * (this->template getBase<g::p2>(j) - this->template getBase<g::p1>(j)));
914 }
915 }
916
917 /*! \brief Shrink moving p2 of sh quantity (on each direction)
918 *
919 * \param sh
920 *
921 */
922 inline void shrinkP2(T sh)
923 {
924 for (size_t j = 0 ; j < dim ; j++)
925 {
926 this->setHigh(j,this->getHigh(j) - sh);
927 }
928 }
929
930 /*! \brief Refine the box to enclose the given box and itself
931 *
932 * \param en Box to enclose
933 *
934 */
935 inline void enclose(const Box<dim,T> & en)
936 {
937 for (size_t j = 0 ; j < dim ; j++)
938 {
939 if (getLow(j) > en.getLow(j))
940 this->setLow(j,en.getLow(j));
941
942 if (getHigh(j) < en.getHigh(j))
943 this->setHigh(j,en.getHigh(j));
944 }
945 }
946
947 /*! \brief Refine the box to be contained in the given box and itself
948 *
949 * All the boxes are considered centered at p1, so it only count its relative size
950 *
951 * \param en Box to be contained
952 * \param reset_p1 if true set p1 to 0
953 *
954 */
955 inline void contained(const Box<dim,T> & en, const bool reset_p1 = true)
956 {
957 for (size_t j = 0 ; j < dim ; j++)
958 {
959 if (getHigh(j) > (en.getHigh(j) - en.getLow(j)))
960 setHigh(j,en.getHigh(j) - en.getLow(j));
961
962 if (reset_p1 == true)
963 setLow(j,0);
964 }
965 }
966
967 /*! \brief Set p1 and p2 to 0
968 *
969 */
970 inline void zero()
971 {
972 for (size_t j = 0 ; j < dim ; j++)
973 {
974 setHigh(j,0);
975 setLow(j,0);
976 }
977 }
978
979 /*! \brief Check if the box is contained
980 *
981 * \param b Box
982 *
983 * \return true if the box is contained
984 *
985 */
986 inline bool isContained(const Box<dim,T> & b) const
987 {
988 bool isc = true;
989
990 isc &= isInside(b.getP1());
991 isc &= isInside(b.getP2());
992
993 return isc;
994 }
995
996 /*! \brief Check if the point is inside the box
997 *
998 * \param p point to check
999 * \return true if the point is inside the space
1000 *
1001 */
1002
1003 inline
1004 __host__ __device__ bool isInside(const Point<dim,T> & p) const
1005 {
1006 // check if bound
1007
1008 for (size_t i = 0 ; i < dim ; i++)
1009 {
1010 // if outside the region return false
1011 if ( boost::fusion::at_c<Point<dim,T>::x>(p.data)[i] < boost::fusion::at_c<Box<dim,T>::p1>(this->data)[i]
1012 || boost::fusion::at_c<Point<dim,T>::x>(p.data)[i] > boost::fusion::at_c<Box<dim,T>::p2>(this->data)[i])
1013 {
1014 // Out of bound
1015
1016 return false;
1017 }
1018
1019 }
1020
1021 // In bound
1022
1023 return true;
1024 }
1025
1026 /*! \brief Check if the point is inside the region excluding the positive part
1027 *
1028 * In periodic boundary conditions the positive border is not included, but match the beginning
1029 *
1030 * \param p point to check
1031 * \return true if the point is inside the space
1032 *
1033 */
1034 __device__ __host__ inline bool isInsideNP(const Point<dim,T> & p) const
1035 {
1036 // check if bound
1037
1038 for (size_t i = 0 ; i < dim ; i++)
1039 {
1040 // if outside the region return false
1041 if ( boost::fusion::at_c<Point<dim,T>::x>(p.data)[i] < boost::fusion::at_c<Box<dim,T>::p1>(this->data)[i]
1042 || boost::fusion::at_c<Point<dim,T>::x>(p.data)[i] >= boost::fusion::at_c<Box<dim,T>::p2>(this->data)[i])
1043 {
1044 // Out of bound
1045
1046
1047
1048 return false;
1049 }
1050
1051 }
1052
1053 // In bound
1054
1055 return true;
1056 }
1057
1058 /*! \brief Check if the point is inside the region excluding the positive part
1059 *
1060 * In periodic boundary conditions the positive border is not included, but match the beginning
1061 *
1062 * \param p point to check
1063 * \return true if the point is inside the space
1064 *
1065 */
1066 template<typename bc_type>
1067 __device__ __host__ inline bool isInsideNP_with_border(const Point<dim,T> & p, const Box<dim,T> & border, const bc_type (& bc)[dim]) const
1068 {
1069 // check if bound
1070
1071 for (size_t i = 0 ; i < dim ; i++)
1072 {
1073 // if outside the region return false
1074 if ( p.get(i) < this->getLow(i)
1075 || (bc[i] == NON_PERIODIC && ((this->getHigh(i) != border.getHigh(i) && p.get(i) >= this->getHigh(i)) || (this->getHigh(i) == border.getHigh(i) && p.get(i) > this->getHigh(i)) ) )
1076 || (bc[i] == PERIODIC && p.get(i) >= this->getHigh(i)))
1077 {
1078 // Out of bound
1079 return false;
1080 }
1081
1082 }
1083
1084 // In bound
1085
1086 return true;
1087 }
1088
1089 /*! \brief Check if the point is inside the region excluding the borders
1090 *
1091 * \param p point to check
1092 * \return true if the point is inside the space
1093 *
1094 */
1095 inline bool isInsideNB(const Point<dim,T> & p) const
1096 {
1097 // check if bound
1098
1099 for (size_t i = 0 ; i < dim ; i++)
1100 {
1101 // if outside the region return false
1102 if ( boost::fusion::at_c<Point<dim,T>::x>(p.data)[i] <= boost::fusion::at_c<Box<dim,T>::p1>(this->data)[i]
1103 || boost::fusion::at_c<Point<dim,T>::x>(p.data)[i] >= boost::fusion::at_c<Box<dim,T>::p2>(this->data)[i])
1104 {
1105 // Out of bound
1106
1107 return false;
1108 }
1109
1110 }
1111
1112 // In bound
1113
1114 return true;
1115 }
1116
1117 /*! \brief Check if the point is inside the region (Border included)
1118 *
1119 * \param p point to check
1120 * \return true if the point is inside the space
1121 *
1122 */
1123
1124 inline bool isInside(const T (&p)[dim]) const
1125 {
1126 // check if bound
1127
1128 for (size_t i = 0 ; i < dim ; i++)
1129 {
1130 // if outside the region return false
1131 if ( p[i] < boost::fusion::at_c<Box<dim,T>::p1>(this->data)[i]
1132 || p[i] > boost::fusion::at_c<Box<dim,T>::p2>(this->data)[i])
1133 {
1134 // Out of bound
1135
1136 return false;
1137 }
1138
1139 }
1140
1141 // In bound
1142
1143 return true;
1144 }
1145
1146 /*! \brief Check if the point is inside the region (Border included)
1147 *
1148 * \param k key to check
1149 * \return true if the point is inside the space
1150 *
1151 */
1152 template<typename KeyType>
1153 __device__ __host__ inline bool isInsideKey(const KeyType & k) const
1154 {
1155 // check if bound
1156
1157 for (size_t i = 0 ; i < dim ; i++)
1158 {
1159 // if outside the region return false
1160 if ( k.get(i) < boost::fusion::at_c<Box<dim,T>::p1>(this->data)[i]
1161 || k.get(i) > boost::fusion::at_c<Box<dim,T>::p2>(this->data)[i])
1162 {
1163 // Out of bound
1164
1165 return false;
1166 }
1167
1168 }
1169
1170 // In bound
1171
1172 return true;
1173 }
1174
1175 /*! \brief Check if the Box is a valid box P2 >= P1
1176 *
1177 * \return true if it is valid
1178 *
1179 */
1180 inline bool isValid() const
1181 {
1182 for (size_t i = 0 ; i < dim ; i++)
1183 {
1184 if (getLow(i) > getHigh(i))
1185 return false;
1186 }
1187
1188 return true;
1189 }
1190
1191 /*! \brief Check if the Box is a valid box P2 > P1
1192 *
1193 * \return true if it is valid
1194 *
1195 */
1196 inline bool isValidN() const
1197 {
1198 for (size_t i = 0 ; i < dim ; i++)
1199 {
1200 if (getLow(i) >= getHigh(i))
1201 return false;
1202 }
1203
1204 return true;
1205 }
1206
1207 /*! \brief Apply the ceil operation to the point P1
1208 *
1209 *
1210 */
1211 inline void floorP1()
1212 {
1213 for (size_t i = 0 ; i < dim ; i++)
1214 {
1215 setLow(i,std::floor(getLow(i)));
1216 }
1217 }
1218
1219 /*! \brief Apply the ceil operation to the point P2
1220 *
1221 *
1222 */
1223 inline void floorP2()
1224 {
1225 for (size_t i = 0 ; i < dim ; i++)
1226 {
1227 setHigh(i,std::floor(getHigh(i)));
1228 }
1229 }
1230
1231 /*! \brief Apply the ceil operation to the point P1
1232 *
1233 *
1234 */
1235 inline void ceilP1()
1236 {
1237 for (size_t i = 0 ; i < dim ; i++)
1238 {
1239 setLow(i,std::ceil(getLow(i)));
1240 }
1241 }
1242
1243 /*! \brief Apply the ceil operation to the point P2
1244 *
1245 *
1246 */
1247 inline void ceilP2()
1248 {
1249 for (size_t i = 0 ; i < dim ; i++)
1250 {
1251 setHigh(i,std::ceil(getHigh(i)));
1252 }
1253 }
1254
1255 /*! \brief Shrink the point P2 by one vector
1256 *
1257 * \param p vector
1258 *
1259 */
1260 inline void shrinkP2(const Point<dim,T> & p)
1261 {
1262 for (size_t i = 0 ; i < dim ; i++)
1263 {
1264 setHigh(i,getHigh(i) - p.get(i));
1265 }
1266 }
1267
1268 /*! \brief exchange the data of two boxes
1269 *
1270 * \param b box to switch
1271 *
1272 */
1273 void swap(Box<dim,T> & b)
1274 {
1275 for (size_t i = 0 ; i < dim ; i++)
1276 {
1277 T tmp_l = getLow(i);
1278 T tmp_h = getHigh(i);
1279
1280 setLow(i,b.getLow(i));
1281 setHigh(i,b.getHigh(i));
1282
1283 b.setLow(i,tmp_l);
1284 b.setHigh(i,tmp_h);
1285 }
1286 }
1287
1288 /*! \brief Check if the point is inside the region
1289 *
1290 * \param p point to check
1291 * \return true if the point is inside the space
1292 *
1293 */
1294
1295 template <typename Mem>
1296 inline bool isInside(const encapc<1,Point<dim,T>,Mem> & p)
1297 {
1298 // check if bound
1299
1300 for (size_t i = 0 ; i < dim ; i++)
1301 {
1302 // if outside the region return false
1303 if ( p.template get<Point<dim,T>::x>()[i] < boost::fusion::at_c<Box<dim,T>::p1>(this->data)[i]
1304 || p.template get<Point<dim,T>::x>()[i] > boost::fusion::at_c<Box<dim,T>::p2>(this->data)[i])
1305 {
1306 // Out of bound
1307
1308 return false;
1309 }
1310
1311 }
1312
1313 // In bound
1314
1315 return true;
1316 }
1317
1318 /*! \brief Get the worst extension
1319 *
1320 * \return the worst extension
1321 */
1322 inline T getRcut() const
1323 {
1324 T r_cut = 0;
1325 for (size_t i = 0 ; i < dim ; i++)
1326 {r_cut = std::max(r_cut,getHigh(i));}
1327
1328 return r_cut;
1329 }
1330
1331 /*! \brief Get the volume of the box
1332 *
1333 * \return the box volume
1334 *
1335 */
1336 inline T getVolume() const
1337 {
1338 T vol = 1.0;
1339
1340 for (size_t i = 0 ; i < dim ; i++)
1341 vol *= (getHigh(i) - getLow(i));
1342
1343 return vol;
1344 }
1345
1346 /*! \brief Get the volume spanned by the Box P1 and P2 interpreted as grid key
1347 *
1348 * \return The volume
1349 *
1350 */
1351 inline T getVolumeKey() const
1352 {
1353 T vol = 1.0;
1354
1355 for (size_t i = 0 ; i < dim ; i++)
1356 vol *= (getHigh(i) - getLow(i) + 1.0);
1357
1358 return vol;
1359 }
1360
1361 /*! \brief Get the volume spanned by the Box as grid_key_dx_iterator_sub
1362 *
1363 * \warning Careful it is not the simple volume calculation there is a +1 on each dimension, consider the case of a subgrid iterator with
1364 * P1 = {5,7} ; P2 = {5,7}, the sub-grid iterator has one point {5,7}, that mean Volume=1, so
1365 * the volume formula is (5 - 5 + 1) * (7 - 7 + 1)
1366 *
1367 * \param p1 point p1
1368 * \param p2 point p2
1369 *
1370 * \return The volume
1371 *
1372 */
1373 inline static T getVolumeKey(const T (&p1)[dim], const T(&p2)[dim])
1374 {
1375 T vol = 1.0;
1376
1377 for (size_t i = 0 ; i < dim ; i++)
1378 vol *= (p2[i] - p1[i] + 1.0);
1379
1380 return vol;
1381 }
1382
1383 //! This structure has no internal pointers
1384 static bool noPointers()
1385 {
1386 return true;
1387 }
1388
1389 /*! \brief Return the middle point of the box
1390 *
1391 * \return the middle point of the box
1392 *
1393 */
1394 inline Point<dim,T> middle() const
1395 {
1396 Point<dim,T> p;
1397
1398 for (size_t i = 0 ; i < dim ; i++)
1399 p.get(i) = (getLow(i) + getHigh(i))/2;
1400
1401 return p;
1402 }
1403
1404 /*! \brief Produce a string from the object
1405 *
1406 * \return string
1407 *
1408 */
1409 std::string toString() const
1410 {
1411 std::stringstream str;
1412
1413 for (size_t i = 0 ; i < dim ; i++)
1414 str << "x[" << i << "]=" << getLow(i) << " ";
1415
1416 str << " | ";
1417
1418 for (size_t i = 0 ; i < dim ; i++)
1419 str << "x[" << i << "]=" << getHigh(i) << " ";
1420
1421 return str.str();
1422 }
1423
1424 /*! \brief Compare two boxes
1425 *
1426 * \param b
1427 *
1428 * \return true if the boxes are equal
1429 *
1430 */
1431 bool operator==(const Box<dim,T> & b) const
1432 {
1433 for (size_t i = 0 ; i < dim ; i++)
1434 {
1435 if (getLow(i) != b.getLow(i))
1436 return false;
1437
1438 if (getHigh(i) != b.getHigh(i))
1439 return false;
1440 }
1441
1442 return true;
1443 }
1444
1445
1446 /*! \brief Compare two boxes
1447 *
1448 * \param b
1449 *
1450 * \return true if the boxes are equal
1451 *
1452 */
1453 bool operator!=(const Box<dim,T> & b) const
1454 {
1455 return ! this->operator==(b);
1456 }
1457};
1458
1459template<typename T, typename Sfinae = void>
1460struct is_Box: std::false_type {};
1461
1462
1463/*! \brief Check if a type T is an aggregate
1464 *
1465 * return true if T is an aggregate
1466 *
1467 */
1468template<typename T>
1469struct is_Box<T, typename Void< typename T::yes_is_box>::type> : std::true_type
1470{};
1471
1472#endif
1473