GCC Code Coverage Report

Directory: Bembel/src/
File: Bembel/src/Helmholtz/SingleLayerOperator.hpp
Date: 2024-03-19 14:38:05
Exec Total Coverage
Lines: 41 41 100.0%
Functions: 5 5 100.0%
Branches: 40 80 50.0%
Line Branch Exec Source
1 // This file is part of Bembel, the higher order C++ boundary element library.
2 //
3 // Copyright (C) 2022 see <http://www.bembel.eu>
4 //
5 // It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz,
6 // M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt,
7 // Universitaet Basel, and Universita della Svizzera italiana, Lugano. This
8 // source code is subject to the GNU General Public License version 3 and
9 // provided WITHOUT ANY WARRANTY, see <http://www.bembel.eu> for further
10 // information.
11 //
12 #ifndef BEMBEL_SRC_HELMHOLTZ_SINGLELAYEROPERATOR_HPP_
13 #define BEMBEL_SRC_HELMHOLTZ_SINGLELAYEROPERATOR_HPP_
14
15 namespace Bembel {
16 // forward declaration of class HelmholtzSingleLayerOperator in order to define
17 // traits
18 class HelmholtzSingleLayerOperator;
19 /**
20 * \brief Specification of the LinerOperatorTraits for Helmholtz.
21 */
22 template <>
23 struct LinearOperatorTraits<HelmholtzSingleLayerOperator> {
24 typedef Eigen::VectorXcd EigenType;
25 typedef Eigen::VectorXcd::Scalar Scalar;
26 enum {
27 OperatorOrder = -1,
28 Form = DifferentialForm::Discontinuous,
29 NumberOfFMMComponents = 1
30 };
31 };
32
33 /**
34 * \ingroup Helmholtz
35 * \brief This class implements the specification of the integration for the
36 * single layer potential for Helmholtz.
37 */
38 class HelmholtzSingleLayerOperator
39 : public LinearOperatorBase<HelmholtzSingleLayerOperator> {
40 // implementation of the kernel evaluation, which may be based on the
41 // information available from the superSpace
42 public:
43 4 HelmholtzSingleLayerOperator() {}
44 template <class T>
45 6607321 void evaluateIntegrand_impl(
46 const T &super_space, const SurfacePoint &p1, const SurfacePoint &p2,
47 Eigen::Matrix<
48 typename LinearOperatorTraits<HelmholtzSingleLayerOperator>::Scalar,
49 Eigen::Dynamic, Eigen::Dynamic> *intval) const {
50 6607321 auto polynomial_degree = super_space.get_polynomial_degree();
51 6607321 auto polynomial_degree_plus_one_squared =
52 6607321 (polynomial_degree + 1) * (polynomial_degree + 1);
53
54 // get evaluation points on unit square
55
1/2
✓ Branch 1 taken 6607321 times.
✗ Branch 2 not taken.
 6607321 auto s = p1.segment<2>(0);
56
1/2
✓ Branch 1 taken 6607321 times.
✗ Branch 2 not taken.
 6607321 auto t = p2.segment<2>(0);
57
58 // get quadrature weights
59
1/2
✓ Branch 1 taken 6607321 times.
✗ Branch 2 not taken.
 6607321 auto ws = p1(2);
60
1/2
✓ Branch 1 taken 6607321 times.
✗ Branch 2 not taken.
 6607321 auto wt = p2(2);
61
62 // get points on geometry and tangential derivatives
63
1/2
✓ Branch 1 taken 6607321 times.
✗ Branch 2 not taken.
 6607321 auto x_f = p1.segment<3>(3);
64
1/2
✓ Branch 1 taken 6607321 times.
✗ Branch 2 not taken.
 6607321 auto x_f_dx = p1.segment<3>(6);
65
1/2
✓ Branch 1 taken 6607321 times.
✗ Branch 2 not taken.
 6607321 auto x_f_dy = p1.segment<3>(9);
66
1/2
✓ Branch 1 taken 6607321 times.
✗ Branch 2 not taken.
 6607321 auto y_f = p2.segment<3>(3);
67
1/2
✓ Branch 1 taken 6607321 times.
✗ Branch 2 not taken.
 6607321 auto y_f_dx = p2.segment<3>(6);
68
1/2
✓ Branch 1 taken 6607321 times.
✗ Branch 2 not taken.
 6607321 auto y_f_dy = p2.segment<3>(9);
69
70 // compute surface measures from tangential derivatives
71
2/4
✓ Branch 1 taken 6607321 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6607321 times.
✗ Branch 5 not taken.
 6607321 auto x_kappa = x_f_dx.cross(x_f_dy).norm();
72
2/4
✓ Branch 1 taken 6607321 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6607321 times.
✗ Branch 5 not taken.
 6607321 auto y_kappa = y_f_dx.cross(y_f_dy).norm();
73
74 // integrand without basis functions
75
3/6
✓ Branch 1 taken 6607321 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6607321 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6607321 times.
✗ Branch 8 not taken.
 6607321 auto integrand = evaluateKernel(x_f, y_f) * x_kappa * y_kappa * ws * wt;
76
77 // multiply basis functions with integrand and add to intval, this is an
78 // efficient implementation of
79 // (*intval) += super_space.BasisInteraction(s, t) * evaluateKernel(x_f,
80 // y_f)
81 // * x_kappa * y_kappa * ws * wt;
82
3/6
✓ Branch 1 taken 6607321 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6607321 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6607321 times.
✗ Branch 8 not taken.
 6607321 super_space.addScaledBasisInteraction(intval, integrand, s, t);
83
84 13214642 return;
85 }
86
87 2519424 Eigen::Matrix<std::complex<double>, 1, 1> evaluateFMMInterpolation_impl(
88 const SurfacePoint &p1, const SurfacePoint &p2) const {
89 // get evaluation points on unit square
90
1/2
✓ Branch 1 taken 2519424 times.
✗ Branch 2 not taken.
 2519424 auto s = p1.segment<2>(0);
91
1/2
✓ Branch 1 taken 2519424 times.
✗ Branch 2 not taken.
 2519424 auto t = p2.segment<2>(0);
92
93 // get points on geometry and tangential derivatives
94
1/2
✓ Branch 1 taken 2519424 times.
✗ Branch 2 not taken.
 2519424 auto x_f = p1.segment<3>(3);
95
1/2
✓ Branch 1 taken 2519424 times.
✗ Branch 2 not taken.
 2519424 auto x_f_dx = p1.segment<3>(6);
96
1/2
✓ Branch 1 taken 2519424 times.
✗ Branch 2 not taken.
 2519424 auto x_f_dy = p1.segment<3>(9);
97
1/2
✓ Branch 1 taken 2519424 times.
✗ Branch 2 not taken.
 2519424 auto y_f = p2.segment<3>(3);
98
1/2
✓ Branch 1 taken 2519424 times.
✗ Branch 2 not taken.
 2519424 auto y_f_dx = p2.segment<3>(6);
99
1/2
✓ Branch 1 taken 2519424 times.
✗ Branch 2 not taken.
 2519424 auto y_f_dy = p2.segment<3>(9);
100
101 // compute surface measures from tangential derivatives
102
2/4
✓ Branch 1 taken 2519424 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2519424 times.
✗ Branch 5 not taken.
 2519424 auto x_kappa = x_f_dx.cross(x_f_dy).norm();
103
2/4
✓ Branch 1 taken 2519424 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2519424 times.
✗ Branch 5 not taken.
 2519424 auto y_kappa = y_f_dx.cross(y_f_dy).norm();
104
105 // interpolation
106
1/2
✓ Branch 1 taken 2519424 times.
✗ Branch 2 not taken.
 2519424 Eigen::Matrix<std::complex<double>, 1, 1> intval;
107
4/8
✓ Branch 1 taken 2519424 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2519424 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2519424 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 2519424 times.
✗ Branch 12 not taken.
 2519424 intval(0) = evaluateKernel(x_f, y_f) * x_kappa * y_kappa;
108
109 5038848 return intval;
110 }
111
112 /**
113 * \brief Fundamental solution of Helmholtz problem
114 */
115 9126745 std::complex<double> evaluateKernel(const Eigen::Vector3d &x,
116 const Eigen::Vector3d &y) const {
117
2/4
✓ Branch 1 taken 9126745 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9126745 times.
✗ Branch 5 not taken.
 9126745 auto r = (x - y).norm();
118
1/2
✓ Branch 2 taken 9126745 times.
✗ Branch 3 not taken.
 9126745 return std::exp(-std::complex<double>(0., 1.) * wavenumber_ * r) / 4. /
119 18253490 BEMBEL_PI / r;
120 }
121 //////////////////////////////////////////////////////////////////////////////
122 // setters
123 //////////////////////////////////////////////////////////////////////////////
124 2 void set_wavenumber(std::complex<double> wavenumber) {
125 2 wavenumber_ = wavenumber;
126 2 }
127 //////////////////////////////////////////////////////////////////////////////
128 // getters
129 //////////////////////////////////////////////////////////////////////////////
130 std::complex<double> get_wavenumber() { return wavenumber_; }
131
132 private:
133 std::complex<double> wavenumber_;
134 };
135
136 } // namespace Bembel
137 #endif // BEMBEL_SRC_HELMHOLTZ_SINGLELAYEROPERATOR_HPP_
138