VerletListM.hpp

/*
 * VerletListM.hpp
 *
 *  Created on: Oct 15, 2016
 *      Author: i-bird
 */

#ifndef OPENFPM_DATA_SRC_NN_VERLETLIST_VERLETLISTM_HPP_
#define OPENFPM_DATA_SRC_NN_VERLETLIST_VERLETLISTM_HPP_

#include "NN/VerletList/VerletNNIteratorM.hpp"
#include "VerletList.hpp"


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

template<unsigned int dim, typename T, typename CellListImpl, typename PartIt>
struct NNTypeM<dim,T,CellListImpl,PartIt,VL_CRS_SYMMETRIC>
{
    static inline auto get(const PartIt & it,
                           const openfpm::vector<Point<dim,T>> & pos,
                           const openfpm::vector<pos_v<dim,T>> & v,
                           Point<dim,T> & xp,
                           size_t pp,
                           size_t p,
                           CellListImpl & cl,
                           T r_cut) -> decltype(it.getNNIteratorCSRM(pos,v))
    {
        return it.getNNIteratorCSRM(pos,v);
    }
};

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



template<unsigned int dim,
         typename T,
         unsigned int sh_byte ,
         typename CellListImpl=CellListM<dim,T,sh_byte>,
         typename transform = shift<dim,T>,
         typename VerletBase=VerletList<dim,T,Mem_fast<>,transform, size_t> >
class VerletListM : public VerletBase
{

    typedef boost::high_bit_mask_t<sh_byte>  mask_high;

    typedef boost::low_bits_mask_t<sizeof(size_t)*8-sh_byte>  mask_low;

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

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

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

/*      this->cl_n.resize(end);
        this->cl_base.resize(end*this->slot);
        this->cl_n.fill(0);*/

        this->init_to_zero(this->slot,end);

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

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

            // Get the neighborhood of the particle
            NN_type NN = NNTypeM<dim,T,CellListImpl,decltype(it),type>::get(it,pos,pos2,xp,pp,i,cli,r_cut);

            while (NN.isNext())
            {
                size_t nnp = NN.getP();
                size_t v = NN.getV();
                size_t nnp_a = NN.get();

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

                if (xp.distance2(xq) < r_cut2)
                    this->addPart(i,nnp_a);

                // Next particle
                ++NN;
            }

            ++it;
        }
    }

public:

    void Initialize(CellListImpl & cli, size_t pp, T r_cut, const openfpm::vector<Point<dim,T>> & pos, const openfpm::vector<struct pos_v<dim,T>> & pos2, size_t g_m, size_t opt = VL_NON_SYMMETRIC)
    {
//      this->cl_n.resize(g_m);
//      this->cl_base.resize(g_m*this->slot);

        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,pp,r_cut,g_m,cli,opt);
        }
        else
        {
            std::cerr << __FILE__ << ":" << __LINE__ << " error iterator with radius is not implemented yet " << std::endl;
        }
    }

    inline size_t getP(size_t part, size_t ele) const
    {
        return VerletBase::get(part,ele) & mask_low::sig_bits_fast;
    }

    inline size_t getV(size_t part, size_t ele) const
    {
        return (VerletBase::get(part,ele)) >> (sizeof(size_t)*8-sh_byte);
    }

    template<unsigned int impl=NO_CHECK> inline VerletNNIteratorM<dim,VerletListM<dim,T,sh_byte,CellListImpl,transform,VerletBase>,sh_byte> getNNIterator(size_t part_id)
    {
        VerletNNIteratorM<dim,VerletListM<dim,T,sh_byte,CellListImpl,transform,VerletBase>,sh_byte> vln(part_id,*this);

        return vln;
    }
};


#endif /* OPENFPM_DATA_SRC_NN_VERLETLIST_VERLETLISTM_HPP_ */