VerletListFast.hpp

/*
 * VerletListFast.hpp
 *
 *  Created on: Aug 16, 2016
 *      Author: i-bird
 */

#ifndef OPENFPM_DATA_SRC_NN_VERLETLIST_VERLETLISTFAST_HPP_
#define OPENFPM_DATA_SRC_NN_VERLETLIST_VERLETLISTFAST_HPP_

#include "VerletNNIterator.hpp"
#include "NN/CellList/CellList_util.hpp"
#include "NN/Mem_type/MemFast.hpp"
#include "NN/Mem_type/MemBalanced.hpp"
#include "NN/Mem_type/MemMemoryWise.hpp"

#define VERLET_STARTING_NSLOT 128


#define WITH_RADIUS 3

template<unsigned int dim, typename T, typename CellListImpl, typename PartIt, int type, typename local_index>
struct NNType
{
    template<typename vector_pos_type>
    static inline auto get(const PartIt & it,
                           const vector_pos_type & v,
                           Point<dim,T> & xp,
                           size_t p,
                           CellListImpl & cl,
                           T r_cut) -> decltype(cl.template getNNIterator<NO_CHECK>(0))
    {
        return cl.template getNNIterator<NO_CHECK>(cl.getCell(xp));
    }

    static inline void add(local_index p, openfpm::vector<local_index> & pc)
    {
    }
};


template<unsigned int dim, typename T, typename CellListImpl, typename PartIt, typename local_index>
struct NNType<dim,T,CellListImpl,PartIt,WITH_RADIUS,local_index>
{
    template<typename vector_pos_type>
    static inline auto get(const PartIt & it,
                           const vector_pos_type & v,
                           Point<dim,T> & xp,
                           size_t p,
                           CellListImpl & cl,
                           T r_cut) -> decltype(cl.template getNNIteratorRadius<NO_CHECK>(0,0.0))
    {
        return cl.template getNNIteratorRadius<NO_CHECK>(cl.getCell(xp),r_cut);
    }

    static inline void add(local_index p, openfpm::vector<local_index> & pc)
    {
    }
};


template<unsigned int dim, typename T, typename CellListImpl, typename PartIt, typename local_index>
struct NNType<dim,T,CellListImpl,PartIt,VL_SYMMETRIC,local_index>
{
    template<typename vector_pos_type>
    static inline auto get(const PartIt & it,
                           const vector_pos_type & v,
                           Point<dim,T> & xp,
                           size_t p,
                           CellListImpl & cl,
                           T r_cut) -> decltype(cl.template getNNIteratorSym<NO_CHECK>(0,0,vector_pos_type()))
    {
        return cl.template getNNIteratorSym<NO_CHECK>(cl.getCell(xp),p,v);
    }

    static inline void add(local_index p, openfpm::vector<local_index> & pc)
    {
    }
};


template<unsigned int dim, typename T, typename CellListImpl, typename PartIt, typename local_index>
struct NNType<dim,T,CellListImpl,PartIt,VL_CRS_SYMMETRIC,local_index>
{
    template<typename vector_pos_type>
    static inline auto get(const PartIt & it,
                           const vector_pos_type & v,
                           Point<dim,T> & xp,
                           size_t p,
                           CellListImpl & cl,
                           T r_cut) -> decltype(it.getNNIteratorCSR(v))
    {
        return it.getNNIteratorCSR(v);
    }


    static inline void add(local_index p, openfpm::vector<local_index> & pc)
    {
        pc.add(p);
    }
};

template<unsigned int type, unsigned int dim, typename vector, typename CellList>
class PartItNN
{
public:

    static inline auto get(const vector & pos, const openfpm::vector<size_t> & dom, const openfpm::vector<subsub_lin<dim>> & anom, CellList & cli, size_t g_m, size_t & end) -> decltype(pos.getIteratorTo(0))
    {
        end = g_m;
        return pos.getIteratorTo(end);
    }
};

template<unsigned int dim, typename vector, typename CellList>
class PartItNN<VL_CRS_SYMMETRIC,dim,vector,CellList>
{
public:

    static inline ParticleItCRS_Cells<dim,CellList,vector> get(const vector & pos, const openfpm::vector<size_t> & dom, const openfpm::vector<subsub_lin<dim>> & anom, CellList & cli, size_t g_m, size_t & end)
    {
        end = pos.size();
        return ParticleItCRS_Cells<dim,CellList,vector>(cli,dom,anom,cli.getNNc_sym());
    }
};

template<unsigned int dim,
         typename T,
         typename Mem_type = Mem_fast<HeapMemory,local_index_>,
         typename transform = no_transform<dim,T>,
         typename vector_pos_type = openfpm::vector<Point<dim,T>>,
         typename CellListImpl = CellList<dim,T,Mem_fast<HeapMemory,typename Mem_type::local_index_type>,transform,vector_pos_type> >
class VerletList: public Mem_type
{
protected:

    typename Mem_type::local_index_type slot;

    openfpm::vector<typename Mem_type::local_index_type> dp;

private:

    size_t n_dec;

    CellListImpl cli;


    void initCl(CellListImpl & cli, vector_pos_type & pos, size_t g_m, size_t opt)
    {
        mgpu::ofp_context_t context(mgpu::gpu_context_opt::dummy);
        if (opt & VL_SYMMETRIC || opt & VL_CRS_SYMMETRIC)
        {populate_cell_list(pos,cli,context,g_m,CL_SYMMETRIC,cl_construct_opt::Full);}
        else
        {populate_cell_list(pos,cli,context,g_m,CL_NON_SYMMETRIC,cl_construct_opt::Full);}
    }

    inline void create(const vector_pos_type & pos, const vector_pos_type & pos2, const openfpm::vector<size_t> & dom, const openfpm::vector<subsub_lin<dim>> & anom, T r_cut, size_t g_m, CellListImpl & cl, size_t opt)
    {
        if (opt == VL_CRS_SYMMETRIC)
        {
            create_<CellNNIteratorSym<dim,CellListImpl,vector_pos_type,RUNTIME,NO_CHECK>,VL_CRS_SYMMETRIC>(pos,pos2,dom,anom,r_cut,g_m,cl,opt);
        }
        else if (opt == VL_SYMMETRIC)
        {
            create_<decltype(cl.template getNNIteratorSym<NO_CHECK>(0,0,pos)),VL_SYMMETRIC>(pos,pos2,dom,anom,r_cut,g_m,cl,opt);
        }
        else
        {
            create_<decltype(cl.template getNNIterator<NO_CHECK>(0)),VL_NON_SYMMETRIC>(pos,pos2,dom,anom,r_cut,g_m,cl,opt);
        }
    }

    template<typename NN_type, int type> inline void create_(const vector_pos_type & pos, const vector_pos_type & pos2 , const openfpm::vector<size_t> & dom, const openfpm::vector<subsub_lin<dim>> & anom, T r_cut, size_t g_m, CellListImpl & cli, size_t opt)
    {
        size_t end;

        auto it = PartItNN<type,dim,vector_pos_type,CellListImpl>::get(pos,dom,anom,cli,g_m,end);

        Mem_type::init_to_zero(slot,end);

        dp.clear();

        // square of the cutting radius
        T r_cut2 = r_cut * r_cut;

        // iterate the particles
        while (it.isNext())
        {
            typename Mem_type::local_index_type i = it.get();
            Point<dim,T> xp = pos.template get<0>(i);

            // Get the neighborhood of the particle
            auto NN = NNType<dim,T,CellListImpl,decltype(it),type,typename Mem_type::local_index_type>::get(it,pos,xp,i,cli,r_cut);
            NNType<dim,T,CellListImpl,decltype(it),type,typename Mem_type::local_index_type>::add(i,dp);

            while (NN.isNext())
            {
                auto nnp = NN.get();

                Point<dim,T> xq = pos2.template get<0>(nnp);

                if (xp.distance2(xq) < r_cut2)
                    addPart(i,nnp);

                // Next particle
                ++NN;
            }

            ++it;
        }
    }

    inline void createR(openfpm::vector<Point<dim,T>> & pos, T r_cut, size_t g_m, CellListImpl & cl)
    {

        Mem_type::init_to_zero(slot,g_m);

        // square of the cutting radius
        T r_cut2 = r_cut * r_cut;

        // iterate the particles
        for (size_t i = 0 ; i < g_m ; i++)
        {
            Point<dim,T> p = pos.template get<0>(i);

            // Get the neighborhood of the particle
            auto NN = cl.template getNNIteratorRadius<NO_CHECK>(cl.getCell(p),r_cut);
            while (NN.isNext())
            {
                auto nnp = NN.get();

                Point<dim,T> q = pos.template get<0>(nnp);

                if (p.distance2(q) < r_cut2)
                    addPart(i,nnp);

                // Next particle
                ++NN;
            }
        }
    }

public:

    typedef Mem_type Mem_type_type;

    typedef size_t value_type;

    typedef CellListImpl CellListImpl_;

    size_t size()
    {
        return Mem_type::size();
    }

    inline void addPart(size_t part_id, size_t ele)
    {
        Mem_type::addCell(part_id,ele);
    }

    void Initialize(const Box<dim,T> & box, const Box<dim,T> & dom, T r_cut, vector_pos_type & pos, size_t g_m, size_t opt = VL_NON_SYMMETRIC)
    {
        // Number of divisions
        size_t div[dim];

        Box<dim,T> bt = box;

        // Calculate the divisions for the Cell-lists
        cl_param_calculate(bt,div,r_cut,Ghost<dim,T>(0.0));

        // Initialize a cell-list
        cli.Initialize(bt,div);

        initCl(cli,pos,g_m,opt);

        // Unuseful empty vector
        openfpm::vector<subsub_lin<dim>> anom_c;
        openfpm::vector<size_t> dom_c;

        // create verlet
        create(pos, pos,dom_c,anom_c,r_cut,g_m,cli,opt);
    }

    void InitializeSym(const Box<dim,T> & box, const Box<dim,T> & dom, const Ghost<dim,T> & g, T r_cut, openfpm::vector<Point<dim,T>> & pos, size_t g_m)
    {
        // Padding
        size_t pad = 0;

        // Cell decomposer
        CellDecomposer_sm<dim,T,shift<dim,T>> cd_sm;

        // Calculate the divisions for the Cell-lists
        cl_param_calculateSym<dim,T>(box,cd_sm,g,r_cut,pad);

        // Initialize a cell-list
        cli.Initialize(cd_sm,dom,pad);
        initCl(cli,pos,g_m,VL_SYMMETRIC);

        // Unused
        openfpm::vector<subsub_lin<dim>> anom_c;
        openfpm::vector<size_t> dom_c;

        // create verlet
        create(pos, pos,dom_c,anom_c,r_cut,g_m,cli,VL_SYMMETRIC);
    }


    void InitializeCrs(const Box<dim,T> & box, const Box<dim,T> & dom, const Ghost<dim,T> & g, T r_cut, openfpm::vector<Point<dim,T>> & pos, size_t g_m)
    {
        // Padding
        size_t pad = 0;

        // Cell decomposer
        CellDecomposer_sm<dim,T,shift<dim,T>> cd_sm;

        // Calculate the divisions for the Cell-lists
        cl_param_calculateSym<dim,T>(box,cd_sm,g,r_cut,pad);

        // Initialize a cell-list
        cli.Initialize(cd_sm,dom,pad);
        initCl(cli,pos,g_m,VL_SYMMETRIC);
    }

    void createVerletCrs(T r_cut, size_t g_m, openfpm::vector<Point<dim,T>> & pos, openfpm::vector<size_t> & dom_c, openfpm::vector<subsub_lin<dim>> & anom_c)
    {
        // create verlet
        create(pos, pos,dom_c,anom_c,r_cut,g_m,cli,VL_CRS_SYMMETRIC);
    }

    void update(const Box<dim,T> & dom, T r_cut, openfpm::vector<Point<dim,T>> & pos, size_t & g_m, size_t opt)
    {
        initCl(cli,pos,g_m,opt);

        // Unused
        openfpm::vector<subsub_lin<dim>> anom_c;
        openfpm::vector<size_t> dom_c;

        create(pos, pos,dom_c,anom_c,r_cut,g_m,cli,opt);
    }

    void updateCrs(const Box<dim,T> & dom, T r_cut, openfpm::vector<Point<dim,T>> & pos, size_t & g_m, const openfpm::vector<size_t> & dom_c, const openfpm::vector<subsub_lin<dim>> & anom_c)
    {
        initCl(cli,pos,g_m,VL_CRS_SYMMETRIC);

        create(pos,pos,dom_c,anom_c,r_cut,g_m,cli,VL_CRS_SYMMETRIC);
    }

    void Initialize(CellListImpl & cli,
                    T r_cut,
                    const vector_pos_type & pos,
                    const vector_pos_type & pos2,
                    size_t g_m, size_t opt = VL_NON_SYMMETRIC)
    {
        Point<dim,T> spacing = cli.getCellBox().getP2();

        // Create with radius or not
        bool wr = true;

        for (size_t i = 0 ; i < dim ; i++)
            wr &= r_cut <= spacing.get(i);

        if (wr == true || opt == VL_SYMMETRIC)
        {
            openfpm::vector<subsub_lin<dim>> anom_c;
            openfpm::vector<size_t> dom_c;

            create(pos,pos2,dom_c,anom_c,r_cut,g_m,cli,opt);
        }
        else
        {
            openfpm::vector<subsub_lin<dim>> anom_c;
            openfpm::vector<size_t> dom_c;

            create_<decltype(cli.template getNNIteratorRadius<NO_CHECK>(0,0.0)),WITH_RADIUS>(pos,pos2,dom_c,anom_c,r_cut,g_m,cli,VL_NON_SYMMETRIC);
        }
    }

    VerletList()
    :Mem_type(VERLET_STARTING_NSLOT),slot(VERLET_STARTING_NSLOT),n_dec(0)
    {};

    VerletList(const VerletList<dim,T,Mem_type,transform,vector_pos_type,CellListImpl> & cell)
    :Mem_type(VERLET_STARTING_NSLOT),slot(VERLET_STARTING_NSLOT)
    {
        this->operator=(cell);
    }

    VerletList(VerletList<dim,T,Mem_type,transform,vector_pos_type,CellListImpl> && cell)
    :Mem_type(VERLET_STARTING_NSLOT),slot(VERLET_STARTING_NSLOT),n_dec(0)
    {
        this->operator=(cell);
    }


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

    VerletList(Box<dim,T> & box, T r_cut, openfpm::vector<Point<dim,T>> & pos, size_t g_m, size_t slot=VERLET_STARTING_NSLOT)
    :slot(slot)
    {
        SpaceBox<dim,T> sbox(box);
        Initialize(sbox,r_cut,pos,g_m);
    }

    VerletList(SpaceBox<dim,T> & box, Box<dim,T> & dom, T r_cut, openfpm::vector<Point<dim,T>> & pos, size_t g_m, size_t slot=VERLET_STARTING_NSLOT)
    :slot(slot)
    {
        Initialize(box,r_cut,pos);
    }

    ~VerletList()
    {}

    VerletList<dim,T,Mem_type,transform,vector_pos_type,CellListImpl> &
    operator=(VerletList<dim,T,Mem_type,transform,vector_pos_type,CellListImpl> && vl)
    {
        slot = vl.slot;

        Mem_type::operator=(vl);
        dp.swap(vl.dp);

        n_dec = vl.n_dec;

        return *this;
    }

    VerletList<dim,T,Mem_type,transform,vector_pos_type,CellListImpl> &
    operator=(const VerletList<dim,T,Mem_type,transform,vector_pos_type,CellListImpl> & vl)
    {
        slot = vl.slot;

        Mem_type::operator=(vl);

        cli = vl.cli;

        dp = vl.dp;
        n_dec = vl.n_dec;

        return *this;
    }

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

    inline size_t get(size_t i, size_t j) const
    {
        return Mem_type::get(i,j);
    }

    inline void swap(VerletList<dim,T,Mem_type,transform,vector_pos_type,CellListImpl> & vl)
    {
        Mem_type::swap(vl);
        dp.swap(vl.dp);

        size_t vl_slot_tmp = vl.slot;
        vl.slot = slot;
        slot = vl_slot_tmp;

        cli.swap(vl.cli);

        size_t n_dec_tmp = vl.n_dec;
        vl.n_dec = n_dec;
        n_dec = n_dec_tmp;
    }

    template<unsigned int impl=NO_CHECK>
    inline VerletNNIterator<dim,VerletList<dim,T,Mem_type,transform,vector_pos_type,CellListImpl>>
    getNNIterator(size_t part_id)
    {
        VerletNNIterator<dim,VerletList<dim,T,Mem_type,transform,vector_pos_type,CellListImpl>> vln(part_id,*this);

        return vln;
    }

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

    inline const typename Mem_type::local_index_type &
    getStart(typename Mem_type::local_index_type part_id)
    {
        return Mem_type::getStartId(part_id);
    }

    inline const typename Mem_type::local_index_type &
    getStop(typename Mem_type::local_index_type part_id)
    {
        return Mem_type::getStopId(part_id);
    }

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

    CellListImpl & getInternalCellList()
    {
        return cli;
    }

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

        cli.set_ndec(n_dec);
    }

    size_t get_ndec()
    {
        return n_dec;
    }

    openfpm::vector<typename Mem_type::local_index_type> & getParticleSeq()
    {
        return dp;
    }
};



#endif /* OPENFPM_DATA_SRC_NN_VERLETLIST_VERLETLISTFAST_HPP_ */