CartesianGraphFactory.hpp

/*
 * CartesianGraphFactory.hpp
 *
 *  Created on: Nov 28, 2014
 *      Author: i-bird
 */

#ifndef CARTESIANGRAPHFACTORY_HPP_
#define CARTESIANGRAPHFACTORY_HPP_

#include "Vector/map_vector.hpp"
#include "Graph/map_graph.hpp"
#include "Grid/grid_sm.hpp"
#include "Space/Shape/Box.hpp"
#include "Space/Shape/HyperCube.hpp"

#define NO_VERTEX_ID -1

template<unsigned int dim, typename G_v, int prp>
struct fill_id
{
    static inline void fill(G_v & g_v, const grid_key_dx<dim> & gk, const grid_sm<dim, void> & gs)
    {
        g_v.template get<prp>() = gs.LinId(gk);
    }
};
template<unsigned int dim, typename G_v>
struct fill_id<dim, G_v, NO_VERTEX_ID>
{
    static inline void fill(G_v & g_v, const grid_key_dx<dim> & gk, const grid_sm<dim, void> & gs)
    {
    }
};

template<unsigned int dim, int lin_id, typename dT, typename G_v, typename v, int impl>
class fill_prop
{
    const Box<dim,dT> & domain;

    const dT (&szd)[dim];

    grid_key_dx<dim> & gk;

    G_v & g_v;

    const grid_sm<dim, void> & gs;

public:

    fill_prop(G_v & g_v, const dT (&szd)[dim], grid_key_dx<dim> & gk, const grid_sm<dim, void> & gs, const Box<dim,dT> & domain)
    :domain(domain), szd(szd), gk(gk), g_v(g_v), gs(gs)
    {
    }

    template<typename T>
    void operator()(T& t) const
    {
        typedef typename boost::fusion::result_of::at<v, boost::mpl::int_<T::value>>::type t_val;

        g_v.template get<t_val::value>() = gk.get(T::value) * szd[T::value] + domain.getLow(T::value);
        fill_id<dim, G_v, lin_id>::fill(g_v, gk, gs);
    }
};

template<unsigned int dim, int lin_id, typename dT, typename G_v, typename v>
class fill_prop<dim, lin_id, dT, G_v, v, 0>
{

public:

    fill_prop(G_v & g_v, const dT (&szd)[dim], grid_key_dx<dim> & gk, const grid_sm<dim, void> & gs, const Box<dim,dT> & domain)
    {
    }

    template<typename T>
    void operator()(T& t) const
    {
    }
};

template<unsigned int dim, int lin_id, typename dT, typename G_v, typename v>
class fill_prop<dim, lin_id, dT, G_v, v, 2>
{
    const Box<dim,dT> & domain;

    const dT (&szd)[dim];

    grid_key_dx<dim> & gk;

    G_v & g_v;

    const grid_sm<dim, void> & gs;

public:

    fill_prop(G_v & g_v, const dT (&szd)[dim], grid_key_dx<dim> & gk, const grid_sm<dim, void> & gs, const Box<dim,dT> & domain)
    :domain(domain), szd(szd), gk(gk), g_v(g_v), gs(gs)
    {
    }

    template<typename T>
    void operator()(T& t) const
    {
        typedef typename boost::fusion::result_of::at<v, boost::mpl::int_<0>>::type t_val;
        typedef typename boost::mpl::at<typename G_v::T_type::type,t_val>::type s_type;

        for (size_t i = 0 ; i < std::extent<s_type>::value ; i++)
            g_v.template get<t_val::value>()[i] = 0.0;

        for (size_t i = 0 ; i < dim ; i++)
            g_v.template get<t_val::value>()[i] = gk.get(i) * static_cast<float>(szd[i]) + domain.getLow(i);

        fill_id<dim, G_v, lin_id>::fill(g_v, gk, gs);
    }
};

template<int i, typename p, typename Graph, int ... pos>
struct fill_prop_by_type
{

    typedef typename boost::mpl::at<p, boost::mpl::int_<0>>::type v_element;

    typedef typename boost::mpl::at<typename Graph::V_type::type, v_element>::type pos_prop_type;

    enum
    {
        value = ((sizeof...(pos) != 0) * (std::is_array<pos_prop_type>::value + 1))
    };

};

template<typename p, typename Graph, int ... pos>
struct fill_prop_by_type<0, p, Graph, pos...>
{
    enum
    {
        value = 0
    };

};

template<unsigned int dim, int lin_id, typename Graph, int se, typename T, unsigned int dim_c, int ... pos>
class Graph_constructor_impl
{
public:

    static Graph construct(const size_t (& sz)[dim], Box<dim,T> dom, const size_t(& bc)[dim])
    {
        // Calculate the size of the hyper-cubes on each dimension
        T szd[dim];

        for (size_t i = 0; i < dim; i++)
        {
            szd[i] = (dom.getHigh(i) - dom.getLow(i)) / sz[i];
        }


        HyperCube<dim> hc;

        // Construct a grid info

        grid_sm<dim, void> g(sz);

        // Create a graph with the number of vertices equal to the number of
        // grid point


        Graph gp(g.size());

        /******************
         *
         * Create the edges and fill spatial
         * information properties
         *
         ******************/

        grid_key_dx_iterator<dim> k_it(g);


        while (k_it.isNext())
        {
            grid_key_dx<dim> key = k_it.get();

            // Vertex object

            auto obj = gp.vertex(g.LinId(key));

            typedef typename to_boost_vmpl<pos...>::type p;

            // vertex spatial properties functor

            fill_prop<dim, lin_id, T, decltype(gp.vertex(g.LinId(key))), typename to_boost_vmpl<pos...>::type, fill_prop_by_type<sizeof...(pos), p, Graph, pos...>::value> flp(obj, szd, key, g, dom);

            // fill properties

            boost::mpl::for_each<boost::mpl::range_c<int, 0, sizeof...(pos)> >(flp);

            // Get the combinations of dimension d

            for (long int d = dim-1 ; d >= dim_c ; d--)
            {
                // create the edges for that dimension

                std::vector<comb<dim>> c = hc.getCombinations_R(d);

                // for each combination calculate a safe linearization and create an edge

                for (size_t j = 0; j < c.size(); j++)
                {
                    // Calculate the element size

                    T ele_sz = 0;

                    // for each dimension multiply and reduce


                    for (size_t s = 0 ; s < dim ; s++)
                        ele_sz += szd[s] * abs(c[j][s]);

                    // Calculate the end point vertex id
                    // Calculate the start point id

                    size_t start_v = g.LinId(key);

                    size_t end_v = g.template LinId<CheckExistence>(key,c[j].getComb(),bc);

                    // Add an edge and set the the edge property to the size of the face (communication weight)
                    gp.template addEdge<CheckExistence>(start_v, end_v).template get<se>() = ele_sz;
                }
            }

            // Fill vertex properties

            ++k_it;
        }

        return gp;
    }
};

template<unsigned int dim, int lin_id, typename Graph, typename T, unsigned int dim_c, int ... pos>
class Graph_constructor_impl<dim, lin_id, Graph, NO_EDGE, T, dim_c, pos...>
{
public:

    static Graph construct(const size_t ( & sz)[dim], Box<dim,T> dom, const size_t(& bc)[dim])
    {
        // Calculate the size of the hyper-cubes on each dimension

        T szd[dim];

        for (size_t i = 0; i < dim; i++)
        {
            szd[i] = (dom.getHigh(i) - dom.getLow(i)) / sz[i];
        }


        HyperCube<dim> hc;

        // Construct a grid info

        grid_sm<dim, void> g(sz);

        // Create a graph with the number of vertices equal to the number of
        // grid point


        Graph gp(g.size());

        /******************
         *
         * Create the edges and fill spatial
         * information properties
         *
         ******************/

        grid_key_dx_iterator<dim> k_it(g);


        while (k_it.isNext())
        {
            grid_key_dx<dim> key = k_it.get();

            // Vertex object
            auto obj = gp.vertex(g.LinId(key));

            typedef typename to_boost_vmpl<pos...>::type p;

            // vertex spatial properties functor

            fill_prop<dim, lin_id, T, decltype(gp.vertex(g.LinId(key))), typename to_boost_vmpl<pos...>::type, fill_prop_by_type<sizeof...(pos), p, Graph, pos...>::value> flp(obj, szd, key, g, dom);

            // fill properties

            boost::mpl::for_each_ref<boost::mpl::range_c<int, 0, sizeof...(pos)> >(flp);

            // Get the combinations of dimension d

            for (long int d = dim-1 ; d >= dim_c ; d--)
            {
                // create the edges for that dimension

                std::vector<comb<dim>> c = hc.getCombinations_R(d);

                // for each combination calculate a safe linearization and create an edge

                for (size_t j = 0; j < c.size(); j++)
                {
                    // Calculate the end point vertex id
                    // Calculate the start point id

                    size_t start_v = g.LinId(key);

                    size_t end_v = g.template LinId<CheckExistence>(key,c[j].getComb(),bc);

                    // Add an edge and set the the edge property to the size of the face (communication weight)
                    gp.template addEdge<CheckExistence>(start_v, end_v);
                }
            }

            // Fill vertex properties

            ++k_it;
        }

        return gp;
    }
};

template<unsigned int dim, typename Graph>
class CartesianGraphFactory
{

public:

    template<int se, int id_prp, typename T, unsigned int dim_c, int ... pos>
    static Graph construct(const size_t (&sz)[dim], Box<dim, T> dom, const size_t (& bc)[dim])
    {
        return Graph_constructor_impl<dim, id_prp, Graph, se, T, dim_c, pos...>::construct(sz, dom, bc);
    }
};

#endif /* CARTESIANGRAPHFACTORY_HPP_ */