CartDecomposition_ext.hpp

/*
 * CartDecomposition_ext.hpp
 *
 *  Created on: Mar 6, 2016
 *      Author: i-bird
 */

#ifndef SRC_DECOMPOSITION_CARTDECOMPOSITION_EXT_HPP_
#define SRC_DECOMPOSITION_CARTDECOMPOSITION_EXT_HPP_

#include "memory/HeapMemory.hpp"
#include "Decomposition/Distribution/ParMetisDistribution.hpp"
#include "Space/Ghost.hpp"
#include "Decomposition/nn_processor.hpp"

template<unsigned int dim, typename T, typename Memory = HeapMemory, typename Distribution = ParMetisDistribution<dim, T>>
class CartDecomposition;

template<unsigned int dim, typename T, typename Memory = HeapMemory, typename Distribution = ParMetisDistribution<dim, T>>
class CartDecomposition_ext: public CartDecomposition<dim,T,Memory,Distribution>
{
private:

    void extend_subdomains(const CartDecomposition<dim,T,Memory,Distribution> & dec, const ::Box<dim,T> & ext_dom)
    {
        // Box
        typedef ::Box<dim,T> b;

        this->bbox.zero();

        // Extend sub-domains
        for (size_t i = 0 ; i < dec.sub_domains.size() ; i++)
        {
            ::Box<dim,T> box;

            // Calculate the extended box
            for (size_t j = 0 ; j < dim ; j++)
            {
                if (dec.sub_domains.template get<b::p1>(i)[j] == dec.domain.getLow(j))
                    box.setLow(j,ext_dom.getLow(j));
                else
                    box.setLow(j,dec.sub_domains.template get<b::p1>(i)[j]);

                if (dec.sub_domains.template get<b::p2>(i)[j] == dec.domain.getHigh(j))
                    box.setHigh(j,ext_dom.getHigh(j));
                else
                    box.setHigh(j,dec.sub_domains.template get<b::p2>(i)[j]);
            }

            // add the subdomain
            this->sub_domains.add(box);

            // Calculate the bound box
            this->bbox.enclose(box);
        }
    }

    void extend_fines(const CartDecomposition<dim,T,Memory,Distribution> & dec)
    {
        // Extension, first we calculate the extensions of the new domain compared
        // to the old one in cell units (each cell unit is a sub-sub-domain)
        ::Box<dim,size_t> ext;
        // Extension of the new fines structure
        ::Box<dim,size_t> n_fines_ext;
        // Extension of the old fines structure
        ::Box<dim,size_t> o_fines_ext;

        size_t sz_new[dim];
        size_t sz_old[dim];

        for (size_t i = 0; i < dim ; i++)
        {
            size_t p1 = (dec.domain.getLow(i) - dec.domain.getLow(i)) / dec.cd.getCellBox().getHigh(i) + 1;
            size_t p2 = (dec.domain.getLow(i) - dec.domain.getLow(i)) / dec.cd.getCellBox().getHigh(i) + 1;

            ext.setLow(i,p1);
            ext.setHigh(i,p2);
            sz_new[i] = p1+p2+dec.cd.getGrid().size(i);
            sz_old[i] = dec.cd.getGrid().size(i);
        }

        grid_sm<dim,void> info_new(sz_new);
        grid_sm<dim,void> info_old(sz_old);

        // resize the new fines
        this->fine_s.resize(info_new.size());

        // we create an iterator that iterate across the full new fines
        grid_key_dx_iterator<dim> fines_t(info_new);

        while (fines_t.isNext())
        {
            auto key = fines_t.get();

            // new_fines is bigger than old_fines structure
            // out of bound key must be adjusted
            // The adjustment produce a natural extension
            // a representation can be seen in the figure of
            // CartDecomposition duplicate function with extended domains

            grid_key_dx<dim> key_old;
            for (size_t i = 0 ; i < dim ; i++)
            {
                key_old.set_d(i,(long int)key.get(i) - ext.getLow(i));
                if (key_old.get(i) < 0)
                    key_old.set_d(i,0);
                else if(key_old.get(i) >= (long int)info_old.size(i) )
                    key_old.set_d(i,info_old.size(i)-1);
            }

            this->fine_s.get(info_new.LinId(key)) = dec.fine_s.get(info_old.LinId(key_old));

            ++fines_t;
        }

        this->gr.setDimensions(sz_new);

        // the new extended CellDecomposer must be consistent with the old cellDecomposer.
        this->cd.setDimensions(dec.cd,ext);
    }

public:

    CartDecomposition_ext(Vcluster & v_cl)
    :CartDecomposition<dim,T,Memory,Distribution>(v_cl)
    {
    }

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

    void setParameters(const CartDecomposition<dim,T,Memory,Distribution> & dec, const Ghost<dim,T> & g, const ::Box<dim,T> & ext_domain)
    {
        this->box_nn_processor = dec.box_nn_processor;

        // Calculate new sub-domains for extended domain
        extend_subdomains(dec,ext_domain);

        // Calculate fine_s structure for the extended domain
        // update the cell decomposer and gr
        extend_fines(dec);

        // Get the old sub-sub-domain grid extension

        this->domain = ext_domain;

        // spacing does not change

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

        this->ghost = g;
        this->dist = dec.dist;

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

        (static_cast<nn_prcs<dim,T> &>(*this)).create(this->box_nn_processor, this->sub_domains);
        (static_cast<nn_prcs<dim,T> &>(*this)).applyBC(ext_domain,g,this->bc);

        this->Initialize_geo_cell_lists();
        this->calculateGhostBoxes();
    }

};



#endif /* SRC_DECOMPOSITION_CARTDECOMPOSITION_EXT_HPP_ */