grid_dist_util.hpp

/*
 * grid_dist_util.hpp
 *
 *  Created on: Jan 28, 2016
 *      Author: i-bird
 */
 
#ifndef SRC_GRID_GRID_DIST_UTIL_HPP_
#define SRC_GRID_GRID_DIST_UTIL_HPP_
 
#include "NN/CellList/CellDecomposer.hpp"
#include "Decomposition/common.hpp"
 
template<unsigned int dim> void getCellDecomposerPar(size_t (& c_g)[dim], const size_t (& g_sz)[dim], const size_t (& bc)[dim])
{
    for (size_t i = 0 ; i < dim ; i++)
    {
        if (bc[i] == NON_PERIODIC)
            c_g[i] = (g_sz[i]-1 > 0)?(g_sz[i]-1):1;
        else
            c_g[i] = g_sz[i];
    }
}
 
 
template<unsigned int dim> periodicity<dim> create_non_periodic()
{
    periodicity<dim> p;
 
    for(size_t i = 0 ; i < dim ; i++)
        p.bc[i] = NON_PERIODIC;
 
    return p;
}
 
template<unsigned int dim>
size_t get_gdb_ext(const openfpm::vector<GBoxes<dim>> & gdb_ext, size_t start , size_t stop, size_t k)
{
    for (size_t i = start ; i < stop ; i++)
    {
        if (gdb_ext.get(i).k == k)
        {return i;}
    }
 
    // not found
    return (size_t)-1;
}
 
template<unsigned int dim>
void add_to_gdb_ext(openfpm::vector<GBoxes<dim>> & gdb_ext, size_t k, Box<dim,long int> & sp_t, Box<dim,long int> & sp_tg)
{
    // Add gdb_ext
    gdb_ext.add();
 
    gdb_ext.last().origin = sp_tg.getP1();
 
    // save information about the local grid: domain box seen inside the domain + ghost box (see GDBoxes for a visual meaning)
    // and where the GDBox start, or the origin of the local grid (+ghost) in global coordinate
    gdb_ext.last().Dbox = sp_t;
    gdb_ext.last().Dbox -= sp_tg.getP1();
 
    gdb_ext.last().GDbox = sp_tg;
    gdb_ext.last().GDbox -= sp_tg.getP1();
 
    gdb_ext.last().k = k;
}
 
template<int dim, typename Decomposition>
inline void create_gdb_ext(openfpm::vector<GBoxes<Decomposition::dims>> & gdb_ext,
                           openfpm::vector<size_t> & gdb_ext_markers,
                           Decomposition & dec,
                           CellDecomposer_sm<Decomposition::dims,typename Decomposition::stype,shift<dim,typename Decomposition::stype>> & cd_sm,
                           openfpm::vector<Box<Decomposition::dims,long int>> & bx_create,
                           const Ghost<Decomposition::dims,long int> & exp,
                           bool use_bx_def)
{
    gdb_ext.clear();
    gdb_ext_markers.clear();
 
    // Get the number of local grid needed
    size_t n_grid = dec.getNSubDomain();
 
    // Allocate the grids
    for (size_t i = 0 ; i < n_grid ; i++)
    {
        gdb_ext_markers.add(gdb_ext.size());
 
        // Get the local sub-domain (Grid conversion must be done with the domain P1 equivalent to 0.0)
        // consider that the sub-domain with point P1 equivalent to the domain P1 is a (0,0,0) in grid unit
        SpaceBox<Decomposition::dims, typename Decomposition::stype> sp = dec.getSubDomain(i);
        SpaceBox<Decomposition::dims, typename Decomposition::stype> sp_g = dec.getSubDomainWithGhost(i);
 
        // Because of round off we expand for safety the ghost area
        // std::nextafter return the next bigger or smaller representable floating
        // point number
        for (size_t i = 0 ; i < Decomposition::dims ; i++)
        {
            sp_g.setLow(i,std::nextafter(sp_g.getLow(i),sp_g.getLow(i) - 1.0));
            sp_g.setHigh(i,std::nextafter(sp_g.getHigh(i),sp_g.getHigh(i) + 1.0));
        }
 
        // Convert from SpaceBox<dim,St> to SpaceBox<dim,long int>
        SpaceBox<Decomposition::dims,long int> sp_t = cd_sm.convertDomainSpaceIntoGridUnits(sp,dec.periodicity());
        SpaceBox<Decomposition::dims,long int> sp_tg = cd_sm.convertDomainSpaceIntoGridUnits(sp_g,dec.periodicity());
 
        for (size_t i = 0 ; i < Decomposition::dims ; i++)
        {
            if (sp_t.getLow(i) < sp_tg.getLow(i))
            {sp_tg.setLow(i,sp_t.getLow(i));}
            if (sp_t.getHigh(i) > sp_tg.getHigh(i))
            {sp_tg.setHigh(i,sp_t.getHigh(i));}
        }
 
        if (use_bx_def == true)
        {
            // intersect the sub-domain with all the boxes
 
            for (size_t k = 0 ; k < bx_create.size() ; k++)
            {
                Box<Decomposition::dims, long int> inte;
 
                if (sp_t.Intersect(bx_create.get(k),inte) == true)
                {
                    // Ok we have a sub-domain now we have to create the ghost part.
                    // The ghost part is created converting the bx_def into a continuous
                    // box expanding this box by the ghost and intersecting this box
                    // with the sub-domain. This is one way to get a ghost area consistent
                    // with the construction of the external and internal ghost boxes,
                    // always calculated in continuous from the decomposition.
                    //
 
                    Box<Decomposition::dims,typename Decomposition::stype> output;
                    Box<Decomposition::dims,typename Decomposition::stype> bx_wg = cd_sm.convertCellUnitsIntoDomainSpaceMiddle(bx_create.get(k));
                    bx_wg.enlarge(dec.getGhost());
                    bx_wg.Intersect(sp_g,output);
 
                    SpaceBox<Decomposition::dims,long int> sp_t2 = inte;
                    SpaceBox<Decomposition::dims,long int> sp_tg2 =  cd_sm.convertDomainSpaceIntoGridUnits(output,dec.periodicity());
 
                    add_to_gdb_ext(gdb_ext,k,sp_t2,sp_tg2);
                }
            }
        }
        else
        {
            add_to_gdb_ext(gdb_ext,0,sp_t,sp_tg);
        }
    }
 
    gdb_ext_markers.add(gdb_ext.size());
}
 
template<int dim, typename Decomposition>
inline void create_gdb_ext(openfpm::vector<GBoxes<dim>> & gdb_ext,
                           Decomposition & dec,
                           const size_t (& sz)[dim],
                           const Box<Decomposition::dims,typename Decomposition::stype> & domain,
                           typename Decomposition::stype (& spacing)[dim])
{
    // Create the cell decomposer
    CellDecomposer_sm<Decomposition::dims,typename Decomposition::stype, shift<Decomposition::dims,typename Decomposition::stype>> cd_sm;
 
    size_t cdp[dim];
 
    // Get the parameters to create a Cell-decomposer
    getCellDecomposerPar<Decomposition::dims>(cdp,sz,dec.periodicity());
 
    // Careful cd_sm require the number of cell
    cd_sm.setDimensions(domain,cdp,0);
 
    // create an empty vector of boxes
    openfpm::vector<Box<Decomposition::dims,long int>> empty;
    Ghost<Decomposition::dims,long int> zero(0);
 
    openfpm::vector<size_t> unused;
    create_gdb_ext<dim,Decomposition>(gdb_ext,unused,dec,cd_sm,empty,zero,false);
 
    // fill the spacing
    for (size_t i = 0 ; i < dim ; i++)
    {spacing[i] = cd_sm.getCellBox().getP2()[i];}
}
 
template <unsigned int dim> struct e_box_id
{
    ::Box<dim,long int> g_e_box;
 
    ::Box<dim,long int> l_e_box;
 
    ::Box<dim,long int> lr_e_box;
 
    comb<dim> cmb;
 
    size_t g_id;
 
    size_t sub;
};
 
template<unsigned int dim, typename T, typename idT>
Box<dim,long int> flip_box(const Box<dim,idT> & box, const comb<dim> & cmb, const grid_sm<dim,T> & ginfo)
{
    Box<dim,long int> flp;
 
    for (size_t i = 0 ; i < dim; i++)
    {
        if (cmb[i] == 0)
        {
            flp.setLow(i,box.getLow(i));
            flp.setHigh(i,box.getHigh(i));
        }
        else if (cmb[i] == 1)
        {
            flp.setLow(i,box.getLow(i) + ginfo.size(i));
            flp.setHigh(i,box.getHigh(i) + ginfo.size(i));
        }
        else if (cmb[i] == -1)
        {
            flp.setLow(i,box.getLow(i) - ginfo.size(i));
            flp.setHigh(i,box.getHigh(i) - ginfo.size(i));
        }
    }
 
    return flp;
}
 
template<unsigned int dim>
inline size_t convert_to_gdb_ext(size_t sub_id,
                                 size_t def_id,
                                 openfpm::vector<GBoxes<dim>> & gdb_ext,
                                 openfpm::vector<size_t> & gdb_ext_markers)
{
    size_t start = gdb_ext_markers.get(sub_id);
    size_t stop = gdb_ext_markers.get(sub_id+1);
    return get_gdb_ext(gdb_ext,start,stop,def_id);
}
 
template <unsigned int dim> struct e_lbox_id
{
    ::Box<dim,long int> ebox;
 
    bool initialized = false;
 
    size_t sub;
 
    size_t sub_gdb_ext;
 
    size_t k;
 
    comb<dim> cmb;
};
 
 
template<unsigned int dim> inline void add_loc_eg_box(size_t le_sub,
                                               size_t se,
                                               size_t j,
                                               size_t k,
                                               openfpm::vector<e_lbox_id<dim>> & bid,
                                               const Box<dim,long int> & ebox,
                                               comb<dim> & cmb)
{
    bid.add();
 
    bid.last().ebox = ebox;
 
    bid.last().sub = se;
    bid.last().sub_gdb_ext = k;
 
    bid.last().cmb = cmb;
    bid.last().k = j;
    bid.last().initialized = true;
}
 
template<unsigned int dim> inline void add_eg_box(size_t k,
                                           const comb<dim> & cmb,
                                           const Box<dim,long int> & output,
                                           size_t g_id,
                                           const Point<dim,long int> & origin,
                                           const Point<dim,long int> & p1,
                                           openfpm::vector<e_box_id<dim>> & bid)
{
    // link
 
    size_t sub_id = k;
 
    e_box_id<dim> bid_t;
    bid_t.sub = sub_id;
    bid_t.cmb = cmb;
    bid_t.cmb.sign_flip();
    ::Box<dim,long int> ib = output;
    bid_t.g_e_box = ib;
    bid_t.g_id = g_id;
 
    // Translate in local coordinate for gdb_ext
    Box<dim,long int> tb = ib;
    tb -= origin;
    bid_t.l_e_box = tb;
 
    // Translate in local coordinates for the received box
    Box<dim,long int> tbr = ib;
    bid_t.lr_e_box = tbr;
 
    bid.add(bid_t);
}
 
template<unsigned int dim>
struct result_box
{
    size_t id;
 
    Box<dim,long int> bx;
};
 
template<unsigned int dim>
void bx_intersect(openfpm::vector<Box<dim,long int>> & bx_def,
                  bool use_bx_def,
                  Box<dim,long int> & bx,
                  openfpm::vector_std<result_box<dim>> & result)
{
    result.clear();
 
    if (use_bx_def == false)
    {
        result_box<dim> tmp;
        tmp.bx = bx;
        tmp.id = 0;
 
        result.add(tmp);
        return;
    }
 
    for (size_t i = 0 ; i < bx_def.size() ; i++)
    {
        result_box<dim> inte;
        if (bx.Intersect(bx_def.get(i),inte.bx))
        {
            inte.id = i;
            result.add(inte);
        }
    }
}
 
template<unsigned int dim> struct i_box_id
{
    ::Box<dim,long int> box;
 
    size_t g_id;
 
    size_t r_sub;
 
    comb<dim> cmb;
 
 
 
    size_t sub;
};
 
template<unsigned int dim>
struct ibox_send
{
    size_t gid;
 
    Box<dim,long int> ibox;
};
 
 
 
template<unsigned int dim> struct i_lbox_id
{
    ::Box<dim,long int> box;
 
    size_t sub;
 
    size_t sub_gdb_ext;
 
    openfpm::vector<size_t> k;
 
    comb<dim> cmb;
};
 
 
template<unsigned int dim>
struct e_box_multi
{
    openfpm::vector<size_t> eb_list;
 
    size_t full_match;
 
    size_t e_id;
};
 
 
template <unsigned int dim> struct ip_box_grid
{
    // ghost in grid units
    openfpm::vector<i_box_id<dim>> bid;
 
    size_t prc;
};
 
template <unsigned int dim> struct i_lbox_grid
{
    // ghost in grid units
    openfpm::vector<i_lbox_id<dim>> bid;
};
 
template <unsigned int dim>struct ep_box_grid
{
    // ghost in grid units
    openfpm::vector<e_box_id<dim>> bid;
 
    size_t prc;
 
    size_t recv_pnt;
 
    size_t n_r_box;
};
 
template <unsigned int dim> struct e_lbox_grid
{
    // ghost in grid units
    openfpm::vector<e_lbox_id<dim>> bid;
};
 
 
#endif /* SRC_GRID_GRID_DIST_UTIL_HPP_ */