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 "Grid/Encap.hpp"
#include <sstream>

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;

    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;

    // Get method
    template <int i> inline auto get() -> decltype(boost::fusion::at_c<i>(data))
    {
        return boost::fusion::at_c<i>(data);
    }

    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> bool Intersect(const encapc<1,Box<dim,T>,Mem> & e_b, Box<dim,T> & b_out) const
    {
        return Intersect(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();
    }

    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);
    }

    Box<dim,T> & operator=(const Box<dim,T> & box)
    {
        for(size_t i = 0 ; i < dim ; i++)
        {setLow(i,box.getLow(i));}

        for(size_t i = 0 ; i < dim ; i++)
        {setHigh(i,box.getHigh(i));}

        // return itself
        return *this;
    }

    public:

    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 * high, T * low)
    {
        // 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];
        }
    }

    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];
        }
    }

    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> 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> 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> inline Box(const Box<dim,S> & b)
    {
        for (size_t d = 0 ; d < dim ; d++)
        {this->setLow(d,b.getLow(d));}

        for (size_t d = 0 ; d < dim ; d++)
        {this->setHigh(d,b.getHigh(d));}
    }

    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;
        }
    }

    inline void setLow(int i, T val)
    {
        boost::fusion::at_c<p1>(data)[i] = val;
    }

    inline void setHigh(int i, T val)
    {
        boost::fusion::at_c<p2>(data)[i] = val;
    }

    inline T getLow(int i) const
    {
        return boost::fusion::at_c<p1>(data)[i];
    }

    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;
    }

    type & getVector()
    {
        return data;
    }

    /* \brief Get the point p1 as grid_key_dx
     *
     * \return the key
     *
     */
    grid_key_dx<dim> getKP1() const
    {
        // grid key to return
        grid_key_dx<dim> ret(boost::fusion::at_c<p1>(data));

        return ret;
    }

    /* \brief Get the point p12 as grid_key_dx
     *
     * \return the key
     *
     */
    grid_key_dx<dim> getKP2() const
    {
        // grid key to return
        grid_key_dx<dim> ret(boost::fusion::at_c<p2>(data));

        return ret;
    }

    /* \brief Get the point p1
     *
     * \return the point p1
     *
     */
    inline Point<dim,T> getP1() const
    {
        // grid key to return
        Point<dim,T> ret(boost::fusion::at_c<p1>(data));

        return ret;
    }

    /* \brief Get the point p2
     *
     * \return the point p2
     *
     */

    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;
    }

    /* \brief expand expand the box by a vector
     *
     * \param vector
     *
     */
    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++)
        {
            this->setLow(j,this->template getBase<g::p1>(j) + gh.template getBase<g::p1>(j));
            this->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++)
        {
            this->setHigh(j,this->template getBase<g::p2>(j) + gh.template getBase<g::p2>(j) - gh.template getBase<g::p1>(j));
        }
    }


    void magnify(T mg)
    {
        typedef ::Box<dim,T> g;

        for (size_t j = 0 ; j < dim ; j++)
        {
            this->setLow(j,mg * this->template getBase<g::p1>(j));
            this->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++)
        {
            this->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++)
        {
            this->setHigh(j,this->getHigh(j) - sh);
        }
    }

    inline void enclose(Box<dim,T> & en)
    {
        for (size_t j = 0 ; j < dim ; j++)
        {
            if (getLow(j) > en.getLow(j))
                this->setLow(j,en.getLow(j));

            if (getHigh(j) < en.getHigh(j))
                this->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);
        }
    }

    void zero()
    {
        for (size_t j = 0 ; j < dim ; j++)
        {
            setHigh(j,0);
            setLow(j,0);
        }
    }

    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;
    }

    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;
    }

    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;
    }


    bool isValid() const
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (getLow(i) > getHigh(i))
                return false;
        }

        return true;
    }

    void ceilP1()
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            setLow(i,std::ceil(getLow(i)));
        }
    }

    void ceilP2()
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            setHigh(i,std::ceil(getHigh(i)));
        }
    }

    void shrinkP2(const Point<dim,T> & p)
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            setHigh(i,getHigh(i) - p.get(i));
        }
    }

    template <typename Mem> 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 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(T (&p1)[dim], 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);
    }
};

#endif