CellListFast.hpp

/*
 * CellListStandard.hpp
 *
 *  Created on: Mar 22, 2015
 *      Author: Pietro Incardona
 */

#ifndef CELLLISTSTANDARD_HPP_
#define CELLLISTSTANDARD_HPP_

#include "CellDecomposer.hpp"
#include "Space/SpaceBox.hpp"
#include "util/mathutil.hpp"
#include "CellNNIterator.hpp"
#include "Space/Shape/HyperCube.hpp"

#include "util/common.hpp"

template<unsigned int dim, typename T, typename transform, typename base>
class CellList<dim,T,FAST,transform,base> : public CellDecomposer_sm<dim,T,transform>
{
    // The array contain the neighborhood of the cell-id in case of asymmetric interaction
    //
    //    * * *
    //    * x *
    //    * * *

    long int NNc_full[openfpm::math::pow(3,dim)];

    // The array contain the neighborhood of the cell-id in case of symmetric interaction
    //
    //   * * *
    //     x *
    //
    long int NNc_sym[openfpm::math::pow(3,dim)/2+1];

    // The array contain the neighborhood of the cell-id in case of symmetric interaction (Optimized)
    //
    //   * *
    //   x *
    //
    long int NNc_cr[openfpm::math::pow(2,dim)];

    // Number of slot for each cell
    size_t slot;

    // number of particle in each cell list
    openfpm::vector<size_t> cl_n;

    // elements that each cell store (each cell can store a number
    // of elements == slot )
    base cl_base;

    //Origin point
    Point<dim,T> orig;

    void realloc()
    {
        // we do not have enough slots reallocate the basic structure with more
        // slots
        base cl_base_(2*slot * cl_n.size());

        // copy cl_base
        for (size_t i = 0 ; i < cl_n.size() ; i++)
        {
            for (size_t j = 0 ; j < cl_n.get(i) ; j++)
                cl_base_.get(2*i*slot + j) = cl_base.get(slot * i + j);
        }

        // Double the number of slots
        slot *= 2;

        // swap the memory
        cl_base.swap(cl_base_);
    }

public:

    // Object type that the structure store
    typedef T value_type;

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

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

        // Initialize point transformation

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

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

        CellDecomposer_sm<dim,T,transform>::setDimensions(box,div, mat, orig, pad);
        this->slot = slot;
        this->orig = orig;

        // create the array that store the number of particle on each cell and se it to 0

        cl_n.resize(this->tot_n_cell);
        cl_n.fill(0);

        // create the array that store the cell id

        cl_base.resize(this->tot_n_cell * slot);

        // 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;
        typedef typename generate_array<size_t,dim, Fill_one>::result NNone;

        // Generate the sub-grid iterator

        grid_key_dx_iterator_sub<dim> gr_sub3(this->gr_cell,NNzero::data,NNtwo::data);

        // Calculate the NNc array

        size_t middle = this->gr_cell.LinId(NNone::data);
        size_t i = 0;
        while (gr_sub3.isNext())
        {
            NNc_full[i] = (long int)this->gr_cell.LinId(gr_sub3.get()) - middle;

            ++gr_sub3;
            i++;
        }

        // Calculate the NNc_sym array

        i = 0;
        gr_sub3.reset();
        while (gr_sub3.isNext())
        {
            auto key = gr_sub3.get();

            size_t lin = this->gr_cell.LinId(key);

            // Only the first half is considered
            if (lin < middle)
            {
                ++gr_sub3;
                continue;
            }

            NNc_sym[i] = lin - middle;

            ++gr_sub3;
            i++;
        }

        // Calculate the NNc_cross array

        i = 0;
        grid_key_dx_iterator_sub<dim> gr_sub2(this->gr_cell,NNzero::data,NNone::data);

        while (gr_sub2.isNext())
        {
            auto key = gr_sub2.get();

            NNc_cr[i] = (long int)this->gr_cell.LinId(key);

            ++gr_sub2;
            i++;
        }
    }

    CellList()
    :slot(16)
    {};

    CellList(const CellList<dim,T,FAST,transform,base> & cell)
    {
        this->operator=(cell);
    }

    CellList(CellList<dim,T,FAST,transform,base> && cell)
    {
        this->operator=(cell);
    }

    CellList<dim,T,FAST,transform,base> & operator=(CellList<dim,T,FAST,transform,base> && 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]);
        std::copy(&cell.NNc_cr[0],&cell.NNc_cr[openfpm::math::pow(2,dim)],&NNc_cr[0]);

        size_t tslot = slot;
        slot = cell.slot;
        cell.slot = tslot;

        cl_n.swap(cell.cl_n);
        cl_base.swap(cell.cl_base);

        Point<dim,T> torig = orig;
        orig = cell.orig;
        cell.orig = torig;

        return this;
    }

    CellList<dim,T,FAST,transform,base> & operator=(const CellList<dim,T,FAST,transform,base> & 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]);
        std::copy(&cell.NNc_cr[0],&cell.NNc_cr[openfpm::math::pow(2,dim)],&NNc_cr[0]);

        slot = cell.slot;

        cl_n = cell.cl_n;
        cl_base = cell.cl_base;

        orig = cell.orig;

        return *this;
    }

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

    CellList(Box<dim,T> & box, const size_t (&div)[dim], Point<dim,T> & orig, const size_t pad = 1, size_t slot=16)
    {
        SpaceBox<dim,T> sbox(box);
        Initialize(sbox,div,orig,pad,slot);
    }

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


    inline void addCell(size_t cell_id, typename base::value_type ele)
    {
        // Get the number of element the cell is storing

        size_t nl = getNelements(cell_id);

        if (nl + 1 >= slot)
        {
            realloc();
        }

        // we have enough slot to store another neighbor element

        cl_base.get(slot * cell_id + cl_n.get(cell_id)) = ele;
        cl_n.get(cell_id)++;
    }

    inline void add(const T (& pos)[dim], typename base::value_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 add(const Point<dim,T> & pos, typename base::value_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)
    {
        cl_n.get(cell)--;
    }

    inline size_t getNelements(const size_t cell_id) const
    {
        return cl_n.get(cell_id);
    }

    inline auto get(size_t cell, size_t ele) -> decltype(cl_base.get(cell * slot + ele))
    {
        return cl_base.get(cell * slot + ele);
    }

    template<unsigned int i> inline auto get(size_t cell, size_t ele) -> decltype(cl_base.get(cell * slot + ele))
    {
        return cl_base.template get<i>(cell * slot + ele);
    }

    inline void swap(CellList<dim,T,FAST,transform,base> & cl)
    {
        cl_n.swap(cl.cl_n);
        cl_base.swap(cl.cl_base);
    }

    CellIterator<CellList<dim,T,FAST,transform,base>> getIterator(size_t cell)
    {
        return CellIterator<CellList<dim,T,FAST,transform,base>>(cell,*this);
    }

    template<unsigned int impl> inline CellNNIterator<dim,CellList<dim,T,FAST,transform,base>,FULL,impl> getNNIterator(size_t cell)
    {
        CellNNIterator<dim,CellList<dim,T,FAST,transform,base>,FULL,impl> cln(cell,NNc_full,*this);

        return cln;
    }


    template<unsigned int impl> inline CellNNIterator<dim,CellList<dim,T,FAST,transform,base>,SYM,impl> getNNIteratorSym(size_t cell)
    {
        CellNNIterator<dim,CellList<dim,T,FAST,transform,base>,SYM,impl> cln(cell,NNc_sym,*this);

        return cln;
    }


    template<unsigned int impl> inline CellNNIterator<dim,CellList<dim,T,FAST,transform,base>,CRS,impl> getNNIteratorCross(size_t cell)
    {
        CellNNIterator<dim,CellList<dim,T,FAST,transform,base>,CRS,impl> cln(cell,NNc_cr,*this);

        return cln;
    }
};


#endif /* CELLLISTSTANDARD_HPP_ */