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"
#include "cuda/CartDecomposition_gpu.cuh"
#include "Domain_icells_cart.hpp"

#define CARTDEC_ERROR 2000lu

enum dec_options
{
    DEC_NONE = 0,
    DEC_SKIP_ICELL = 1
};

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, template <typename> class layout_base, typename Distribution>
class CartDecomposition: public ie_loc_ghost<dim, T,layout_base, Memory>,
                         public nn_prcs<dim, T,layout_base,Memory>,
                         public ie_ghost<dim,T,Memory,layout_base>,
                         public domain_nn_calculator_cart<dim>,
                         public domain_icell_calculator<dim,T,layout_base,Memory>
{
public:

    typedef T domain_type;

    typedef SpaceBox<dim, T> Box;

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

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

    typedef ie_loc_ghost<dim, T,layout_base, Memory> ie_loc_ghost_type;

    typedef nn_prcs<dim, T,layout_base,Memory> nn_prcs_type;

    typedef ie_ghost<dim,T,Memory,layout_base> ie_ghost_type;

protected:

    bool host_dev_transfer = false;

    bool commCostSet = false;

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

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

    mutable openfpm::vector<Box_map<dim, T>,Memory,layout_base> sub_domains_global;

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

    CellList<dim,T,Mem_fast<Memory,int>,shift<dim,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;
    }

    void collect_all_sub_domains(openfpm::vector<Box_map<dim,T>,Memory,layout_base> & sub_domains_global)
    {
#ifdef SE_CLASS2
        check_valid(this,8);
#endif

        sub_domains_global.clear();
        openfpm::vector<Box_map<dim,T>,Memory,layout_base> bm;

        for (size_t i = 0 ; i < sub_domains.size() ; i++)
        {
            bm.add();

            bm.template get<0>(bm.size()-1) = ::SpaceBox<dim,T>(sub_domains.get(i));
            bm.template get<1>(bm.size()-1) = v_cl.rank();
        }

        v_cl.SGather<decltype(bm),decltype(sub_domains_global),layout_base>(bm,sub_domains_global,0);

        size_t size = sub_domains_global.size();

        v_cl.max(size);
        v_cl.execute();

        sub_domains_global.resize(size);

        v_cl.Bcast(sub_domains_global,0);
        v_cl.execute();
    }

public:

    void initialize_fine_s(const ::Box<dim,T> & domain)
    {
        fine_s.clear();
        size_t div_g[dim];

        // We reduce the size of the cells by a factor 8 in 3d 4 in 2d
        for (size_t i = 0 ; i < dim ; i++)
        {div_g[i] = (gr.size(i) == 1)?1:gr.size(i)/2;}

        fine_s.Initialize(domain,div_g);
    }

    void construct_fine_s()
    {
        collect_all_sub_domains(sub_domains_global);

        // now draw all sub-domains in fine-s

        for (size_t i = 0 ; i < sub_domains_global.size() ; i++)
        {

            // get the cells this box span
            const grid_key_dx<dim> p1 = fine_s.getCellGrid_me(sub_domains_global.template get<0>(i).getP1());
            const grid_key_dx<dim> p2 = fine_s.getCellGrid_pe(sub_domains_global.template get<0>(i).getP2());

            // Get the grid and the sub-iterator
            auto & gi = fine_s.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();
                fine_s.addCell(gi.LinId(key),i);

                ++g_sub;
            }
        }

        host_dev_transfer = false;
    }

    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
        initialize_fine_s(domain);

        // Optimize the decomposition creating bigger spaces
        // And reducing Ghost over-stress
        dec_optimizer<dim, Graph_CSR<nm_v<dim>, 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 or merge sub-sub-domain
        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,layout_base,Memory>::create(box_nn_processor, sub_domains);
        nn_prcs<dim,T,layout_base,Memory>::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)


        construct_fine_s();

        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)
        {
            // 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,Memory,layout_base>::Initialize_geo_cell(bound,div);

            // Initialize shift vectors
            ie_ghost<dim,T,Memory,layout_base>::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,Memory,layout_base>::create_box_nn_processor_ext(v_cl, ghost, sub_domains, box_nn_processor, *this);
        ie_ghost<dim, T,Memory,layout_base>::create_box_nn_processor_int(v_cl, ghost, sub_domains, box_nn_processor, *this);

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

        // Ghost box information must be re-offloaded
        host_dev_transfer = false;
        ie_ghost<dim, T,Memory,layout_base>::reset_host_dev_transfer();
    }


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, layout_base,Memory>(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,layout_base,Distribution> & cart)
    :nn_prcs<dim,T,layout_base,Memory>(cart.v_cl),v_cl(cart.v_cl),dist(v_cl),ref_cnt(0)
    {
        this->operator=(cart);
    }

    CartDecomposition(CartDecomposition<dim,T,Memory,layout_base,Distribution> && cart)
    :nn_prcs<dim,T,layout_base,Memory>(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:
        template<typename encap_type> inline static size_t id(const encap_type & p, size_t b_id)
        {
            return p.template get<proc_>();
        }
    };

    class lc_processor_id
    {
    public:
        template<typename encap_type> inline static size_t id(const encap_type & p, size_t b_id)
        {
            return p.template get<lc_proc_>();
        }
    };

    class shift_id
    {
    public:
        template<typename encap_type> inline static size_t id(const encap_type & p, size_t b_id)
        {
            return p.template get<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,layout_base,Distribution> duplicate(const Ghost<dim,T> & g) const
    {
        CartDecomposition<dim,T,Memory,layout_base,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.magn[i] = magn[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,layout_base,Memory> &>(cart)).create(box_nn_processor, sub_domains);
        (static_cast<nn_prcs<dim,T,layout_base,Memory> &>(cart)).applyBC(domain,ghost,bc);

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

        cart.collect_all_sub_domains(cart.sub_domains_global);

        return cart;
    }

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

        (static_cast<ie_loc_ghost<dim,T, layout_base,Memory>*>(&cart))->operator=(static_cast<ie_loc_ghost<dim,T,layout_base,Memory>>(*this));
        (static_cast<nn_prcs<dim,T,layout_base,Memory>*>(&cart))->operator=(static_cast<nn_prcs<dim,T,layout_base,Memory>>(*this));
        (static_cast<ie_ghost<dim,T,Memory,layout_base>*>(&cart))->operator=(static_cast<ie_ghost<dim,T,Memory,layout_base>>(*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;
        cart.sub_domains_global = sub_domains_global;
        for (size_t i = 0 ; i < dim ; i++)
        {
            cart.spacing[i] = spacing[i];
            cart.magn[i] = magn[i];
        };

        cart.ghost = ghost;

        cart.bbox = bbox;

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

        return cart;
    }

    template<typename Memory2, template <typename> class layout_base2>
    CartDecomposition<dim,T,Memory2,layout_base2,Distribution> duplicate_convert() const
    {
        CartDecomposition<dim,T> cart(v_cl);

        (static_cast<ie_loc_ghost<dim,T,layout_base2,Memory2>*>(&cart))->operator=(static_cast<ie_loc_ghost<dim,T,layout_base,Memory>>(*this));
        (static_cast<nn_prcs<dim,T,layout_base2,Memory2>*>(&cart))->operator=(static_cast<nn_prcs<dim,T,layout_base,Memory>>(*this));
        ie_ghost<dim,T,Memory,layout_base> * ptr = static_cast<ie_ghost<dim,T,Memory,layout_base> *>((CartDecomposition<dim,T,Memory,layout_base,Distribution> *)this);
        (static_cast<ie_ghost<dim,T,Memory2,layout_base2>*>(&cart))->operator=(ptr->template duplicate<Memory2,layout_base2>());

        cart.private_get_sub_domains() = sub_domains;
        cart.private_get_box_nn_processor() = box_nn_processor;
        cart.private_get_fine_s() = fine_s;
        cart.private_get_gr() = gr;
        cart.private_get_gr_dist() = gr_dist;
        cart.private_get_dist() = dist;
        cart.private_get_commCostSet() = commCostSet;
        cart.private_get_cd() = cd;
        cart.private_get_domain() = domain;
        cart.private_get_sub_domains_global() = sub_domains_global;
        for (size_t i = 0 ; i < dim ; i++)
        {cart.private_get_spacing(i) = spacing[i];};

        cart.private_get_ghost() = ghost;

        cart.private_get_bbox() = bbox;

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

        return cart;
    }

    CartDecomposition<dim,T,Memory, layout_base, Distribution> & operator=(const CartDecomposition & cart)
    {
        static_cast<ie_loc_ghost<dim,T,layout_base,Memory>*>(this)->operator=(static_cast<ie_loc_ghost<dim,T,layout_base,Memory>>(cart));
        static_cast<nn_prcs<dim,T,layout_base,Memory>*>(this)->operator=(static_cast<nn_prcs<dim,T,layout_base,Memory>>(cart));
        static_cast<ie_ghost<dim,T,Memory,layout_base>*>(this)->operator=(static_cast<ie_ghost<dim,T,Memory,layout_base>>(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;
        sub_domains_global = cart.sub_domains_global;

        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;
    }

    template<typename CartDecomposition2>
    CartDecomposition<dim,T,Memory, layout_base, Distribution> & operator=(const CartDecomposition2 & cart)
    {
        static_cast<ie_loc_ghost<dim,T,layout_base,Memory>*>(this)->operator=(static_cast<typename CartDecomposition2::ie_loc_ghost_type>(cart));
        static_cast<nn_prcs<dim,T,layout_base,Memory>*>(this)->operator=(static_cast<typename CartDecomposition2::nn_prcs_type>(cart));
        static_cast<ie_ghost<dim,T,Memory,layout_base>*>(this)->operator=(static_cast<typename CartDecomposition2::ie_ghost_type>(cart));

        sub_domains = cart.private_get_sub_domains();
        box_nn_processor = cart.private_get_box_nn_processor();
        fine_s = cart.private_get_fine_s();
        gr = cart.private_get_gr();
        gr_dist = cart.private_get_gr_dist();
        dist = cart.private_get_dist();
        commCostSet = cart.private_get_commCostSet();
        cd = cart.private_get_cd();
        domain = cart.private_get_domain();
        sub_domains_global = cart.private_get_sub_domains_global();

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

        ghost = cart.private_get_ghost();

        bbox = cart.private_get_bbox();

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

        return *this;
    }

    CartDecomposition<dim,T,Memory,layout_base, Distribution> & operator=(CartDecomposition && cart)
    {
        static_cast<ie_loc_ghost<dim,T,layout_base,Memory>*>(this)->operator=(static_cast<ie_loc_ghost<dim,T,layout_base,Memory>>(cart));
        static_cast<nn_prcs<dim,T,layout_base,Memory>*>(this)->operator=(static_cast<nn_prcs<dim,T,layout_base,Memory>>(cart));
        static_cast<ie_ghost<dim,T,Memory,layout_base>*>(this)->operator=(static_cast<ie_ghost<dim,T,Memory,layout_base>>(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;
        gr_dist = cart.gr_dist;
        dist = cart.dist;

        domain = cart.domain;
        sub_domains_global.swap(cart.sub_domains_global);

        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 processorID_impl(p,fine_s,sub_domains_global,getDomain(),bc);
    }

    size_t inline processorID(const Point<dim,T> &p) const
    {
        return processorID_impl(p,fine_s,sub_domains_global,getDomain(),bc);
    }

    size_t inline processorID(const T (&p)[dim]) const
    {
        return processorID_impl(p,fine_s,sub_domains_global,getDomain(),bc);
    }

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


        return processorID_impl(pt,fine_s,sub_domains_global,getDomain(),bc);
    }

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

        // Get the number of elements in the cell

        return processorID_impl(pt,fine_s,sub_domains_global,getDomain(),bc);
    }

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

        return processorID_impl(pt,fine_s,sub_domains_global,getDomain(),bc);
    }

    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 > (unsigned int long)(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,layout_base,Memory>::reset();
        ie_ghost<dim,T,Memory,layout_base>::reset();
        ie_loc_ghost<dim, T,layout_base,Memory>::reset();
    }

    void decompose(dec_options opt = dec_options::DEC_NONE)
    {
        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);

        if (opt != dec_options::DEC_SKIP_ICELL)
        {

            domain_icell_calculator<dim,T,layout_base,Memory>
			::CalculateInternalCells(v_cl,
                                 ie_ghost<dim, T,Memory,layout_base>::private_get_vb_int_box(),
                                 sub_domains,
                                 this->getProcessorBounds(),
                                 this->getGhost().getRcut(),
                                 this->getGhost());
        }
    }

    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;
    }

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

    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();
    }

    bool isLocalBC(const Point<dim,T> & 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[i],domain.getHigh(i),domain.getLow(i));
        }

        return processorID(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>,Memory,layout_base>, 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,layout_base,Memory>::write(output);
        ie_ghost<dim,T,Memory,layout_base>::write(output, v_cl.getProcessUnitID());
        ie_loc_ghost<dim, T,layout_base,Memory>::write(output, v_cl.getProcessUnitID());

        return true;
    }

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

    void free_geo_cell()
    {
        ie_ghost<dim,T,Memory,layout_base>::free_geo_cell();
    }

    void free_fines()
    {
        fine_s.clear();
        fine_s.destroy();
    }

    bool check_consistency()
    {
        if (ie_loc_ghost<dim, T,layout_base,Memory>::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 << "Subdomains global\n";
        for (size_t p = 0; p < sub_domains_global.size(); p++)
        {
            std::cout << ::SpaceBox<dim, T>(sub_domains_global.template get<0>(p)).toString() << " proc:" << sub_domains_global.template get<1>(p) << "\n";
        }

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

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

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

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

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

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

        if (static_cast<ie_ghost<dim,T,Memory,layout_base>*>(this)->is_equal(static_cast<ie_ghost<dim,T,Memory,layout_base>&>(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,layout_base,Memory>*>(this)->is_equal_ng(static_cast<ie_loc_ghost<dim,T,layout_base,Memory>&>(cart)) == false)
            return false;

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

        if (static_cast<ie_ghost<dim,T,Memory,layout_base>*>(this)->is_equal_ng(static_cast<ie_ghost<dim,T,Memory,layout_base>&>(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;
    }


    CartDecomposition_gpu<dim,T,Memory,layout_base> toKernel()
    {
        if (host_dev_transfer == false)
        {
            fine_s.hostToDevice();
            sub_domains_global.template hostToDevice<0,1>();
            host_dev_transfer = true;
        }

        int bc_[dim];

        for (int i = 0 ; i < dim ; i++) {bc_[i] = this->periodicity(i);}

        CartDecomposition_gpu<dim,T,Memory,layout_base> cdg(fine_s.toKernel(),
                                                            ie_ghost<dim,T,Memory,layout_base>::toKernel(),
                                                sub_domains_global.toKernel(),
                                                getDomain(),
                                                bc_);

        return cdg;
    }

    friend extended_type;

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

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

    openfpm::vector<openfpm::vector<long unsigned int> > & private_get_box_nn_processor()
    {
        return box_nn_processor;
    }

    const openfpm::vector<openfpm::vector<long unsigned int> > & private_get_box_nn_processor() const
    {
        return box_nn_processor;
    }

    CellList<dim,T,Mem_fast<Memory,int>,shift<dim,T>> & private_get_fine_s()
    {
        return fine_s;
    }

    const CellList<dim,T,Mem_fast<Memory,int>,shift<dim,T>> & private_get_fine_s() const
    {
        return fine_s;
    }

    grid_sm<dim, void> & private_get_gr()
    {
        return gr;
    }

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

    grid_sm<dim, void> & private_get_gr_dist()
    {
        return gr_dist;
    }

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

    Distribution & private_get_dist()
    {
        return dist;
    }

    const Distribution & private_get_dist() const
    {
        return dist;
    }

    bool & private_get_commCostSet()
    {
        return commCostSet;
    }

    const bool & private_get_commCostSet() const
    {
        return commCostSet;
    }

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

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

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

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

    openfpm::vector<Box_map<dim, T>,Memory,layout_base> & private_get_sub_domains_global()
    {
        return sub_domains_global;
    }

    const openfpm::vector<Box_map<dim, T>,Memory,layout_base> & private_get_sub_domains_global() const
    {
        return sub_domains_global;
    }

    T & private_get_spacing(int i)
    {
        return spacing[i];
    }

    const T & private_get_spacing(int i) const
    {
        return spacing[i];
    }

    T & private_get_magn(int i)
    {
        return spacing[i];
    }

    const T & private_get_magn(int i) const
    {
        return spacing[i];
    }


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

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

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

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

    size_t & private_get_bc(int i)
    {
        return bc[i];
    }

    const size_t & private_get_bc(int i) const
    {
        return bc[i];
    }
};


#endif