7 #ifndef NIHU_HELMHOLTZ_BURTON_MILLER_SOLVER_HPP_INCLUDED
8 #define NIHU_HELMHOLTZ_BURTON_MILLER_SOLVER_HPP_INCLUDED
18 #include "preconditioner.hpp"
25 #include <Eigen/IterativeLinearSolvers>
28 #ifdef NIHU_FMM_PARALLEL
40 template <
class Fmm,
class TrialSpace>
45 typedef TrialSpace trial_space_t;
48 typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1> excitation_t;
49 typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1> response_t;
56 typename trial_space_t::field_type_vector_t
58 >::type trial_field_t;
65 : m_trial_space(trial_space)
66 , m_test_space(
dirac(m_trial_space))
69 , m_far_field_order(6)
76 void set_restart(
size_t restart)
81 void set_tolerance(
double tolerance)
83 m_tolerance = tolerance;
86 void set_wave_number(
double k)
91 double get_wave_number()
const
96 response_t
const &get_response()
const
101 void set_excitation(excitation_t
const &xct)
106 void set_accuracy(
double accuracy)
108 m_accuracy = accuracy;
111 void set_far_field_order(
size_t order)
113 m_far_field_order = order;
119 template <
class Div
ideDerived>
123 std::cout <<
"Building cluster tree ..." << std::endl;
130 std::cout <<
"Tree:\n" << tree << std::endl;
133 std::cout <<
"Instantiating fmm object ..." << std::endl;
134 fmm_t fmm(m_wave_number);
136 fmm.set_accuracy(m_accuracy);
139 std::cout <<
"Initializing tree data ..." << std::endl;
140 fmm.init_level_data(tree);
142 tree[c].set_p_level_data(&fmm.get_level_data(tree[c].get_level()));
144 #if NIHU_FMM_PARALLEL
145 auto max_num_threads = omp_get_max_threads();
146 std::cout <<
"Expanding to " << max_num_threads <<
" threads" << std::endl;
148 fmm.get_level_data(i).set_num_threads(max_num_threads);
152 std::cout <<
"Computing interaction lists..." << std::endl;
157 m_far_field_order,
true);
159 auto idx_fctr = create_indexed_functor(
160 m_test_space.template field_begin<test_field_t>(),
161 m_test_space.template field_end<test_field_t>(),
162 m_trial_space.template field_begin<trial_field_t>(),
163 m_trial_space.template field_end<trial_field_t>(),
166 auto pre_fctr = create_precompute_functor(tree, lists);
169 std::complex<double> alpha(0.0, -1.0 / m_wave_number);
174 int_fctr(fmm.template create_p2p<0, 1>())
175 + alpha * int_fctr(fmm.template create_p2p<1, 1>())
177 int_fctr(fmm.template create_p2m<1>()),
178 int_fctr(fmm.template create_p2l<1>()),
179 int_fctr(fmm.template create_m2p<0>())
180 + alpha * int_fctr(fmm.template create_m2p<1>()),
181 int_fctr(fmm.template create_l2p<0>())
182 + alpha * int_fctr(fmm.template create_l2p<1>()),
189 int_fctr(fmm.template create_p2p<0, 0>())
190 + alpha * int_fctr(fmm.template create_p2p<1, 0>())
192 int_fctr(fmm.template create_p2m<0>()),
193 int_fctr(fmm.template create_p2l<0>())
197 auto lhs_cix_collection = lhs_collection.transform(idx_fctr);
198 auto rhs_cix_collection = rhs_collection.transform(idx_fctr);
201 std::cout <<
"Precomputing fmm operators ..." << std::endl;
202 auto lhs_pre_collection = lhs_cix_collection.transform(pre_fctr);
205 std::cout <<
"Assembling rhs matrix ..." << std::endl;
207 pre_fctr(rhs_cix_collection.get(
p2p_tag())),
208 rhs_cix_collection.get(
p2m_tag()),
209 rhs_cix_collection.get(
p2l_tag()),
210 lhs_pre_collection.get(
m2p_tag()),
211 lhs_pre_collection.get(
l2p_tag()),
212 lhs_pre_collection.get(
m2m_tag()),
213 lhs_pre_collection.get(
l2l_tag()),
214 lhs_pre_collection.get(
m2l_tag()),
218 std::cout <<
"Computing rhs ..." << std::endl;
219 response_t rhs = slp_matrix * m_excitation;
222 std::cout <<
"Assembling lhs matrix ..." << std::endl;
228 std::cout <<
"Starting iterative solution ..." << std::endl;
231 Eigen::GMRES<decltype(M), Eigen::IdentityPreconditioner > solver(M);
232 solver.setTolerance(m_tolerance);
233 solver.set_restart(m_restart);
234 m_response = solver.solve(rhs);
235 m_iters = solver.iterations();
240 size_t get_iterations()
const
246 trial_space_t
const &m_trial_space;
247 test_space_t
const &m_test_space;
249 double m_wave_number;
251 size_t m_far_field_order;
252 excitation_t m_excitation;
253 response_t m_response;
259 template <
class FmmTag,
class TrialSpace>
260 auto create_helmholtz_burton_miller_solver(FmmTag, TrialSpace
const& trial_space)
262 return helmholtz_burton_miller_solver<typename tag2type<FmmTag>::type, TrialSpace>(trial_space);