laplace_2d_fmm.hpp

Go to the documentation of this file.

 
#ifndef NIHU_LAPLACE_2D_FMM_HPP_INCLUDED
#define NIHU_LAPLACE_2D_FMM_HPP_INCLUDED
 
#include "cluster.hpp"
#include "m2l_indices.hpp"
#include "p2p.hpp"
 
#include "library/laplace_2d.hpp"
 
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/binomial.hpp>
 
#include <complex>
 
namespace NiHu
{
namespace fmm
{
 
class laplace_2d_cluster;
 
template <>
struct cluster_traits<laplace_2d_cluster>
{
    static size_t const dimension = 2;
    typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1> multipole_t;
    typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1> local_t;
};
 
class laplace_2d_cluster
    : public cluster_base<laplace_2d_cluster>
{
    typedef cluster_base<laplace_2d_cluster> base_t;
    
public:
    typedef typename base_t::multipole_t multipole_t;
    typedef typename base_t::local_t local_t;
 
    void set_expansion_length(size_t L);
 
    size_t get_expansion_length() const;
 
    multipole_t zero_multipole() const;
 
    local_t zero_local() const;
 
private:
    size_t m_expansion_length;
};
 
class laplace_2d_fmm
{
public:
    typedef laplace_2d_cluster cluster_t;
    typedef cluster_t::bounding_box_t bounding_box_t;
    typedef cluster_t::location_t location_t;
 
private:
    static std::complex<double> center2complex(cluster_t const &c);
 
public:
    template <int Nx, int Ny>
    fmm::p2p<NiHu::normal_derivative_kernel<
        NiHu::laplace_kernel<NiHu::space_2d<double> >, Nx, Ny
    > >
        create_p2p() const
    {
        typedef NiHu::normal_derivative_kernel<
            NiHu::laplace_kernel<NiHu::space_2d<double> >, Nx, Ny
        > kernel_t;
        return fmm::p2p<kernel_t>(kernel_t());
    }
 
 
    template <unsigned Ny>
    class p2m
    {
    public:
        typedef cluster_t test_input_t;
 
        typedef typename NiHu::normal_derivative_kernel<
            NiHu::laplace_kernel<NiHu::space_2d<double> >,
            0, Ny
        >::trial_input_t trial_input_t;
 
        typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1> result_t;
 
        size_t rows(test_input_t const &to) const
        {
            return to.get_expansion_length() + 1;
        }
 
        result_t operator()(test_input_t const &to, trial_input_t const & from) const
        {
            return eval(to, from, std::integral_constant<unsigned, Ny>());
        }
 
        result_t eval(test_input_t const &to,
            trial_input_t const &from,
            std::integral_constant<unsigned, 0>) const
        {
            using namespace boost::math::double_constants;
            
            auto const& y = from.get_x();
            std::complex<double> zy(y(0), y(1));
            zy -= center2complex(to);
 
            result_t res(to.get_expansion_length() + 1, 1);
            res(0) = 1.0;
            for (size_t k = 1; k <= to.get_expansion_length(); ++k)
                res(k) = -std::pow(zy, k) / double(k);
            return -res / two_pi;
        }
 
        result_t eval(test_input_t const &to,
            trial_input_t const &from,
            std::integral_constant<unsigned, 1>) const
        {
            using namespace boost::math::double_constants;
            
            auto const& y = from.get_x();
            std::complex<double> zy(y(0), y(1));
            zy -= center2complex(to);
            std::complex<double> zydny(from.get_unit_normal()(0), from.get_unit_normal()(1));
            result_t res(to.get_expansion_length() + 1, 1);
            res(0) = 0.0;
            for (size_t k = 1; k <= to.get_expansion_length(); ++k)
                res(k) = -std::pow(zy, k-1);
            res *= zydny;
            return -res / two_pi;
        }
    };
 
 
    template <unsigned Ny>
    class p2l
    {
    public:
        typedef cluster_t test_input_t;
        typedef typename NiHu::normal_derivative_kernel<
            NiHu::laplace_kernel<NiHu::space_2d<double> >,
            0, Ny
        >::trial_input_t trial_input_t;
        typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1> result_t;
 
        size_t rows(test_input_t const &to) const
        {
            return to.get_expansion_length() + 1;
        }
 
        result_t operator()(test_input_t const &to, trial_input_t const &from) const
        {
            return eval(to, from, std::integral_constant<unsigned, Ny>());
        }
 
    private:
        result_t eval(test_input_t const &to, trial_input_t const & from,
            std::integral_constant<unsigned, 0>) const
        {
            using namespace boost::math::double_constants;
 
            std::complex<double> z = 
                std::complex<double>(from.get_x()(0), from.get_x()(1)) 
                - center2complex(to);
 
            result_t res(to.get_expansion_length() + 1, 1);
            res(0) = std::log(-z);
            for (size_t k = 1; k <= to.get_expansion_length(); ++k)
                res(k) = -1. / (double(k) * std::pow(z, k));
            return -res / two_pi;
        }
 
        result_t eval(test_input_t const &to, trial_input_t const & from,
            std::integral_constant<unsigned, 1>) const
        {
            using namespace boost::math::double_constants;
            
            std::complex<double> z =
                std::complex<double>(from.get_x()(0), from.get_x()(1))
                - center2complex(to);
 
            std::complex<double> zdny(from.get_unit_normal()(0), from.get_unit_normal()(1));
 
            result_t res(to.get_expansion_length() + 1, 1);
            res(0) = 1. / z;
            for (size_t k = 1; k <= to.get_expansion_length(); ++k)
                res(k) = 1. / std::pow(z, k + 1);
            res *= zdny;
            return -res / two_pi;
        }
    };
    
 
    template <unsigned Nx>
    class l2p
    {
    public:
        typedef typename NiHu::normal_derivative_kernel<
            NiHu::laplace_kernel<NiHu::space_2d<double> >,
            Nx, 0
        >::test_input_t test_input_t;
 
        typedef cluster_t trial_input_t;
        typedef Eigen::Matrix<std::complex<double>, 1, Eigen::Dynamic> complex_result_t;
        typedef complex_result_t result_t;
 
        size_t cols(trial_input_t const &from) const
        {
            return from.get_expansion_length() + 1;
        }
 
        result_t operator()(test_input_t const &to, trial_input_t const & from) const
        {
            return eval(to, from, std::integral_constant<unsigned, Nx>());
        }
 
    private:
        result_t eval(test_input_t const &to,
            trial_input_t const &from,
            std::integral_constant<unsigned, 0>) const
        {
            std::complex<double> z = std::complex<double>(to.get_x()(0), to.get_x()(1)) - center2complex(from);
            complex_result_t res(1, from.get_expansion_length() + 1);
            for (size_t k = 0; k <= from.get_expansion_length(); ++k)
                res(k) = std::pow(z, k);
            return res;
        }
 
        result_t eval(test_input_t const &to,
            trial_input_t const &from,
            std::integral_constant<unsigned, 1>) const
        {
            std::complex<double> z = std::complex<double>(to.get_x()(0), to.get_x()(1)) - center2complex(from);
            std::complex<double> zdnx(to.get_unit_normal()(0), to.get_unit_normal()(1));
            complex_result_t res(1, from.get_expansion_length() + 1);
            for (size_t k = 0; k <= from.get_expansion_length(); ++k)
                res(k) = double(k) * std::pow(z, k-1);
            res *= zdnx;
            return res;
        }
    };
    
 
    template <unsigned Nx>
    class m2p
    {
    public:
        typedef typename NiHu::normal_derivative_kernel<
            NiHu::laplace_kernel<NiHu::space_2d<double> >,
            Nx, 0
        >::test_input_t test_input_t;
        typedef cluster_t trial_input_t;
        typedef Eigen::Matrix<std::complex<double>, 1, Eigen::Dynamic> complex_result_t;
        typedef complex_result_t result_t;
 
        size_t cols(trial_input_t const &from) const
        {
            return from.get_expansion_length() + 1;
        }
 
        result_t operator()(test_input_t const &to, trial_input_t const &from) const
        {
            return eval(to, from, std::integral_constant<unsigned, Nx>());
        }
 
    private:
        result_t eval(test_input_t const &to,
            trial_input_t const &from,
            std::integral_constant<unsigned, 0>) const
        {
            complex_result_t res(1, from.get_expansion_length() + 1);
            std::complex<double> z = std::complex<double>(to.get_x()(0), to.get_x()(1)) - center2complex(from);
            res(0) = std::log(z);
            for (size_t k = 1; k <= from.get_expansion_length(); ++k)
                res(k) = 1. / std::pow(z, k);
            return res;
        }
 
        result_t eval(test_input_t const &to,
            trial_input_t const &from,
            std::integral_constant<unsigned, 1>) const
        {
            complex_result_t res(1, from.get_expansion_length() + 1);
            std::complex<double> z = std::complex<double>(to.get_x()(0), to.get_x()(1)) - center2complex(from);
            std::complex<double> zdnx(to.get_unit_normal()(0), to.get_unit_normal()(1));
            res(0) = 1. / z;
            for (size_t k = 1; k <= from.get_expansion_length(); ++k)
                res(k) = -double(k) / std::pow(z, k + 1);
            res *= zdnx;
            return res;
        }
    };
 
 
    class m2m
    {
    public:
        typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, Eigen::Dynamic> result_t;
        typedef laplace_2d_fmm::cluster_t cluster_t;
 
        static size_t unique_idx(cluster_t const &to, cluster_t const &from);
 
        result_t operator()(cluster_t to, cluster_t from) const;
    };
 
 
    class l2l
    {
    public:
        typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, Eigen::Dynamic> result_t;
        typedef laplace_2d_fmm::cluster_t cluster_t;
 
        static size_t unique_idx(cluster_t const &to, cluster_t const &from);
 
        result_t operator()(cluster_t to, cluster_t from) const;
    };
    
 
    class m2l
    {
    public:
        typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, Eigen::Dynamic> result_t;
        typedef laplace_2d_fmm::cluster_t cluster_t;
 
        static size_t unique_idx(cluster_t const &to, cluster_t const &from);
 
        result_t operator()(cluster_t to, cluster_t from) const;
    };
    
 
    template <unsigned Ny>
    p2m<Ny> create_p2m() const
    {
        return p2m<Ny>();
    }
    
 
    template <unsigned Ny>
    p2l<Ny> create_p2l() const
    {
        return p2l<Ny>();
    }
    
 
    template <unsigned Nx>
    m2p<Nx> create_m2p() const
    {
        return m2p<Nx>();
    }
    
 
    template <unsigned Nx>
    l2p<Nx> create_l2p() const
    {
        return l2p<Nx>(); 
    }
    
 
    m2m create_m2m() const;
    
 
    l2l create_l2l() const;
    
 
    m2l create_m2l() const;
};
 
} // end of namespace fmm
} // namespace NiHu
#endif /* NIHU_LAPLACE_2D_FMM_HPP_INCLUDED */