source: ThirdParty/vmg/src/units/particle/bspline.cpp@ 1ca493a

Candidate_v1.7.0 stable
Last change on this file since 1ca493a was 7faa5c, checked in by Frederik Heber <heber@…>, 9 years ago

Merge commit 'de061d9d851257a04e924d4472df4523d33bb08b' as 'ThirdParty/vmg'
  • Property mode set to 100644
File size: 6.9 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 <libvmg_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 int& near_field_cells, const vmg_float& h) :
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(near_field_cells*h),
57 near_field_cells(near_field_cells)
58{
59 for (unsigned int i=0; i<intervals.size(); ++i)
60 intervals[i] = R * ( -1.0 + 2.0 / static_cast<vmg_float>(BSPLINE_DEGREE) * (i + BSPLINE_DEGREE/2 + 1));
61
62#if BSPLINE_DEGREE == 3
63
64 spline_nom[0] = Polynomial(2, 81.0*POW(R,2), 0.0, -243.0);
65 spline_nom[1] = Polynomial(2, 243.0*POW(R,2), -486.0*R, 243.0);
66
67 spline_denom[0] = Polynomial(0, 16.0 * Math::pi * POW(R,5));
68 spline_denom[1] = Polynomial(0, 32.0 * Math::pi * POW(R,5));
69
70 potential_nom[0] = Polynomial(5,
71 0.0,
72 -195.0*POW(R,4),
73 0.0,
74 270.0*POW(R,2),
75 0.0,
76 -243.0);
77
78 potential_nom[1] = Polynomial(5,
79 2.0 * POW(R,5),
80 -405.0 * POW(R,4),
81 0.0,
82 810.0 * POW(R,2),
83 -810.0 * R,
84 243.0);
85
86 potential_nom[2] = Polynomial(0, -1.0);
87
88 potential_denom[0] = Polynomial(0, 80.0*POW(R,5));
89 potential_denom[1] = Polynomial(0, 160.0*POW(R,5));
90 potential_denom[2] = Polynomial(0, 1.0);
91
92 field_nom[0] = Polynomial(5,
93 20.0 * POW(R,5),
94 0.0,
95 0.0,
96 -135.0 * POW(R,2),
97 0.0,
98 243.0);
99
100 field_nom[1] = Polynomial(5,
101 81.0 * POW(R,5),
102 0.0,
103 0.0,
104 -810.0 * POW(R,2),
105 1215.0 * R,
106 -486.0);
107
108 field_denom[0] = Polynomial(3,
109 0.0,
110 0.0,
111 0.0,
112 20.0 * POW(R,5));
113
114 field_denom[1] = Polynomial(3,
115 0.0,
116 0.0,
117 0.0,
118 80.0 * POW(R,5));
119
120 antid = 39.0 / (16.0 * R);
121
122#elif BSPLINE_DEGREE == 6
123
124 spline_nom[0] = Polynomial(5,
125 297.0 * POW(R,5),
126 0.0,
127 -2430.0 * POW(R,3),
128 0.0,
129 10935.0 * R,
130 -10935.0);
131
132 spline_nom[1] = Polynomial(5,
133 459.0 * POW(R,5),
134 2025.0 * POW(R,4),
135 -17010.0 * POW(R,3),
136 36450.0 * POW(R,2),
137 -32805.0 * R,
138 10935.0);
139
140 spline_nom[2] = Polynomial(5,
141 2187.0 * POW(R,5),
142 -10935.0 * POW(R,4),
143 21870.0 * POW(R,3),
144 -21870.0 * POW(R,2),
145 10935.0 * R,
146 -2187.0);
147
148 spline_denom[0] = Polynomial(0, 20.0 * Math::pi * POW(R,8));
149 spline_denom[1] = Polynomial(0, 40.0 * Math::pi * POW(R,8));
150 spline_denom[2] = Polynomial(0, 40.0 * Math::pi * POW(R,8));
151
152 potential_nom[0] = Polynomial(8,
153 0.0,
154 -956.0 * POW(R,7),
155 0.0,
156 2772.0 * POW(R,5),
157 0.0,
158 -6804.0 * POW(R,3),
159 0.0,
160 14580.0 * R,
161 -10935.0);
162
163 potential_nom[1] = Polynomial(8,
164 -5.0 * POW(R,8),
165 -5676.0 * POW(R,7),
166 0.0,
167 12852.0 * POW(R,5),
168 28350.0 * POW(R,4),
169 -142884.0 * POW(R,3),
170 204120.0 * POW(R,2),
171 -131220.0 * R,
172 32805.0);
173
174 potential_nom[2] = Polynomial(8,
175 169.0 * POW(R,8),
176 -2916.0 * POW(R,7),
177 0.0,
178 20412.0 * POW(R,5),
179 -51030.0 * POW(R,4),
180 61236.0 * POW(R,3),
181 -40824.0 * POW(R,2),
182 14580.0 * R,
183 -2187.0);
184
185 potential_nom[3] = Polynomial(0, -1.0);
186
187 potential_denom[0] = Polynomial(0, 280.0 * POW(R,8));
188 potential_denom[1] = Polynomial(0, 1680.0 * POW(R,8));
189 potential_denom[2] = Polynomial(0, 560.0 * POW(R,8));
190 potential_denom[3] = Polynomial(0, 0.0);
191
192 field_nom[0] = Polynomial(8,
193 280.0 * POW(R,8),
194 0.0,
195 0.0,
196 -5544.0 * POW(R,5),
197 0.0,
198 27216.0 * POW(R,3),
199 0.0,
200 -87480.0 * R,
201 76545.0);
202
203 field_nom[1] = Polynomial(8,
204 1675.0 * POW(R,8),
205 0.0,
206 0.0,
207 -25704.0 * POW(R,5),
208 -85050.0 * POW(R,4),
209 571536.0 * POW(R,3),
210 -1020600.0 * POW(R,2),
211 787320.0 * R,
212 -229635.0);
213
214 field_nom[2] = Polynomial(8,
215 729.0 * POW(R,8),
216 0.0,
217 0.0,
218 -40824.0 * POW(R,5),
219 153090.0 * POW(R,4),
220 -244944.0 * POW(R,3),
221 204120.0 * POW(R,2),
222 -87480.0 * R,
223 15309.0);
224
225 field_denom[0] = Polynomial(3,
226 0.0,
227 0.0,
228 0.0,
229 280.0 * POW(R,8));
230 field_denom[1] = Polynomial(3,
231 0.0,
232 0.0,
233 0.0,
234 1680.0 * POW(R,8));
235
236 field_denom[2] = Polynomial(3,
237 0.0,
238 0.0,
239 0.0,
240 560.0 * POW(R,8));
241
242 antid = 239.0 / (70.0 * R);
243
244#else
245#error B-Spline degree not supported. Choose 3 or 6.
246#endif /* BSPLINE_DEGREE */
247}
Note: See TracBrowser for help on using the repository browser.