GCC Code Coverage Report

Directory: Bembel/src/
File: Bembel/src/Spline/Bezierextraction.hpp
Date: 2024-09-30 07:01:38
Exec Total Coverage
Lines: 52 52 100.0%
Functions: 2 2 100.0%
Branches: 31 40 77.5%
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_SPLINE_BEZIEREXTRACTION_HPP_
12 #define BEMBEL_SRC_SPLINE_BEZIEREXTRACTION_HPP_
13
14 namespace Bembel {
15 namespace Spl {
16
17 /**
18 * \ingroup Spline
19 * \brief implements Bezier extraction via the solution of an interpolation
20 * problem.
21 */
22 template <typename T>
23 1 Eigen::SparseMatrix<T> MakeProjection(const std::vector<T> &x_knots,
24 const std::vector<T> &y_knots,
25 const std::vector<T> &x_unique_knots,
26 const std::vector<T> &y_unique_knots,
27 const int polynomial_degree_x,
28 const int polynomial_degree_y) noexcept {
29 // In here phi refers to the BIG bezier basis, and psi refers to the small
30 // B-Spline basis given by x_knots,y_knots and xp,yp.
31 // Sparse does weird things. We need not be that accurate, since the right
32 // coeefs will be >=.1
33 1 constexpr double sparseTol = 0.001;
34 // Number of functions in X,Y and total
35 1 const int size_phi_x = (x_unique_knots.size() - 1) * polynomial_degree_x;
36 1 const int size_phi_y = (y_unique_knots.size() - 1) * polynomial_degree_y;
37 1 const int size_phi = size_phi_x * size_phi_y;
38 1 const int size_psi_x = x_knots.size() - polynomial_degree_x;
39 1 const int size_psi_y = y_knots.size() - polynomial_degree_y;
40 1 const int size_psi = size_psi_x * size_psi_y;
41
42 // Interpolation points
43 1 const auto xmask = MakeInterpolationMask(polynomial_degree_x);
44 1 const auto ymaks = MakeInterpolationMask(polynomial_degree_y);
45
46 1 const auto xpoint = MakeInterpolationPoints(x_unique_knots, xmask);
47 1 const auto ypoint = MakeInterpolationPoints(y_unique_knots, ymaks);
48
49
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
 1 assert(((int)(xmask.size() * (x_unique_knots.size() - 1))) == size_phi_x);
50
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
 1 assert(((int)(xpoint.size())) == size_phi_x);
51
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
 1 assert(((int)(ypoint.size())) == size_phi_y);
52
53 1 Eigen::Matrix<T, -1, -1> tempsolver =
54 Eigen::Matrix<T, -1, -1>::Zero(size_phi, size_phi);
55
56
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 double *coefficients_x = new double[polynomial_degree_x];
57
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 double *coefficients_y = new double[polynomial_degree_y];
58
59
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
 3 for (int i = 0; i < polynomial_degree_x; i++) coefficients_x[i] = 0;
60
61
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
 3 for (int i = 0; i < polynomial_degree_y; i++) coefficients_y[i] = 0;
62
63 512 auto BezBasis = [](std::vector<double> uniq, int deg, int pos, double pt,
64 double *coef) {
65 512 std::div_t loc = std::div(pos, deg);
66
6/6
✓ Branch 1 taken 384 times.
✓ Branch 2 taken 128 times.
✓ Branch 4 taken 128 times.
✓ Branch 5 taken 256 times.
✓ Branch 6 taken 256 times.
✓ Branch 7 taken 256 times.
 512 if (pt > uniq[loc.quot + 1] || pt < uniq[loc.quot]) {
67 256 return 0.;
68 } else {
69 256 coef[loc.rem] = 1;
70
1/2
✓ Branch 2 taken 256 times.
✗ Branch 3 not taken.
 256 double out = Bembel::Basis::ShapeFunctionHandler::evalCoef(
71 256 deg - 1, coef, Rescale(pt, uniq[loc.quot], uniq[loc.quot + 1]));
72 256 coef[loc.rem] = 0;
73 256 return out;
74 }
75 };
76
77 // A matrix, each row corresponds to one tp-basisfunction, each col to one
78 // tp-interpolation point
79
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
 5 for (int iy = 0; iy < size_phi_y; iy++) {
80
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 4 times.
 20 for (int ix = 0; ix < size_phi_x; ix++) {
81
2/2
✓ Branch 0 taken 64 times.
✓ Branch 1 taken 16 times.
 80 for (int y = 0; y < size_phi_y; y++) {
82
2/2
✓ Branch 0 taken 256 times.
✓ Branch 1 taken 64 times.
 320 for (int x = 0; x < size_phi_x; x++) {
83 512 tempsolver(y * size_phi_x + x, ix * size_phi_y + iy) =
84 512 BezBasis(x_unique_knots, polynomial_degree_x, ix, xpoint[x],
85 256 coefficients_x) *
86 256 BezBasis(y_unique_knots, polynomial_degree_y, iy, ypoint[y],
87 coefficients_y);
88 }
89 }
90 }
91 }
92
93
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 delete[] coefficients_x;
94
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 delete[] coefficients_y;
95
96 1 Eigen::FullPivLU<Eigen::Matrix<double, -1, -1>> lu_decomp(tempsolver);
97
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
 1 assert(lu_decomp.rank() == size_phi);
98 /// The inverse matrix, for the solition of the linear system to come.
99 1 Eigen::Matrix<T, -1, -1> solve = tempsolver.inverse();
100
101 1 Eigen::Matrix<T, -1, -1> Proj(size_phi, size_psi);
102 1 Eigen::Matrix<T, -1, -1> unit =
103 Eigen::Matrix<T, -1, -1>::Zero(size_psi_y, size_psi_x);
104
105 /// For each basisfunction phi_j solves for the right coeef such that \Sum
106
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
 4 for (int i = 0; i < size_psi_x; i++)
107
2/2
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 3 times.
 12 for (int j = 0; j < size_psi_y; j++) {
108 9 unit(j, i) = 1.;
109 9 Eigen::Matrix<T, -1, 1> smallBase =
110 (Unroll(DeBoorTP(unit, x_knots, y_knots, xpoint, ypoint)));
111 9 Proj.col(i * size_psi_y + j) = (solve * smallBase).transpose();
112 9 unit(j, i) = 0.;
113 }
114
115 1 return (Proj.sparseView()).pruned(sparseTol);
116 1 }
117 } // namespace Spl
118 } // namespace Bembel
119 #endif // BEMBEL_SRC_SPLINE_BEZIEREXTRACTION_HPP_
120