CartDecomposition.hpp

/*
 * CartDecomposition.hpp
 *
 *  Created on: Oct 07, 2015
 *      Author: Pietro Incardona, Antonio Leo
 */

#ifndef CARTDECOMPOSITION_HPP
#define CARTDECOMPOSITION_HPP

#include "config.h"
#include <cmath>
#include "VCluster/VCluster.hpp"
#include "Graph/CartesianGraphFactory.hpp"
#include "Decomposition.hpp"
#include "Vector/map_vector.hpp"
#include <vector>
#include <initializer_list>
#include "SubdomainGraphNodes.hpp"
#include "dec_optimizer.hpp"
#include "Space/Shape/Box.hpp"
#include "Space/Shape/Point.hpp"
#include "NN/CellList/CellDecomposer.hpp"
#include <unordered_map>
#include "NN/CellList/CellList.hpp"
#include "Space/Ghost.hpp"
#include "common.hpp"
#include "ie_loc_ghost.hpp"
#include "ie_ghost.hpp"
#include "nn_processor.hpp"
#include "GraphMLWriter/GraphMLWriter.hpp"
#include "Distribution/ParMetisDistribution.hpp"
#include "Distribution/DistParMetisDistribution.hpp"
#include "Distribution/MetisDistribution.hpp"
#include "DLB/DLB.hpp"
#include "util/se_util.hpp"
#include "util/mathutil.hpp"
#include "CartDecomposition_ext.hpp"
#include "data_type/aggregate.hpp"
#include "Domain_NN_calculator_cart.hpp"

#define CARTDEC_ERROR 2000lu

template<unsigned int dim> static void nsub_to_div2(size_t (& div)[dim], size_t n_sub, size_t dim_r)
{
    for (size_t i = 0; i < dim; i++)
    {
        if (i < dim_r)
        {div[i] = openfpm::math::round_big_2(pow(n_sub, 1.0 / dim_r));}
        else
        {div[i] = 1;}
    }
}

template<unsigned int dim> static void nsub_to_div(size_t (& div)[dim], size_t n_sub, size_t dim_r)
{
    for (size_t i = 0; i < dim; i++)
    {
        if (i < dim_r)
        {div[i] = std::floor(pow(n_sub, 1.0 / dim_r));}
        else
        {div[i] = 1;}
    }
}

#define COMPUTE_SKIN_SUB 1

template<unsigned int dim, typename T, typename Memory, typename Distribution>
class CartDecomposition: public ie_loc_ghost<dim, T>, public nn_prcs<dim, T>, public ie_ghost<dim, T>, public domain_nn_calculator_cart<dim>
{

public:

    typedef T domain_type;

    typedef SpaceBox<dim, T> Box;

    typedef CartDecomposition<dim,T,Memory,Distribution> base_type;

    typedef CartDecomposition_ext<dim,T,Memory,Distribution> extended_type;

protected:

    bool commCostSet = false;

    typedef typename openfpm::vector<SpaceBox<dim, T>,
            Memory,
            typename memory_traits_lin<SpaceBox<dim, T>>::type,
            memory_traits_lin,
            openfpm::vector_grow_policy_default,
            openfpm::vect_isel<SpaceBox<dim, T>>::value>::access_key acc_key;

    openfpm::vector<SpaceBox<dim, T>> sub_domains;

    mutable openfpm::vector<openfpm::vector<SpaceBox<dim, T>>> sub_domains_global;

    openfpm::vector<openfpm::vector<long unsigned int> > box_nn_processor;

    openfpm::vector<size_t> fine_s;

    grid_sm<dim, void> gr;

    grid_sm<dim, void> gr_dist;

    CellDecomposer_sm<dim, T, shift<dim,T>> cd;

    ::Box<dim,T> domain;

    T spacing[dim];

    size_t magn[dim];

    Vcluster & v_cl;

    Distribution dist;

    ::Box<dim,T> bbox;

    long int ref_cnt;

    Ghost<dim,T> ghost;

    size_t bc[dim];

    ::Box<dim,size_t> proc_box;

    openfpm::vector<::Box<dim, size_t>> loc_box;

    template<typename Memory_bx> SpaceBox<dim,T> convertDecBoxIntoSubDomain(encapc<1,::Box<dim,size_t>,Memory_bx> loc_box)
    {
        // A point with all coordinate to one
        size_t one[dim];
        for (size_t i = 0 ; i < dim ; i++)  {one[i] = 1;}

        SpaceBox<dim, size_t> sub_dc = loc_box;
        SpaceBox<dim, size_t> sub_dce = sub_dc;
        sub_dce.expand(one);
        sub_dce.mul(magn);

        // shrink by one
        for (size_t i = 0 ; i < dim ; i++)
        {
            loc_box.template get<Box::p1>()[i] = sub_dce.getLow(i);
            loc_box.template get<Box::p2>()[i] = sub_dce.getHigh(i) - 1;
        }

        SpaceBox<dim, T> sub_d(sub_dce);
        sub_d.mul(spacing);
        sub_d += domain.getP1();

        // we add the

        // Fixing sub-domains to cover all the domain

        // Fixing sub_d
        // if (loc_box) is at the boundary we have to ensure that the box span the full
        // domain (avoiding rounding off error)
        for (size_t i = 0; i < dim; i++)
        {
            if (sub_dc.getHigh(i) == gr.size(i) - 1)
            {sub_d.setHigh(i, domain.getHigh(i));}

            if (sub_dc.getLow(i) == 0)
            {sub_d.setLow(i,domain.getLow(i));}
        }

        return sub_d;
    }


public:

    void createSubdomains(Vcluster & v_cl, const size_t (& bc)[dim], size_t opt = 0)
    {
        int p_id = v_cl.getProcessUnitID();

        // Calculate the total number of box and and the spacing
        // on each direction
        // Get the box containing the domain
        SpaceBox<dim, T> bs = domain.getBox();

        for (unsigned int i = 0; i < dim; i++)
        {
            // Calculate the spacing
            spacing[i] = (bs.getHigh(i) - bs.getLow(i)) / gr.size(i);
        }

        // fill the structure that store the processor id for each sub-domain
        fine_s.resize(gr.size());

        // Optimize the decomposition creating bigger spaces
        // And reducing Ghost over-stress
        dec_optimizer<dim, Graph_CSR<nm_v, nm_e>> d_o(dist.getGraph(), gr_dist.getSize());

        // Ghost
        Ghost<dim,long int> ghe;

        // Set the ghost
        for (size_t i = 0 ; i < dim ; i++)
        {
            ghe.setLow(i,static_cast<long int>(ghost.getLow(i)/spacing[i]) - 1);
            ghe.setHigh(i,static_cast<long int>(ghost.getHigh(i)/spacing[i]) + 1);
        }

        // optimize the decomposition
        d_o.template optimize<nm_v::sub_id, nm_v::proc_id>(dist.getGraph(), p_id, loc_box, box_nn_processor,ghe,bc);

        // Initialize
        if (loc_box.size() > 0)
        {
            bbox = convertDecBoxIntoSubDomain(loc_box.get(0));
            proc_box = loc_box.get(0);
            sub_domains.add(bbox);
        }
        else
        {
            // invalidate all the boxes
            for (size_t i = 0 ; i < dim ; i++)
            {
                proc_box.setLow(i,0.0);
                proc_box.setHigh(i,0);

                bbox.setLow(i,0.0);
                bbox.setHigh(i,0);
            }
        }

        // convert into sub-domain
        for (size_t s = 1; s < loc_box.size(); s++)
        {
            SpaceBox<dim,T> sub_d = convertDecBoxIntoSubDomain(loc_box.get(s));

            // add the sub-domain
            sub_domains.add(sub_d);

            // Calculate the bound box
            bbox.enclose(sub_d);
            proc_box.enclose(loc_box.get(s));
        }

        nn_prcs<dim,T>::create(box_nn_processor, sub_domains);
        nn_prcs<dim,T>::applyBC(domain,ghost,bc);

        // fill fine_s structure
        // fine_s structure contain the processor id for each sub-sub-domain
        // with sub-sub-domain we mean the sub-domain decomposition before
        // running dec_optimizer (before merging sub-domains)


        grid_key_dx_iterator<dim> git(gr);

        while (git.isNext())
        {
            auto key = git.get();
            grid_key_dx<dim> key2;

            for (size_t i = 0 ; i < dim ; i++)
                key2.set_d(i,key.get(i) / magn[i]);

            size_t lin = gr_dist.LinId(key2);
            size_t lin2 = gr.LinId(key);

            fine_s.get(lin2) = dist.getGraph().template vertex_p<nm_v::proc_id>(lin);

            ++git;
        }

        Initialize_geo_cell_lists();
    }

    void Initialize_geo_cell_lists()
    {
        // Get the processor bounding Box
        ::Box<dim,T> bound = getProcessorBounds();

        // Check if the box is valid
        if (bound.isValidN() == true)
        {
            // Not necessary, but I prefer
            bound.enlarge(ghost);

            // calculate the sub-divisions
            size_t div[dim];
            for (size_t i = 0; i < dim; i++)
            {div[i] = (size_t) ((bound.getHigh(i) - bound.getLow(i)) / cd.getCellBox().getP2()[i]);}

            // Initialize the geo_cell structure
            ie_ghost<dim,T>::Initialize_geo_cell(bound,div);

            // Initialize shift vectors
            ie_ghost<dim,T>::generateShiftVectors(domain,bc);
        }
    }

    void computeCommunicationAndMigrationCosts(size_t ts)
    {
        float migration = 0;

        SpaceBox<dim, T> cellBox = cd.getCellBox();
        float b_s = static_cast<float>(cellBox.getHigh(0));
        float gh_s = static_cast<float>(ghost.getHigh(0));

        // compute the gh_area for 2 dim case
        float gh_v = (gh_s * b_s);

        // multiply for sub-sub-domain side for each domain
        for (size_t i = 2; i < dim; i++)
        {
            /* coverity[dead_error_line] */
            gh_v *= b_s;
        }

        size_t norm = (size_t) (1.0 / gh_v);

        migration = pow(b_s, dim);

        size_t prev = 0;

        for (size_t i = 0; i < dist.getNSubSubDomains(); i++)
        {
            dist.setMigrationCost(i, norm * migration /* * dist.getSubSubDomainComputationCost(i)*/ );

            for (size_t s = 0; s < dist.getNSubSubDomainNeighbors(i); s++)
            {
                // We have to remove dist.getSubSubDomainComputationCost(i) otherwise the graph is
                // not directed
                dist.setCommunicationCost(i, s, 1   *  ts);
            }
            prev += dist.getNSubSubDomainNeighbors(i);
        }

        commCostSet = true;
    }

    void CreateSubspaces()
    {
        // Create a grid where each point is a space
        grid_sm<dim, void> g(div);

        // create a grid_key_dx iterator
        grid_key_dx_iterator<dim> gk_it(g);

        // Divide the space into subspaces
        while (gk_it.isNext())
        {
            grid_key_dx<dim> key = gk_it.get();

            SpaceBox<dim, T> tmp;

            for (int i = 0; i < dim; i++)
            {
                tmp.setHigh(i, (key.get(i) + 1) * spacing[i]);
                tmp.setLow(i, key.get(i) * spacing[i]);
            }

            sub_domains.add(tmp);

            // Next sub-domain
            ++gk_it;
        }
    }


    void calculateGhostBoxes()
    {
        // Intersect all the local sub-domains with the sub-domains of the contiguous processors

        // create the internal structures that store ghost information
        ie_ghost<dim, T>::create_box_nn_processor_ext(v_cl, ghost, sub_domains, box_nn_processor, *this);
        ie_ghost<dim, T>::create_box_nn_processor_int(v_cl, ghost, sub_domains, box_nn_processor, *this);

        ie_loc_ghost<dim,T>::create(sub_domains,domain,ghost,bc);
    }

public:

    static constexpr int dims = dim;

    typedef T stype;

    void incRef()
    {ref_cnt++;}

    void decRef()
    {ref_cnt--;}

    long int ref()
    {
        return ref_cnt;
    }

    CartDecomposition(Vcluster & v_cl)
    :nn_prcs<dim, T>(v_cl), v_cl(v_cl), dist(v_cl),ref_cnt(0)
    {
        // Reset the box to zero
        bbox.zero();
    }

    CartDecomposition(const CartDecomposition<dim,T,Memory,Distribution> & cart)
    :nn_prcs<dim,T>(cart.v_cl),v_cl(cart.v_cl),dist(v_cl),ref_cnt(0)
    {
        this->operator=(cart);
    }

    CartDecomposition(CartDecomposition<dim,T,Memory,Distribution> && cart)
    :nn_prcs<dim,T>(cart.v_cl),v_cl(cart.v_cl),dist(v_cl),ref_cnt(0)
    {
        this->operator=(cart);
    }

    ~CartDecomposition()
    {
    }

    class box_id
    {
    public:
        inline static size_t id(p_box<dim, T> & p, size_t b_id)
        {
            return b_id;
        }
    };

    class processor_id
    {
    public:
        inline static size_t id(p_box<dim, T> & p, size_t b_id)
        {
            return p.proc;
        }
    };

    class lc_processor_id
    {
    public:
        inline static size_t id(p_box<dim, T> & p, size_t b_id)
        {
            return p.lc_proc;
        }
    };

    class shift_id
    {
    public:
        inline static size_t id(p_box<dim,T> & p, size_t b_id)
        {
            return p.shift_id;
        }
    };

    void applyPointBC(float (& pt)[dim]) const
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (bc[i] == PERIODIC)
                pt[i] = openfpm::math::periodic_l(pt[i],domain.getHigh(i),domain.getLow(i));
        }
    }

    void applyPointBC(Point<dim,T> & pt) const
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (bc[i] == PERIODIC)
                pt.get(i) = openfpm::math::periodic_l(pt.get(i),domain.getHigh(i),domain.getLow(i));
        }
    }

    template<typename Mem> void applyPointBC(encapc<1,Point<dim,T>,Mem> && pt) const
    {
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (bc[i] == PERIODIC)
                pt.template get<0>()[i] = openfpm::math::periodic_l(pt.template get<0>()[i],domain.getHigh(i),domain.getLow(i));
        }
    }

    CartDecomposition<dim,T,Memory,Distribution> duplicate(const Ghost<dim,T> & g) const
    {
        CartDecomposition<dim,T,Memory,Distribution> cart(v_cl);

        cart.box_nn_processor = box_nn_processor;
        cart.sub_domains = sub_domains;
        cart.fine_s = fine_s;

        cart.gr = gr;
        cart.cd = cd;
        cart.domain = domain;
        for (size_t i = 0 ; i < dim ; i++)
        {cart.spacing[i] = spacing[i];};

        cart.bbox = bbox;
        cart.ghost = g;

        cart.dist = dist;

        for (size_t i = 0 ; i < dim ; i++)
            cart.bc[i] = bc[i];

        (static_cast<nn_prcs<dim,T> &>(cart)).create(box_nn_processor, sub_domains);
        (static_cast<nn_prcs<dim,T> &>(cart)).applyBC(domain,ghost,bc);

        cart.Initialize_geo_cell_lists();
        cart.calculateGhostBoxes();

        return cart;
    }

    CartDecomposition<dim,T,Memory,Distribution> duplicate() const
    {
        CartDecomposition<dim,T,Memory,Distribution> cart(v_cl);

        (static_cast<ie_loc_ghost<dim,T>*>(&cart))->operator=(static_cast<ie_loc_ghost<dim,T>>(*this));
        (static_cast<nn_prcs<dim,T>*>(&cart))->operator=(static_cast<nn_prcs<dim,T>>(*this));
        (static_cast<ie_ghost<dim,T>*>(&cart))->operator=(static_cast<ie_ghost<dim,T>>(*this));

        cart.sub_domains = sub_domains;
        cart.box_nn_processor = box_nn_processor;
        cart.fine_s = fine_s;
        cart.gr = gr;
        cart.gr_dist = gr_dist;
        cart.dist = dist;
        cart.commCostSet = commCostSet;
        cart.cd = cd;
        cart.domain = domain;
        for (size_t i = 0 ; i < dim ; i++)
        {cart.spacing[i] = spacing[i];};

        cart.ghost = ghost;

        cart.bbox = bbox;

        for (size_t i = 0 ; i < dim ; i++)
            cart.bc[i] = this->bc[i];

        return cart;
    }

    CartDecomposition<dim,T,Memory, Distribution> & operator=(const CartDecomposition & cart)
    {
        static_cast<ie_loc_ghost<dim,T>*>(this)->operator=(static_cast<ie_loc_ghost<dim,T>>(cart));
        static_cast<nn_prcs<dim,T>*>(this)->operator=(static_cast<nn_prcs<dim,T>>(cart));
        static_cast<ie_ghost<dim,T>*>(this)->operator=(static_cast<ie_ghost<dim,T>>(cart));

        sub_domains = cart.sub_domains;
        box_nn_processor = cart.box_nn_processor;
        fine_s = cart.fine_s;
        gr = cart.gr;
        gr_dist = cart.gr_dist;
        dist = cart.dist;
        commCostSet = cart.commCostSet;
        cd = cart.cd;
        domain = cart.domain;

        for (size_t i = 0 ; i < dim ; i++)
        {
            spacing[i] = cart.spacing[i];
            magn[i] = cart.magn[i];
        };

        ghost = cart.ghost;

        bbox = cart.bbox;

        for (size_t i = 0 ; i < dim ; i++)
            bc[i] = cart.bc[i];

        return *this;
    }

    CartDecomposition<dim,T,Memory,Distribution> & operator=(CartDecomposition && cart)
    {
        static_cast<ie_loc_ghost<dim,T>*>(this)->operator=(static_cast<ie_loc_ghost<dim,T>>(cart));
        static_cast<nn_prcs<dim,T>*>(this)->operator=(static_cast<nn_prcs<dim,T>>(cart));
        static_cast<ie_ghost<dim,T>*>(this)->operator=(static_cast<ie_ghost<dim,T>>(cart));

        sub_domains.swap(cart.sub_domains);
        box_nn_processor.swap(cart.box_nn_processor);
        fine_s.swap(cart.fine_s);
        gr = cart.gr;
        gr_dist = cart.gr_dist;
        dist = cart.dist;
        commCostSet = cart.commCostSet;
        cd = cart.cd;
        domain = cart.domain;
        for (size_t i = 0 ; i < dim ; i++)
        {
            spacing[i] = cart.spacing[i];
            magn[i] = cart.magn[i];
        };

        ghost = cart.ghost;

        bbox = cart.bbox;

        for (size_t i = 0 ; i < dim ; i++)
            bc[i] = cart.bc[i];

        return *this;
    }

    static size_t getDefaultGrid(size_t n_sub)
    {
        // Calculate the number of sub-sub-domain on
        // each dimension
        return openfpm::math::round_big_2(pow(n_sub, 1.0 / dim));
    }

    template<typename Mem> size_t inline processorID(const encapc<1, Point<dim,T>, Mem> & p) const
    {
        return fine_s.get(cd.template getCell(p));
    }

    size_t inline processorID(const Point<dim,T> &p) const
    {
        return fine_s.get(cd.getCell(p));
    }

    size_t inline processorID(const T (&p)[dim]) const
    {
        return fine_s.get(cd.getCell(p));
    }

    template<typename Mem> size_t inline processorIDBC(encapc<1, Point<dim,T>, Mem> p)
    {
        Point<dim,T> pt = p;
        applyPointBC(pt);

        return fine_s.get(cd.getCell(pt));
    }

    template<typename ofb> size_t inline processorIDBC(const Point<dim,T> &p) const
    {
        Point<dim,T> pt = p;
        applyPointBC(pt);

        return fine_s.get(cd.getCell(p));
    }

    template<typename ofb> size_t inline processorIDBC(const T (&p)[dim]) const
    {
        Point<dim,T> pt = p;
        applyPointBC(pt);

        return fine_s.get(cd.getCell(p));
    }

    inline size_t periodicity(size_t i) const
    {
        return bc[i];
    }

    inline const size_t (& periodicity() const) [dim]
    {
        return bc;
    }

    void calculate_magn(const grid_sm<dim,void> & gm)
    {
        if (gm.size() == 0)
        {
            for (size_t i = 0 ; i < dim ; i++)
                magn[i] = 1;
        }
        else
        {
            for (size_t i = 0 ; i < dim ; i++)
            {
                if (gr.size(i) % gm.size(i) != 0)
                    std::cerr << __FILE__ << ":" << __LINE__ << ".Error the decomposition grid specified as gr.size(" << i << ")=" << gr.size(i) << " is not multiple of the distribution grid gm.size(" << i << ")=" << gm.size(i) << std::endl;

                magn[i] = gr.size(i) / gm.size(i);
            }
        }
    }

    void setGoodParameters(::Box<dim,T> domain_,
                           const size_t (& bc)[dim],
                           const Ghost<dim,T> & ghost,
                           size_t dec_gran,
                           const grid_sm<dim,void> & sec_dist = grid_sm<dim,void>())
    {
        size_t div[dim];

        // Create a valid decomposition of the space
        // Get the number of processor and calculate the number of sub-domain
        // for decomposition
        size_t n_proc = v_cl.getProcessingUnits();
        size_t n_sub = n_proc * dec_gran;

        // Calculate the maximum number (before merging) of sub-domain on
        // each dimension

        nsub_to_div2(div,n_sub,dim);

/*      for (size_t i = 0; i < dim; i++)
        {
            div[i] = openfpm::math::round_big_2(pow(n_sub, 1.0 / dim));
        }*/

        if (dim > 3)
        {
            long int dim_r = dim-1;
            do
            {
                // Check for adjustment
                size_t tot_size = 1;
                for (size_t i = 0 ; i < dim ; i++)
                {tot_size *= div[i];}

                // the granularity is too coarse increase the divisions
                if (tot_size / n_proc > 0.75*dec_gran )
                {break;}

                nsub_to_div(div,n_sub,dim_r);

                dim_r--;
            } while(dim_r > 0);
        }

        setParameters(div,domain_,bc,ghost,sec_dist);
    }

    void getParameters(size_t (& div_)[dim])
    {
        for (size_t i = 0 ; i < dim ; i++)
        {div_[i] = this->gr.size(i);}
    }

    void setParameters(const size_t (& div_)[dim],
                       ::Box<dim,T> domain_,
                        const size_t (& bc)[dim],
                        const Ghost<dim,T> & ghost,
                        const grid_sm<dim,void> & sec_dist = grid_sm<dim,void>())
    {
        // set the boundary conditions
        for (size_t i = 0 ; i < dim ; i++)
            this->bc[i] = bc[i];

        // set the ghost
        this->ghost = ghost;

        // Set the decomposition parameters
        gr.setDimensions(div_);
        domain = domain_;
        cd.setDimensions(domain, div_, 0);

        // We we have a secondary grid costruct a reduced graph
        if (sec_dist.size(0) != 0)
        {
            calculate_magn(sec_dist);
            gr_dist.setDimensions(sec_dist.getSize());
        }
        else
        {
            calculate_magn(sec_dist);
            gr_dist = gr;
        }

        // init distribution
        dist.createCartGraph(gr_dist, domain);

    }

    void reset()
    {
        sub_domains.clear();
        box_nn_processor.clear();
        fine_s.clear();
        loc_box.clear();
        nn_prcs<dim, T>::reset();
        ie_ghost<dim, T>::reset();
        ie_loc_ghost<dim, T>::reset();
    }

    void decompose()
    {
        reset();

        if (commCostSet == false)
        {computeCommunicationAndMigrationCosts(1);}

        dist.decompose();

        createSubdomains(v_cl,bc);

        calculateGhostBoxes();

        domain_nn_calculator_cart<dim>::reset();
        domain_nn_calculator_cart<dim>::setParameters(proc_box);
    }

    void refine(size_t ts)
    {
        reset();

        if (commCostSet == false)
        {computeCommunicationAndMigrationCosts(ts);}

        dist.refine();

        createSubdomains(v_cl,bc);

        calculateGhostBoxes();

        domain_nn_calculator_cart<dim>::reset();
        domain_nn_calculator_cart<dim>::setParameters(proc_box);
    }

    void redecompose(size_t ts)
    {
        reset();

        if (commCostSet == false)
        {computeCommunicationAndMigrationCosts(ts);}

        dist.redecompose();

        createSubdomains(v_cl,bc);

        calculateGhostBoxes();

        domain_nn_calculator_cart<dim>::reset();
        domain_nn_calculator_cart<dim>::setParameters(proc_box);
    }

    bool refine(DLB & dlb)
    {
        // if the DLB heuristic to use is the "Unbalance Threshold" get unbalance percentage
        if (dlb.getHeurisitc() == DLB::Heuristic::UNBALANCE_THRLD)
        {
            float unbalance = dist.getUnbalance();
            dlb.setUnbalance(unbalance);
            if (v_cl.getProcessUnitID() == 0)
            {
                std::cout << std::setprecision(3) << unbalance << "\n";
            }
        }

        if (dlb.rebalanceNeeded())
        {
            refine(dlb.getNTimeStepSinceDLB());

            return true;
        }
        return false;
    }

//  size_t n_step = 0;

    float getUnbalance()
    {
        return dist.getUnbalance();
    }

    size_t getProcessorLoad()
    {
        return dist.getProcessorLoad();
    }

    inline void getSubSubDomainPosition(size_t id, T (&pos)[dim])
    {
        dist.getSubSubDomainPosition(id, pos);
    }

    //TODO fix in Parmetis distribution to get only the right amount of vertices
    size_t getNSubSubDomains()
    {
        return dist.getNSubSubDomains();
    }

    inline void setSubSubDomainComputationCost(size_t id, size_t weight)
    {
        dist.setComputationCost(id, weight);
    }

    inline size_t getSubSubDomainComputationCost(size_t id)
    {
        return dist.getSubSubDomainComputationCost(id);
    }

    size_t subSize()
    {
        return dist.subSize();
    }

    size_t getNSubDomain()
    {
        return sub_domains.size();
    }

    SpaceBox<dim, T> getSubDomain(size_t lc)
    {
        // Create a space box
        SpaceBox<dim, T> sp;

        // fill the space box

        for (size_t k = 0; k < dim; k++)
        {
            // create the SpaceBox Low and High
            sp.setLow(k, sub_domains.template get<Box::p1>(lc)[k]);
            sp.setHigh(k, sub_domains.template get<Box::p2>(lc)[k]);
        }

        return sp;
    }

    SpaceBox<dim, T> getSubDomainWithGhost(size_t lc)
    {
        // Create a space box
        SpaceBox<dim, T> sp = sub_domains.get(lc);

        // enlarge with ghost
        sp.enlarge(ghost);

        return sp;
    }

    const ::Box<dim,T> & getDomain() const
    {
        return domain;
    }

    openfpm::vector<SpaceBox<dim, T>> getSubDomains() const
    {
        return sub_domains;
    }

    openfpm::vector<openfpm::vector<SpaceBox<dim, T>>> & getSubDomainsGlobal()
    {
        return sub_domains_global;
    }

    template<typename Mem> bool isLocal(const encapc<1, Point<dim, T>, Mem> p) const
    {
        return processorID<Mem>(p) == v_cl.getProcessUnitID();
    }

    bool isLocal(const T (&pos)[dim]) const
    {
        return processorID(pos) == v_cl.getProcessUnitID();
    }

    bool isLocal(const Point<dim,T> & pos) const
    {
        return processorID(pos) == v_cl.getProcessUnitID();
    }

    template<typename Mem> bool isLocalBC(const encapc<1, Point<dim,T>, Mem> p, const size_t (& bc)[dim]) const
    {
        Point<dim,T> pt = p;

        for (size_t i = 0 ; i < dim ; i++)
        {
            if (bc[i] == PERIODIC)
                pt.get(i) = openfpm::math::periodic_l(p.template get<0>()[i],domain.getHigh(i),domain.getLow(i));
        }

        return processorID<Mem>(pt) == v_cl.getProcessUnitID();
    }

    openfpm::vector<size_t> & getDomainCells()
    {
        return domain_nn_calculator_cart<dim>::getDomainCells();
    }

    openfpm::vector<size_t> & getCRSDomainCells()
    {
        return domain_nn_calculator_cart<dim>::getCRSDomainCells();
    }

    void setNNParameters(grid_key_dx<dim> & shift, grid_sm<dim,void> & gs)
    {
        domain_nn_calculator_cart<dim>::setNNParameters(loc_box, shift, gs);
    }

    openfpm::vector<subsub_lin<dim>> & getCRSAnomDomainCells()
    {
        return domain_nn_calculator_cart<dim>::getCRSAnomDomainCells();
    }

    bool isLocalBC(const T (&p)[dim], const size_t (& bc)[dim]) const
    {
        Point<dim,T> pt = p;

        for (size_t i = 0 ; i < dim ; i++)
        {
            if (bc[i] == PERIODIC)
                pt.get(i) = openfpm::math::periodic_l(p[i],domain.getHigh(i),domain.getLow(i));
        }

        return processorID(pt) == v_cl.getProcessUnitID();
    }


    ::Box<dim, T> & getProcessorBounds()
    {
        return bbox;
    }


    const Ghost<dim,T> & getGhost() const
    {
        return ghost;
    }

    const grid_sm<dim,void> getGrid()
    {
        return gr;
    }

    const grid_sm<dim,void> getDistGrid()
    {
        return gr_dist;
    }


    bool write(std::string output) const
    {
        VTKWriter<openfpm::vector<::SpaceBox<dim, T>>, VECTOR_BOX> vtk_box1;
        vtk_box1.add(sub_domains);
        vtk_box1.write(output + std::string("subdomains_") + std::to_string(v_cl.getProcessUnitID()) + std::string(".vtk"));

        nn_prcs<dim, T>::write(output);
        ie_ghost<dim, T>::write(output, v_cl.getProcessUnitID());
        ie_loc_ghost<dim, T>::write(output, v_cl.getProcessUnitID());

        return true;
    }

    Vcluster & getVC() const
    {
#ifdef SE_CLASS2
        check_valid(this,8);
#endif
        return v_cl;
    }

    bool check_consistency()
    {
        if (ie_loc_ghost<dim, T>::check_consistency(getNSubDomain()) == false)
            return false;

        return true;
    }

    void debugPrint()
    {
        std::cout << "Subdomains\n";
        for (size_t p = 0; p < sub_domains.size(); p++)
        {
            std::cout << ::SpaceBox<dim, T>(sub_domains.get(p)).toString() << "\n";
        }

        std::cout << "External ghost box\n";

        for (size_t p = 0; p<nn_prcs < dim, T>::getNNProcessors(); p++)
        {
            for (size_t i = 0; i<ie_ghost < dim, T>::getProcessorNEGhost(p); i++)
            {
                std::cout << ie_ghost<dim, T>::getProcessorEGhostBox(p, i).toString() << "   prc=" << nn_prcs<dim, T>::IDtoProc(p) << "   id=" << ie_ghost<dim, T>::getProcessorEGhostId(p, i) << "\n";
            }
        }

        std::cout << "Internal ghost box\n";

        for (size_t p = 0; p<nn_prcs < dim, T>::getNNProcessors(); p++)
        {
            for (size_t i = 0; i<ie_ghost < dim, T>::getProcessorNIGhost(p); i++)
            {
                std::cout << ie_ghost<dim, T>::getProcessorIGhostBox(p, i).toString() << "   prc=" << nn_prcs<dim, T>::IDtoProc(p) << "   id=" << ie_ghost<dim, T>::getProcessorIGhostId(p, i) << "\n";
            }
        }
    }

    bool is_equal(CartDecomposition<dim,T,Memory> & cart)
    {
        if (static_cast<ie_loc_ghost<dim,T>*>(this)->is_equal(static_cast<ie_loc_ghost<dim,T>&>(cart)) == false)
            return false;

        if (static_cast<nn_prcs<dim,T>*>(this)->is_equal(static_cast<nn_prcs<dim,T>&>(cart)) == false)
            return false;

        if (static_cast<ie_ghost<dim,T>*>(this)->is_equal(static_cast<ie_ghost<dim,T>&>(cart)) == false)
            return false;

        if (sub_domains != cart.sub_domains)
            return false;

        if (box_nn_processor != cart.box_nn_processor)
            return false;

        if (fine_s != cart.fine_s)
            return false;

        if (gr != cart.gr)
            return false;

        if (cd != cart.cd)
            return false;

        if (domain != cart.domain)
            return false;

        if (meta_compare<T[dim]>::meta_compare_f(cart.spacing,spacing) == false)
            return false;

        if (ghost != cart.ghost)
            return false;

        return true;
    }

    bool is_equal_ng(CartDecomposition<dim,T,Memory> & cart)
    {
        if (static_cast<ie_loc_ghost<dim,T>*>(this)->is_equal_ng(static_cast<ie_loc_ghost<dim,T>&>(cart)) == false)
            return false;

        if (static_cast<nn_prcs<dim,T>*>(this)->is_equal(static_cast<nn_prcs<dim,T>&>(cart)) == false)
            return false;

        if (static_cast<ie_ghost<dim,T>*>(this)->is_equal_ng(static_cast<ie_ghost<dim,T>&>(cart)) == false)
            return false;

        if (sub_domains != cart.sub_domains)
            return false;

        if (box_nn_processor != cart.box_nn_processor)
            return false;

        if (fine_s != cart.fine_s)
            return false;

        if (gr != cart.gr)
            return false;

        if (cd != cart.cd)
            return false;

        if (domain != cart.domain)
            return false;

        if (meta_compare<T[dim]>::meta_compare_f(cart.spacing,spacing) == false)
            return false;

        return true;
    }

    Distribution & getDistribution()
    {
        return dist;
    }


    inline void addComputationCost(size_t gid, size_t i)
    {
        size_t c = dist.getSubSubDomainComputationCost(gid);

        dist.setComputationCost(gid, c + i);
    }

    size_t get_ndec()
    {
        return dist.get_ndec();
    }

    const CellDecomposer_sm<dim, T, shift<dim,T>> & getCellDecomposer()
    {
        return cd;
    }

    friend extended_type;

};


#endif