source: src/units/particle/bspline.cpp@ f003a9

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

Refactored vmg in order to separate the core library and the particle simulation part properly.

git-svn-id: ​https://svn.version.fz-juelich.de/scafacos/trunk@1798 5161e1c8-67bf-11de-9fd5-51895aff932f
  • Property mode set to 100644
File size: 6.0 KB
Line 
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#include "base/helper.hpp"
23#include "base/math.hpp"
24#include "units/particle/bspline.hpp"
25
26#define POW(x,y) Helper::pow(x,y)
27
28using namespace VMG;
29
30Particle::BSpline::BSpline(const vmg_float& width) :
31 spline_nom((BSPLINE_DEGREE+1)/2),
32 spline_denom((BSPLINE_DEGREE+1)/2),
33 potential_nom((BSPLINE_DEGREE+1)/2+1),
34 potential_denom((BSPLINE_DEGREE+1)/2+1),
35 field_nom((BSPLINE_DEGREE+1)/2),
36 field_denom((BSPLINE_DEGREE+1)/2),
37 intervals((BSPLINE_DEGREE+1)/2),
38 R(width)
39{
40 for (unsigned int i=0; i<intervals.size(); ++i)
41 intervals[i] = R * ( -1.0 + 2.0 / static_cast<vmg_float>(BSPLINE_DEGREE) * (i + BSPLINE_DEGREE/2 + 1));
42
43#if BSPLINE_DEGREE == 3
44
45 spline_nom[0] = Polynomial(2, 81.0*POW(R,2), 0.0, -243.0);
46 spline_nom[1] = Polynomial(2, 243.0*POW(R,2), -486.0*R, 243.0);
47
48 spline_denom[0] = Polynomial(0, 16.0 * Math::pi * POW(R,5));
49 spline_denom[1] = Polynomial(0, 32.0 * Math::pi * POW(R,5));
50
51 potential_nom[0] = Polynomial(5,
52 0.0,
53 -195.0*POW(R,4),
54 0.0,
55 270.0*POW(R,2),
56 0.0,
57 -243.0);
58
59 potential_nom[1] = Polynomial(5,
60 2.0 * POW(R,5),
61 -405.0 * POW(R,4),
62 0.0,
63 810.0 * POW(R,2),
64 -810.0 * R,
65 243.0);
66
67 potential_nom[2] = Polynomial(0, -1.0);
68
69 potential_denom[0] = Polynomial(0, 80.0*POW(R,5));
70 potential_denom[1] = Polynomial(0, 160.0*POW(R,5));
71 potential_denom[2] = Polynomial(0, 1.0);
72
73 field_nom[0] = Polynomial(5,
74 20.0 * POW(R,5),
75 0.0,
76 0.0,
77 -135.0 * POW(R,2),
78 0.0,
79 243.0);
80
81 field_nom[1] = Polynomial(5,
82 81.0 * POW(R,5),
83 0.0,
84 0.0,
85 -810.0 * POW(R,2),
86 1215.0 * R,
87 -486.0);
88
89 field_denom[0] = Polynomial(3,
90 0.0,
91 0.0,
92 0.0,
93 20.0 * POW(R,5));
94
95 field_denom[1] = Polynomial(3,
96 0.0,
97 0.0,
98 0.0,
99 80.0 * POW(R,5));
100
101 antid = 39.0 / (16.0 * R);
102
103#elif BSPLINE_DEGREE == 6
104
105 spline_nom[0] = Polynomial(5,
106 297.0 * POW(R,5),
107 0.0,
108 -2430.0 * POW(R,3),
109 0.0,
110 10935.0 * R,
111 -10935.0);
112
113 spline_nom[1] = Polynomial(5,
114 459.0 * POW(R,5),
115 2025.0 * POW(R,4),
116 -17010.0 * POW(R,3),
117 36450.0 * POW(R,2),
118 -32805.0 * R,
119 10935.0);
120
121 spline_nom[2] = Polynomial(5,
122 2187.0 * POW(R,5),
123 -10935.0 * POW(R,4),
124 21870.0 * POW(R,3),
125 -21870.0 * POW(R,2),
126 10935.0 * R,
127 -2187.0);
128
129 spline_denom[0] = Polynomial(0, 20.0 * Math::pi * POW(R,8));
130 spline_denom[1] = Polynomial(0, 40.0 * Math::pi * POW(R,8));
131 spline_denom[2] = Polynomial(0, 40.0 * Math::pi * POW(R,8));
132
133 potential_nom[0] = Polynomial(8,
134 0.0,
135 -956.0 * POW(R,7),
136 0.0,
137 2772.0 * POW(R,5),
138 0.0,
139 -6804.0 * POW(R,3),
140 0.0,
141 14580.0 * R,
142 -10935.0);
143
144 potential_nom[1] = Polynomial(8,
145 -5.0 * POW(R,8),
146 -5676.0 * POW(R,7),
147 0.0,
148 12852.0 * POW(R,5),
149 28350.0 * POW(R,4),
150 -142884.0 * POW(R,3),
151 204120.0 * POW(R,2),
152 -131220.0 * R,
153 32805.0);
154
155 potential_nom[2] = Polynomial(8,
156 169.0 * POW(R,8),
157 -2916.0 * POW(R,7),
158 0.0,
159 20412.0 * POW(R,5),
160 -51030.0 * POW(R,4),
161 61236.0 * POW(R,3),
162 -40824.0 * POW(R,2),
163 14580.0 * R,
164 -2187.0);
165
166 potential_nom[3] = Polynomial(0, -1.0);
167
168 potential_denom[0] = Polynomial(0, 280.0 * POW(R,8));
169 potential_denom[1] = Polynomial(0, 1680.0 * POW(R,8));
170 potential_denom[2] = Polynomial(0, 560.0 * POW(R,8));
171 potential_denom[3] = Polynomial(0, 0.0);
172
173 field_nom[0] = Polynomial(8,
174 280.0 * POW(R,8),
175 0.0,
176 0.0,
177 -5544.0 * POW(R,5),
178 0.0,
179 27216.0 * POW(R,3),
180 0.0,
181 -87480.0 * R,
182 76545.0);
183
184 field_nom[1] = Polynomial(8,
185 1675.0 * POW(R,8),
186 0.0,
187 0.0,
188 -25704.0 * POW(R,5),
189 -85050.0 * POW(R,4),
190 571536.0 * POW(R,3),
191 -1020600.0 * POW(R,2),
192 787320.0 * R,
193 -229635.0);
194
195 field_nom[2] = Polynomial(8,
196 729.0 * POW(R,8),
197 0.0,
198 0.0,
199 -40824.0 * POW(R,5),
200 153090.0 * POW(R,4),
201 -244944.0 * POW(R,3),
202 204120.0 * POW(R,2),
203 -87480.0 * R,
204 15309.0);
205
206 field_denom[0] = Polynomial(3,
207 0.0,
208 0.0,
209 0.0,
210 280.0 * POW(R,8));
211 field_denom[1] = Polynomial(3,
212 0.0,
213 0.0,
214 0.0,
215 1680.0 * POW(R,8));
216
217 field_denom[2] = Polynomial(3,
218 0.0,
219 0.0,
220 0.0,
221 560.0 * POW(R,8));
222
223 antid = 239.0 / (70.0 * R);
224
225#else
226#error B-Spline degree not supported. Choose 3 or 6.
227#endif /* BSPLINE_DEGREE */
228}
Note: See TracBrowser for help on using the repository browser.