GCC Code Coverage Report

Directory: Bembel/src/
File: HomogenisedLaplace/Coefficients.hpp
Date: 2024-12-18 07:36:36
Exec Total Coverage
Lines: 174 176 98.9%
Functions: 8 8 100.0%
Branches: 335 636 52.7%
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_HOMOGENISEDLAPLACE_COEFFICIENTS_HPP_
13 #define BEMBEL_SRC_HOMOGENISEDLAPLACE_COEFFICIENTS_HPP_
14
15 #ifndef POINT_DEGREE
16 #define POINT_DEGREE 20 /** the number of points on the surface of the cube */
17 #endif
18
19 #include <Eigen/Dense>
20 #include <Bembel/HomogenisedLaplace>
21 #include <Bembel/Quadrature>
22
23 #include <functional>
24
25 namespace Bembel {
26
27 Eigen::VectorXd getCoefficients(double precision);
28
29 inline unsigned int getDegree(double precision);
30
31 Eigen::VectorXd getDisplacement(Eigen::MatrixXd ps_l, Eigen::MatrixXd ps_f,
32 Eigen::MatrixXd ps_b);
33
34 inline double k_mod(Eigen::Vector3d in);
35
36 inline Eigen::Vector3d Dk_mod(Eigen::Vector3d in);
37
38 double calculateFirstCoefficient(Eigen::VectorXd cs, unsigned int deg,
39 Eigen::MatrixXd ps_l, Eigen::MatrixXd ps_f, Eigen::MatrixXd ps_b);
40
41 /**
42 * \brief Calculates the coefficients for the solid harmonics expansion of
43 * the periodic kernel \f$ k(x) = \sum_{m \in \{-1, 0, 1\}^3} \frac{1}{4 \pi |x-m|}
44 * + \frac{|x|^2}{6} + \sum_{n = 0}^{N} \sum_{m = -n}^n c_m^n \phi_n^m(x) \f$.
45 *
46 * @param precision: The desired mean error from periodicity.
47 *
48 * \see A. Barnett and L. Greengard. A new integral representation for quasi-periodic
49 * fields and its application to two-dimensional band structure calculations.
50 * _Journal of Computational Physics_, 229(19):6898--6914, 2010.
51 *
52 * \see P. Cazeaux and O. Zahm. A fast boundary element method
53 * for the solution of periodic many-inclusion problems via
54 * hierarchical matrix techniques. _ESAIM: Proceedings and Surveys_,
55 * 48:156--168, 2015.
56 *
57 * \see R. von Rickenbach. _Boundary Element Methods for Shape
58 * Optimisation in Homogenisation_. MSc thesis, Universität Basel, 2021.
59 */
60 5 Eigen::VectorXd getCoefficients(double precision) {
61
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 Eigen::VectorXd diff;
62
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 Eigen::RowVectorXd difft;
63
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 Eigen::Vector3d v;
64
3/6
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
 5 Eigen::MatrixXd spherical_values_pre, spherical_values_pst, Dsolid_values_pre,
65
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 Dsolid_values_pst;
66
67 5 unsigned int deg = getDegree(precision);
68 5 unsigned int Msquare = (POINT_DEGREE + 1) * (POINT_DEGREE + 1);
69
70 unsigned int m, n, k;
71 double scale, fac, norm;
72
73
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 GaussSquare<POINT_DEGREE> GS;
74
2/4
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
 5 Eigen::MatrixXd xs = GS[POINT_DEGREE].xi_;
75
3/6
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 5 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 5 times.
✗ Branch 10 not taken.
 5 xs -= 0.5 * Eigen::MatrixXd::Ones(xs.rows(), xs.cols());
76
77
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 Eigen::Vector3d ex(1.0, 0.0, 0.0);
78
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 Eigen::Vector3d ey(0.0, 1.0, 0.0);
79
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 Eigen::Vector3d ez(0.0, 0.0, 1.0);
80
81
1/2
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 5 Eigen::MatrixXd ps_left(3, xs.cols());
82
4/8
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 5 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 5 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 5 times.
✗ Branch 12 not taken.
 5 ps_left.row(0) = -0.5 * Eigen::VectorXd::Ones(xs.cols());
83
3/6
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 5 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 5 times.
✗ Branch 10 not taken.
 5 ps_left.block(1, 0, 2, xs.cols()) = xs.block(0, 0, 2, xs.cols());
84
85
1/2
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 5 Eigen::MatrixXd ps_front(3, xs.cols());
86
3/6
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
 5 ps_front.row(0) = xs.row(0);
87
4/8
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 5 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 5 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 5 times.
✗ Branch 12 not taken.
 5 ps_front.row(1) = -0.5 * Eigen::VectorXd::Ones(xs.cols());
88
3/6
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
 5 ps_front.row(2) = xs.row(1);
89
90
1/2
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 5 Eigen::MatrixXd ps_bottom(3, xs.cols());
91
3/6
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 5 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 5 times.
✗ Branch 10 not taken.
 5 ps_bottom.block(0, 0, 2, xs.cols()) = xs.block(0, 0, 2, xs.cols());
92
4/8
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 5 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 5 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 5 times.
✗ Branch 12 not taken.
 5 ps_bottom.row(2) = -0.5 * Eigen::VectorXd::Ones(xs.cols());
93
94
4/8
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 5 times.
✗ Branch 11 not taken.
 10 Eigen::VectorXd displacement = getDisplacement(ps_left, ps_front, ps_bottom);
95
96
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 Eigen::MatrixXd systemMatrix(6 * Msquare, ((deg + 1) * (deg + 2)) / 2 - 1);
97
98
2/2
✓ Branch 0 taken 2205 times.
✓ Branch 1 taken 5 times.
 2210 for (k = 0; k < Msquare; k++) {
99 /* left - right difference */
100
2/4
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
 2205 v = ps_left.col(k);
101
1/2
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
 2205 norm = v.norm();
102
103
3/6
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
 2205 spherical_values_pre = spherical_harmonics_full(v, deg);
104
4/8
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2205 times.
✗ Branch 11 not taken.
 2205 spherical_values_pst = spherical_harmonics_full(v + ex, deg);
105
106
4/8
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2205 times.
✗ Branch 11 not taken.
 2205 Dsolid_values_pre = Dsolid_harmonics_full(v, deg, spherical_values_pre);
107
4/8
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2205 times.
✗ Branch 11 not taken.
 4410 Dsolid_values_pst = Dsolid_harmonics_full(v + ex, deg,
108
1/2
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
 2205 spherical_values_pst);
109
110
2/2
✓ Branch 0 taken 26460 times.
✓ Branch 1 taken 2205 times.
 28665 for (n = 1; n <= deg; n++) {
111 26460 scale = pow(norm, n);
112 diff = scale
113
2/4
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26460 times.
✗ Branch 5 not taken.
 52920 * (spherical_values_pre.block(n * n + n, 0, n + 1, 1)
114
3/6
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26460 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 26460 times.
✗ Branch 8 not taken.
 79380 - spherical_values_pst.block(n * n + n, 0, n + 1, 1));
115
1/2
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
 26460 difft = Dsolid_values_pre.block(0, n * n + n, 1, n + 1)
116
3/6
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26460 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 26460 times.
✗ Branch 8 not taken.
 52920 - Dsolid_values_pst.block(0, n * n + n, 1, n + 1);
117
118 /* adapt the scaling */
119
2/4
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26460 times.
✗ Branch 5 not taken.
 26460 diff.segment(1, n) *= 2.0;
120
2/4
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26460 times.
✗ Branch 5 not taken.
 26460 difft.segment(1, n) *= 2.0;
121
122
3/6
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26460 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 26460 times.
✗ Branch 8 not taken.
 26460 systemMatrix.block(k, (n * (n + 1)) / 2 - 1, 1, n + 1) = diff.transpose();
123
2/4
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26460 times.
✗ Branch 5 not taken.
 26460 systemMatrix.block(k + Msquare, (n * (n + 1)) / 2 - 1, 1, n + 1) = difft;
124 }
125
126 /* front - back difference */
127
2/4
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
 2205 v = ps_front.col(k);
128
1/2
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
 2205 norm = v.norm();
129
130
3/6
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
 2205 spherical_values_pre = spherical_harmonics_full(v, deg);
131
4/8
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2205 times.
✗ Branch 11 not taken.
 2205 spherical_values_pst = spherical_harmonics_full(v + ey, deg);
132
133
4/8
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2205 times.
✗ Branch 11 not taken.
 2205 Dsolid_values_pre = Dsolid_harmonics_full(v, deg, spherical_values_pre);
134
4/8
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2205 times.
✗ Branch 11 not taken.
 4410 Dsolid_values_pst = Dsolid_harmonics_full(v + ey, deg,
135
1/2
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
 2205 spherical_values_pst);
136
137
2/2
✓ Branch 0 taken 26460 times.
✓ Branch 1 taken 2205 times.
 28665 for (n = 1; n <= deg; n++) {
138 26460 scale = pow(norm, n);
139 diff = scale
140
2/4
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26460 times.
✗ Branch 5 not taken.
 52920 * (spherical_values_pre.block(n * n + n, 0, n + 1, 1)
141
3/6
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26460 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 26460 times.
✗ Branch 8 not taken.
 79380 - spherical_values_pst.block(n * n + n, 0, n + 1, 1));
142
1/2
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
 26460 difft = Dsolid_values_pre.block(1, n * n + n, 1, n + 1)
143
3/6
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26460 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 26460 times.
✗ Branch 8 not taken.
 52920 - Dsolid_values_pst.block(1, n * n + n, 1, n + 1);
144
145 /* adapt the scaling */
146
2/4
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26460 times.
✗ Branch 5 not taken.
 26460 diff.segment(1, n) *= 2.0;
147
2/4
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26460 times.
✗ Branch 5 not taken.
 26460 difft.segment(1, n) *= 2.0;
148
149
1/2
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
 26460 systemMatrix.block(k + 2 * Msquare, (n * (n + 1)) / 2 - 1, 1, n + 1) =
150
2/4
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26460 times.
✗ Branch 5 not taken.
 52920 diff.transpose();
151
1/2
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
 52920 systemMatrix.block(k + 3 * Msquare, (n * (n + 1)) / 2 - 1, 1, n + 1) =
152
1/2
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
 26460 difft;
153 }
154
155 /* bottom - top difference */
156
2/4
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
 2205 v = ps_bottom.col(k);
157
1/2
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
 2205 norm = v.norm();
158
159
3/6
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
 2205 spherical_values_pre = spherical_harmonics_full(v, deg);
160
4/8
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2205 times.
✗ Branch 11 not taken.
 2205 spherical_values_pst = spherical_harmonics_full(v + ez, deg);
161
162
4/8
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2205 times.
✗ Branch 11 not taken.
 2205 Dsolid_values_pre = Dsolid_harmonics_full(v, deg, spherical_values_pre);
163
4/8
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2205 times.
✗ Branch 11 not taken.
 4410 Dsolid_values_pst = Dsolid_harmonics_full(v + ez, deg,
164
1/2
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
 2205 spherical_values_pst);
165
166
2/2
✓ Branch 0 taken 26460 times.
✓ Branch 1 taken 2205 times.
 28665 for (n = 1; n <= deg; n++) {
167 26460 scale = pow(norm, n);
168 diff = scale
169
2/4
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26460 times.
✗ Branch 5 not taken.
 52920 * (spherical_values_pre.block(n * n + n, 0, n + 1, 1)
170
3/6
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26460 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 26460 times.
✗ Branch 8 not taken.
 79380 - spherical_values_pst.block(n * n + n, 0, n + 1, 1));
171
1/2
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
 26460 difft = Dsolid_values_pre.block(2, n * n + n, 1, n + 1)
172
3/6
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26460 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 26460 times.
✗ Branch 8 not taken.
 52920 - Dsolid_values_pst.block(2, n * n + n, 1, n + 1);
173
174 /* adapt the scaling */
175
2/4
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26460 times.
✗ Branch 5 not taken.
 26460 diff.segment(1, n) *= 2.0;
176
2/4
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26460 times.
✗ Branch 5 not taken.
 26460 difft.segment(1, n) *= 2.0;
177
178
1/2
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
 26460 systemMatrix.block(k + 4 * Msquare, (n * (n + 1)) / 2 - 1, 1, n + 1) =
179
2/4
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26460 times.
✗ Branch 5 not taken.
 52920 diff.transpose();
180
1/2
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
 52920 systemMatrix.block(k + 5 * Msquare, (n * (n + 1)) / 2 - 1, 1, n + 1) =
181
1/2
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
 26460 difft;
182 }
183 }
184
185 /* solve the system */
186
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 Eigen::VectorXd coeffs(((deg + 1) * (deg + 2)) / 2);
187
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 coeffs.setZero();
188
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 coeffs.segment(1, ((deg + 1) * (deg + 2)) / 2 - 1) =
189
4/8
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 5 times.
✗ Branch 11 not taken.
 10 systemMatrix.colPivHouseholderQr().solve(-displacement);
190
191 /* Copy the stuff into the full Coefficient list */
192
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 Eigen::VectorXd coeffs_full((deg + 1) * (deg + 1));
193
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 coeffs_full(0) = 0;
194
2/2
✓ Branch 0 taken 60 times.
✓ Branch 1 taken 5 times.
 65 for (n = 1; n <= deg; n++) {
195
2/4
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 60 times.
✗ Branch 5 not taken.
 60 coeffs_full(n * n + n) = coeffs((n * (n + 1)) / 2);
196
2/2
✓ Branch 0 taken 390 times.
✓ Branch 1 taken 60 times.
 450 for (m = 1; m <= n; m++) {
197
2/4
✓ Branch 1 taken 390 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 390 times.
✗ Branch 5 not taken.
 390 coeffs_full(n * n + n + m) = coeffs((n * (n + 1)) / 2 + m);
198
2/4
✓ Branch 1 taken 390 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 390 times.
✗ Branch 5 not taken.
 390 coeffs_full(n * n + n - m) = coeffs((n * (n + 1)) / 2 + m);
199 }
200 }
201
202 /* calculate the first coefficient */
203
6/12
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 5 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 5 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 5 times.
✗ Branch 17 not taken.
 5 coeffs_full(0) = calculateFirstCoefficient(coeffs_full, deg, ps_left,
204 ps_front, ps_bottom);
205
206 10 return coeffs_full;
207 5 }
208
209 /**
210 * \brief Returns the degree of the sphericals expansion
211 * given a precision. Can be extended, use even numbers only!
212 *
213 * \see P. Cazeaux and O. Zahm. A fast boundary element method
214 * for the solution of periodic many-inclusion problems via
215 * hierarchical matrix techniques. _ESAIM: Proceedings and Surveys_,
216 * 48:156--168, 2015.
217 *
218 * \see R. von Rickenbach. _Boundary Element Methods for Shape
219 * Optimisation in Homogenisation_. MSc thesis, Universität Basel, 2021.
220 */
221 10 inline unsigned int getDegree(double precision) {
222
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 10 if (precision > 1e-4) {
223 ✗ return 4;
224
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 10 } else if (precision > 1e-6) {
225 ✗ return 8;
226 } else {
227 10 return 12;
228 }
229 }
230
231 /**
232 * \brief Returns the right-hand side for the homogenised coefficient calculation.
233 *
234 * @param ps_l: Points on the left side, i.e. \f$ \subseteq \{-0.5\} \times [-0.5, 0.5] \times [-0.5, 0.5] \f$.
235 * @param ps_f: Points on the front side, i.e. \f$ \subseteq [-0.5, 0.5] \times \{-0.5\} \times [-0.5, 0.5] \f$.
236 * @param ps_b: Points on the bottom side, i.e. \f$ \subseteq [-0.5, 0.5] \times [-0.5, 0.5] \times \{-0.5\} \f$.
237 */
238 5 Eigen::VectorXd getDisplacement(Eigen::MatrixXd ps_l, Eigen::MatrixXd ps_f,
239 Eigen::MatrixXd ps_b) {
240
241 13230 std::function<double(Eigen::Vector3d)> u = [](Eigen::Vector3d in) {
242
1/2
✓ Branch 2 taken 13230 times.
✗ Branch 3 not taken.
 13230 return k_mod(in);
243 5 };
244 13230 std::function<Eigen::Vector3d(Eigen::Vector3d)> Du = [](Eigen::Vector3d in) {
245
1/2
✓ Branch 2 taken 13230 times.
✗ Branch 3 not taken.
 13230 return Dk_mod(in);
246 5 };
247
248 5 unsigned int Msquare = (POINT_DEGREE + 1) * (POINT_DEGREE + 1);
249 unsigned int k;
250
251
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 Eigen::Vector3d ex(1.0, 0.0, 0.0);
252
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 Eigen::Vector3d ey(0.0, 1.0, 0.0);
253
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 Eigen::Vector3d ez(0.0, 0.0, 1.0);
254
255
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 Eigen::Vector3d tmp;
256
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 Eigen::VectorXd d(6 * Msquare); /* the return vector */
257
258 /* left - right displacement */
259
2/2
✓ Branch 0 taken 2205 times.
✓ Branch 1 taken 5 times.
 2210 for (k = 0; k < Msquare; k++) {
260
9/18
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2205 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2205 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 2205 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 2205 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 2205 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 2205 times.
✗ Branch 26 not taken.
 2205 tmp = Du(ps_l.col(k)) - Du(ps_l.col(k) + ex);
261
8/16
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2205 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2205 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 2205 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 2205 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 2205 times.
✗ Branch 23 not taken.
 2205 d(k) = u(ps_l.col(k)) - u(ps_l.col(k) + ex);
262
2/4
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
 2205 d(k + Msquare) = tmp(0);
263 }
264
265 /* front - back displacement */
266
2/2
✓ Branch 0 taken 2205 times.
✓ Branch 1 taken 5 times.
 2210 for (k = 0; k < Msquare; k++) {
267
9/18
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2205 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2205 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 2205 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 2205 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 2205 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 2205 times.
✗ Branch 26 not taken.
 2205 tmp = Du(ps_f.col(k)) - Du(ps_f.col(k) + ey);
268
8/16
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2205 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2205 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 2205 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 2205 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 2205 times.
✗ Branch 23 not taken.
 2205 d(k + 2 * Msquare) = u(ps_f.col(k)) - u(ps_f.col(k) + ey);
269
2/4
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
 2205 d(k + 3 * Msquare) = tmp(1);
270 }
271
272 /* bottom - top displacement */
273
2/2
✓ Branch 0 taken 2205 times.
✓ Branch 1 taken 5 times.
 2210 for (k = 0; k < Msquare; k++) {
274
9/18
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2205 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2205 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 2205 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 2205 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 2205 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 2205 times.
✗ Branch 26 not taken.
 2205 tmp = Du(ps_b.col(k)) - Du(ps_b.col(k) + ez);
275
8/16
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2205 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2205 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 2205 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 2205 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 2205 times.
✗ Branch 23 not taken.
 2205 d(k + 4 * Msquare) = u(ps_b.col(k)) - u(ps_b.col(k) + ez);
276
2/4
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
 2205 d(k + 5 * Msquare) = tmp(2);
277 }
278
279 10 return d;
280 5 }
281
282 /**
283 * \brief Returns the modified kernel \f$ \sum_{m \in \{-1, 0, 1 \}^3} \frac{1}{4 \pi |x-m|}
284 * + \frac{|x|^2}{6} \f$.
285 */
286 3237941 inline double k_mod(Eigen::Vector3d in) {
287 3237941 double r = 0.0;
288 int i, j, k;
289
1/2
✓ Branch 1 taken 3237941 times.
✗ Branch 2 not taken.
 3237941 Eigen::Vector3d m;
290
291
2/2
✓ Branch 0 taken 9713823 times.
✓ Branch 1 taken 3237941 times.
 12951764 for (i = -1; i <= 1; i++) {
292
2/2
✓ Branch 0 taken 29141469 times.
✓ Branch 1 taken 9713823 times.
 38855292 for (j = -1; j <= 1; j++) {
293
2/2
✓ Branch 0 taken 87424407 times.
✓ Branch 1 taken 29141469 times.
 116565876 for (k = -1; k <= 1; k++) {
294
1/2
✓ Branch 1 taken 87424407 times.
✗ Branch 2 not taken.
 87424407 m = Eigen::Vector3d(i, j, k);
295
2/4
✓ Branch 1 taken 87424407 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 87424407 times.
✗ Branch 5 not taken.
 87424407 r += 1.0 / ((in - m).norm());
296 }
297 }
298 }
299
300 3237941 r /= (4 * M_PI);
301
302 /* the part to ensure the vanishing mean on the Laplacian */
303
1/2
✓ Branch 1 taken 3237941 times.
✗ Branch 2 not taken.
 3237941 r += (in.dot(in)) / 6.0;
304
305 3237941 return r;
306 }
307
308 /**
309 * \brief Returns the gradient of the modified kernel \f$ - \sum_{m \in \{-1, 0, 1\}^3}
310 * \frac{x-m}{|x-m|^3} + \frac{x}{3} \f$.
311 */
312 74437 inline Eigen::Vector3d Dk_mod(Eigen::Vector3d in) {
313
2/4
✓ Branch 1 taken 74437 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 74437 times.
✗ Branch 5 not taken.
 74437 Eigen::Vector3d r, s;
314 double snorm;
315
1/2
✓ Branch 1 taken 74437 times.
✗ Branch 2 not taken.
 74437 r.setZero();
316
317 int i, j, k;
318
2/2
✓ Branch 0 taken 223311 times.
✓ Branch 1 taken 74437 times.
 297748 for (i = -1; i <= 1; i++) {
319
2/2
✓ Branch 0 taken 669933 times.
✓ Branch 1 taken 223311 times.
 893244 for (j = -1; j <= 1; j++) {
320
2/2
✓ Branch 0 taken 2009799 times.
✓ Branch 1 taken 669933 times.
 2679732 for (k = -1; k <= 1; k++) {
321
3/6
✓ Branch 1 taken 2009799 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2009799 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2009799 times.
✗ Branch 8 not taken.
 2009799 s = in - Eigen::Vector3d(i, j, k);
322
1/2
✓ Branch 1 taken 2009799 times.
✗ Branch 2 not taken.
 2009799 snorm = s.norm();
323
2/4
✓ Branch 1 taken 2009799 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2009799 times.
✗ Branch 5 not taken.
 2009799 r -= s / (snorm * snorm * snorm);
324 }
325 }
326 }
327
328
1/2
✓ Branch 1 taken 74437 times.
✗ Branch 2 not taken.
 74437 r /= (4.0 * M_PI);
329
330 /* the part to ensure the vanishing mean on the Laplacian */
331
2/4
✓ Branch 1 taken 74437 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 74437 times.
✗ Branch 5 not taken.
 74437 r += in / 3.0;
332
333 148874 return r;
334 }
335
336 /**
337 * \brief Calculates the first coefficient such that the mean of the
338 * kernel vanishes.
339 *
340 * @param cs: The already calculated other coefficients
341 * @param deg: The degree
342 * @param ps_l: Gauss quadrature points on the left side of the cube
343 * @param ps_f: Gauss quadrature points on the front side of the cube
344 * @param ps_b: Gauss quadrature points on the bottom side of the cube
345 */
346 5 double calculateFirstCoefficient(Eigen::VectorXd cs, unsigned int deg,
347 Eigen::MatrixXd ps_l, Eigen::MatrixXd ps_f, Eigen::MatrixXd ps_b) {
348 5 double res = 0.0;
349
350
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 GaussSquare<POINT_DEGREE> GS;
351
2/4
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
 5 Eigen::VectorXd ws = GS[POINT_DEGREE].w_;
352
353
1/2
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 5 Eigen::VectorXd cs_tmp(cs.rows());
354
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 cs_tmp.setZero();
355
356 unsigned int k, n;
357 double norm;
358
359
2/2
✓ Branch 1 taken 2205 times.
✓ Branch 2 taken 5 times.
 2210 for (k = 0; k < ps_l.cols(); k++) {
360
2/4
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
 2205 norm = ps_l.col(k).norm(); /* equal for ps_f and ps_b */
361
2/2
✓ Branch 0 taken 26460 times.
✓ Branch 1 taken 2205 times.
 28665 for (n = 1; n <= deg; n++) {
362
1/2
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
 26460 cs_tmp.segment(n * n, 2 * n + 1) = pow(norm, n)
363
3/6
✓ Branch 1 taken 26460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26460 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 26460 times.
✗ Branch 8 not taken.
 52920 * cs.segment(n * n, 2 * n + 1);
364 }
365
366
1/2
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
 2205 res += ws(k)
367
3/6
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
 2205 * (k_mod(ps_l.col(k))
368
5/10
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2205 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2205 times.
✗ Branch 14 not taken.
 2205 + evaluate_sphericals(ps_l.col(k) / norm, cs_tmp, deg));
369
1/2
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
 2205 res += ws(k)
370
3/6
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
 2205 * (k_mod(ps_f.col(k))
371
5/10
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2205 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2205 times.
✗ Branch 14 not taken.
 2205 + evaluate_sphericals(ps_f.col(k) / norm, cs_tmp, deg));
372
1/2
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
 2205 res += ws(k)
373
3/6
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
 2205 * (k_mod(ps_b.col(k))
374
5/10
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2205 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2205 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2205 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2205 times.
✗ Branch 14 not taken.
 2205 + evaluate_sphericals(ps_b.col(k) / norm, cs_tmp, deg));
375 }
376
377 5 res /= 3.0;
378
379 5 res += (1.0 / 24);
380
381 5 return -2.0 * sqrt(M_PI) * res;
382 5 }
383
384 } /* namespace Bembel */
385
386 #endif // BEMBEL_SRC_HOMOGENISEDLAPLACE_COEFFICIENTS_HPP_
387