NiHu  2.0
helmholtz_field_point.hpp
Go to the documentation of this file.
1 
7 #ifndef NIHU_HELMHOLTZ_FIELD_POINT_HPP_INCLUDED
8 #define NIHU_HELMHOLTZ_FIELD_POINT_HPP_INCLUDED
9 
10 #include "cluster_tree.hpp"
11 #include "divide.hpp"
12 #include "elem_center_iterator.hpp"
13 #include "fmm_matrix.hpp"
14 #include "fmm_precompute.hpp"
15 #include "fmm_indexed.hpp"
17 #include "fmm_integrated.hpp"
18 #include "matrix_free.hpp"
19 
20 #include "core/field.hpp"
21 #include "core/function_space.hpp"
22 
23 #include "util/timer.h"
24 #include "util/type2tag.hpp"
25 
26 #ifdef NIHU_FMM_PARALLEL
27 #include <omp.h>
28 #endif
29 
30 namespace NiHu
31 {
32 namespace fmm
33 {
34 
35 template <class Fmm, class TestSpace, class TrialSpace>
37 {
38 public:
39  typedef Fmm fmm_t;
40  typedef TestSpace test_space_t;
41  typedef TrialSpace trial_space_t;
42 
43  typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1> cvector_t;
44  typedef cvector_t excitation_t;
45  typedef cvector_t response_t;
46 
47  typedef typename fmm_t::cluster_t cluster_t;
49 
50  typedef typename tmp::deref<
51  typename tmp::begin<
52  typename trial_space_t::field_type_vector_t
53  >::type
54  >::type trial_field_t;
55 
56  enum { num_trial_dofs = trial_field_t::num_dofs };
57 
58  typedef typename tmp::deref<
59  typename tmp::begin<
60  typename test_space_t::field_type_vector_t
61  >::type
62  >::type test_field_t;
63 
66 
67  helmholtz_field_point(test_space_t const &test_space,
68  trial_space_t const &trial_space)
69  : m_test_space(test_space)
70  , m_trial_space(trial_space)
71  , m_wave_number(0.0)
72  , m_accuracy(3.0)
73  , m_far_field_order(6)
74  {
75  }
76 
77  void set_wave_number(double k)
78  {
79  m_wave_number = k;
80  }
81 
82  double get_wave_number() const
83  {
84  return m_wave_number;
85  }
86 
87  void set_accuracy(double accuracy)
88  {
89  m_accuracy = accuracy;
90  }
91 
92  void set_far_field_order(size_t far_field_order)
93  {
94  m_far_field_order = far_field_order;
95  }
96 
97 
98  response_t const &get_response() const
99  {
100  return m_response;
101  }
102 
103  void set_psurf(excitation_t const &psurf)
104  {
105  m_psurf = psurf;
106  }
107 
108  void set_qsurf(excitation_t const &qsurf)
109  {
110  m_qsurf = qsurf;
111  }
112 
113  template <class DivideDerived>
114  response_t const &eval(divide_base<DivideDerived> const &divide)
115  {
116  // create cluster tree
117  std::cout << "Building cluster tree ..." << std::endl;
118  cluster_tree_t tree(
119  create_field_center_iterator(m_trial_space.template field_begin<trial_field_t>()),
120  create_field_center_iterator(m_trial_space.template field_end<trial_field_t>()),
121  create_field_center_iterator(m_test_space.template field_begin<test_field_t>()),
122  create_field_center_iterator(m_test_space.template field_end<test_field_t>()),
123  divide.derived());
124  std::cout << tree << std::endl;
125 
126  // initialize tree data
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);
131  for (size_t c = 0; c < tree.get_n_clusters(); ++c)
132  tree[c].set_p_level_data(&fmm.get_level_data(tree[c].get_level()));
133 
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;
137  for (size_t i = 0; i < tree.get_n_levels(); ++i)
138  fmm.get_level_data(i).set_num_threads(max_num_threads);
139 #endif
140 
141  // create interaction lists
142  std::cout << "Computing interaction lists ..." << std::endl;
143  interaction_lists lists(tree);
144 
145  // create functors
147  m_far_field_order, false);
148 
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>(),
154  tree);
155 
156  auto pre_fctr = create_precompute_functor(tree, lists);
157 
158  // operator manipulations
159  std::cout << "Precomputing fmm operators ..." << std::endl;
160  auto pre_collection = create_fmm_operator_collection(
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>()),
166  fmm.create_m2m(),
167  fmm.create_l2l(),
168  fmm.create_m2l()
169  ).transform(idx_fctr).transform(pre_fctr);
170 
171  // create matrix objects
172  std::cout << "Starting assembling | SLP | DLP | " << std::endl;
173  auto combined_matrix = create_fmm_matrix(pre_collection, tree, lists);
174  std::cout << "Combined matrix assembled" << std::endl;
175 
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)
179  {
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);
184  }
185 
186  auto wc_t0 = NiHu::wc_time::tic();
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;
191 
192  combined_matrix.get_timer().print(std::cout);
193 
194  return m_response;
195  }
196 
197 private:
198  test_space_t const &m_test_space;
199  trial_space_t const &m_trial_space;
200 
201  double m_wave_number;
202  excitation_t m_psurf;
203  excitation_t m_qsurf;
204  response_t m_response;
205  double m_accuracy;
206  size_t m_far_field_order;
207 };
208 
209 template <class FmmTag, class TestSpace, class TrialSpace>
210 auto create_helmholtz_field_point(FmmTag, TestSpace const &test_space, TrialSpace const &trial_space)
211 {
213  test_space, trial_space);
214 }
215 
216 } // end of namespace fmm
217 } // end of namespace NiHu
218 
219 #endif /* NIHU_HELMHOLTZ_FIELD_POINT_HPP_INCLUDED */
type2tag.hpp
Assign a tag to a type.
fmm_matrix.hpp
Class NiHu::fmm::fmm_matrix.
NiHu::fmm::chebyshev_cluster
Cluster class of the Black Box FMM.
Definition: chebyshev_cluster.hpp:28
NiHu::fmm::cluster_tree< cluster_t >
tmp::begin
return begin iterator of a sequence
Definition: sequence.hpp:79
tmp::deref
metafunction to dereference an iterator
Definition: sequence.hpp:114
NiHu::type2tag
Metafunction assigning a tag to a type.
Definition: type2tag.hpp:17
NiHu::fmm::divide_base
Base CRTP class for cluster division.
Definition: divide.hpp:28
NiHu::fmm::cluster_tree::get_n_clusters
size_t get_n_clusters() const
return number of clusters
Definition: cluster_tree.hpp:282
NiHu::fmm::cluster_tree::get_n_levels
size_t get_n_levels() const
return number of levels
Definition: cluster_tree.hpp:253
fmm_operator_collection.hpp
Implementation of NiHu::fmm::fmm_operator_collection.
field.hpp
implementation of fields and field views
matrix_free.hpp
An Eigen::Matrix adaptor for the NiHu::fmm::fmm_matrix class.
elem_center_iterator.hpp
Classes for iterating through elements and fields.
cluster_tree.hpp
Implementation of class NiHu::fmm::cluster_tree.
NiHu::wc_time::tic
static time_point_t tic()
returns the current time point
Definition: timer.h:34
fmm_indexed.hpp
Interface for indexing FMM operators.
fmm_precompute.hpp
Operator pre-computation interface.
NiHu::fmm::create_field_center_iterator
auto create_field_center_iterator(It it)
factory function to create a field center iterator
Definition: elem_center_iterator.hpp:232
NiHu::fmm::create_fmm_matrix
fmm_matrix< P2P, P2M, P2L, M2P, L2P, M2M, L2L, M2L > create_fmm_matrix(P2P &&p2p, P2M &&p2m, P2L &&p2l, M2P &&m2p, L2P &&l2p, M2M &&m2m, L2L &&l2l, M2L &&m2l, cluster_tree< Cluster > const &tree, interaction_lists const &lists)
factory function to create an fmm_matrix object
Definition: fmm_matrix.hpp:626
NiHu::fmm::create_integrated_functor
auto create_integrated_functor(TestTag test_tag, TrialTag trial_tag, size_t quadrature_order, bool sing_check)
factory function to create an integrated functor
Definition: fmm_integrated.hpp:131
timer.h
portable implementation of wall clock and cpu timers
divide.hpp
Cluster division strategies.
fmm_integrated.hpp
Integrated operators (P2P, P2X, X2P)
function_space.hpp
declaration of class function_space
NiHu::fmm::create_fmm_operator_collection
fmm_operator_collection< Ops... > create_fmm_operator_collection(Ops &&... ops)
Create function for the collection.
Definition: fmm_operator_collection.hpp:166
NiHu::wc_time::toc
static double toc(time_point_t const &t0)
returns the time elapsed since t0
Definition: timer.h:42
NiHu::fmm::helmholtz_field_point
Definition: helmholtz_field_point.hpp:36
NiHu::fmm::interaction_lists
class storing the different interaction lists of a tree
Definition: lists.hpp:22