Box.hpp

 
#ifndef BOX_HPP_
#define BOX_HPP_
 
 
#include "Space/Shape/Sphere.hpp"
#include <boost/fusion/sequence/intrinsic/at_c.hpp>
#include "Grid/grid_key.hpp"
#include "memory_ly/Encap.hpp"
#include <sstream>
 
#define PERIODIC 1
#define NON_PERIODIC 0
 
enum Base
{
    UP,
    DOWN
};
 
template<unsigned int dim , typename T>
class Box
{
public:
 
    typedef boost::fusion::vector<T[dim],T[dim]> type;
    typedef T btype;
 
    typedef int yes_is_box;
 
    type data;
 
    static const unsigned int p1 = 0;
    static const unsigned int p2 = 1;
    static const unsigned int max_prop = 2;
 
    static const unsigned int dims = dim;
 
    __device__ __host__ bool Intersect(const Box<dim,T> & b, Box<dim,T> & b_out) const
    {
        // check if p1 of b is smaller than
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (getLow(i) <= b.getLow(i))
                b_out.setLow(i,b.getLow(i));
            else if (getLow(i) <= b.getHigh(i))
                b_out.setLow(i,getLow(i));
            else
                return false;
 
            if (getHigh(i) >= b.getHigh(i))
                b_out.setHigh(i,b.getHigh(i));
            else if (getHigh(i) >= b.getLow(i))
                b_out.setHigh(i,getHigh(i));
            else
                return false;
        }
        return true;
    }
 
 
 
    template<typename Mem>
    __device__ __host__
    bool Intersect(const encapc<1,Box<dim,T>,Mem> & e_b, Box<dim,T> & b_out) const
    {
        return Intersect(Box<dim,T>(e_b),b_out);
    }
 
    template <typename distance> bool Intersect(Sphere<dim,T> & sphere)
    {
        // distance functor
        distance dist;
 
        // Get the nearest point of the box from the center of the sphere
        typename distance::ResultType distance_r = 0;
 
        for (size_t i = 0 ; i < dim ; i++)
        {
 
            // if the center of the sphere on dimension i is not in the i box interval
            // do not accumulate, otherwise accumulate from the nearest point on that
            // dimension
            if (boost::fusion::at_c<p1>(data)[i] < sphere.center(i))
            {
                // accumulate the distance from p1
                distance_r += dist.accum_dist(sphere.center(i),boost::fusion::at_c<p1>(data)[i],i);
            }
            else if ( boost::fusion::at_c<p2>(data)[i] <= sphere.center(i))
            {
                // accumulate the distance from p2
                distance_r += dist.accum_dist(sphere.center(i),boost::fusion::at_c<p2>(data)[i],i);
            }
        }
 
        // return if there is intersection
        return distance_r < sphere.radius();
    }
 
    /* \brief return the minimum distance between two boxes
     *
     * Consider two box in space like in the figure below. Given two points P and Q in the two boxes 1 and 2.
     * The minimum distance, is the distance that satify min( dist(P,Q) ). There P and Q is the standard euclidean distance
     *
     * \verbatim b box 2
     *
     * \verbatim
 
 
 
                       _____ p2
                      |    |
                      |    |   Box 2
                      |p1__|
                     /
                    /
                   /    <----------- min distance
                  /
            _____/
           |  p2|
           |    |  Box 1
           |____|
          p1
 
     \endverbatim
     *
     *
     */
    T min_distance(Box<dim,T> & b) const
    {
        Point<dim,T> dist;
 
        for (unsigned int i = 0 ; i < dim ; i++)
        {
            if (getHigh(i) <= b.getLow(i))
            {dist.get(i) = b.getLow(i) - getHigh(i);}
            else if (b.getHigh(i) <= getLow(i))
            {dist.get(i) = getLow(i) - b.getHigh(i);}
            else
            {
                T d1 = fabs(getHigh(i) - b.getLow(i) );
                T d2 = fabs(getLow(i) - b.getLow(i));
 
                if (d1 <= d2)
                {
                    T d3 = fabs(getHigh(i) - b.getHigh(i));
                    dist.get(i) = (d3 < d1)?d3:d1;
                }
                else
                {
                    T d3 = fabs(getHigh(i) - b.getHigh(i));
                    dist.get(i) = (d3 < d1)?d3:d2;
                }
            }
        }
 
        return dist.norm();
    }
 
    template<unsigned int b> T getBase(const unsigned int i) const
    {
        return boost::fusion::at_c<b>(data)[i];
    }
 
    Box<dim-1,T> getSubBox()
    {
        return Box<dim-1,T>(data);
    }
 
    __device__ __host__ Box<dim,T> & operator=(const Box<dim,T> & box)
    {
        for(size_t i = 0 ; i < dim ; i++)
        {
            setLow(i,box.getLow(i));
            setHigh(i,box.getHigh(i));
        }
 
        // return itself
        return *this;
    }
 
    public:
 
    __device__ __host__ Box()
    {}
 
    Box(const Point<dim,T> & p1, const Point<dim,T> & p2)
    {
        setP1(p1);
        setP2(p2);
    }
 
    Box(std::initializer_list<T> p1, std::initializer_list<T> p2)
    {
        set(p1,p2);
    }
 
    inline Box(T * low, T * high)
    {
        // copy all the data
 
        for (int i = 0 ; i < dim ; i++)
        {
            // p1 is low p2 is high
 
            boost::fusion::at_c<Box::p1>(data)[i] = low[i];
            boost::fusion::at_c<Box::p2>(data)[i] = high[i];
        }
    }
 
    __device__ __host__ inline Box(const Box<dim,T> & box)
    {
        // we copy the data
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            boost::fusion::at_c<p1>(data)[i] = boost::fusion::at_c<p1>(box.data)[i];
            boost::fusion::at_c<p2>(data)[i] = boost::fusion::at_c<p2>(box.data)[i];
        }
    }
 
    explicit inline Box(type box_data)
    {
        // we copy the data
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            boost::fusion::at_c<p1>(data)[i] = boost::fusion::at_c<p1>(box_data)[i];
            boost::fusion::at_c<p2>(data)[i] = boost::fusion::at_c<p2>(box_data)[i];
        }
    }
 
    inline Box(T (& box_data)[dim])
    {
        // we copy the data
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            boost::fusion::at_c<p1>(data)[i] = 0;
            boost::fusion::at_c<p2>(data)[i] = box_data[i];
        }
    }
 
    inline Box(const grid_key_dx<dim> & key1, const grid_key_dx<dim> & key2)
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            setLow(i,key1.get(i));
            setHigh(i,key2.get(i));
        }
    }
 
    template<unsigned int dimS>
    __device__ __host__ inline Box(boost::fusion::vector<T[dimS],T[dimS]> & box_data)
    {
        // we copy the data
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            boost::fusion::at_c<p1>(data)[i] = boost::fusion::at_c<p1>(box_data)[i];
            boost::fusion::at_c<p2>(data)[i] = boost::fusion::at_c<p2>(box_data)[i];
        }
    }
 
    template<typename Mem>
    __device__ __host__
    inline Box(const encapc<1,Box<dim,T>,Mem> & b)
    {
        // we copy the data
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            boost::fusion::at_c<p1>(data)[i] = b.template get<p1>()[i];
            boost::fusion::at_c<p2>(data)[i] = b.template get<p2>()[i];
        }
    }
 
    template <typename S>
    __device__ __host__
    inline Box(const Box<dim,S> & b)
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            setLow(i,b.getLow(i));
            setHigh(i,b.getHigh(i));
        }
    }
 
    inline Box<dim,T> & operator/=(const Point<dim,T> & p)
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            setLow(i, getLow(i)/p.get(i));
            setHigh(i, getHigh(i)/p.get(i));
        }
        return *this;
    }
 
    inline Box<dim,T> operator*(const Point<dim,T> & p)
    {
        Box<dim,T> ret;
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            ret.setLow(i, getLow(i)*p.get(i));
            ret.setHigh(i, getHigh(i)*p.get(i));
        }
        return ret;
    }
 
    inline void set(std::initializer_list<T> p1, std::initializer_list<T> p2)
    {
        size_t i = 0;
        for(T x : p1)
        {
            setLow(i,x);
            i++;
            if (i >= dim)
                break;
        }
 
        i = 0;
        for(T x : p2)
        {
            setHigh(i,x);
            i++;
            if (i >= dim)
                break;
        }
    }
 
    __device__ __host__  inline void setLow(int i, T val)
    {
        boost::fusion::at_c<p1>(data)[i] = val;
    }
 
    __device__ __host__  inline void setHigh(int i, T val)
    {
        boost::fusion::at_c<p2>(data)[i] = val;
    }
 
    __device__ __host__  inline T getLow(int i) const
    {
        return boost::fusion::at_c<p1>(data)[i];
    }
 
    __device__ __host__ inline T getHigh(int i) const
    {
        return boost::fusion::at_c<p2>(data)[i];
    }
 
    inline void setP1(const grid_key_dx<dim> & p1)
    {
        for (size_t i = 0 ; i < dim ; i++)
            setLow(i,p1.get(i));
    }
 
    inline void setP2(const grid_key_dx<dim> & p2)
    {
        for (size_t i = 0 ; i < dim ; i++)
            setHigh(i,p2.get(i));
    }
 
    inline void setP1(const Point<dim,T> & p1)
    {
        for (size_t i = 0 ; i < dim ; i++)
            setLow(i,p1.get(i));
    }
 
    inline void setP2(const Point<dim,T> & p2)
    {
        for (size_t i = 0 ; i < dim ; i++)
            setHigh(i,p2.get(i));
    }
 
    Box<dim,T> & getBox()
    {
        return *this;
    }
 
    const Box<dim,T> & getBox() const
    {
        return *this;
    }
 
    type & getVector()
    {
        return data;
    }
 
    grid_key_dx<dim> getKP1() const
    {
        // grid key to return
        grid_key_dx<dim> ret(boost::fusion::at_c<p1>(data));
 
        return ret;
    }
 
    grid_key_dx<dim> getKP2() const
    {
        // grid key to return
        grid_key_dx<dim> ret(boost::fusion::at_c<p2>(data));
 
        return ret;
    }
 
    grid_key_dx<dim,int> getKP1int() const
    {
        // grid key to return
        grid_key_dx<dim,int> ret(boost::fusion::at_c<p1>(data));
 
        return ret;
    }
 
    grid_key_dx<dim,int> getKP2int() const
    {
        // grid key to return
        grid_key_dx<dim,int> ret(boost::fusion::at_c<p2>(data));
 
        return ret;
    }
 
    inline Point<dim,T> getP1() const
    {
        // grid key to return
        Point<dim,T> ret(boost::fusion::at_c<p1>(data));
 
        return ret;
    }
 
    inline Point<dim,T> getP2() const
    {
        // grid key to return
        Point<dim,T> ret(boost::fusion::at_c<p2>(data));
 
        return ret;
    }
 
    inline Box<dim,T> & operator-=(const Point<dim,T> & p)
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            boost::fusion::at_c<p2>(data)[i] -= p.get(i);
            boost::fusion::at_c<p1>(data)[i] -= p.get(i);
        }
 
        return *this;
    }
 
    inline Box<dim,T> & operator+=(const Point<dim,T> & p)
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            boost::fusion::at_c<p2>(data)[i] += p.get(i);
            boost::fusion::at_c<p1>(data)[i] += p.get(i);
        }
 
        return *this;
    }
 
    inline Box<dim,T> operator+(const Point<dim,T> & p) const
    {
        Box<dim,T> b;
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            b.setHigh(i,boost::fusion::at_c<p2>(data)[i] + p.get(i));
            b.setLow(i,boost::fusion::at_c<p1>(data)[i] + p.get(i));
        }
 
        return b;
    }
 
    inline void expand(T (& exp)[dim])
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            boost::fusion::at_c<p2>(data)[i] = boost::fusion::at_c<p2>(data)[i] + exp[i];
        }
    }
 
    void enlarge(const Box<dim,T> & gh)
    {
        typedef ::Box<dim,T> g;
 
        for (size_t j = 0 ; j < dim ; j++)
        {
            setLow(j,this->template getBase<g::p1>(j) + gh.template getBase<g::p1>(j));
            setHigh(j,this->template getBase<g::p2>(j) + gh.template getBase<g::p2>(j));
        }
    }
 
    template<typename S> inline void enlarge_fix_P1(Box<dim,S> & gh)
    {
        typedef ::Box<dim,T> g;
 
        for (size_t j = 0 ; j < dim ; j++)
        {
            setHigh(j,this->template getBase<g::p2>(j) + gh.template getBase<g::p2>(j) - gh.template getBase<g::p1>(j));
        }
    }
 
    void invalidate()
    {
        for (size_t j = 0 ; j < dim ; j++)
        {
            setLow(j,this->getHigh(j)+1);
        }
    }
 
 
    void magnify(T mg)
    {
        typedef ::Box<dim,T> g;
 
        for (size_t j = 0 ; j < dim ; j++)
        {
            setLow(j,mg * this->template getBase<g::p1>(j));
            setHigh(j,mg * this->template getBase<g::p2>(j));
        }
    }
 
    inline void magnify_fix_P1(T mg)
    {
        typedef ::Box<dim,T> g;
 
        for (size_t j = 0 ; j < dim ; j++)
        {
            setHigh(j,this->template getBase<g::p1>(j) + mg * (this->template getBase<g::p2>(j) - this->template getBase<g::p1>(j)));
        }
    }
 
    inline void shrinkP2(T sh)
    {
        for (size_t j = 0 ; j < dim ; j++)
        {
            setHigh(j,this->getHigh(j) - sh);
        }
    }
 
    inline void enlargeP2(T sh)
    {
        for (size_t j = 0 ; j < dim ; j++)
        {
            setHigh(j,this->getHigh(j) + sh);
        }
    }
 
    inline void enclose(const Box<dim,T> & en)
    {
        for (size_t j = 0 ; j < dim ; j++)
        {
            if (getLow(j) > en.getLow(j))
                setLow(j,en.getLow(j));
 
            if (getHigh(j) < en.getHigh(j))
                setHigh(j,en.getHigh(j));
        }
    }
 
    inline void contained(const Box<dim,T> & en, const bool reset_p1 = true)
    {
        for (size_t j = 0 ; j < dim ; j++)
        {
            if (getHigh(j) > (en.getHigh(j) - en.getLow(j)))
                setHigh(j,en.getHigh(j) - en.getLow(j));
 
            if (reset_p1 == true)
                setLow(j,0);
        }
    }
 
    inline void zero()
    {
        for (size_t j = 0 ; j < dim ; j++)
        {
            setHigh(j,0);
            setLow(j,0);
        }
    }
 
    inline bool isContained(const Box<dim,T> & b) const
    {
        bool isc = true;
 
        isc &= isInside(b.getP1());
        isc &= isInside(b.getP2());
 
        return isc;
    }
 
    inline
    __host__ __device__ bool isInside(const Point<dim,T> & p) const
    {
        // check if bound
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            // if outside the region return false
            if (   boost::fusion::at_c<Point<dim,T>::x>(p.data)[i] < boost::fusion::at_c<Box<dim,T>::p1>(this->data)[i]
                || boost::fusion::at_c<Point<dim,T>::x>(p.data)[i] > boost::fusion::at_c<Box<dim,T>::p2>(this->data)[i])
            {
                // Out of bound
 
                return false;
            }
        }
 
        // In bound
 
        return true;
    }
 
    __device__ __host__ inline bool isInsideNP(const Point<dim,T> & p) const
    {
        // check if bound
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            // if outside the region return false
            if (   boost::fusion::at_c<Point<dim,T>::x>(p.data)[i] < boost::fusion::at_c<Box<dim,T>::p1>(this->data)[i]
                || boost::fusion::at_c<Point<dim,T>::x>(p.data)[i] >= boost::fusion::at_c<Box<dim,T>::p2>(this->data)[i])
            {
                // Out of bound
 
 
 
                return false;
            }
        }
 
        // In bound
 
        return true;
    }
 
    template<typename bc_type>
    __device__ __host__ inline bool isInsideNP_with_border(const Point<dim,T> & p, const Box<dim,T> & border, const bc_type (& bc)[dim]) const
    {
        // check if bound
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            // if outside the region return false
            if (   p.get(i) < this->getLow(i)
                || (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)) ) )
                || (bc[i] == PERIODIC && p.get(i) >= this->getHigh(i)))
            {
                // Out of bound
                return false;
            }
        }
 
        // In bound
 
        return true;
    }
 
    inline bool isInsideNB(const Point<dim,T> & p) const
    {
        // check if bound
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            // if outside the region return false
            if (   boost::fusion::at_c<Point<dim,T>::x>(p.data)[i] <= boost::fusion::at_c<Box<dim,T>::p1>(this->data)[i]
                || boost::fusion::at_c<Point<dim,T>::x>(p.data)[i] >= boost::fusion::at_c<Box<dim,T>::p2>(this->data)[i])
            {
                // Out of bound
 
                return false;
            }
        }
 
        // In bound
 
        return true;
    }
 
    inline bool isInside(const T (&p)[dim]) const
    {
        // check if bound
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            // if outside the region return false
            if (   p[i] < boost::fusion::at_c<Box<dim,T>::p1>(this->data)[i]
                || p[i] > boost::fusion::at_c<Box<dim,T>::p2>(this->data)[i])
            {
                // Out of bound
 
                return false;
            }
        }
 
        // In bound
 
        return true;
    }
 
    template<typename KeyType>
    __device__ __host__ inline bool isInsideKey(const KeyType & k) const
    {
        // check if bound
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            // if outside the region return false
            if (   k.get(i) < boost::fusion::at_c<Box<dim,T>::p1>(this->data)[i]
                || k.get(i) > boost::fusion::at_c<Box<dim,T>::p2>(this->data)[i])
            {
                // Out of bound
 
                return false;
            }
        }
 
        // In bound
 
        return true;
    }
 
    inline bool isValid() const
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (getLow(i) > getHigh(i))
                return false;
        }
 
        return true;
    }
 
    inline bool isValidN() const
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (getLow(i) >= getHigh(i))
                return false;
        }
 
        return true;
    }
 
    inline void floorP1()
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            setLow(i,std::floor(getLow(i)));
        }
    }
 
    inline void floorP2()
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            setHigh(i,std::floor(getHigh(i)));
        }
    }
 
    inline void ceilP1()
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            setLow(i,std::ceil(getLow(i)));
        }
    }
 
    inline void ceilP2()
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            setHigh(i,std::ceil(getHigh(i)));
        }
    }
 
    inline void shrinkP2(const Point<dim,T> & p)
    {
        for (size_t i = 0 ; i < dim ; i++)
            setHigh(i,getHigh(i) - p.get(i));
    }
 
    void swap(Box<dim,T> & b)
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            T tmp_l = getLow(i);
            T tmp_h = getHigh(i);
 
            setLow(i,b.getLow(i));
            setHigh(i,b.getHigh(i));
 
            b.setLow(i,tmp_l);
            b.setHigh(i,tmp_h);
        }
    }
 
    template <typename Mem>
    inline bool isInside(const encapc<1,Point<dim,T>,Mem> & p)
    {
        // check if bound
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            // if outside the region return false
            if (   p.template get<Point<dim,T>::x>()[i] < boost::fusion::at_c<Box<dim,T>::p1>(this->data)[i]
                || p.template get<Point<dim,T>::x>()[i] > boost::fusion::at_c<Box<dim,T>::p2>(this->data)[i])
            {
                // Out of bound
 
                return false;
            }
        }
 
        // In bound
 
        return true;
    }
 
    inline T getRcut() const
    {
        T r_cut = 0;
        for (size_t i = 0 ; i < dim ; i++)
            r_cut = std::max(r_cut, getHigh(i));
 
        return r_cut;
    }
 
    inline T getVolume() const
    {
        T vol = 1.0;
 
        for (size_t i = 0 ; i < dim ; i++)
            vol *= (getHigh(i) - getLow(i));
 
        return vol;
    }
 
    inline T getVolumeKey() const
    {
        T vol = 1.0;
 
        for (size_t i = 0 ; i < dim ; i++)
            vol *= (getHigh(i) - getLow(i) + 1.0);
 
        return vol;
    }
 
    inline static T getVolumeKey(const T (&p1)[dim], const T(&p2)[dim])
    {
        T vol = 1.0;
 
        for (size_t i = 0 ; i < dim ; i++)
            vol *= (p2[i] - p1[i] + 1.0);
 
        return vol;
    }
 
    static bool noPointers()
    {
        return true;
    }
 
    inline Point<dim,T> middle() const
    {
        Point<dim,T> p;
 
        for (size_t i = 0 ; i < dim ; i++)
            p.get(i) = (getLow(i) + getHigh(i))/2;
 
        return p;
    }
 
    std::string toString() const
    {
        std::stringstream str;
 
        for (size_t i = 0 ; i < dim ; i++)
            str << "x[" << i << "]=" << getLow(i) << " ";
 
        str << "   |  ";
 
        for (size_t i = 0 ; i < dim ; i++)
            str << "x[" << i << "]=" << getHigh(i) << " ";
 
        return str.str();
    }
 
    bool operator==(const Box<dim,T> & b) const
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (getLow(i) != b.getLow(i))
                return false;
 
            if (getHigh(i) != b.getHigh(i))
                return false;
        }
 
        return true;
    }
 
    bool operator!=(const Box<dim,T> & b) const
    {
        return ! this->operator==(b);
    }
 
    void rescale(T (& sp)[dim])
    {
        for (size_t d = 0 ; d < dim ; d++)
            setHigh(d,getLow(d) + (getHigh(d) -getLow(d)) * sp[d]);
    }
 
    template<typename S> void rescale(S (& sp)[dim])
    {
        for (size_t d = 0 ; d < dim ; d++)
            setHigh(d,getLow(d) + (getHigh(d) -getLow(d)) * sp[d]);
    }
 
    void mul(T (& sp)[dim])
    {
        for (size_t d = 0 ; d < dim ; d++)
        {
            setLow(d,getLow(d) * sp[d]);
            setHigh(d,getHigh(d) * sp[d]);
        }
    }
 
    template<typename S> void mul(S (& sp)[dim])
    {
        for (size_t d = 0 ; d < dim ; d++)
        {
            setLow(d,getLow(d) * sp[d]);
            setHigh(d,getHigh(d) * sp[d]);
        }
    }
 
    Point<dim,T> rnd()
    {
        Point<dim,T> p;
 
        for (size_t i = 0 ; i < dim ; i++)
            p.get(i) = ((T)rand())/(T)RAND_MAX * (getHigh(i) - getLow(i)) + getLow(i);
 
        return p;
    }
};
 
#endif