ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
bspline.hpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 The ESPResSo project
3 *
4 * This file is part of ESPResSo.
5 *
6 * ESPResSo is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * ESPResSo is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20#pragma once
21
22#include "sqr.hpp"
23#include "utils/device_qualifier.hpp"
24
25#if defined(__CUDACC__)
26#include <cfloat>
27#else
28#include <limits>
29#endif
30#include <stdexcept>
31
32#if defined(__clang__)
33#pragma clang diagnostic push
34#pragma clang diagnostic ignored "-Wimplicit-fallthrough"
35#elif defined(__GNUC__) or defined(__GNUG__)
36#pragma GCC diagnostic push
37#pragma GCC diagnostic ignored "-Wimplicit-fallthrough"
38#endif
39
40namespace Utils {
41/** @brief Formula of the B-spline. */
42template <int order, typename T>
43DEVICE_QUALIFIER auto bspline(int i, T x) -> T
44 requires((order > 0) and (order <= 7))
45{
46#if defined(__CUDACC__)
47 // LCOV_EXCL_START
48 [[maybe_unused]] auto constexpr tolerance = T(6) * (T)(FLT_EPSILON);
49 // LCOV_EXCL_STOP
50#else
51 [[maybe_unused]] auto constexpr tolerance =
52 T(6) * std::numeric_limits<T>::epsilon();
53#endif
54 DEVICE_ASSERT(i < order);
55 DEVICE_ASSERT(x >= T(-0.5) or (T(-0.5) - x) < tolerance);
56 DEVICE_ASSERT(x <= T(0.5) or (x - T(0.5)) < tolerance);
57
58 switch (order) {
59 case 1:
60 return T(1.);
61 case 2:
62 switch (i) {
63 case 0:
64 return T(0.5) - x;
65 case 1:
66 return T(0.5) + x;
67 }
68 case 3:
69 switch (i) {
70 case 0:
71 return T(0.5) * sqr(T(0.5) - x);
72 case 1:
73 return T(0.75) - sqr(x);
74 case 2:
75 return T(0.5) * sqr(T(0.5) + x);
76 }
77 case 4:
78 switch (i) {
79 case 0:
80 return (T(1.) + x * (T(-6.) + x * (T(12.) - x * T(8.)))) / T(48.);
81 case 1:
82 return (T(23.) + x * (T(-30.) + x * (T(-12.) + x * T(24.)))) / T(48.);
83 case 2:
84 return (T(23.) + x * (T(30.) + x * (T(-12.) - x * T(24.)))) / T(48.);
85 case 3:
86 return (T(1.) + x * (T(6.) + x * (T(12.) + x * T(8.)))) / T(48.);
87 }
88 case 5:
89 switch (i) {
90 case 0:
91 return (T(1.) +
92 x * (T(-8.) + x * (T(24.) + x * (T(-32.) + x * T(16.))))) /
93 T(384.);
94 case 1:
95 return (T(19.) +
96 x * (T(-44.) + x * (T(24.) + x * (T(16.) - x * T(16.))))) /
97 T(96.);
98 case 2:
99 return (T(115.) + x * x * (T(-120.) + x * x * T(48.))) / T(192.);
100 case 3:
101 return (T(19.) +
102 x * (T(44.) + x * (T(24.) + x * (T(-16.) - x * T(16.))))) /
103 T(96.);
104 case 4:
105 return (T(1.) + x * (T(8.) + x * (T(24.) + x * (T(32.) + x * T(16.))))) /
106 T(384.);
107 }
108 case 6:
109 switch (i) {
110 case 0:
111 return (T(1.) +
112 x * (T(-10.) +
113 x * (T(40.) + x * (T(-80.) + x * (T(80.) - x * T(32.)))))) /
114 T(3840.);
115 case 1:
116 return (T(237.) +
117 x * (T(-750.) +
118 x * (T(840.) +
119 x * (T(-240.) + x * (T(-240.) + x * T(160.)))))) /
120 T(3840.);
121 case 2:
122 return (T(841.) +
123 x * (T(-770.) +
124 x * (T(-440.) +
125 x * (T(560.) + x * (T(80.) - x * T(160.)))))) /
126 T(1920.);
127 case 3:
128 return (T(841.) +
129 x * (+T(770.) +
130 x * (T(-440.) +
131 x * (T(-560.) + x * (T(80.) + x * T(160.)))))) /
132 T(1920.);
133 case 4:
134 return (T(237.) +
135 x * (T(750.) +
136 x * (T(840.) +
137 x * (T(240.) + x * (T(-240.) - x * T(160.)))))) /
138 T(3840.);
139 case 5:
140 return (T(1.) +
141 x * (T(10.) +
142 x * (T(40.) + x * (T(80.) + x * (T(80.) + x * T(32.)))))) /
143 T(3840.);
144 }
145 case 7:
146 switch (i) {
147 case 0:
148 return (T(1.) +
149 x * (T(-12.) +
150 x * (T(60.) +
151 x * (T(-160.) +
152 x * (T(240.) + x * (T(-192.) + x * T(64.))))))) /
153 T(46080.);
154 case 1:
155 return (T(361.) +
156 x * (T(-1416.) +
157 x * (T(2220.) +
158 x * (T(-1600.) +
159 x * (T(240.) + x * (T(384.) - x * T(192.))))))) /
160 T(23040.);
161 case 2:
162 return (T(10543.) +
163 x * (T(-17340.) +
164 x * (T(4740.) +
165 x * (T(6880.) + x * (T(-4080.) +
166 x * (T(-960.) + x * T(960.))))))) /
167 T(46080.);
168 case 3:
169 return (T(5887.) +
170 x * x * (T(-4620.) + x * x * (T(1680.) - x * x * T(320.)))) /
171 T(11520.);
172 case 4:
173 return (T(10543.) +
174 x * (T(17340.) +
175 x * (T(4740.) +
176 x * (T(-6880.) +
177 x * (T(-4080.) + x * (T(960.) + x * T(960.))))))) /
178 T(46080.);
179 case 5:
180 return (T(361.) +
181 x * (T(1416.) +
182 x * (T(2220.) +
183 x * (T(1600.) +
184 x * (T(240.) + x * (T(-384.) - x * T(192.))))))) /
185 T(23040.);
186 case 6:
187 return (T(1.) +
188 x * (T(12.) +
189 x * (T(60.) +
190 x * (T(160.) +
191 x * (T(240.) + x * (T(192.) + x * T(64.))))))) /
192 T(46080.);
193 }
194 }
195
196 DEVICE_THROW(std::runtime_error("Internal interpolation error."));
197 return T{};
198}
199
200/** @brief Derivative of the B-spline. */
201template <int order, typename T = double>
202DEVICE_QUALIFIER auto bspline_d(int i, T x) -> T
203 requires((order > 0) and (order <= 7))
204{
205#if defined(__CUDACC__)
206 // LCOV_EXCL_START
207 [[maybe_unused]] auto constexpr tolerance = T(6) * (T)(FLT_EPSILON);
208 // LCOV_EXCL_STOP
209#else
210 [[maybe_unused]] auto constexpr tolerance =
211 T(6) * std::numeric_limits<T>::epsilon();
212#endif
213 DEVICE_ASSERT(i < order);
214 DEVICE_ASSERT(x >= T(-0.5) or (T(-0.5) - x) < tolerance);
215 DEVICE_ASSERT(x <= T(0.5) or (x - T(0.5)) < tolerance);
216
217 switch (order - 1) {
218 case 0:
219 return 0.;
220 case 1:
221 switch (i) {
222 case 0:
223 return T(-1.);
224 case 1:
225 return T(1.);
226 }
227 case 2:
228 switch (i) {
229 case 0:
230 return x - T(0.5);
231 case 1:
232 return T(-2.) * x;
233 case 2:
234 return x + T(0.5);
235 }
236 case 3:
237 switch (i) {
238 case 0:
239 return (T(-1.) + x * (T(4.) + x * T(-4.))) / T(8.);
240 case 1:
241 return (T(-5.) + x * (T(-4.) + x * T(12.))) / T(8.);
242 case 2:
243 return (T(5.) + x * (T(-4.) + x * T(-12.))) / T(8.);
244 case 3:
245 return (T(1.) + x * (T(4.) + x * T(4.))) / T(8.);
246 }
247 case 4:
248 switch (i) {
249 case 0:
250 return (T(-1.) + x * (T(6.) + x * (T(-12.) + x * T(8.)))) / T(48.);
251 case 1:
252 return (T(-11.) + x * (T(12.) + x * (T(12.) + x * T(-16.)))) / T(24.);
253 case 2:
254 return (x * (T(-5.) + x * x * T(4.))) / T(4.);
255 case 3:
256 return (T(11.) + x * (T(12.) + x * (T(-12.) + x * T(-16.)))) / T(24.);
257 case 4:
258 return (T(1.) + x * (T(6.) + x * (T(12.) + x * T(8.)))) / T(48.);
259 }
260 case 5:
261 switch (i) {
262 case 0:
263 return (T(-1.) +
264 x * (T(8.) + x * (T(-24.) + x * (T(32.) + x * T(-16.))))) /
265 T(384.);
266 case 1:
267 return (T(-75.) +
268 x * (T(168.) + x * (T(-72.) + x * (T(-96.) + x * (T(80.)))))) /
269 T(384.);
270 case 2:
271 return (T(-77.) +
272 x * (T(-88.) + x * (T(168.) + x * (T(32.) + x * T(-80.))))) /
273 T(192.);
274 case 3:
275 return (T(77.) +
276 x * (T(-88.) + x * (T(-168.) + x * (T(32.) + x * T(80.))))) /
277 T(192.);
278 case 4:
279 return (T(75.) +
280 x * (T(168.) + x * (T(72.) + x * (T(-96.) + x * T(-80.))))) /
281 T(384.);
282 case 5:
283 return (T(1.) + x * (T(8.) + x * (T(24.) + x * (T(32.) + x * T(16.))))) /
284 T(384.);
285 }
286 case 6:
287 switch (i) {
288 case 0:
289 return (T(-1.) +
290 x * (T(10.) +
291 x * (T(-40.) + x * (T(80.) + x * (T(-80.) + x * T(32.)))))) /
292 T(3840.);
293 case 1:
294 return (T(-59.) +
295 x * (T(185.) +
296 x * (T(-200.) + x * (T(40.) + x * (T(80.) - x * T(48.)))))) /
297 T(960.);
298 case 2:
299 return (T(-289.) +
300 x * (T(158.) +
301 x * (T(344.) +
302 x * (T(-272.) + x * (T(-80.) + x * T(96.)))))) /
303 T(768.);
304 case 3:
305 return (x * (T(-77.) + x * x * (T(56.) - x * x * T(16.)))) / T(96.);
306 case 4:
307 return (T(289.) + x * (T(158.) + x * (T(-344.) +
308 x * (T(-272.) +
309 x * (T(80.) + x * T(96.)))))) /
310 T(768.);
311 case 5:
312 return (T(59.) +
313 x * (T(185.) +
314 x * (T(200.) + x * (T(40.) + x * (T(-80.) - x * T(48.)))))) /
315 T(960.);
316 case 6:
317 return (T(1.) +
318 x * (T(10.) +
319 x * (T(40.) + x * (T(80.) + x * (T(80.) + x * T(32.)))))) /
320 T(3840.);
321 }
322 }
323
324 DEVICE_THROW(std::runtime_error("Internal interpolation error."));
325 return T{};
326}
327} // namespace Utils
328
329#if defined(__clang__)
330#pragma clang diagnostic pop
331#elif defined(__GNUC__) or defined(__GNUG__)
332#pragma GCC diagnostic pop
333#endif
stream
cudaStream_t stream[1]
CUDA streams for parallel computing on CPU and GPU.
Definition common_cuda.cu:34
device_qualifier.hpp
DEVICE_ASSERT
#define DEVICE_ASSERT(A)
Definition device_qualifier.hpp:30
DEVICE_QUALIFIER
#define DEVICE_QUALIFIER
Definition device_qualifier.hpp:29
DEVICE_THROW
#define DEVICE_THROW(E)
Definition device_qualifier.hpp:28
Utils::bspline_d
DEVICE_QUALIFIER auto bspline_d(int i, T x) -> T requires((order > 0) and(order<=7))
Derivative of the B-spline.
Definition bspline.hpp:202
Utils::sqr
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Definition sqr.hpp:28
Utils::bspline
DEVICE_QUALIFIER auto bspline(int i, T x) -> T requires((order > 0) and(order<=7))
Formula of the B-spline.
Definition bspline.hpp:43