GCC Code Coverage Report

Directory: Bembel/src/
File: Bembel/src/Maxwell/SingleLayerOperator.hpp
Date: 2024-03-19 14:38:05
Exec Total Coverage
Lines: 64 64 100.0%
Functions: 6 6 100.0%
Branches: 67 134 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_MAXWELL_SINGLELAYEROPERATOR_HPP_
13 #define BEMBEL_SRC_MAXWELL_SINGLELAYEROPERATOR_HPP_
14
15 namespace Bembel {
16 // forward declaration of class MaxwellSingleLayerOperator in order to define
17 // traits
18 class MaxwellSingleLayerOperator;
19 /**
20 * \brief Specification of the LinerOperatorTraits for Maxwell.
21 */
22 template <>
23 struct LinearOperatorTraits<MaxwellSingleLayerOperator> {
24 typedef Eigen::VectorXcd EigenType;
25 typedef Eigen::VectorXcd::Scalar Scalar;
26 enum {
27 OperatorOrder = -1,
28 Form = DifferentialForm::DivConforming,
29 NumberOfFMMComponents = 2
30 };
31 };
32
33 /**
34 * \ingroup Maxwell
35 * \brief This class implements the specification of the integration for the
36 * Electric Field Integral Equation.
37 */
38 class MaxwellSingleLayerOperator
39 : public LinearOperatorBase<MaxwellSingleLayerOperator> {
40 // implementation of the kernel evaluation, which may be based on the
41 // information available from the superSpace
42 public:
43 2 MaxwellSingleLayerOperator() {}
44 template <class T>
45 2563584 void evaluateIntegrand_impl(
46 const T &super_space, const SurfacePoint &p1, const SurfacePoint &p2,
47 Eigen::Matrix<
48 typename LinearOperatorTraits<MaxwellSingleLayerOperator>::Scalar,
49 Eigen::Dynamic, Eigen::Dynamic> *intval) const {
50 2563584 auto polynomial_degree = super_space.get_polynomial_degree();
51 2563584 auto polynomial_degree_plus_one_squared =
52 2563584 (polynomial_degree + 1) * (polynomial_degree + 1);
53
54 // get evaluation points on unit square
55
1/2
✓ Branch 1 taken 2563584 times.
✗ Branch 2 not taken.
 2563584 auto s = p1.segment<2>(0);
56
1/2
✓ Branch 1 taken 2563584 times.
✗ Branch 2 not taken.
 2563584 auto t = p2.segment<2>(0);
57
58 // get quadrature weights
59
1/2
✓ Branch 1 taken 2563584 times.
✗ Branch 2 not taken.
 2563584 auto ws = p1(2);
60
1/2
✓ Branch 1 taken 2563584 times.
✗ Branch 2 not taken.
 2563584 auto wt = p2(2);
61
62 // get points on geometry and tangential derivatives
63
1/2
✓ Branch 1 taken 2563584 times.
✗ Branch 2 not taken.
 2563584 auto x_f = p1.segment<3>(3);
64
1/2
✓ Branch 1 taken 2563584 times.
✗ Branch 2 not taken.
 2563584 auto x_f_dx = p1.segment<3>(6);
65
1/2
✓ Branch 1 taken 2563584 times.
✗ Branch 2 not taken.
 2563584 auto x_f_dy = p1.segment<3>(9);
66
1/2
✓ Branch 1 taken 2563584 times.
✗ Branch 2 not taken.
 2563584 auto y_f = p2.segment<3>(3);
67
1/2
✓ Branch 1 taken 2563584 times.
✗ Branch 2 not taken.
 2563584 auto y_f_dx = p2.segment<3>(6);
68
1/2
✓ Branch 1 taken 2563584 times.
✗ Branch 2 not taken.
 2563584 auto y_f_dy = p2.segment<3>(9);
69
70 // compute surface measures from tangential derivatives
71 2563584 auto h = 1. / (1 << super_space.get_refinement_level()); // h = 1 ./ (2^M)
72
73 // integrand without basis functions
74
3/6
✓ Branch 1 taken 2563584 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2563584 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2563584 times.
✗ Branch 8 not taken.
 2563584 auto kernel_evaluation = evaluateKernel(x_f, y_f) * ws * wt;
75 2563584 auto integrand_vector = kernel_evaluation;
76
1/2
✓ Branch 2 taken 2563584 times.
✗ Branch 3 not taken.
 2563584 auto integrand_divergence = -kernel_evaluation / wavenumber2_ / h / h;
77
78 // vector part: mulitply shape functions with integrand and add to buffer
79
7/14
✓ Branch 1 taken 2563584 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2563584 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2563584 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2563584 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2563584 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 2563584 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 2563584 times.
✗ Branch 20 not taken.
 2563584 super_space.addScaledVectorBasisInteraction(intval, integrand_vector, s, t,
80 x_f_dx, x_f_dy, y_f_dx, y_f_dy);
81
82 // divergence part: multiply shape functions with integrand and add to
83 // buffer
84
3/6
✓ Branch 1 taken 2563584 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2563584 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2563584 times.
✗ Branch 8 not taken.
 2563584 super_space.addScaledVectorBasisDivergenceInteraction(
85 intval, integrand_divergence, s, t);
86
87 5127168 return;
88 }
89
90 1259712 Eigen::Matrix<std::complex<double>, 4, 4> evaluateFMMInterpolation_impl(
91 const SurfacePoint &p1, const SurfacePoint &p2) const {
92 // get evaluation points on unit square
93
1/2
✓ Branch 1 taken 1259712 times.
✗ Branch 2 not taken.
 1259712 auto s = p1.segment<2>(0);
94
1/2
✓ Branch 1 taken 1259712 times.
✗ Branch 2 not taken.
 1259712 auto t = p2.segment<2>(0);
95
96 // get points on geometry and tangential derivatives
97
1/2
✓ Branch 1 taken 1259712 times.
✗ Branch 2 not taken.
 1259712 auto x_f = p1.segment<3>(3);
98
1/2
✓ Branch 1 taken 1259712 times.
✗ Branch 2 not taken.
 1259712 auto x_f_dx = p1.segment<3>(6);
99
1/2
✓ Branch 1 taken 1259712 times.
✗ Branch 2 not taken.
 1259712 auto x_f_dy = p1.segment<3>(9);
100
1/2
✓ Branch 1 taken 1259712 times.
✗ Branch 2 not taken.
 1259712 auto y_f = p2.segment<3>(3);
101
1/2
✓ Branch 1 taken 1259712 times.
✗ Branch 2 not taken.
 1259712 auto y_f_dx = p2.segment<3>(6);
102
1/2
✓ Branch 1 taken 1259712 times.
✗ Branch 2 not taken.
 1259712 auto y_f_dy = p2.segment<3>(9);
103
104 // evaluate kernel
105
3/6
✓ Branch 1 taken 1259712 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1259712 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1259712 times.
✗ Branch 8 not taken.
 1259712 auto kernel = evaluateKernel(x_f, y_f);
106
107 // interpolation
108
1/2
✓ Branch 1 taken 1259712 times.
✗ Branch 2 not taken.
 1259712 Eigen::Matrix<std::complex<double>, 4, 4> intval;
109
1/2
✓ Branch 1 taken 1259712 times.
✗ Branch 2 not taken.
 1259712 intval.setZero();
110
2/4
✓ Branch 1 taken 1259712 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1259712 times.
✗ Branch 5 not taken.
 1259712 intval(0, 0) = kernel * x_f_dx.dot(y_f_dx);
111
2/4
✓ Branch 1 taken 1259712 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1259712 times.
✗ Branch 5 not taken.
 1259712 intval(0, 2) = kernel * x_f_dx.dot(y_f_dy);
112
2/4
✓ Branch 1 taken 1259712 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1259712 times.
✗ Branch 5 not taken.
 1259712 intval(2, 0) = kernel * x_f_dy.dot(y_f_dx);
113
2/4
✓ Branch 1 taken 1259712 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1259712 times.
✗ Branch 5 not taken.
 1259712 intval(2, 2) = kernel * x_f_dy.dot(y_f_dy);
114
2/4
✓ Branch 2 taken 1259712 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1259712 times.
✗ Branch 6 not taken.
 1259712 intval(1, 1) = -kernel / wavenumber2_;
115
2/4
✓ Branch 2 taken 1259712 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1259712 times.
✗ Branch 6 not taken.
 1259712 intval(1, 3) = -kernel / wavenumber2_;
116
2/4
✓ Branch 2 taken 1259712 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1259712 times.
✗ Branch 6 not taken.
 1259712 intval(3, 1) = -kernel / wavenumber2_;
117
2/4
✓ Branch 2 taken 1259712 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1259712 times.
✗ Branch 6 not taken.
 1259712 intval(3, 3) = -kernel / wavenumber2_;
118
119 2519424 return intval;
120 }
121
122 /**
123 * \brief Fundamental solution of Helmholtz/Maxwell problem
124 */
125 3823296 std::complex<double> evaluateKernel(const Eigen::Vector3d &x,
126 const Eigen::Vector3d &y) const {
127
2/4
✓ Branch 1 taken 3823296 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3823296 times.
✗ Branch 5 not taken.
 3823296 auto r = (x - y).norm();
128
1/2
✓ Branch 2 taken 3823296 times.
✗ Branch 3 not taken.
 3823296 return std::exp(-std::complex<double>(0., 1.) * wavenumber_ * r) / 4. /
129 7646592 BEMBEL_PI / r;
130 }
131 //////////////////////////////////////////////////////////////////////////////
132 // setters
133 //////////////////////////////////////////////////////////////////////////////
134 2 void set_wavenumber(std::complex<double> wavenumber) {
135 2 wavenumber_ = wavenumber;
136 2 wavenumber2_ = wavenumber_ * wavenumber_;
137 2 }
138 //////////////////////////////////////////////////////////////////////////////
139 // getters
140 //////////////////////////////////////////////////////////////////////////////
141 std::complex<double> get_wavenumber() { return wavenumber_; }
142
143 private:
144 std::complex<double> wavenumber_;
145 std::complex<double> wavenumber2_;
146 };
147
148 /**
149 * \brief The Maxwell single layer operator requires a special treatment of the
150 * moment matrices of the FMM due to the involved derivatives on the ansatz
151 * functions.
152 */
153 template <typename InterpolationPoints>
154 struct H2Multipole::Moment2D<InterpolationPoints, MaxwellSingleLayerOperator> {
155 1 static std::vector<Eigen::MatrixXd> compute2DMoment(
156 const SuperSpace<MaxwellSingleLayerOperator> &super_space,
157 const int cluster_level, const int cluster_refinements,
158 const int number_of_points) {
159
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 Eigen::MatrixXd moment = moment2DComputer<
160 Moment1D<InterpolationPoints, MaxwellSingleLayerOperator>,
161 Moment1D<InterpolationPoints, MaxwellSingleLayerOperator>>(
162 super_space, cluster_level, cluster_refinements, number_of_points);
163
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 Eigen::MatrixXd moment_dx = moment2DComputer<
164 Moment1DDerivative<InterpolationPoints, MaxwellSingleLayerOperator>,
165 Moment1D<InterpolationPoints, MaxwellSingleLayerOperator>>(
166 super_space, cluster_level, cluster_refinements, number_of_points);
167
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 Eigen::MatrixXd moment_dy = moment2DComputer<
168 Moment1D<InterpolationPoints, MaxwellSingleLayerOperator>,
169 Moment1DDerivative<InterpolationPoints, MaxwellSingleLayerOperator>>(
170 super_space, cluster_level, cluster_refinements, number_of_points);
171
172
1/2
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 Eigen::MatrixXd moment1(moment.rows() + moment_dx.rows(), moment.cols());
173
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 moment1 << moment, moment_dx;
174
1/2
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 Eigen::MatrixXd moment2(moment.rows() + moment_dy.rows(), moment.cols());
175
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 moment2 << moment, moment_dy;
176
177 1 std::vector<Eigen::MatrixXd> vector_of_moments;
178
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 vector_of_moments.push_back(moment1);
179
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 vector_of_moments.push_back(moment2);
180
181 2 return vector_of_moments;
182 1 }
183 };
184
185 } // namespace Bembel
186 #endif // BEMBEL_SRC_MAXWELL_SINGLELAYEROPERATOR_HPP_
187