helmholtz_3d_hf_fmm.hpp

Go to the documentation of this file.

 
#ifndef NIHU_HELMHOLTZ_3D_HF_FMM_HPP_INCLUDED
#define NIHU_HELMHOLTZ_3D_HF_FMM_HPP_INCLUDED
 
#include "cluster.hpp"
#include "cluster_tree.hpp"
#include "fmm_operator.hpp"
#include "helmholtz_3d_hf_cluster.h"
#include "helmholtz_3d_hf_level_data.h"
#include "m2l_indices.hpp"
#include "operator_with_wave_number.hpp"
#include "p2p.hpp"
#include "unit_sphere.h"
 
#include "library/helmholtz_3d.hpp"
 
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/hankel.hpp>
#include <boost/math/special_functions/legendre.hpp>
 
#define L2L_SHIFT_FIRST 1
 
#include <Eigen/Dense>
#include <complex>
 
namespace NiHu
{
namespace fmm
{
 
class up_shift
{
    using cvector_t = Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1>;
 
public:
    up_shift()
        : m_level_data(nullptr)
    {
    }
 
    up_shift(cvector_t const &shift, helmholtz_3d_hf_level_data const &ld)
        : m_shift(shift)
        , m_level_data(&ld)
    {
    }
 
    cvector_t operator*(cvector_t const &rhs) const
    {
        return m_shift.array() * m_level_data->interp_up(rhs).array();
    }
 
private:
    cvector_t m_shift;
    helmholtz_3d_hf_level_data const *m_level_data;
};
 
 
class down_shift
{
    using cvector_t = Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1>;
 
public:
    down_shift()
        : m_level_data(nullptr)
    {
    }
 
    down_shift(cvector_t const &shift, helmholtz_3d_hf_level_data const &ld)
        : m_shift(shift)
        , m_level_data(&ld)
    {
    }
 
    cvector_t operator*(cvector_t const &rhs) const
    {
#if L2L_SHIFT_FIRST
        return m_level_data->interp_down(m_shift.array() * rhs.array());
#else
        return m_shift.array() * m_level_data->interp_down(rhs).array();
#endif
    }
 
private:
    cvector_t m_shift;
    helmholtz_3d_hf_level_data const *m_level_data;
};
 
 
template <class WaveNumber>
class helmholtz_3d_hf_fmm
{
public:
    using wave_number_t = WaveNumber;
    using cvector_t = Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1>;
    using cluster_t = helmholtz_3d_hf_cluster;
    using bounding_box_t = typename cluster_t::bounding_box_t;
    using location_t = typename bounding_box_t::location_t;
    using cluster_tree_t = cluster_tree<cluster_t>;
 
    using distance_dependent_kernel_t = NiHu::helmholtz_kernel<NiHu::space_3d<>, wave_number_t>;
 
    helmholtz_3d_hf_fmm(wave_number_t const &k)
        : m_wave_number(k)
        , m_accuracy(3.0)
    {
    }
 
    static size_t compute_expansion_length(double drel, double C)
    {
        using namespace boost::math::double_constants;
        auto kd = two_pi * drel;
        return size_t(std::ceil(kd + C * std::log(kd + pi)));
    }
 
    void set_accuracy(double accuracy)
    {
        m_accuracy = accuracy;
    }
 
    void init_level_data(cluster_tree_t const &tree)
    {
        using namespace boost::math::double_constants;
        double lambda = two_pi / std::real(m_wave_number);
        m_level_data_vector.clear();
        m_level_data_vector.resize(tree.get_n_levels());
 
        // set the expansion length for each level
        for (size_t i = 0; i < tree.get_n_levels(); ++i)
        {
            auto &ld = m_level_data_vector[i];
            // get the diameter from the first cluster on the level
            auto idx = tree.level_begin(i);
            auto d = tree[idx].get_bounding_box().get_diameter();
            auto L = compute_expansion_length(d / lambda, m_accuracy);
            ld.set_expansion_length(L);
        }
 
        // compute interpolation matrices
        for (size_t i = 0; i < tree.get_n_levels(); ++i)
        {
            auto &ld = m_level_data_vector[i];
            auto const &Sto = ld.get_unit_sphere();
            if (i >= 3)
                ld.set_interp_dn(interpolator(m_level_data_vector[i - 1].get_unit_sphere(), Sto));
            if (i < tree.get_n_levels() - 1 && i >= 2)
                ld.set_interp_up(interpolator(m_level_data_vector[i + 1].get_unit_sphere(), Sto));
        }
    } // end of function init_level_data
 
    helmholtz_3d_hf_level_data const &get_level_data(size_t idx) const
    {
        return m_level_data_vector[idx];
    }
 
    helmholtz_3d_hf_level_data &get_level_data(size_t idx)
    {
        return m_level_data_vector[idx];
    }
 
    class m2m
        : public operator_with_wave_number<wave_number_t>
        , public fmm_operator<m2m_tag>
    {
    public:
        using wn_base_t = operator_with_wave_number<wave_number_t>;
        using result_t = up_shift;
        using cluster_t = helmholtz_3d_hf_fmm::cluster_t;
 
        m2m(wave_number_t const &wave_number)
            : wn_base_t(wave_number)
        {
        }
 
        static size_t unique_idx(cluster_t const &to, cluster_t const &from)
        {
            return bounding_box_t::dist2idx(
                from.get_bounding_box().get_center(),
                to.get_bounding_box().get_center());
        }
 
        result_t operator()(cluster_t const &to, cluster_t const &from) const
        {
            using namespace std::complex_literals;
            auto const &Sto = to.get_level_data().get_unit_sphere();
            location_t D = to.get_bounding_box().get_center()
                - from.get_bounding_box().get_center();
            auto const &k = this->get_wave_number();
            cvector_t shift = exp((-1.i * k * Sto.get_s().transpose() * D).array());
            return result_t(shift, to.get_level_data());
        }
    };
 
 
    class l2l
        : public operator_with_wave_number<wave_number_t>
        , public fmm_operator<l2l_tag>
    {
    public:
        using wn_base_t = operator_with_wave_number<WaveNumber>;
        using result_t = down_shift;
        using cluster_t = helmholtz_3d_hf_fmm::cluster_t;
 
        l2l(wave_number_t const &wave_number)
            : wn_base_t(wave_number)
        {
        }
 
        static size_t unique_idx(cluster_t const &to, cluster_t const &from)
        {
            return bounding_box_t::dist2idx(
                to.get_bounding_box().get_center(),
                from.get_bounding_box().get_center());
        }
 
        result_t operator()(cluster_t const &to, cluster_t const &from) const
        {
            using namespace std::complex_literals;
#if L2L_SHIFT_FIRST
            auto const &Sfrom = from.get_level_data().get_unit_sphere();
#else
            auto const &Sto = to.get_level_data().get_unit_sphere();
#endif
            location_t D = to.get_bounding_box().get_center()
                - from.get_bounding_box().get_center();
            auto const &k = this->get_wave_number();
#if L2L_SHIFT_FIRST
            cvector_t shift = exp((-1.i * k * Sfrom.get_s().transpose() * D).array());
#else
            cvector_t shift = exp((-1.i * k * Sto.get_s().transpose() * D).array());
#endif
            return result_t(shift, to.get_level_data());
        }
    };
 
 
    template <int Ny>
    class p2m
        : public operator_with_wave_number<wave_number_t>
        , public fmm_operator<p2m_tag>
    {
    public:
        using wn_base_t = operator_with_wave_number<wave_number_t>;
 
        using test_input_t = cluster_t;
        using trial_input_t = typename NiHu::normal_derivative_kernel<
            distance_dependent_kernel_t, 0, Ny
        >::trial_input_t;
        using result_t = cvector_t;
 
        p2m(wave_number_t const &wave_number)
            : wn_base_t(wave_number)
        {
        }
 
        size_t rows(test_input_t const &to) const
        {
            return to.get_level_data().get_unit_sphere().get_s().cols();
        }
 
        result_t operator()(test_input_t const &to, trial_input_t const &tri) const
        {
            return eval(to, tri, std::integral_constant<int, Ny>());
        }
 
    private:
        result_t eval(test_input_t const &to, trial_input_t const &tri,
            std::integral_constant<int, 0>) const
        {
            using namespace std::complex_literals;
            auto const &Y = to.get_bounding_box().get_center();
            auto const &y = tri.get_x();
            auto const &s = to.get_level_data().get_unit_sphere().get_s();
            auto const &k = this->get_wave_number();
            location_t d = Y - y;
            return Eigen::exp(-1.i * k * (s.transpose() * d).array());
        }
 
        result_t eval(test_input_t const &to, trial_input_t const &tri,
            std::integral_constant<int, 1>) const
        {
            using namespace std::complex_literals;
            auto const &Y = to.get_bounding_box().get_center();
            auto const &y = tri.get_x();
            auto const &s = to.get_level_data().get_unit_sphere().get_s();
            auto const &k = this->get_wave_number();
            location_t d = Y - y;
            return Eigen::exp(-1.i * k * (s.transpose() * d).array())
                * ((1.i * k) * (s.transpose() * tri.get_unit_normal()).array());
        }
    };
 
    template <int Ny>
    class p2l
        : public operator_with_wave_number<wave_number_t>
        , public fmm_operator<p2l_tag>
    {
    public:
        using wn_base_t =operator_with_wave_number<wave_number_t>;
 
        using test_input_t = cluster_t;
        using trial_input_t = typename NiHu::normal_derivative_kernel<
            distance_dependent_kernel_t, 0, Ny
        >::trial_input_t;
        using result_t = cvector_t;
 
        p2l(wave_number_t const &wave_number)
            : wn_base_t(wave_number)
        {
        }
 
        size_t rows(test_input_t const &to) const
        {
            return to.get_level_data().get_unit_sphere().get_s().cols();
        }
 
        result_t operator()(test_input_t const &to, trial_input_t const &tri) const
        {
            return eval(to, tri, std::integral_constant<int, Ny>());
        }
 
    private:
        result_t eval(test_input_t const &to, trial_input_t const &tri,
            std::integral_constant<int, 0>) const
        {
            throw std::logic_error("Unimplemented p2l operator");
        }
 
        result_t eval(test_input_t const &to, trial_input_t const &tri,
            std::integral_constant<int, 1>) const
        {
            throw std::logic_error("Unimplemented p2l operator");
        }
    };
 
 
    template <int Nx>
    class l2p
        : public operator_with_wave_number<wave_number_t>
        , public fmm_operator<l2p_tag>
    {
    public:
        using wn_base_t = operator_with_wave_number<wave_number_t>;
 
        using trial_input_t = cluster_t;
        using test_input_t = typename NiHu::normal_derivative_kernel<
            distance_dependent_kernel_t, Nx, 0
        >::test_input_t;
        using result_t = Eigen::Matrix<std::complex<double>, 1, Eigen::Dynamic>;
 
        l2p(wave_number_t const &wave_number)
            : wn_base_t(wave_number)
        {
        }
 
        size_t cols(trial_input_t const &from) const
        {
            return from.get_level_data().get_unit_sphere().get_s().cols();
        }
 
        result_t operator()(test_input_t const &tsi, trial_input_t const &from) const
        {
            return eval(tsi, from, std::integral_constant<int, Nx>());
        }
 
    private:
        result_t eval(test_input_t const &tsi, trial_input_t const &from,
            std::integral_constant<int, 0>) const
        {
            using namespace std::complex_literals;
            auto const &X = from.get_bounding_box().get_center();
            auto const &x = tsi.get_x();
            auto const &S = from.get_level_data().get_unit_sphere();
            auto const &s = S.get_s();
            auto const &w = S.get_w();
            auto const &k = this->get_wave_number();
            location_t d = x - X;
            return Eigen::exp(-1.i * k * (s.transpose() * d).array()) * w.array();
        }
 
        result_t eval(test_input_t const &tsi, trial_input_t const &from,
            std::integral_constant<int, 1>) const
        {
            using namespace std::complex_literals;
            auto const &X = from.get_bounding_box().get_center();
            auto const &x = tsi.get_x();
            auto const &S = from.get_level_data().get_unit_sphere();
            auto const &s = S.get_s();
            auto const &w = S.get_w();
            auto const &k = this->get_wave_number();
            location_t d = x - X;
            return Eigen::exp(-1.i * k * (s.transpose() * d).array())
                * (-1.i * k) * (s.transpose() * tsi.get_unit_normal()).array()
                * w.array();
        }
    };
 
 
    template <int Nx>
    class m2p
        : public operator_with_wave_number<wave_number_t>
        , public fmm_operator<m2p_tag>
    {
    public:
        using wn_base_t = operator_with_wave_number<wave_number_t>;
 
        using trial_input_t = cluster_t;
        using test_input_t = typename NiHu::normal_derivative_kernel<
            distance_dependent_kernel_t, Nx, 0
        >::test_input_t;
        using result_t = Eigen::Matrix<std::complex<double>, 1, Eigen::Dynamic>;
 
        m2p(wave_number_t const &wave_number)
            : wn_base_t(wave_number)
        {
        }
 
        size_t cols(trial_input_t const &from) const
        {
            return from.get_level_data().get_unit_sphere().get_s().cols();
        }
 
        result_t operator()(test_input_t const &tsi, trial_input_t const &from) const
        {
            return eval(tsi, from, std::integral_constant<int, Nx>());
        }
 
    private:
        result_t eval(test_input_t const &tsi, trial_input_t const &from,
            std::integral_constant<int, 0>) const
        {
            throw std::logic_error("Unimplemented m2p operator");
        }
 
        result_t eval(test_input_t const &tsi, trial_input_t const &from,
            std::integral_constant<int, 1>) const
        {
            throw std::logic_error("Unimplemented m2p operator");
        }
    };
 
    class m2l
        : public operator_with_wave_number<wave_number_t>
        , public fmm_operator<m2l_tag>
    {
    public:
        using wn_base_t = operator_with_wave_number<wave_number_t>;
        using result_t = Eigen::DiagonalMatrix<std::complex<double>, Eigen::Dynamic>;
        using cluster_t = helmholtz_3d_hf_fmm::cluster_t;
 
        m2l(wave_number_t const &wave_number)
            : wn_base_t(wave_number)
        {
        }
 
        static size_t unique_idx(cluster_t const &to, cluster_t const &from)
        {
            return m2l_indices<3>::eval(to.get_bounding_box(), from.get_bounding_box());
        }
 
        result_t operator()(cluster_t const &to, cluster_t const &from) const
        {
            auto const &k = this->get_wave_number();
 
            auto const &X = to.get_bounding_box().get_center();
            auto const &Y = from.get_bounding_box().get_center();
            location_t Dvec = X - Y;
 
            auto L = to.get_level_data().get_expansion_length();
            auto const &s = to.get_level_data().get_unit_sphere().get_s();
 
            return m2l_matrix_impl(Dvec, s, k, L).asDiagonal();
        }
 
    private:
        static cvector_t m2l_matrix_impl(location_t const &Dvec,
            Eigen::Matrix<double, 3, Eigen::Dynamic> const &s,
            wave_number_t const &k,
            size_t L)
        {
            using namespace boost::math::double_constants;
            using namespace std::complex_literals;
    
            double D = Dvec.norm();
            location_t Uvec = Dvec / D;
            auto N = s.cols();
 
            Eigen::Matrix<double, 1, Eigen::Dynamic> x = Uvec.transpose() * s;
            auto z = -k * D;
 
            cvector_t M2L = cvector_t::Zero(N, 1);
 
            for (size_t l = 0; l <= L; ++l)
            {
                auto h = boost::math::sph_hankel_1(l, z);
                auto c = double(2 * l + 1) * std::pow(1.i, l) * h;
                for (int i = 0; i < s.cols(); ++i)
                    M2L(i) += c * boost::math::legendre_p(int(l), x(0, i));
            }
            M2L *= -1.i * k / (4.0 * pi) / (4.0 * pi);
 
            return M2L;
        }
    };
 
    template <int Ny>
    using p2m_type = p2m<Ny>;
 
    template <int Ny>
    using p2l_type = p2l<Ny>;
 
    template <int Nx>
    using m2p_type = m2p<Nx>;
 
    template <int Nx>
    using l2p_type = l2p<Nx>;
 
    template <int Nx, int Ny>
    using p2p_type = fmm::p2p<
            NiHu::normal_derivative_kernel<
            distance_dependent_kernel_t, Nx, Ny
            >
        >;
 
    template <int Ny>
    p2m_type<Ny> create_p2m() const
    {
        return p2m<Ny>(m_wave_number);
    }
 
    template <int Ny>
    p2l_type<Ny> create_p2l() const
    {
        return p2l<Ny>(m_wave_number);
    }
 
    template <int Nx>
    l2p_type<Nx> create_l2p() const
    {
        return l2p<Nx>(m_wave_number);
    }
 
    template <int Nx>
    m2p_type<Nx> create_m2p() const
    {
        return m2p<Nx>(m_wave_number);
    }
 
    template <int Nx, int Ny>
    p2p_type<Nx, Ny> create_p2p() const
    {
        typedef NiHu::normal_derivative_kernel<
            distance_dependent_kernel_t, Nx, Ny
        > kernel_t;
        return p2p_type<Nx, Ny>(kernel_t(m_wave_number));
    }
 
    m2m create_m2m() const
    {
        return m2m(m_wave_number);
    }
 
    l2l create_l2l() const
    {
        return l2l(m_wave_number);
    }
 
    m2l create_m2l() const
    {
        return m2l(m_wave_number);
    }
 
private:
    wave_number_t m_wave_number;
    std::vector<helmholtz_3d_hf_level_data> m_level_data_vector;
    double m_accuracy; // C
};
 
 
} // end of namespace fmm
} // namespace NiHu
 
#endif /* NIHU_HELMHOLTZ_3D_HF_FMM_HPP_INCLUDED */