NiHu  2.0
convected_helmholtz_3d_kernel.hpp
1 #ifndef NIHU_CONVECTED_HELMHOLTZ_3D_KERNEL_HPP_INCLUDED
2 #define NIHU_CONVECTED_HELMHOLTZ_3D_KERNEL_HPP_INCLUDED
3 
4 #include "../core/gaussian_quadrature.hpp"
5 #include "../core/kernel.hpp"
6 #include "guiggiani_1992.hpp"
7 #include "location_normal.hpp"
8 #include "../core/space.hpp"
9 
10 #include "wave_number_kernel.hpp"
11 
12 #include "Eigen/Geometry"
13 
14 namespace NiHu
15 {
16 
23 template <class WaveNumber, class Space>
24 class convected_helmholtz_kernel
25  : public wave_number_kernel<WaveNumber>
26 {
27 public:
28  typedef Space space_t;
29  typedef typename space_t::location_t x_t;
30  typedef WaveNumber wave_number_t;
31  typedef Eigen::Matrix3d rotation_matrix_t;
32  typedef Eigen::DiagonalMatrix<double, 3> scaling_matrix_t;
33 
40  convected_helmholtz_kernel(wave_number_t const &k, x_t const &mach_number_vector)
41  : wave_number_kernel<WaveNumber>(k)
42  , m_mach_number_vector(mach_number_vector)
43  , m_normalized_mach_number_vector(m_mach_number_vector.normalized())
44  , m_mach_number(m_mach_number_vector.norm())
45  , m_kappa(this->get_wave_number() / (1 - m_mach_number*m_mach_number))
46  {
47  // Compute the rotation matrix
48  Eigen::Vector3d axis = m_normalized_mach_number_vector.cross(Eigen::Vector3d::UnitX());
49  if (m_mach_number > 1e-10 && axis.norm() > 1e-10) {
50  double angle = std::acos(m_normalized_mach_number_vector(0));
51  m_rot = Eigen::AngleAxis<double>(angle, axis.normalized()).matrix();
52  } else {
53  // Mach number very small, or already aligned with X axis (axis.norm very small)
54  m_rot = rotation_matrix_t::Identity();
55  }
56 
57  // Set up scaling matrix
58  double Q = 1.0 - m_mach_number * m_mach_number;
59  m_scal.setZero();
60  m_scal.diagonal() << 1.0, std::sqrt(Q), std::sqrt(Q);
61  }
62 
63  wave_number_t const &get_kappa() const
64  {
65  return m_kappa;
66  }
67 
68  x_t const &get_mach_number_vector() const
69  {
70  return m_mach_number_vector;
71  }
72 
73  double const &get_mach_number() const
74  {
75  return m_mach_number;
76  }
77 
78  x_t const &get_normalized_mach_number_vector() const
79  {
80  return m_normalized_mach_number_vector;
81  }
82 
83  rotation_matrix_t const &get_rotation_matrix() const
84  {
85  return m_rot;
86  }
87 
88  scaling_matrix_t const &get_scaling_matrix() const
89  {
90  return m_scal;
91  }
92 
93 private:
94  x_t const m_mach_number_vector;
95  x_t const m_normalized_mach_number_vector;
96  double const m_mach_number;
97  wave_number_t const m_kappa;
98  rotation_matrix_t m_rot;
99  scaling_matrix_t m_scal;
100 };
101 
102 
108 template <class WaveNumber>
109 class convected_helmholtz_3d_SLP_kernel;
110 
111 template <class WaveNumber>
112 struct kernel_traits<convected_helmholtz_3d_SLP_kernel<WaveNumber>>
113 {
114  typedef location_input_3d test_input_t;
115  typedef location_normal_input_3d trial_input_t;
116  typedef std::complex<double> result_t;
117  typedef gauss_family_tag quadrature_family_t;
118  static bool const is_symmetric = false;
119 
120  typedef asymptotic::inverse<1> far_field_behaviour_t;
121 
122  static bool const is_singular = true;
123 };
124 
125 template <class WaveNumber>
126 struct singular_kernel_traits<convected_helmholtz_3d_SLP_kernel<WaveNumber>>
127 {
128  typedef asymptotic::inverse<1> singularity_type_t;
129  static unsigned const singular_quadrature_order = 7;
130 };
131 
132 
133 template <class WaveNumber>
134 class convected_helmholtz_3d_SLP_kernel
135  : public kernel_base<convected_helmholtz_3d_SLP_kernel<WaveNumber>>
136  , public convected_helmholtz_kernel<WaveNumber, space_3d<double>>
137 {
138 public:
139  typedef kernel_traits<convected_helmholtz_3d_SLP_kernel<WaveNumber>> traits_t;
140  typedef typename wave_number_kernel<WaveNumber>::wave_number_t wave_number_t;
141  typedef typename traits_t::test_input_t test_input_t;
142  typedef typename traits_t::trial_input_t trial_input_t;
143  typedef typename test_input_t::x_t x_t;
144  typedef typename traits_t::result_t result_t;
145 
146  convected_helmholtz_3d_SLP_kernel(wave_number_t const &k, x_t const &mach_number_vector)
147  : convected_helmholtz_kernel<WaveNumber, space_3d<double> >(k, mach_number_vector)
148  {
149  }
150 
151  result_t operator()(x_t const &x, x_t const &y, x_t const &ny) const
152  {
153  using boost::math::double_constants::pi;
154  using namespace std::complex_literals;
155 
156  double const &Ma = this->get_mach_number();
157  x_t const &Ma_vec = this->get_mach_number_vector();
158  wave_number_t const &kappa = this->get_kappa();
159 
160  x_t r_vec = y - x;
161  double M_dot_r = Ma_vec.dot(r_vec);
162  double R = std::sqrt((1.0 - Ma*Ma) * r_vec.squaredNorm() + (M_dot_r * M_dot_r));
163  result_t G = std::exp(-1i*kappa*(R+M_dot_r)) / (4 * pi * R);
164  double M_dot_ny = Ma_vec.dot(ny);
165  return G * (1.0 - M_dot_ny*M_dot_ny);
166  }
167 
168  result_t operator()(test_input_t const &x, trial_input_t const &y) const
169  {
170  return (*this)(x.get_x(), y.get_x(), y.get_unit_normal());
171  }
172 };
173 
174 
175 
176 template <class WaveNumber>
177 class convected_helmholtz_3d_DLP_kernel;
178 
179 template <class WaveNumber>
180 struct kernel_traits<convected_helmholtz_3d_DLP_kernel<WaveNumber>>
181 {
182  typedef location_input_3d test_input_t;
183  typedef location_normal_input_3d trial_input_t;
184  typedef std::complex<double> result_t;
185  typedef gauss_family_tag quadrature_family_t;
186  static bool const is_symmetric = false;
187 
188  typedef asymptotic::inverse<2> far_field_behaviour_t;
189 
190  static bool const is_singular = true;
191 };
192 
193 template <class WaveNumber>
194 struct singular_kernel_traits<convected_helmholtz_3d_DLP_kernel<WaveNumber>>
195 {
196  typedef asymptotic::inverse<1> singularity_type_t; // remains weakly singular because of normal derivative
197  static unsigned const singular_quadrature_order = 7;
198 };
199 
200 template <class WaveNumber>
201 class convected_helmholtz_3d_DLP_kernel
202  : public kernel_base<convected_helmholtz_3d_DLP_kernel<WaveNumber>>
203  , public convected_helmholtz_kernel<WaveNumber, space_3d<double> >
204 {
205 public:
206  typedef kernel_traits<convected_helmholtz_3d_DLP_kernel<WaveNumber>> traits_t;
207  typedef typename wave_number_kernel<WaveNumber>::wave_number_t wave_number_t;
208  typedef typename traits_t::test_input_t test_intput_t;
209  typedef typename traits_t::trial_input_t trial_input_t;
210  typedef typename trial_input_t::x_t x_t;
211  typedef typename traits_t::result_t result_t;
212 
213  convected_helmholtz_3d_DLP_kernel(wave_number_t const &k, x_t const &mach_number_vector)
214  : convected_helmholtz_kernel<WaveNumber, space_3d<double> >(k, mach_number_vector)
215  {
216  }
217 
218  result_t operator()(x_t const &x, x_t const &y, x_t const &ny) const
219  {
220  using boost::math::double_constants::pi;
221  using namespace std::complex_literals;
222 
223  wave_number_t const &k = this->get_wave_number();
224  wave_number_t const& kappa = this->get_kappa();
225  double const &M = this->get_mach_number();
226  x_t const &Ma_vec = this->get_mach_number_vector();
227 
228  x_t r_vec = y - x;
229 
230  double M_dot_r = Ma_vec.dot(r_vec);
231  double M_dot_ny = Ma_vec.dot(ny);
232  double R = std::sqrt((1-M*M) * r_vec.squaredNorm() + (M_dot_r * M_dot_r));
233 
234  std::complex<double> G = std::exp(-1i*kappa*(R + M_dot_r)) / (4. * pi * R);
235  Eigen::Matrix<std::complex<double>, 3, 1> gradG = -G * ((1. + 1i*kappa*R) / (R * R) * ((1-M*M)*r_vec + Ma_vec * M_dot_r) + 1i*kappa*Ma_vec);
236  std::complex<double> dGn = ny.dot(gradG); // order of arguments is important, dot conjugates 1st arg!
237  std::complex<double> dGM = Ma_vec.dot(gradG); //this->get_normalized_mach_number_vector().dot(gradG);
238 
239  return dGn + (2i*k*G - dGM) * M_dot_ny; // ny.dot(this->get_normalized_mach_number_vector());
240  //return dGn - dGM * M_dot_ny;
241  }
242 
243  result_t operator()(test_intput_t const &x, trial_input_t const &y) const
244  {
245  return (*this)(x.get_x(), y.get_x(), y.get_unit_normal());
246  }
247 };
248 
249 
250 template <class WaveNumber>
251 class convected_helmholtz_3d_DLP_tan_kernel;
252 
253 template <class WaveNumber>
254 struct kernel_traits<convected_helmholtz_3d_DLP_tan_kernel<WaveNumber>>
255 {
256  typedef location_input_3d test_input_t;
257  typedef location_normal_input_3d trial_input_t;
258  typedef std::complex<double> result_t;
259  typedef gauss_family_tag quadrature_family_t;
260  static bool const is_symmetric = false;
261 
262  typedef asymptotic::inverse<1> far_field_behaviour_t;
263 
264  static bool const is_singular = true;
265 };
266 
267 template <class WaveNumber>
268 struct singular_kernel_traits<convected_helmholtz_3d_DLP_tan_kernel<WaveNumber>>
269 {
270  typedef asymptotic::inverse<1> singularity_type_t; // remains weakly singular because of normal derivative
271  static unsigned const singular_quadrature_order = 7;
272 };
273 
274 template <class WaveNumber>
275 class convected_helmholtz_3d_DLP_tan_kernel
276  : public kernel_base<convected_helmholtz_3d_DLP_tan_kernel<WaveNumber>>
277  , public convected_helmholtz_kernel<WaveNumber, space_3d<double> >
278 {
279 public:
280  typedef kernel_traits<convected_helmholtz_3d_DLP_tan_kernel<WaveNumber>> traits_t;
281  typedef typename wave_number_kernel<WaveNumber>::wave_number_t wave_number_t;
282  typedef typename traits_t::test_input_t test_intput_t;
283  typedef typename traits_t::trial_input_t trial_input_t;
284  typedef typename trial_input_t::x_t x_t;
285  typedef typename traits_t::result_t result_t;
286 
287  convected_helmholtz_3d_DLP_tan_kernel(wave_number_t const &k, x_t const &mach_number_vector)
288  : convected_helmholtz_kernel<WaveNumber, space_3d<double> >(k, mach_number_vector)
289  {
290  }
291 
292  result_t operator()(x_t const &x, x_t const &y, x_t const &ny) const
293  {
294  using boost::math::double_constants::pi;
295  using namespace std::complex_literals;
296 
297  wave_number_t const &kappa = this->get_kappa();
298  double const &M = this->get_mach_number();
299  x_t const &Ma_vec = this->get_mach_number_vector();
300 
301  x_t rvec = y - x;
302 
303  double M_dot_r = Ma_vec.dot(rvec);
304  double M_dot_ny = Ma_vec.dot(ny);
305  double R = std::sqrt((1-M*M) * rvec.squaredNorm() + (M_dot_r * M_dot_r));
306  std::complex<double> G = std::exp(-1i*kappa*(R+M_dot_r)) / (4. * pi * R);
307  return G * M_dot_ny; // ny.dot(this->get_normalized_mach_number_vector());
308  }
309 
310  result_t operator()(test_intput_t const &x, trial_input_t const &y) const
311  {
312  return (*this)(x.get_x(), y.get_x(), y.get_unit_normal());
313  }
314 };
315 
316 
317 
318 template <class WaveNumber>
319 class convected_helmholtz_3d_DLPt_kernel;
320 
321 template <class WaveNumber>
322 struct kernel_traits<convected_helmholtz_3d_DLPt_kernel<WaveNumber>>
323 {
324  typedef location_normal_input_3d test_input_t;
325  //typedef location_input_3d trial_input_t;
326  typedef location_normal_input_3d trial_input_t;
327  typedef std::complex<double> result_t;
328  typedef gauss_family_tag quadrature_family_t;
329  static bool const is_symmetric = false;
330 
331  typedef asymptotic::inverse<2> far_field_behaviour_t;
332 
333  static bool const is_singular = true;
334 };
335 
336 template <class WaveNumber>
337 struct singular_kernel_traits<convected_helmholtz_3d_DLPt_kernel<WaveNumber>>
338 {
339  typedef asymptotic::inverse<2> singularity_type_t;
340  static unsigned const singular_quadrature_order = 7;
341  typedef convected_helmholtz_3d_DLPt_kernel<WaveNumber> singular_core_t;
342 };
343 
344 template <class WaveNumber>
345 class convected_helmholtz_3d_DLPt_kernel
346  : public kernel_base<convected_helmholtz_3d_DLPt_kernel<WaveNumber>>
347  , public convected_helmholtz_kernel<WaveNumber, space_3d<double> >
348 {
349 public:
350  typedef kernel_traits<convected_helmholtz_3d_DLPt_kernel<WaveNumber>> traits_t;
351  typedef typename wave_number_kernel<WaveNumber>::wave_number_t wave_number_t;
352  typedef typename traits_t::test_input_t test_intput_t;
353  typedef typename traits_t::trial_input_t trial_input_t;
354  typedef typename trial_input_t::x_t x_t;
355  typedef typename traits_t::result_t result_t;
356 
357  convected_helmholtz_3d_DLPt_kernel(wave_number_t const &k, x_t const &mach_number_vector)
358  : convected_helmholtz_kernel<WaveNumber, space_3d<double> >(k, mach_number_vector)
359  {
360  }
361 
362  result_t operator()(x_t const &x, x_t const& nx, x_t const &y, x_t const& ny) const
363  {
364  using boost::math::double_constants::pi;
365  using namespace std::complex_literals;
366 
367  wave_number_t const& kappa = this->get_kappa();
368  double const &M = this->get_mach_number();
369  x_t const &Ma_vec = this->get_mach_number_vector();
370 
371  x_t r_vec = y - x;
372 
373  double M_dot_r = Ma_vec.dot(r_vec);
374  double M_dot_nx = Ma_vec.dot(nx);
375 
376  double R = std::sqrt((1-M*M) * r_vec.squaredNorm() + (M_dot_r * M_dot_r));
377 
378  std::complex<double> G = std::exp(-1i*kappa*(R + M_dot_r)) / (4. * pi * R);
379  Eigen::Matrix<std::complex<double>, 3, 1> gradG = -G * ((1. + 1i*kappa*R) / (R * R) * ((1-M*M)*r_vec + Ma_vec * M_dot_r) + 1i*kappa*Ma_vec);
380  std::complex<double> dGn = -nx.dot(gradG); // order of arguments is important, dot conjugates 1st arg!
381  std::complex<double> dGM = -Ma_vec.dot(gradG);
382 
383  double M_dot_ny = Ma_vec.dot(ny);
384 
385  return (dGn - dGM*M_dot_nx) * (1 - M_dot_ny*M_dot_ny);
386  }
387 
388  result_t operator()(test_intput_t const &x, trial_input_t const &y) const
389  {
390  return (*this)(x.get_x(), x.get_unit_normal(), y.get_x(), y.get_unit_normal());
391  }
392 };
393 
394 
395 template <class WaveNumber>
396 class convected_helmholtz_3d_DLPt_n_kernel;
397 
398 template <class WaveNumber>
399 struct kernel_traits<convected_helmholtz_3d_DLPt_n_kernel<WaveNumber>>
400 {
401  typedef location_normal_input_3d test_input_t;
402  //typedef location_input_3d trial_input_t;
403  typedef location_normal_input_3d trial_input_t;
404  typedef std::complex<double> result_t;
405  typedef gauss_family_tag quadrature_family_t;
406  static bool const is_symmetric = false;
407 
408  typedef asymptotic::inverse<2> far_field_behaviour_t;
409 
410  static bool const is_singular = true;
411 };
412 
413 
414 
415 template <class WaveNumber>
416 struct singular_kernel_traits<convected_helmholtz_3d_DLPt_n_kernel<WaveNumber>>
417 {
418  typedef asymptotic::inverse<2> singularity_type_t;
419  static unsigned const singular_quadrature_order = 7;
420  typedef convected_helmholtz_3d_DLPt_n_kernel<WaveNumber> singular_core_t;
421 };
422 
423 template <class WaveNumber>
424 class convected_helmholtz_3d_DLPt_n_kernel
425  : public kernel_base<convected_helmholtz_3d_DLPt_n_kernel<WaveNumber>>
426  , public convected_helmholtz_kernel<WaveNumber, space_3d<double> >
427 {
428 public:
429  typedef kernel_traits<convected_helmholtz_3d_DLPt_n_kernel<WaveNumber>> traits_t;
430  typedef typename wave_number_kernel<WaveNumber>::wave_number_t wave_number_t;
431  typedef typename traits_t::test_input_t test_intput_t;
432  typedef typename traits_t::trial_input_t trial_input_t;
433  typedef typename trial_input_t::x_t x_t;
434  typedef typename traits_t::result_t result_t;
435 
436  convected_helmholtz_3d_DLPt_n_kernel(wave_number_t const &k, x_t const &mach_number_vector)
437  : convected_helmholtz_kernel<WaveNumber, space_3d<double> >(k, mach_number_vector)
438  {
439  }
440 
441  result_t operator()(x_t const &x, x_t const& nx, x_t const &y, x_t const& ny) const
442  {
443  using boost::math::double_constants::pi;
444  using namespace std::complex_literals;
445 
446  wave_number_t const& kappa = this->get_kappa();
447  double const &M = this->get_mach_number();
448  x_t const &Ma_vec = this->get_mach_number_vector();
449 
450  x_t r_vec = y - x;
451 
452  double M_dot_r = Ma_vec.dot(r_vec);
453  double M_dot_nx = Ma_vec.dot(nx);
454 
455  double R = std::sqrt((1-M*M) * r_vec.squaredNorm() + (M_dot_r * M_dot_r));
456 
457  std::complex<double> G = std::exp(-1i*kappa*(R + M_dot_r)) / (4. * pi * R);
458  Eigen::Matrix<std::complex<double>, 3, 1> gradG = -G * ((1. + 1i*kappa*R) / (R * R) * ((1-M*M)*r_vec + Ma_vec * M_dot_r) + 1i*kappa*Ma_vec);
459  std::complex<double> dGn = -nx.dot(gradG); // order of arguments is important, dot conjugates 1st arg!
460  std::complex<double> dGM = -Ma_vec.dot(gradG);
461 
462  double M_dot_ny = Ma_vec.dot(ny);
463 
464  return (dGn - dGM*M_dot_nx) * M_dot_ny;
465  }
466 
467  result_t operator()(test_intput_t const &x, trial_input_t const &y) const
468  {
469  return (*this)(x.get_x(), x.get_unit_normal(), y.get_x(), y.get_unit_normal());
470  }
471 };
472 
473 
477 template <class WaveNumber>
478 class convected_helmholtz_3d_HSP_kernel;
479 
480 template <class WaveNumber>
481 struct kernel_traits<convected_helmholtz_3d_HSP_kernel<WaveNumber>>
482 {
483  typedef location_normal_input_3d test_input_t;
484  typedef location_normal_input_3d trial_input_t;
485  typedef std::complex<double> result_t;
486  typedef gauss_family_tag quadrature_family_t;
487  static bool const is_symmetric = false;
488 
489  typedef asymptotic::inverse<3> far_field_behaviour_t;
490 
491  static bool const is_singular = true;
492 };
493 
494 template <class WaveNumber>
495 struct singular_kernel_traits<convected_helmholtz_3d_HSP_kernel<WaveNumber>>
496 {
497  typedef asymptotic::inverse<3> singularity_type_t;
498  static unsigned const singular_quadrature_order = 7;
499  typedef convected_helmholtz_3d_HSP_kernel<WaveNumber> singular_core_t;
500 };
501 
502 
503 template <class WaveNumber>
504 class convected_helmholtz_3d_HSP_kernel
505  : public kernel_base<convected_helmholtz_3d_HSP_kernel<WaveNumber>>
506  , public convected_helmholtz_kernel<WaveNumber, space_3d<double> >
507 {
508 public:
509  typedef kernel_traits<convected_helmholtz_3d_HSP_kernel<WaveNumber>> traits_t;
510  typedef typename wave_number_kernel<WaveNumber>::wave_number_t wave_number_t;
511  typedef typename traits_t::test_input_t test_intput_t;
512  typedef typename traits_t::trial_input_t trial_input_t;
513  typedef typename trial_input_t::x_t x_t;
514  typedef typename traits_t::result_t result_t;
515 
516  convected_helmholtz_3d_HSP_kernel(wave_number_t const &k, x_t const &mach_number_vector)
517  : convected_helmholtz_kernel<WaveNumber, space_3d<double> >(k, mach_number_vector)
518  {
519  }
520 
521  result_t operator()(x_t const &x, x_t const& nx, x_t const &y, x_t const& ny) const
522  {
523  using boost::math::double_constants::pi;
524  using namespace std::complex_literals;
525 
526  wave_number_t const& kappa = this->get_kappa();
527  double const &M = this->get_mach_number();
528  x_t const &Ma_vec = this->get_mach_number_vector();
529 
530  x_t r_vec = y - x;
531  double r = r_vec.norm();
532 
533  double Q = 1.0 - M*M;
534  double M_dot_r = Ma_vec.dot(r_vec);
535  double M_dot_nx = Ma_vec.dot(nx);
536  double M_dot_ny = Ma_vec.dot(ny);
537  double R = std::sqrt((1-M*M) * r_vec.squaredNorm() + (M_dot_r * M_dot_r));
538 
539  double r_per_R = r/R;
540  double nx_r_per_r = r_vec.dot(nx)/r;
541  double ny_r_per_r = r_vec.dot(ny)/r;
542 
543  std::complex<double> G = std::exp(-1i*kappa*(R + M_dot_r)) / (4. * pi * R);
544 
545  return G*Q*(
546  - Q * r_per_R*r_per_R * (3.0 + 3i*kappa*R - (kappa*R)*(kappa*R)) / (R*R) * nx_r_per_r * ny_r_per_r
547  - Q / R * r_per_R * 1i*kappa* (1.0 + 1i*kappa*R) * (M_dot_nx * ny_r_per_r + M_dot_ny * nx_r_per_r)
548  + 1.0 / (R*R) * (1.0 + 1i*kappa*R) * (ny.dot(nx) - M_dot_nx * M_dot_ny)
549  + Q * kappa*kappa * M_dot_nx * M_dot_ny
550  );
551 
552  }
553 
554  result_t operator()(test_intput_t const &x, trial_input_t const &y) const
555  {
556  return (*this)(x.get_x(), x.get_unit_normal(), y.get_x(), y.get_unit_normal());
557  }
558 };
559 
560 
561 template <class WaveNumber>
562 class polar_laurent_coeffs<convected_helmholtz_3d_DLPt_kernel<WaveNumber> >
563 {
564 public:
565  typedef convected_helmholtz_3d_DLPt_kernel<WaveNumber> kernel_t;
566 
567  template <class guiggiani>
568  static void eval(guiggiani &obj)
569  {
570  using boost::math::double_constants::pi;
571  using namespace std::complex_literals;
572 
573  kernel_t const &kernel = obj.get_kernel();
574 
575  double const M = kernel.get_mach_number();
576  double const Q = 1 - M*M;
577 
578  auto const& Mav = kernel.get_mach_number_vector();
579 
580  auto const& D = kernel.get_scaling_matrix();
581  auto const& T = kernel.get_rotation_matrix();
582 
583  auto const& A_vec = obj.get_rvec_series(_1());
584  auto A_hat_vec = D * T * A_vec;
585  double A = A_vec.norm(); // should be 1 because of Rong's improvement
586  double A_hat = A_hat_vec.norm();
587 
588  auto const &J0_vec = obj.get_Jvec_series(_0());
589  double J0 = std::sqrt(J0_vec.dot(J0_vec));
590 
591  auto d0 = A_vec.normalized();
592  auto const &nx = obj.get_n0();
593 
594  auto const &N0 = obj.get_shape_series(_0());
595 
596  double const M_n = nx.dot(Mav);
597 
598  obj.set_laurent_coeff(_m1(),
599  Q / (4.0 * pi) * A / std::pow(A_hat, 3) * nx.dot(d0) * J0 * N0 * (1 - M_n*M_n)
600  );
601  }
602 };
603 
604 
605 template <class WaveNumber>
606 class polar_laurent_coeffs<convected_helmholtz_3d_DLPt_n_kernel<WaveNumber> >
607 {
608 public:
609  typedef convected_helmholtz_3d_DLPt_n_kernel<WaveNumber> kernel_t;
610 
611  template <class guiggiani>
612  static void eval(guiggiani &obj)
613  {
614  using boost::math::double_constants::pi;
615  using namespace std::complex_literals;
616 
617  kernel_t const &kernel = obj.get_kernel();
618 
619  double const M = kernel.get_mach_number();
620  double const Q = 1 - M*M;
621 
622  auto const& Mav = kernel.get_mach_number_vector();
623 
624  auto const& D = kernel.get_scaling_matrix();
625  auto const& T = kernel.get_rotation_matrix();
626 
627  auto const& A_vec = obj.get_rvec_series(_1());
628  auto A_hat_vec = D * T * A_vec;
629  double A = A_vec.norm(); // should be 1 because of Rong's improvement
630  double A_hat = A_hat_vec.norm();
631 
632  auto const &J0_vec = obj.get_Jvec_series(_0());
633  double J0 = std::sqrt(J0_vec.dot(J0_vec));
634 
635  auto d0 = A_vec.normalized();
636  auto const &nx = obj.get_n0();
637 
638  auto const &N0 = obj.get_shape_series(_0());
639 
640  double const M_n = nx.dot(Mav);
641 
642  obj.set_laurent_coeff(_m1(),
643  Q / (4.0 * pi) * A / std::pow(A_hat, 3) * nx.dot(d0) * J0 * N0 * M_n
644  );
645  }
646 };
647 
648 
649 template <class WaveNumber>
650 class polar_laurent_coeffs<convected_helmholtz_3d_HSP_kernel<WaveNumber> >
651 {
652 public:
653  typedef WaveNumber wave_number_t;
654  typedef convected_helmholtz_3d_HSP_kernel<wave_number_t> kernel_t;
655 
656  template <class guiggiani>
657  static void eval(guiggiani &obj)
658  {
659  // Could be like this:
660  // typedef typename guiggiani::kernel_t kernel_t
661  using boost::math::double_constants::pi;
662  using namespace std::complex_literals;
663 
664  kernel_t const &kernel = obj.get_kernel();
665 
666  double const M = kernel.get_mach_number();
667  double const Q = 1. - M*M;
668  wave_number_t const &kappa = kernel.get_kappa();
669  auto const& Mav = kernel.get_mach_number_vector();
670 
671  auto const& D = kernel.get_scaling_matrix();
672  auto const& T = kernel.get_rotation_matrix();
673 
674  auto const& A_vec = obj.get_rvec_series(_1());
675  auto const& B_vec = obj.get_rvec_series(_2());
676 
677  auto A_hat_vec = D * T * A_vec;
678  auto B_hat_vec = D * T * B_vec;
679 
680  double A = A_vec.norm(); // should be 1 because of Rong's improvement
681  double A_hat = A_hat_vec.norm();
682 
683  auto d0 = A_vec.normalized();
684  auto d1 = (B_vec - d0 * d0.dot(B_vec)) / A;
685 
686  auto const &nx = obj.get_n0();
687 
688  auto const &J0 = obj.get_Jvec_series(_0());
689  auto const &J1 = obj.get_Jvec_series(_1());
690 
691  auto const &N0 = obj.get_shape_series(_0());
692  auto const &N1 = obj.get_shape_series(_1());
693 
694 
695  obj.set_laurent_coeff(_m1(),
696  -Q*Q / (4.0 * pi) * 3.0 * A*A / std::pow(A_hat, 5) * (d0.dot(nx)) * (d1.dot(J0) + d0.dot(J1)) * N0
697  -Q*Q / (4.0 * pi) * A / std::pow(A_hat, 3) * (1i*kappa) * (J0.dot(Mav)) * (d0.dot(nx)) * N0
698  +Q / (4.0 * pi) / std::pow(A_hat, 3) * (J1.dot(nx) - nx.dot(Mav) * J1.dot(Mav)) * N0
699  +Q / (4.0 * pi) / std::pow(A_hat, 3) * (J0.dot(nx) - nx.dot(Mav) * J0.dot(Mav)) * N1
700  +Q / (4.0 * pi) * (-3. * A_hat_vec.dot(B_hat_vec) / std::pow(A_hat, 5) - 1i*kappa*A_vec.dot(Mav)/std::pow(A_hat, 3)) * (J0.dot(nx) - nx.dot(Mav) * J0.dot(Mav)) * N0
701  );
702 
703  obj.set_laurent_coeff(_m2(),
704  Q / (4.0 * pi) / std::pow(A_hat, 3) * (J0.dot(nx) - nx.dot(Mav) * J0.dot(Mav)) * N0
705  );
706 
707  }
708 };
709 
710 
711 }
712 
713 #endif // NIHU_CONVECTED_HELMHOLTZ_3D_KERNEL_HPP_INCLUDED
NiHu::convected_helmholtz_kernel
Base class of kernels of the convected Helmholtz equation.
Definition: convected_helmholtz_3d_kernel.hpp:24
NiHu::convected_helmholtz_kernel::convected_helmholtz_kernel
convected_helmholtz_kernel(wave_number_t const &k, x_t const &mach_number_vector)
Construct a new convected helmholtz kernel object.
Definition: convected_helmholtz_3d_kernel.hpp:40
NiHu::kernel_traits::quadrature_family_t
kernel_traits_ns::quadrature_family< Derived >::type quadrature_family_t
the far field asymptotic behaviour of the kernel
Definition: kernel.hpp:84
NiHu::convected_helmholtz_3d_DLP_tan_kernel
Definition: convected_helmholtz_3d_kernel.hpp:251
NiHu::convected_helmholtz_3d_DLP_kernel
Definition: convected_helmholtz_3d_kernel.hpp:177
NiHu::_m2
std::integral_constant< int, -2 > _m2
shorthand for constant minus two
Definition: guiggiani_1992.hpp:64
NiHu::kernel_traits::far_field_behaviour_t
kernel_traits_ns::far_field_behaviour< Derived >::type far_field_behaviour_t
the far field asymptotic behaviour of the kernel
Definition: kernel.hpp:82
NiHu::singular_kernel_traits
singular traits class of a kernel
Definition: kernel.hpp:101
NiHu::space
class representing a coordinate space with a scalar and a dimension
Definition: space.hpp:36
NiHu::_2
std::integral_constant< int, 2 > _2
shorthand for constant two
Definition: guiggiani_1992.hpp:75
NiHu::kernel_traits::test_input_t
kernel_traits_ns::test_input< Derived >::type test_input_t
kernel test input type
Definition: kernel.hpp:76
NiHu::location_input::x_t
space_t::location_t x_t
the location type
Definition: location_normal.hpp:43
NiHu::location_normal_jacobian_input
a class representing a normal + Jacobian input
Definition: location_normal.hpp:79
NiHu::kernel_traits::is_singular
@ is_singular
indicates if the kernel is singular
Definition: kernel.hpp:91
guiggiani_1992.hpp
Guiggiani's method for CPV and HPF collocational integrals.
NiHu::kernel_base
CRTP base class of all BEM kernels.
Definition: kernel.hpp:133
NiHu::_m1
std::integral_constant< int, -1 > _m1
shorthand for constant minus one
Definition: guiggiani_1992.hpp:69
NiHu::kernel_traits
traits class of a kernel
Definition: kernel.hpp:71
NiHu::singular_kernel_traits::singular_quadrature_order
@ singular_quadrature_order
the quadrature order singular integrals are handled with
Definition: kernel.hpp:111
NiHu::convected_helmholtz_3d_SLP_kernel
SLP kernel of the convected Helmholtz equation in 3D.
Definition: convected_helmholtz_3d_kernel.hpp:109
NiHu::singular_kernel_traits::singularity_type_t
kernel_traits_ns::singularity_type< Derived >::type singularity_type_t
the kernel's singularity type
Definition: kernel.hpp:104
NiHu::kernel_traits::result_t
kernel_traits_ns::result< Derived >::type result_t
the kernel result type
Definition: kernel.hpp:80
NiHu::kernel_traits::is_symmetric
@ is_symmetric
indicates if the kernel is symmetric
Definition: kernel.hpp:89
NiHu::_1
std::integral_constant< int, 1 > _1
shorthand for constant one
Definition: guiggiani_1992.hpp:73
NiHu::asymptotic::inverse< 1 >
NiHu::convected_helmholtz_3d_DLPt_kernel
Definition: convected_helmholtz_3d_kernel.hpp:319
location_normal.hpp
implementation of location and normal kernel inputs
NiHu::location_normal_jacobian_input::x_t
space_t::location_t x_t
the location type
Definition: location_normal.hpp:86
NiHu::polar_laurent_coeffs
definition of Laurent coefficients of singularities
Definition: guiggiani_1992.hpp:64
NiHu::_0
std::integral_constant< int, 0 > _0
shorthand for constant zero
Definition: guiggiani_1992.hpp:71
NiHu::wave_number_kernel
Definition: wave_number_kernel.hpp:8
NiHu::convected_helmholtz_3d_HSP_kernel
Definition: convected_helmholtz_3d_kernel.hpp:478
NiHu::kernel_traits::trial_input_t
kernel_traits_ns::trial_input< Derived >::type trial_input_t
kernel trial input type
Definition: kernel.hpp:78
NiHu::convected_helmholtz_3d_DLPt_n_kernel
Definition: convected_helmholtz_3d_kernel.hpp:396
NiHu::guiggiani
Implementation of Guiggiani's method.
Definition: guiggiani_1992.hpp:84
NiHu::location_input
a class representing a simple location
Definition: location_normal.hpp:37
NiHu::gauss_family_tag
tag for the family of Gaussian quadratures
Definition: gaussian_quadrature.hpp:517