staggered_dist_grid_util.hpp

/*
 * staggered_util.hpp
 *
 *  Created on: Aug 19, 2015
 *      Author: i-bird
 */

#ifndef SRC_GRID_STAGGERED_DIST_GRID_UTIL_HPP_
#define SRC_GRID_STAGGERED_DIST_GRID_UTIL_HPP_

#include "util/common.hpp"
#include "VTKWriter/VTKWriter.hpp"
#include "util/convert.hpp"


template<typename ele, typename vtk, bool has_attributes=has_attributes<ele>::value>
struct vtk_write
{
    vtk_write(vtk vtk_w, const std::string output, const size_t i)
    {
        vtk_w.write(output + "_" + ele::attributes::name[i] + ".vtk",ele::attributes::name[i]);
    }
};

template<typename ele, typename vtk>
struct vtk_write<ele,vtk,false>
{
    vtk_write(vtk vtk_w, const std::string output, const size_t i)
    {
        vtk_w.write(output + "_" + std::to_string(i) + ".vtk","attr" + std::to_string(i));
    }
};


template<typename T>
struct extends
{
    static inline size_t mul()
    {
        return 1;
    }

    static inline size_t dim()
    {
        return 0;
    }
};

template<typename T,size_t N1>
struct extends<T[N1]>
{
    static inline size_t mul()
    {
        return N1;
    }

    static inline size_t dim()
    {
        return 1;
    }
};

template<typename T,size_t N1,size_t N2>
struct extends<T[N1][N2]>
{
    static inline size_t mul()
    {
        return N1 * N2;
    }

    static inline size_t dim()
    {
        return 2;
    }
};

template<typename T,size_t N1,size_t N2,size_t N3>
struct extends<T[N1][N2][N3]>
{
    static inline size_t mul()
    {
        return N1 * N2 * N3;
    }

    static inline size_t dim()
    {
        return 3;
    }
};

template<typename T,size_t N1,size_t N2,size_t N3,size_t N4>
struct extends<T[N1][N2][N3][N4]>
{
    static inline size_t mul()
    {
        return N1 * N2 * N3 * N4;
    }

    static inline size_t dim()
    {
        return 4;
    }
};

template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5>
struct extends<T[N1][N2][N3][N4][N5]>
{
    static inline size_t mul()
    {
        return N1 * N2 * N3 * N4 * N5;
    }

    static inline size_t dim()
    {
        return 5;
    }
};

template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6>
struct extends<T[N1][N2][N3][N4][N5][N6]>
{
    static inline size_t mul()
    {
        return N1 * N2 * N3 * N4 * N5 * N6;
    }

    static inline size_t dim()
    {
        return 6;
    }
};

template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6, size_t N7>
struct extends<T[N1][N2][N3][N4][N5][N6][N7]>
{
    static inline size_t mul()
    {
        return N1 * N2 * N3 * N4 * N5 * N6 * N7;
    }

    static inline size_t dim()
    {
        return 7;
    }
};

template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6, size_t N7, size_t N8>
struct extends<T[N1][N2][N3][N4][N5][N6][N7][N8]>
{
    static inline size_t mul()
    {
        return N1 * N2 * N3 * N4 * N5 * N6 * N7 * N8;
    }

    static inline size_t dim()
    {
        return 8;
    }
};

template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6, size_t N7, size_t N8, size_t N9>
struct extends<T[N1][N2][N3][N4][N5][N6][N7][N8][N9]>
{
    static inline size_t mul()
    {
        return N1 * N2 * N3 * N4 * N5 * N6 * N7 * N8 * N9;
    }

    static inline size_t dim()
    {
        return 9;
    }
};

template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6, size_t N7, size_t N8, size_t N9, size_t N10>
struct extends<T[N1][N2][N3][N4][N5][N6][N7][N8][N9][N10]>
{
    static inline size_t mul()
    {
        return N1 * N2 * N3 * N4 * N5 * N6 * N7 * N8 * N9 * N10;
    }

    static inline size_t dim()
    {
        return 10;
    }
};


template<typename T>
struct write_stag
{
    template<unsigned int p_val, typename sg, typename v_g> static inline void write(sg & st_g, v_g & vg,size_t lg)
    {
        // Add a grid;
        vg.add();
        size_t k = vg.size() - 1;

        // Get the source and destination grid
        auto & g_src = st_g.get_loc_grid(lg);
        auto & g_dst = vg.get(k);

        // Set dimensions and memory
        g_dst.resize(g_src.getGrid().getSize());

        // copy

        auto it = vg.get(k).getIterator();

        while(it.isNext())
        {
            g_dst.template get<0>(it.get()) = g_src.template get<p_val>(it.get());

            ++it;
        }
    }
};

template<typename T,size_t N1>
struct write_stag<T[N1]>
{
    template<unsigned int p_val, typename sg, typename v_g> static inline void write(sg & st_g, v_g & vg,size_t lg)
    {
        for (size_t i = 0 ; i < N1 ; i++)
        {
            // Add a grid;
            vg.add();
            size_t k = vg.size() - 1;

            // Get the source and destination grid
            auto & g_src = st_g.get_loc_grid(lg);
            auto & g_dst = vg.get(k);

            // Set dimensions and memory
            g_dst.resize(g_src.getGrid().getSize());

            auto it = vg.get(k).getIterator();

            while(it.isNext())
            {
                g_dst.template get<0>(it.get()) = g_src.template get<p_val>(it.get())[i];

                ++it;
            }
        }
    }
};

template<typename T,size_t N1,size_t N2>
struct write_stag<T[N1][N2]>
{
    template<unsigned int p_val, typename sg, typename v_g> static inline void write(sg & st_g, v_g & vg,size_t lg)
    {
        for (size_t i = 0 ; i < N1 ; i++)
        {
            for (size_t j = 0 ; j < N2 ; j++)
            {
                // Add a grid;
                vg.add();
                size_t k = vg.size() - 1;

                // Set dimensions and memory
                vg.get(k).resize(st_g.get_loc_grid(lg).getGrid().getSize());

                // copy
                auto & g_src = st_g.get_loc_grid(lg);
                auto & g_dst = vg.get(k);
                auto it = vg.get(k).getIterator();

                while(it.isNext())
                {
                    g_dst.template get<0>(it.get()) = g_src.template get<p_val>(it.get())[i][j];

                    ++it;
                }
            }
        }
    }
};


template<unsigned int dim, typename v, bool has_pM = has_posMask<v>::value>
class stag_set_position
{
    // within the cell
    openfpm::vector<comb<dim>> (& pos_prp)[boost::fusion::result_of::size<v>::type::value];

public:

    stag_set_position( openfpm::vector<comb<dim>> (& pos_prp)[boost::fusion::result_of::size<v>::type::value])
    :pos_prp(pos_prp)
    {}

    template<typename T>
    void operator()(T& t) const
    {
        // This is the type of the object we have to copy
        typedef typename boost::mpl::at<v,typename boost::mpl::int_<T::value>>::type prop;

        bool val = prop::stag_pos_mask[T::value];

        if (val == false)
            return;

        // Dimension of the object
        size_t dim_prp = extends<prop>::dim();

        // It is a scalar
        if (dim_prp == 0)
        {
            comb<dim> c;
            c.zero();

            // It stay in the center
            pos_prp[T::value].add(c);
        }
        else if (dim_prp == 1)
        {
            // It stay on the object of dimension dim-1 (Negative part)
            for (size_t i = 0 ; i < dim ; i++)
            {
                comb<dim> c;
                c.zero();
                c.value(i) = -1;

                pos_prp[T::value].add(c);
            }
        }
        else if (dim_prp == 2)
        {
            // Create an hypercube object
            HyperCube<dim> hyp;

            // Diagonal part live in
            for (size_t i = 0 ; i < dim ; i++)
            {
                comb<dim> c1 = pos_prp[T::value-1].get(i);
                for (size_t j = 0 ; j < dim ; j++)
                {
                    comb<dim> c2;
                    c2.zero();
                    c2.value(i) = 1;

                    comb<dim> c_res = (c1 + c2) & 0x1;

                    pos_prp[T::value].add(c_res);
                }
            }
        }
        else if (dim_prp > 2)
        {
            std::cerr << __FILE__ << ":" << __LINE__ << " Tensor of rank bigger than 2 are not supported";
        }
    }
};


template<unsigned int dim, typename v>
class stag_set_position<dim,v,false>
{
private:

    // within the cell
    openfpm::vector<comb<dim>> (& pos_prp)[boost::fusion::result_of::size<v>::type::value];


public:

    stag_set_position( openfpm::vector<comb<dim>> (& pos_prp)[boost::fusion::result_of::size<v>::type::value])
    :pos_prp(pos_prp)
    {}

    template<typename T>
    void operator()(T& t) const
    {
        // This is the type of the object we have to copy
        typedef typename boost::mpl::at<v,typename boost::mpl::int_<T::value>>::type prop;

        // Dimension of the object
        size_t dim_prp = extends<prop>::dim();

        // It is a scalar
        if (dim_prp == 0)
        {
            comb<dim> c;
            c.zero();

            // It stay in the center
            pos_prp[T::value].add(c);
        }
        else if (dim_prp == 1)
        {
            // It stay on the object of dimension dim-1 (Negative part)
            for (size_t i = 0 ; i < dim ; i++)
            {
                comb<dim> c;
                c.zero();
                c.getComb()[i] = -1;

                pos_prp[T::value].add(c);
            }
        }
        else if (dim_prp == 2)
        {
            // Diagonal part live in
            for (size_t i = 0 ; i < dim ; i++)
            {
                comb<dim> c1 = pos_prp[T::value-1].get(i);
                for (size_t j = 0 ; j < dim ; j++)
                {
                    comb<dim> c2;
                    c2.zero();
                    c2.getComb()[j] = 1;

                    comb<dim> c_res = (c2 + c1).flip();

                    pos_prp[T::value].add(c_res);
                }
            }
        }
        else if (dim_prp > 2)
        {
            std::cerr << __FILE__ << ":" << __LINE__ << " Tensor of rank bigger than 2 are not supported";
        }
    }
};

template<unsigned int dim, typename st_grid, typename St>
class stag_create_and_add_grid
{

    size_t p_id;

    // staggered grid to write
    st_grid & st_g;

public:

    stag_create_and_add_grid(st_grid & st_g, size_t p_id)
    :p_id(p_id),st_g(st_g)
    {}

    template<unsigned int p_val> void out_normal()
    {
        // property type
        typedef typename boost::mpl::at< typename st_grid::value_type::type , typename boost::mpl::int_<p_val> >::type ele;

        // create an openfpm format object from the property type
        typedef object<typename boost::fusion::vector<ele>> d_object;

        VTKWriter<boost::mpl::pair<grid_cpu<dim, d_object >,St>,VECTOR_GRIDS> vtk_w;

        // Create a vector of grids

        openfpm::vector< grid_cpu<dim, d_object > > vg(st_g.getN_loc_grid());

        // for each domain grid
        for (size_t i = 0 ; i < vg.size() ; i++)
        {
            // Set dimansions and memory
            vg.get(i).resize(st_g.get_loc_grid(i).getGrid().getSize());

            auto & g_src = st_g.get_loc_grid(i);
            auto & g_dst = vg.get(i);

            auto it = vg.get(i).getIterator();

            while(it.isNext())
            {
                object_si_d< decltype(g_src.get_o(it.get())),decltype(g_dst.get_o(it.get())) ,OBJ_ENCAP,p_val>(g_src.get_o(it.get()),g_dst.get_o(it.get()));

                ++it;
            }

            Point<dim,St> offset = st_g.getOffset(i);
            Point<dim,St> spacing = st_g.getSpacing();
            Box<dim,size_t> dom = st_g.getDomain(i);

            vtk_w.add(g_dst,offset,spacing,dom);
        }

        vtk_w.write("vtk_grids_st_" + std::to_string(p_id) + "_" + std::to_string(p_val) + ".vtk");
    }

    template<unsigned int p_val> void out_staggered()
    {
        // property type
        typedef typename boost::mpl::at< typename st_grid::value_type::type , typename boost::mpl::int_<p_val> >::type ele;

        // Eliminate the extends
        typedef typename std::remove_all_extents<ele>::type r_ele;

        // create an openfpm format object from the property type
        typedef object<typename boost::fusion::vector<r_ele>> d_object;

        VTKWriter<boost::mpl::pair<grid_cpu<dim, d_object >,St>,VECTOR_ST_GRIDS> vtk_w;

        // Create a vector of grids
        openfpm::vector< grid_cpu<dim, d_object > > vg;
        vg.reserve(st_g.getN_loc_grid() * extends<ele>::mul());

        size_t k = 0;

        // for each domain grid
        for (size_t i = 0 ; i < st_g.getN_loc_grid() ; i++)
        {
            write_stag<ele>::template write<p_val, st_grid,openfpm::vector< grid_cpu<dim, d_object > > >(st_g,vg,i);

            // for each component
            for ( ; k < vg.size() ; k++)
            {
                Point<dim,St> offset = st_g.getOffset(i);
                Point<dim,St> spacing = st_g.getSpacing();
                Box<dim,size_t> dom = st_g.getDomain(i);

                vtk_w.add(i,vg.get(k),offset,spacing,dom,st_g.c_prp[p_val].get(k));
            }

            k = vg.size();
        }

        vtk_write<typename st_grid::value_type,VTKWriter<boost::mpl::pair<grid_cpu<dim, d_object >,St>,VECTOR_ST_GRIDS>> v(vtk_w,"vtk_grids_st_" + std::to_string(p_id),p_val);
    }

    template<typename T>
    void operator()(T& t)
    {
        if (st_g.is_staggered_prop(T::value) == false)
            out_normal<T::value>();
        else
            out_staggered<T::value>();
    }
};

template<typename Eqs_sys, typename it1_type, typename it2_type> bool checkIterator(const it1_type & it1, const it2_type & it2)
{
#ifdef SE_CLASS1

    grid_key_dx<Eqs_sys::dims> it1_k = it1.getStop() - it1.getStart();
    grid_key_dx<Eqs_sys::dims> it2_k = it2.getStop() - it2.getStart();

    for (size_t i = 0 ; i < Eqs_sys::dims ; i++)
    {
        if (it1_k.get(i) !=  it2_k.get(i))
        {
            std::cerr << __FILE__ << ":" << __LINE__ << " error src iterator and destination iterator does not match in size\n";
            return false;
        }
    }

    return true;
#else

    return true;

#endif
}

template<unsigned int dim, typename v_prp_id, typename v_prp_type>
struct interp_points
{
/*#ifdef SE_CLASS3

    // number of properties we are processing
    typedef boost::mpl::size<v_prp_id> v_size_old;

    typedef boost::mpl::int_<v_size_old::value-3> v_size;

#else*/

    // number of properties we are processing
    typedef boost::mpl::size<v_prp_id> v_size;

//#endif

    // interpolation points for each property
    openfpm::vector<std::vector<comb<dim>>> (& interp_pts)[v_size::value];

    // staggered position for each property
    const openfpm::vector<comb<dim>> (&stag_pos)[v_size::value];

    inline interp_points(openfpm::vector<std::vector<comb<dim>>> (& interp_pts)[v_size::value],const openfpm::vector<comb<dim>> (&stag_pos)[v_size::value])
    :interp_pts(interp_pts),stag_pos(stag_pos){};

    template<typename T>
    inline void operator()(T& t)
    {
        // This is the type of the object we have to copy
        typedef typename boost::mpl::at_c<v_prp_type,T::value>::type prp_type;

        interp_pts[T::value].resize(stag_pos[T::value].size());

        for (size_t i = 0 ; i < stag_pos[T::value].size() ; i++)
        {
            // Create the interpolation points
            interp_pts[T::value].get(i) = SubHyperCube<dim,dim - std::rank<prp_type>::value>::getCombinations_R(stag_pos[T::value].get(i),0);

            // interp_point are -1,0,1, map the -1 to 0 and 1 to -1
            for (size_t j = 0 ; j < interp_pts[T::value].get(i).size() ; j++)
            {
                for (size_t k = 0 ; k < dim ; k++)
                    interp_pts[T::value].get(i)[j].getComb()[k] = - ((interp_pts[T::value].get(i)[j].getComb()[k] == -1)?0:interp_pts[T::value].get(i)[j].getComb()[k]);
            }
        }
    }
};

#endif /* SRC_GRID_STAGGERED_DIST_GRID_UTIL_HPP_ */