CellList.hpp

/*
 * CellList.hpp
 *
 *  Created on: Mar 21, 2015
 *      Author: Pietro Incardona
 */

#ifndef CELLLIST_HPP_
#define CELLLIST_HPP_

#include "CellList_def.hpp"
#include "Vector/map_vector.hpp"
#include "CellDecomposer.hpp"
#include "Space/SpaceBox.hpp"
#include "util/mathutil.hpp"
#include "CellNNIterator.hpp"
#include "Space/Shape/HyperCube.hpp"
#include "CellListNNIteratorRadius.hpp"
#include <unordered_map>

#include "CellListIterator.hpp"
#include "ParticleIt_Cells.hpp"
#include "ParticleItCRS_Cells.hpp"
#include "util/common.hpp"

#include "NN/Mem_type/MemFast.hpp"
#include "NN/Mem_type/MemBalanced.hpp"
#include "NN/Mem_type/MemMemoryWise.hpp"
#include "NN/CellList/NNc_array.hpp"
#include "cuda/CellList_cpu_ker.cuh"

template<typename key,typename val>
class wrap_unordered_map: public std::unordered_map<key,val>
{
};

#ifdef HAVE_LIBQUADMATH

#include <boost/multiprecision/float128.hpp>


template<typename val>
class wrap_unordered_map<boost::multiprecision::float128,val>
{
};

#endif

#define CELL_REALLOC 16ul

#define STARTING_NSLOT 16

template<unsigned int dim> void NNcalc_csr(openfpm::vector<std::pair<grid_key_dx<dim>,grid_key_dx<dim>>> & cNN)
{
    // Calculate the NNc_full array, it is a structure to get the neighborhood array

    // compile-time array {0,0,0,....}  {2,2,2,...} {1,1,1,...}

    typedef typename generate_array<size_t,dim, Fill_zero>::result NNzero;
    typedef typename generate_array<size_t,dim, Fill_two>::result NNtwo;

    // Generate the sub-grid iterator

    size_t div[dim];

    // Calculate the divisions

    for (size_t i = 0 ; i < dim ; i++)
        div[i] = 4;

    grid_sm<dim,void> gs(div);
    grid_key_dx_iterator_sub<dim> gr_sub3(gs,NNzero::data,NNtwo::data);

    grid_key_dx<dim> src_; // Source cell
    for (size_t i = 0; i < dim; i++)
        src_.set_d(i,1);

    size_t middle = gs.LinId(src_);

    // Calculate the symmetric crs array
    while (gr_sub3.isNext())
    {
        auto dst = gr_sub3.get();
        grid_key_dx<dim> src = src_;

        if ((long int)middle > gs.LinId(dst))
        {
            ++gr_sub3;
            continue;
        }

        // Here we adjust src and dst to be in the positive quadrant

        for (size_t i = 0 ; i < dim; i++)
        {
            if (dst.get(i) == 0)
            {
                src.set_d(i,src.get(i) + 1);
                dst.set_d(i,dst.get(i) + 1);
            }
        }

        src -= src_;
        dst -= src_;

        cNN.add(std::pair<grid_key_dx<dim>,grid_key_dx<dim>>(src,dst));

        ++gr_sub3;

    }
};

template<unsigned int dim> void NNcalc_sym(openfpm::vector<grid_key_dx<dim>> & cNN)
{
    // Calculate the NNc_full array, it is a structure to get the neighborhood array

    // compile-time array {0,0,0,....}  {2,2,2,...} {1,1,1,...}

    typedef typename generate_array<size_t,dim, Fill_zero>::result NNzero;
    typedef typename generate_array<size_t,dim, Fill_two>::result NNtwo;

    // Generate the sub-grid iterator

    size_t div[dim];

    // Calculate the divisions

    for (size_t i = 0 ; i < dim ; i++)
        div[i] = 4;

    grid_sm<dim,void> gs(div);
    grid_key_dx_iterator_sub<dim> gr_sub3(gs,NNzero::data,NNtwo::data);

    grid_key_dx<dim> src_; // Source cell
    for (size_t i = 0; i < dim; i++)
        src_.set_d(i,1);

    size_t middle = gs.LinId(src_);

    // Calculate the symmetric array
    while (gr_sub3.isNext())
    {
        auto dst = gr_sub3.get();

        if ((long int)middle > gs.LinId(dst))
        {
            ++gr_sub3;
            continue;
        }

        cNN.add(dst - src_);

        ++gr_sub3;

    }
};

template<unsigned int dim> void NNcalc_full(openfpm::vector<grid_key_dx<dim>> & cNN)
{
    // Calculate the NNc_full array, it is a structure to get the neighborhood array

    // compile-time array {0,0,0,....}  {2,2,2,...} {1,1,1,...}

    typedef typename generate_array<size_t,dim, Fill_zero>::result NNzero;
    typedef typename generate_array<size_t,dim, Fill_two>::result NNtwo;

    // Generate the sub-grid iterator

    size_t div[dim];

    // Calculate the divisions

    for (size_t i = 0 ; i < dim ; i++)
    {div[i] = 4;}

    grid_sm<dim,void> gs(div);
    grid_key_dx_iterator_sub<dim> gr_sub3(gs,NNzero::data,NNtwo::data);

    grid_key_dx<dim> src_; // Source cell
    for (size_t i = 0; i < dim; i++)
        src_.set_d(i,1);

    // Calculate the symmetric crs array
    while (gr_sub3.isNext())
    {
        auto dst = gr_sub3.get();

        cNN.add(dst - src_);

        ++gr_sub3;
    }
};


template<unsigned int dim, typename T>
void NNcalc_rad(T r_cut, openfpm::vector<long int> & NNcell, const Box<dim,T> & cell_box, const grid_sm<dim,void> & gs)
{
    size_t n_cell[dim];
    size_t n_cell_mid[dim];

    Point<dim,T> spacing = cell_box.getP2();

    for (size_t i = 0 ; i < dim ; i++)
    {
        n_cell[i] = 2*(std::ceil(r_cut / spacing.get(i)))+1;
        n_cell_mid[i] = n_cell[i] / 2;
    }

    grid_sm<dim,void> gsc(n_cell);
    grid_key_dx_iterator<dim> gkdi(gsc);

    Box<dim,T> cell_zero;

    for (unsigned int i = 0 ; i < dim ; i++)
    {
        cell_zero.setLow(i,n_cell_mid[i]*spacing.get(i));
        cell_zero.setHigh(i,(n_cell_mid[i]+1)*spacing.get(i));
    }

    NNcell.clear();
    while (gkdi.isNext())
    {
        auto key = gkdi.get();

        Box<dim,T> cell;

        for (unsigned int i = 0 ; i < dim ; i++)
        {
            cell.setLow(i,key.get(i)*spacing.get(i));
            cell.setHigh(i,(key.get(i)+1)*spacing.get(i));
        }

        // here we check if the cell is in the radius.
        T min_distance = cell.min_distance(cell_zero);
        if (min_distance > r_cut)
        {++gkdi;continue;}

        for (size_t i = 0 ; i < dim ; i++)
        {key.set_d(i,key.get(i) - n_cell_mid[i]);}

        NNcell.add(gs.LinId(key));

        ++gkdi;
    }
}

/* NOTE all the implementations
 *
 * has complexity O(1) in getting the cell id and the elements in a cell
 * but with different constants
 *
 */

template<unsigned int dim, typename T,  typename Mem_type, typename transform = no_transform<dim,T>, typename vector_pos_type = openfpm::vector<Point<dim,T>>>
class CellList : public CellDecomposer_sm<dim,T,transform>, public Mem_type
{
protected:
    //
    //    * * *
    //    * x *
    //    * * *


    NNc_array<dim,(unsigned int)openfpm::math::pow(3,dim)> NNc_full;

    //
    //   * * *
    //     x *
    //
//  long int NNc_sym[openfpm::math::pow(3,dim)/2+1];
    NNc_array<dim,(unsigned int)openfpm::math::pow(3,dim)/2+1> NNc_sym;

private:

    wrap_unordered_map<T,openfpm::vector<long int>> rcache;

    bool from_cd;

    size_t n_dec;

    openfpm::vector<long int> nnc_rad;



    void InitializeStructures(const size_t (& div)[dim], size_t tot_n_cell, size_t slot=STARTING_NSLOT)
    {
        Mem_type::init_to_zero(slot,tot_n_cell);

        NNc_full.set_size(div);
        NNc_full.init_full();

        NNc_sym.set_size(div);
        NNc_sym.init_sym();
    }

    void setCellDecomposer(CellDecomposer_sm<dim,T,transform> & cd, const CellDecomposer_sm<dim,T,transform> & cd_sm, const Box<dim,T> & dom_box, size_t pad) const
    {
        size_t bc[dim];

        for (size_t i = 0 ; i < dim ; i++)
        {bc[i] = NON_PERIODIC;}

        Box<dim,long int> bx = cd_sm.convertDomainSpaceIntoCellUnits(dom_box,bc);

        size_t div[dim];
        size_t div_big[dim];

        for (size_t i = 0 ; i < dim ; i++)
        {
            div[i] = bx.getHigh(i) - bx.getLow(i);
            div_big[i] = cd_sm.getDiv()[i] - 2*cd_sm.getPadding(i);
        }

        cd.setDimensions(cd_sm.getDomain(),div_big,div, pad, bx.getP1());
    }

public:

    typedef Mem_type Mem_type_type;

    typedef CellNNIteratorSym<dim,CellList<dim,T,Mem_type,transform,vector_pos_type>,vector_pos_type,RUNTIME,NO_CHECK> SymNNIterator;

    typedef typename Mem_type::local_index_type value_type;

    typedef T stype;

    typedef vector_pos_type internal_vector_pos_type;

    const grid_sm<dim,void> & getGrid()
    {
        return CellDecomposer_sm<dim,T,transform>::getGrid();
    }

    void Initialize(CellDecomposer_sm<dim,T,transform> & cd_sm, const Box<dim,T> & dom_box,const size_t pad = 1, size_t slot=STARTING_NSLOT)
    {
        size_t bc[dim];
        for (size_t i = 0 ; i < dim ; i++)  {bc[i] = NON_PERIODIC;}

        Box<dim,long int> bx = cd_sm.convertDomainSpaceIntoCellUnits(dom_box,bc);

        setCellDecomposer(*this,cd_sm,dom_box,pad);

        size_t div_w_pad[dim];
        size_t tot_cell = 1;

        for (size_t i = 0 ; i < dim ; i++)
        {
            div_w_pad[i] = bx.getHigh(i) - bx.getLow(i) + 2*pad;
            tot_cell *= div_w_pad[i];
        }

        // here we set the cell-shift for the CellDecomposer

        InitializeStructures(div_w_pad,tot_cell);

        // Initialized from CellDecomposer
        from_cd = true;
    }

    void Initialize(const Box<dim,T> & box, const size_t (&div)[dim], const size_t pad = 1, size_t slot=STARTING_NSLOT)
    {
        SpaceBox<dim,T> sbox(box);

        // Initialize point transformation

        Initialize(sbox,div,pad,slot);
    }

    void Initialize(const SpaceBox<dim,T> & box, const size_t (&div)[dim], const size_t pad = 1, size_t slot=STARTING_NSLOT)
    {
        Matrix<dim,T> mat;

        CellDecomposer_sm<dim,T,transform>::setDimensions(box,div, mat, pad);
        Mem_type::set_slot(slot);

        // create the array that store the number of particle on each cell and se it to 0
        InitializeStructures(this->gr_cell.getSize(),this->gr_cell.size());

        from_cd = false;
    }

    CellList()

    :Mem_type(STARTING_NSLOT)
    {};

    CellList(const CellList<dim,T,Mem_type,transform,vector_pos_type> & cell)
    :Mem_type(STARTING_NSLOT)
    {
        this->operator=(cell);
    }

    CellList(CellList<dim,T,Mem_type,transform,vector_pos_type> && cell)
    :Mem_type(STARTING_NSLOT)
    {
        this->operator=(cell);
    }


    CellList(Box<dim,T> & box, const size_t (&div)[dim], Matrix<dim,T> mat, const size_t pad = 1, size_t slot=STARTING_NSLOT)
    :Mem_type(slot),CellDecomposer_sm<dim,T,transform>(box,div,mat,box.getP1(),pad)
    {
        SpaceBox<dim,T> sbox(box);
        Initialize(sbox,div,pad,slot);
    }

    CellList(Box<dim,T> & box, const size_t (&div)[dim], const size_t pad = 1, size_t slot=STARTING_NSLOT)
    :Mem_type(slot),n_dec(0)
    {
        SpaceBox<dim,T> sbox(box);
        Initialize(sbox,div,pad,slot);
    }

    CellList(SpaceBox<dim,T> & box, const size_t (&div)[dim], const size_t pad = 1, size_t slot=STARTING_NSLOT)
    :Mem_type(slot)
    {
        Initialize(box,div,pad,slot);
    }

    CellList(CellDecomposer_sm<dim,T,transform> & cd_sm, const Box<dim,T> & box, const size_t pad = 1, size_t slot=STARTING_NSLOT)
    :Mem_type(slot),n_dec(0)
    {
        Initialize(cd_sm,box,pad,slot);
    }

    ~CellList()
    {}

    CellList<dim,T,Mem_type,transform> & operator=(CellList<dim,T,Mem_type,transform> && cell)
    {
        std::copy(&cell.NNc_full[0],&cell.NNc_full[openfpm::math::pow(3,dim)],&NNc_full[0]);
        std::copy(&cell.NNc_sym[0],&cell.NNc_sym[openfpm::math::pow(3,dim)/2+1],&NNc_sym[0]);

        Mem_type::swap(static_cast<Mem_type &&>(cell));

        static_cast<CellDecomposer_sm<dim,T,transform> &>(*this).swap(cell);

        n_dec = cell.n_dec;
        from_cd = cell.from_cd;

        return *this;
    }

    CellList<dim,T,Mem_type,transform,vector_pos_type> & operator=(const CellList<dim,T,Mem_type,transform,vector_pos_type> & cell)
    {
        NNc_full = cell.NNc_full;
        NNc_sym = cell.NNc_sym;

        Mem_type::operator=(static_cast<const Mem_type &>(cell));

        static_cast<CellDecomposer_sm<dim,T,transform> &>(*this) = static_cast<const CellDecomposer_sm<dim,T,transform> &>(cell);

        n_dec = cell.n_dec;
        from_cd = cell.from_cd;

        return *this;
    }

    template<typename Mem_type2>
    CellList<dim,T,Mem_type,transform,vector_pos_type> & operator=(const CellList<dim,T,Mem_type2,transform,vector_pos_type> & cell)
    {
        NNc_full = cell.private_get_NNc_full();
        NNc_sym = cell.private_get_NNc_sym();

        Mem_type::copy_general(static_cast<const Mem_type2 &>(cell));

        static_cast<CellDecomposer_sm<dim,T,transform> &>(*this) = static_cast<const CellDecomposer_sm<dim,T,transform> &>(cell);

        n_dec = cell.get_ndec();
        from_cd = cell.private_get_from_cd();

        return *this;
    }

    void setRadius(T radius)
    {
        NNcalc_rad(radius,nnc_rad,this->getCellBox(),this->getGrid());
    }

    ParticleIt_Cells<dim,CellList<dim,T,Mem_fast<>,transform>> getDomainIterator(openfpm::vector<size_t> & dom_cells)
    {
        ParticleIt_Cells<dim,CellList<dim,T,Mem_fast<>,transform>> it(*this,dom_cells);

        return it;
    }

    inline void addCell(size_t cell_id, typename Mem_type::local_index_type ele)
    {
        Mem_type::addCell(cell_id,ele);
    }

    inline void add(const T (& pos)[dim], typename Mem_type::local_index_type ele)
    {
        // calculate the Cell id
        size_t cell_id = this->getCell(pos);

        Mem_type::add(cell_id,ele);
    }

    inline void add(const Point<dim,T> & pos, typename Mem_type::local_index_type ele)
    {
        // calculate the Cell id
        size_t cell_id = this->getCell(pos);

        Mem_type::add(cell_id,ele);
    }


    inline void addDom(const T (& pos)[dim], typename Mem_type::local_index_type ele)
    {
        // calculate the Cell id

        size_t cell_id = this->getCellDom(pos);

        // add the element to the cell

        addCell(cell_id,ele);
    }

    inline void addDom(const Point<dim,T> & pos, typename Mem_type::local_index_type ele)
    {
        // calculate the Cell id

        size_t cell_id = this->getCellDom(pos);

        // add the element to the cell

        addCell(cell_id,ele);
    }

    inline void addPad(const T (& pos)[dim], typename Mem_type::local_index_type ele)
    {
        // calculate the Cell id

        size_t cell_id = this->getCellPad(pos);

        // add the element to the cell

        addCell(cell_id,ele);
    }

    inline void addPad(const Point<dim,T> & pos, typename Mem_type::local_index_type ele)
    {
        // calculate the Cell id

        size_t cell_id = this->getCell(pos);

        // add the element to the cell

        addCell(cell_id,ele);
    }

    inline void remove(size_t cell, size_t ele)
    {
        Mem_type::remove(cell,ele);
    }

    inline size_t getNCells() const
    {
        return Mem_type::size();
    }

    inline size_t getNelements(const size_t cell_id) const
    {
        return Mem_type::getNelements(cell_id);
    }

    inline auto get(size_t cell, size_t ele) -> decltype(this->Mem_type::get(cell,ele))
    {
        return Mem_type::get(cell,ele);
    }

    inline auto get(size_t cell, size_t ele) const -> decltype(this->Mem_type::get(cell,ele))
    {
        return Mem_type::get(cell,ele);
    }

    inline void swap(CellList<dim,T,Mem_type,transform,vector_pos_type> & cl)
    {
        NNc_full.swap(cl.NNc_full);
        NNc_sym.swap(cl.NNc_sym);

        Mem_type::swap(static_cast<Mem_type &>(cl));

        static_cast<CellDecomposer_sm<dim,T,transform> &>(*this).swap(static_cast<CellDecomposer_sm<dim,T,transform> &>(cl));

        n_dec = cl.n_dec;
        from_cd = cl.from_cd;
    }

    CellIterator<CellList<dim,T,Mem_type,transform>> getCellIterator(size_t cell)
    {
        return CellIterator<CellList<dim,T,Mem_type,transform>>(cell,*this);
    }

    template<unsigned int impl=NO_CHECK> inline CellNNIteratorRadius<dim,CellList<dim,T,Mem_type,transform>,impl> getNNIteratorRadius(size_t cell)
    {
        CellNNIteratorRadius<dim,CellList<dim,T,Mem_type,transform>,impl> cln(cell,nnc_rad,*this);
        return cln;
    }

    template<unsigned int impl=NO_CHECK>
    __attribute__((always_inline)) inline CellNNIterator<dim,CellList<dim,T,Mem_type,transform,vector_pos_type>,(int)FULL,impl> getNNIterator(size_t cell)
    {
        CellNNIterator<dim,CellList<dim,T,Mem_type,transform,vector_pos_type>,(int)FULL,impl> cln(cell,NNc_full,*this);
        return cln;

    }

    template<unsigned int impl=NO_CHECK>
    __attribute__((always_inline)) inline CellNNIteratorRadius<dim,CellList<dim,T,Mem_type,transform,vector_pos_type>,impl> getNNIteratorRadius(size_t cell, T r_cut)
    {
        openfpm::vector<long int> & NNc = rcache[r_cut];

        if (NNc.size() == 0)
        {NNcalc_rad(r_cut,NNc,this->getCellBox(),this->getGrid());}

        CellNNIteratorRadius<dim,CellList<dim,T,Mem_type,transform,vector_pos_type>,impl> cln(cell,NNc,*this);

        return cln;
    }



    template<unsigned int impl>
    __attribute__((always_inline)) inline CellNNIteratorSym<dim,CellList<dim,T,Mem_type,transform,vector_pos_type>,vector_pos_type,(unsigned int)SYM,impl>
    getNNIteratorSym(size_t cell, size_t p, const vector_pos_type & v)
    {
#ifdef SE_CLASS1
        if (from_cd == false)
        {std::cerr << __FILE__ << ":" << __LINE__ << " Warning when you try to get a symmetric neighborhood iterator, you must construct the Cell-list in a symmetric way" << std::endl;}
#endif

        CellNNIteratorSym<dim,CellList<dim,T,Mem_type,transform,vector_pos_type>,vector_pos_type,SYM,impl> cln(cell,p,NNc_sym,*this,v);
        return cln;
    }

    template<unsigned int impl, typename vector_pos_type2>
    __attribute__((always_inline)) inline CellNNIteratorSymMP<dim,CellList<dim,T,Mem_type,transform,vector_pos_type>,vector_pos_type2,(unsigned int)SYM,impl>
    getNNIteratorSymMP(size_t cell, size_t p, const vector_pos_type2 & v_p1, const vector_pos_type2 & v_p2)
    {
#ifdef SE_CLASS1
        if (from_cd == false)
            std::cerr << __FILE__ << ":" << __LINE__ << " Warning when you try to get a symmetric neighborhood iterator, you must construct the Cell-list in a symmetric way" << std::endl;
#endif

        CellNNIteratorSymMP<dim,CellList<dim,T,Mem_type,transform,vector_pos_type>,vector_pos_type2,SYM,impl> cln(cell,p,NNc_sym,*this,v_p1,v_p2);
        return cln;
    }

    const NNc_array<dim,(unsigned int)openfpm::math::pow(3,dim)/2+1> & getNNc_sym() const
    {
        return NNc_sym;
    }

    size_t getPadding(size_t i) const
    {
        return CellDecomposer_sm<dim,T,transform>::getPadding(i);
    }

    size_t (& getPadding())[dim]
    {
        return CellDecomposer_sm<dim,T,transform>::getPadding();
    }

    void clear()
    {
        Mem_type::clear();
    }

    void destroy()
    {
        Mem_type::destroy();
    }

    __attribute__((always_inline)) inline const typename Mem_type::local_index_type & getStartId(typename Mem_type::local_index_type cell_id) const
    {
        return Mem_type::getStartId(cell_id);
    }

    __attribute__((always_inline)) inline const typename Mem_type::local_index_type & getStopId(typename Mem_type::local_index_type cell_id) const
    {
        return Mem_type::getStopId(cell_id);
    }

    __attribute__((always_inline)) inline const typename Mem_type::local_index_type & get_lin(const typename Mem_type::local_index_type * part_id) const
    {
        return Mem_type::get_lin(part_id);
    }


    size_t g_m = 0;

    inline size_t get_gm()
    {
        return g_m;
    }

    inline void set_gm(size_t g_m)
    {
        this->g_m = g_m;
    }

#ifdef CUDA_GPU

    CellList_cpu_ker<dim,T,typename Mem_type::toKernel_type,transform> toKernel()
    {
        typedef typename Mem_type::local_index_type ids_type;
        typedef typename Mem_type::local_index_type cnt_type;

        openfpm::array<T,dim,cnt_type> spacing_c;
        openfpm::array<ids_type,dim,cnt_type> div_c;
        openfpm::array<ids_type,dim,cnt_type> off;

        for (size_t i = 0 ; i < dim ; i++)
        {
            spacing_c[i] = CellDecomposer_sm<dim,T,transform>::getCellBox().getHigh(i);
            div_c[i] = CellDecomposer_sm<dim,T,transform>::getDiv()[i];
            off[i] = CellDecomposer_sm<dim,T,transform>::getPadding(i);
        }

        CellList_cpu_ker<dim,T,typename Mem_type::toKernel_type,transform> cl(Mem_type::toKernel(),
                                                                              spacing_c,
                                                                              div_c,
                                                                              off,
                                                                              CellDecomposer_sm<dim,T,transform>::getTransform());

        return cl;
    }


    void hostToDevice()
    {
        Mem_type::hostToDevice();
    }

#endif

    const NNc_array<dim,(unsigned int)openfpm::math::pow(3,dim)> & private_get_NNc_full () const
    {
        return NNc_full;
    }

    const NNc_array<dim,(unsigned int)openfpm::math::pow(3,dim)/2+1> & private_get_NNc_sym () const
    {
        return NNc_sym;
    }

    bool private_get_from_cd() const
    {
        return from_cd;
    }


    void re_setBoxNN()
    {}


    void set_ndec(size_t n_dec)
    {
        this->n_dec = n_dec;
    }

    size_t get_ndec() const
    {
        return n_dec;
    }

};

template<unsigned int dim, typename St> static inline void cl_param_calculate(Box<dim, St> & pbox, size_t (&div)[dim], St r_cut, const Ghost<dim, St> & enlarge)
{
    // calculate the parameters of the cell list

    // extend by the ghost
    pbox.enlarge(enlarge);

    // Calculate the division array and the cell box
    for (size_t i = 0; i < dim; i++)
    {
        div[i] = static_cast<size_t>((pbox.getP2().get(i) - pbox.getP1().get(i)) / r_cut);
        div[i]++;
        pbox.setHigh(i,pbox.getLow(i) + div[i]*r_cut);
    }
}

template<unsigned int dim, typename St> static
inline void cl_param_calculateSym(const Box<dim,St> & dom,
                                  CellDecomposer_sm<dim,St,shift<dim,St>> & cd_sm,
                                  Ghost<dim,St> g,
                                  St r_cut,
                                  size_t & pad)
{
    size_t div[dim];

    for (size_t i = 0 ; i < dim ; i++)
        div[i] = (dom.getHigh(i) - dom.getLow(i)) / r_cut;

    g.magnify(1.013);

    // Calculate the maximum padding
    for (size_t i = 0 ; i < dim ; i++)
    {
        size_t tmp = std::ceil(fabs(g.getLow(i)) / r_cut);
        pad = (pad > tmp)?pad:tmp;

        tmp = std::ceil(fabs(g.getHigh(i)) / r_cut);
        pad = (pad > tmp)?pad:tmp;
    }

    cd_sm.setDimensions(dom,div,pad);
}


template<unsigned int dim, typename St>
static inline void cl_param_calculateSym(const Box<dim,St> & dom,
                                         CellDecomposer_sm<dim,St,shift<dim,St>> & cd_sm,
                                         Ghost<dim,St> g,
                                         size_t & pad)
{
    size_t div[dim];

    for (size_t i = 0 ; i < dim ; i++)
        div[i] = (dom.getHigh(i) - dom.getLow(i)) / g.getHigh(i);

    g.magnify(1.013);

    pad = 1;

    cd_sm.setDimensions(dom,div,pad);
}


#endif /* CELLLIST_HPP_ */