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#include <stdexcept>
26#include <type_traits>
27
28namespace Utils {
29/** @brief Formula of the B-spline. */
30template <int order, typename T>
31DEVICE_QUALIFIER auto
32bspline(int i, T x) -> std::enable_if_t<(order > 0) && (order <= 7), T> {
33 DEVICE_ASSERT(i < order);
34 DEVICE_ASSERT(x >= T(-0.5));
35 DEVICE_ASSERT(x <= T(0.5));
36
37 switch (order) {
38 case 1:
39 return T(1.);
40 case 2:
41 switch (i) {
42 case 0:
43 return T(0.5) - x;
44 case 1:
45 return T(0.5) + x;
46 }
47 case 3:
48 switch (i) {
49 case 0:
50 return T(0.5) * sqr(T(0.5) - x);
51 case 1:
52 return T(0.75) - sqr(x);
53 case 2:
54 return T(0.5) * sqr(T(0.5) + x);
55 }
56 case 4:
57 switch (i) {
58 case 0:
59 return (T(1.) + x * (T(-6.) + x * (T(12.) - x * T(8.)))) / T(48.);
60 case 1:
61 return (T(23.) + x * (T(-30.) + x * (T(-12.) + x * T(24.)))) / T(48.);
62 case 2:
63 return (T(23.) + x * (T(30.) + x * (T(-12.) - x * T(24.)))) / T(48.);
64 case 3:
65 return (T(1.) + x * (T(6.) + x * (T(12.) + x * T(8.)))) / T(48.);
66 }
67 case 5:
68 switch (i) {
69 case 0:
70 return (T(1.) +
71 x * (T(-8.) + x * (T(24.) + x * (T(-32.) + x * T(16.))))) /
72 T(384.);
73 case 1:
74 return (T(19.) +
75 x * (T(-44.) + x * (T(24.) + x * (T(16.) - x * T(16.))))) /
76 T(96.);
77 case 2:
78 return (T(115.) + x * x * (T(-120.) + x * x * T(48.))) / T(192.);
79 case 3:
80 return (T(19.) +
81 x * (T(44.) + x * (T(24.) + x * (T(-16.) - x * T(16.))))) /
82 T(96.);
83 case 4:
84 return (T(1.) + x * (T(8.) + x * (T(24.) + x * (T(32.) + x * T(16.))))) /
85 T(384.);
86 }
87 case 6:
88 switch (i) {
89 case 0:
90 return (T(1.) +
91 x * (T(-10.) +
92 x * (T(40.) + x * (T(-80.) + x * (T(80.) - x * T(32.)))))) /
93 T(3840.);
94 case 1:
95 return (T(237.) +
96 x * (T(-750.) +
97 x * (T(840.) +
98 x * (T(-240.) + x * (T(-240.) + x * T(160.)))))) /
99 T(3840.);
100 case 2:
101 return (T(841.) +
102 x * (T(-770.) +
103 x * (T(-440.) +
104 x * (T(560.) + x * (T(80.) - x * T(160.)))))) /
105 T(1920.);
106 case 3:
107 return (T(841.) +
108 x * (+T(770.) +
109 x * (T(-440.) +
110 x * (T(-560.) + x * (T(80.) + x * T(160.)))))) /
111 T(1920.);
112 case 4:
113 return (T(237.) +
114 x * (T(750.) +
115 x * (T(840.) +
116 x * (T(240.) + x * (T(-240.) - x * T(160.)))))) /
117 T(3840.);
118 case 5:
119 return (T(1.) +
120 x * (T(10.) +
121 x * (T(40.) + x * (T(80.) + x * (T(80.) + x * T(32.)))))) /
122 T(3840.);
123 }
124 case 7:
125 switch (i) {
126 case 0:
127 return (T(1.) +
128 x * (T(-12.) +
129 x * (T(60.) +
130 x * (T(-160.) +
131 x * (T(240.) + x * (T(-192.) + x * T(64.))))))) /
132 T(46080.);
133 case 1:
134 return (T(361.) +
135 x * (T(-1416.) +
136 x * (T(2220.) +
137 x * (T(-1600.) +
138 x * (T(240.) + x * (T(384.) - x * T(192.))))))) /
139 T(23040.);
140 case 2:
141 return (T(10543.) +
142 x * (T(-17340.) +
143 x * (T(4740.) +
144 x * (T(6880.) + x * (T(-4080.) +
145 x * (T(-960.) + x * T(960.))))))) /
146 T(46080.);
147 case 3:
148 return (T(5887.) +
149 x * x * (T(-4620.) + x * x * (T(1680.) - x * x * T(320.)))) /
150 T(11520.);
151 case 4:
152 return (T(10543.) +
153 x * (T(17340.) +
154 x * (T(4740.) +
155 x * (T(-6880.) +
156 x * (T(-4080.) + x * (T(960.) + x * T(960.))))))) /
157 T(46080.);
158 case 5:
159 return (T(361.) +
160 x * (T(1416.) +
161 x * (T(2220.) +
162 x * (T(1600.) +
163 x * (T(240.) + x * (T(-384.) - x * T(192.))))))) /
164 T(23040.);
165 case 6:
166 return (T(1.) +
167 x * (T(12.) +
168 x * (T(60.) +
169 x * (T(160.) +
170 x * (T(240.) + x * (T(192.) + x * T(64.))))))) /
171 T(46080.);
172 }
173 }
174
175 DEVICE_THROW(std::runtime_error("Internal interpolation error."));
176 return T{};
177}
178
179/**
180 * @brief Calculate B-splines.
181 * @param i knot number, using 0-based indexing
182 * @param x position in the range (-0.5, 0.5)
183 * @param k order of the B-spline, using 1-based indexing, i.e. a
184 * B-spline of order @p k is a polynomial of degree <tt>k-1</tt>
185 */
186template <class T> auto bspline(int i, T x, int k) {
187 switch (k) {
188 case 1:
189 return bspline<1>(i, x);
190 case 2:
191 return bspline<2>(i, x);
192 case 3:
193 return bspline<3>(i, x);
194 case 4:
195 return bspline<4>(i, x);
196 case 5:
197 return bspline<5>(i, x);
198 case 6:
199 return bspline<6>(i, x);
200 case 7:
201 return bspline<7>(i, x);
202 }
203
204 return T(0.);
205}
206
207/** @brief Derivative of the B-spline. */
208template <int order, typename T = double>
209DEVICE_QUALIFIER auto
210bspline_d(int i, T x) -> std::enable_if_t<(order > 0) && (order <= 7), T> {
211 DEVICE_ASSERT(i < order);
212 DEVICE_ASSERT(x >= T(-0.5));
213 DEVICE_ASSERT(x <= T(0.5));
214
215 switch (order - 1) {
216 case 0:
217 return 0.;
218 case 1:
219 switch (i) {
220 case 0:
221 return T(-1.);
222 case 1:
223 return T(1.);
224 }
225 case 2:
226 switch (i) {
227 case 0:
228 return x - T(0.5);
229 case 1:
230 return T(-2.) * x;
231 case 2:
232 return x + T(0.5);
233 }
234 case 3:
235 switch (i) {
236 case 0:
237 return (T(-1.) + x * (T(4.) + x * T(-4.))) / T(8.);
238 case 1:
239 return (T(-5.) + x * (T(-4.) + x * T(12.))) / T(8.);
240 case 2:
241 return (T(5.) + x * (T(-4.) + x * T(-12.))) / T(8.);
242 case 3:
243 return (T(1.) + x * (T(4.) + x * T(4.))) / T(8.);
244 }
245 case 4:
246 switch (i) {
247 case 0:
248 return (T(-1.) + x * (T(6.) + x * (T(-12.) + x * T(8.)))) / T(48.);
249 case 1:
250 return (T(-11.) + x * (T(12.) + x * (T(12.) + x * T(-16.)))) / T(24.);
251 case 2:
252 return (x * (T(-5.) + x * x * T(4.))) / T(4.);
253 case 3:
254 return (T(11.) + x * (T(12.) + x * (T(-12.) + x * T(-16.)))) / T(24.);
255 case 4:
256 return (T(1.) + x * (T(6.) + x * (T(12.) + x * T(8.)))) / T(48.);
257 }
258 case 5:
259 switch (i) {
260 case 0:
261 return (T(-1.) +
262 x * (T(8.) + x * (T(-24.) + x * (T(32.) + x * T(-16.))))) /
263 T(384.);
264 case 1:
265 return (T(-75.) +
266 x * (T(168.) + x * (T(-72.) + x * (T(-96.) + x * (T(80.)))))) /
267 T(384.);
268 case 2:
269 return (T(-77.) +
270 x * (T(-88.) + x * (T(168.) + x * (T(32.) + x * T(-80.))))) /
271 T(192.);
272 case 3:
273 return (T(77.) +
274 x * (T(-88.) + x * (T(-168.) + x * (T(32.) + x * T(80.))))) /
275 T(192.);
276 case 4:
277 return (T(75.) +
278 x * (T(168.) + x * (T(72.) + x * (T(-96.) + x * T(-80.))))) /
279 T(384.);
280 case 5:
281 return (T(1.) + x * (T(8.) + x * (T(24.) + x * (T(32.) + x * T(16.))))) /
282 T(384.);
283 }
284 case 6:
285 switch (i) {
286 case 0:
287 return (T(-1.) +
288 x * (T(10.) +
289 x * (T(-40.) + x * (T(80.) + x * (T(-80.) + x * T(32.)))))) /
290 T(3840.);
291 case 1:
292 return (T(-59.) +
293 x * (T(185.) +
294 x * (T(-200.) + x * (T(40.) + x * (T(80.) - x * T(48.)))))) /
295 T(960.);
296 case 2:
297 return (T(-289.) +
298 x * (T(158.) +
299 x * (T(344.) +
300 x * (T(-272.) + x * (T(-80.) + x * T(96.)))))) /
301 T(768.);
302 case 3:
303 return (x * (T(-77.) + x * x * (T(56.) - x * x * T(16.)))) / T(96.);
304 case 4:
305 return (T(289.) + x * (T(158.) + x * (T(-344.) +
306 x * (T(-272.) +
307 x * (T(80.) + x * T(96.)))))) /
308 T(768.);
309 case 5:
310 return (T(59.) +
311 x * (T(185.) +
312 x * (T(200.) + x * (T(40.) + x * (T(-80.) - x * T(48.)))))) /
313 T(960.);
314 case 6:
315 return (T(1.) +
316 x * (T(10.) +
317 x * (T(40.) + x * (T(80.) + x * (T(80.) + x * T(32.)))))) /
318 T(3840.);
319 }
320 }
321
322 DEVICE_THROW(std::runtime_error("Internal interpolation error."));
323 return T{};
324}
325} // namespace Utils
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::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) -> std::enable_if_t<(order > 0) &&(order<=7), T >
Formula of the B-spline.
Definition bspline.hpp:32
Utils::bspline_d
DEVICE_QUALIFIER auto bspline_d(int i, T x) -> std::enable_if_t<(order > 0) &&(order<=7), T >
Derivative of the B-spline.
Definition bspline.hpp:210