Laplacian.hpp

/*
 * Laplacian.hpp
 *
 *  Created on: Oct 5, 2015
 *      Author: i-bird
 */

#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_LAPLACIAN_HPP_
#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_LAPLACIAN_HPP_

#include "FD_util_include.hpp"
#include "util/util_num.hpp"
#include "FiniteDifference/eq.hpp"

template<typename arg, typename Sys_eqs, unsigned int impl=CENTRAL>
class Lap
{
    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,
                             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";
    }
};

template<typename arg, typename Sys_eqs>
class Lap<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)
    {
        // for each dimension
        for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
        {
            long int old_val = kmap.getKeyRef().get(i);
            kmap.getKeyRef().set_d(i, kmap.getKeyRef().get(i) + 1);
            arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[i]/spacing[i]);
            kmap.getKeyRef().set_d(i,old_val);

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

            arg::value(g_map,kmap,gs,spacing,cols, - 2.0 * coeff/spacing[i]/spacing[i]);
        }
    }
};


template<typename arg, typename Sys_eqs>
class Lap<arg,Sys_eqs,CENTRAL_SYM>
{
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)
    {
        // for each dimension
        for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
        {
            long int old_val = kmap.getKeyRef().get(i);
            kmap.getKeyRef().set_d(i, kmap.getKeyRef().get(i) + 2);
            arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[i]/spacing[i]/4.0);
            kmap.getKeyRef().set_d(i,old_val);

            old_val = kmap.getKeyRef().get(i);
            kmap.getKeyRef().set_d(i, kmap.getKeyRef().get(i) - 2);
            arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[i]/spacing[i]/4.0);
            kmap.getKeyRef().set_d(i,old_val);

            arg::value(g_map,kmap,gs,spacing,cols, - 2.0 * coeff/spacing[i]/spacing[i]/4.0);
        }
    }
};

#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_LAPLACIAN_HPP_ */