GCC Code Coverage Report

Directory: Bembel/src/
File: Bembel/src/Spline/DeBoorTP.hpp
Date: 2024-03-19 14:38:05
Exec Total Coverage
Lines: 16 16 100.0%
Functions: 2 2 100.0%
Branches: 2 2 100.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 #ifndef BEMBEL_SRC_SPLINE_DEBOORTP_HPP_
12 #define BEMBEL_SRC_SPLINE_DEBOORTP_HPP_
13
14 namespace Bembel {
15 namespace Spl {
16 /**
17 * \ingroup Spline
18 * \brief A "by the book" implementation of the derivatives and TP-algos based
19 * on the DeBoor Recursion.
20 **/
21 template <typename T>
22 3670 Eigen::Matrix<T, -1, -1> DeBoorDer(
23 Eigen::Matrix<T, -1, -1> const &control_points,
24 std::vector<double> const &knot,
25 std::vector<double> const &evaluation_points) noexcept {
26 3670 const int control_points_cols = control_points.cols();
27 3670 const int polynomial_degree = knot.size() - control_points_cols - 1;
28 3670 Eigen::Matrix<T, -1, -1> temp(control_points.rows(), control_points_cols - 1);
29
2/2
✓ Branch 0 taken 48510 times.
✓ Branch 1 taken 3670 times.
 52180 for (int i = control_points_cols - 1; 0 <= --i;) {
30 97020 temp.col(i) = (polynomial_degree) /
31 48510 (knot[i + polynomial_degree + 1] - knot[i + 1]) *
32 97020 (control_points.col(i + 1) - control_points.col(i));
33 }
34
35 3670 std::vector<double> tempknt(knot.begin() + 1, knot.end() - 1);
36 3670 return DeBoor(temp, tempknt, evaluation_points);
37 3670 }
38
39 template <typename T>
40 std::vector<Eigen::Matrix<T, -1, -1>> DeBoorDer(
41 std::vector<Eigen::Matrix<T, -1, -1>> const &control_points,
42 std::vector<double> const &knot,
43 std::vector<double> const &evaluation_points) noexcept {
44 const int control_points_cols = control_points[0].cols();
45 const int dimension = control_points.size();
46 for (int ll = 0; ll < dimension; ll++) {
47 assert(control_points[0].cols() == control_points[ll].cols() &&
48 control_points[0].rows() == control_points[ll].rows());
49 }
50 const int polynomial_degree = knot.size() - control_points_cols - 1;
51 std::vector<Eigen::Matrix<T, -1, -1>> temp(dimension);
52 for (int ll = 0; ll < dimension; ll++) {
53 temp[ll].resize(control_points[ll].rows(), control_points_cols - 1);
54 }
55
56 for (int i = control_points_cols - 1; 0 <= --i;) {
57 double factor =
58 (polynomial_degree) / (knot[i + polynomial_degree + 1] - knot[i + 1]);
59 for (int ll = 0; ll < dimension; ll++)
60 temp[ll].col(i) =
61 factor * (control_points[ll].col(i + 1) - control_points[ll].col(i));
62 }
63
64 std::vector<double> tempknt(knot.begin() + 1, knot.end() - 1);
65
66 return DeBoor(temp, tempknt, evaluation_points);
67 }
68
69 template <typename T>
70 std::vector<Eigen::Matrix<T, -1, -1>> deBoorDerGiveData(
71 std::vector<Eigen::Matrix<T, -1, -1>> const &control_points,
72 std::vector<double> const &knot) noexcept {
73 const int control_points_cols = control_points[0].cols();
74 const int dimension = control_points.size();
75 for (int ll = 0; ll < dimension; ll++) {
76 assert(control_points[0].cols() == control_points[ll].cols() &&
77 control_points[0].rows() == control_points[ll].rows());
78 }
79 const int polynomial_degree = knot.size() - control_points_cols - 1;
80 std::vector<Eigen::Matrix<T, -1, -1>> temp(dimension);
81 for (int ll = 0; ll < dimension; ll++) {
82 temp[ll].resize(control_points[ll].rows(), control_points_cols - 1);
83 }
84
85 for (int i = control_points_cols - 1; 0 <= --i;) {
86 double factor =
87 (polynomial_degree) / (knot[i + polynomial_degree + 1] - knot[i + 1]);
88 for (int ll = 0; ll < dimension; ll++)
89 temp[ll].col(i) =
90 factor * (control_points[ll].col(i + 1) - control_points[ll].col(i));
91 }
92
93 return temp;
94 }
95 ////////////////////////////////////////////////////////////////////////////////
96 /// Simple TP Bsplines
97 ////////////////////////////////////////////////////////////////////////////////
98 template <typename T>
99 33 Eigen::Matrix<T, -1, -1> DeBoorTP(
100 Eigen::Matrix<T, -1, -1> const &control_points,
101 std::vector<double> const &knots_x, std::vector<double> const &knots_y,
102 std::vector<double> const &evaluation_points_x,
103 std::vector<double> const &evaluation_points_y) noexcept {
104 33 Eigen::Matrix<T, -1, -1> tmp =
105 66 DeBoor(control_points, knots_x, evaluation_points_x).transpose();
106 33 return DeBoor(tmp, knots_y, evaluation_points_y);
107 33 }
108
109 template <typename T>
110 std::vector<Eigen::Matrix<T, -1, -1>> DeBoorTP(
111 std::vector<Eigen::Matrix<T, -1, -1>> const &control_points,
112 std::vector<double> const &knots_x, std::vector<double> const &knots_y,
113 std::vector<double> const &evaluation_points_x,
114 std::vector<double> const &evaluation_points_y) noexcept {
115 std::vector<Eigen::Matrix<T, -1, -1>> tmp =
116 DeBoor(control_points, knots_x, evaluation_points_x);
117 for (int ll = tmp.size() - 1; ll >= 0; ll--)
118 tmp[ll] = tmp[ll].transpose().eval();
119 return DeBoor(tmp, knots_y, evaluation_points_y);
120 }
121 ////////////////////////////////////////////////////////////////////////////////
122 // TP Der with direction declaration
123 ////////////////////////////////////////////////////////////////////////////////
124 template <typename T>
125 Eigen::Matrix<T, -1, -1> DeBoorTPDer(
126 Eigen::Matrix<T, -1, -1> const &control_points,
127 std::vector<double> const &knots_x, std::vector<double> const &knots_y,
128 std::vector<double> const &evaluation_points_x,
129 std::vector<double> const &evaluation_points_y, bool const &x_to_be_derived,
130 bool const &y_to_be_derived) noexcept {
131 Eigen::Matrix<T, -1, -1> tmp =
132 x_to_be_derived
133 ? DeBoorDer(control_points, knots_x, evaluation_points_x).transpose()
134 : DeBoor(control_points, knots_x, evaluation_points_x).transpose();
135 return y_to_be_derived ? DeBoorDer(tmp, knots_y, evaluation_points_y)
136 : DeBoor(tmp, knots_y, evaluation_points_y);
137 }
138
139 template <typename T>
140 std::vector<Eigen::Matrix<T, -1, -1>> DeBoorTPDer(
141 std::vector<Eigen::Matrix<T, -1, -1>> const &control_points,
142 std::vector<double> const &knots_x, std::vector<double> const &knots_y,
143 std::vector<double> const &evaluation_points_x,
144 std::vector<double> const &evaluation_points_y, bool const &x_to_be_derived,
145 bool const &y_to_be_derived) noexcept {
146 std::vector<Eigen::Matrix<T, -1, -1>> tmp =
147 x_to_be_derived ? DeBoorDer(control_points, knots_x, evaluation_points_x)
148 : DeBoor(control_points, knots_x, evaluation_points_x);
149 for (int ll = tmp.size() - 1; ll >= 0; ll--)
150 tmp[ll] = tmp[ll].transpose().eval();
151 return y_to_be_derived ? DeBoorDer(tmp, knots_y, evaluation_points_y)
152 : DeBoor(tmp, knots_y, evaluation_points_y);
153 }
154 } // namespace Spl
155 } // namespace Bembel
156 #endif // BEMBEL_SRC_SPLINE_DEBOORTP_HPP_
157