GCC Code Coverage Report

Directory: Bembel/src/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 99.8% 4044 / 0 / 4053
Functions: 99.7% 1176 / 0 / 1180
Branches: -% 0 / 0 / 0
Spline/Bernstein.hpp
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_SPLINE_BERNSTEIN_HPP_
13 #define BEMBEL_SRC_SPLINE_BERNSTEIN_HPP_
14
15 namespace Bembel {
16 namespace Basis {
17 /**
18 * \ingroup Spline
19 * \brief Template recursion to produce Bernstein polynomials. This is only
20 * limited by the binomial coefficient, see Pascal.hpp
21 */
22 template <int N>
23 527891264 inline constexpr double BernsteinX(double evaluation_point) noexcept {
24 #ifdef _spline_debug_flag_
25 assert((evaluation_point > -.0000001) && (evaluation_point < 1.0000001) &&
26 ("Function only valid for 0 <= x <= 1!"));
27 #endif
28 527891264 return evaluation_point * BernsteinX<N - 1>(evaluation_point);
29 }
30
31 template <>
32 882863776 inline constexpr double BernsteinX<1>(double evaluation_point) noexcept {
33 882863776 return evaluation_point;
34 }
35 template <>
36 1350950626 inline constexpr double BernsteinX<0>(double evaluation_point) noexcept {
37 1350950626 return 1.;
38 }
39 template <>
40 ✗ inline constexpr double BernsteinX<-1>(double evaluation_point) noexcept {
41 ✗ return 0.;
42 }
43
44 template <int N, int P>
45 1116907201 inline constexpr double Bernstein(double evaluation_point) noexcept {
46 1116907201 return Binomial<N, P>::value * BernsteinX<N>(evaluation_point) *
47 1116907201 BernsteinX<P - N>(1. - evaluation_point);
48 }
49 ////////////////////////////////////////////////////////////////////////////////
50 /// Hidden Classes
51 ////////////////////////////////////////////////////////////////////////////////
52 template <typename T, int N, int P>
53 class HiddenBernsteinClass {
54 public:
55 2566 static inline T EvalCoefs(T *in, double evaluation_point) noexcept {
56 2566 return in[N] * Bernstein<N, P>(evaluation_point) +
57 2566 HiddenBernsteinClass<T, N - 1, P>::EvalCoefs(in, evaluation_point);
58 }
59 2090 static inline T EvalDerCoefs(T *in, double evaluation_point) noexcept {
60 return (
61 2090 (in[N + 1] - in[N]) * Bernstein<N, P>(evaluation_point) +
62 2090 HiddenBernsteinClass<T, N - 1, P>::EvalDerCoefs(in, evaluation_point));
63 }
64 static inline void EvalBasisPEQ(T *in, double evaluation_point) noexcept {
65 in[N] += Bernstein<N, P>(evaluation_point);
66 HiddenBernsteinClass<T, N - 1, P>::EvalBasisPEQ(in, evaluation_point);
67 return;
68 }
69 static inline void EvalDerBasisPEQ(T *in, double evaluation_point) noexcept {
70 in[N] += (P + 1) * (Bernstein<N - 1, P>(evaluation_point) -
71 Bernstein<N, P>(evaluation_point));
72 HiddenBernsteinClass<T, N - 1, P>::EvalDerBasisPEQ(in, evaluation_point);
73 return;
74 }
75 309869034 static inline void EvalBasis(T *in, double evaluation_point) noexcept {
76 309869034 in[N] = Bernstein<N, P>(evaluation_point);
77 309869034 HiddenBernsteinClass<T, N - 1, P>::EvalBasis(in, evaluation_point);
78 309869034 return;
79 }
80 65779099 static inline void EvalDerBasis(T *in, double evaluation_point) noexcept {
81 65779099 in[N] = (P + 1) * (Bernstein<N - 1, P>(evaluation_point) -
82 65779099 Bernstein<N, P>(evaluation_point));
83 65779099 HiddenBernsteinClass<T, N - 1, P>::EvalDerBasis(in, evaluation_point);
84 65779099 return;
85 }
86 };
87
88 template <typename T, int P>
89 class HiddenBernsteinClass<T, 0, P> {
90 public:
91 487 static inline T EvalCoefs(T *in, double evaluation_point) noexcept {
92 487 return in[0] * Bernstein<0, P>(evaluation_point);
93 }
94 220 static inline T EvalDerCoefs(T *in, double evaluation_point) noexcept {
95 // P needs to be passed lower to avoid infinite recursion
96 220 return (in[1] - in[0]) * Bernstein<0, P>(evaluation_point);
97 }
98 static inline void EvalBasisPEQ(T *in, double evaluation_point) noexcept {
99 in[0] += Bernstein<0, P>(evaluation_point);
100 return;
101 }
102 static inline void EvalDerBasisPEQ(T *in, double evaluation_point) noexcept {
103 // P needs to be passed lower to avoid infinite recursion
104 in[0] += (-P - 1) * Bernstein<0, P>(evaluation_point);
105 return;
106 }
107
108 255293626 static inline void EvalBasis(T *in, double evaluation_point) noexcept {
109 255293626 in[0] = Bernstein<0, P>(evaluation_point);
110 255293626 return;
111 }
112 210090490 static inline void EvalDerBasis(T *in, double evaluation_point) noexcept {
113 // P needs to be passed lower to avoid infinite recursion
114 210090490 in[0] = (-P - 1) * Bernstein<0, P>(evaluation_point);
115 210090490 return;
116 }
117 };
118
119 // This specialization is needed to get a specialized recursion anchor for the
120 // case P = 0.
121 template <typename T, int P>
122 class HiddenBernsteinClass<T, -1, P> {
123 public:
124 static inline T EvalCoefs(T *in, double evaluation_point) noexcept {
125 (void)in;
126 (void)evaluation_point;
127 assert(
128 false &&
129 "Pos.A This should not happen. Something is wrong with the recursion");
130 };
131 ✗ static inline T EvalDerCoefs(T *in, double evaluation_point) noexcept {
132 // P needs to be passed lower to avoid infinite recursion
133 (void)in;
134 (void)evaluation_point;
135 ✗ return 0;
136 };
137 static inline void EvalBasis(T *in, double evaluation_point) noexcept {
138 (void)in;
139 (void)evaluation_point;
140 assert(
141 false &&
142 "Pos.C This should not happen. Something is wrong with the recursion");
143 };
144 ✗ static inline void EvalDerBasis(T *in, double evaluation_point) noexcept {
145 (void)in;
146 (void)evaluation_point;
147 // P needs to be passed lower to avoid infinite recursion
148 ✗ return;
149 };
150 static inline void EvalBasisPEQ(T *in, double evaluation_point) noexcept {
151 (void)in;
152 (void)evaluation_point;
153 assert(
154 false &&
155 "Pos.C This should not happen. Something is wrong with the recursion");
156 };
157 static inline void EvalDerBasisPEQ(T *in, double evaluation_point) noexcept {
158 (void)in;
159 (void)evaluation_point;
160 // P needs to be passed lower to avoid infinite recursion
161 return;
162 };
163 };
164 ////////////////////////////////////////////////////////////////////////////////
165 /// Evaluation Routines
166 ////////////////////////////////////////////////////////////////////////////////
167 template <typename T, int P>
168 487 T EvalBernstein(T *in, double evaluation_point) noexcept {
169 487 return HiddenBernsteinClass<T, P, P>::EvalCoefs(in, evaluation_point);
170 }
171
172 template <typename T, int P>
173 void EvalBernstein(T *in, const std::vector<double> &evaluation_points,
174 T *out) noexcept {
175 const int N = evaluation_points.size();
176 for (int i = 0; i < N; i++)
177 out[i] = HiddenBernsteinClass<T, P, P>::EvalCoefs(in, evaluation_points[i]);
178 return;
179 }
180
181 template <typename T, int P>
182 std::vector<T> EvalBernstein(
183 T *in, const std::vector<double> &evaluation_points) noexcept {
184 const int N = evaluation_points.size();
185 std::vector<double> out(N);
186 for (int i = 0; i < N; i++)
187 out[i] = HiddenBernsteinClass<T, P, P>::EvalCoefs(in, evaluation_points[i]);
188 return out;
189 }
190
191 template <typename T, int P>
192 void EvalBernsteinBasisPEQ(T *in, double evaluation_point) noexcept {
193 HiddenBernsteinClass<T, P, P>::EvalBasisPEQ(in, evaluation_point);
194 return;
195 }
196
197 template <typename T, int P>
198 255293626 void EvalBernsteinBasis(T *in, double evaluation_point) noexcept {
199 255293626 HiddenBernsteinClass<T, P, P>::EvalBasis(in, evaluation_point);
200 255293626 return;
201 }
202 ////////////////////////////////////////////////////////////////////////////////
203 /// Evaluation of the Derivatives
204 ////////////////////////////////////////////////////////////////////////////////
205 template <typename T, int P>
206 220 T EvalBernsteinDer(T *in, double evaluation_point) noexcept {
207 220 return P * HiddenBernsteinClass<T, P - 1, P - 1>::EvalDerCoefs(
208 220 in, evaluation_point);
209 }
210
211 template <typename T, int P>
212 void EvalBernsteinDer(T *in, const std::vector<double> &evaluation_points,
213 T *out) noexcept {
214 const int N = evaluation_points.size();
215 for (int i = 0; i < N; i++)
216 out[i] = P * HiddenBernsteinClass<T, P - 1, P - 1>::EvalDerCoefs(
217 in, evaluation_points[i]);
218 return;
219 }
220
221 template <typename T, int P>
222 std::vector<T> EvalBernsteinDer(
223 T *in, const std::vector<double> &evaluation_points) noexcept {
224 const int N = evaluation_points.size();
225 std::vector<double> out(N);
226 for (int i = 0; i < N; i++)
227 out[i] = P * HiddenBernsteinClass<T, P - 1, P - 1>::EvalDerCoefs(
228 in, evaluation_points[i]);
229 return out;
230 }
231
232 template <typename T, int P>
233 void EvalBernsteinDerBasisPEQ(T *in, double evaluation_point) noexcept {
234 in[P] += P * Bernstein<P - 1, P - 1>(evaluation_point);
235 HiddenBernsteinClass<T, P - 1, P - 1>::EvalDerBasisPEQ(in, evaluation_point);
236 return;
237 }
238
239 template <typename T, int P>
240 210090490 void EvalBernsteinDerBasis(T *in, double evaluation_point) noexcept {
241 210090490 in[P] = P * Bernstein<P - 1, P - 1>(evaluation_point);
242 210090490 HiddenBernsteinClass<T, P - 1, P - 1>::EvalDerBasis(in, evaluation_point);
243 210090490 return;
244 }
245
246 } // namespace Basis
247 } // namespace Bembel
248 #endif // BEMBEL_SRC_SPLINE_BERNSTEIN_HPP_
249