SpaceDistributionWeight.hpp

/*
 * SpaceDistributionWeight.hpp
 *
 *  Created on: Apr 5, 2017
 *      Author: i-bird
 */
 
#ifndef SRC_DECOMPOSITION_DISTRIBUTION_SPACEDISTRIBUTIONWEIGHT_HPP_
#define SRC_DECOMPOSITION_DISTRIBUTION_SPACEDISTRIBUTIONWEIGHT_HPP_
 
#if 0
 
#include "util/mathutil.hpp"
#include "NN/CellList/CellDecomposer.hpp"
 
template<unsigned int dim, typename T>
class SpaceDistributionWeight
{
    Vcluster & v_cl;
 
    grid_sm<dim, void> gr;
 
    Box<dim, T> domain;
 
    Graph_CSR<nm_v, nm_e> gp;
 
 
    openfpm::vector<unsigned int> dec;
 
public:
 
    SpaceDistributionWeight(Vcluster & v_cl)
    :v_cl(v_cl)
    {
    }
 
    SpaceDistributionWeight(SpaceDistributionWeight<dim,T> && pm)
    {
        this->operator=(pm);
    }
 
    void createCartGraph(grid_sm<dim, void> & grid, Box<dim, T> dom)
    {
        size_t bc[dim];
 
        for (size_t i = 0 ; i < dim ; i++)
            bc[i] = NON_PERIODIC;
 
        // Set grid and domain
        gr = grid;
        domain = dom;
 
        // Create a cartesian grid graph
        CartesianGraphFactory<dim, Graph_CSR<nm_v, nm_e>> g_factory_part;
        gp = g_factory_part.template construct<NO_EDGE, nm_v::id, T, dim - 1, 0>(gr.getSize(), domain, bc);
 
        // Init to 0.0 axis z (to fix in graphFactory)
        if (dim < 3)
        {
            for (size_t i = 0; i < gp.getNVertex(); i++)
                gp.vertex(i).template get<nm_v::x>()[2] = 0.0;
        }
        for (size_t i = 0; i < gp.getNVertex(); i++)
            gp.vertex(i).template get<nm_v::global_id>() = i;
 
    }
 
    Graph_CSR<nm_v, nm_e> & getGraph()
    {
        return gp;
    }
 
    void decompose()
    {
        // collect all the weight on the graph index + weight on master
 
        openfpm::vector<id_w> id_w;
        v_cl.SGather(id_w,0);
 
        // Calculate the total weight
        v_cl.sum();
        v_cl.execute();
 
        tot /= v_cl.getProcessingUnits();
 
        // only master compute
        if (v_cl.getProcessUnitID() == 0)
        {
            for (size_t i = 0 ; i < id_w.size() ; i++)
                gp.template vertex_p<nm_v::computation>(id_w.get(i).id) = id_w.get(i).w;
 
            // Get the maximum along dimensions and take the smallest n number
            // such that 2^n < m. n it will be order of the hilbert curve
 
            size_t max = 0;
 
            for (size_t i = 0; i < dim ; i++)
            {
                if (max < gr.size(i))
                    max = gr.size(i);
            }
 
            // Get the order of the hilbert-curve
            size_t order = openfpm::math::log2_64(max);
            if (1ul << order < max)
                order += 1;
 
            size_t n = 1 << order;
 
            // Create the CellDecomoser
 
            CellDecomposer_sm<dim,T> cd_sm;
            cd_sm.setDimensions(domain, gr.getSize(), 0);
 
            // create the hilbert curve
 
            //hilbert curve iterator
            grid_key_dx_iterator_hilbert<dim> h_it(order);
 
            T spacing[dim];
 
            // Calculate the hilbert curve spacing
            for (size_t i = 0 ; i < dim ; i++)
                spacing[i] = (domain.getHigh(i) - domain.getLow(i)) / n;
 
            // Small grid to detect already assigned sub-sub-domains
            grid_cpu<dim,aggregate<long int>> g(gr.getSize());
            g.setMemory();
 
            // Reset the grid to -1
            grid_key_dx_iterator<dim> it(gr);
            while (it.isNext())
            {
                auto key = it.get();
 
                g.template get<0>(key) = -1;
 
                ++it;
            }
 
            // Go along the hilbert-curve and divide the space
            size_t proc_d = 0;
            size_t ele_d = 0;
            while (h_it.isNext())
            {
              auto key = h_it.get();
 
              // Point p
              Point<dim,T> p;
 
              for (size_t i = 0 ; i < dim ; i++)
                  p.get(i) = key.get(i) * spacing[i] + spacing[i] / 2;
 
              grid_key_dx<dim> sp = cd_sm.getCellGrid(p);
 
              if (g.template get<0>(sp) == -1)
              {
                  g.template get<0>(sp) = proc_d;
                  ele_d++;
 
                  if (ele_d >= avg)
                      proc_d++;
              }
 
              ++h_it;
            }
 
            // Fill from the grid to the graph
 
            // Reset the grid to -1
            grid_key_dx_iterator<dim> it2(gr);
            while (it2.isNext())
            {
                auto key = it2.get();
 
                gp.template vertex_p<nm_v::proc_id>(gr.LinId(key)) = g.template get<0>(key);
 
                id_w.get(g.template get<1>(key)) = g.template get<1>(key);
 
                ++it2;
            }
        }
 
        // Broadcast the result to all the processors.
 
        v_cl.Bcast(id_w);
        v_cl.execute();
 
        //
 
        return;
    }
 
    void refine()
    {
        decompose();
    }
 
    float getUnbalance()
    {
        return gr.size() % v_cl.getProcessingUnits();
    }
 
    void getSubSubDomainPosition(size_t id, T (&pos)[dim])
    {
#ifdef SE_CLASS1
        if (id >= gp.getNVertex())
            std::cerr << __FILE__ << ":" << __LINE__ << "Such vertex doesn't exist (id = " << id << ", " << "total size = " << gp.getNVertex() << ")\n";
#endif
 
        // Copy the geometrical informations inside the pos vector
        pos[0] = gp.vertex(id).template get<nm_v::x>()[0];
        pos[1] = gp.vertex(id).template get<nm_v::x>()[1];
        if (dim == 3)
            pos[2] = gp.vertex(id).template get<nm_v::x>()[2];
    }
 
    inline void setComputationCost(size_t id, size_t weight)
    {
        std::cout << __FILE__ << ":" << __LINE__ << " You are trying to set the computation cost on a fixed decomposition, this operation has no effect" << std::endl;
    }
 
    bool weightsAreUsed()
    {
        return false;
    }
 
    size_t getSubSubDomainComputationCost(size_t id)
    {
        return 1.0;
    }
 
    size_t getProcessorLoad()
    {
        // Get the number of processing units
        size_t Np = v_cl.getProcessingUnits();
 
        // Calculate the best number of sub-domains for each
        // processor
        size_t N_tot = gr.size();
        size_t N_best_each = N_tot / Np;
        size_t N_rest = N_tot % Np;
 
        if (v_cl.getProcessUnitID() < N_rest)
            N_best_each += 1;
 
        return N_best_each;
    }
 
    void setMigrationCost(size_t id, size_t migration)
    {
    }
 
    void setCommunicationCost(size_t v_id, size_t e, size_t communication)
    {
    }
 
    size_t getNSubSubDomains()
    {
        return gp.getNVertex();
    }
 
    size_t getNSubSubDomainNeighbors(size_t id)
    {
        return gp.getNChilds(id);
    }
 
    void write(const std::string & file)
    {
        VTKWriter<Graph_CSR<nm_v, nm_e>, VTK_GRAPH> gv2(gp);
        gv2.write(std::to_string(v_cl.getProcessUnitID()) + "_" + file + ".vtk");
    }
 
    const SpaceDistribution<dim,T> & operator=(const SpaceDistribution<dim,T> & dist)
    {
        gr = dist.gr;
        domain = dist.domain;
        gp = dist.gp;
 
        return *this;
    }
 
    const SpaceDistribution<dim,T> & operator=(SpaceDistribution<dim,T> && dist)
    {
        v_cl = dist.v_cl;
        gr = dist.gr;
        domain = dist.domain;
        gp.swap(dist.gp);
 
        return *this;
    }
 
    size_t get_ndec()
    {
        return 0;
    }
};
 
#endif
 
#endif /* SRC_DECOMPOSITION_DISTRIBUTION_SPACEDISTRIBUTIONWEIGHT_HPP_ */