Domain_NN_calculator_cart.hpp

/*
 * Domain_NN_calculator.hpp
 *
 *  Created on: Nov 26, 2016
 *      Author: i-bird
 */
 
#ifndef SRC_DECOMPOSITION_DOMAIN_NN_CALCULATOR_CART_HPP_
#define SRC_DECOMPOSITION_DOMAIN_NN_CALCULATOR_CART_HPP_
 
#include <Vector/map_vector.hpp>
#include "NN/CellList/ParticleItCRS_Cells.hpp"
 
template<unsigned int dim>
class domain_nn_calculator_cart
{
    bool are_domain_anom_computed;
 
 
    openfpm::vector<subsub<dim>> anom;
 
    openfpm::vector<subsub_lin<dim>> anom_lin;
 
    openfpm::vector<grid_key_dx<dim>> dom;
 
    openfpm::vector<size_t> dom_lin;
 
 
    openfpm::vector<grid_key_dx<dim>> dom_cells;
 
    openfpm::vector<size_t> dom_cells_lin;
 
 
    Box<dim,long int> proc_box;
 
    grid_sm<dim,void> gs;
 
    grid_key_dx<dim> one;
 
    void CalculateDomAndAnomCells(openfpm::vector<subsub<dim>> & sub_keys,
                                  openfpm::vector<grid_key_dx<dim>> & dom_subsub,
                                  openfpm::vector<grid_key_dx<dim>> & dom_cells,
                                  const ::Box<dim,long int> & proc_box,
                                  const openfpm::vector<::Box<dim, size_t>> & loc_box)
    {
        // Reset dom and dom_subsub
        sub_keys.clear();
        dom_subsub.clear();
 
        size_t sz[dim];
 
        // ----Grid size = proc_box.getHigh(j) - proc_box.getLow(j)
        // +2 is padding
 
        for (size_t j = 0 ; j < dim ; j++)
        {sz[j] = proc_box.getHigh(j) - proc_box.getLow(j) + 2 + 1;}
 
        gs.setDimensions(sz);
 
        // Set the grid
        grid_cpu<dim, aggregate<openfpm::vector<grid_key_dx<dim>> >> g(sz);
        g.setMemory();
 
        for (size_t i = 0 ; i < dim ; i++)
        {one.set_d(i,1);}
 
        // Calculate the csr neighborhood
        openfpm::vector<std::pair<grid_key_dx<dim>,grid_key_dx<dim>>> csr;
        NNcalc_csr(csr);
 
        // Draw the domain on this grid
        for (size_t i = 0 ; i < loc_box.size() ; i++)
        {
            grid_key_dx<dim> start;
            grid_key_dx<dim> stop;
 
            for (size_t j = 0 ; j < dim ; j++)
            {
                start.set_d(j,loc_box.template get<0>(i)[j] - proc_box.getLow(j) + 1);
                stop.set_d(j,loc_box.template get<1>(i)[j] - proc_box.getLow(j) + 1);
            }
 
            grid_key_dx_iterator_sub<dim> sub(g.getGrid(),start,stop);
 
            while (sub.isNext())
            {
                auto key = sub.get();
 
                for (size_t j = 0 ; j < csr.size() ; j++)
                {
                    grid_key_dx<dim> src = key + csr.get(j).first;
                    grid_key_dx<dim> dst = key + csr.get(j).second;
                    g.template get<0>(src).add(dst);
                }
 
                dom_cells.add(key - one);
 
                ++sub;
            }
        }
 
        // Span all the grid point and take all the sub-sub-domains that has a
        // neighborhood non-consistent to the normal symmetric NN, and the one that
        // are consistent to NN symmetric
        grid_key_dx_iterator<dim> it(g.getGrid());
 
        while (it.isNext())
        {
            auto key = it.get();
 
            // Adding in the list of the non-normal neighborhood cells
            if (g.template get<0>(key).size() == openfpm::math::pow(3,dim)/2+1)
            {
                // Add in the list of the normal neighborhood list
                dom_subsub.add(key - one);
            }
            else if (g.template get<0>(key).size() != 0)
            {
                sub_keys.add();
                sub_keys.last().subsub = key - one;
                // Adding the neighborhood of the cell
 
                sub_keys.last().NN_subsub.resize(g.template get<0>(key).size());
 
                for (size_t i = 0 ; i < g.template get<0>(key).size() ; i++)
                {sub_keys.last().NN_subsub.get(i) = g.template get<0>(key).get(i) - one;}
            }
 
            ++it;
        }
    }
 
    void linearize_subsub(const openfpm::vector<subsub<dim>> & anom,
                          openfpm::vector<subsub_lin<dim>> & anom_lin,
                          const grid_key_dx<dim> & shift,
                          const grid_sm<dim,void> & gs)
    {
        anom_lin.clear();
        for (size_t i = 0 ; i < anom.size() ; i++)
        {
            grid_key_dx<dim> tmp = anom.get(i).subsub + shift;
            anom_lin.add();
            anom_lin.last().subsub = gs.LinId(tmp);
 
            long int self_cell = -1;
 
            for (size_t j = 0 ; j < anom.get(i).NN_subsub.size() ; j++)
            {
                grid_key_dx<dim> tmp = anom.get(i).NN_subsub.get(j) + shift;
                anom_lin.get(i).NN_subsub.add((long int)gs.LinId(tmp) - anom_lin.get(i).subsub);
 
                // This indicate that for example in the neighborhood of one cell it-self is included in the list
                // For example the cell 100 is in the neighborhood of the cell 100
                if (anom_lin.get(i).NN_subsub.last() == 0)
                    self_cell = anom_lin.get(i).NN_subsub.size() - 1;
            }
 
            // if exist the self interacting cell (Example cell 100 neighborhood of cell 100), this cell MUST BE ALWAYS at the beginning
            if (self_cell != -1)
            {
                // bring the self-cell into the beginning
                size_t tmp = anom_lin.get(i).NN_subsub.get(0);
                anom_lin.get(i).NN_subsub.get(0) = 0;
                anom_lin.get(i).NN_subsub.get(self_cell) = tmp;
            }
        }
    }
 
public:
 
    domain_nn_calculator_cart()
    :are_domain_anom_computed(false)
    {
    }
 
    void setParameters(const Box<dim,long int> & proc_box)
    {
        this->proc_box = proc_box;
    }
 
    void setNNParameters(openfpm::vector<::Box<dim, size_t>> & loc_box,
                         const grid_key_dx<dim> & shift,
                         const grid_sm<dim,void> & gs)
    {
        if (are_domain_anom_computed == false)
        {
            CalculateDomAndAnomCells(anom,dom,dom_cells,proc_box,loc_box);
            are_domain_anom_computed = true;
 
            dom_cells_lin.clear();
            for (size_t i = 0 ; i < dom_cells.size() ; i++)
            {
                grid_key_dx<dim> tmp =  dom_cells.get(i) + shift;
                dom_cells_lin.add(gs.LinId(tmp));
            }
 
 
            dom_lin.clear();
            for (size_t i = 0 ; i < dom.size() ; i++)
            {
                grid_key_dx<dim> tmp = dom.get(i) + shift;
                dom_lin.add(gs.LinId(tmp));
            }
 
            linearize_subsub(anom,anom_lin,shift,gs);
        }
    }
 
    openfpm::vector<size_t> & getDomainCells()
    {
        return dom_cells_lin;
    }
 
    openfpm::vector<size_t> & getCRSDomainCells()
    {
        return dom_lin;
    }
 
    openfpm::vector<subsub_lin<dim>> & getCRSAnomDomainCells()
    {
        return anom_lin;
    }
 
 
    void reset()
    {
        are_domain_anom_computed = false;
    }
};
 
#endif /* SRC_DECOMPOSITION_DOMAIN_NN_CALCULATOR_CART_HPP_ */