ie_ghost.hpp

/*
 * ie_ghost.hpp
 *
 *  Created on: Aug 8, 2015
 *      Author: i-bird
 */
 
#ifndef SRC_DECOMPOSITION_IE_GHOST_HPP_
#define SRC_DECOMPOSITION_IE_GHOST_HPP_
 
#include "common.hpp"
#include "nn_processor.hpp"
#include "Decomposition/shift_vect_converter.hpp"
#include "Decomposition/cuda/ie_ghost_gpu.cuh"
 
struct proc_box_id
{
    size_t proc_id;
    size_t box_id;
    size_t shift_id;
 
    bool operator<(const proc_box_id & pbi) const
    {
        if (proc_id < pbi.proc_id)
        {return true;}
        else if (proc_id == pbi.proc_id)
        {
            return shift_id < pbi.shift_id;
        }
 
        return false;
    }
};
 
 
template<unsigned int dim, typename T, typename Memory, template<typename> class layout_base >
class ie_ghost
{
    openfpm::vector< openfpm::vector< Box_proc<dim,T> > > box_nn_processor_int;
 
    openfpm::vector< Box_dom<dim,T> > proc_int_box;
 
    openfpm::vector<p_box<dim,T> > vb_ext;
 
    openfpm::vector<aggregate<unsigned int,unsigned int,unsigned int>,Memory,layout_base> vb_int;
 
    openfpm::vector<Box<dim,T>,Memory,layout_base> vb_int_box;
 
    CellList<dim,T,Mem_fast<Memory,int>,shift<dim,T>> geo_cell;
 
    typedef openfpm::vector<Box<dim,T>,Memory,layout_base> proc_boxes;
 
    openfpm::vector<Point<dim,T>,Memory,layout_base> shifts;
 
    openfpm::vector<std::pair<size_t,size_t>> ids_p;
 
    openfpm::vector<size_t> ids;
 
    shift_vect_converter<dim,T,Memory,layout_base> sc_convert;
 
    bool host_dev_transfer = false;
 
    Box<dim,T> domain;
 
    size_t bc[dim];
 
    inline size_t link_ebx_ibx(const nn_prcs<dim,T,layout_base,Memory> & nn_p, size_t p_id, size_t i)
    {
        // Search for the correct id
        size_t k = 0;
        size_t p_idp = nn_p.ProctoID(p_id);
        for (k = 0 ; k < nn_p.getSentSubdomains(p_idp).size() ; k++)
        {
            if (nn_p.getSentSubdomains(p_idp).get(k) == i)
                break;
        }
        if (k == nn_p.getSentSubdomains(p_idp).size())
            std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " sub-domain not found\n";
 
        return k;
    }
 
    inline size_t ebx_ibx_form(size_t k, size_t b, size_t p_id, const comb<dim> & c ,size_t N_b, Vcluster<> & v_cl, const bool ei)
    {
        comb<dim> cext = c;
 
        if (ei == true)
            cext.sign_flip();
 
        return ((k * N_b + b) * v_cl.getProcessingUnits() + p_id) * openfpm::math::pow(3,dim) + cext.lin();
    }
 
protected:
 
    void generateShiftVectors(const Box<dim,T> & domain, size_t (& bc)[dim])
    {
        sc_convert.generateShiftVectors(domain,bc,shifts);
 
        this->domain = domain;
 
        for (int i = 0 ; i < dim ; i++)
        {
            this->bc[i] = bc[i];
        }
    }
 
    void Initialize_geo_cell(const Box<dim,T> & domain, const size_t (&div)[dim])
    {
        // Initialize the geo_cell structure
        geo_cell.Initialize(domain,div,1);
    }
 
    void free_geo_cell()
    {
        geo_cell.destroy();
    }
 
    void create_box_nn_processor_ext(Vcluster<> & v_cl,
                                     Ghost<dim,T> & ghost,
                                     openfpm::vector<Box<dim,T>,Memory,layout_base> & sub_domains,
                                     const openfpm::vector<openfpm::vector<long unsigned int> > & box_nn_processor,
                                     const nn_prcs<dim,T,layout_base,Memory> & nn_p)
    {
        box_nn_processor_int.resize(sub_domains.size());
        proc_int_box.resize(nn_p.getNNProcessors());
 
        // For each sub-domain
        for (size_t i = 0 ; i < sub_domains.size() ; i++)
        {
            Box<dim,T> sub_with_ghost = sub_domains.get(i);
 
            // enlarge the sub-domain with the ghost
            sub_with_ghost.enlarge(ghost);
 
            // resize based on the number of near processors
            box_nn_processor_int.get(i).resize(box_nn_processor.get(i).size());
 
            // For each processor near to this sub-domain
            for (size_t j = 0 ; j < box_nn_processor.get(i).size() ; j++)
            {
                // near processor
                size_t p_id = box_nn_processor.get(i).get(j);
 
                // used later
                Box_dom<dim,T> & proc_int_box_g = proc_int_box.get(nn_p.ProctoID(p_id));
 
                // Number of received sub-domains
                size_t n_r_sub = nn_p.getNRealSubdomains(p_id);
 
                // get the set of sub-domains, sector position, and real sub-domain id of the near processor p_id
                const openfpm::vector< ::Box<dim,T> > & nn_processor_subdomains_g = nn_p.getNearSubdomains(p_id);
                const openfpm::vector< comb<dim> > & nnpsg_pos = nn_p.getNearSubdomainsPos(p_id);
                const openfpm::vector< size_t > & r_sub = nn_p.getNearSubdomainsRealId(p_id);
 
                // used later
                openfpm::vector< ::Box<dim,T> > & box_nn_processor_int_gg = box_nn_processor_int.get(i).get(j).bx;
 
                // for each near processor sub-domain intersect with the enlarged local sub-domain and store it
                for (size_t b = 0 ; b < nn_processor_subdomains_g.size() ; b++)
                {
                    ::Box<dim,T> bi;
                    ::Box<dim,T> sub_bb(nn_processor_subdomains_g.get(b));
 
                    bool intersect = sub_with_ghost.Intersect(sub_bb,bi);
 
                    if (intersect == true)
                    {
                        struct p_box<dim,T> pb;
 
                        pb.box = bi;
                        pb.proc = p_id;
                        pb.lc_proc = nn_p.ProctoID(p_id);
                        pb.shift_id = (size_t)-1;
 
                        //
                        // Updating
                        //
                        // vb_ext
                        // box_nn_processor_int
                        // proc_int_box
                        //
                        // They all store the same information but organized in different ways
                        // read the description of each for more information
                        //
 
                        vb_ext.add(pb);
                        box_nn_processor_int_gg.add(bi);
                        proc_int_box_g.ebx.add();
                        proc_int_box_g.ebx.last().bx = bi;
                        proc_int_box_g.ebx.last().sub = i;
                        proc_int_box_g.ebx.last().cmb = nnpsg_pos.get(b);
 
                        // Search where the sub-domain i is in the sent list for processor p_id
                        size_t k = link_ebx_ibx(nn_p,p_id,i);
 
                        proc_int_box_g.ebx.last().id = ebx_ibx_form(k,r_sub.get(b),p_id,nnpsg_pos.get(b),n_r_sub,v_cl,true);
                    }
                }
            }
        }
    }
 
 
    void create_box_nn_processor_int(Vcluster<> & v_cl,
                                     Ghost<dim,T> & ghost,
                                     openfpm::vector<Box<dim,T>,Memory,layout_base> & sub_domains,
                                     const openfpm::vector<openfpm::vector<long unsigned int> > & box_nn_processor,
                                     const nn_prcs<dim,T,layout_base,Memory> & nn_p)
    {
        box_nn_processor_int.resize(sub_domains.size());
        proc_int_box.resize(nn_p.getNNProcessors());
 
        // For each sub-domain
        for (size_t i = 0 ; i < sub_domains.size() ; i++)
        {
            // For each processor contiguous to this sub-domain
            for (size_t j = 0 ; j < box_nn_processor.get(i).size() ; j++)
            {
                // Near processor
                size_t p_id = box_nn_processor.get(i).get(j);
 
                // get the set of sub-domains of the near processor p_id
                const openfpm::vector< ::Box<dim,T> > & nn_p_box = nn_p.getNearSubdomains(p_id);
 
                // get the sector position for each sub-domain in the list
                const openfpm::vector< comb<dim> > nn_p_box_pos = nn_p.getNearSubdomainsPos(p_id);
 
                // get the real sub-domain id for each sub-domain
                const openfpm::vector<size_t> r_sub = nn_p.getNearSubdomainsRealId(p_id);
 
                // get the local processor id
                size_t lc_proc = nn_p.getNearProcessor(p_id);
 
                // For each near processor sub-domains enlarge and intersect with the local sub-domain and store the result
                for (size_t k = 0 ; k < nn_p_box.size() ; k++)
                {
                    // enlarge the near-processor sub-domain
                    ::Box<dim,T> n_sub = nn_p_box.get(k);
 
                    // local sub-domain
                    ::Box<dim,T> l_sub = sub_domains.get(i);
 
                    // Create a margin of ghost size around the near processor sub-domain
                    n_sub.enlarge(ghost);
 
                    // Intersect with the local sub-domain
                    Box<dim,T> b_int;
                    bool intersect = n_sub.Intersect(l_sub,b_int);
 
                    // store if it intersect
                    if (intersect == true)
                    {
                        vb_int.add();
 
                        size_t last = vb_int.size() - 1;
 
                        // the box fill with the processor id
                        vb_int.template get<proc_>(last) = p_id;
 
                        // fill the local processor id
                        vb_int.template get<lc_proc_>(last) = lc_proc;
 
                        // fill the shift id
                        vb_int.template get<shift_id_>(last) = convertShift(nn_p_box_pos.get(k));
 
                        //
                        // Updating
                        //
                        // vb_int
                        // box_nn_processor_int
                        // proc_int_box
                        //
                        // They all store the same information but organized in different ways
                        // read the description of each for more information
                        //
 
                        // add the box to the near processor sub-domain intersections
                        openfpm::vector< ::Box<dim,T> > & p_box_int = box_nn_processor_int.get(i).get(j).nbx;
                        p_box_int.add(b_int);
                        vb_int_box.add(b_int);
 
                        // store the box in proc_int_box storing from which sub-domain they come from
                        Box_sub<dim,T> sb;
                        sb.bx = b_int;
                        sb.sub = i;
                        sb.r_sub = r_sub.get(k);
                        sb.cmb = nn_p_box_pos.get(k);
 
                        size_t p_idp = nn_p.ProctoID(p_id);
 
                        // Search where the sub-domain i is in the sent list for processor p_id
                        size_t s = link_ebx_ibx(nn_p,p_id,i);
 
                        // calculate the id of the internal box
                        sb.id = ebx_ibx_form(r_sub.get(k),s,v_cl.getProcessUnitID(),nn_p_box_pos.get(k),nn_p.getSentSubdomains(p_idp).size(),v_cl,false);
 
                        Box_dom<dim,T> & pr_box_int = proc_int_box.get(nn_p.ProctoID(p_id));
                        pr_box_int.ibx.add(sb);
 
                        // update the geo_cell list
 
                        // get the cells this box span
                        const grid_key_dx<dim> p1 = geo_cell.getCellGrid_me(b_int.getP1());
                        const grid_key_dx<dim> p2 = geo_cell.getCellGrid_pe(b_int.getP2());
 
                        // Get the grid and the sub-iterator
                        auto & gi = geo_cell.getGrid();
                        grid_key_dx_iterator_sub<dim> g_sub(gi,p1,p2);
 
                        // add the box-id to the cell list
                        while (g_sub.isNext())
                        {
                            auto key = g_sub.get();
                            size_t cell = gi.LinId(key);
 
                            geo_cell.addCell(cell,vb_int.size()-1);
 
                            ++g_sub;
                        }
                    }
                }
            }
        }
 
        reorder_geo_cell();
    }
 
    void reorder_geo_cell()
    {
        openfpm::vector<proc_box_id> tmp_sort;
 
        size_t div[dim];
 
        for (size_t i = 0 ; i < dim ; i++)  {div[i] = geo_cell.getDiv()[i];}
 
        grid_sm<dim,void> gs(div);
 
        grid_key_dx_iterator<dim> it(gs);
 
        while (it.isNext())
        {
            size_t cell = gs.LinId(it.get());
 
            size_t sz = geo_cell.getNelements(cell);
            tmp_sort.resize(sz);
 
            for (size_t i = 0 ; i < sz ; i++)
            {
                tmp_sort.get(i).box_id = geo_cell.get(cell,i);
                tmp_sort.get(i).proc_id = vb_int.template get<proc_>(tmp_sort.get(i).box_id);
                tmp_sort.get(i).shift_id = vb_int.template get<shift_id_>(tmp_sort.get(i).box_id);
            }
 
            tmp_sort.sort();
 
            // now we set again the cell in an ordered way
 
            for (size_t i = 0 ; i < sz ; i++)
            {geo_cell.get(cell,i) = tmp_sort.get(i).box_id;}
 
            ++it;
        }
    }
 
public:
 
    ie_ghost() {};
 
    ie_ghost(const ie_ghost<dim,T,Memory,layout_base> & ie)
    {
        this->operator=(ie);
    }
 
    ie_ghost(ie_ghost<dim,T,Memory,layout_base> && ie)
    {
        this->operator=(ie);
    }
 
    inline ie_ghost<dim,T,Memory,layout_base> & operator=(ie_ghost<dim,T,Memory,layout_base> && ie)
    {
        box_nn_processor_int.swap(ie.box_nn_processor_int);
        proc_int_box.swap(ie.proc_int_box);
        vb_ext.swap(ie.vb_ext);
        vb_int.swap(ie.vb_int);
        vb_int_box.swap(ie.vb_int_box);
        geo_cell.swap(ie.geo_cell);
        shifts.swap(ie.shifts);
        ids_p.swap(ie.ids_p);
        ids.swap(ie.ids);
 
        // it die anyway we can avoid to swap
        domain = ie.domain;
 
        for (int i = 0 ; i < dim ; i++)
        {bc[i] = ie.bc[i];}
 
        return *this;
    }
 
    inline ie_ghost<dim,T,Memory,layout_base> & operator=(const ie_ghost<dim,T,Memory,layout_base> & ie)
    {
        box_nn_processor_int = ie.box_nn_processor_int;
        proc_int_box = ie.proc_int_box;
        vb_ext = ie.vb_ext;
        vb_int = ie.vb_int;
        vb_int_box = ie.vb_int_box;
        geo_cell = ie.geo_cell;
        shifts = ie.shifts;
        ids_p = ie.ids_p;
        ids = ie.ids;
        domain = ie.domain;
 
        domain = ie.domain;
 
        for (int i = 0 ; i < dim ; i++)
        {bc[i] = ie.bc[i];}
 
        return *this;
    }
 
    template<typename Memory2, template<typename> class layout_base2>
    inline ie_ghost<dim,T,Memory,layout_base> & operator=(const ie_ghost<dim,T,Memory2,layout_base2> & ie)
    {
        box_nn_processor_int = ie.private_get_box_nn_processor_int();
        proc_int_box = ie.private_get_proc_int_box();
        vb_ext = ie.private_get_vb_ext();
        vb_int = ie.private_get_vb_int();
        vb_int_box = ie.private_get_vb_int_box();
        geo_cell = ie.private_geo_cell();
        shifts = ie.private_get_shifts();
        ids_p = ie.private_get_ids_p();
        ids = ie.private_get_ids();
        domain = ie.private_get_domain();
 
        for (int i = 0 ; i < dim ; i++)
        {bc[i] = ie.private_get_bc(i);}
 
        return *this;
    }
 
 
 
    template<typename Memory2, template <typename> class layout_base2>
    inline ie_ghost<dim,T,Memory2,layout_base2> duplicate()
    {
        ie_ghost<dim,T,Memory2,layout_base2> tmp;
 
        tmp.private_get_box_nn_processor_int() = box_nn_processor_int;
        tmp.private_get_proc_int_box() = proc_int_box;
        tmp.private_get_vb_ext() = vb_ext;
        tmp.private_get_vb_int() = vb_int;
        tmp.private_get_vb_int_box() = vb_int_box;
        tmp.private_geo_cell() = geo_cell;
        tmp.private_get_shifts() = shifts;
        tmp.private_get_ids_p() = ids_p;
        tmp.private_get_ids() = ids;
 
        tmp.private_get_domain() = domain;
 
        for (int i = 0 ; i < dim ; i++)
        {tmp.private_get_bc()[i] = bc[i];}
 
        return tmp;
    }
 
    const openfpm::vector<Point<dim,T>,Memory,layout_base> & getShiftVectors()
    {
        if (host_dev_transfer == false)
        {
            shifts.template hostToDevice<0>();
        }
 
        return shifts;
    }
 
    size_t convertShift(const comb<dim> & cmb)
    {
        return sc_convert.linId(cmb);
    }
 
    inline size_t getProcessorNIGhost(size_t id) const
    {
        return proc_int_box.get(id).ibx.size();
    }
 
    inline size_t getProcessorNEGhost(size_t id) const
    {
        return proc_int_box.get(id).ebx.size();
    }
 
    inline const ::Box<dim,T> & getProcessorIGhostBox(size_t id, size_t j) const
    {
        return proc_int_box.get(id).ibx.get(j).bx;
    }
 
    inline const ::Box<dim,T> & getProcessorEGhostBox(size_t id, size_t j) const
    {
        return proc_int_box.get(id).ebx.get(j).bx;
    }
 
    inline const comb<dim> & getProcessorEGhostPos(size_t id, size_t j) const
    {
        return proc_int_box.get(id).ebx.get(j).cmb;
    }
 
    inline const comb<dim> & getProcessorIGhostPos(size_t id, size_t j) const
    {
        return proc_int_box.get(id).ibx.get(j).cmb;
    }
 
    inline size_t getProcessorIGhostId(size_t id, size_t j) const
    {
        return proc_int_box.get(id).ibx.get(j).id;
    }
 
    inline size_t getProcessorEGhostId(size_t id, size_t j) const
    {
        return proc_int_box.get(id).ebx.get(j).id;
    }
 
    inline size_t getProcessorIGhostSSub(size_t id, size_t j) const
    {
        return proc_int_box.get(id).ibx.get(j).r_sub;
    }
 
    inline size_t getProcessorIGhostSub(size_t id, size_t j) const
    {
        return proc_int_box.get(id).ibx.get(j).sub;
    }
 
    inline size_t getProcessorEGhostSub(size_t id, size_t j) const
    {
        return proc_int_box.get(id).ebx.get(j).sub;
    }
 
    inline size_t getNIGhostBox() const
    {
        return vb_int.size();
    }
 
    inline ::Box<dim,T> getIGhostBox(size_t b_id) const
    {
        return vb_int_box.get(b_id);
    }
 
    inline size_t getIGhostBoxProcessor(size_t b_id) const
    {
        return vb_int.template get<proc_>(b_id);
    }
 
    inline size_t getNEGhostBox() const
    {
        return vb_ext.size();
    }
 
    inline ::Box<dim,T> getEGhostBox(size_t b_id) const
    {
        return vb_ext.get(b_id).box;
    }
 
    inline size_t getEGhostBoxProcessor(size_t b_id) const
    {
        return vb_ext.get(b_id).proc;
    }
 
    auto getInternalIDBoxes(Point<dim,T> & p) -> decltype(geo_cell.getParticleInCellIterator(geo_cell.getCell(p)))
    {
        return geo_cell.getParticleInCellIterator(geo_cell.getCell(p));
    }
 
    inline auto labelPoint(Point<dim,T> & p) -> decltype(geo_cell.getParticleInCellIterator(geo_cell.getCell(p)))
    {
        return geo_cell.getParticleInCellIterator(geo_cell.getCell(p));
    }
 
    template<typename output_type> inline void ghost_processor_ID(const Point<dim,T> & p, output_type & output, unsigned int base, unsigned int pi)
    {
        ID_operation<output_type> op(output);
 
        ghost_processorID_general_impl(p,base,pi,geo_cell,vb_int_box,vb_int,op);
    }
 
    inline unsigned int ghost_processorID_N(const Point<dim,T> & p)
    {
        return ghost_processorID_N_impl(p,geo_cell,vb_int_box,vb_int);
    }
 
    template <typename id1, typename id2> inline const openfpm::vector<std::pair<size_t,size_t>> ghost_processorID_pair(Point<dim,T> & p, const int opt = MULTIPLE)
    {
        ids_p.clear();
 
        // Check with geo-cell if a particle is inside one Cell containing boxes
        int cellID = geo_cell.getCell(p);
 
        if (cellID == -1)
            return ids_p;
 
        auto cell_it = geo_cell.getParticleInCellIterator(cellID);
 
        // For each element in the cell, check if the point is inside the box
        // if it is, store the processor id
        while (cell_it.isNext())
        {
            size_t bid = cell_it.get();
 
            if (Box<dim,T>(vb_int_box.get(bid)).isInsideNP_with_border(p,domain,bc) == true)
            {
                ids_p.add(std::pair<size_t,size_t>(id1::id(vb_int.get(bid),bid),id2::id(vb_int.get(bid),bid)));
            }
 
            ++cell_it;
        }
 
        // Make the id unique
        if (opt == UNIQUE)
        {
            ids_p.sort();
            ids_p.unique();
        }
 
        return ids_p;
    }
 
    template <typename id> inline const openfpm::vector<size_t> ghost_processorID(const Point<dim,T> & p, const int opt = MULTIPLE)
    {
        ids.clear();
 
        // Check with geo-cell if a particle is inside one Cell containing boxes
 
        int cellID = geo_cell.getCell(p);
 
        if (cellID == -1)
            return ids;
 
        auto cell_it = geo_cell.getParticleInCellIterator(cellID);
 
        // For each element in the cell, check if the point is inside the box
        // if it is, store the processor id
        while (cell_it.isNext())
        {
            size_t bid = cell_it.get();
 
            if (Box<dim,T>(vb_int_box.get(bid)).isInsideNP_with_border(p,domain,bc) == true)
            {
                ids.add(id::id(vb_int.get(bid),bid));
            }
 
            ++cell_it;
        }
 
        // Make the id unique
        if (opt == UNIQUE)
        {
            ids_p.sort();
            ids_p.unique();
        }
 
        return ids;
    }
 
    template<typename id1, typename id2, typename Mem> inline const openfpm::vector<std::pair<size_t,size_t>> & ghost_processorID_pair(const encapc<1,Point<dim,T>,Mem> & p, const int opt = MULTIPLE)
    {
        ids_p.clear();
 
        // Check with geo-cell if a particle is inside one Cell containing boxes
 
        int cellID = geo_cell.getCell(p);
 
        if (cellID == -1)
            return ids_p;
 
        auto cell_it = geo_cell.getParticleInCellIterator(cellID);
 
        // For each element in the cell, check if the point is inside the box
        // if it is, store the processor id
        while (cell_it.isNext())
        {
            size_t bid = cell_it.get();
 
            if (Box<dim,T>(vb_int_box.get(bid)).isInsideNP_with_border(p,domain,bc) == true)
            {
                ids_p.add(std::pair<size_t,size_t>(id1::id(vb_int.get(bid),bid),id2::id(vb_int.get(bid),bid)));
            }
 
            ++cell_it;
        }
 
        // Make the id unique
        if (opt == UNIQUE)
        {
            ids_p.sort();
            ids_p.unique();
        }
 
        return ids_p;
    }
 
    template<typename id, typename Mem> inline const openfpm::vector<size_t> & ghost_processorID(const encapc<1,Point<dim,T>,Mem> & p, const int opt = MULTIPLE)
    {
        ids.clear();
 
        // Check with geo-cell if a particle is inside one Cell containing boxes
 
        int cellID = geo_cell.getCell(p);
 
        if (cellID == -1)
            return ids;
 
        auto cell_it = geo_cell.getParticleInCellIterator(cellID);
 
        // For each element in the cell, check if the point is inside the box
        // if it is, store the processor id
        while (cell_it.isNext())
        {
            size_t bid = cell_it.get();
 
            if (vb_int.get(bid).box.isInsideNP_with_border(p,bc,domain) == true)
            {
                ids.add(id::id(vb_int.get(bid),bid));
            }
 
            ++cell_it;
        }
 
        // Make the id unique
        if (opt == UNIQUE)
        {
            ids_p.sort();
            ids_p.unique();
        }
 
        return ids;
    }
 
    bool write(std::string output, size_t p_id) const
    {
        VTKWriter<openfpm::vector<::Box<dim,T>>,VECTOR_BOX> vtk_box3;
        for (size_t p = 0 ; p < box_nn_processor_int.size() ; p++)
        {
            for (size_t s = 0 ; s < box_nn_processor_int.get(p).size() ; s++)
            {
                vtk_box3.add(box_nn_processor_int.get(p).get(s).nbx);
            }
        }
        vtk_box3.write(output + std::string("internal_ghost_") + std::to_string(p_id) + std::string(".vtk"));
 
        VTKWriter<openfpm::vector<::Box<dim,T>>,VECTOR_BOX> vtk_box4;
        for (size_t p = 0 ; p < box_nn_processor_int.size() ; p++)
        {
            for (size_t s = 0 ; s < box_nn_processor_int.get(p).size() ; s++)
            {
                vtk_box4.add(box_nn_processor_int.get(p).get(s).bx);
            }
        }
        vtk_box4.write(output + std::string("external_ghost_") + std::to_string(p_id) + std::string(".vtk"));
 
        return true;
    }
 
    bool is_equal(ie_ghost<dim,T,Memory,layout_base> & ig)
    {
        if (getNEGhostBox() != ig.getNEGhostBox())
            return false;
 
        if (getNIGhostBox() != ig.getNIGhostBox())
            return false;
 
        for (size_t i = 0 ; i < getNIGhostBox() ; i++)
        {
            if (getIGhostBox(i) != ig.getIGhostBox(i))
                return false;
            if (getIGhostBoxProcessor(i) != ig.getIGhostBoxProcessor(i))
                return false;
        }
 
        for (size_t i = 0 ; i < proc_int_box.size() ; i++)
        {
            if (getProcessorNIGhost(i) != ig.getProcessorNIGhost(i))
                return false;
            for (size_t j = 0 ; j < getProcessorNIGhost(i) ; j++)
            {
                if (getProcessorIGhostBox(i,j) != ig.getProcessorIGhostBox(i,j))
                    return false;
                if (getProcessorIGhostId(i,j) != ig.getProcessorIGhostId(i,j))
                    return false;
                if (getProcessorIGhostSub(i,j) != ig.getProcessorIGhostSub(i,j))
                    return false;
            }
        }
 
        for (size_t i = 0 ; i < getNEGhostBox() ; i++)
        {
            if (getEGhostBox(i) != ig.getEGhostBox(i))
                return false;
            if (getEGhostBoxProcessor(i) != ig.getEGhostBoxProcessor(i))
                return false;
        }
 
        for (size_t i = 0 ; i < proc_int_box.size() ; i++)
        {
            if (getProcessorNEGhost(i) != ig.getProcessorNEGhost(i))
                return false;
            for (size_t j = 0 ; j < getProcessorNEGhost(i) ; j++)
            {
                if (getProcessorEGhostBox(i,j) != ig.getProcessorEGhostBox(i,j))
                    return false;
                if (getProcessorEGhostId(i,j) != ig.getProcessorEGhostId(i,j))
                    return false;
                if (getProcessorEGhostSub(i,j) != ig.getProcessorEGhostSub(i,j))
                    return false;
            }
        }
 
        if (domain != ig.domain)
        {return false;}
 
        for (int i = 0 ; i < dim ; i++)
        {
            if (bc[i] != ig.bc[i])
            {return false;}
        }
 
        return true;
    }
 
    bool is_equal_ng(ie_ghost<dim,T,Memory,layout_base> & ig)
    {
        return true;
    }
 
    void reset()
    {
        box_nn_processor_int.clear();
        proc_int_box.clear();
        vb_ext.clear();
        vb_int.clear();
        vb_int_box.clear();
        geo_cell.clear();
        shifts.clear();
        ids_p.clear();
        ids.clear();
    }
 
    inline const openfpm::vector< openfpm::vector< Box_proc<dim,T> > > & private_get_box_nn_processor_int() const
    {
        return box_nn_processor_int;
    }
 
    inline openfpm::vector< openfpm::vector< Box_proc<dim,T> > > & private_get_box_nn_processor_int()
    {
        return box_nn_processor_int;
    }
 
    inline openfpm::vector< Box_dom<dim,T> > & private_get_proc_int_box()
    {
        return proc_int_box;
    }
 
    inline const openfpm::vector< Box_dom<dim,T> > & private_get_proc_int_box() const
    {
        return proc_int_box;
    }
 
    inline openfpm::vector<p_box<dim,T> > & private_get_vb_ext()
    {
        return vb_ext;
    }
 
    inline const openfpm::vector<p_box<dim,T> > & private_get_vb_ext() const
    {
        return vb_ext;
    }
 
    inline openfpm::vector<aggregate<unsigned int,unsigned int,unsigned int>,Memory,layout_base> &
    private_get_vb_int()
    {
        return vb_int;
    }
 
    inline const openfpm::vector<aggregate<unsigned int,unsigned int,unsigned int>,Memory,layout_base> &
    private_get_vb_int() const
    {
        return vb_int;
    }
 
    inline openfpm::vector<Box<dim,T>,Memory,layout_base> &
    private_get_vb_int_box()
    {
        return vb_int_box;
    }
 
    inline const openfpm::vector<Box<dim,T>,Memory,layout_base> &
    private_get_vb_int_box() const
    {
        return vb_int_box;
    }
 
    inline CellList<dim,T,Mem_fast<Memory,int>,shift<dim,T>> &
    private_geo_cell()
    {
        return geo_cell;
    }
 
    inline const CellList<dim,T,Mem_fast<Memory,int>,shift<dim,T>> &
    private_geo_cell() const
    {
        return geo_cell;
    }
 
    inline openfpm::vector<Point<dim,T>,Memory,layout_base> &
    private_get_shifts()
    {
        return shifts;
    }
 
    inline const openfpm::vector<Point<dim,T>,Memory,layout_base> &
    private_get_shifts() const
    {
        return shifts;
    }
 
    inline openfpm::vector<std::pair<size_t,size_t>> &
    private_get_ids_p()
    {
        return ids_p;
    }
 
    inline const openfpm::vector<std::pair<size_t,size_t>> &
    private_get_ids_p() const
    {
        return ids_p;
    }
 
    inline openfpm::vector<size_t> &
    private_get_ids()
    {
        return ids;
    }
 
    inline const openfpm::vector<size_t> &
    private_get_ids() const
    {
        return ids;
    }
 
    inline Box<dim,T> &
    private_get_domain()
    {
        return domain;
    }
 
    inline const Box<dim,T> &
    private_get_domain() const
    {
        return domain;
    }
 
    size_t private_get_bc(int i) const
    {
        return bc[i];
    }
 
    inline  size_t (& private_get_bc())[dim]
    {
        return bc;
    }
 
    ie_ghost_gpu<dim,T,Memory,layout_base> toKernel()
    {
        if (host_dev_transfer == false)
        {
            geo_cell.hostToDevice();
            vb_int_box.template hostToDevice<0,1>();
            vb_int.template hostToDevice<0,1,2>();
            shifts.template hostToDevice<0>();
 
            host_dev_transfer = true;
        }
 
        ie_ghost_gpu<dim,T,Memory,layout_base> igg(geo_cell.toKernel(),
                                                    vb_int_box.toKernel(),
                                                    vb_int.toKernel());
 
        return igg;
    }
 
    void reset_host_dev_transfer()
    {
        host_dev_transfer = false;
    }
};
 
 
#endif /* SRC_DECOMPOSITION_IE_GHOST_HPP_ */