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_ */