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

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

vmg: Added license files and headers (GPLv3).

git-svn-id: ​https://svn.version.fz-juelich.de/scafacos/trunk@1812 5161e1c8-67bf-11de-9fd5-51895aff932f
  • Property mode set to 100644
File size: 6.8 KB
Line 
1/*
2 * vmg - a versatile multigrid solver
3 * Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
4 *
5 * vmg is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * vmg is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19/**
20 * @file bspline.cpp
21 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
22 * @date Mon Nov 21 13:27:22 2011
23 *
24 * @brief B-Splines for molecular dynamics.
25 *
26 */
27
28#ifdef HAVE_CONFIG_H
29#include <config.h>
30#endif
31
32#ifndef BSPLINE_DEGREE
33#error BSPLINE_DEGREE not defined.
34#endif
35
36#if BSPLINE_DEGREE < 3 || BSPLINE_DEGREE > 6
37#error Please choose a B-Spline degree between 3 and 6
38#endif
39
40#include "base/helper.hpp"
41#include "base/math.hpp"
42#include "units/particle/bspline.hpp"
43
44#define POW(x,y) Helper::pow(x,y)
45
46using namespace VMG;
47
48Particle::BSpline::BSpline(const vmg_float& width) :
49 spline_nom((BSPLINE_DEGREE+1)/2),
50 spline_denom((BSPLINE_DEGREE+1)/2),
51 potential_nom((BSPLINE_DEGREE+1)/2+1),
52 potential_denom((BSPLINE_DEGREE+1)/2+1),
53 field_nom((BSPLINE_DEGREE+1)/2),
54 field_denom((BSPLINE_DEGREE+1)/2),
55 intervals((BSPLINE_DEGREE+1)/2),
56 R(width)
57{
58 for (unsigned int i=0; i<intervals.size(); ++i)
59 intervals[i] = R * ( -1.0 + 2.0 / static_cast<vmg_float>(BSPLINE_DEGREE) * (i + BSPLINE_DEGREE/2 + 1));
60
61#if BSPLINE_DEGREE == 3
62
63 spline_nom[0] = Polynomial(2, 81.0*POW(R,2), 0.0, -243.0);
64 spline_nom[1] = Polynomial(2, 243.0*POW(R,2), -486.0*R, 243.0);
65
66 spline_denom[0] = Polynomial(0, 16.0 * Math::pi * POW(R,5));
67 spline_denom[1] = Polynomial(0, 32.0 * Math::pi * POW(R,5));
68
69 potential_nom[0] = Polynomial(5,
70 0.0,
71 -195.0*POW(R,4),
72 0.0,
73 270.0*POW(R,2),
74 0.0,
75 -243.0);
76
77 potential_nom[1] = Polynomial(5,
78 2.0 * POW(R,5),
79 -405.0 * POW(R,4),
80 0.0,
81 810.0 * POW(R,2),
82 -810.0 * R,
83 243.0);
84
85 potential_nom[2] = Polynomial(0, -1.0);
86
87 potential_denom[0] = Polynomial(0, 80.0*POW(R,5));
88 potential_denom[1] = Polynomial(0, 160.0*POW(R,5));
89 potential_denom[2] = Polynomial(0, 1.0);
90
91 field_nom[0] = Polynomial(5,
92 20.0 * POW(R,5),
93 0.0,
94 0.0,
95 -135.0 * POW(R,2),
96 0.0,
97 243.0);
98
99 field_nom[1] = Polynomial(5,
100 81.0 * POW(R,5),
101 0.0,
102 0.0,
103 -810.0 * POW(R,2),
104 1215.0 * R,
105 -486.0);
106
107 field_denom[0] = Polynomial(3,
108 0.0,
109 0.0,
110 0.0,
111 20.0 * POW(R,5));
112
113 field_denom[1] = Polynomial(3,
114 0.0,
115 0.0,
116 0.0,
117 80.0 * POW(R,5));
118
119 antid = 39.0 / (16.0 * R);
120
121#elif BSPLINE_DEGREE == 6
122
123 spline_nom[0] = Polynomial(5,
124 297.0 * POW(R,5),
125 0.0,
126 -2430.0 * POW(R,3),
127 0.0,
128 10935.0 * R,
129 -10935.0);
130
131 spline_nom[1] = Polynomial(5,
132 459.0 * POW(R,5),
133 2025.0 * POW(R,4),
134 -17010.0 * POW(R,3),
135 36450.0 * POW(R,2),
136 -32805.0 * R,
137 10935.0);
138
139 spline_nom[2] = Polynomial(5,
140 2187.0 * POW(R,5),
141 -10935.0 * POW(R,4),
142 21870.0 * POW(R,3),
143 -21870.0 * POW(R,2),
144 10935.0 * R,
145 -2187.0);
146
147 spline_denom[0] = Polynomial(0, 20.0 * Math::pi * POW(R,8));
148 spline_denom[1] = Polynomial(0, 40.0 * Math::pi * POW(R,8));
149 spline_denom[2] = Polynomial(0, 40.0 * Math::pi * POW(R,8));
150
151 potential_nom[0] = Polynomial(8,
152 0.0,
153 -956.0 * POW(R,7),
154 0.0,
155 2772.0 * POW(R,5),
156 0.0,
157 -6804.0 * POW(R,3),
158 0.0,
159 14580.0 * R,
160 -10935.0);
161
162 potential_nom[1] = Polynomial(8,
163 -5.0 * POW(R,8),
164 -5676.0 * POW(R,7),
165 0.0,
166 12852.0 * POW(R,5),
167 28350.0 * POW(R,4),
168 -142884.0 * POW(R,3),
169 204120.0 * POW(R,2),
170 -131220.0 * R,
171 32805.0);
172
173 potential_nom[2] = Polynomial(8,
174 169.0 * POW(R,8),
175 -2916.0 * POW(R,7),
176 0.0,
177 20412.0 * POW(R,5),
178 -51030.0 * POW(R,4),
179 61236.0 * POW(R,3),
180 -40824.0 * POW(R,2),
181 14580.0 * R,
182 -2187.0);
183
184 potential_nom[3] = Polynomial(0, -1.0);
185
186 potential_denom[0] = Polynomial(0, 280.0 * POW(R,8));
187 potential_denom[1] = Polynomial(0, 1680.0 * POW(R,8));
188 potential_denom[2] = Polynomial(0, 560.0 * POW(R,8));
189 potential_denom[3] = Polynomial(0, 0.0);
190
191 field_nom[0] = Polynomial(8,
192 280.0 * POW(R,8),
193 0.0,
194 0.0,
195 -5544.0 * POW(R,5),
196 0.0,
197 27216.0 * POW(R,3),
198 0.0,
199 -87480.0 * R,
200 76545.0);
201
202 field_nom[1] = Polynomial(8,
203 1675.0 * POW(R,8),
204 0.0,
205 0.0,
206 -25704.0 * POW(R,5),
207 -85050.0 * POW(R,4),
208 571536.0 * POW(R,3),
209 -1020600.0 * POW(R,2),
210 787320.0 * R,
211 -229635.0);
212
213 field_nom[2] = Polynomial(8,
214 729.0 * POW(R,8),
215 0.0,
216 0.0,
217 -40824.0 * POW(R,5),
218 153090.0 * POW(R,4),
219 -244944.0 * POW(R,3),
220 204120.0 * POW(R,2),
221 -87480.0 * R,
222 15309.0);
223
224 field_denom[0] = Polynomial(3,
225 0.0,
226 0.0,
227 0.0,
228 280.0 * POW(R,8));
229 field_denom[1] = Polynomial(3,
230 0.0,
231 0.0,
232 0.0,
233 1680.0 * POW(R,8));
234
235 field_denom[2] = Polynomial(3,
236 0.0,
237 0.0,
238 0.0,
239 560.0 * POW(R,8));
240
241 antid = 239.0 / (70.0 * R);
242
243#else
244#error B-Spline degree not supported. Choose 3 or 6.
245#endif /* BSPLINE_DEGREE */
246}
Note: See TracBrowser for help on using the repository browser.