CellDecomposer.hpp

/*
 * CellDecomposer.hpp
 *
 *  Created on: Apr 1, 2015
 *      Author: Pietro Incardona
 */

#ifndef CELLDECOMPOSER_HPP_
#define CELLDECOMPOSER_HPP_

#include "Space/SpaceBox.hpp"
#include "Space/Matrix.hpp"
#include "util/copy_compare/meta_compare.hpp"

// Shift transformation
template<unsigned int dim, typename T>
class shift
{
    Point<dim,T> sh;

public:

    shift(const Matrix<dim,T> & t, const Point<dim,T> & s)
    :sh(s)
    {
    }

    inline T transform(const T(&s)[dim], const size_t i) const
    {
        return s[i] - sh.get(i);
    }

    inline T transform(const Point<dim,T> & s, const size_t i) const
    {
        return s.get(i) - sh.get(i);
    }

    template<typename Mem> inline T transform(const encapc<1,Point<dim,T>,Mem> & s, const size_t i) const
    {
        return s.get(i) - sh.get(i);
    }

    inline void setTransform(const Matrix<dim,T> & mat, const Point<dim,T> & orig)
    {
        for (size_t i = 0 ; i < dim ; i++)
            sh.get(i) = orig.get(i);
    }
};

// No transformation
template<unsigned int dim, typename T>
class no_transform
{
public:

    no_transform(const Matrix<dim,T> & t, const Point<dim,T> & s)
    {
    }

    inline T transform(const T(&s)[dim], const size_t i) const
    {
        return s[i];
    }

    inline T transform(const Point<dim,T> & s, const size_t i) const
    {
        return s.get(i);
    }

    template<typename Mem> inline T transform(const encapc<1,Point<dim,T>,Mem> & s, const size_t i) const
    {
        return s.template get<Point<dim,T>::x>()[i];
    }

    inline void setTransform(const Matrix<dim,T> & mat, const Point<dim,T> & orig)
    {

    }

    inline bool operator==(const no_transform<dim,T> & nt)
    {
        return true;
    }

    inline bool operator!=(const no_transform<dim,T> & nt)
    {
        return false;
    }
};


template<unsigned int dim,typename T, typename transform = no_transform<dim,T>>
class CellDecomposer_sm
{
    // Point in the middle of a cell
    //
    // \verbatim
    //
    //         C (0.1,0.1)
    // +-----+
    // |     |
    // |  .P |
    // |     |
    // +-----+
    //
    // \endverbatim
    //
    //    C is the cell and P is the point inside the middle of the cell
    //    for example if the cell is (0.1,0.1) P is (0.05,0.05)
    //
    //
    Point<dim,T> p_middle;

    // Point transformation before get the Cell object (useful for example to shift the cell list)
    transform t;

    inline size_t ConvertToID(const T (&x)[dim] ,size_t s) const
    {
        size_t id = (size_t)(t.transform(x,s) / box_unit.getHigh(s)) + off[s];
        id = (id >= (gr_cell.size(s) + off[0]))?(gr_cell.size(s)-1):id;
        return id;
    }

    inline size_t ConvertToID(const Point<dim,T> & x ,size_t s) const
    {
        size_t id = (size_t)(t.transform(x,s) / box_unit.getHigh(s)) + off[s];
        id = (id >= (gr_cell.size(s) + off[0]))?(gr_cell.size(s)-1):id;
        return id;
    }

    template <typename Mem> inline size_t ConvertToID_(const encapc<1,Point<dim,T>,Mem> & x ,size_t s) const
    {
        size_t id = (size_t)(t.transform(x,s) / box_unit.getHigh(s)) + off[s];
        id = (id >= (gr_cell.size(s) + off[0]))?(gr_cell.size(s)-1):id;
        return id;
    }

protected:

    // Total number of cell
    size_t tot_n_cell;

    // Domain of the cell list
    SpaceBox<dim,T> box;

    // Unit box of the Cell list
    SpaceBox<dim,T> box_unit;

    // Grid structure of the Cell list
    grid_sm<dim,void> gr_cell;

    // cell padding on each dimension
    size_t off[dim];


    void Initialize(const size_t pad, const size_t (& div)[dim])
    {
#ifdef DEBUG
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (div[i] == 0)
            {
                std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " the number of cells on each dimension must be different from zero\n";
            }
        }
#endif

        // created a padded div
        size_t div_p[dim];

        for (size_t i = 0 ; i < dim ; i++)
            div_p[i] = div[i] + 2*pad;

        gr_cell.setDimensions(div_p);

        tot_n_cell = 1;

        // Total number of cells and calculate the unit cell size

        for (size_t i = 0 ; i < dim ; i++)
        {
            tot_n_cell *= gr_cell.size(i);

            // Cell are padded by
            box_unit.setHigh(i,box.getHigh(i) / (gr_cell.size(i)- 2*pad) );
        }

        for (size_t i = 0; i < dim ; i++)
            off[i] = pad;

        // Initialize p_middle

        p_middle = box_unit.getP2();
        p_middle = p_middle / 2;
    }

public:

    const grid_sm<dim,void> & getGrid() const
    {
#ifdef DEBUG
        if (tot_n_cell == 0)
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " using an uninitialized CellDecomposer";
#endif

        return gr_cell;
    }

    inline grid_key_dx<dim> getCellGrid(const T (& pos)[dim]) const
    {
#ifdef SE_CLASS1
        if (tot_n_cell == 0)
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " using an uninitialized CellDecomposer";
#endif

        grid_key_dx<dim> key;
        key.set_d(0,ConvertToID(pos,0));

        for (size_t s = 1 ; s < dim ; s++)
        {
#ifdef SE_CLASS1
            if ((size_t)(t.transform(pos,s) / box_unit.getHigh(s)) + off[s] < 0)
                std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point is not inside the cell space\n";
#endif
            key.set_d(s,ConvertToID(pos,s));

        }

        return key;
    }

    grid_key_dx<dim> getCellGrid(const Point<dim,T> pos) const
    {
#ifdef SE_CLASS1
        if (tot_n_cell == 0)
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " using an uninitialized CellDecomposer\n";
#endif

        grid_key_dx<dim> key;
        key.set_d(0,ConvertToID(pos,0));

        for (size_t s = 1 ; s < dim ; s++)
        {
#ifdef SE_CLASS1
            if ((size_t)(t.transform(pos,s) / box_unit.getHigh(s)) + off[s] < 0)
                std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point is not inside the cell space\n";
#endif
            key.set_d(s,ConvertToID(pos,s));
        }

        return key;
    }

    size_t getCell(const T (& pos)[dim]) const
    {
#ifdef SE_CLASS1
        if (tot_n_cell == 0)
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " using an uninitialized CellDecomposer";

        if (t.transform(pos,0) < box.getLow(0) || t.transform(pos,0) > box.getHigh(0))
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point " << toPointString(pos) << " is not inside the cell space";
#endif

        size_t cell_id = ConvertToID(pos,0);

        for (size_t s = 1 ; s < dim ; s++)
        {
#ifdef SE_CLASS1
            if (t.transform(pos,s) < box.getLow(s) || t.transform(pos,s) > box.getHigh(s))
                std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point " << toPointString(pos) << " is not inside the cell space";
#endif
            cell_id += gr_cell.size_s(s-1) * ConvertToID(pos,s);
        }

        return cell_id;
    }

    size_t getCell(const Point<dim,T> & pos) const
    {
#ifdef SE_CLASS1
        if (tot_n_cell == 0)
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " using an uninitialized CellDecomposer";

        if (t.transform(pos,0) < box.getLow(0) || t.transform(pos,0) > box.getHigh(0))
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point " << pos.toPointString() << " is not inside the cell space";
#endif

        size_t cell_id = ConvertToID(pos,0);

        for (size_t s = 1 ; s < dim ; s++)
        {
#ifdef DEBUG
            if (t.transform(pos,s) < box.getLow(s) || t.transform(pos,s) > box.getHigh(s))
                std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point " << pos.toPointString() << " is not inside the cell space";
#endif
            cell_id += gr_cell.size_s(s-1) * ConvertToID(pos,s);
        }

        return cell_id;
    }

    template<typename Mem> size_t getCell(const encapc<1,Point<dim,T>,Mem> & pos) const
    {

#ifdef SE_CLASS1
        if (tot_n_cell == 0)
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " using an uninitialized CellDecomposer";

        if (t.transform(pos,0) < box.getLow(0) || t.transform(pos,0) > box.getHigh(0))
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point " << toPointString(pos) << " is not inside the cell space";
#endif

        size_t cell_id = ConvertToID_(pos,0);

        for (size_t s = 1 ; s < dim ; s++)
        {
#ifdef SE_CLASS1
            if (t.transform(pos,s) < box.getLow(s) || t.transform(pos,s) > box.getHigh(s))
                std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point " << toPointString(pos) << " is not inside the cell space";
#endif
            cell_id += gr_cell.size_s(s-1) * ConvertToID_(pos,s);
        }

        return cell_id;
    }

    Box<dim,size_t> getGridPoints(const Box<dim,T> & s_box) const
    {
        // Box with inside grid
        Box<dim,size_t> bx;

        // Point p2
        Point<dim,T> p2 = s_box.getP2();
        p2 = p2 - p_middle;

        // Point p1
        Point<dim,T> p1 = s_box.getP1();
        p1 = p1 + p_middle;

        bx.setP2(getCellGrid(p2));
        bx.setP1(getCellGrid(p1));

        return bx;
    }

    void setDimensions(const Box<dim,T> & box, const size_t (&div)[dim], const size_t pad)
    {
        this->box = box;
        this->gr_cell.setDimensions(div);
        Initialize(pad,div);
    }

    void setDimensions(const Box<dim,T> & box, const size_t (&div)[dim], const Matrix<dim,T> & mat, const Point<dim,T> & orig, const size_t pad)
    {
        t.setTransform(mat,orig);
        this->box = box;
        this->gr_cell.setDimensions(div);
        Initialize(pad,div);
    }

    CellDecomposer_sm(const SpaceBox<dim,T> & box, const size_t (&div)[dim], Matrix<dim,T> & mat, Point<dim,T> & orig, const size_t pad)
    :t(Matrix<dim,T>::identity(),Point<dim,T>::zero()),box(box),gr_cell()
    {
        Initialize(pad);
    }

    CellDecomposer_sm(const SpaceBox<dim,T> & box, const size_t (&div)[dim], Point<dim,T> & orig, const size_t pad)
    :t(Matrix<dim,T>::identity(),orig),box(box),gr_cell(div)
    {
        Initialize(pad,div);
    }

    CellDecomposer_sm(const SpaceBox<dim,T> & box, const size_t (&div)[dim], const size_t pad)
    :t(Matrix<dim,T>::identity(),Point<dim,T>::zero_p()),box(box),gr_cell()
    {
        Initialize(pad,div);
    }


    CellDecomposer_sm()
    :t(Matrix<dim,T>::identity(),Point<dim,T>::zero_p()),tot_n_cell(0)
    {

    }

    const Box<dim,T> & getCellBox() const
    {
        return box_unit;
    }

/*  Box<dim,T> BoundaryFixation()
    {

    }*/

    Box<dim,size_t> convertDomainSpaceIntoGridUnits(const Box<dim,T> & b_d) const
    {
        Box<dim,size_t> g_box;
        Box<dim,T> b = b_d;

        // Convert b into grid units
        b /= getCellBox().getP2();

        // Considering that we are interested in a box decomposition of the space
        // where basically all the point are uniquely assigned we include the positive
        // countour and exclude the negative one. So ceilP1 do the job for P1 while ceilP2 - 1
        // do the job for P2

        b.ceilP1();
        b.ceilP2();

        g_box = b;

        // on the other hand if we are at the positive (with non periodic boundary condition)
        // we have to include also the positive border

        Point<dim,size_t> p_move;

        for (size_t i = 0 ; i < dim ; i++)
        {
            // we are at the positive border (We are assuming that there are not rounding error, check boundary Fixation)
            if (b_d.getHigh(i) == box.getHigh(i))
            {
                p_move.get(i) = 0;
                g_box.setHigh(i,gr_cell.size(i));
            }
            else
                p_move.get(i) = 1;
        }

        g_box.shrinkP2(p_move);

        return g_box;
    }

    bool operator==(const CellDecomposer_sm<dim,T,transform> & cd)
    {
        if (meta_compare<Point<dim,T>>::meta_compare_f(p_middle,cd.p_middle) == false)
            return false;

        if (t != cd.t)
            return false;

        if (tot_n_cell != cd.tot_n_cell)
            return false;

        if (box != cd.box)
            return false;

        if (box_unit != cd.box_unit)
            return false;

        if (gr_cell != cd.gr_cell)
            return false;

        for (size_t i = 0 ; i < dim ; i++)
        {
            if (off[i] != cd.off[i])
                return false;
        }

        return true;
    }

    bool operator!=(const CellDecomposer_sm<dim,T,transform> & cd)
    {
        return ! this->operator==(cd);
    }
};


#endif /* CELLDECOMPOSER_HPP_ */