Derivative.hpp

/*
 * Derivative.hpp
 *
 *  Created on: Oct 5, 2015
 *      Author: Pietro Incardona
 */

#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_DERIVATIVE_HPP_
#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_DERIVATIVE_HPP_


#include "util/mathutil.hpp"
#include "Vector/map_vector.hpp"
#include "Grid/comb.hpp"
#include "FiniteDifference/util/common.hpp"
#include "util/util_num.hpp"
#include <unordered_map>
#include "FD_util_include.hpp"

template<unsigned int d, typename Field, typename Sys_eqs, unsigned int impl=CENTRAL>
class D
{
    inline static void value(const grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
    {
        std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
    }

    inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos,
                                                      const grid_sm<Sys_eqs::dims,void> & gs,
                                                      const comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
    {
        std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";

        return pos;
    }
};

template<unsigned int d, typename arg, typename Sys_eqs>
class D<d,arg,Sys_eqs,CENTRAL>
{
    public:

    inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map,
                             grid_dist_key_dx<Sys_eqs::dims> & kmap,
                             const grid_sm<Sys_eqs::dims,void> & gs,
                             typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
                             std::unordered_map<long int,typename Sys_eqs::stype > & cols,
                             typename Sys_eqs::stype coeff)
    {
        // if the system is staggered the CENTRAL derivative is equivalent to a forward derivative
        if (is_grid_staggered<Sys_eqs>::value())
        {
            D<d,arg,Sys_eqs,BACKWARD>::value(g_map,kmap,gs,spacing,cols,coeff);
            return;
        }

        long int old_val = kmap.getKeyRef().get(d);
        kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 1);
        arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[d]/2.0 );
        kmap.getKeyRef().set_d(d,old_val);


        old_val = kmap.getKeyRef().get(d);
        kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 1);
        arg::value(g_map,kmap,gs,spacing,cols,-coeff/spacing[d]/2.0 );
        kmap.getKeyRef().set_d(d,old_val);
    }


    inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos,
                                                      const grid_sm<Sys_eqs::dims,void> & gs,
                                                      const comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
    {
        auto arg_pos = arg::position(pos,gs,s_pos);
        if (is_grid_staggered<Sys_eqs>::value())
        {
            if (arg_pos.get(d) == -1)
            {
                arg_pos.set_d(d,0);
                return arg_pos;
            }
            else
            {
                arg_pos.set_d(d,-1);
                return arg_pos;
            }
        }

        return arg_pos;
    }
};


template<unsigned int d, typename arg, typename Sys_eqs>
class D<d,arg,Sys_eqs,CENTRAL_B_ONE_SIDE>
{
public:

    static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map,
                      grid_dist_key_dx<Sys_eqs::dims> & kmap,
                      const grid_sm<Sys_eqs::dims,void> & gs,
                      typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
                      std::unordered_map<long int,typename Sys_eqs::stype > & cols,
                      typename Sys_eqs::stype coeff)
    {
#ifdef SE_CLASS1
        if (Sys_eqs::boundary[d] == PERIODIC)
            std::cerr << __FILE__ << ":" << __LINE__ << " error, it make no sense use one sided derivation with periodic boundary, please use CENTRAL\n";
#endif

        grid_key_dx<Sys_eqs::dims> pos = g_map.getGKey(kmap);

        if (pos.get(d) == (long int)gs.size(d)-1 )
        {
            arg::value(g_map,kmap,gs,spacing,cols,1.5*coeff/spacing[d]);

            long int old_val = kmap.getKeyRef().get(d);
            kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 1);
            arg::value(g_map,kmap,gs,spacing,cols,-2.0*coeff/spacing[d]);
            kmap.getKeyRef().set_d(d,old_val);

            old_val = kmap.getKeyRef().get(d);
            kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 2);
            arg::value(g_map,kmap,gs,spacing,cols,0.5*coeff/spacing[d]);
            kmap.getKeyRef().set_d(d,old_val);
        }
        else if (pos.get(d) == 0)
        {
            arg::value(g_map,kmap,gs,spacing,cols,-1.5*coeff/spacing[d]);

            long int old_val = kmap.getKeyRef().get(d);
            kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 1);
            arg::value(g_map,kmap,gs,spacing,cols,2.0*coeff/spacing[d]);
            kmap.getKeyRef().set_d(d,old_val);

            old_val = kmap.getKeyRef().get(d);
            kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 2);
            arg::value(g_map,kmap,gs,spacing,cols,-0.5*coeff/spacing[d]);
            kmap.getKeyRef().set_d(d,old_val);
        }
        else
        {
            long int old_val = kmap.getKeyRef().get(d);
            kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 1);
            arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[d]);
            kmap.getKeyRef().set_d(d,old_val);

            old_val = kmap.getKeyRef().get(d);
            kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 1);
            arg::value(g_map,kmap,gs,spacing,cols,-coeff/spacing[d]);
            kmap.getKeyRef().set_d(d,old_val);
        }
    }

    inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, const comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
    {
        return arg::position(pos,gs,s_pos);
    }
};


template<unsigned int d, typename arg, typename Sys_eqs>
class D<d,arg,Sys_eqs,FORWARD>
{
    public:

    inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map,
                             grid_dist_key_dx<Sys_eqs::dims> & kmap,
                             const grid_sm<Sys_eqs::dims,void> & gs,
                             typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
                             std::unordered_map<long int,typename Sys_eqs::stype > & cols,
                             typename Sys_eqs::stype coeff)
    {

        long int old_val = kmap.getKeyRef().get(d);
        kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 1);
        arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[d]);
        kmap.getKeyRef().set_d(d,old_val);

        // backward
        arg::value(g_map,kmap,gs,spacing,cols,-coeff/spacing[d]);
    }


    inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos,
                                                      const grid_sm<Sys_eqs::dims,void> & gs,
                                                      const comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
    {
        return arg::position(pos,gs,s_pos);
    }
};

template<unsigned int d, typename arg, typename Sys_eqs>
class D<d,arg,Sys_eqs,BACKWARD>
{
    public:

    inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map,
                             grid_dist_key_dx<Sys_eqs::dims> & kmap,
                             const grid_sm<Sys_eqs::dims,void> & gs,
                             typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
                             std::unordered_map<long int,typename Sys_eqs::stype > & cols,
                             typename Sys_eqs::stype coeff)
    {

        long int old_val = kmap.getKeyRef().get(d);
        kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 1);
        arg::value(g_map,kmap,gs,spacing,cols,-coeff/spacing[d]);
        kmap.getKeyRef().set_d(d,old_val);

        // forward
        arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[d]);
    }


    inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, const comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
    {
        return arg::position(pos,gs,s_pos);
    }
};

#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_DERIVATIVE_HPP_ */