GCC Code Coverage Report

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		1055782528	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		1055782528	return evaluation_point * BernsteinX<N - 1>(evaluation_point);
29			}
30
31			template <>
32		882863328	inline constexpr double BernsteinX<1>(double evaluation_point) noexcept {
33		882863328	return evaluation_point;
34			}
35			template <>
36		1350949858	inline constexpr double BernsteinX<0>(double evaluation_point) noexcept {
37		1350949858	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		2233813186	inline constexpr double Bernstein(double evaluation_point) noexcept {
46		2233813186	return Binomial<N, P>::value * BernsteinX<N>(evaluation_point) *
47		2233813186	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		5132	static inline T EvalCoefs(T *in, double evaluation_point) noexcept {
56		5132	return in[N] * Bernstein<N, P>(evaluation_point) +
57		5132	HiddenBernsteinClass<T, N - 1, P>::EvalCoefs(in, evaluation_point);
58			}
59		4180	static inline T EvalDerCoefs(T *in, double evaluation_point) noexcept {
60			return (
61		4180	(in[N + 1] - in[N]) * Bernstein<N, P>(evaluation_point) +
62		4180	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		619737620	static inline void EvalBasis(T *in, double evaluation_point) noexcept {
76		619737620	in[N] = Bernstein<N, P>(evaluation_point);
77		619737620	HiddenBernsteinClass<T, N - 1, P>::EvalBasis(in, evaluation_point);
78		619737620	return;
79			}
80		131558198	static inline void EvalDerBasis(T *in, double evaluation_point) noexcept {
81		131558198	in[N] = (P + 1) * (Bernstein<N - 1, P>(evaluation_point) -
82		131558198	Bernstein<N, P>(evaluation_point));
83		131558198	HiddenBernsteinClass<T, N - 1, P>::EvalDerBasis(in, evaluation_point);
84		131558198	return;
85			}
86			};
87
88			template <typename T, int P>
89			class HiddenBernsteinClass<T, 0, P> {
90			public:
91		974	static inline T EvalCoefs(T *in, double evaluation_point) noexcept {
92		974	return in[0] * Bernstein<0, P>(evaluation_point);
93			}
94		440	static inline T EvalDerCoefs(T *in, double evaluation_point) noexcept {
95			// P needs to be passed lower to avoid infinite recursion
96		440	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		510586804	static inline void EvalBasis(T *in, double evaluation_point) noexcept {
109		510586804	in[0] = Bernstein<0, P>(evaluation_point);
110		510586804	return;
111			}
112		420180820	static inline void EvalDerBasis(T *in, double evaluation_point) noexcept {
113			// P needs to be passed lower to avoid infinite recursion
114		420180820	in[0] = (-P - 1) * Bernstein<0, P>(evaluation_point);
115		420180820	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		718	T EvalBernstein(T *in, double evaluation_point) noexcept {
169		718	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		510586804	void EvalBernsteinBasis(T *in, double evaluation_point) noexcept {
199		510586804	HiddenBernsteinClass<T, P, P>::EvalBasis(in, evaluation_point);
200		510586804	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		210090410	void EvalBernsteinDerBasis(T *in, double evaluation_point) noexcept {
241		210090410	in[P] = P * Bernstein<P - 1, P - 1>(evaluation_point);
242		210090410	HiddenBernsteinClass<T, P - 1, P - 1>::EvalDerBasis(in, evaluation_point);
243		210090410	return;
244			}
245
246			} // namespace Basis
247			} // namespace Bembel
248			#endif // BEMBEL_SRC_SPLINE_BERNSTEIN_HPP_
249