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"


template<unsigned int dim, typename T>
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<p_box<dim,T> > vb_int;

    CellList<dim,T,Mem_fast<>,shift<dim,T>> geo_cell;

    openfpm::vector<Point<dim,T>> shifts;

    openfpm::vector<std::pair<size_t,size_t>> ids_p;

    openfpm::vector<size_t> ids;

    shift_vect_converter<dim,T> sc_convert;

    inline size_t link_ebx_ibx(const nn_prcs<dim,T> & 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);
    }

    void Initialize_geo_cell(const Box<dim,T> & domain, const size_t (&div)[dim])
    {
        // Initialize the geo_cell structure
        geo_cell.Initialize(domain,div,0);
    }

    void create_box_nn_processor_ext(Vcluster & v_cl,
                                     Ghost<dim,T> & ghost,
                                     openfpm::vector<SpaceBox<dim,T>> & sub_domains,
                                     const openfpm::vector<openfpm::vector<long unsigned int> > & box_nn_processor,
                                     const nn_prcs<dim,T> & 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++)
        {
            SpaceBox<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<SpaceBox<dim,T>> & sub_domains,
                                     const openfpm::vector<openfpm::vector<long unsigned int> > & box_nn_processor,
                                     const nn_prcs<dim,T> & 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
                    ::SpaceBox<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
                    p_box<dim,T> b_int;
                    bool intersect = n_sub.Intersect(l_sub,b_int.box);

                    // store if it intersect
                    if (intersect == true)
                    {
                        // the box fill with the processor id
                        b_int.proc = p_id;

                        // fill the local processor id
                        b_int.lc_proc = lc_proc;

                        // fill the shift id
                        b_int.shift_id = 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.box);
                        vb_int.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.box;
                        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(b_int.box.getP1());
                        const grid_key_dx<dim> p2 = geo_cell.getCellGrid(b_int.box.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();
                            geo_cell.addCell(gi.LinId(key),vb_int.size()-1);
                            ++g_sub;
                        }
                    }
                }
            }
        }
    }


public:

    ie_ghost() {};

    ie_ghost(const ie_ghost<dim,T> & ie)
    {
        this->operator =(ie);
    }

    ie_ghost(ie_ghost<dim,T> && ie)
    {
        this->operator=(ie);
    }

    inline ie_ghost<dim,T> & operator=(ie_ghost<dim,T> && 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);
        geo_cell.swap(ie.geo_cell);
        shifts.swap(ie.shifts);
        ids_p.swap(ie.ids_p);
        ids.swap(ie.ids);

        return *this;
    }

    inline ie_ghost<dim,T> & operator=(const ie_ghost<dim,T> & 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;
        geo_cell = ie.geo_cell;
        shifts = ie.shifts;
        ids_p = ie.ids_p;
        ids = ie.ids;

        return *this;
    }

    const openfpm::vector<Point<dim,T>> & getShiftVectors()
    {
        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 const ::Box<dim,T> & getIGhostBox(size_t b_id) const
    {
        return vb_int.get(b_id).box;
    }

    inline size_t getIGhostBoxProcessor(size_t b_id) const
    {
        return vb_int.get(b_id).proc;
    }

    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.getCellIterator(geo_cell.getCell(p)))
    {
        return geo_cell.getCellIterator(geo_cell.getCell(p));
    }

    inline auto labelPoint(Point<dim,T> & p) -> decltype(geo_cell.getCellIterator(geo_cell.getCell(p)))
    {
        return geo_cell.getCellIterator(geo_cell.getCell(p));
    }

    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

        auto cell_it = geo_cell.getCellIterator(geo_cell.getCell(p));

        // 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.isInside(p) == 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

        auto cell_it = geo_cell.getCellIterator(geo_cell.getCell(p));

        // 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.isInside(p) == 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

        auto cell_it = geo_cell.getCellIterator(geo_cell.getCell(p));

        // 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.isInside(p) == 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

        auto cell_it = geo_cell.getCellIterator(geo_cell.getCell(p));

        // 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.isInside(p) == 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> & 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;
            }
        }

        return true;
    }

    bool is_equal_ng(ie_ghost<dim,T> & ig)
    {
        return true;
    }

    void reset()
    {
        box_nn_processor_int.clear();
        proc_int_box.clear();
        vb_ext.clear();
        vb_int.clear();
        geo_cell.clear();
        shifts.clear();
        ids_p.clear();
        ids.clear();
    }
};


#endif /* SRC_DECOMPOSITION_IE_GHOST_HPP_ */