FD_expressions.hpp

/*
 * FD_expressions.hpp
 *
 *  Created on: May 7, 2020
 *      Author: i-bird
 */

#ifndef FD_EXPRESSIONS_HPP_
#define FD_EXPRESSIONS_HPP_

namespace FD
{

    constexpr int NORM_EXPRESSION = 0;
    constexpr int STAG_EXPRESSION = 1;
    constexpr int GRID_COMP = 2;

    template <typename exp1,typename exp2, typename impl>
    class grid_dist_expression_op
    {
    };

    struct g_comp {};

    template<unsigned int i>
    struct grid_dist_expression_value_impl_func_scal
    {
        template<unsigned int prp, typename base_type, typename gtype>
        static void inte(gtype & g, grid_dist_key_dx<gtype::dims> & k, comb<gtype::dims> & c_where, comb<gtype::dims> & c_o1, base_type & inte_out, int & c, int comp)
        {
            if (c_where[i] != c_o1[i])
            {
                int sign = (c_where[i] > c_o1[i])?1:-1;

                grid_dist_expression_value_impl_func_scal<i-1>::template inte<prp,base_type>(g,k,c_where,c_o1,inte_out,c,comp);
                long int x0 = k.getKeyRef().get(i);

                k.getKeyRef().set_d(i, x0 + sign);
                grid_dist_expression_value_impl_func_scal<i-1>::template inte<prp,base_type>(g,k,c_where,c_o1,inte_out,c,comp);
                k.getKeyRef().set_d(i, x0);
            }
            else
            {
                grid_dist_expression_value_impl_func_scal<i-1>::template inte<prp,base_type>(g,k,c_where,c_o1,inte_out,c,comp);
            }
        }
    };

    template<>
    struct grid_dist_expression_value_impl_func_scal<0>
    {
        template<unsigned int prp, typename base_type, typename gtype>
        static void inte(gtype & g, grid_dist_key_dx<gtype::dims> & k, comb<gtype::dims> & c_where, comb<gtype::dims> & c_o1, base_type & inte_out , int & c , int comp)
        {
            if (c_where[0] != c_o1[0])
            {
                int sign = (c_where[0] > c_o1[0])?1:-1;

                inte_out += g.template getProp<prp>(k);

                long int x0 = k.getKeyRef().get(0);

                k.getKeyRef().set_d(0, x0 + sign);
                inte_out += g.template getProp<prp>(k);
                k.getKeyRef().set_d(0, x0);
                c += 2;
            }
            else
            {
                inte_out += g.template getProp<prp>(k);
                c += 1;
            }
        }
    };

    template<typename base_type>
    struct grid_dist_expression_value_impl
    {
        typedef base_type type;

        template<unsigned int prp, typename gtype>
        static base_type inte(gtype & g, grid_dist_key_dx<gtype::dims> & k, comb<gtype::dims> & c_where, comb<gtype::dims> & c_o1, int comp)
        {
            int c = 0;
            base_type inte = 0;

            grid_dist_expression_value_impl_func_scal<gtype::dims-1>::template inte<prp,base_type>(g,k,c_where,c_o1,inte,c,comp);

            inte /= c;

            return inte;
        }

        template<unsigned int prp, typename gtype>
        static base_type value_n(gtype & g, const grid_dist_key_dx<gtype::dims> & k, int comp)
        {
            return g.template getProp<prp>(k);
        }

        template<unsigned int prp, typename gtype>
        static auto value_ref(gtype & g, const grid_dist_key_dx<gtype::dims> & k, int comp) -> decltype(g.template getProp<prp>(k))
        {
            return g.template getProp<prp>(k);
        }
    };


    template<unsigned int i>
    struct grid_dist_expression_value_impl_func_vec
    {
        template<unsigned int prp, typename base_type, typename gtype>
        static void inte(gtype & g, grid_dist_key_dx<gtype::dims> & k, comb<gtype::dims> & c_where, comb<gtype::dims> & c_o1, base_type & inte_out, int & c, int comp)
        {
            if (c_where[i] != c_o1[i])
            {
                grid_dist_expression_value_impl_func_vec<i-1>::template inte<prp,base_type>(g,k,c_where,c_o1,inte_out,c,comp);
                long int x0 = k.getKeyRef().get(i);

                int sign = (c_where[i] > c_o1[i])?1:-1;

                k.getKeyRef().set_d(i, x0 + sign);
                grid_dist_expression_value_impl_func_vec<i-1>::template inte<prp,base_type>(g,k,c_where,c_o1,inte_out,c,comp);
                k.getKeyRef().set_d(i, x0);
            }
            else
            {
                grid_dist_expression_value_impl_func_vec<i-1>::template inte<prp,base_type>(g,k,c_where,c_o1,inte_out,c,comp);
            }
        }
    };

    template<>
    struct grid_dist_expression_value_impl_func_vec<0>
    {
        template<unsigned int prp, typename base_type, typename gtype>
        static void inte(gtype & g, grid_dist_key_dx<gtype::dims> & k, comb<gtype::dims> & c_where, comb<gtype::dims> & c_o1, base_type & inte_out , int & c , int comp)
        {
            if (c_where[0] != c_o1[0])
            {
                int sign = (c_where[0] > c_o1[0])?1:-1;

                inte_out += g.template getProp<prp>(k)[comp];

                long int x0 = k.getKeyRef().get(0);

                k.getKeyRef().set_d(0, x0 + sign);
                inte_out += g.template getProp<prp>(k)[comp];
                k.getKeyRef().set_d(0, x0);
                c += 2;
            }
            else
            {
                inte_out += g.template getProp<prp>(k)[comp];
                c += 1;
            }
        }
    };

    template<typename base_type, unsigned int N1>
    struct grid_dist_expression_value_impl<base_type[N1]>
    {
        typedef base_type type;

        template<unsigned int prp, typename gtype>
        static base_type inte(gtype & g, grid_dist_key_dx<gtype::dims> & k, comb<gtype::dims> & c_where, comb<gtype::dims> & c_o1, int comp)
        {
            int c = 0;
            base_type inte = 0;

            grid_dist_expression_value_impl_func_vec<gtype::dims-1>::template inte<prp,base_type>(g,k,c_where,c_o1,inte,c,comp);

            if (c != 0)
            {inte /= c;}
            else
            {inte = g.template getProp<prp>(k)[comp];}

            return inte;
        }

        template<unsigned int prp, typename gtype>
        static base_type value_n(gtype & g, const grid_dist_key_dx<gtype::dims> & k, int comp)
        {
            return g.template getProp<prp>(k)[comp];
        }

        template<unsigned int prp, typename gtype>
        static auto value_ref(gtype & g, const grid_dist_key_dx<gtype::dims> & k, int comp) -> decltype(g.template getProp<prp>(k)[comp])
        {
            return g.template getProp<prp>(k)[comp];
        }
    };

    template<typename base_type, unsigned int N1>
    struct grid_dist_expression_value_impl<Point<N1,base_type>>
    {
        typedef base_type type;

        template<unsigned int prp, typename gtype>
        static base_type inte(gtype & g, grid_dist_key_dx<gtype::dims> & k, comb<gtype::dims> & c_where, comb<gtype::dims> & c_o1, int comp)
        {
            int c = 0;
            base_type inte = 0;

            grid_dist_expression_value_impl_func_vec<gtype::dims-1>::template inte<prp,base_type>(g,k,c_where,c_o1,inte,c,comp);

            if (c != 0)
            {inte /= c;}
            else
            {inte = g.template getProp<prp>(k)[comp];}

            return inte;
        }

        template<unsigned int prp, typename gtype>
        static base_type value_n(gtype & g, const grid_dist_key_dx<gtype::dims> & k, int comp)
        {
            return g.template getProp<prp>(k)[comp];
        }

        template<unsigned int prp, typename gtype>
        static auto value_ref(gtype & g, const grid_dist_key_dx<gtype::dims> & k, int comp) -> decltype(g.template getProp<prp>(k)[comp])
        {
            return g.template getProp<prp>(k)[comp];
        }
    };

    template<unsigned int i>
    struct grid_dist_expression_value_impl_vnz
    {
        template<typename Sys_eqs, typename gmap_type, typename unordered_map_type ,typename gtype>
        static void value_nz(const gmap_type & g_map,
                             unordered_map_type & cols,
                             gtype & g, 
                             grid_dist_key_dx<gtype::dims> & key, 
                             comb<gtype::dims> & c_where, 
                             comb<gtype::dims> & c_o1, 
                             typename Sys_eqs::stype coeff, 
                             int & c, 
                             int comp,
                             int var_id)
        {
            if (c_where[i] != c_o1[i])
            {
                int sign = (c_where[i] > c_o1[i])?1:-1;

                grid_dist_expression_value_impl_vnz<i-1>::template value_nz<Sys_eqs>(g_map,cols,g,key,c_where,c_o1,coeff,c,comp,var_id);
                //cols[g_map.template getProp<0>(key)*Sys_eqs::nvar + var_id + comp] += coeff / c;

                long int x0 = key.getKeyRef().get(i);

                key.getKeyRef().set_d(i, x0 + sign);
                grid_dist_expression_value_impl_vnz<i-1>::template value_nz<Sys_eqs>(g_map,cols,g,key,c_where,c_o1,coeff,c,comp,var_id);
                //cols[g_map.template getProp<0>(key)*Sys_eqs::nvar + var_id + comp] += coeff / c;
                key.getKeyRef().set_d(i, x0);
            }
            else
            {grid_dist_expression_value_impl_vnz<i-1>::template value_nz<Sys_eqs>(g_map,cols,g,key,c_where,c_o1,coeff,c,comp,var_id);}
        }
    };

    template<>
    struct grid_dist_expression_value_impl_vnz<0>
    {
        template<typename Sys_eqs, typename gmap_type, typename unordered_map_type ,typename gtype>
        static void value_nz(const gmap_type & g_map,
                             unordered_map_type & cols,
                             gtype & g, 
                             grid_dist_key_dx<gtype::dims> & key, 
                             comb<gtype::dims> & c_where, 
                             comb<gtype::dims> & c_o1, 
                             typename Sys_eqs::stype coeff, 
                             int & c , 
                             int comp,
                             int var_id)
        {
            if (c_where[0] != c_o1[0])
            {
                int sign = (c_where[0] > c_o1[0])?1:-1;

                cols[g_map.template getProp<0>(key)*Sys_eqs::nvar + var_id + comp] += coeff / c;

                long int x0 = key.getKeyRef().get(0);

                key.getKeyRef().set_d(0, x0 + sign);
                cols[g_map.template getProp<0>(key)*Sys_eqs::nvar + var_id + comp] += coeff / c;
                key.getKeyRef().set_d(0, x0);
            }
            else
            {
                cols[g_map.template getProp<0>(key)*Sys_eqs::nvar + var_id + comp] += coeff;
            }
        }
    };

    template<unsigned int prp, typename grid, unsigned int impl>
    class grid_dist_expression
    {};

    template<unsigned int prp, typename grid>
    class grid_dist_expression<prp,grid,NORM_EXPRESSION>
    {
        grid & g;

        typedef typename boost::mpl::at<typename grid::value_type::type,boost::mpl::int_<prp>>::type type_proc;

    public:

        typedef grid gtype;

        static const unsigned int prop = prp;

        int var_id = 0;

        void setVarId(int var_id)
        {
            this->var_id = var_id;
        }

        grid_dist_expression(grid & g)
        :g(g)
        {}

        grid & getGrid()
        {
            return g;
        }

        const grid & getGrid() const
        {
            return g;
        }

        inline void init() const
        {}

        inline auto value(const grid_dist_key_dx<grid::dims> & k, comb<grid::dims> & c_where, int comp = 0) const -> decltype(grid_dist_expression_value_impl<type_proc>::template value_n<prp>(g,k,comp))
        {
            return grid_dist_expression_value_impl<type_proc>::template value_n<prp>(g,k,comp);
//          return g.template getProp<prp>(k);
        }

        inline auto value_ref(const grid_dist_key_dx<grid::dims> & k, comb<grid::dims> & c_where, int comp = 0) const -> decltype(grid_dist_expression_value_impl<type_proc>::template value_ref<prp>(g,k,comp))
        {
            return grid_dist_expression_value_impl<type_proc>::template value_ref<prp>(g,k,comp);
//          return g.template getProp<prp>(k);
        }

        template<unsigned int prp2> grid & operator=(const grid_dist_expression<prp2,grid,NORM_EXPRESSION> & g_exp)
        {
            g_exp.init();

            comb<grid::dims> s_pos;
            s_pos.zero;

            auto it = g.getDomainIterator();

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

                g.template getProp<prp>(key) = g_exp.value(key,s_pos);

                ++it;
            }

            return g;
        }

        template<typename exp1, typename exp2, typename op> grid & operator=(const grid_dist_expression_op<exp1,exp2,op> & g_exp)
        {
            g_exp.init();

            comb<grid::dims> s_pos;
            s_pos.mone();

            auto it = g.getDomainIterator();

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

                g.template getProp<prp>(key) = g_exp.value(key,s_pos);

                ++it;
            }

            return g;
        }

        grid & operator=(double d)
        {
            auto it = g.getDomainIterator();

            //comb<grid::dims> s_pos;
            //s_pos.mone();

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

                g.template getProp<prp>(key) = d;

                ++it;
            }

            return g;
        }


        template<typename Sys_eqs, typename gmap_type, typename unordered_map_type>
        inline void value_nz(const gmap_type & g_map,
                grid_dist_key_dx<Sys_eqs::dims> & key,
                 const grid_sm<Sys_eqs::dims,void> & gs,
                 typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
                 unordered_map_type & cols,
                 typename Sys_eqs::stype coeff,
                 unsigned int comp,
                 comb<Sys_eqs::dims> & c_where) const
        {
            cols[g_map.template getProp<0>(key)*Sys_eqs::nvar + var_id + comp] += coeff;
        }

        inline grid_dist_expression_op<grid_dist_expression<prp,grid,NORM_EXPRESSION>,boost::mpl::int_<1>,g_comp> operator[](int comp)
        {
            int comp_n[1];

            comp_n[0] = comp;

            grid_dist_expression_op<grid_dist_expression<prp,grid,NORM_EXPRESSION>,boost::mpl::int_<1>,g_comp> v_exp(*this,comp_n,var_id);

            return v_exp;
        }

        //Need more treatment for staggered (c_where based on exp)
        inline typename gtype::stype get(grid_dist_key_dx<gtype::dims> & key)
        {
            comb<gtype::dims> c_where;
            c_where.zero();
            return this->value(key,c_where);
        }

        int isConstant(){
            return false;
        }
    };

    template<unsigned int prp, typename grid>
    class grid_dist_expression<prp,grid,STAG_EXPRESSION>
    {
        grid & g;

        typedef typename boost::mpl::at<typename grid::value_type::type,boost::mpl::int_<prp>>::type type_proc;

    public:

        typedef grid gtype;

        static const unsigned int prop = prp;

        int var_id = 0;

        void setVarId(int var_id)
        {
            this->var_id = var_id;
        }

        grid_dist_expression(grid & g)
        :g(g)
        {}

        grid & getGrid()
        {
            return g;
        }

        const grid & getGrid() const
        {
            return g;
        }

        inline void init() const
        {}

        inline auto value_ref(const grid_dist_key_dx<grid::dims> & k, comb<grid::dims> & c_where,int comp = 0) const -> decltype(grid_dist_expression_value_impl<type_proc>::template value_ref<prp>(g,k,comp))
        {
            return grid_dist_expression_value_impl<type_proc>::template value_ref<prp>(g,k,comp);
            //return g.template getProp<prp>(k);
        }

        inline auto value(grid_dist_key_dx<grid::dims> & k, comb<grid::dims> & c_where, int comp = 0) const -> decltype(grid_dist_expression_value_impl<type_proc>::template inte<prp>(g,k,c_where,c_where,comp))
        {
            comb<grid::dims> c_o1 = g.getStagPositions()[prp].get(comp);

            return grid_dist_expression_value_impl<type_proc>::template inte<prp>(g,k,c_where,c_o1,comp);
        }

        template<unsigned int prp2> grid & operator=(const grid_dist_expression<prp2,grid,STAG_EXPRESSION> & g_exp)
        {
            g_exp.init();

            auto it = g.getDomainIterator();

            comb<grid::dims> s_pos = g.getStagPosition()[prp].get(0);

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

                g.template getProp<prp>(key) = g_exp.value(key,s_pos);

                ++it;
            }

            return g;
        }

        template<typename exp1, typename exp2, typename op> grid & operator=(const grid_dist_expression_op<exp1,exp2,op> & g_exp)
        {
            g_exp.init();

            auto it = g.getDomainIterator();

            comb<grid::dims> s_pos = g.getStagPositions()[prp].get(0);

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

                g.template getProp<prp>(key) = g_exp.value(key,s_pos);

                ++it;
            }

            return g;
        }

        grid & operator=(double d)
        {
            auto it = g.getDomainIterator();

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

                g.template getProp<prp>(key) = d;

                ++it;
            }

            return g;
        }

        template<typename Sys_eqs, typename gmap_type, typename unordered_map_type>
        inline void value_nz(const gmap_type & g_map,
                grid_dist_key_dx<Sys_eqs::dims> & key,
                 const grid_sm<Sys_eqs::dims,void> & gs,
                 typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
                 unordered_map_type & cols,
                 typename Sys_eqs::stype coeff,
                 unsigned int comp,
                 comb<Sys_eqs::dims> & c_where) const
        {
            int c = 1;
            comb<grid::dims> c_o1 = g.getStagPositions()[prp].get(comp);

            // x0, dx are defined in proper dir є(x, y, z)

            for (int i = 0 ; i < grid::dims ; i++)
            {
                if (c_where[i] != c_o1[i])
                {c *= 2;}
            }

            grid_dist_expression_value_impl_vnz<Sys_eqs::dims-1>::template value_nz<Sys_eqs>(g_map,cols,g,key,c_where,c_o1,coeff,c,comp,var_id);

        }

        inline grid_dist_expression_op<grid_dist_expression<prp,grid,STAG_EXPRESSION>,boost::mpl::int_<1>,g_comp> operator[](int comp)
        {
            int comp_n[1];

            comp_n[0] = comp;

            grid_dist_expression_op<grid_dist_expression<prp,grid,STAG_EXPRESSION>,boost::mpl::int_<1>,g_comp> v_exp(*this,comp_n,var_id);

            return v_exp;
        }

    };

    template<unsigned int dim>
    class grid_dist_expression<dim,double,NORM_EXPRESSION>
    {
        double d;

    public:

        inline grid_dist_expression(const double & d)
        :d(d)
        {}

        inline void init() const
        {}

        inline double value(const grid_dist_key_dx<dim> & k, comb<dim> & c_where) const
        {
            return d;
        }
    };

    template<unsigned int dim>
    class grid_dist_expression<dim,double,STAG_EXPRESSION>
    {
        double d;

    public:

        inline grid_dist_expression(const double & d)
        :d(d)
        {}

        inline void init() const
        {}

        inline double value(const grid_dist_key_dx<dim> & k, comb<dim> & c_where) const
        {
            return d;
        }
    };

    template<unsigned int dim, unsigned int impl>
    class grid_dist_expression<dim,float,impl>
    {
        float d;

    public:

        typedef float vtype;

        inline grid_dist_expression(const float & d)
        :d(d)
        {}

        inline void init() const
        {}

        inline float value(const grid_dist_key_dx<dim> & k) const
        {
            return d;
        }
    };

    struct sum {};

    template <typename exp1, typename exp2>
    class grid_dist_expression_op<exp1,exp2,sum>
    {
        const exp1 o1;

        const exp2 o2;

    public:

        typedef typename exp1::gtype gtype;

        inline grid_dist_expression_op(const exp1 & o1,const exp2 & o2)
                :o1(o1),o2(o2)
        {}

        inline void init() const
        {
            o1.init();
            o2.init();
        }

        inline auto value(grid_dist_key_dx<gtype::dims> & key, comb<gtype::dims> & c_where) const -> typename std::remove_reference<decltype(o1.value(key,c_where))>::type
        {
            typename std::remove_reference<decltype(o1.value(key,c_where))>::type val;

            return o1.value(key,c_where) + o2.value(key,c_where);
        }

        gtype & getGrid()
        {
            return o1.getGrid();
        }

        const gtype & getGrid() const
        {
            return o1.getGrid();
        }

        template<typename Sys_eqs, typename gmap_type, typename unordered_map_type>
        inline void value_nz(const gmap_type & g_map,
                             grid_dist_key_dx<Sys_eqs::dims> & key,
                             const grid_sm<Sys_eqs::dims,void> & gs,
                             typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
                             unordered_map_type & cols,
                             typename Sys_eqs::stype coeff,
                             unsigned int comp,
                             comb<Sys_eqs::dims> & c_where) const
        {
            o1.template value_nz<Sys_eqs>(g_map,key,gs,spacing,cols,coeff,comp,c_where);
            o2.template value_nz<Sys_eqs>(g_map,key,gs,spacing,cols,coeff,comp,c_where);
        }
    };

    struct sub {};


    template <typename exp1, typename exp2>
    class grid_dist_expression_op<exp1,exp2,sub>
    {
        const exp1 o1;

        const exp2 o2;

    public:

        typedef typename exp1::gtype gtype;

        inline grid_dist_expression_op(const exp1 & o1,const exp2 & o2)
                :o1(o1),o2(o2)
        {}

        inline void init() const
        {
            o1.init();
            o2.init();
        }

        inline auto value(grid_dist_key_dx<gtype::dims> & key, comb<gtype::dims> & c_where) const -> typename std::remove_reference<decltype(o1.value(key,c_where))>::type
        {
            typename std::remove_reference<decltype(o1.value(key,c_where))>::type val;

            return o1.value(key,c_where) - o2.value(key,c_where);
        }

        gtype & getGrid()
        {
            return o1.getGrid();
        }

        const gtype & getGrid() const
        {
            return o1.getGrid();
        }

        template<typename Sys_eqs, typename gmap_type, typename unordered_map_type>
        inline void value_nz(const gmap_type & g_map,
                             grid_dist_key_dx<Sys_eqs::dims> & key,
                             const grid_sm<Sys_eqs::dims,void> & gs,
                             typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
                             unordered_map_type & cols,
                             typename Sys_eqs::stype coeff,
                             unsigned int comp,
                             comb<Sys_eqs::dims> & c_where) const
        {
            o1.template value_nz<Sys_eqs>(g_map,key,gs,spacing,cols,coeff,comp,c_where);
            o2.template value_nz<Sys_eqs>(g_map,key,gs,spacing,cols,-coeff,comp,c_where);
        }
    };

    struct subuni{};

    template <typename exp1>
    class grid_dist_expression_op<exp1,void,subuni>
    {
        const exp1 o1;

    public:

        typedef typename exp1::gtype gtype;

        inline grid_dist_expression_op(const exp1 & o1)
                :o1(o1)
        {}

        inline void init() const
        {
            o1.init();
        }

        inline auto value(grid_dist_key_dx<gtype::dims> & key, comb<gtype::dims> & c_where) const -> typename std::remove_reference<decltype(o1.value(key,c_where))>::type
        {
            typename std::remove_reference<decltype(o1.value(key,c_where))>::type val;

            return -o1.value(key,c_where);
        }

        gtype & getGrid()
        {
            return o1.getGrid();
        }

        const gtype & getGrid() const
        {
            return o1.getGrid();
        }

        template<typename Sys_eqs, typename gmap_type, typename unordered_map_type>
        inline void value_nz(const gmap_type & g_map,
                             grid_dist_key_dx<Sys_eqs::dims> & key,
                             const grid_sm<Sys_eqs::dims,void> & gs,
                             typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
                             unordered_map_type & cols,
                             typename Sys_eqs::stype coeff,
                             unsigned int comp,
                             comb<Sys_eqs::dims> & c_where) const
        {
            o1.template value_nz<Sys_eqs>(g_map,key,gs,spacing,cols,-coeff,comp,c_where);
        }
    };

    struct mul {};

    template <typename exp1, typename exp2>
    class grid_dist_expression_op<exp1,exp2,mul>
    {
        const exp1 o1;

        const exp2 o2;

    public:

        typedef typename exp2::gtype gtype;

        inline grid_dist_expression_op(const exp1 & o1,const exp2 & o2)
                :o1(o1),o2(o2)
        {}

        inline void init() const
        {
            o1.init();
            o2.init();
        }

        inline auto value(grid_dist_key_dx<gtype::dims> & key, comb<gtype::dims> & c_where) const -> typename std::remove_reference<decltype(o1.value(key,c_where))>::type
        {
            typename std::remove_reference<decltype(o1.value(key,c_where))>::type val;

            return o1.value(key,c_where) * o2.value(key,c_where);
        }

        gtype & getGrid()
        {
            return o1.getGrid();
        }

        const gtype & getGrid() const
        {
            return o1.getGrid();
        }

        template<typename Sys_eqs, typename gmap_type, typename unordered_map_type>
        inline void value_nz(const gmap_type & g_map,
                grid_dist_key_dx<Sys_eqs::dims> & key,
                 const grid_sm<Sys_eqs::dims,void> & gs,
                 typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
                 unordered_map_type & cols,
                 typename Sys_eqs::stype coeff,
                 unsigned int comp,
                 comb<Sys_eqs::dims> & c_where) const
        {
            typename Sys_eqs::stype coeff_tmp = o1.value(key,c_where) * coeff;

            o2.template value_nz<Sys_eqs>(g_map,key,gs,spacing,cols,coeff_tmp,comp,c_where);
        }
    };

    template <typename grid_type, unsigned int prp>
    struct pos_or_propL
    {
        __device__ __host__ static inline auto value(grid_type & v, const grid_dist_key_dx<grid_type::dims> & k) -> decltype(v.template getProp<prp>(k))
        {
            return v.template getProp<prp>(k);
        }
    };

    template<unsigned int, bool is_valid>
    struct get_grid_dist_expression_op
    {
        template<typename exp_type>
        inline static auto get(exp_type & o1, grid_dist_key_dx<exp_type::gtype::dims> & key, comb<exp_type::gtype::dims> & c_where) -> decltype(o1.value(key,c_where))
        {
            return o1.value(key,c_where);
        }

        template<typename exp_type>
        inline static auto get_ref(exp_type & o1, grid_dist_key_dx<exp_type::gtype::dims> & key, comb<exp_type::gtype::dims> & c_where) -> decltype(o1.value_ref(key,c_where))
        {
            return o1.value_ref(key,c_where);
        }

        template<unsigned int prop, typename exp_type, typename grid_type>
        inline static void assign(exp_type & o1, grid_type & g, const grid_dist_key_dx<exp_type::gtype::dims> & key)
        {
            pos_or_propL<grid_type,exp_type::prop>::value(g,key) = o1.value(key);
        }

        template<unsigned int prop, typename grid_type>
        inline static void assign_double(double d, grid_type & g, const grid_dist_key_dx<grid_type::dims> & key)
        {
            pos_or_propL<grid_type,prop>::value(g,key) = d;
        }
    };

    template<>
    struct get_grid_dist_expression_op<1,false>
    {
        template<typename exp_type>
        static int get(exp_type & o1, grid_dist_key_dx<exp_type::gtype::dims> & key, comb<exp_type::gtype::dims> & c_where, const int (& comp)[1])
        {
            return 0;
        }

        template<unsigned int prop, typename exp_type, typename grid_type>
        inline static void assign(exp_type & o1, grid_type & g, const grid_dist_key_dx<exp_type::gtype::dims> & key)
        {
        }

        template<unsigned int prop, typename grid_type>
        inline static void assign_double(double d, grid_type & g, const grid_dist_key_dx<grid_type::dims> & key)
        {
        }
    };

    template<>
    struct get_grid_dist_expression_op<1,true>
    {
        template<typename exp_type>
        static auto get(exp_type & o1, grid_dist_key_dx<exp_type::gtype::dims> & key, comb<exp_type::gtype::dims> & c_where, const int (& comp)[1]) -> decltype(o1.value(key,c_where,comp[0]))
        {
            return o1.value(key,c_where,comp[0]);
        }

        template<typename exp_type>
        static auto get_ref(exp_type & o1, grid_dist_key_dx<exp_type::gtype::dims> & key, comb<exp_type::gtype::dims> & c_where, const int (& comp)[1]) -> decltype(o1.value_ref(key,c_where,comp[0]))
        {
            return o1.value_ref(key,c_where,comp[0]);
        }

        template<unsigned int prop,typename exp_type, typename grid_type>
        inline static void assign(exp_type & o1, grid_type & g, grid_dist_key_dx<exp_type::gtype::dims> & key,comb<exp_type::gtype::dims> & c_where, const int (& comp)[1])
        {
            pos_or_propL<grid_type,prop>::value(g,key)[comp[0]] = o1.value(key,c_where);
        }

        template<unsigned int prop, typename grid_type>
        inline static void assign_double(double d, grid_type & g, const grid_dist_key_dx<grid_type::dims> & key, const int (& comp)[1])
        {
            pos_or_propL<grid_type,prop>::value(g,key)[comp[0]] = d;
        }
    };

    template<>
    struct get_grid_dist_expression_op<2,true>
    {
        template<typename exp_type>
        static auto get(exp_type & o1, grid_dist_key_dx<exp_type::gtype::dims> & key, comb<exp_type::gtype::dims> & c_where, const int (& comp)[2]) -> decltype(o1.value(key,c_where)[0][0])
        {
            return o1.value(key)[comp[0]][comp[1]];
        }

        template<unsigned int prop,typename exp_type, typename grid_type>
        inline static void assign(exp_type & o1, grid_type & g, const grid_dist_key_dx<grid_type::dims> & key, const int (& comp)[2])
        {
            pos_or_propL<grid_type,prop>::value(g,key)[comp[0]][comp[1]] = o1.value(key);
        }

        template<unsigned int prop, typename grid_type>
        inline static void assign_double(double d, grid_type & g, const grid_dist_key_dx<grid_type::dims> & key, const int (& comp)[2])
        {
            pos_or_propL<grid_type,prop>::value(g,key)[comp[0]][comp[1]] = d;
        }
    };

    template <typename exp1,int n>
    class grid_dist_expression_op<exp1,boost::mpl::int_<n>,g_comp>
    {
        exp1 o1;

        int comp[n];

        int var_id = 0;
        void setVarId(int var_id)
        {
            this->var_id = var_id;
        }

        typedef grid_dist_expression_op<exp1,boost::mpl::int_<n>,g_comp> myself;

    public:

            typedef std::false_type is_ker;

        typedef typename exp1::gtype gtype;

        static const unsigned int prop = exp1::prop;


        grid_dist_expression_op(const exp1 & o1, int (& comp)[n], int var_id)
        :o1(o1),var_id(var_id)
        {
            for (int i = 0 ; i < n ; i++)
            {this->comp[i] = comp[i];}
        }

        const gtype & getGrid() const
        {
            return o1.getGrid();
        }

        gtype & getGrid()
        {
            return o1.getGrid();
        }

        inline void init() const
        {
            o1.init();
        }

        typedef typename boost::mpl::at<typename gtype::value_type::type,boost::mpl::int_<exp1::prop>>::type property_act;

        inline auto value(grid_dist_key_dx<gtype::dims> & key, comb<gtype::dims> & c_where) const -> decltype(get_grid_dist_expression_op<n,n == rank_gen<property_act>::type::value>::get(o1,key,c_where,comp))
        {
            return get_grid_dist_expression_op<n,n == rank_gen<property_act>::type::value>::get(o1,key,c_where,comp);
        }

        inline auto value_ref(grid_dist_key_dx<gtype::dims> & key, comb<gtype::dims> & c_where) const -> decltype(get_grid_dist_expression_op<n,n == rank_gen<property_act>::type::value>::get_ref(o1,key,c_where,comp))
        {
            return get_grid_dist_expression_op<n,n == rank_gen<property_act>::type::value>::get_ref(o1,key,c_where,comp);
        }

        inline auto get(grid_dist_key_dx<gtype::dims> & key, comb<gtype::dims> & c_where) const -> decltype(value(key,c_where))
        {
            return this->value(key,c_where);
        }

        template<typename Sys_eqs, typename gmap_type, typename unordered_map_type>
        inline void value_nz(const gmap_type & g_map,
                grid_dist_key_dx<Sys_eqs::dims> & key,
                 const grid_sm<Sys_eqs::dims,void> & gs,
                 typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
                 unordered_map_type & cols,
                 typename Sys_eqs::stype coeff,
                 unsigned int comp_,
                 comb<Sys_eqs::dims> & c_where) const
        {
    #ifdef SE_CLASS1

            if (n != 1)
            {
                std::cout << __FILE__ << ":" << __LINE__ << " Error it only work for tensore of rank 1 ... like vectors " << std::endl;
            }

    #endif

            o1.template value_nz<Sys_eqs>(g_map,key,gs,spacing,cols,coeff,comp_ + var_id + comp[0],c_where);
        }

        inline grid_dist_expression_op<exp1,boost::mpl::int_<2>,g_comp> operator[](int comp_)
        {
            int comp_n[n+1];

            for (int i = 0 ; i < n ; i++)
            {comp_n[i] = comp[i];}
            comp_n[n] = comp_;

            grid_dist_expression_op<exp1,boost::mpl::int_<2>,g_comp> v_exp(o1,comp_n,var_id);

            return v_exp;
        }

        //Need more treatment for staggered (c_where based on exp)
        inline typename gtype::stype get(grid_dist_key_dx<gtype::dims> & key)
        {
            comb<gtype::dims> c_where;
            c_where.zero();
            return this->value(key,c_where);
        }


        template<unsigned int prp2, unsigned int impl> gtype & operator=(const grid_dist_expression<prp2,gtype,impl> & v_exp)
        {
            v_exp.init();

            auto & g = getGrid();

            auto it = g.getDomainIterator();

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

                get_grid_dist_expression_op<n,n == rank_gen<property_act>::type::value>::template assign<exp1::prop>(v_exp,g,key,comp);

                ++it;
            }

            return g;
        }

        template<typename exp1_, typename exp2_, typename op> gtype & operator=(const grid_dist_expression_op<exp1_,exp2_,op> & v_exp)
        {
            v_exp.init();

            auto & g = getGrid();

            auto it = g.getDomainIterator();

            comb<gtype::dims> c_where;
            c_where.zero();

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

                get_grid_dist_expression_op<n,n == rank_gen<property_act>::type::value>::template assign<exp1::prop>(v_exp,g,key,c_where,comp);

                ++it;
            }

            return g;
        }

        gtype & operator=(double d)
        {
            auto & v = getGrid();

            auto it = v.getDomainIterator();

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

                //pos_or_propL<vtype,exp1::prp>::value(v,key) = d;
                get_grid_dist_expression_op<n,n == rank_gen<property_act>::type::value>::template assign_double<exp1::prop>(d,v,key,comp);


                ++it;
            }

            return v;
        }

        int isConstant(){
            return false;
        }
    };


    template <unsigned int prp,typename grid> inline grid_dist_expression<prp,grid,NORM_EXPRESSION> getV(grid & g)
    {
        grid_dist_expression<prp,grid,NORM_EXPRESSION> exp_g(g);

        return exp_g;
    }

    template <unsigned int prp, typename grid> inline grid_dist_expression<prp,grid,STAG_EXPRESSION> getV_stag(grid & g)
    {
        grid_dist_expression<prp,grid,STAG_EXPRESSION> exp_g(g);

        return exp_g;
    }

};

/* \brief sum two distributed grid expression
 *
 * \param ga grid expression one
 * \param gb grid expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int p1, unsigned int p2, typename g1, typename g2, unsigned int impl_p1, unsigned int impl_p2>
inline FD::grid_dist_expression_op<FD::grid_dist_expression<p1,g1,impl_p1>,FD::grid_dist_expression<p2,g2,impl_p2>,FD::sum>
operator+(const FD::grid_dist_expression<p1,g1,impl_p1> & ga, const FD::grid_dist_expression<p2,g2,impl_p2> & gb)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression<p1,g1,impl_p1>,FD::grid_dist_expression<p2,g2,impl_p2>,FD::sum> exp_sum(ga,gb);

    return exp_sum;
}

/* \brief sum two distributed grid expression
 *
 * \param ga grid expression one
 * \param gb grid expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1 , typename exp2, typename op1, unsigned int prp1, typename g1, unsigned int impl_p1>
inline FD::grid_dist_expression_op<FD::grid_dist_expression_op<exp1,exp2,op1>,FD::grid_dist_expression<prp1,g1,impl_p1>,FD::sum>
operator+(const FD::grid_dist_expression_op<exp1,exp2,op1> & ga, const FD::grid_dist_expression<prp1,g1,impl_p1> & gb)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression_op<exp1,exp2,op1>,FD::grid_dist_expression<prp1,g1,impl_p1>,FD::sum> exp_sum(ga,gb);

    return exp_sum;
}

/* \brief sum two distributed grid expression
 *
 * \param ga grid expression one
 * \param gb grid expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1 , typename exp2, typename op1, unsigned int prp1, typename g1, unsigned int impl_p1>
inline FD::grid_dist_expression_op<FD::grid_dist_expression<prp1,g1,impl_p1>,FD::grid_dist_expression_op<exp1,exp2,op1>,FD::sum>
operator+(const FD::grid_dist_expression<prp1,g1,impl_p1> & ga, const FD::grid_dist_expression_op<exp1,exp2,op1> & gb)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression<prp1,g1,impl_p1>,FD::grid_dist_expression_op<exp1,exp2,op1>,FD::sum> exp_sum(ga,gb);

    return exp_sum;
}

/* \brief sum two distributed grid expression
 *
 * \param ga grid expression one
 * \param gb grid expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1 , typename exp2, typename op1, typename exp3 , typename exp4, typename op2>
inline FD::grid_dist_expression_op<FD::grid_dist_expression_op<exp1,exp2,op1>,FD::grid_dist_expression_op<exp3,exp4,op2>,FD::sum>
operator+(const FD::grid_dist_expression_op<exp1,exp2,op1> & ga, const FD::grid_dist_expression_op<exp3,exp4,op2> & gb)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression_op<exp1,exp2,op1>,FD::grid_dist_expression_op<exp3,exp4,op2>,FD::sum> exp_sum(ga,gb);

    return exp_sum;
}

/* \brief sum two distributed grid expression
 *
 * \param ga grid expression one
 * \param gb grid expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int prp1 , typename g1, unsigned int impl_p1>
inline FD::grid_dist_expression_op<FD::grid_dist_expression<prp1,g1,impl_p1>,FD::grid_dist_expression<g1::dims,double,impl_p1>,FD::sum>
operator+(const FD::grid_dist_expression<prp1,g1,impl_p1> & ga, double d)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression<prp1,g1,impl_p1>,FD::grid_dist_expression<g1::dims,double,impl_p1>,FD::sum> exp_sum(ga,FD::grid_dist_expression<g1::dims,double,impl_p1>(d));

    return exp_sum;
}

/* \brief sum two distributed grid expression
 *
 * \param ga grid expression one
 * \param gb grid expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int prp1 , typename g1, unsigned int impl_p1>
inline FD::grid_dist_expression_op<FD::grid_dist_expression<g1::dims,double,impl_p1>,FD::grid_dist_expression<prp1,g1,impl_p1>,FD::sum>
operator+(double d, const FD::grid_dist_expression<prp1,g1,impl_p1> & gb)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression<g1::dims,double,impl_p1>,FD::grid_dist_expression<prp1,g1,impl_p1>,FD::sum> exp_sum(FD::grid_dist_expression<g1::dims,double,impl_p1>(d),gb);

    return exp_sum;
}

/* \brief sum two distributed grid expression
 *
 * \param ga grid expression one
 * \param gb grid expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1 , typename exp2, typename op1>
inline FD::grid_dist_expression_op<FD::grid_dist_expression_op<exp1,exp2,op1>,FD::grid_dist_expression<exp1::gtype::dims,double,FD::NORM_EXPRESSION>,FD::sum>
operator+(const FD::grid_dist_expression_op<exp1,exp2,op1> & ga, double d)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression_op<exp1,exp2,op1>,FD::grid_dist_expression<exp1::gtype::dims,double,FD::NORM_EXPRESSION>,FD::sum> exp_sum(ga,FD::grid_dist_expression<exp1::gtype::dims,double,FD::NORM_EXPRESSION>(d));

    return exp_sum;
}

/* \brief sum two distributed grid expression
 *
 * \param ga grid expression one
 * \param gb grid expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int p1, unsigned int p2, typename g1, typename g2, unsigned int impl_p1, unsigned int impl_p2>
inline FD::grid_dist_expression_op<FD::grid_dist_expression<p1,g1,impl_p1>,FD::grid_dist_expression<p2,g2,impl_p2>,FD::sub>
operator-(const FD::grid_dist_expression<p1,g1,impl_p1> & ga, const FD::grid_dist_expression<p2,g2,impl_p2> & gb)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression<p1,g1,impl_p1>,FD::grid_dist_expression<p2,g2,impl_p2>,FD::sub> exp_sum(ga,gb);

    return exp_sum;
}

/* \brief sum two distributed grid expression
 *
 * \param ga grid expression one
 * \param gb grid expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1 , typename exp2, typename op1, unsigned int prp1, typename g1, unsigned int impl_p1>
inline FD::grid_dist_expression_op<FD::grid_dist_expression_op<exp1,exp2,op1>,FD::grid_dist_expression<prp1,g1,impl_p1>,FD::sub>
operator-(const FD::grid_dist_expression_op<exp1,exp2,op1> & ga, const FD::grid_dist_expression<prp1,g1,impl_p1> & gb)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression_op<exp1,exp2,op1>,FD::grid_dist_expression<prp1,g1,impl_p1>,FD::sub> exp_sum(ga,gb);

    return exp_sum;
}

/* \brief sum two distributed grid expression
 *
 * \param ga grid expression one
 * \param gb grid expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1 , typename exp2, typename op1, unsigned int prp1, typename g1, unsigned int impl_p1>
inline FD::grid_dist_expression_op<FD::grid_dist_expression<prp1,g1,impl_p1>,FD::grid_dist_expression_op<exp1,exp2,op1>,FD::sub>
operator-(const FD::grid_dist_expression<prp1,g1,impl_p1> & ga, const FD::grid_dist_expression_op<exp1,exp2,op1> & gb)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression<prp1,g1,impl_p1>,FD::grid_dist_expression_op<exp1,exp2,op1>,FD::sub> exp_sum(ga,gb);

    return exp_sum;
}

/* \brief minus of a distributed grid expression
 *
 * \param ga grid expression one
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1 , typename exp2_, typename op1>
inline FD::grid_dist_expression_op<FD::grid_dist_expression_op<exp1,exp2_,op1>,void,FD::subuni>
operator-(const FD::grid_dist_expression_op<exp1,exp2_,op1> & ga)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression_op<exp1,exp2_,op1>,void,FD::subuni> exp_sum(ga);

    return exp_sum;
}

template<unsigned int prp1 , typename g1, unsigned int impl_p1>
inline FD::grid_dist_expression_op<FD::grid_dist_expression<prp1,g1,impl_p1>,void,FD::subuni>
operator-(const FD::grid_dist_expression<prp1,g1,impl_p1> & ga)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression<prp1,g1,impl_p1>,void,FD::subuni> exp_sum(ga);

    return exp_sum;
}


/* \brief sum two distributed grid expression
 *
 * \param ga grid expression one
 * \param gb grid expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1 , typename exp2, typename op1, typename exp3 , typename exp4, typename op2>
inline FD::grid_dist_expression_op<FD::grid_dist_expression_op<exp1,exp2,op1>,FD::grid_dist_expression_op<exp3,exp4,op2>,FD::sub>
operator-(const FD::grid_dist_expression_op<exp1,exp2,op1> & ga, const FD::grid_dist_expression_op<exp3,exp4,op2> & gb)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression_op<exp1,exp2,op1>,FD::grid_dist_expression_op<exp3,exp4,op2>,FD::sub> exp_sum(ga,gb);

    return exp_sum;
}

/* \brief sum two distributed grid expression
 *
 * \param ga grid expression one
 * \param gb grid expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int prp1 , typename g1, unsigned int impl_p1>
inline FD::grid_dist_expression_op<FD::grid_dist_expression<prp1,g1,impl_p1>,FD::grid_dist_expression<g1::dims,double,impl_p1>,FD::sub>
operator-(const FD::grid_dist_expression<prp1,g1,impl_p1> & ga, double d)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression<prp1,g1,impl_p1>,FD::grid_dist_expression<g1::dims,double,impl_p1>,FD::sub> exp_sum(ga,FD::grid_dist_expression<g1::dims,double,impl_p1>(d));

    return exp_sum;
}

/* \brief sum two distributed grid expression
 *
 * \param ga grid expression one
 * \param gb grid expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int prp1 , typename g1, unsigned int impl_p1>
inline FD::grid_dist_expression_op<FD::grid_dist_expression<g1::dims,double,impl_p1>,FD::grid_dist_expression<prp1,g1,impl_p1>,FD::sub>
operator-(double d, const FD::grid_dist_expression<prp1,g1,impl_p1> & gb)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression<g1::dims,double,impl_p1>,FD::grid_dist_expression<prp1,g1,impl_p1>,FD::sub> exp_sum(FD::grid_dist_expression<g1::dims,double,impl_p1>(d),gb);

    return exp_sum;
}

/* \brief sum two distributed grid expression
 *
 * \param ga grid expression one
 * \param gb grid expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1 , typename exp2, typename op1>
inline FD::grid_dist_expression_op<FD::grid_dist_expression_op<exp1,exp2,op1>,FD::grid_dist_expression<exp1::gtype::dims,double,FD::NORM_EXPRESSION>,FD::sub>
operator-(const FD::grid_dist_expression_op<exp1,exp2,op1> & ga, double d)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression_op<exp1,exp2,op1>,FD::grid_dist_expression<exp1::gtype::dims,double,FD::NORM_EXPRESSION>,FD::sub> exp_sum(ga,FD::grid_dist_expression<exp1::gtype::dims,double,FD::NORM_EXPRESSION>(d));

    return exp_sum;
}

/* \brief Multiply two distributed grid expression
 *
 * \param va grid expression one
 * \param vb grid expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int p2, typename g2, unsigned int impl_p2>
inline FD::grid_dist_expression_op<FD::grid_dist_expression<g2::dims,double,impl_p2>,FD::grid_dist_expression<p2,g2,impl_p2>,FD::mul>
operator*(double d, const FD::grid_dist_expression<p2,g2,impl_p2> & vb)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression<g2::dims,double,impl_p2>,FD::grid_dist_expression<p2,g2,impl_p2>,FD::mul> exp_sum(FD::grid_dist_expression<g2::dims,double,impl_p2>(d),vb);

    return exp_sum;
}

/* \brief Multiply two distributed grid expression
 *
 * \param va grid expression one
 * \param vb grid expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int p2, typename g2, unsigned int impl_p2>
inline FD::grid_dist_expression_op<FD::grid_dist_expression<p2,g2,impl_p2>,FD::grid_dist_expression<g2::dims,double,impl_p2>,FD::mul>
operator*(const FD::grid_dist_expression<p2,g2,impl_p2> & va, double d)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression<p2,g2,impl_p2>,FD::grid_dist_expression<g2::dims,double,impl_p2>,FD::mul> exp_sum(va,FD::grid_dist_expression<g2::dims,double,impl_p2>(d));

    return exp_sum;
}

/* \brief Multiply two distributed grid expression
 *
 * \param va grid expression one
 * \param vb grid expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int p1, typename v1,unsigned int p2, typename v2, unsigned int impl_p1, unsigned int impl_p2>
inline FD::grid_dist_expression_op<FD::grid_dist_expression<p1,v1,impl_p1>,FD::grid_dist_expression<p2,v2,impl_p2>,FD::mul>
operator*(const FD::grid_dist_expression<p1,v1, impl_p1> & va, const FD::grid_dist_expression<p2,v2,impl_p2> & vb)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression<p1,v1,impl_p1>,FD::grid_dist_expression<p2,v2,impl_p2>,FD::mul> exp_sum(va,vb);

    return exp_sum;
}

/* \brief Multiply two distributed grid expression
 *
 * \param va grid expression one
 * \param vb grid expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int p1, typename v1, typename exp1, typename exp2, typename op1, unsigned int impl_p1>
inline FD::grid_dist_expression_op<FD::grid_dist_expression<p1,v1,impl_p1>,FD::grid_dist_expression_op<exp1,exp2,op1>,FD::mul>
operator*(const FD::grid_dist_expression<p1,v1,impl_p1> & va, const FD::grid_dist_expression_op<exp1,exp2,op1> & vb)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression<p1,v1,impl_p1>,FD::grid_dist_expression_op<exp1,exp2,op1>,FD::mul> exp_sum(va,vb);

    return exp_sum;
}

/* \brief Multiply two distributed grid expression
 *
 * \param va grid expression one
 * \param vb grid expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int p1, typename v1, typename exp1, typename exp2, typename op1, unsigned int impl_p1>
inline FD::grid_dist_expression_op<FD::grid_dist_expression_op<exp1,exp2,op1>,FD::grid_dist_expression<p1,v1,impl_p1>,FD::mul>
operator*(const FD::grid_dist_expression_op<exp1,exp2,op1> & va, const FD::grid_dist_expression<p1,v1,impl_p1> & vb)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression_op<exp1,exp2,op1>,FD::grid_dist_expression<p1,v1,impl_p1>,FD::mul> exp_sum(va,vb);

    return exp_sum;
}

/* \brief Multiply two distributed grid expression
 *
 * \param va grid expression one
 * \param vb grid expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1, typename exp2, typename op1, typename exp3 , typename exp4, typename op2>
inline FD::grid_dist_expression_op<FD::grid_dist_expression_op<exp1,exp2,op1>,FD::grid_dist_expression_op<exp3,exp4,op2>,FD::mul>
operator*(const FD::grid_dist_expression_op<exp1,exp2,op1> & va, const FD::grid_dist_expression_op<exp3,exp4,op2> & vb)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression_op<exp1,exp2,op1>,FD::grid_dist_expression_op<exp3,exp4,op2>,FD::mul> exp_sum(va,vb);

    return exp_sum;
}

/* \brief Multiply a distributed grid expression by a number
 *
 * \param va grid expression
 * \param d number
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1 , typename exp2, typename op1>
inline FD::grid_dist_expression_op<FD::grid_dist_expression_op<exp1,exp2,op1>,FD::grid_dist_expression<exp1::gtype::dims,double,FD::NORM_EXPRESSION>,FD::mul>
operator*(const FD::grid_dist_expression_op<exp1,exp2,op1> & va, double d)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression_op<exp1,exp2,op1>,FD::grid_dist_expression<exp1::gtype::dims,double,FD::NORM_EXPRESSION>,FD::mul> exp_sum(va,FD::grid_dist_expression<exp1::gtype::dims,double,FD::NORM_EXPRESSION>(d));

    return exp_sum;
}

/* \brief Multiply a distributed grid expression by a number
 *
 * \param d number
 * \param vb grid expression
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1 , typename exp2, typename op1>
inline FD::grid_dist_expression_op<FD::grid_dist_expression<exp1::gtype::dims,double,FD::NORM_EXPRESSION>,FD::grid_dist_expression_op<exp1,exp2,op1>,FD::mul>
operator*(double d, const FD::grid_dist_expression_op<exp1,exp2,op1> & vb)
{
    FD::grid_dist_expression_op<FD::grid_dist_expression<exp1::gtype::dims,double,FD::NORM_EXPRESSION>,FD::grid_dist_expression_op<exp1,exp2,op1>,FD::mul> exp_sum(FD::grid_dist_expression<exp1::gtype::dims,double,FD::NORM_EXPRESSION>(d),vb);

    return exp_sum;
}

#endif /* FD_EXPRESSIONS_HPP_ */