source: src/samples/bspline.cpp@ facba0

Last change on this file since facba0 was dfed1c, checked in by Julian Iseringhausen <isering@…>, 14 years ago

Major vmg update.

git-svn-id: ​https://svn.version.fz-juelich.de/scafacos/trunk@1136 5161e1c8-67bf-11de-9fd5-51895aff932f
  • Property mode set to 100644
File size: 7.8 KB
RevLine 
[dfed1c]1/**
2 * @file bspline.cpp
3 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
4 * @date Mon Nov 21 13:27:22 2011
5 *
6 * @brief B-Splines for molecular dynamics.
7 *
8 */
9
10#ifdef HAVE_CONFIG_H
11#include <config.h>
12#endif
13
14#ifndef BSPLINE_DEGREE
15#error BSPLINE_DEGREE not defined.
16#endif
17
18#if BSPLINE_DEGREE < 3 || BSPLINE_DEGREE > 6
19#error Please choose a B-Spline degree between 3 and 6
20#endif
21
22#define _USE_MATH_DEFINES
23#include <cmath>
24
25#include "base/helper.hpp"
26#include "samples/bspline.hpp"
27
28#define POW(x,y) Helper::pow(x,y)
29
30using namespace VMG;
31
32BSpline::BSpline(const vmg_float& width) :
33 R(width)
34{
35 spline_nom.resize(BSPLINE_DEGREE/2+1);
36 spline_denom.resize(BSPLINE_DEGREE/2+1);
37 potential_nom.resize(BSPLINE_DEGREE/2+2);
38 potential_denom.resize(BSPLINE_DEGREE/2+2);
39 intervals.resize(BSPLINE_DEGREE/2+1);
40
41 for (unsigned int i=0; i<intervals.size(); ++i)
42 intervals[i] = R * ( -1.0 + 2.0 / static_cast<vmg_float>(BSPLINE_DEGREE) * (i + BSPLINE_DEGREE/2 + 1));
43
44#if BSPLINE_DEGREE == 3
45
46 spline_nom[0] = Polynomial(2, 81.0*POW(R,2), 0.0, -243.0);
47 spline_nom[1] = Polynomial(2, 243.0*POW(R,2), -486.0*R, 243.0);
48
49 spline_denom[0] = Polynomial(0, 16.0*M_PI*POW(R,5));
50 spline_denom[1] = Polynomial(0, 32.0*M_PI*POW(R,5));
51
52 potential_nom[0] = Polynomial(5,
53 0.0,
54 -195.0*POW(R,4),
55 0.0,
56 270.0*POW(R,2),
57 0.0,
58 -243.0);
59
60 potential_nom[1] = Polynomial(5,
61 2.0*POW(R,5),
62 -405.0*POW(R,4),
63 0.0,
64 810.0*POW(R,2),
65 -810.0*R,
66 243.0);
67
68 potential_nom[2] = Polynomial(0, -1.0);
69
70 potential_denom[0] = Polynomial(0, 320.0*M_PI*POW(R,5));
71 potential_denom[1] = Polynomial(0, 640.0*M_PI*POW(R,5));
72 potential_denom[2] = Polynomial(0, 4.0*M_PI);
73
74 antid = 39.0 / (64.0 * M_PI * R);
75
76#elif BSPLINE_DEGREE == 4
77
78 spline_nom[0] = Polynomial(3,
79 8.0*POW(R,3),
80 0.0,
81 -48.0*R,
82 48.0);
83
84 spline_nom[1] = Polynomial(3,
85 16.0*POW(R,3),
86 -48.0*POW(R,2),
87 48.0*R,
88 -16.0);
89
90 spline_denom[0] = Polynomial(0, M_PI*POW(R,6));
91 spline_denom[1] = Polynomial(0, M_PI*POW(R,6));
92
93 potential_nom[0] = Polynomial(6,
94 0.0,
95 -21.0*POW(R,5),
96 0.0,
97 40.0*POW(R,3),
98 0.0,
99 -72.0*R,
100 48.0);
101
102 potential_nom[1] = Polynomial(6,
103 1.0*POW(R,6),
104 -48.0*POW(R,5),
105 0.0,
106 160.0*POW(R,3),
107 -240.0*POW(R,2),
108 144.0*R,
109 -32.0);
110
111 potential_nom[2] = Polynomial(0, -1.0);
112
113 potential_denom[0] = Polynomial(0, 30.0*M_PI*POW(R,6));
114 potential_denom[1] = Polynomial(0, 60.0*M_PI*POW(R,6));
115 potential_denom[2] = Polynomial(0, 4.0*M_PI);
116
117 antid = 7.0 / (10.0 * M_PI * R);
118
119#elif BSPLINE_DEGREE == 5
120
121 spline_nom[0] = Polynomial(4,
122 2875.0*POW(R, 4),
123 0.0,
124 -18750.0*POW(R,2),
125 0.0,
126 46875.0);
127
128 spline_nom[1] = Polynomial(4,
129 1375.0*POW(R,4),
130 1250.0*POW(R,3),
131 -18750.0*POW(R,2),
132 31250.0*R,
133 -15625.0);
134
135 spline_nom[2] = Polynomial(4,
136 15625.0*POW(R,4),
137 -62500.0*POW(R,3),
138 93750.0*POW(R,2),
139 -62500.0*R,
140 15625.0);
141
142 spline_denom[0] = Polynomial(0, 256.0*M_PI*POW(R,7));
143 spline_denom[1] = Polynomial(0, 128.0*M_PI*POW(R,7));
144 spline_denom[2] = Polynomial(0, 512.0*M_PI*POW(R,7));
145
146 potential_nom[0] = Polynomial(7,
147 0.0,
148 -8393.0*POW(R,6),
149 0.0,
150 20125.0*POW(R,4),
151 0.0,
152 -39375.0*POW(R,2),
153 0.0,
154 46875.0);
155
156 potential_nom[1] = Polynomial(7,
157 -POW(R,7),
158 -20965.0*POW(R,6),
159 0.0,
160 48125.0*POW(R,4),
161 21875.0*POW(R,3),
162 -196875.0*POW(R,2),
163 218750.0*R,
164 78125.0);
165
166 potential_nom[2] = Polynomial(7,
167 874.0*POW(R,7),
168 -21875.0*POW(R,6),
169 0.0,
170 109375.0*POW(R,4),
171 -218750.0*POW(R,3),
172 196875.0*POW(R,2),
173 -87500.0*R,
174 15625.0);
175
176 potential_nom[3] = Polynomial(0, -1.0);
177
178 potential_denom[0] = Polynomial(0, 10752.0*M_PI*POW(R,7));
179 potential_denom[1] = Polynomial(0, 26880.0*M_PI*POW(R,7));
180 potential_denom[2] = Polynomial(0, 21504.0*M_PI*POW(R,7));
181 potential_denom[3] = Polynomial(0, 4.0*M_PI);
182
183 antid = 1199.0 / (1536.0 * M_PI * R);
184
185#elif BSPLINE_DEGREE == 6
186
187 spline_nom[0] = Polynomial(5,
188 297.0*POW(R,5),
189 0.0,
190 -2430.0*POW(R,3),
191 0.0,
192 10935.0*R,
193 -10935.0);
194
195 spline_nom[1] = Polynomial(5,
196 459.0*POW(R,5),
197 2025.0*POW(R,4),
198 -17010.0*POW(R,3),
199 36450.0*POW(R,2),
200 -32805.0*R,
201 10935.0);
202
203 spline_nom[2] = Polynomial(5,
204 2187.0*POW(R,5),
205 -10935.0*POW(R,4),
206 21870.0*POW(R,3),
207 -21870.0*POW(R,2),
208 10935.0*R,
209 -2187.0);
210
211 spline_denom[0] = Polynomial(0, 20.0*M_PI*POW(R,8));
212 spline_denom[1] = Polynomial(0, 40.0*M_PI*POW(R,8));
213 spline_denom[2] = Polynomial(0, 40.0*M_PI*POW(R,8));
214
215 potential_nom[0] = Polynomial(8,
216 0.0,
217 -956.0*POW(R,7),
218 0.0,
219 2772.0*POW(R,5),
220 0.0,
221 -6804.0*POW(R,3),
222 0.0,
223 14580.0*R,
224 -10935.0);
225
226 potential_nom[1] = Polynomial(8,
227 -5.0*POW(R,8),
228 -5676.0*POW(R,7),
229 0.0,
230 12852.0*POW(R,5),
231 28350.0*POW(R,4),
232 -142884.0*POW(R,3),
233 204120.0*POW(R,2),
234 -131220.0*R,
235 32805.0);
236
237 potential_nom[2] = Polynomial(8,
238 169.0*POW(R,8),
239 -2916.0*POW(R,7),
240 0.0,
241 20412.0*POW(R,5),
242 -51030.0*POW(R,4),
243 61236.0*POW(R,3),
244 -40824.0*POW(R,2),
245 14580.0*R,
246 -2187.0);
247
248 potential_nom[3] = Polynomial(0, -1.0);
249
250 potential_denom[0] = Polynomial(0, 1120.0*M_PI*POW(R,8));
251 potential_denom[1] = Polynomial(0, 6720.0*M_PI*POW(R,8));
252 potential_denom[2] = Polynomial(0, 2240.0*M_PI*POW(R,8));
253 potential_denom[3] = Polynomial(0, 4.0*M_PI);
254
255 antid = 239.0 / (280.0 * M_PI * R);
256
257#endif /* BSPLINE_DEGREE */
258}
259
260vmg_float BSpline::EvaluateDerivative(const Vector& v, const int& dv)
261{
262 const vmg_float& v1 = v[dv];
263 const vmg_float L = v.Length();
264 vmg_float result = 0.0;
265
266#if BSPLINE_DEGREE == 3
267
268 if (L < intervals[0])
269 result = v1*(243.0*L*L-135.0*R*R) / (80.0*M_PI*POW(R,5));
270 else if (L < intervals[1])
271 result = (81.0 * v1 * (POW(R,4)-6.0*POW(L,2)*POW(R,2)+8.0*POW(L,3)*R-8.0*POW(L,4))) / (128.0*M_PI*L*POW(R,5));
272 else
273 result = -1.0*v1 / (4*M_PI*POW(L,3));
274
275#elif BSPLINE_DEGREE == 4
276
277 if (L < intervals[0])
278 result = 8.0*v1 * (-5.0*POW(R,3) + 18.0*R*POW(L,2) - 15.0*POW(L,3)) / (15.0*M_PI*POW(R,6));
279 else if (L<intervals[1])
280 result = v1 * (POW(R,6) - 320.0*POW(L,3)*POW(R,3) + 720.0*POW(L,4)*POW(R,2) - 576*POW(L,5)*R + 160.0*POW(L,6)) / (60.0*M_PI*POW(L,3)*POW(R,6));
281 else
282 result = -1.0*v1 / (4*M_PI*POW(L,3));
283
284#elif BSPLINE_DEGREE == 5
285
286 if (L < intervals[0])
287 result = -125.0*v1 * (161.0*POW(R,4) - 630.0*POW(L,2)*POW(R,2) + 1125.0*POW(L,4)) / (5376.0*M_PI*POW(R,7));
288 else if (L < intervals[1])
289 result = v1 * (POW(R,7) + 96250.0*POW(L,3)*POW(R,4) + 65625.0*POW(L,4)*POW(R,3) - 787500.0*POW(L,5)*POW(R,2) + 1093750.0*POW(L,6)*R - 468750.0*POW(L,7)) / (26880.0*M_PI*POW(L,3)*POW(R,7));
290 else if (L < intervals[2])
291 result = v1 * (437.0*POW(R,7) - 109375.0*POW(L,3)*POW(R,4) + 328125.0*POW(L,4)*POW(R,3) - 393750.0*POW(L,5)*POW(R,2) + 218750.0*POW(L,6)*R - 46875.0*POW(L,7)) / (10752.0*M_PI*POW(L,3)*POW(R,7));
292 else
293 result = -1.0*v1 / (4*M_PI*POW(L,3));
294
295#elif BSPLINE_DEGREE == 6
296
297 if (L < intervals[0])
298 result = 9.0*v1 * (616.0*POW(R,5) - 3024.0*POW(L,2)*POW(R,3) + 9720.0*POW(L,4)*R - 8505*POW(L,5)) / (1120.0*M_PI*POW(R,8));
299 else if (L < intervals[1])
300 result = -1.0*v1 * (5.0*POW(R,8) + 25704*POW(L,3)*POW(R,8) + 85050.0*POW(L,4)*POW(R,4) - 571536.0*POW(L,5)*POW(R,3) + 1020600*POW(L,6)*POW(R,2) - 787320.0*POW(L,7)*R + 229635.0*POW(L,8)) / (6720.0*M_PI*POW(L,3)*POW(R,8));
301 else if (L < intervals[2])
302 result = v1 * (169.0*POW(R,8) - 40824.0*POW(L,3)*POW(R,5) + 153090.0*POW(L,4)*POW(R,4) - 244944.0*POW(L,5)*POW(R,3) + 204120.0*POW(L,6)*POW(R,2) - 87480.0*POW(L,7)*R + 15309.0*POW(L,8)) / (2240.0*M_PI*POW(L,3)*POW(R,8));
303 else
304 result = -1.0*v1 / (4*M_PI*POW(L,3));
305
306#endif
307
308 return result;
309}
Note: See TracBrowser for help on using the repository browser.