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

#define CELL_DECOMPOSER 8001lu




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

    Matrix<dim,T> mat;

public:

    __device__ __host__ shift(const Matrix<dim,T> & t, const Point<dim,T> & s)
    :sh(s)
    {
        mat.identity();
    }

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

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

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

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

    inline const Point<dim,T> & getOrig() const
    {
        return sh;
    }

    inline const Matrix<dim,T> & getMat() const
    {
        return mat;
    }

    inline bool operator==(const shift<dim,T> & s)
    {
        return sh == s.sh;
    }

    inline bool operator!=(const shift<dim,T> & s)
    {
        return !this->operator==(s);
    }
};

template<unsigned int dim, typename T>
class shift_only
{
    Point<dim,T> sh;

public:

    shift_only()
    {
        sh.zero();
    }

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

    __device__ __host__  inline T transform(const T * s, const int i) const
    {
        return s[i] - sh.get(i);
    }

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

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

    template<typename Mem> __device__ __host__   inline T transform(const encapc<1,Point<dim,T>,Mem> & s, const int i) const
    {
        return s.template get<0>()[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);
    }

    __device__ __host__ inline const Point<dim,T> & getOrig() const
    {
        return sh;
    }


    inline bool operator==(const shift<dim,T> & s)
    {
        return sh == s.sh;
    }

    inline bool operator!=(const shift<dim,T> & s)
    {
        return !this->operator==(s);
    }
};

template<unsigned int dim, typename T>
class no_transform
{
    Point<dim,T> orig;

    Matrix<dim,T> mat;

public:

    __device__ __host__ no_transform(const Matrix<dim,T> & t, const Point<dim,T> & s)
    {
        orig.zero();
        mat.identity();
    }

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

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

    template<typename Mem> __device__ __host__ 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;
    }

    inline const Point<dim,T> & getOrig() const
    {
        return orig;
    }

    inline const Matrix<dim,T> & getMat() const
    {
        return mat;
    }
};

template<unsigned int dim, typename T>
class no_transform_only
{

public:

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

    __device__ __host__ inline T transform(const T * s, const int i) const
    {
        return s[i];
    }

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

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

    template<typename Mem> __device__ __host__ inline T transform(const encapc<1,Point<dim,T>,Mem> & s, const int 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 = openfpm::math::size_t_floor(t.transform(x,s) / box_unit.getHigh(s)) + off[s];
        id = (id >= gr_cell.size(s))?(gr_cell.size(s)-1-cell_shift.get(s)):id-cell_shift.get(s);
        return id;
    }

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

    inline size_t ConvertToID_me(const Point<dim,T> & x ,size_t s, size_t sc = 0) const
    {
        T cc = t.transform(x,s) / box_unit.getHigh(s) - 0.015625;
        size_t id = openfpm::math::size_t_floor(cc) + off[s];
        id = (id >= gr_cell.size(s))?(gr_cell.size(s)-1-cell_shift.get(s)):id-cell_shift.get(s);
        return id;
    }

    inline size_t ConvertToID_pe(const Point<dim,T> & x ,size_t s, size_t sc = 0) const
    {
        T cc = t.transform(x,s) / box_unit.getHigh(s) + 0.015625;
        size_t id = openfpm::math::size_t_floor(cc) + off[s];
        id = (id >= gr_cell.size(s))?(gr_cell.size(s)-1-cell_shift.get(s)):id-cell_shift.get(s);
        return id;
    }

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

    inline size_t ConvertToID_ns(const Point<dim,T> & x ,size_t s, size_t sc = 0) const
    {
        size_t id = openfpm::math::size_t_floor(t.transform(x,s) / box_unit.getHigh(s)) + off[s];
        id = (id >= gr_cell.size(s))?(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, size_t sc = 0) const
    {
        size_t id = (size_t)(t.transform(x,s) / box_unit.getHigh(s)) + off[s];
        id = (id >= gr_cell.size(s))?(gr_cell.size(s)-1-cell_shift.get(s)):id-cell_shift.get(s);
        return id;
    }

    template<typename Ele> inline void check_and_print_error(const Ele & pos ,size_t s) const
    {
#ifdef SE_CLASS1
        if (tot_n_cell == 0)
        {
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " using an uninitialized CellDecomposer";
            ACTION_ON_ERROR(CELL_DECOMPOSER);
        }

        if (pos[s] < box.getLow(s) - off[s]*box_unit.getP2()[s] || pos[s] > box.getHigh(s) + off[s]*box_unit.getP2()[s])
        {
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point " << Point<dim,T>(pos).toString() << " is not inside the cell space";
            ACTION_ON_ERROR(CELL_DECOMPOSER);
        }
#endif
    }


    template<typename Ele> inline size_t getCellDom_impl(const Ele & pos) const
    {
        check_and_print_error(pos,0);

        size_t cell_id = ConvertToID_ns(pos,0);

        cell_id = (cell_id == gr_cell.size(0) - off[0])?gr_cell.size(0) - off[0] - 1:cell_id;
        cell_id = (cell_id == off[0]-1)?off[0]:cell_id;

        cell_id -= cell_shift.get(0);

        for (size_t s = 1 ; s < dim ; s++)
        {
            /* coverity[dead_error_begin] */
            check_and_print_error(pos,s);

            size_t cell_idt = ConvertToID_ns(pos,s);

            cell_idt = (cell_idt == gr_cell.size(s) - off[s])?gr_cell.size(s) - off[s] - 1:cell_idt;
            cell_idt = (cell_idt == off[s]-1)?off[s]:cell_idt;

            cell_idt -= cell_shift.get(s);

            cell_id += gr_cell2.size_s(s-1) * cell_idt;
        }

        return cell_id;
    }

    template<typename Ele> inline size_t getCellPad_impl(const Ele & pos) const
    {
        check_and_print_error(pos,0);

        size_t cell_id = ConvertToID_ns(pos,0);
        cell_id = (cell_id == off[0])?off[0]-1:cell_id;
        cell_id = (cell_id == gr_cell.size(0) - off[0] - 1)?gr_cell.size(0) - off[0]:cell_id;

        cell_id -= cell_shift.get(0);

        for (size_t s = 1 ; s < dim ; s++)
        {
            check_and_print_error(pos,s);

            size_t cell_idt = ConvertToID_ns(pos,s);
            cell_idt = (cell_idt == off[s])?off[s]-1:cell_idt;
            cell_idt = (cell_idt == gr_cell.size(s) - off[s] - 1)?gr_cell.size(s) - off[s]:cell_idt;

            cell_idt -= cell_shift.get(s);

            cell_id += gr_cell2.size_s(s-1) * cell_idt;
        }

        return cell_id;
    }


protected:

    size_t tot_n_cell;

    SpaceBox<dim,T> box;

    SpaceBox<dim,T> box_unit;

    grid_sm<dim,void> gr_cell;

    grid_sm<dim,void> gr_cell2;

    Box<dim,T> box_gr_cell2;

    size_t off[dim];

    Point<dim,long int> cell_shift;

    // temporal buffer
    mutable size_t div_wp[dim];

    void Initialize(const size_t pad, const size_t (& div)[dim])
    {
#ifdef SE_CLASS1

        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";
                ACTION_ON_ERROR(CELL_DECOMPOSER)
            }
        }

#endif

        // set div_wp to zero
        for (size_t i = 0 ; i < dim ; i++)
        {div_wp[i] = 0;}

        // 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);
        gr_cell2.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) - box.getLow(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 transform & getTransform()
    {
        return t;
    }

    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";
            ACTION_ON_ERROR(CELL_DECOMPOSER);
        }
#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";
                ACTION_ON_ERROR(CELL_DECOMPOSER);
            }
#endif
            key.set_d(s,ConvertToID(pos,s));

        }

        return key;
    }

    inline size_t getCellLin(grid_key_dx<dim> && k) const
    {
#ifdef SE_CLASS1
        if (tot_n_cell == 0)
        {
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " using an uninitialized CellDecomposer";
            ACTION_ON_ERROR(CELL_DECOMPOSER);
        }

        if (gr_cell2.size(0) < k.get(0) + off[0] - cell_shift.get(0))
        {
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " cell coordinate 0 = " << k.get(0) + off[0] - cell_shift.get(0) << " bigger than cell space " << gr_cell2.size(0) << std::endl;
            ACTION_ON_ERROR(CELL_DECOMPOSER);
        }
#endif

        size_t cell_id = k.get(0) + off[0] - cell_shift.get(0);

        for (size_t s = 1 ; s < dim ; s++)
        {
#ifdef SE_CLASS1
            if (gr_cell2.size(s) < k.get(s) + off[s] - cell_shift.get(s))
            {
                std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " cell coordinate 0 = " << k.get(0) + off[0] - cell_shift.get(0) << " bigger than cell space " << gr_cell2.size(s) << std::endl;
                ACTION_ON_ERROR(CELL_DECOMPOSER);
            }
#endif
            cell_id += gr_cell2.size_s(s-1) * (k.get(s) + off[s] -cell_shift.get(s));
        }

        return cell_id;
    }

    inline size_t getCellLin(const grid_key_dx<dim> & k) const
    {
#ifdef SE_CLASS1
        if (tot_n_cell == 0)
        {
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " using an uninitialized CellDecomposer";
            ACTION_ON_ERROR(CELL_DECOMPOSER);
        }

        if (gr_cell2.size(0) < k.get(0) + off[0] - cell_shift.get(0))
        {
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " cell coordinate 0 = " << k.get(0) + off[0] - cell_shift.get(0) << " bigger than cell space " << gr_cell2.size(0) << std::endl;
            ACTION_ON_ERROR(CELL_DECOMPOSER);
        }
#endif

        size_t cell_id = k.get(0) + off[0] -cell_shift.get(0);

        for (size_t s = 1 ; s < dim ; s++)
        {
#ifdef SE_CLASS1
            if (gr_cell2.size(s) < k.get(s) + off[s] - cell_shift.get(s))
            {
                std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " cell coordinate 0 = " << k.get(0) + off[0] - cell_shift.get(0) << " bigger than cell space " << gr_cell2.size(s) << std::endl;
                ACTION_ON_ERROR(CELL_DECOMPOSER);
            }
#endif
            cell_id += gr_cell2.size_s(s-1) * (k.get(s) + off[s] -cell_shift.get(s));
        }

        return cell_id;
    }

    inline grid_key_dx<dim> getCellGrid_me(const Point<dim,T> & pos) const
    {
#ifdef SE_CLASS1
        if (tot_n_cell == 0)
        {
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " using an uninitialized CellDecomposer" << std::endl;
            ACTION_ON_ERROR(CELL_DECOMPOSER);
        }
#endif

        grid_key_dx<dim> key;
        key.set_d(0,ConvertToID_me(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" << std::endl;
                ACTION_ON_ERROR(CELL_DECOMPOSER);
            }
#endif
            /* coverity[dead_error_line] */
            key.set_d(s,ConvertToID_me(pos,s));
        }

        return key;
    }

    inline grid_key_dx<dim> getCellGrid_pe(const Point<dim,T> & pos) const
    {
#ifdef SE_CLASS1
        if (tot_n_cell == 0)
        {
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " using an uninitialized CellDecomposer" << std::endl;
            ACTION_ON_ERROR(CELL_DECOMPOSER);
        }
#endif

        grid_key_dx<dim> key;
        key.set_d(0,ConvertToID_pe(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" << std::endl;
                ACTION_ON_ERROR(CELL_DECOMPOSER);
            }
#endif
            /* coverity[dead_error_line] */
            key.set_d(s,ConvertToID_pe(pos,s));
        }

        return key;
    }

    inline 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" << std::endl;
            ACTION_ON_ERROR(CELL_DECOMPOSER);
        }
#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" << std::endl;
                ACTION_ON_ERROR(CELL_DECOMPOSER);
            }
#endif
            /* coverity[dead_error_line] */
            key.set_d(s,ConvertToID(pos,s));
        }

        return key;
    }

    inline size_t getCellDom(const Point<dim,T> & pos) const
    {
        return getCellDom_impl<Point<dim,T>>(pos);
    }


    inline size_t getCellDom(const T (& pos)[dim]) const
    {
        return getCellDom_impl<T[dim]>(pos);
    }

    inline size_t getCellPad(const Point<dim,T> & pos) const
    {
        return getCellPad_impl<Point<dim,T>>(pos);
    }

    inline size_t getCellPad(const T (& pos)[dim]) const
    {
        return getCellPad_impl<T[dim]>(pos);
    }

    inline 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";
            ACTION_ON_ERROR(CELL_DECOMPOSER);
        }

        if (pos[0] < box.getLow(0) - off[0]*box_unit.getP2()[0] || pos[0] > box.getHigh(0) + off[0]*box_unit.getP2()[0])
        {
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point " << toPointString(pos) << " is not inside the cell space";
            ACTION_ON_ERROR(CELL_DECOMPOSER);
        }
#endif

        size_t cell_id = ConvertToID(pos,0);

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

        return cell_id;
    }

    inline 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";
            ACTION_ON_ERROR(CELL_DECOMPOSER);
        }

        if (pos.get(0) < box.getLow(0) - off[0]*box_unit.getP2()[0] || pos.get(0) > box.getHigh(0) + off[0]*box_unit.getP2()[0])
        {
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point " << pos.toPointString() << " is not inside the cell space";
            ACTION_ON_ERROR(CELL_DECOMPOSER);
        }
#endif

        size_t cell_id = ConvertToID(pos,0);

        for (size_t s = 1 ; s < dim ; s++)
        {
#ifdef SE_CLASS1
            if (pos.get(s) < box.getLow(s) - off[s]*box_unit.getP2()[s] || pos.get(s) > box.getHigh(s) + off[s]*box_unit.getP2()[s])
            {
                std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point " << pos.toPointString() << " is not inside the cell space";
                ACTION_ON_ERROR(CELL_DECOMPOSER);
            }
#endif
            /* coverity[dead_error_line] */
            cell_id += gr_cell2.size_s(s-1) * ConvertToID(pos,s);
        }

        return cell_id;
    }

    template<typename Mem> inline 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";
            ACTION_ON_ERROR(CELL_DECOMPOSER);
        }

        if (pos.template get<0>()[0] < box.getLow(0) - off[0]*box_unit.getP2()[0] || pos.template get<0>()[0] > box.getHigh(0) + off[0]*box_unit.getP2()[0])
        {
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point " << toPointString(pos) << " is not inside the cell space " << box.toString() << std::endl;
            ACTION_ON_ERROR(CELL_DECOMPOSER);
        }
#endif

        size_t cell_id = ConvertToID_(pos,0);

        for (size_t s = 1 ; s < dim ; s++)
        {
#ifdef SE_CLASS1

            if (pos.template get<0>()[s] < box.getLow(s) - off[s]*box_unit.getP2()[s] || pos.template get<0>()[s] > box.getHigh(s) + off[s]*box_unit.getP2()[s])
            {
                std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " point " << toPointString(pos) << " is not inside the cell space";
                ACTION_ON_ERROR(CELL_DECOMPOSER);
            }
#endif
            cell_id += gr_cell2.size_s(s-1) * ConvertToID_(pos,s);
        }

        return cell_id;
    }

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

    inline void setDimensions(const Box<dim,T> & box, const size_t (&div)[dim], const size_t (&div2)[dim], const size_t pad, Point<dim,long int> cell_shift)
    {
        Matrix<dim,T> mat;
        mat.identity();
        t.setTransform(mat,box.getP1());
        this->box = box;

        Initialize(pad,div);

        size_t cells[dim];

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

        gr_cell2.setDimensions(cells);

        for (size_t i = 0 ; i < dim ; i++)
            this->cell_shift.get(i) = cell_shift.get(i) - off[i];
    }

    inline void setDimensions(const Box<dim,T> & box, const size_t (&div)[dim], const size_t pad)
    {
        Matrix<dim,T> mat;
        mat.identity();
        t.setTransform(mat,box.getP1());
        this->box = box;
        Initialize(pad,div);
        this->cell_shift = 0;
    }

    inline void setDimensions(const Box<dim,T> & box, const size_t (&div)[dim], const size_t (&div2)[dim], const Matrix<dim,T> & mat, const size_t pad, Point<dim,long int> cell_shift)
    {
        t.setTransform(mat,box.getP1());
        this->box = box;

        Initialize(pad,div);

        // The nested cell is big div2 + 2*off

        size_t div_with_off[dim];

        for(size_t i = 0 ; i < dim ; i++)
            div_with_off[i] = div2[i] + 2*off[i];

        gr_cell2.setDimensions(div_with_off);

        for (size_t i = 0 ; i < dim ; i++)
            this->cell_shift.get(i) = cell_shift.get(i) - off[i];
    }

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


    inline void setDimensions(const CellDecomposer_sm<dim,T,transform> & cd, const Box<dim,size_t> & cell_extension)
    {
        this->cell_shift = 0;

        // Get the space transformation

        t.setTransform(cd.getMat(),cd.getOrig());

        // The domain is equivalent to the old one
        this->box = cd.box;

        // The padding must be calculated

        size_t pad = 0;

        for (size_t i = 0 ; i < dim ; i++)
        {
            if (pad < cell_extension.getLow(i))
                pad = cell_extension.getLow(i);
            else if (pad > cell_extension.getHigh(i))
                pad = cell_extension.getHigh(i);
        }

        // We have to give the old division

        size_t sz_div[dim];

        for (size_t i = 0 ; i < dim ; i++)
            sz_div[i] = cd.gr_cell.size(i) - 2*cd.off[i];

        Initialize(pad,sz_div);
    }

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

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

    CellDecomposer_sm(const CellDecomposer_sm<dim,T,transform> & cd, Box<dim,size_t> & ext)
    :t(Matrix<dim,T>::identity(),cd.getOrig())
    {
        setDimensions(cd,ext);
    }


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

    }

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

    inline const Matrix<dim,T> & getMat() const
    {
        return t.getMat();
    }

    inline const Point<dim,T> & getOrig() const
    {
        return t.getOrig();
    }

    inline Box<dim,long int> convertDomainSpaceIntoCellUnits(const Box<dim,T> & b_d, const size_t (& bc)[dim]) const
    {
        Box<dim,long int> g_box;
        Box<dim,T> b = b_d;
        b -= getOrig();

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

        // Considering that we are interested in a box decomposition of the space
        // where each box does not intersect any other boxes in the decomposition we include the negative
        // countour and exclude the positive one. So ceilP1 do the job for P1 while ceilP2 - 1
        // do the job for P2

        b.floorP1();
        b.ceilP2();

        g_box = b;

        // Translate the box by the offset

        for (size_t i = 0 ; i < dim ; i++)
        {
            g_box.setLow(i,g_box.getLow(i) + off[i]);
            g_box.setHigh(i,g_box.getHigh(i) + off[i]);
        }

        // on the other hand with non periodic boundary condition, the positive border of the
        // sub-domain at the edge of the domain must be included

        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 in the decomposition)
            if (b_d.getHigh(i) == box.getHigh(i))
            {
                if (bc[i] == NON_PERIODIC)
                {
                    // Here the positive boundary is included
                    g_box.setHigh(i,gr_cell.size(i) - off[i]);
                }
                else
                {
                    // Carefull in periodic gr_cell is one bigger than the non-periodic
                    // and the positive boundary is excluded
                    g_box.setHigh(i,gr_cell.size(i)-1 - off[i]);
                }
            }


            if (b_d.getLow(i) == box.getHigh(i))
            {
                if (bc[i] == NON_PERIODIC)
                {
                    // The instruction is the same but the meaning is different
                    // for this reason there is anyway a branch
                    // Here the border is not included
                    g_box.setLow(i,gr_cell.size(i) - off[i]);
                }
                else
                {
                    // Carefull in periodic gr_cell is one bigger than the non-periodic
                    // Here the border is included
                    g_box.setLow(i,gr_cell.size(i) - off[i]);
                }
            }


            // we are at the positive border (We are assuming that there are not rounding error in the decomposition)
            /* coverity[copy_paste_error] */
            if (b_d.getHigh(i) == box.getLow(i))
                    g_box.setHigh(i,off[i]);


            if (b_d.getLow(i) == box.getLow(i))
                g_box.setLow(i,off[i]);
        }

        return g_box;
    }


    inline Box<dim,long int> convertDomainSpaceIntoGridUnits(const Box<dim,T> & b_d, const size_t (& bc)[dim]) const
    {
        Box<dim,long int> g_box;
        Box<dim,T> b = b_d;
        b -= getOrig();

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

        // Considering that we are interested in a box decomposition of the space
        // where each box does not intersect any other boxes in the decomposition we include the negative
        // countour and exclude the positive one. So ceilP1 do the job for P1 while ceilP2 - 1
        // do the job for P2

        b.ceilP1();

        // (we do -1 later)
        b.ceilP2();
        for (size_t i = 0 ; i < dim ; i++)  {b.setHigh(i,b.getHigh(i)-1);}

        g_box = b;

        // on the other hand with non periodic boundary condition, the positive border of the
        // sub-domain at the edge of the domain must be included

        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 in the decomposition)
            if (b_d.getHigh(i) == box.getHigh(i))
            {
                if (bc[i] == NON_PERIODIC)
                {
                    // Here the positive boundary is included
                    g_box.setHigh(i,gr_cell.size(i) - off[i]);
                }
                else
                {
                    // Carefull in periodic gr_cell is one bigger than the non-periodic
                    // and the positive boundary is excluded
                    g_box.setHigh(i,gr_cell.size(i)-1 - off[i]);
                }
            }


            if (b_d.getLow(i) == box.getHigh(i))
            {
                if (bc[i] == NON_PERIODIC)
                {
                    // The instruction is the same but the meaning is different
                    // for this reason there is anyway a branch
                    // Here the border is not included
                    g_box.setLow(i,gr_cell.size(i) - off[i]);
                }
                else
                {
                    // Carefull in periodic gr_cell is one bigger than the non-periodic
                    // Here the border is included
                    g_box.setLow(i,gr_cell.size(i) - off[i]);
                }
            }


            // we are at the positive border (We are assuming that there are not rounding error in the decomposition)
            /* coverity[copy_paste_error] */
            if (b_d.getHigh(i) == box.getLow(i))
                    g_box.setHigh(i,off[i]);


            if (b_d.getLow(i) == box.getLow(i))
                g_box.setLow(i,off[i]);
        }

        return g_box;
    }


    inline Box<dim,T> convertCellUnitsIntoDomainSpace(const Box<dim,long int> & b_d) const
    {
        Box<dim,T> be;

        for (size_t i = 0 ; i < dim ; i++)
        {
            if ((long int)gr_cell.size(i) - (long int)off[i] == b_d.getLow(i))
                be.setLow(i,box.getHigh(i));
            else if ((long int)off[i] == b_d.getLow(i))
                be.setLow(i,box.getLow(i));
            else
                be.setLow(i,(b_d.getLow(i) - off[i]) * box_unit.getP2()[i] + box.getLow(i));

            if ((long int)gr_cell.size(i) - (long int)off[i] == b_d.getHigh(i))
                be.setHigh(i,box.getHigh(i));
            else if ((long int)off[i] == b_d.getHigh(i))
                be.setHigh(i,box.getLow(i));
            else
                be.setHigh(i,(b_d.getHigh(i) - off[i]) * box_unit.getP2()[i] + box.getLow(i));
        }

        return be;
    }


    inline Box<dim,T> convertCellUnitsIntoDomainSpaceMiddle(const Box<dim,long int> & b_d) const
    {
        Box<dim,T> be;

        for (size_t i = 0 ; i < dim ; i++)
        {
            if ((long int)gr_cell.size(i) - (long int)off[i] == b_d.getLow(i))
            {be.setLow(i,box.getHigh(i));}
            else if ((long int)off[i] == b_d.getLow(i))
            {be.setLow(i,box.getLow(i));}
            else
            {be.setLow(i,(b_d.getLow(i) - off[i]) * box_unit.getP2()[i] + box.getLow(i) - box_unit.getP2()[i] / 2.0);}

            if ((long int)gr_cell.size(i) - (long int)off[i] == b_d.getHigh(i))
            {be.setHigh(i,box.getHigh(i));}
            else if ((long int)off[i] == b_d.getHigh(i))
            {be.setHigh(i,box.getLow(i));}
            else
            {be.setHigh(i,(b_d.getHigh(i) - off[i]) * box_unit.getP2()[i] + box.getLow(i) + box_unit.getP2()[i] / 2.0);}
        }

        return be;
    }

    const size_t (& getDiv() const)[dim]
    {
        return gr_cell.getSize();
    }


    const size_t (& getDivWP() const)[dim]
    {
        for (size_t i = 0 ; i < dim ; i++)
        {div_wp[i] = gr_cell.getSize()[i] - 2*getPadding(i);}

        return div_wp;
    }

    const Box<dim,T> & getDomain() const
    {
        return box;
    }

    inline void swap(CellDecomposer_sm<dim,T,transform> & cd)
    {
        // swap all the members
        p_middle.swap(cd.p_middle);

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

        // Total number of cell
        size_t tot_n_cell_t = tot_n_cell;
        tot_n_cell = cd.tot_n_cell;
        cd.tot_n_cell = tot_n_cell_t;

        box.swap(cd.box);
        box_unit.swap(cd.box_unit);
        gr_cell.swap(cd.gr_cell);
        gr_cell2.swap(cd.gr_cell2);

        for (size_t i = 0 ; i < dim ; i++)
        {
            size_t off_t = off[i];
            off[i] = cd.off[i];
            cd.off[i] = off_t;

            size_t cs_t = cell_shift.get(i);
            cell_shift.get(i) = cd.cell_shift.get(i);
            cd.cell_shift.get(i) = cs_t;
        }
    }

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

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

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

            if (cell_shift.get(i) != cd.cell_shift.get(i))
                return false;
        }

        return true;
    }

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

    size_t getPadding(size_t i) const
    {
        return off[i];
    }

    size_t (& getPadding())[dim]
    {
        return off;
    }

    grid_key_dx<dim> getStartDomainCell() const
    {
        grid_key_dx<dim> key;

        for (size_t i = 0 ; i < dim ; i++)
        {
            key.set_d(i, cell_shift.get(i));
        }

        return key;
    }

    grid_key_dx<dim> getStopDomainCell() const
    {
        grid_key_dx<dim> key;

        for (size_t i = 0 ; i < dim ; i++)
        {
            key.set_d(i,cell_shift.get(i) + gr_cell2.size(i) - 2*getPadding(i) - 1);
        }

        return key;
    }

    grid_key_dx<dim> getShift() const
    {
        grid_key_dx<dim> k;

        for (size_t i = 0 ; i < dim ; i++)
            k.set_d(i,cell_shift.get(i));

        return k;
    }

    const grid_sm<dim,void> & getInternalGrid() const
    {
        return gr_cell2;
    }
};


#endif /* CELLDECOMPOSER_HPP_ */