GCC Code Coverage Report

Directory: Bembel/src/
File: H2Matrix/H2Multipole.hpp
Date: 2024-12-18 07:36:36
Exec Total Coverage
Lines: 237 242 97.9%
Functions: 31 31 100.0%
Branches: 258 440 58.6%
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 #ifndef BEMBEL_SRC_H2MATRIX_H2MULTIPOLE_HPP_
12 #define BEMBEL_SRC_H2MATRIX_H2MULTIPOLE_HPP_
13
14 namespace Bembel {
15 namespace H2Multipole {
16 /**
17 * \ingroup H2Matrix
18 * \brief Computes the number_of_points Chebychev points.
19 */
20 struct ChebychevRoots {
21 ChebychevRoots() {}
22 14873 explicit ChebychevRoots(int number_of_points) {
23
1/2
✓ Branch 1 taken 14873 times.
✗ Branch 2 not taken.
 14873 init_ChebyshevRoots(number_of_points);
24 14873 }
25 14873 void init_ChebyshevRoots(int n_pts) {
26
2/4
✓ Branch 1 taken 14873 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14873 times.
✗ Branch 5 not taken.
 14873 auto grid_pts = Eigen::ArrayXd::LinSpaced(n_pts, 1, n_pts).reverse();
27 14873 double alpha = BEMBEL_PI / (2. * double(n_pts));
28
7/14
✓ Branch 1 taken 14873 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14873 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 14873 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 14873 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 14873 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 14873 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 14873 times.
✗ Branch 20 not taken.
 14873 points_ = 0.5 * (alpha * (2 * grid_pts - 1)).cos() + 0.5;
29 29746 return;
30 }
31 Eigen::VectorXd points_;
32 };
33 /**
34 * \ingroup H2Matrix
35 * \brief computes the Lagrange polynomials wrt the interpolation
36 * points given by the InterpolationPoints struct wrt.
37 * Newton basis
38 **/
39 template <typename InterpolationPoints>
40 36 Eigen::MatrixXd computeLagrangePolynomials(int number_of_points) {
41
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
 36 Eigen::MatrixXd retval =
42 Eigen::MatrixXd::Identity(number_of_points, number_of_points);
43
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 Eigen::VectorXd x = InterpolationPoints(number_of_points).points_;
44
2/2
✓ Branch 0 taken 324 times.
✓ Branch 1 taken 36 times.
 360 for (auto i = 0; i < number_of_points; ++i)
45
2/2
✓ Branch 0 taken 2592 times.
✓ Branch 1 taken 324 times.
 2916 for (auto j = 1; j < number_of_points; ++j)
46
2/2
✓ Branch 0 taken 11664 times.
✓ Branch 1 taken 2592 times.
 14256 for (auto k = number_of_points - 1; k >= j; --k)
47
5/10
✓ Branch 1 taken 11664 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 11664 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 11664 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 11664 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 11664 times.
✗ Branch 14 not taken.
 11664 retval(k, i) = (retval(k, i) - retval(k - 1, i)) / (x(k) - x(k - j));
48 72 return retval;
49 36 }
50 /**
51 * \ingroup H2Matrix
52 * \brief evaluates a given polynomial in the Newton basis wrt the
53 * interpolation points at a given location
54 **/
55 template <typename InterpolationPoints>
56 14796 double evaluatePolynomial(const Eigen::VectorXd &L, double xi) {
57 14796 int number_of_points = L.size();
58
1/2
✓ Branch 1 taken 14796 times.
✗ Branch 2 not taken.
 14796 Eigen::VectorXd x = InterpolationPoints(number_of_points).points_;
59 14796 double retval = 0;
60
61
1/2
✓ Branch 1 taken 14796 times.
✗ Branch 2 not taken.
 14796 retval = L(number_of_points - 1);
62
63
2/2
✓ Branch 0 taken 118368 times.
✓ Branch 1 taken 14796 times.
 133164 for (auto i = number_of_points - 2; i >= 0; --i)
64
2/4
✓ Branch 1 taken 118368 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 118368 times.
✗ Branch 5 not taken.
 118368 retval = retval * (xi - x(i)) + L(i);
65
66 14796 return retval;
67 14796 }
68 /**
69 * \ingroup H2Matrix
70 * \brief Computes transfer matrices in required order to apply them
71 * all-at-once in a matrix-vector-product, i.e. the order
72 * in the output is [T0 T2 T3 T1].
73 *
74 *
75 */
76 template <typename InterpolationPoints>
77 8 Eigen::MatrixXd computeTransferMatrices(int number_of_points) {
78 8 int np2 = number_of_points * number_of_points;
79
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 Eigen::MatrixXd E(number_of_points, 2 * number_of_points);
80
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 Eigen::MatrixXd T(np2, 4 * np2);
81 /// initialize Lagrange polynomials and interpolation nodes
82
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 Eigen::VectorXd x = InterpolationPoints(number_of_points).points_;
83
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 Eigen::MatrixXd L =
84 computeLagrangePolynomials<InterpolationPoints>(number_of_points);
85 /**
86 * initialize values of Lagrange functions on the interpolation points
87 * as follows:
88 * E(:,0:npts-1) containes the values of the Lagrange polynomials on the
89 * interpolation points scaled to [0, 0.5].
90 * E(:,npts:2*npts-1) containes the values of the Lagrange polynomials on
91 * the interpolation points scaled to [0.5, 1].
92 * E(i,j) containes the value of the j.th polynomial on the ith point
93 **/
94
2/2
✓ Branch 0 taken 72 times.
✓ Branch 1 taken 8 times.
 80 for (auto j = 0; j < number_of_points; ++j)
95
2/2
✓ Branch 0 taken 648 times.
✓ Branch 1 taken 72 times.
 720 for (auto i = 0; i < number_of_points; ++i) {
96
5/10
✓ Branch 1 taken 648 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 648 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 648 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 648 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 648 times.
✗ Branch 14 not taken.
 648 E(i, j) = evaluatePolynomial<InterpolationPoints>(L.col(j), 0.5 * x(i));
97
1/2
✓ Branch 1 taken 648 times.
✗ Branch 2 not taken.
 648 E(i, j + number_of_points) =
98
4/8
✓ Branch 1 taken 648 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 648 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 648 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 648 times.
✗ Branch 11 not taken.
 1296 evaluatePolynomial<InterpolationPoints>(L.col(j), 0.5 * x(i) + 0.5);
99 }
100 // This construction, regard permuation vector, results in the
101 // order T0 T2 T3 T1 to apply to the moments
102
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 Eigen::Vector4i permutation;
103
4/8
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 8 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 8 times.
✗ Branch 11 not taken.
 8 permutation << 0, 3, 1, 2;
104
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 8 times.
 40 for (auto k = 0; k < 4; ++k)
105
2/2
✓ Branch 0 taken 288 times.
✓ Branch 1 taken 32 times.
 320 for (auto i = 0; i < number_of_points; ++i)
106
2/2
✓ Branch 0 taken 2592 times.
✓ Branch 1 taken 288 times.
 2880 for (auto ii = 0; ii < number_of_points; ++ii)
107
2/2
✓ Branch 0 taken 23328 times.
✓ Branch 1 taken 2592 times.
 25920 for (auto j = 0; j < number_of_points; ++j)
108
2/2
✓ Branch 0 taken 209952 times.
✓ Branch 1 taken 23328 times.
 233280 for (auto jj = 0; jj < number_of_points; ++jj)
109 419904 T(j * number_of_points + jj,
110
1/2
✓ Branch 1 taken 209952 times.
✗ Branch 2 not taken.
 209952 i * number_of_points + ii + np2 * permutation(k)) =
111
1/2
✓ Branch 1 taken 209952 times.
✗ Branch 2 not taken.
 419904 E(i, j + (k / 2) * number_of_points) *
112
2/4
✓ Branch 1 taken 209952 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 209952 times.
✗ Branch 5 not taken.
 209952 E(ii, jj + (k % 2) * number_of_points);
113
114 16 return T;
115 8 }
116 /**
117 * \ingroup H2Matrix
118 * \brief Computes 1D moment for FMM. All calculations ar performed on [0,1]
119 */
120 template <typename InterpolationPoints, typename LinOp>
121 struct Moment1D {
122 26 static Eigen::MatrixXd computeMoment1D(const SuperSpace<LinOp> &super_space,
123 const int cluster_level,
124 const int cluster_refinements,
125 const int number_of_points) {
126 26 int n = 1 << cluster_refinements;
127 26 double h = 1. / n;
128 26 int N = 1 << cluster_level;
129 26 double H = 1. / N;
130 26 int polynomial_degree = super_space.get_polynomial_degree();
131 26 int polynomial_degree_plus_one = polynomial_degree + 1;
132 26 int polynomial_degree_plus_one_squared =
133 26 polynomial_degree_plus_one * polynomial_degree_plus_one;
134
1/2
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
 26 GaussLegendre<Constants::maximum_quadrature_degree> GL;
135
1/2
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
 26 auto Q =
136
1/2
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
 26 GL[(int)std::ceil(0.5 * (number_of_points + polynomial_degree - 2))];
137
138
1/2
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
 26 Eigen::VectorXd x = InterpolationPoints(number_of_points).points_;
139
1/2
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
 26 Eigen::MatrixXd L =
140 computeLagrangePolynomials<InterpolationPoints>(number_of_points);
141
142
1/2
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
 26 Eigen::MatrixXd moment(number_of_points, n * polynomial_degree_plus_one);
143
144
2/2
✓ Branch 1 taken 234 times.
✓ Branch 2 taken 26 times.
 260 for (auto i = 0; i < number_of_points; ++i) {
145 Eigen::Matrix<typename LinearOperatorTraits<LinOp>::Scalar,
146 Eigen::Dynamic, 1>
147
1/2
✓ Branch 1 taken 234 times.
✗ Branch 2 not taken.
 234 intval(polynomial_degree_plus_one, 1);
148
2/2
✓ Branch 0 taken 2664 times.
✓ Branch 1 taken 234 times.
 2898 for (auto j = 0; j < n; ++j) {
149
1/2
✓ Branch 1 taken 2664 times.
✗ Branch 2 not taken.
 2664 intval.setZero();
150
2/2
✓ Branch 1 taken 13320 times.
✓ Branch 2 taken 2664 times.
 15984 for (auto k = 0; k < Q.w_.size(); ++k) {
151
1/2
✓ Branch 1 taken 720 times.
✗ Branch 2 not taken.
 25920 super_space.addScaledBasis1D(
152 &intval,
153 26640 Q.w_(k) * std::sqrt(h * H) *
154
3/6
✓ Branch 1 taken 13320 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13320 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 12600 times.
✗ Branch 8 not taken.
 14040 evaluatePolynomial<InterpolationPoints>(L.col(i),
155
2/4
✓ Branch 1 taken 13320 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13320 times.
✗ Branch 5 not taken.
 13320 h * (j + Q.xi_(k))),
156
2/4
✓ Branch 1 taken 13320 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13320 times.
✗ Branch 5 not taken.
 13320 Q.xi_(k));
157 }
158 // we are integrating real stuff, so take the real part
159
1/2
✓ Branch 1 taken 2664 times.
✗ Branch 2 not taken.
 2664 moment.block(i, j * polynomial_degree_plus_one, 1,
160
5/8
✓ Branch 1 taken 144 times.
✓ Branch 2 taken 2520 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 144 times.
✓ Branch 5 taken 2520 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 144 times.
✗ Branch 8 not taken.
 5328 polynomial_degree_plus_one) = intval.real().transpose();
161 }
162 }
163
164 52 return moment;
165 26 }
166 };
167 /**
168 * \ingroup H2Matrix
169 * \brief Computes 1D moment for FMM using derivatives of the basis functions.
170 * All calculations ar performed on [0,1].
171 */
172 template <typename InterpolationPoints, typename LinOp>
173 struct Moment1DDerivative {
174 2 static Eigen::MatrixXd computeMoment1D(const SuperSpace<LinOp> &super_space,
175 const int cluster_level,
176 const int cluster_refinements,
177 const int number_of_points) {
178 2 int n = 1 << cluster_refinements;
179 2 double h = 1. / n;
180 2 int N = 1 << cluster_level;
181 2 double H = 1. / N;
182 2 int polynomial_degree = super_space.get_polynomial_degree();
183 2 int polynomial_degree_plus_one = polynomial_degree + 1;
184 2 int polynomial_degree_plus_one_squared =
185 2 polynomial_degree_plus_one * polynomial_degree_plus_one;
186
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 GaussLegendre<Constants::maximum_quadrature_degree> GL;
187
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 auto Q =
188
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 GL[(int)std::ceil(0.5 * (number_of_points + polynomial_degree - 2))];
189
190
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 Eigen::VectorXd x = InterpolationPoints(number_of_points).points_;
191
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 Eigen::MatrixXd L =
192 computeLagrangePolynomials<InterpolationPoints>(number_of_points);
193
194
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 Eigen::MatrixXd moment(number_of_points, n * polynomial_degree_plus_one);
195
196
2/2
✓ Branch 1 taken 18 times.
✓ Branch 2 taken 2 times.
 20 for (auto i = 0; i < number_of_points; ++i) {
197 Eigen::Matrix<typename LinearOperatorTraits<LinOp>::Scalar,
198 Eigen::Dynamic, 1>
199
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 18 intval(polynomial_degree_plus_one, 1);
200
2/2
✓ Branch 0 taken 36 times.
✓ Branch 1 taken 18 times.
 54 for (auto j = 0; j < n; ++j) {
201
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 intval.setZero();
202
2/2
✓ Branch 1 taken 180 times.
✓ Branch 2 taken 36 times.
 216 for (auto k = 0; k < Q.w_.size(); ++k) {
203
1/2
✓ Branch 1 taken 180 times.
✗ Branch 2 not taken.
 180 super_space.addScaledBasis1DDx(
204 &intval,
205 360 Q.w_(k) / std::sqrt(h * H) *
206
2/4
✓ Branch 1 taken 180 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 180 times.
✗ Branch 5 not taken.
 360 evaluatePolynomial<InterpolationPoints>(L.col(i),
207
2/4
✓ Branch 1 taken 180 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 180 times.
✗ Branch 5 not taken.
 180 h * (j + Q.xi_(k))),
208
2/4
✓ Branch 1 taken 180 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 180 times.
✗ Branch 5 not taken.
 180 Q.xi_(k));
209 }
210 // we are integrating real stuff, so take the real part
211
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 moment.block(i, j * polynomial_degree_plus_one, 1,
212
3/6
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 36 times.
✗ Branch 8 not taken.
 72 polynomial_degree_plus_one) = intval.real().transpose();
213 }
214 }
215
216 4 return moment;
217 2 }
218 };
219 /**
220 * \ingroup H2Matrix
221 * \brief Computes a single 2D moment for the FMM by tensorisation of the 1D
222 * moments
223 */
224 template <typename Mom1D_1, typename Mom1D_2, typename LinOp>
225 17 Eigen::MatrixXd moment2DComputer(const SuperSpace<LinOp> &super_space,
226 const int cluster_level,
227 const int cluster_refinements,
228 const int number_of_points) {
229 17 int n = 1 << cluster_refinements;
230 17 auto n2 = n * n;
231 17 int polynomial_degree = super_space.get_polynomial_degree();
232 17 int polynomial_degree_plus_one = polynomial_degree + 1;
233 17 int polynomial_degree_plus_one_squared =
234 polynomial_degree_plus_one * polynomial_degree_plus_one;
235
236 // compute 1D moments
237
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 17 Eigen::MatrixXd moment1D_1 = Mom1D_1::computeMoment1D(
238 super_space, cluster_level, cluster_refinements, number_of_points);
239
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 17 Eigen::MatrixXd moment1D_2 = Mom1D_2::computeMoment1D(
240 super_space, cluster_level, cluster_refinements, number_of_points);
241 /**
242 * Throughout this code we face the problem of memory serialisation for
243 * the traversel of elements in the element tree. the canonical orders would
244 * either be row major or column major traversal of a given patch
245 * resulting in e.g. element(i,j) = *(elementRoot + i * n + j)
246 * for n = 1 << j. However, to achieve a localisation of matrix blocks if
247 * their corresponding shape functions are geometrically close, we typically
248 * use a hierarchical ordering of the elements. The mapping from hierarchical
249 * to row major ordering is easily facilitated by computing the elements
250 * position in its patch using its llc_. For our special situation of a
251 * balanced quadtree, this mapping is the same for every patch.
252 * In particular, the moments in our construction of the FMM only depend on
253 * the current level of uniform refinement and not on a particular patch.
254 * Next, we explicitly set up this mapping, where index_s determines the
255 * location of a given element along the first coordinate in the reference
256 * domain and index_t its second coordinate in the reference domain.
257 * Note that this piece of code only has to be used once in the entire setup
258 * of the FMM, such that the overhead is quite small.
259 **/
260
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 17 Eigen::VectorXi index_s(n2);
261
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 17 Eigen::VectorXi index_t(n2);
262
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 17 index_s.setZero();
263
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 17 index_t.setZero();
264
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 14 times.
 52 for (auto j = 0; j < cluster_refinements; ++j) {
265 35 auto inc = 1 << 2 * j;
266
2/2
✓ Branch 0 taken 1850 times.
✓ Branch 1 taken 32 times.
 1888 for (auto i = 0; i < inc; ++i) {
267
2/4
✓ Branch 1 taken 1850 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1850 times.
✗ Branch 5 not taken.
 1853 index_s(4 * i) = 2 * index_s(i);
268
2/4
✓ Branch 1 taken 1850 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1850 times.
✗ Branch 5 not taken.
 1853 index_s(4 * i + 1) = 2 * index_s(i) + 1;
269
2/4
✓ Branch 1 taken 1850 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1850 times.
✗ Branch 5 not taken.
 1853 index_s(4 * i + 2) = 2 * index_s(i) + 1;
270
2/4
✓ Branch 1 taken 1850 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1850 times.
✗ Branch 5 not taken.
 1853 index_s(4 * i + 3) = 2 * index_s(i);
271
2/4
✓ Branch 1 taken 1850 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1850 times.
✗ Branch 5 not taken.
 1853 index_t(4 * i) = 2 * index_t(i);
272
2/4
✓ Branch 1 taken 1850 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1850 times.
✗ Branch 5 not taken.
 1853 index_t(4 * i + 1) = 2 * index_t(i);
273
2/4
✓ Branch 1 taken 1850 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1850 times.
✗ Branch 5 not taken.
 1853 index_t(4 * i + 2) = 2 * index_t(i) + 1;
274
2/4
✓ Branch 1 taken 1850 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1850 times.
✗ Branch 5 not taken.
 1853 index_t(4 * i + 3) = 2 * index_t(i) + 1;
275 }
276 }
277
278 // assemble 2D tensor-product moments
279 ✗ Eigen::MatrixXd moment2D(number_of_points * number_of_points,
280
1/2
✓ Branch 3 taken 14 times.
✗ Branch 4 not taken.
 17 moment1D_1.cols() * moment1D_2.cols());
281
2/2
✓ Branch 0 taken 126 times.
✓ Branch 1 taken 14 times.
 170 for (auto i = 0; i < number_of_points; ++i)
282
2/2
✓ Branch 0 taken 1134 times.
✓ Branch 1 taken 126 times.
 1530 for (auto j = 0; j < number_of_points; ++j)
283
2/2
✓ Branch 0 taken 450684 times.
✓ Branch 1 taken 1134 times.
 453033 for (auto k = 0; k < n2; ++k)
284
2/2
✓ Branch 0 taken 451980 times.
✓ Branch 1 taken 450684 times.
 905580 for (auto m1 = 0; m1 < polynomial_degree_plus_one; ++m1)
285
2/2
✓ Branch 0 taken 454572 times.
✓ Branch 1 taken 451980 times.
 912384 for (auto m2 = 0; m2 < polynomial_degree_plus_one; ++m2)
286 916920 moment2D(i * number_of_points + j,
287 458460 polynomial_degree_plus_one_squared * k +
288 916920 m1 * polynomial_degree_plus_one + m2) =
289
2/4
✓ Branch 1 taken 454572 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 454572 times.
✗ Branch 5 not taken.
 916920 moment1D_1(i, index_s(k) * polynomial_degree_plus_one + m2) *
290
3/6
✓ Branch 1 taken 454572 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 454572 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 454572 times.
✗ Branch 8 not taken.
 458460 moment1D_2(j, index_t(k) * polynomial_degree_plus_one + m1);
291
292 34 return moment2D;
293 17 }
294 /**
295 * \ingroup H2Matrix
296 * \brief Computes all 2D moment for the FMM by tensorisation of the 1D
297 * moments. Specialice this for your linear operator if you need derivatives on
298 * your local shape functions. See e.g. MaxwellSingleLayerOperator for an
299 * example.
300 */
301 template <typename InterpolationPoints, typename LinOp>
302 struct Moment2D {
303 11 static std::vector<Eigen::MatrixXd> compute2DMoment(
304 const SuperSpace<LinOp> &super_space, const int cluster_level,
305 const int cluster_refinements, const int number_of_points) {
306 11 std::vector<Eigen::MatrixXd> vector_of_moments;
307 11 for (int i = 0;
308
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 11 times.
 22 i <
309 22 getFunctionSpaceVectorDimension<LinearOperatorTraits<LinOp>::Form>();
310 ++i)
311
2/4
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 11 times.
✗ Branch 5 not taken.
 11 vector_of_moments.push_back(
312 moment2DComputer<Moment1D<InterpolationPoints, LinOp>,
313 Moment1D<InterpolationPoints, LinOp>>(
314 super_space, cluster_level, cluster_refinements,
315 number_of_points));
316 11 return vector_of_moments;
317 ✗ }
318 };
319 /**
320 * \ingroup H2Matrix
321 * \brief Compute tensor interpolation points on unit square from 1D
322 * interpolation points.
323 */
324 5 Eigen::MatrixXd interpolationPoints2D(const Eigen::VectorXd &x) {
325
1/2
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
 5 Eigen::MatrixXd x2(x.size() * x.size(), 2);
326
2/2
✓ Branch 1 taken 45 times.
✓ Branch 2 taken 5 times.
 50 for (auto i = 0; i < x.size(); ++i)
327
2/2
✓ Branch 1 taken 405 times.
✓ Branch 2 taken 45 times.
 450 for (auto j = 0; j < x.size(); ++j) {
328
5/10
✓ Branch 1 taken 405 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 405 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 405 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 405 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 405 times.
✗ Branch 15 not taken.
 405 x2.row(i * x.size() + j) = Eigen::Vector2d(x(i), x(j));
329 }
330 5 return x2;
331 ✗ }
332 /**
333 * \ingroup H2Matrix
334 * \brief Interpolate kernel function on reference domain for FMM.
335 */
336 template <typename LinOp>
337 Eigen::Matrix<typename LinearOperatorTraits<LinOp>::Scalar, Eigen::Dynamic,
338 Eigen::Dynamic>
339 1056 interpolateKernel(const LinOp &linOp, const SuperSpace<LinOp> &super_space,
340 const Eigen::MatrixXd &x,
341 const Bembel::ElementTreeNode &cluster1,
342 const Bembel::ElementTreeNode &cluster2) {
343
2/4
✓ Branch 1 taken 1056 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1056 times.
✗ Branch 5 not taken.
 1056 SurfacePoint qp1, qp2;
344 Eigen::Matrix<
345 typename LinearOperatorTraits<LinOp>::Scalar,
346 getFunctionSpaceVectorDimension<LinearOperatorTraits<LinOp>::Form>() *
347 LinearOperatorTraits<LinOp>::NumberOfFMMComponents,
348 getFunctionSpaceVectorDimension<LinearOperatorTraits<LinOp>::Form>() *
349 LinearOperatorTraits<LinOp>::NumberOfFMMComponents>
350
1/2
✓ Branch 1 taken 1056 times.
✗ Branch 2 not taken.
 1056 interpval;
351 1056 const int vector_dimension = Bembel::getFunctionSpaceVectorDimension<
352 Bembel::LinearOperatorTraits<LinOp>::Form>();
353 1056 const int number_of_FMM_components =
354 Bembel::LinearOperatorTraits<LinOp>::NumberOfFMMComponents;
355 Eigen::Matrix<typename LinearOperatorTraits<LinOp>::Scalar, Eigen::Dynamic,
356 Eigen::Dynamic>
357
1/2
✓ Branch 1 taken 1056 times.
✗ Branch 2 not taken.
 1056 F(vector_dimension * number_of_FMM_components * x.rows(),
358 1056 vector_dimension * number_of_FMM_components * x.rows());
359
2/2
✓ Branch 1 taken 85536 times.
✓ Branch 2 taken 1056 times.
 86592 for (int i = 0; i < x.rows(); ++i) {
360
3/6
✓ Branch 1 taken 85536 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 85536 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 85536 times.
✗ Branch 8 not taken.
 85536 super_space.map2surface(cluster1, x.row(i), 1., &qp1);
361
2/2
✓ Branch 1 taken 6928416 times.
✓ Branch 2 taken 85536 times.
 7013952 for (int j = 0; j < x.rows(); ++j) {
362
3/6
✓ Branch 1 taken 6928416 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6928416 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6928416 times.
✗ Branch 8 not taken.
 6928416 super_space.map2surface(cluster2, x.row(j), 1., &qp2);
363
1/2
✓ Branch 1 taken 6928416 times.
✗ Branch 2 not taken.
 6928416 auto FMM_output = linOp.evaluateFMMInterpolation(qp1, qp2);
364
2/2
✓ Branch 0 taken 8188128 times.
✓ Branch 1 taken 6928416 times.
 15116544 for (int ii = 0; ii < vector_dimension; ++ii)
365
2/2
✓ Branch 0 taken 10707552 times.
✓ Branch 1 taken 8188128 times.
 18895680 for (int jj = 0; jj < vector_dimension; ++jj)
366
2/2
✓ Branch 0 taken 15746400 times.
✓ Branch 1 taken 10707552 times.
 26453952 for (int iii = 0; iii < number_of_FMM_components; ++iii)
367
2/2
✓ Branch 0 taken 25824096 times.
✓ Branch 1 taken 15746400 times.
 41570496 for (int jjj = 0; jjj < number_of_FMM_components; ++jjj)
368
1/2
✓ Branch 1 taken 25824096 times.
✗ Branch 2 not taken.
 25824096 F(i + (ii * number_of_FMM_components + iii) * x.rows(),
369 28973376 j + (jj * number_of_FMM_components + jjj) * x.rows()) =
370 25824096 FMM_output(iii + ii * number_of_FMM_components,
371
1/2
✓ Branch 1 taken 25824096 times.
✗ Branch 2 not taken.
 25824096 jjj + jj * number_of_FMM_components);
372 }
373 }
374 2112 return F;
375 ✗ }
376 /**
377 * \ingroup H2Matrix
378 * \brief Forward transformation for FMM.
379 */
380 template <typename Scalar>
381 std::vector<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>>
382 775 forwardTransformation(const Eigen::MatrixXd &moment_matrices,
383 const Eigen::MatrixXd &transfer_matrices, const int steps,
384 const Eigen::Matrix<Scalar, Eigen::Dynamic,
385 Eigen::Dynamic> &long_rhs_matrix) {
386 // get numbers
387 775 int number_of_points = transfer_matrices.rows();
388 775 int number_of_FMM_components = moment_matrices.rows() / number_of_points;
389 // apply moment matrices
390 std::vector<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>>
391 775 long_rhs_forward;
392
3/6
✓ Branch 1 taken 775 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 775 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 775 times.
✗ Branch 8 not taken.
 775 long_rhs_forward.push_back(moment_matrices * long_rhs_matrix);
393 // apply transfer matrices
394
2/2
✓ Branch 2 taken 778 times.
✓ Branch 3 taken 775 times.
 1553 for (int i = 0; i < steps; ++i) {
395
1/2
✓ Branch 1 taken 778 times.
✗ Branch 2 not taken.
 778 Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> reshaped =
396
3/6
✓ Branch 1 taken 778 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 778 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 778 times.
✗ Branch 8 not taken.
 1556 Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>>(
397 778 long_rhs_forward.back().data(), 4 * long_rhs_forward.back().rows(),
398 778 long_rhs_forward.back().cols() / 4);
399 778 Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> transferred(
400
1/2
✓ Branch 3 taken 778 times.
✗ Branch 4 not taken.
 778 reshaped.rows() / 4, reshaped.cols());
401
402 // iterate over FMM components
403
2/2
✓ Branch 0 taken 1410 times.
✓ Branch 1 taken 778 times.
 2188 for (int j = 0; j < number_of_FMM_components; ++j) {
404 1410 Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> to_transfer(
405
1/2
✓ Branch 2 taken 1410 times.
✗ Branch 3 not taken.
 1410 4 * number_of_points, reshaped.cols());
406 1410 to_transfer << reshaped.block(
407
1/2
✓ Branch 1 taken 1410 times.
✗ Branch 2 not taken.
 1410 (j + 0 * number_of_FMM_components) * number_of_points, 0,
408
2/4
✓ Branch 1 taken 1410 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1410 times.
✗ Branch 5 not taken.
 2820 number_of_points, reshaped.cols()),
409
1/2
✓ Branch 2 taken 1410 times.
✗ Branch 3 not taken.
 1410 reshaped.block((j + 1 * number_of_FMM_components) * number_of_points,
410
1/2
✓ Branch 1 taken 1410 times.
✗ Branch 2 not taken.
 1410 0, number_of_points, reshaped.cols()),
411
1/2
✓ Branch 2 taken 1410 times.
✗ Branch 3 not taken.
 1410 reshaped.block((j + 2 * number_of_FMM_components) * number_of_points,
412
1/2
✓ Branch 1 taken 1410 times.
✗ Branch 2 not taken.
 1410 0, number_of_points, reshaped.cols()),
413
1/2
✓ Branch 2 taken 1410 times.
✗ Branch 3 not taken.
 2820 reshaped.block((j + 3 * number_of_FMM_components) * number_of_points,
414 0, number_of_points, reshaped.cols());
415
1/2
✓ Branch 2 taken 1410 times.
✗ Branch 3 not taken.
 1410 transferred.block(j * number_of_points, 0, number_of_points,
416
2/4
✓ Branch 1 taken 1410 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1410 times.
✗ Branch 5 not taken.
 2820 reshaped.cols()) = transfer_matrices * to_transfer;
417 }
418
419
1/2
✓ Branch 1 taken 778 times.
✗ Branch 2 not taken.
 778 long_rhs_forward.push_back(transferred);
420 }
421 775 return long_rhs_forward;
422 ✗ }
423 /**
424 * \ingroup H2Matrix
425 * \brief Backward transformation for FMM. The content of long_dst_backward is
426 * destroyed.
427 */
428 template <typename Scalar>
429 775 Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> backwardTransformation(
430 const Eigen::MatrixXd &moment_matrices,
431 const Eigen::MatrixXd &transfer_matrices, const int steps,
432 std::vector<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>>
433 &long_dst_backward) {
434 // get numbers
435 775 int number_of_points = transfer_matrices.rows();
436 775 int number_of_FMM_components = moment_matrices.rows() / number_of_points;
437 // apply transfer matrices
438
2/2
✓ Branch 0 taken 775 times.
✓ Branch 1 taken 775 times.
 1550 for (int i = steps; i > 0; --i) {
439 775 Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> transferred(
440
1/2
✓ Branch 5 taken 775 times.
✗ Branch 6 not taken.
 775 4 * long_dst_backward[i].rows(), long_dst_backward[i].cols());
441
2/2
✓ Branch 0 taken 1404 times.
✓ Branch 1 taken 775 times.
 2179 for (int j = 0; j < number_of_FMM_components; ++j) {
442
1/2
✓ Branch 1 taken 1404 times.
✗ Branch 2 not taken.
 1404 Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> prod =
443
2/4
✓ Branch 1 taken 1404 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1404 times.
✗ Branch 5 not taken.
 1404 transfer_matrices.transpose() *
444
1/2
✓ Branch 3 taken 1404 times.
✗ Branch 4 not taken.
 2808 long_dst_backward[i].block(j * number_of_points, 0, number_of_points,
445 1404 long_dst_backward[i].cols());
446
1/2
✓ Branch 2 taken 1404 times.
✗ Branch 3 not taken.
 1404 transferred.block((j + 0 * number_of_FMM_components) * number_of_points,
447
1/2
✓ Branch 1 taken 1404 times.
✗ Branch 2 not taken.
 1404 0, number_of_points, prod.cols()) =
448
1/2
✓ Branch 2 taken 1404 times.
✗ Branch 3 not taken.
 1404 prod.block(0 * number_of_points, 0, number_of_points, prod.cols());
449
1/2
✓ Branch 2 taken 1404 times.
✗ Branch 3 not taken.
 1404 transferred.block((j + 1 * number_of_FMM_components) * number_of_points,
450
1/2
✓ Branch 1 taken 1404 times.
✗ Branch 2 not taken.
 1404 0, number_of_points, prod.cols()) =
451
1/2
✓ Branch 2 taken 1404 times.
✗ Branch 3 not taken.
 1404 prod.block(1 * number_of_points, 0, number_of_points, prod.cols());
452
1/2
✓ Branch 2 taken 1404 times.
✗ Branch 3 not taken.
 1404 transferred.block((j + 2 * number_of_FMM_components) * number_of_points,
453
1/2
✓ Branch 1 taken 1404 times.
✗ Branch 2 not taken.
 1404 0, number_of_points, prod.cols()) =
454
1/2
✓ Branch 2 taken 1404 times.
✗ Branch 3 not taken.
 1404 prod.block(2 * number_of_points, 0, number_of_points, prod.cols());
455
1/2
✓ Branch 2 taken 1404 times.
✗ Branch 3 not taken.
 1404 transferred.block((j + 3 * number_of_FMM_components) * number_of_points,
456
1/2
✓ Branch 1 taken 1404 times.
✗ Branch 2 not taken.
 1404 0, number_of_points, prod.cols()) =
457
1/2
✓ Branch 2 taken 1404 times.
✗ Branch 3 not taken.
 2808 prod.block(3 * number_of_points, 0, number_of_points, prod.cols());
458 }
459
1/2
✓ Branch 2 taken 775 times.
✗ Branch 3 not taken.
 775 long_dst_backward[i - 1] +=
460
3/6
✓ Branch 1 taken 775 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 775 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 775 times.
✗ Branch 10 not taken.
 3100 Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>>(
461 775 transferred.data(), long_dst_backward[i - 1].rows(),
462 775 long_dst_backward[i - 1].cols());
463 }
464 // apply moment matrices
465
1/2
✓ Branch 1 taken 775 times.
✗ Branch 2 not taken.
 775 Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> long_dst_matrix =
466
2/4
✓ Branch 2 taken 775 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 775 times.
✗ Branch 6 not taken.
 775 moment_matrices.transpose() * long_dst_backward[0];
467
3/6
✓ Branch 1 taken 775 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 775 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 775 times.
✗ Branch 9 not taken.
 775 return Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, 1>>(
468
1/2
✓ Branch 1 taken 775 times.
✗ Branch 2 not taken.
 1550 long_dst_matrix.data(), long_dst_matrix.size());
469 775 }
470
471 } // namespace H2Multipole
472 } // namespace Bembel
473 #endif // BEMBEL_SRC_H2MATRIX_H2MULTIPOLE_HPP_
474