laplace_3d_fmm.hpp

Go to the documentation of this file.

 
#ifndef NIHU_LAPLACE_3D_FMM_HPP_INCLUDED
#define NIHU_LAPLACE_3D_FMM_HPP_INCLUDED
 
#include "cluster.hpp"
#include "m2l_indices.hpp"
#include "p2p.hpp"
 
#include "library/laplace_3d.hpp"
 
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/spherical_harmonic.hpp>
#include <boost/math/special_functions/factorials.hpp>
 
#include <Eigen/Dense>
 
#include <complex>
 
namespace NiHu
{
namespace fmm
{
 
class laplace_3d_cluster;
 
template <>
struct cluster_traits<laplace_3d_cluster>
{
private:
    typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1> cvector_t;
 
public:
    static size_t const dimension = 3;
    typedef cvector_t multipole_t;
    typedef cvector_t local_t;
};
 
class laplace_3d_cluster
    : public cluster_base<laplace_3d_cluster>
{
    typedef cluster_base<laplace_3d_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;
 
    size_t data_size() const;
 
    multipole_t zero_multipole() const;
 
    local_t zero_local() const;
 
    Eigen::Index linear_index(Eigen::Index n, Eigen::Index m) const;
 
private:
    size_t m_expansion_length;
};
 
class laplace_3d_fmm
{
public:
    typedef laplace_3d_cluster cluster_t;
    typedef cluster_t::bounding_box_t bounding_box_t;
    typedef cluster_t::location_t location_t;
 
public:
    template <int Nx, int Ny>
    fmm::p2p<NiHu::normal_derivative_kernel<
        NiHu::laplace_kernel<NiHu::space_3d<double> >, Nx, Ny
    > >
        create_p2p() const
    {
        typedef NiHu::normal_derivative_kernel<
            NiHu::laplace_kernel<NiHu::space_3d<double> >, Nx, Ny
        > kernel_t;
        return fmm::p2p<kernel_t>(kernel_t());
    }
 
    static std::complex<double> Y(size_t n, int m, double theta, double phi)
    {
        using boost::math::spherical_harmonic;
        using namespace boost::math::double_constants;
        
        std::complex<double> res = 
            spherical_harmonic(unsigned(n), std::abs(m), theta, phi) 
            / std::sqrt((2 * n + 1.0) / (4.0 * pi));
        return (m < 0) ? std::conj(res) : res;
    }
 
    static double A(int n, int m)
    {
        using boost::math::factorial;
        double res = std::pow(-1, n) / std::sqrt(
            factorial<double>(n - m) * factorial<double>(n + m));
        return res;
    }
 
 
    template <unsigned Ny>
    class p2m
    {
    public:
        typedef cluster_t test_input_t;
        typedef cluster_t::bounding_box_t::location_t x_t;
 
        typedef typename NiHu::normal_derivative_kernel<
            NiHu::laplace_kernel<NiHu::space_3d<double> >,
            0, Ny
        >::trial_input_t trial_input_t;
 
        typedef cluster_t::multipole_t result_t;
 
        size_t rows(test_input_t const &to) const
        {
            return to.data_size();
        }
 
        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;
            
            result_t res(to.data_size(), 1);
 
            x_t x = from.get_x() - to.get_bounding_box().get_center();
            double r = x.norm();
            double phi = std::atan2(x(1), x(0));
            double theta = std::acos(x(2) / r);
 
            for (int n = 0; n <= to.get_expansion_length(); ++n)
            {
                double rn = std::pow(r, n);
                for (int m = -n; m <= n; ++m)
                {
                    auto idx = to.linear_index(n, m);
                    res(idx, 0) = rn * Y(n, -m, theta, phi);
                }
            }
 
            return res / (4. * 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;
            
            result_t res = to.zero_multipole();
            return res / (4. * pi);
        }
    };
 
 
    template <unsigned Ny>
    class p2l
    {
    public:
        typedef cluster_t::bounding_box_t::location_t x_t;
        typedef cluster_t test_input_t;
        typedef typename NiHu::normal_derivative_kernel<
            NiHu::laplace_kernel<NiHu::space_3d<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.data_size();
        }
 
        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;
            
            result_t res(rows(to), 1);
 
            x_t x = from.get_x() - to.get_bounding_box().get_center();
            double r = x.norm();
            double phi = std::atan2(x(1), x(0));
            double theta = std::acos(x(2) / r);
 
            for (int n = 0; n <= to.get_expansion_length(); ++n)
            {
                double rn = std::pow(r, -(n + 1));
                for (int m = -n; m <= n; ++m)
                {
                    auto idx = to.linear_index(n, m);
                    res(idx, 0) = rn * Y(n, -m, theta, phi);
                }
            }
 
            return res / (4.0 * pi);
        }
 
        result_t eval(test_input_t const &to, trial_input_t const & from,
            std::integral_constant<unsigned, 1>) const
        {
            throw std::logic_error("Unimplemented");
            return result_t();
        }
    };
    
 
    template <unsigned Nx>
    class l2p
    {
    public:
        typedef typename NiHu::normal_derivative_kernel<
            NiHu::laplace_kernel<NiHu::space_3d<double> >,
            Nx, 0
        >::test_input_t test_input_t;
        typedef cluster_t::bounding_box_t::location_t x_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.data_size();
        }
 
        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
        {
            result_t res(1, from.data_size());
 
            x_t x = to.get_x() - from.get_bounding_box().get_center();
            double r = x.norm();
            double phi = std::atan2(x(1), x(0));
            double theta = std::acos(x(2) / r);
 
            for (int n = 0; n <= from.get_expansion_length(); ++n)
            {
                double rn = std::pow(r, n);
                for (int m = -n; m <= n; ++m)
                {
                    auto idx = from.linear_index(n, m);
                    res(0, idx) = rn * Y(n, m, theta, phi);
                }
            }
 
            return res;
        }
 
        result_t eval(test_input_t const &to,
            trial_input_t const &from,
            std::integral_constant<unsigned, 1>) const
        {
            throw std::logic_error("Unimplemented");
            return result_t();
        }
    };
    
 
    template <unsigned Nx>
    class m2p
    {
    public:
        typedef typename NiHu::normal_derivative_kernel<
            NiHu::laplace_kernel<NiHu::space_3d<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;
        typedef cluster_t::bounding_box_t::location_t x_t;
 
        size_t cols(trial_input_t const &from) const
        {
            return from.data_size();
        }
 
        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
        {
            result_t res(1, from.data_size());
 
            x_t x = to.get_x() - from.get_bounding_box().get_center();
            double r = x.norm();
            double phi = std::atan2(x(1), x(0));
            double theta = std::acos(x(2) / r);
 
            for (size_t n = 0; n <= from.get_expansion_length(); ++n)
            {
                double rn = std::pow(r, -int(n+1));
                for (int m = -int(n); m <= int(n); ++m)
                {
                    auto idx = from.linear_index(n, m);
                    res(0, idx) = rn * Y(n, m, theta, phi);
                }
            }
 
            return res;
        }
 
        result_t eval(test_input_t const &to,
            trial_input_t const &from,
            std::integral_constant<unsigned, 1>) const
        {
            throw std::logic_error("Unimplemented");
            return result_t();
        }
    };
 
 
    class m2m
    {
    public:
        typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, Eigen::Dynamic> result_t;
        typedef laplace_3d_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_3d_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_3d_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_3D_FMM_HPP_INCLUDED */