vector_dist.hpp

/*
 * Vector.hpp
 *
 *  Created on: Mar 5, 2015
 *      Author: Pietro Incardona
 */

#ifndef VECTOR_HPP_
#define VECTOR_HPP_

#include "HDF5_wr/HDF5_wr.hpp"
#include "VCluster/VCluster.hpp"
#include "Space/Shape/Point.hpp"
#include "Vector/Iterators/vector_dist_iterator.hpp"
#include "Space/Shape/Box.hpp"
#include "Vector/vector_dist_key.hpp"
#include "memory/PtrMemory.hpp"
#include "NN/CellList/CellList.hpp"
#include "NN/CellList/CellListFast_gen.hpp"
#include "util/common.hpp"
#include "util/object_util.hpp"
#include "memory/ExtPreAlloc.hpp"
#include "CSVWriter/CSVWriter.hpp"
#include "VTKWriter/VTKWriter.hpp"
#include "Decomposition/common.hpp"
#include "Grid/Iterators/grid_dist_id_iterator_dec.hpp"
#include "Grid/grid_key_dx_iterator_hilbert.hpp"
#include "Vector/vector_dist_ofb.hpp"
#include "Decomposition/CartDecomposition.hpp"
#include "data_type/aggregate.hpp"
#include "NN/VerletList/VerletList.hpp"
#include "vector_dist_comm.hpp"
#include "DLB/LB_Model.hpp"
#include "Vector/vector_map_iterator.hpp"
#include "NN/CellList/ParticleIt_Cells.hpp"
#include "NN/CellList/ProcKeys.hpp"

#define DEC_GRAN(gr) ((size_t)gr << 32)

#define VECTOR_DIST_ERROR_OBJECT std::runtime_error("Runtime vector distributed error");

#ifdef SE_CLASS3
#include "se_class3_vector.hpp"
#endif

#ifdef SE_CLASS3
    #define SE_CLASS3_VDIST_CONSTRUCTOR ,se3(getDecomposition(),*this)
#else
    #define SE_CLASS3_VDIST_CONSTRUCTOR
#endif


#define NO_ID false
#define ID true

// Perform a ghost get or a ghost put
#define GET 1
#define PUT 2

// Write the particles with ghost
#define NO_GHOST 0
#define WITH_GHOST 2

#define GCL_NON_SYMMETRIC 0
#define GCL_SYMMETRIC 1
#define GCL_HILBERT 2


template<unsigned int dim, typename St, typename CellL, typename Vector, unsigned int impl>
struct gcl
{
    static inline CellL get(Vector & vd, const St & r_cut, const Ghost<dim,St> & g)
    {
        return vd.template getCellList<CellL>(r_cut);
    }
};

template<unsigned int dim, typename St, typename CellL, typename Vector>
struct gcl<dim,St,CellL,Vector,GCL_HILBERT>
{
    static inline CellL get(Vector & vd, const St & r_cut, const Ghost<dim,St> & g)
    {
        return vd.getCellList_hilb(r_cut,g);
    }
};

template<unsigned int dim, typename St, typename CellL, typename Vector>
struct gcl<dim,St,CellL,Vector,GCL_SYMMETRIC>
{
    static inline CellL get(Vector & vd, const St & r_cut, const Ghost<dim,St> & g)
    {
        return vd.getCellListSym(r_cut);
    }
};


template<unsigned int dim, typename St, typename CellL, typename Vector, unsigned int impl>
struct gcl_An
{
    static inline CellL get(Vector & vd, const size_t (& div)[dim], const size_t (& pad)[dim], const Ghost<dim,St> & g)
    {
        return vd.template getCellListSym<CellL>(div,pad);
    }
};

#define CELL_MEMFAST(dim,St) CellList_gen<dim, St, Process_keys_lin, Mem_fast<>, shift<dim, St> >
#define CELL_MEMBAL(dim,St) CellList_gen<dim, St, Process_keys_lin, Mem_bal<>, shift<dim, St> >
#define CELL_MEMMW(dim,St) CellList_gen<dim, St, Process_keys_lin, Mem_mw<>, shift<dim, St> >

#define CELL_MEMFAST_HILB(dim,St) CellList_gen<dim, St, Process_keys_hilb, Mem_fast<>, shift<dim, St> >
#define CELL_MEMBAL_HILB(dim,St) CellList_gen<dim, St, Process_keys_hilb, Mem_bal<>, shift<dim, St> >
#define CELL_MEMMW_HILB(dim,St) CellList_gen<dim, St, Process_keys_hilb, Mem_mw<>, shift<dim, St> >

#define VERLET_MEMFAST(dim,St) VerletList<dim,St,Mem_fast<>,shift<dim,St> >
#define VERLET_MEMBAL(dim,St)  VerletList<dim,St,Mem_bal<>,shift<dim,St> >
#define VERLET_MEMMW(dim,St)   VerletList<dim,St,Mem_mw<>,shift<dim,St> >

#define VERLET_MEMFAST_INT(dim,St) VerletList<dim,St,Mem_fast<unsigned int>,shift<dim,St> >
#define VERLET_MEMBAL_INT(dim,St)  VerletList<dim,St,Mem_bal<unsigned int>,shift<dim,St> >
#define VERLET_MEMMW_INT(dim,St)   VerletList<dim,St,Mem_mw<unsigned int>,shift<dim,St> >

enum reorder_opt
{
    NO_REORDER = 0,
    HILBERT = 1,
    LINEAR = 2
};

template<unsigned int dim,
         typename St,
         typename prop,
         typename layout = typename memory_traits_lin<prop>::type,
         template <typename> class layout_base = memory_traits_lin,
         typename Decomposition = CartDecomposition<dim,St>,
         typename Memory = HeapMemory>
class vector_dist : public vector_dist_comm<dim,St,prop,layout,layout_base,Decomposition,Memory>
{
public:

    typedef vector_dist<dim,St,prop,layout,layout_base,Decomposition,Memory> self;

    typedef prop value_type;

private:

    size_t g_m = 0;

    openfpm::vector<Point<dim, St>,Memory> v_pos;

    openfpm::vector<prop,Memory,layout,layout_base> v_prp;

    Vcluster & v_cl;

    size_t opt = 0;

    openfpm::vector<std::string> prp_names;

#ifdef SE_CLASS3

    se_class3_vector<prop::max_prop,dim,St,Decomposition,self> se3;

#endif

    void init_structures(size_t np)
    {
        // convert to a local number of elements
        size_t p_np = np / v_cl.getProcessingUnits();

        // Get non divisible part
        size_t r = np % v_cl.getProcessingUnits();

        // Distribute the remain particles
        if (v_cl.getProcessUnitID() < r)
            p_np++;

        // resize the position vector
        v_pos.resize(p_np);

        // resize the properties vector
        v_prp.resize(p_np);

        g_m = p_np;
    }

    void check_parameters(Box<dim,St> & box)
    {
        // if the box is not valid return an error
        if (box.isValid() == false)
        {
            std::cerr << __FILE__ << ":" << __LINE__ << "  Error the domain is not valid " << box.toString() << std::endl;
            ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
        }

    }

    void check_ghost_compatible_rcut(St r_cut)
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (fabs(getDecomposition().getGhost().getLow(i)) < r_cut)
            {
                std::cerr << __FILE__ << ":" << __LINE__ << " Error the cut off radius " << r_cut << " is bigger that the ghost layer on the dimension " << i << " lower=" << getDecomposition().getGhost().getLow(i) << std::endl;
                ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
            }
        }
    }

    template<typename CellL, typename sfc_it>
    void reorder_sfc(openfpm::vector<Point<dim,St>> & v_pos_dest,
                         openfpm::vector<prop> & v_prp_dest,
                         sfc_it & h_it,
                         CellL & cell_list)
    {
        v_pos_dest.resize(v_pos.size());
        v_prp_dest.resize(v_prp.size());

        //Index for v_pos_dest
        size_t count = 0;

        grid_key_dx<dim> ksum;

        for (size_t i = 0; i < dim ; i++)
        {ksum.set_d(i,cell_list.getPadding(i));}

        while (h_it.isNext())
        {
          auto key = h_it.get();
          key += ksum;

          size_t lin = cell_list.getGrid().LinId(key);

          // for each particle in the Cell "lin"
          for (size_t i = 0; i < cell_list.getNelements(lin); i++)
          {
              //reorder
              auto v = cell_list.get(lin,i);
              v_pos_dest.get(count) = v_pos.get(v);
              v_prp_dest.get(count) = v_prp.get(v);

              count++;
          }
          ++h_it;
        }
    }

public:

    typedef St stype;

    static const unsigned int dims = dim;

    vector_dist<dim,St,prop,layout,layout_base,Decomposition,Memory> & operator=(const vector_dist<dim,St,prop,layout,layout_base,Decomposition,Memory> & v)
    {
        static_cast<vector_dist_comm<dim,St,prop,layout,layout_base,Decomposition,Memory> *>(this)->operator=(static_cast<vector_dist_comm<dim,St,prop,layout,layout_base,Decomposition,Memory>>(v));

        g_m = v.g_m;
        v_pos = v.v_pos;
        v_prp = v.v_prp;

#ifdef SE_CLASS3
        se3 = v.se3;
#endif

        opt = v.opt;

        return *this;
    }

    vector_dist<dim,St,prop,layout,layout_base,Decomposition,Memory> & operator=(vector_dist<dim,St,prop,layout,layout_base,Decomposition,Memory> && v)
    {
        static_cast<vector_dist_comm<dim,St,prop,layout,layout_base,Decomposition,Memory> *>(this)->operator=(static_cast<vector_dist_comm<dim,St,prop,layout,layout_base,Decomposition,Memory> >(v));

        g_m = v.g_m;
        v_pos.swap(v.v_pos);
        v_prp.swap(v.v_prp);

#ifdef SE_CLASS3
        se3 = v.se3;
#endif

        opt = v.opt;

        return *this;
    }


    vector_dist(const vector_dist<dim,St,prop,layout,layout_base,Decomposition,Memory> & v)
    :vector_dist_comm<dim,St,prop,layout,layout_base,Decomposition,Memory>(v.getDecomposition()),v_cl(v.v_cl) SE_CLASS3_VDIST_CONSTRUCTOR
    {
#ifdef SE_CLASS2
        check_new(this,8,VECTOR_DIST_EVENT,4);
#endif

        this->operator=(v);
    }

    vector_dist(vector_dist<dim,St,prop,layout,layout_base,Decomposition,Memory> && v) noexcept
    :v_cl(v.v_cl) SE_CLASS3_VDIST_CONSTRUCTOR
    {
#ifdef SE_CLASS2
        check_new(this,8,VECTOR_DIST_EVENT,4);
#endif

        this->operator=(v);

#ifdef SE_CLASS3
        se3.Initialize();
#endif
    }

    vector_dist(const Decomposition & dec, size_t np) :
    vector_dist_comm<dim,St,prop,layout,layout_base,Decomposition,Memory>(dec), v_cl(create_vcluster()) SE_CLASS3_VDIST_CONSTRUCTOR
    {
#ifdef SE_CLASS2
        check_new(this,8,VECTOR_DIST_EVENT,4);
#endif

        init_structures(np);

#ifdef SE_CLASS3
        se3.Initialize();
#endif
    }


    vector_dist(size_t np, Box<dim, St> box, const size_t (&bc)[dim], const Ghost<dim, St> & g, size_t opt = 0, const grid_sm<dim,void> & gdist = grid_sm<dim,void>())
    :v_cl(create_vcluster()),opt(opt) SE_CLASS3_VDIST_CONSTRUCTOR
    {
#ifdef SE_CLASS2
        check_new(this,8,VECTOR_DIST_EVENT,4);
#endif

        if (opt >> 32 != 0)
            this->setDecompositionGranularity(opt >> 32);

        check_parameters(box);

        init_structures(np);
        this->init_decomposition(box,bc,g,opt,gdist);

#ifdef SE_CLASS3
        se3.Initialize();
#endif
    }

    ~vector_dist()
    {
#ifdef SE_CLASS2
        check_delete(this);
#endif
    }

    void clear()
    {
        resize(0);
    }

    size_t size_local() const
    {
        return g_m;
    }

    size_t size_local_with_ghost() const
    {
        return v_pos.size();
    }

#ifndef ONLY_READWRITE_GETTER

    inline auto getPos(vect_dist_key_dx vec_key) -> decltype(v_pos.template get<0>(vec_key.getKey()))
    {
#ifdef SE_CLASS3
        check_for_pos_nan_inf<prop::max_prop_real,prop::max_prop>(*this,vec_key.getKey());
#endif

        return v_pos.template get<0>(vec_key.getKey());
    }

    inline auto getPos(vect_dist_key_dx vec_key) const -> decltype(v_pos.template get<0>(vec_key.getKey()))
    {
#ifdef SE_CLASS3
        check_for_pos_nan_inf<prop::max_prop_real,prop::max_prop>(*this,vec_key.getKey());
#endif
        return v_pos.template get<0>(vec_key.getKey());
    }

    inline auto getPos(size_t vec_key) -> decltype(v_pos.template get<0>(vec_key))
    {
#ifdef SE_CLASS3
        check_for_pos_nan_inf<prop::max_prop_real,prop::max_prop>(*this,vec_key);
#endif
        return v_pos.template get<0>(vec_key);
    }

    inline auto getPos(size_t vec_key) const -> decltype(v_pos.template get<0>(vec_key))
    {
#ifdef SE_CLASS3
        check_for_pos_nan_inf<prop::max_prop_real,prop::max_prop>(*this,vec_key);
#endif
        return v_pos.template get<0>(vec_key);
    }

    template<unsigned int id> inline auto getProp(vect_dist_key_dx vec_key) -> decltype(v_prp.template get<id>(vec_key.getKey()))
    {
#ifdef SE_CLASS3
        check_for_prop_nan_inf<id,prop::max_prop+SE3_STATUS>(*this,vec_key.getKey());
#endif
        return v_prp.template get<id>(vec_key.getKey());
    }

    template<unsigned int id> inline auto getProp(vect_dist_key_dx vec_key) const -> decltype(v_prp.template get<id>(vec_key.getKey()))
    {
#ifdef SE_CLASS3
        check_for_prop_nan_inf<id,prop::max_prop+SE3_STATUS>(*this,vec_key.getKey());
#endif
        return v_prp.template get<id>(vec_key.getKey());
    }

    template<unsigned int id> inline auto getProp(size_t vec_key) -> decltype(v_prp.template get<id>(vec_key))
    {
#ifdef SE_CLASS3
        check_for_prop_nan_inf<id,prop::max_prop+SE3_STATUS>(*this,vec_key);
#endif
        return v_prp.template get<id>(vec_key);
    }

    template<unsigned int id> inline auto getProp(size_t vec_key) const -> decltype(v_prp.template get<id>(vec_key))
    {
#ifdef SE_CLASS3
        check_for_prop_nan_inf<id,prop::max_prop+SE3_STATUS>(*this,vec_key);
#endif
        return v_prp.template get<id>(vec_key);
    }

#endif


    inline auto getPosNC(vect_dist_key_dx vec_key) -> decltype(v_pos.template get<0>(vec_key.getKey()))
    {
        return v_pos.template get<0>(vec_key.getKey());
    }

    inline auto getPosNC(vect_dist_key_dx vec_key) const -> decltype(v_pos.template get<0>(vec_key.getKey()))
    {
        return v_pos.template get<0>(vec_key.getKey());
    }

    inline auto getPosNC(size_t vec_key) -> decltype(v_pos.template get<0>(vec_key))
    {
        return v_pos.template get<0>(vec_key);
    }

    inline auto getPosNC(size_t vec_key) const -> decltype(v_pos.template get<0>(vec_key))
    {
        return v_pos.template get<0>(vec_key);
    }

    template<unsigned int id> inline auto getPropNC(vect_dist_key_dx vec_key) -> decltype(v_prp.template get<id>(vec_key.getKey()))
    {
        return v_prp.template get<id>(vec_key.getKey());
    }

    template<unsigned int id> inline auto getPropNC(vect_dist_key_dx vec_key) const -> decltype(v_prp.template get<id>(vec_key.getKey()))
    {
        return v_prp.template get<id>(vec_key.getKey());
    }

    template<unsigned int id> inline auto getPropNC(size_t vec_key) -> decltype(v_prp.template get<id>(vec_key))
    {
        return v_prp.template get<id>(vec_key);
    }

    template<unsigned int id> inline auto getPropNC(size_t vec_key) const -> decltype(v_prp.template get<id>(vec_key))
    {
        return v_prp.template get<id>(vec_key);
    }


    inline auto getPosWrite(vect_dist_key_dx vec_key) -> decltype(v_pos.template get<0>(vec_key.getKey()))
    {
#ifdef SE_CLASS3
        se3.template write<prop::max_prop_real>(*this,vec_key.getKey());
#endif

        return v_pos.template get<0>(vec_key.getKey());
    }

    inline auto getPosRead(vect_dist_key_dx vec_key) const -> decltype(v_pos.template get<0>(vec_key.getKey()))
    {
#ifdef SE_CLASS3
        se3.template read<prop::max_prop_real>(*this,vec_key.getKey());
#endif

        return v_pos.template get<0>(vec_key.getKey());
    }

    template<unsigned int id> inline auto getPropWrite(vect_dist_key_dx vec_key) -> decltype(v_prp.template get<id>(vec_key.getKey()))
    {
#ifdef SE_CLASS3
        se3.template write<id>(*this,vec_key.getKey());
#endif

        return v_prp.template get<id>(vec_key.getKey());
    }

    template<unsigned int id> inline auto getPropRead(vect_dist_key_dx vec_key) const -> decltype(v_prp.template get<id>(vec_key.getKey()))
    {
#ifdef SE_CLASS3
        se3.template read<id>(*this,vec_key.getKey());
#endif

        return v_prp.template get<id>(vec_key.getKey());
    }


    void add()
    {
        v_prp.insert(g_m);
        v_pos.insert(g_m);

        g_m++;

#ifdef SE_CLASS3
        for (size_t i = 0 ; i < prop::max_prop_real+1 ; i++)
            v_prp.template get<prop::max_prop_real>(g_m-1)[i] = UNINITIALIZED;
#endif
    }

#ifndef ONLY_READWRITE_GETTER

    inline auto getLastPos() -> decltype(v_pos.template get<0>(0))
    {
        return v_pos.template get<0>(g_m - 1);
    }

    template<unsigned int id> inline auto getLastProp() -> decltype(v_prp.template get<id>(0))
    {
        return v_prp.template get<id>(g_m - 1);
    }

#endif


    inline auto getLastPosRead() -> decltype(v_pos.template get<0>(0))
    {
#ifdef SE_CLASS3
        se3.template read<prop::max_prop_real>(*this,g_m-1);
#endif

        return v_pos.template get<0>(g_m - 1);
    }

    template<unsigned int id> inline auto getLastPropRead() -> decltype(v_prp.template get<id>(0))
    {
#ifdef SE_CLASS3
        se3.read<id>(*this,g_m-1);
#endif

        return v_prp.template get<id>(g_m - 1);
    }


    inline auto getLastPosWrite() -> decltype(v_pos.template get<0>(0))
    {
#ifdef SE_CLASS3
        se3.template write<prop::max_prop_real>(*this,g_m-1);
#endif

        return v_pos.template get<0>(g_m - 1);
    }

    template<unsigned int id> inline auto getLastPropWrite() -> decltype(v_prp.template get<id>(0))
    {
#ifdef SE_CLASS3
        se3.template write<id>(*this,g_m-1);
#endif

        return v_prp.template get<id>(g_m - 1);
    }


    template<typename CellL = CellList<dim, St, Mem_fast<>, shift<dim, St> > > CellL getCellListSym(St r_cut)
    {
#ifdef SE_CLASS1
        if (!(opt & BIND_DEC_TO_GHOST))
        {
            if (getDecomposition().getGhost().getLow(dim-1) == 0.0)
            {
                std::cerr << __FILE__ << ":" << __LINE__ << " Error the vector has been constructed without BIND_DEC_TO_GHOST, If you construct a vector without BIND_DEC_TO_GHOST the ghost must be full without reductions " << std::endl;
                ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
            }
        }
#endif

        // Cell list
        CellL cell_list;

        size_t pad = 0;
        CellDecomposer_sm<dim,St,shift<dim,St>> cd_sm;
        cl_param_calculateSym(getDecomposition().getDomain(),cd_sm,getDecomposition().getGhost(),r_cut,pad);

        // Processor bounding box
        Box<dim, St> pbox = getDecomposition().getProcessorBounds();

        // Ghost padding extension
        Ghost<dim,size_t> g_ext(0);
        cell_list.Initialize(cd_sm,pbox,pad);
        cell_list.set_ndec(getDecomposition().get_ndec());

        updateCellListSym(cell_list);

        return cell_list;
    }

    template<typename CellL = CellList<dim, St, Mem_fast<>, shift<dim, St> > >
    CellL getCellListSym(const size_t (& div)[dim],
                         const size_t (& pad)[dim])
    {
#ifdef SE_CLASS1
        if (!(opt & BIND_DEC_TO_GHOST))
        {
            if (getDecomposition().getGhost().getLow(dim-1) == 0.0)
            {
                std::cerr << __FILE__ << ":" << __LINE__ << " Error the vector has been constructed without BIND_DEC_TO_GHOST, If you construct a vector without BIND_DEC_TO_GHOST the ghost must be full without reductions " << std::endl;
                ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
            }
        }
#endif

        size_t pad_max = pad[0];
        for (size_t i = 1 ; i < dim ; i++)
        {if (pad[i] > pad_max)  {pad_max = pad[i];}}

        // Cell list
        CellL cell_list;

        CellDecomposer_sm<dim,St,shift<dim,St>> cd_sm;
        cd_sm.setDimensions(getDecomposition().getDomain(),div,pad_max);

        // Processor bounding box
        Box<dim, St> pbox = getDecomposition().getProcessorBounds();

        // Ghost padding extension
        Ghost<dim,size_t> g_ext(0);
        cell_list.Initialize(cd_sm,pbox,pad_max);
        cell_list.set_ndec(getDecomposition().get_ndec());

        updateCellListSym(cell_list);

        return cell_list;
    }

    template<typename CellL = CellList_gen<dim, St, Process_keys_lin, Mem_fast<>, shift<dim, St> > >
    CellL getCellList(St r_cut, bool no_se3 = false)
    {
#ifdef SE_CLASS3
        if (no_se3 == false)
        {se3.getNN();}
#endif
#ifdef SE_CLASS1
        check_ghost_compatible_rcut(r_cut);
#endif

        // Get ghost and anlarge by 1%
        Ghost<dim,St> g = getDecomposition().getGhost();
        g.magnify(1.013);

        return getCellList<CellL>(r_cut, g,no_se3);
    }

    template<typename CellL = CellList_gen<dim, St, Process_keys_hilb, Mem_fast<>, shift<dim, St> > >
    CellL getCellList_hilb(St r_cut)
    {
#ifdef SE_CLASS3
        se3.getNN();
#endif
#ifdef SE_CLASS1
        check_ghost_compatible_rcut(r_cut);
#endif

        // Get ghost and anlarge by 1%
        Ghost<dim,St> g = getDecomposition().getGhost();
        g.magnify(1.013);

        return getCellList_hilb(r_cut, g);
    }

    template<typename CellL> void updateCellList(CellL & cell_list, bool no_se3 = false)
    {
#ifdef SE_CLASS3
        if (no_se3 == false)
        {se3.getNN();}
#endif

        // This function assume equal spacing in all directions
        // but in the worst case we take the maximum
        St r_cut = 0;
        for (size_t i = 0 ; i < dim ; i++)
        {r_cut = std::max(r_cut,cell_list.getCellBox().getHigh(i));}

        // Here we have to check that the Cell-list has been constructed
        // from the same decomposition
        bool to_reconstruct = cell_list.get_ndec() != getDecomposition().get_ndec();

        if (to_reconstruct == false)
        {
            populate_cell_list(v_pos,cell_list,g_m,CL_NON_SYMMETRIC);

            cell_list.set_gm(g_m);
        }
        else
        {
            CellL cli_tmp = gcl<dim,St,CellL,self,GCL_NON_SYMMETRIC>::get(*this,r_cut,getDecomposition().getGhost());

            cell_list.swap(cli_tmp);
        }
    }

    template<typename CellL = CellList<dim, St, Mem_fast<>, shift<dim, St> > >
    void updateCellListSym(CellL & cell_list)
    {
#ifdef SE_CLASS3
        se3.getNN();
#endif

        // Here we have to check that the Cell-list has been constructed
        // from the same decomposition
        bool to_reconstruct = cell_list.get_ndec() != getDecomposition().get_ndec();

        if (to_reconstruct == false)
        {
            populate_cell_list(v_pos,cell_list,g_m,CL_SYMMETRIC);

            cell_list.set_gm(g_m);
        }
        else
        {
            CellL cli_tmp = gcl_An<dim,St,CellL,self,GCL_SYMMETRIC>::get(*this,
                                                                         cell_list.getDivWP(),
                                                                         cell_list.getPadding(),
                                                                         getDecomposition().getGhost());

            cell_list.swap(cli_tmp);
        }
    }

    template<typename CellL = CellList_gen<dim, St, Process_keys_lin, Mem_fast<>, shift<dim, St> > >
    CellL getCellList(St r_cut, const Ghost<dim, St> & enlarge, bool no_se3 = false)
    {
#ifdef SE_CLASS3
        if (no_se3 == false)
        {se3.getNN();}
#endif

        CellL cell_list;

        // Division array
        size_t div[dim];

        // get the processor bounding box
        Box<dim, St> pbox = getDecomposition().getProcessorBounds();

        // Processor bounding box
        cl_param_calculate(pbox, div, r_cut, enlarge);

        cell_list.Initialize(pbox, div);
        cell_list.set_gm(g_m);
        cell_list.set_ndec(getDecomposition().get_ndec());

        updateCellList(cell_list,no_se3);

        return cell_list;
    }

    template<typename CellL = CellList_gen<dim, St, Process_keys_hilb, Mem_fast<>, shift<dim, St> > > CellL getCellList_hilb(St r_cut, const Ghost<dim, St> & enlarge)
    {
#ifdef SE_CLASS3
        se3.getNN();
#endif

        CellL cell_list;

        // Division array
        size_t div[dim];

        // get the processor bounding box
        Box<dim, St> pbox = getDecomposition().getProcessorBounds();

        // Processor bounding box
        cl_param_calculate(pbox,div, r_cut, enlarge);

        cell_list.Initialize(pbox, div);
        cell_list.set_gm(g_m);
        cell_list.set_ndec(getDecomposition().get_ndec());

        updateCellList(cell_list);

        return cell_list;
    }

    template <typename VerletL = VerletList<dim,St,Mem_fast<>,shift<dim,St> >>
    VerletL getVerletSym(St r_cut)
    {
#ifdef SE_CLASS3
        se3.getNN();
#endif

        VerletL ver;

        // Processor bounding box
        Box<dim, St> pbox = getDecomposition().getProcessorBounds();

        ver.InitializeSym(getDecomposition().getDomain(),pbox,getDecomposition().getGhost(),r_cut,v_pos,g_m);

        ver.set_ndec(getDecomposition().get_ndec());

        return ver;
    }

    template <typename VerletL = VerletList<dim,St,Mem_fast<>,shift<dim,St> >>
    VerletL getVerletCrs(St r_cut)
    {
#ifdef SE_CLASS1
        if (!(opt & BIND_DEC_TO_GHOST))
        {
            std::cerr << __FILE__ << ":" << __LINE__ << " Error the vector has been constructed without BIND_DEC_TO_GHOST, getVerletCrs require the vector to be constructed with BIND_DEC_TO_GHOST option " << std::endl;
            ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
        }
#endif

#ifdef SE_CLASS3
        se3.getNN();
#endif

        VerletL ver;

        // Processor bounding box
        Box<dim, St> pbox = getDecomposition().getProcessorBounds();

        // Initialize the verlet list
        ver.InitializeCrs(getDecomposition().getDomain(),pbox,getDecomposition().getGhost(),r_cut,v_pos,g_m);

        // Get the internal cell list
        auto & NN = ver.getInternalCellList();

        // Shift
        grid_key_dx<dim> shift;

        // Add padding
        for (size_t i = 0 ; i < dim ; i++)
            shift.set_d(i,NN.getPadding(i));

        grid_sm<dim,void> gs = NN.getInternalGrid();

        getDecomposition().setNNParameters(shift,gs);

        ver.createVerletCrs(r_cut,g_m,v_pos,
                            getDecomposition().getCRSDomainCells(),
                            getDecomposition().getCRSAnomDomainCells());

        ver.set_ndec(getDecomposition().get_ndec());

        return ver;
    }

    template <typename VerletL = VerletList<dim,St,Mem_fast<>,shift<dim,St> >>
    VerletL getVerlet(St r_cut)
    {
#ifdef SE_CLASS3
        se3.getNN();
#endif

        VerletL ver;

        // get the processor bounding box
        Box<dim, St> bt = getDecomposition().getProcessorBounds();

        // Get the ghost
        Ghost<dim,St> g = getDecomposition().getGhost();
        g.magnify(1.013);

        // enlarge the box where the Verlet is defined
        bt.enlarge(g);

        ver.Initialize(bt,getDecomposition().getProcessorBounds(),r_cut,v_pos,g_m,VL_NON_SYMMETRIC);

        ver.set_ndec(getDecomposition().get_ndec());

        return ver;
    }

    template<typename Mem_type> void updateVerlet(VerletList<dim,St,Mem_type,shift<dim,St> > & ver, St r_cut, size_t opt = VL_NON_SYMMETRIC)
    {
#ifdef SE_CLASS3
        se3.getNN();
#endif

        if (opt == VL_SYMMETRIC)
        {
            auto & NN = ver.getInternalCellList();

            // Here we have to check that the Box defined by the Cell-list is the same as the domain box of this
            // processor. if it is not like that we have to completely reconstruct from stratch
            bool to_reconstruct = NN.get_ndec() != getDecomposition().get_ndec();

            if (to_reconstruct == false)
                ver.update(getDecomposition().getDomain(),r_cut,v_pos,g_m, opt);
            else
            {
                VerletList<dim,St,Mem_type,shift<dim,St> > ver_tmp;

                ver_tmp = getVerlet<VerletList<dim,St,Mem_type,shift<dim,St> >>(r_cut);
                ver.swap(ver);
            }
        }
        else if (opt == VL_CRS_SYMMETRIC)
        {
#ifdef SE_CLASS1
            if ((opt & BIND_DEC_TO_GHOST))
            {
                std::cerr << __FILE__ << ":" << __LINE__ << " Error the vector has been constructed without BIND_DEC_TO_GHOST, updateVerlet with the option VL_CRS_SYMMETRIC require the vector to be constructed with BIND_DEC_TO_GHOST option " << std::endl;
                ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
            }
#endif

            auto & NN = ver.getInternalCellList();

            // Here we have to check that the Box defined by the Cell-list is the same as the domain box of this
            // processor. if it is not like that we have to completely reconstruct from stratch
            bool to_reconstruct = NN.get_ndec() != getDecomposition().get_ndec();

            if (to_reconstruct == false)
            {
                // Shift
                grid_key_dx<dim> shift;

                // Add padding
                for (size_t i = 0 ; i < dim ; i++)
                    shift.set_d(i,NN.getPadding(i));

                grid_sm<dim,void> gs = NN.getInternalGrid();

                getDecomposition().setNNParameters(shift,gs);

                ver.updateCrs(getDecomposition().getDomain(),r_cut,v_pos,g_m,
                              getDecomposition().getCRSDomainCells(),
                              getDecomposition().getCRSAnomDomainCells());
            }
            else
            {
                VerletList<dim,St,Mem_type,shift<dim,St> > ver_tmp;

                ver_tmp = getVerletCrs<VerletList<dim,St,Mem_type,shift<dim,St> >>(r_cut);
                ver.swap(ver_tmp);
            }
        }
        else
        {
            auto & NN = ver.getInternalCellList();

            // Here we have to check that the Box defined by the Cell-list is the same as the domain box of this
            // processor. if it is not like that we have to completely reconstruct from stratch
            bool to_reconstruct = NN.get_ndec() != getDecomposition().get_ndec();

            if (to_reconstruct == false)
                ver.update(getDecomposition().getDomain(),r_cut,v_pos,g_m, opt);
            else
            {
                VerletList<dim,St,Mem_type,shift<dim,St> > ver_tmp;

                ver_tmp = getVerlet<VerletList<dim,St,Mem_type,shift<dim,St> >>(r_cut);
                ver.swap(ver_tmp);
            }
        }
    }


    template<typename CellL=CellList_gen<dim,St,Process_keys_lin,Mem_bal<>,shift<dim,St> > >
    void reorder (int32_t m, reorder_opt opt = reorder_opt::HILBERT)
    {
        reorder<CellL>(m,getDecomposition().getGhost(),opt);
    }


    template<typename CellL=CellList_gen<dim,St,Process_keys_lin,Mem_bal<>,shift<dim,St> > >
    void reorder(int32_t m, const Ghost<dim,St> & enlarge, reorder_opt opt = reorder_opt::HILBERT)
    {
        // reset the ghost part
        v_pos.resize(g_m);
        v_prp.resize(g_m);


        CellL cell_list;

        // calculate the parameters of the cell list

        // get the processor bounding box
        Box<dim,St> pbox = getDecomposition().getProcessorBounds();
        // extend by the ghost
        pbox.enlarge(enlarge);

        size_t div[dim];

        // Calculate the division array and the cell box
        for (size_t i = 0 ; i < dim ; i++)
        {
            div[i] = 1 << m;
        }

        cell_list.Initialize(pbox,div);
        cell_list.set_gm(g_m);

        // for each particle add the particle to the cell list

        auto it = getIterator();

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

            Point<dim,St> xp = this->getPos(key);

            cell_list.add(xp,key.getKey());

            ++it;
        }

        // Use cell_list to reorder v_pos

        //destination vector
        openfpm::vector<Point<dim,St>> v_pos_dest;
        openfpm::vector<prop> v_prp_dest;

        if (opt == reorder_opt::HILBERT)
        {
            grid_key_dx_iterator_hilbert<dim> h_it(m);

            reorder_sfc<CellL,grid_key_dx_iterator_hilbert<dim>>(v_pos_dest,v_prp_dest,h_it,cell_list);
        }
        else if (opt == reorder_opt::LINEAR)
        {
            grid_sm<dim,void> gs(div);
            grid_key_dx_iterator<dim> h_it(gs);

            reorder_sfc<CellL,grid_key_dx_iterator<dim>>(v_pos_dest,v_prp_dest,h_it,cell_list);
        }
        else
        {
            // We do nothing, we second swap nullify the first
            v_pos.swap(v_pos_dest);
            v_prp.swap(v_prp_dest);
        }

        v_pos.swap(v_pos_dest);
        v_prp.swap(v_prp_dest);
    }

    size_t init_size_accum(size_t np)
    {
        size_t accum = 0;

        // convert to a local number of elements
        size_t p_np = np / v_cl.getProcessingUnits();

        // Get non divisible part
        size_t r = np % v_cl.getProcessingUnits();

        accum = p_np * v_cl.getProcessUnitID();

        // Distribute the remain particles
        if (v_cl.getProcessUnitID() <= r)
            accum += v_cl.getProcessUnitID();
        else
            accum += r;

        return accum;
    }

    vector_dist_iterator getIterator()
    {
#ifdef SE_CLASS3
        se3.getIterator();
#endif
        return vector_dist_iterator(0, v_pos.size());
    }

    inline grid_dist_id_iterator_dec<Decomposition> getGridIterator(const size_t (&sz)[dim])
    {
        grid_key_dx<dim> start;
        grid_key_dx<dim> stop;
        for (size_t i = 0; i < dim; i++)
        {
            start.set_d(i, 0);
            stop.set_d(i, sz[i] - 1);
        }

        grid_dist_id_iterator_dec<Decomposition> it_dec(getDecomposition(), sz, start, stop);
        return it_dec;
    }

    vector_dist_iterator getGhostIterator() const
    {
#ifdef SE_CLASS3
        se3.getIterator();
#endif

        return vector_dist_iterator(g_m, v_pos.size());
    }

    vector_dist_iterator getGhostIterator_no_se3() const
    {
        return vector_dist_iterator(g_m, v_pos.size());
    }

    template<typename CellList> ParticleIt_Cells<dim,CellList>
    getDomainIteratorCells(CellList & NN)
    {
#ifdef SE_CLASS3
        se3.getIterator();
#endif

        // Shift
        grid_key_dx<dim> shift;

        // Add padding
        for (size_t i = 0 ; i < dim ; i++)
            shift.set_d(i,NN.getPadding(i));

        grid_sm<dim,void> gs = NN.getInternalGrid();

        getDecomposition().setNNParameters(shift,gs);

        return ParticleIt_Cells<dim,CellList>(NN,getDecomposition().getDomainCells(),g_m);
    }

    vector_dist_iterator getDomainIterator() const
    {
#ifdef SE_CLASS3
        se3.getIterator();
#endif

        return vector_dist_iterator(0, g_m);
    }

    vector_dist_iterator getDomainIterator_no_se3() const
    {
        return vector_dist_iterator(0, g_m);
    }

    vector_dist_iterator getDomainAndGhostIterator() const
    {
#ifdef SE_CLASS3
        se3.getIterator();
#endif

        return vector_dist_iterator(0, v_pos.size());
    }

    vector_dist_iterator getDomainAndGhostIterator_no_se3() const
    {
        return vector_dist_iterator(0, v_pos.size());
    }

    inline Decomposition & getDecomposition()
    {
        return vector_dist_comm<dim,St,prop,layout,layout_base,Decomposition,Memory>::getDecomposition();
    }

    inline const Decomposition & getDecomposition() const
    {
        return vector_dist_comm<dim,St,prop,layout,layout_base,Decomposition,Memory>::getDecomposition();
    }

    template<unsigned int ... prp> void map_list(size_t opt = NONE)
    {
#ifdef SE_CLASS3
        se3.map_pre();
#endif

        this->template map_list_<prp...>(v_pos,v_prp,g_m,opt);

#ifdef SE_CLASS3
        se3.map_post();
#endif
    }


    template<typename obp = KillParticle> void map(size_t opt = NONE)
    {
#ifdef SE_CLASS3
        se3.map_pre();
#endif

        this->template map_<obp>(v_pos,v_prp,g_m,opt);

#ifdef SE_CLASS3
        se3.map_post();
#endif
    }

    template<int ... prp> inline void ghost_get(size_t opt = WITH_POSITION)
    {
#ifdef SE_CLASS1
        if (getDecomposition().getProcessorBounds().isValid() == false && size_local() != 0)
        {
            std::cerr << __FILE__ << ":" << __LINE__ << " Error the processor " << v_cl.getProcessUnitID() << " has particles, but is supposed to be unloaded" << std::endl;
            ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
        }
#endif

#ifdef SE_CLASS3
        se3.template ghost_get_pre<prp...>(opt);
#endif

        this->template ghost_get_<prp...>(v_pos,v_prp,g_m,opt);

#ifdef SE_CLASS3

        this->template ghost_get_<prop::max_prop_real>(v_pos,v_prp,g_m,opt | KEEP_PROPERTIES);

        se3.template ghost_get_post<prp...>(opt);
#endif
    }

    template<template<typename,typename> class op, int ... prp> inline void ghost_put(size_t opt_ = NONE)
    {
#ifdef SE_CLASS3
        se3.template ghost_put<prp...>();
#endif
        this->template ghost_put_<op,prp...>(v_pos,v_prp,g_m,opt_);
    }

    void remove(openfpm::vector<size_t> & keys, size_t start = 0)
    {
        v_pos.remove(keys, start);
        v_prp.remove(keys, start);

        g_m -= keys.size();
    }

    void remove(size_t key)
    {
        v_pos.remove(key);
        v_prp.remove(key);

        g_m--;
    }

    template <typename Model=ModelLin>inline void addComputationCosts(const self & vd, Model md=Model())
    {
        CellDecomposer_sm<dim, St, shift<dim,St>> cdsm;

        Decomposition & dec = getDecomposition();

        cdsm.setDimensions(dec.getDomain(), dec.getDistGrid().getSize(), 0);

        auto it = vd.getDomainIterator();

        while (it.isNext())
        {
            size_t v = cdsm.getCell(vd.getPos(it.get()));

            md.addComputation(dec,vd,v,it.get().getKey());

            ++it;
        }
    }

    template <typename Model=ModelLin> void finalizeComputationCosts(Model md=Model(), size_t ts = 1)
    {
        Decomposition & dec = getDecomposition();
        auto & dist = getDecomposition().getDistribution();

        dec.computeCommunicationAndMigrationCosts(ts);

        // Go throught all the sub-sub-domains and apply the model

        for (size_t i = 0 ; i < dist.getNOwnerSubSubDomains(); i++)
        {md.applyModel(dec,dist.getOwnerSubSubDomain(i));}

        dist.setDistTol(md.distributionTol());
    }

    void initializeComputationCosts()
    {
        Decomposition & dec = getDecomposition();
        auto & dist = getDecomposition().getDistribution();

        for (size_t i = 0; i < dist.getNOwnerSubSubDomains() ; i++)
        {dec.setSubSubDomainComputationCost(dist.getOwnerSubSubDomain(i) , 1);}
    }

    template <typename Model=ModelLin>inline void addComputationCosts(Model md=Model(), size_t ts = 1)
    {
        initializeComputationCosts();

        addComputationCosts(*this,md);

        finalizeComputationCosts(md,ts);
    }

    inline void save(const std::string & filename) const
    {
        HDF5_writer<VECTOR_DIST> h5s;

        h5s.save(filename,v_pos,v_prp);
    }

    inline void load(const std::string & filename)
    {
        HDF5_reader<VECTOR_DIST> h5l;

        h5l.load(filename,v_pos,v_prp,g_m);
    }

    inline bool write(std::string out, int opt = VTK_WRITER)
    {

        if ((opt & 0x0FFF0000) == CSV_WRITER)
        {
            // CSVWriter test
            CSVWriter<openfpm::vector<Point<dim,St>>, openfpm::vector<prop> > csv_writer;

            std::string output = std::to_string(out + "_" + std::to_string(v_cl.getProcessUnitID()) + std::to_string(".csv"));

            // Write the CSV
            return csv_writer.write(output,v_pos,v_prp);
        }
        else
        {
            file_type ft = file_type::ASCII;

            if (opt & FORMAT_BINARY)
                ft = file_type::BINARY;

            // VTKWriter for a set of points
            VTKWriter<boost::mpl::pair<openfpm::vector<Point<dim,St>>, openfpm::vector<prop>>, VECTOR_POINTS> vtk_writer;
            vtk_writer.add(v_pos,v_prp,g_m);

            std::string output = std::to_string(out + "_" + std::to_string(v_cl.getProcessUnitID()) + std::to_string(".vtk"));

            // Write the VTK file
            return vtk_writer.write(output,prp_names,"particles",ft);
        }

        return false;
    }

    void deleteGhost()
    {
        v_pos.resize(g_m);
        v_prp.resize(g_m);
    }

    void resize(size_t rs)
    {
        deleteGhost();

        v_pos.resize(rs);
        v_prp.resize(rs);

        g_m = rs;
    }

    inline bool write_frame(std::string out, size_t iteration, int opt = VTK_WRITER)
    {
        if ((opt & 0x0FFF0000) == CSV_WRITER)
        {
            // CSVWriter test
            CSVWriter<openfpm::vector<Point<dim, St>>, openfpm::vector<prop> > csv_writer;

            std::string output = std::to_string(out + "_" + std::to_string(v_cl.getProcessUnitID()) + "_" + std::to_string(iteration) + std::to_string(".csv"));

            // Write the CSV
            return csv_writer.write(output, v_pos, v_prp);
        }
        else
        {
            file_type ft = file_type::ASCII;

            if (opt & FORMAT_BINARY)
                ft = file_type::BINARY;

            // VTKWriter for a set of points
            VTKWriter<boost::mpl::pair<openfpm::vector<Point<dim,St>>, openfpm::vector<prop>>, VECTOR_POINTS> vtk_writer;
            vtk_writer.add(v_pos,v_prp,g_m);

            std::string output = std::to_string(out + "_" + std::to_string(v_cl.getProcessUnitID()) + "_" + std::to_string(iteration) + std::to_string(".vtk"));

            // Write the VTK file
            return vtk_writer.write(output,prp_names,"particles",ft);
        }
    }

    void getCellListParams(St r_cut, size_t (&div)[dim],Box<dim, St> & box, Ghost<dim,St> enlarge = Ghost<dim,St>(0.0))
    {
        // get the processor bounding box
        Box<dim, St> pbox = getDecomposition().getProcessorBounds();

        // enlarge the processor bounding box by the ghost
        Ghost<dim,St> g = getDecomposition().getGhost();
        pbox.enlarge(g);

        cl_param_calculate(pbox, div,r_cut,enlarge);

        // output the fixed domain
        box = pbox;
    }

    long int who()
    {
#ifdef SE_CLASS2
        return check_whoami(this,8);
#else
        return -1;
#endif
    }

    Vcluster & getVC()
    {
#ifdef SE_CLASS2
        check_valid(this,8);
#endif
        return v_cl;
    }

    const openfpm::vector<Point<dim,St>> & getPosVector() const
    {
        return v_pos;
    }

    openfpm::vector<Point<dim,St>> & getPosVector()
    {
        return v_pos;
    }

    const openfpm::vector<prop> & getPropVector() const
    {
        return v_prp;
    }

    openfpm::vector<prop> & getPropVector()
    {
        return v_prp;
    }

    size_t accum()
    {
        openfpm::vector<size_t> accu;

        size_t sz = size_local();

        v_cl.allGather(sz,accu);
        v_cl.execute();

        sz = 0;

        for (size_t i = 0 ; i < v_cl.getProcessUnitID() ; i++)
        {sz += accu.get(i);}

        return sz;
    }

    template<typename cli> ParticleItCRS_Cells<dim,cli> getParticleIteratorCRS_Cell(cli & NN)
    {
#ifdef SE_CLASS3
        se3.getIterator();
#endif

#ifdef SE_CLASS1
        if (!(opt & BIND_DEC_TO_GHOST))
        {
            std::cerr << __FILE__ << ":" << __LINE__ << " Error the vector has been constructed without BIND_DEC_TO_GHOST, getParticleIteratorCRS_Cell require the vector to be constructed with BIND_DEC_TO_GHOST option " << std::endl;
            ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
        }
#endif

        // Shift
        grid_key_dx<dim> shift;

        // Add padding
        for (size_t i = 0 ; i < dim ; i++)
            shift.set_d(i,NN.getPadding(i));

        grid_sm<dim,void> gs = NN.getInternalGrid();

        getDecomposition().setNNParameters(shift,gs);

        // First we check that
        return ParticleItCRS_Cells<dim,cli>(NN,getDecomposition().getCRSDomainCells(),
                                               getDecomposition().getCRSAnomDomainCells(),
                                               NN.getNNc_sym());
    }

    void setPropNames(const openfpm::vector<std::string> & names)
    {
        prp_names = names;
    }

    template<typename vrl> openfpm::vector_key_iterator_seq<typename vrl::Mem_type_type::loc_index> getParticleIteratorCRS(vrl & NN)
    {
#ifdef SE_CLASS1
        if (!(opt & BIND_DEC_TO_GHOST))
        {
            std::cerr << __FILE__ << ":" << __LINE__ << " Error the vector has been constructed without BIND_DEC_TO_GHOST, getParticleIteratorCRS_Cell require the vector to be constructed with BIND_DEC_TO_GHOST option " << std::endl;
            ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
        }
#endif

        // First we check that
        return openfpm::vector_key_iterator_seq<typename vrl::Mem_type_type::loc_index>(NN.getParticleSeq());
    }

    template<typename Celllist> grid_key_dx<dim> getCRSStart(Celllist & NN)
    {
        return NN.getStartDomainCell();
    }

    template<typename Celllist> grid_key_dx<dim> getCRSStop(Celllist & NN)
    {
        grid_key_dx<dim> key = NN.getStopDomainCell();

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

#ifdef SE_CLASS3

    se_class3_vector<prop::max_prop,dim,St,Decomposition,self> & get_se_class3()
    {
        return se3;
    }

#endif
};


#endif /* VECTOR_HPP_ */