7 #ifndef NIHU_LAPLACE_2D_FMM_HPP_INCLUDED
8 #define NIHU_LAPLACE_2D_FMM_HPP_INCLUDED
14 #include "library/laplace_2d.hpp"
16 #include <boost/math/constants/constants.hpp>
17 #include <boost/math/special_functions/binomial.hpp>
27 class laplace_2d_cluster;
34 static size_t const dimension = 2;
36 typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1>
multipole_t;
38 typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1>
local_t;
67 size_t m_expansion_length;
85 static std::complex<double> center2complex(
cluster_t const &c);
92 template <
int Nx,
int Ny>
107 template <
unsigned Ny>
116 >::trial_input_t trial_input_t;
118 typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1> result_t;
125 result_t operator()(
test_input_t const &to, trial_input_t
const & from)
const
127 return eval(to, from, std::integral_constant<unsigned, Ny>());
131 trial_input_t
const &from,
132 std::integral_constant<unsigned, 0>)
const
134 using namespace boost::math::double_constants;
136 auto const& y = from.get_x();
137 std::complex<double> zy(y(0), y(1));
138 zy -= center2complex(to);
143 res(k) = -std::pow(zy, k) / double(k);
144 return -res / two_pi;
148 trial_input_t
const &from,
149 std::integral_constant<unsigned, 1>)
const
151 using namespace boost::math::double_constants;
153 auto const& y = from.get_x();
154 std::complex<double> zy(y(0), y(1));
155 zy -= center2complex(to);
156 std::complex<double> zydny(from.get_unit_normal()(0), from.get_unit_normal()(1));
160 res(k) = -std::pow(zy, k-1);
162 return -res / two_pi;
169 template <
unsigned Ny>
177 >::trial_input_t trial_input_t;
178 typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1> result_t;
185 result_t operator()(
test_input_t const &to, trial_input_t
const &from)
const
187 return eval(to, from, std::integral_constant<unsigned, Ny>());
191 result_t eval(
test_input_t const &to, trial_input_t
const & from,
192 std::integral_constant<unsigned, 0>)
const
194 using namespace boost::math::double_constants;
196 std::complex<double> z =
197 std::complex<double>(from.get_x()(0), from.get_x()(1))
198 - center2complex(to);
201 res(0) = std::log(-z);
203 res(k) = -1. / (double(k) * std::pow(z, k));
204 return -res / two_pi;
207 result_t eval(
test_input_t const &to, trial_input_t
const & from,
208 std::integral_constant<unsigned, 1>)
const
210 using namespace boost::math::double_constants;
212 std::complex<double> z =
213 std::complex<double>(from.get_x()(0), from.get_x()(1))
214 - center2complex(to);
216 std::complex<double> zdny(from.get_unit_normal()(0), from.get_unit_normal()(1));
221 res(k) = 1. / std::pow(z, k + 1);
223 return -res / two_pi;
230 template <
unsigned Nx>
237 >::test_input_t test_input_t;
240 typedef Eigen::Matrix<std::complex<double>, 1, Eigen::Dynamic> complex_result_t;
241 typedef complex_result_t result_t;
248 result_t operator()(test_input_t
const &to,
trial_input_t const & from)
const
250 return eval(to, from, std::integral_constant<unsigned, Nx>());
254 result_t eval(test_input_t
const &to,
256 std::integral_constant<unsigned, 0>)
const
258 std::complex<double> z = std::complex<double>(to.get_x()(0), to.get_x()(1)) - center2complex(from);
261 res(k) = std::pow(z, k);
265 result_t eval(test_input_t
const &to,
267 std::integral_constant<unsigned, 1>)
const
269 std::complex<double> z = std::complex<double>(to.get_x()(0), to.get_x()(1)) - center2complex(from);
270 std::complex<double> zdnx(to.get_unit_normal()(0), to.get_unit_normal()(1));
273 res(k) = double(k) * std::pow(z, k-1);
282 template <
unsigned Nx>
289 >::test_input_t test_input_t;
291 typedef Eigen::Matrix<std::complex<double>, 1, Eigen::Dynamic> complex_result_t;
292 typedef complex_result_t result_t;
299 result_t operator()(test_input_t
const &to,
trial_input_t const &from)
const
301 return eval(to, from, std::integral_constant<unsigned, Nx>());
305 result_t eval(test_input_t
const &to,
307 std::integral_constant<unsigned, 0>)
const
310 std::complex<double> z = std::complex<double>(to.get_x()(0), to.get_x()(1)) - center2complex(from);
311 res(0) = std::log(z);
313 res(k) = 1. / std::pow(z, k);
317 result_t eval(test_input_t
const &to,
319 std::integral_constant<unsigned, 1>)
const
322 std::complex<double> z = std::complex<double>(to.get_x()(0), to.get_x()(1)) - center2complex(from);
323 std::complex<double> zdnx(to.get_unit_normal()(0), to.get_unit_normal()(1));
326 res(k) = -double(k) / std::pow(z, k + 1);
337 typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, Eigen::Dynamic> result_t;
350 typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, Eigen::Dynamic> result_t;
363 typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, Eigen::Dynamic> result_t;
375 template <
unsigned Ny>
385 template <
unsigned Ny>
395 template <
unsigned Nx>
405 template <
unsigned Nx>