7 #ifndef NIHU_HELMHOLTZ_FIELD_POINT_HPP_INCLUDED
8 #define NIHU_HELMHOLTZ_FIELD_POINT_HPP_INCLUDED
26 #ifdef NIHU_FMM_PARALLEL
35 template <
class Fmm,
class TestSpace,
class TrialSpace>
40 typedef TestSpace test_space_t;
41 typedef TrialSpace trial_space_t;
43 typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1> cvector_t;
44 typedef cvector_t excitation_t;
45 typedef cvector_t response_t;
52 typename trial_space_t::field_type_vector_t
54 >::type trial_field_t;
56 enum { num_trial_dofs = trial_field_t::num_dofs };
60 typename test_space_t::field_type_vector_t
68 trial_space_t
const &trial_space)
69 : m_test_space(test_space)
70 , m_trial_space(trial_space)
73 , m_far_field_order(6)
77 void set_wave_number(
double k)
82 double get_wave_number()
const
87 void set_accuracy(
double accuracy)
89 m_accuracy = accuracy;
92 void set_far_field_order(
size_t far_field_order)
94 m_far_field_order = far_field_order;
98 response_t
const &get_response()
const
103 void set_psurf(excitation_t
const &psurf)
108 void set_qsurf(excitation_t
const &qsurf)
113 template <
class Div
ideDerived>
117 std::cout <<
"Building cluster tree ..." << std::endl;
124 std::cout << tree << std::endl;
127 std::cout <<
"Initializing fmm object and level data ..." << std::endl;
128 fmm_t fmm(m_wave_number);
129 fmm.set_accuracy(m_accuracy);
130 fmm.init_level_data(tree);
132 tree[c].set_p_level_data(&fmm.get_level_data(tree[c].get_level()));
134 #if NIHU_FMM_PARALLEL
135 auto max_num_threads = omp_get_max_threads();
136 std::cout <<
"Expanding to " << max_num_threads <<
" threads" << std::endl;
138 fmm.get_level_data(i).set_num_threads(max_num_threads);
142 std::cout <<
"Computing interaction lists ..." << std::endl;
147 m_far_field_order,
false);
149 auto idx_fctr = create_indexed_functor(
150 m_test_space.template field_begin<test_field_t>(),
151 m_test_space.template field_end<test_field_t>(),
152 m_trial_space.template field_begin<trial_field_t>(),
153 m_trial_space.template field_end<trial_field_t>(),
156 auto pre_fctr = create_precompute_functor(tree, lists);
159 std::cout <<
"Precomputing fmm operators ..." << std::endl;
161 src_concatenate(int_fctr(fmm.template create_p2p<0, 0>()), int_fctr(fmm.template create_p2p<0, 1>())),
162 src_concatenate(int_fctr(fmm.template create_p2m<0>()), int_fctr(fmm.template create_p2m<1>())),
163 src_concatenate(int_fctr(fmm.template create_p2l<0>()), int_fctr(fmm.template create_p2l<1>())),
164 int_fctr(fmm.template create_m2p<0>()),
165 int_fctr(fmm.template create_l2p<0>()),
169 ).transform(idx_fctr).transform(pre_fctr);
172 std::cout <<
"Starting assembling | SLP | DLP | " << std::endl;
174 std::cout <<
"Combined matrix assembled" << std::endl;
176 std::cout <<
"Computing MVP " << std::endl;
177 cvector_t xct(m_psurf.rows() + m_qsurf.rows(), 1);
178 for (
int i = 0; i < m_psurf.rows()/num_trial_dofs; ++i)
180 xct.segment(2 * i * num_trial_dofs, num_trial_dofs) =
181 -m_qsurf.segment(i * num_trial_dofs, num_trial_dofs);
182 xct.segment((2 * i + 1) * num_trial_dofs, num_trial_dofs) =
183 m_psurf.segment(i * num_trial_dofs, num_trial_dofs);
187 auto cpu_t0 = NiHu::cpu_time::tic();
188 m_response = combined_matrix * xct;
189 std::cout <<
"MVP wall clock time: " <<
NiHu::wc_time::toc(wc_t0) <<
" s" << std::endl;
190 std::cout <<
"MVP CPU time: " << NiHu::cpu_time::toc(cpu_t0) <<
" s" << std::endl;
192 combined_matrix.get_timer().print(std::cout);
198 test_space_t
const &m_test_space;
199 trial_space_t
const &m_trial_space;
201 double m_wave_number;
202 excitation_t m_psurf;
203 excitation_t m_qsurf;
204 response_t m_response;
206 size_t m_far_field_order;
209 template <
class FmmTag,
class TestSpace,
class TrialSpace>
210 auto create_helmholtz_field_point(FmmTag, TestSpace
const &test_space, TrialSpace
const &trial_space)
213 test_space, trial_space);