OPAL (Object Oriented Parallel Accelerator Library) 2024.2
OPAL
FM3DMagnetoStaticExtended.cpp
Go to the documentation of this file.
1#include "Fields/FM3DMagnetoStaticExtended.h"
2#include "Fields/Fieldmap.hpp"
3#include "Physics/Units.h"
4#include "Utilities/GeneralClassicException.h"
5#include "Utilities/Util.h"
6
7#include <cmath>
8#include <fstream>
9#include <ios>
10#include <algorithm>
11
12extern Inform *gmsg;
13
14_FM3DMagnetoStaticExtended::_FM3DMagnetoStaticExtended(const std::string& filename):
15 _Fieldmap(filename),
16 FieldstrengthBx_m(nullptr),
17 FieldstrengthBy_m(nullptr),
18 FieldstrengthBz_m(nullptr) {
19 std::ifstream file;
20 std::string tmpString;
21 double tmpDouble;
22
23 Type = T2DMagnetoStatic;
24
25 // open field map, parse it and disable element on error
26 file.open(Filename_m.c_str());
27 if(file.good()) {
28 bool parsing_passed = true;
29 try {
30 parsing_passed = interpretLine<std::string>(file, tmpString);
31 } catch (GeneralClassicException &e) {
32 parsing_passed = interpretLine<std::string, std::string>(file, tmpString, tmpString);
33
34 tmpString = Util::toUpper(tmpString);
35 if (tmpString != "TRUE" &&
36 tmpString != "FALSE")
37 throw GeneralClassicException("_FM3DMagnetoStaticExtended::_FM3DMagnetoStaticExtended",
38 "The second string on the first line of 3D field "
39 "maps has to be either TRUE or FALSE");
40
41 normalize_m = (tmpString == "TRUE");
42 }
43 parsing_passed = (parsing_passed &&
44 interpretLine<double, double, unsigned int>(file, xbegin_m, xend_m, num_gridpx_m));
45 parsing_passed = (parsing_passed &&
46 interpretLine<double, double, unsigned int>(file, ybegin_m, yend_m, num_gridpy_m));
47 parsing_passed = (parsing_passed &&
48 interpretLine<double, double, unsigned int>(file, zbegin_m, zend_m, num_gridpz_m));
49
50 for(unsigned long i = 0; (i < (num_gridpz_m + 1) * (num_gridpx_m + 1)) && parsing_passed; ++ i) {
51 parsing_passed = parsing_passed && interpretLine<double>(file, tmpDouble);
52 }
53
54 parsing_passed = parsing_passed &&
55 interpreteEOF(file);
56
57 file.close();
58 lines_read_m = 0;
59
60 if(!parsing_passed) {
61 disableFieldmapWarning();
62 zend_m = zbegin_m - 1e-3;
63 throw GeneralClassicException("_FM3DMagnetoStaticExtended::_FM3DMagnetoStaticExtended(std::string)",
64 "Format of fieldmap '" + Filename_m + "' didn't pass basic test");
65 } else {
66 // conversion from cm to m
67 xbegin_m *= Units::cm2m;
68 xend_m *= Units::cm2m;
69 yend_m = std::max(std::abs(ybegin_m), yend_m) * Units::cm2m;
70 ybegin_m = 0.0;
71 zbegin_m *= Units::cm2m;
72 zend_m *= Units::cm2m;
73 length_m = zend_m - zbegin_m;
74
75 hx_m = (xend_m - xbegin_m) / num_gridpx_m;
76 hy_m = (yend_m - ybegin_m) / num_gridpy_m;
77 hz_m = (zend_m - zbegin_m) / num_gridpz_m;
78
79 // num spacings -> num grid points
80 num_gridpx_m++;
81 num_gridpy_m++;
82 num_gridpz_m++;
83 }
84 } else {
85 throw GeneralClassicException("_FM3DMagnetoStaticExtended::_FM3DMagnetoStaticExtended(std::string)",
86 "Couldn't read fieldmap '" + Filename_m + "'");
87 }
88}
89
90_FM3DMagnetoStaticExtended::~_FM3DMagnetoStaticExtended() {
91 freeMap();
92}
93
94FM3DMagnetoStaticExtended _FM3DMagnetoStaticExtended::create(const std::string& filename)
95{
96 return FM3DMagnetoStaticExtended(new _FM3DMagnetoStaticExtended(filename));
97}
98
99void _FM3DMagnetoStaticExtended::readMap() {
100 if(FieldstrengthBz_m == nullptr) {
101 // declare variables and allocate memory
102 std::ifstream in;
103 std::string tmpString;
104 const size_t totalSize = num_gridpx_m * num_gridpy_m * num_gridpz_m;
105 FieldstrengthBx_m = new double[totalSize];
106 FieldstrengthBy_m = new double[totalSize];
107 FieldstrengthBz_m = new double[totalSize];
108
109 // read in and parse field map
110 in.open(Filename_m.c_str());
111 getLine(in, tmpString);
112 getLine(in, tmpString);
113 getLine(in, tmpString);
114 getLine(in, tmpString);
115
116 for(unsigned int i = 0; i < num_gridpx_m; i++) {
117 for(unsigned int k = 0; k < num_gridpz_m; k++) {
118 unsigned long index = getIndex(i,0,k);
119 interpretLine<double>(in, FieldstrengthBy_m[index]);
120 }
121 }
122 in.close();
123
124 std::fill(FieldstrengthBx_m, FieldstrengthBx_m + totalSize, 0.0);
125 std::fill(FieldstrengthBz_m, FieldstrengthBz_m + totalSize, 0.0);
126
127 if (normalize_m) {
128 double Bymax = 0.0;
129
130 // find maximum field
131 unsigned int centerX = static_cast<unsigned int>(std::round(-xbegin_m / hx_m));
132 for(unsigned int k = 0; k < num_gridpz_m; ++ k) {
133 double By = FieldstrengthBy_m[getIndex(centerX, 0, k)];
134 if(std::abs(By) > std::abs(Bymax)) {
135 Bymax = By;
136 }
137 }
138
139 // normalize field
140 for(unsigned int i = 0; i < num_gridpx_m; i ++) {
141 for(unsigned int k = 0; k < num_gridpz_m; k ++) {
142 unsigned long index = getIndex(i,0,k);
143 FieldstrengthBy_m[index] /= Bymax;
144 }
145 }
146 }
147
148 smoothData(FieldstrengthBy_m, 0);
149 for (unsigned int j = 1; j < num_gridpy_m; j ++) {
150 integrateBx(j);
151 integrateBz(j);
152 integrateBy(j);
153
154 smoothData(FieldstrengthBx_m, j);
155 smoothData(FieldstrengthBz_m, j);
156 smoothData(FieldstrengthBy_m, j);
157 }
158
159 INFOMSG(level3 << typeset_msg("read in fieldmap '" + Filename_m + "'", "info") << endl);
160 }
161}
162
163void _FM3DMagnetoStaticExtended::integrateBx(unsigned j) {
164 if (j == 1) {
165 {
166 unsigned int i = 0;
167 // treat cells at i = 0, j = 1, k = 0:num_gridpz_m - 1;
168 for (unsigned int k = 0; k < num_gridpz_m; k ++) {
169 unsigned long index = getIndex(i,j,k);
170 FieldstrengthBx_m[index] = 0.5 * hy_m / hx_m * (-3 * FieldstrengthBy_m[getIndex(i, j - 1, k)] +
171 4 * FieldstrengthBy_m[getIndex(i + 1, j - 1, k)] -
172 FieldstrengthBy_m[getIndex(i + 2, j - 1, k)]);
173 }
174 }
175 // treat cells at i = 1:num_gridpx_m - 2, j = 1, k = 0:num_gridpz_m - 1;
176 for(unsigned int i = 1; i < num_gridpx_m - 1; i ++) {
177 for (unsigned int k = 0; k < num_gridpz_m; k ++) {
178 unsigned long index = getIndex(i,j,k);
179 FieldstrengthBx_m[index] = 0.5 * hy_m / hx_m * (FieldstrengthBy_m[getIndex(i + 1, j - 1, k)] -
180 FieldstrengthBy_m[getIndex(i - 1, j - 1, k)]);
181 }
182 }
183 {
184 unsigned int i = num_gridpx_m - 1;
185 // treat cells at i = num_gridpx_m - 1, j = 1, k = 0:num_gridpz_m - 1;
186 for (unsigned int k = 0; k < num_gridpz_m; k ++) {
187 unsigned long index = getIndex(i,j,k);
188 FieldstrengthBx_m[index] = 0.5 * hy_m / hx_m * (FieldstrengthBy_m[getIndex(i - 2, j - 1, k)] -
189 4 * FieldstrengthBy_m[getIndex(i - 1, j - 1, k)] +
190 3 * FieldstrengthBy_m[getIndex(i, j - 1, k)]);
191 }
192 }
193
194 } else {
195 { // treat cells at i = 0, j > 1, k = 0:num_gridpz_m - 1;
196 unsigned int i = 0;
197 for (unsigned int k = 0; k < num_gridpz_m; k ++) {
198 unsigned long index = getIndex(i,j,k);
199 FieldstrengthBx_m[index] = (4 * FieldstrengthBx_m[getIndex(i, j - 1, k)] -
200 FieldstrengthBx_m[getIndex(i, j - 2, k)] +
201 hy_m / hx_m * (-3 * FieldstrengthBy_m[getIndex(i, j - 1, k)] +
202 4 * FieldstrengthBy_m[getIndex(i + 1, j - 1, k)] -
203 FieldstrengthBy_m[getIndex(i + 2, j - 1, k)])) / 3;
204 }
205 }
206 // treat cells at i = 1:num_gridpx_m - 2, j > 1, k = 0:num_gridpz_m - 1;
207 for(unsigned int i = 1; i < num_gridpx_m - 1; i ++) {
208 for (unsigned int k = 0; k < num_gridpz_m; k ++) {
209 unsigned long index = getIndex(i,j,k);
210 FieldstrengthBx_m[index] = (4 * FieldstrengthBx_m[getIndex(i, j - 1, k)] -
211 FieldstrengthBx_m[getIndex(i, j - 2, k)] +
212 hy_m / hx_m * (FieldstrengthBy_m[getIndex(i + 1, j - 1, k)] -
213 FieldstrengthBy_m[getIndex(i - 1, j - 1, k)])) / 3;
214 }
215 }
216 { // treat cells at i = num_gridpx_m - 1, j > 1, k = 0:num_gridpz_m - 1;
217 unsigned int i = num_gridpx_m - 1;
218 for (unsigned int k = 1; k < num_gridpz_m - 1; k ++) {
219 unsigned long index = getIndex(i,j,k);
220 FieldstrengthBx_m[index] = (4 * FieldstrengthBx_m[getIndex(i, j - 1, k)] -
221 FieldstrengthBx_m[getIndex(i, j - 2, k)] +
222 hy_m / hx_m * (FieldstrengthBy_m[getIndex(i - 2, j - 1, k)] -
223 4 * FieldstrengthBy_m[getIndex(i - 1, j - 1, k)] +
224 3 * FieldstrengthBy_m[getIndex(i, j - 1, k)])) / 3;
225 }
226 }
227 }
228}
229
230void _FM3DMagnetoStaticExtended::integrateBz(unsigned j) {
231 if (j == 1) {
232 // treat cells at i = 0:num_gridpx_m - 1, j = 1, k = 1:num_gridpz_m - 2;
233 for(unsigned int i = 0; i < num_gridpx_m; i ++) {
234 for (unsigned int k = 1; k < num_gridpz_m - 1; k ++) {
235 unsigned long index = getIndex(i,j,k);
236 FieldstrengthBz_m[index] = 0.5 * hy_m / hz_m * (FieldstrengthBy_m[getIndex(i, j - 1, k + 1)] -
237 FieldstrengthBy_m[getIndex(i, j - 1, k - 1)]);
238 }
239
240 { // treat cells at i = 0:num_gridpx_m - 1, j = 1, k = 0;
241 unsigned int k = 0;
242 unsigned int index = getIndex(i,j,k);
243 FieldstrengthBz_m[index] = 0.5 * hy_m / hz_m * (-3 * FieldstrengthBy_m[getIndex(i, j - 1, k)] +
244 4 * FieldstrengthBy_m[getIndex(i, j - 1, k + 1)] -
245 FieldstrengthBy_m[getIndex(i, j - 1, k + 2)]);
246 }
247 { // treat cells at i = 0:num_gridpx_m - 1, j = 1, k = num_gridpz_m - 1;
248 unsigned int k = num_gridpz_m - 1;
249 unsigned int index = getIndex(i,j,k);
250 FieldstrengthBz_m[index] = 0.5 * hy_m / hz_m * (FieldstrengthBy_m[getIndex(i, j - 1, k - 2)] -
251 4 * FieldstrengthBy_m[getIndex(i, j - 1, k - 1)] +
252 3 * FieldstrengthBy_m[getIndex(i,j - 1, k)]);
253 }
254 }
255 } else {
256 // treat cells at i = 0:num_gridpx_m - 1, j > 1, k = 1:num_gridpz_m - 2;
257 for(unsigned int i = 0; i < num_gridpx_m; i ++) {
258 for (unsigned int k = 1; k < num_gridpz_m - 1; k ++) {
259 unsigned long index = getIndex(i,j,k);
260 FieldstrengthBz_m[index] = (4 * FieldstrengthBz_m[getIndex(i, j - 1, k)] -
261 FieldstrengthBz_m[getIndex(i, j - 2, k)] +
262 hy_m / hz_m * (FieldstrengthBy_m[getIndex(i, j - 1, k + 1)] -
263 FieldstrengthBy_m[getIndex(i, j - 1, k - 1)])) / 3;
264 }
265
266 { // treat cells at i = 0:num_gridpx_m - 1, j > 1, k = 0;
267 unsigned int k = 0;
268 unsigned int index = getIndex(i,j,k);
269 FieldstrengthBz_m[index] = (4 * FieldstrengthBz_m[getIndex(i, j - 1, k)] -
270 FieldstrengthBz_m[getIndex(i, j - 2, k)] +
271 hy_m / hz_m * (-3 * FieldstrengthBy_m[getIndex(i, j - 1, k)] +
272 4 * FieldstrengthBy_m[getIndex(i, j - 1, k + 1)] -
273 FieldstrengthBy_m[getIndex(i, j - 1, k + 2)])) / 3;
274 }
275 { // treat cells at i = 0:num_gridpx_m - 1, j > 1, k = num_gridpz_m - 1;
276 unsigned int k = num_gridpz_m - 1;
277 unsigned int index = getIndex(i,j,k);
278 FieldstrengthBz_m[index] = (4 * FieldstrengthBz_m[getIndex(i, j - 1, k)] -
279 FieldstrengthBz_m[getIndex(i, j - 2, k)] +
280 hy_m / hz_m * (FieldstrengthBy_m[getIndex(i, j - 1, k - 2)] -
281 4 * FieldstrengthBy_m[getIndex(i, j - 1, k - 1)] +
282 3 * FieldstrengthBy_m[getIndex(i,j - 1, k)])) / 3;
283 }
284 }
285 }
286}
287
288void _FM3DMagnetoStaticExtended::integrateBy(unsigned j) {
289 if (j == 1) {
290 for (unsigned int i = 1; i < num_gridpx_m - 1; i ++) {
291 { // treat cells at i = 1:num_gridpx_m - 2, j = 1, k = 0
292 unsigned int k = 0;
293 unsigned long index = getIndex(i, j, k);
294 FieldstrengthBy_m[index] = (FieldstrengthBy_m[getIndex(i, j - 1, k)] +
295 hy_m * ((FieldstrengthBx_m[getIndex(i + 1, j, k)] -
296 FieldstrengthBx_m[getIndex(i - 1, j, k)]) / hx_m -
297 (-3 * FieldstrengthBz_m[getIndex(i, j, k)] +
298 4 * FieldstrengthBz_m[getIndex(i, j, k + 1)] -
299 FieldstrengthBz_m[getIndex(i, j, k + 2)]) / hz_m));
300 }
301 // treat cells at i = 1:num_gridpx_m - 2, j = 1, k = 1:num_gridpz_m - 2
302 for (unsigned int k = 1; k < num_gridpz_m - 1; k ++) {
303 unsigned long index = getIndex(i,j,k);
304 FieldstrengthBy_m[index] = (FieldstrengthBy_m[getIndex(i, j - 1, k)] +
305 hy_m * ((FieldstrengthBx_m[getIndex(i + 1, j, k)] -
306 FieldstrengthBx_m[getIndex(i - 1, j, k)]) / hx_m -
307 (FieldstrengthBz_m[getIndex(i, j, k + 1)] -
308 FieldstrengthBz_m[getIndex(i, j, k - 1)]) / hz_m));
309 }
310 { // treat cells at i = 1:num_gridpx_m - 2, j = 1, k = num_gridpz_m - 1
311 unsigned int k = num_gridpz_m - 1;
312 unsigned long index = getIndex(i, j, k);
313 FieldstrengthBy_m[index] = (FieldstrengthBy_m[getIndex(i, j - 1, k)] +
314 hy_m * ((FieldstrengthBx_m[getIndex(i + 1, j, k)] -
315 FieldstrengthBx_m[getIndex(i - 1, j, k)]) / hx_m -
316 (FieldstrengthBz_m[getIndex(i, j, k - 2)] -
317 4 * FieldstrengthBz_m[getIndex(i, j, k - 1)] +
318 3 * FieldstrengthBz_m[getIndex(i, j, k)]) / hz_m));
319 }
320 }
321 {
322 unsigned int i = 0;
323 { // treat cells at i = 0, j = 1, k = 0
324 unsigned int k = 0;
325 unsigned long index = getIndex(i, j, k);
326 FieldstrengthBy_m[index] = (FieldstrengthBy_m[getIndex(i, j - 1, k)] +
327 hy_m * ((-3 * FieldstrengthBx_m[getIndex(i, j, k)] +
328 4 * FieldstrengthBx_m[getIndex(i + 1, j, k)] -
329 FieldstrengthBx_m[getIndex(i + 2, j, k)]) / hx_m -
330 (-3 * FieldstrengthBz_m[getIndex(i, j, k)] +
331 4 * FieldstrengthBz_m[getIndex(i, j, k + 1)] -
332 FieldstrengthBz_m[getIndex(i, j, k + 2)]) / hz_m));
333 }
334 // treat cells at i = 0, j = 1, k = 1:num_gridpz_m - 2
335 for (unsigned int k = 1; k < num_gridpz_m - 1; k ++) {
336 unsigned long index = getIndex(i,j,k);
337 FieldstrengthBy_m[index] = (FieldstrengthBy_m[getIndex(i, j - 1, k)] +
338 hy_m * ((-3 * FieldstrengthBx_m[getIndex(i, j, k)] +
339 4 * FieldstrengthBx_m[getIndex(i + 1, j, k)] -
340 FieldstrengthBx_m[getIndex(i + 2, j, k)]) / hx_m -
341 (FieldstrengthBz_m[getIndex(i, j, k + 1)] -
342 FieldstrengthBz_m[getIndex(i, j, k - 1)]) / hz_m));
343 }
344 { // treat cells at i = 0, j = 1, k = num_gridpz_m - 1
345 unsigned int k = num_gridpz_m - 1;
346 unsigned long index = getIndex(i, j, k);
347 FieldstrengthBy_m[index] = (FieldstrengthBy_m[getIndex(i, j - 1, k)] +
348 hy_m * ((-3 * FieldstrengthBx_m[getIndex(i, j, k)] +
349 4 * FieldstrengthBx_m[getIndex(i + 1, j, k)] -
350 FieldstrengthBx_m[getIndex(i + 2, j, k)]) / hx_m -
351 (FieldstrengthBz_m[getIndex(i, j, k - 2)] -
352 4 * FieldstrengthBz_m[getIndex(i, j, k - 1)] +
353 3 * FieldstrengthBz_m[getIndex(i, j, k)]) / hz_m));
354 }
355 }
356 {
357 unsigned int i = num_gridpx_m - 1;
358 { // treat cells at i = 0, j = 1, k = 0
359 unsigned int k = 0;
360 unsigned long index = getIndex(i, j, k);
361 FieldstrengthBy_m[index] = (FieldstrengthBy_m[getIndex(i, j - 1, k)] +
362 hy_m * ((FieldstrengthBx_m[getIndex(i - 2, j, k)] -
363 4 * FieldstrengthBx_m[getIndex(i - 1, j, k)] +
364 3 * FieldstrengthBx_m[getIndex(i, j, k)]) / hx_m -
365 (-3 * FieldstrengthBz_m[getIndex(i, j, k)] +
366 4 * FieldstrengthBz_m[getIndex(i, j, k + 1)] -
367 FieldstrengthBz_m[getIndex(i, j, k + 2)]) / hz_m));
368 }
369 // treat cells at i = 0, j = 1, k = 1:num_gridpz_m - 2
370 for (unsigned int k = 1; k < num_gridpz_m - 1; k ++) {
371 unsigned long index = getIndex(i,j,k);
372 FieldstrengthBy_m[index] = (FieldstrengthBy_m[getIndex(i, j - 1, k)] +
373 hy_m * ((FieldstrengthBx_m[getIndex(i - 2, j, k)] -
374 4 * FieldstrengthBx_m[getIndex(i - 1, j, k)] +
375 3 * FieldstrengthBx_m[getIndex(i, j, k)]) / hx_m -
376 (FieldstrengthBz_m[getIndex(i, j, k + 1)] -
377 FieldstrengthBz_m[getIndex(i, j, k - 1)]) / hz_m));
378 }
379 { // treat cells at i = 0, j = 1, k = num_gridpz_m - 1
380 unsigned int k = num_gridpz_m - 1;
381 unsigned long index = getIndex(i, j, k);
382 FieldstrengthBy_m[index] = (FieldstrengthBy_m[getIndex(i, j - 1, k)] +
383 hy_m * ((FieldstrengthBx_m[getIndex(i - 2, j, k)] -
384 4 * FieldstrengthBx_m[getIndex(i - 1, j, k)] +
385 3 * FieldstrengthBx_m[getIndex(i, j, k)]) / hx_m -
386 (FieldstrengthBz_m[getIndex(i, j, k - 2)] -
387 4 * FieldstrengthBz_m[getIndex(i, j, k - 1)] +
388 3 * FieldstrengthBz_m[getIndex(i, j, k)]) / hz_m));
389 }
390 }
391 } else {
392 for (unsigned int i = 1; i < num_gridpx_m - 1; i ++) {
393 { // treat cells at i = 1:num_gridpx_m - 2, j > 1, k = 0
394 unsigned int k = 0;
395 unsigned long index = getIndex(i, j, k);
396 FieldstrengthBy_m[index] = (4 * FieldstrengthBy_m[getIndex(i, j - 1, k)] -
397 FieldstrengthBy_m[getIndex(i, j - 2, k)] +
398 hy_m * ((FieldstrengthBx_m[getIndex(i + 1, j, k)] -
399 FieldstrengthBx_m[getIndex(i - 1, j, k)]) / hx_m -
400 (-3 * FieldstrengthBz_m[getIndex(i, j, k)] +
401 4 * FieldstrengthBz_m[getIndex(i, j, k + 1)] -
402 FieldstrengthBz_m[getIndex(i, j, k + 2)]) / hz_m)) / 3;
403 }
404 // treat cells at i = 1:num_gridpx_m - 2, j > 1, k = 1:num_gridpz_m - 2
405 for (unsigned int k = 1; k < num_gridpz_m - 1; k ++) {
406 unsigned long index = getIndex(i,j,k);
407 FieldstrengthBy_m[index] = (4 * FieldstrengthBy_m[getIndex(i, j - 1, k)] -
408 FieldstrengthBy_m[getIndex(i, j - 2, k)] +
409 hy_m * ((FieldstrengthBx_m[getIndex(i + 1, j, k)] -
410 FieldstrengthBx_m[getIndex(i - 1, j, k)]) / hx_m -
411 (FieldstrengthBz_m[getIndex(i, j, k + 1)] -
412 FieldstrengthBz_m[getIndex(i, j, k - 1)]) / hz_m)) / 3;
413 }
414 { // treat cells at i = 1:num_gridpx_m - 2, j > 1, k = num_gridpz_m - 1
415 unsigned int k = num_gridpz_m - 1;
416 unsigned long index = getIndex(i, j, k);
417 FieldstrengthBy_m[index] = (4 * FieldstrengthBy_m[getIndex(i, j - 1, k)] -
418 FieldstrengthBy_m[getIndex(i, j - 2, k)] +
419 hy_m * ((FieldstrengthBx_m[getIndex(i + 1, j, k)] -
420 FieldstrengthBx_m[getIndex(i - 1, j, k)]) / hx_m -
421 (FieldstrengthBz_m[getIndex(i, j, k - 2)] -
422 4 * FieldstrengthBz_m[getIndex(i, j, k - 1)] +
423 3 * FieldstrengthBz_m[getIndex(i, j, k)]) / hz_m)) / 3;
424 }
425 }
426 {
427 unsigned int i = 0;
428 { // treat cells at i = 0, j > 1, k = 0
429 unsigned int k = 0;
430 unsigned long index = getIndex(i, j, k);
431 FieldstrengthBy_m[index] = (4 * FieldstrengthBy_m[getIndex(i, j - 1, k)] -
432 FieldstrengthBy_m[getIndex(i, j - 2, k)] +
433 hy_m * ((-3 * FieldstrengthBx_m[getIndex(i, j, k)] +
434 4 * FieldstrengthBx_m[getIndex(i + 1, j, k)] -
435 FieldstrengthBx_m[getIndex(i + 2, j, k)]) / hx_m -
436 (-3 * FieldstrengthBz_m[getIndex(i, j, k)] +
437 4 * FieldstrengthBz_m[getIndex(i, j, k + 1)] -
438 FieldstrengthBz_m[getIndex(i, j, k + 2)]) / hz_m)) / 3;
439 }
440 // treat cells at i = 0, j > 1, k = 1:num_gridpz_m - 2
441 for (unsigned int k = 1; k < num_gridpz_m - 1; k ++) {
442 unsigned long index = getIndex(i,j,k);
443 FieldstrengthBy_m[index] = (4 * FieldstrengthBy_m[getIndex(i, j - 1, k)] -
444 FieldstrengthBy_m[getIndex(i, j - 2, k)] +
445 hy_m * ((-3 * FieldstrengthBx_m[getIndex(i, j, k)] +
446 4 * FieldstrengthBx_m[getIndex(i + 1, j, k)] -
447 FieldstrengthBx_m[getIndex(i + 2, j, k)]) / hx_m -
448 (FieldstrengthBz_m[getIndex(i, j, k + 1)] -
449 FieldstrengthBz_m[getIndex(i, j, k - 1)]) / hz_m)) / 3;
450 }
451 { // treat cells at i = 0, j > 1, k = num_gridpz_m - 1
452 unsigned int k = num_gridpz_m - 1;
453 unsigned long index = getIndex(i, j, k);
454 FieldstrengthBy_m[index] = (4 * FieldstrengthBy_m[getIndex(i, j - 1, k)] -
455 FieldstrengthBy_m[getIndex(i, j - 2, k)] +
456 hy_m * ((-3 * FieldstrengthBx_m[getIndex(i, j, k)] +
457 4 * FieldstrengthBx_m[getIndex(i + 1, j, k)] -
458 FieldstrengthBx_m[getIndex(i + 2, j, k)]) / hx_m -
459 (FieldstrengthBz_m[getIndex(i, j, k - 2)] -
460 4 * FieldstrengthBz_m[getIndex(i, j, k - 1)] +
461 3 * FieldstrengthBz_m[getIndex(i, j, k)]) / hz_m)) / 3;
462 }
463 }
464 {
465 unsigned int i = num_gridpx_m - 1;
466 { // treat cells at i = 0, j > 1, k = 0
467 unsigned int k = 0;
468 unsigned long index = getIndex(i, j, k);
469 FieldstrengthBy_m[index] = (4 * FieldstrengthBy_m[getIndex(i, j - 1, k)] -
470 FieldstrengthBy_m[getIndex(i, j - 2, k)] +
471 hy_m * ((FieldstrengthBx_m[getIndex(i - 2, j, k)] -
472 4 * FieldstrengthBx_m[getIndex(i - 1, j, k)] +
473 3 * FieldstrengthBx_m[getIndex(i, j, k)]) / hx_m -
474 (-3 * FieldstrengthBz_m[getIndex(i, j, k)] +
475 4 * FieldstrengthBz_m[getIndex(i, j, k + 1)] -
476 FieldstrengthBz_m[getIndex(i, j, k + 2)]) / hz_m)) / 3;
477 }
478 // treat cells at i = 0, j > 1, k = 1:num_gridpz_m - 2
479 for (unsigned int k = 1; k < num_gridpz_m - 1; k ++) {
480 unsigned long index = getIndex(i,j,k);
481 FieldstrengthBy_m[index] = (4 * FieldstrengthBy_m[getIndex(i, j - 1, k)] -
482 FieldstrengthBy_m[getIndex(i, j - 2, k)] +
483 hy_m * ((FieldstrengthBx_m[getIndex(i - 2, j, k)] -
484 4 * FieldstrengthBx_m[getIndex(i - 1, j, k)] +
485 3 * FieldstrengthBx_m[getIndex(i, j, k)]) / hx_m -
486 (FieldstrengthBz_m[getIndex(i, j, k + 1)] -
487 FieldstrengthBz_m[getIndex(i, j, k - 1)]) / hz_m)) / 3;
488 }
489 { // treat cells at i = 0, j > 1, k = num_gridpz_m - 1
490 unsigned int k = num_gridpz_m - 1;
491 unsigned long index = getIndex(i, j, k);
492 FieldstrengthBy_m[index] = (4 * FieldstrengthBy_m[getIndex(i, j - 1, k)] -
493 FieldstrengthBy_m[getIndex(i, j - 2, k)] +
494 hy_m * ((FieldstrengthBx_m[getIndex(i - 2, j, k)] -
495 4 * FieldstrengthBx_m[getIndex(i - 1, j, k)] +
496 3 * FieldstrengthBx_m[getIndex(i, j, k)]) / hx_m -
497 (FieldstrengthBz_m[getIndex(i, j, k - 2)] -
498 4 * FieldstrengthBz_m[getIndex(i, j, k - 1)] +
499 3 * FieldstrengthBz_m[getIndex(i, j, k)]) / hz_m)) / 3;
500 }
501 }
502 }
503}
504
505void _FM3DMagnetoStaticExtended::smoothData(double * data, unsigned j)
506{
507 const double offWeight = 0.1, sumWeightInv = 1.0 / (1.0 + 4 * (1 + offWeight) * offWeight);
508 double *tmp = new double[num_gridpx_m * num_gridpy_m * num_gridpz_m];
509
510 for (unsigned int i = 1; i < num_gridpx_m - 1; ++ i) {
511 for (unsigned int k = 1; k < num_gridpz_m - 1; ++ k) {
512 double sum = 0.0;
513 for (int i2 = -1; i2 < 2; ++ i2) {
514 for (int k2 = -1; k2 < 2; ++ k2) {
515 double weight = std::pow(offWeight, std::abs(i2) + std::abs(k2));
516 sum += weight * data[getIndex(i + i2, j, k + k2)];
517 }
518 }
519
520 tmp[getIndex(i, j, k)] = sum * sumWeightInv;
521 }
522 }
523 {
524 const double sumWeightInv = 1.0 / (1.0 + (3 + 2 * offWeight) * offWeight);
525 unsigned int i = 0;
526 for (unsigned int k = 1; k < num_gridpz_m - 1; ++ k) {
527 double sum = 0.0;
528 for (int i2 = 0; i2 < 2; ++ i2) {
529 for (int k2 = -1; k2 < 2; ++ k2) {
530 double weight = std::pow(offWeight, std::abs(i2) + std::abs(k2));
531 sum += weight * data[getIndex(i + i2, j, k + k2)];
532 }
533 }
534
535 tmp[getIndex(i, j, k)] = sum * sumWeightInv;
536 }
537 }
538 {
539 const double sumWeightInv = 1.0 / (1.0 + (3 + 2 * offWeight) * offWeight);
540 unsigned int i = num_gridpx_m - 1;
541 for (unsigned int k = 1; k < num_gridpz_m - 1; ++ k) {
542 double sum = 0.0;
543 for (int i2 = -1; i2 < 1; ++ i2) {
544 for (int k2 = -1; k2 < 2; ++ k2) {
545 double weight = std::pow(offWeight, std::abs(i2) + std::abs(k2));
546 sum += weight * data[getIndex(i + i2, j, k + k2)];
547 }
548 }
549
550 tmp[getIndex(i, j, k)] = sum * sumWeightInv;
551 }
552 }
553 {
554 const double sumWeightInv = 1.0 / (1.0 + (3 + 2 * offWeight) * offWeight);
555 unsigned int k = 0;
556 for (unsigned int i = 1; i < num_gridpx_m - 1; ++ i) {
557 double sum = 0.0;
558 for (int i2 = -1; i2 < 2; ++ i2) {
559 for (int k2 = 0; k2 < 2; ++ k2) {
560 double weight = std::pow(offWeight, std::abs(i2) + std::abs(k2));
561 sum += weight * data[getIndex(i + i2, j, k + k2)];
562 }
563 }
564
565 tmp[getIndex(i, j, k)] = sum * sumWeightInv;
566 }
567 }
568 {
569 const double sumWeightInv = 1.0 / (1.0 + (3 + 2 * offWeight) * offWeight);
570 unsigned int k = num_gridpz_m - 1;
571 for (unsigned int i = 1; i < num_gridpx_m - 1; ++ i) {
572 double sum = 0.0;
573 for (int i2 = -1; i2 < 2; ++ i2) {
574 for (int k2 = -1; k2 < 1; ++ k2) {
575 double weight = std::pow(offWeight, std::abs(i2) + std::abs(k2));
576 sum += weight * data[getIndex(i + i2, j, k + k2)];
577 }
578 }
579
580 tmp[getIndex(i, j, k)] = sum * sumWeightInv;
581 }
582 }
583 {
584 unsigned long index = getIndex(0, j, 0);
585 tmp[index] = data[index];
586 }
587 {
588 unsigned long index = getIndex(num_gridpx_m - 1, j, 0);
589 tmp[index] = data[index];
590 }
591 {
592 unsigned long index = getIndex(num_gridpx_m - 1, j, num_gridpz_m - 1);
593 tmp[index] = data[index];
594 }
595 {
596 unsigned long index = getIndex(0, j, num_gridpz_m - 1);
597 tmp[index] = data[index];
598 }
599
600 for (unsigned int i = 0; i < num_gridpx_m; ++ i) {
601 for (unsigned int k = 0; k < num_gridpz_m; ++ k) {
602 unsigned long index = getIndex(i, j, k);
603 data[index] = tmp[index];
604 }
605 }
606
607 delete[] tmp;
608}
609
610void _FM3DMagnetoStaticExtended::saveField(const std::string &fname, unsigned int j) const
611{
612 std::ofstream out(fname);
613 out.precision(6);
614 double x = xbegin_m, y = j * hy_m;
615 for (unsigned int i = 0; i < num_gridpx_m; ++ i, x += hx_m) {
616 double z = zbegin_m;
617 for (unsigned int k = 0; k < num_gridpz_m; ++ k, z += hz_m) {
618 unsigned long index = getIndex(i, j, k);
619 out << std::setw(14) << x
620 << std::setw(14) << y
621 << std::setw(14) << z
622 << std::setw(14) << FieldstrengthBx_m[index]
623 << std::setw(14) << FieldstrengthBy_m[index]
624 << std::setw(14) << FieldstrengthBz_m[index]
625 << "\n";
626 }
627 }
628
629 out.close();
630}
631
632void _FM3DMagnetoStaticExtended::freeMap() {
633 if(FieldstrengthBz_m != nullptr) {
634 delete[] FieldstrengthBx_m;
635 delete[] FieldstrengthBy_m;
636 delete[] FieldstrengthBz_m;
637
638 FieldstrengthBx_m = nullptr;
639 FieldstrengthBy_m = nullptr;
640 FieldstrengthBz_m = nullptr;
641 }
642}
643
644Vector_t _FM3DMagnetoStaticExtended::interpolateTrilinearly(const Vector_t &X) const {
645 IndexTriplet idx = getIndex(X);
646 Vector_t B(0.0);
647
648 B(0) = (getWeightedData(FieldstrengthBx_m, idx, LX|LY|LZ) +
649 getWeightedData(FieldstrengthBx_m, idx, LX|LY|HZ) +
650 getWeightedData(FieldstrengthBx_m, idx, LX|HY|LZ) +
651 getWeightedData(FieldstrengthBx_m, idx, LX|HY|HZ) +
652 getWeightedData(FieldstrengthBx_m, idx, HX|LY|LZ) +
653 getWeightedData(FieldstrengthBx_m, idx, HX|LY|HZ) +
654 getWeightedData(FieldstrengthBx_m, idx, HX|HY|LZ) +
655 getWeightedData(FieldstrengthBx_m, idx, HX|HY|HZ));
656
657 B(1) = (getWeightedData(FieldstrengthBy_m, idx, LX|LY|LZ) +
658 getWeightedData(FieldstrengthBy_m, idx, LX|LY|HZ) +
659 getWeightedData(FieldstrengthBy_m, idx, LX|HY|LZ) +
660 getWeightedData(FieldstrengthBy_m, idx, LX|HY|HZ) +
661 getWeightedData(FieldstrengthBy_m, idx, HX|LY|LZ) +
662 getWeightedData(FieldstrengthBy_m, idx, HX|LY|HZ) +
663 getWeightedData(FieldstrengthBy_m, idx, HX|HY|LZ) +
664 getWeightedData(FieldstrengthBy_m, idx, HX|HY|HZ));
665
666 B(2) = (getWeightedData(FieldstrengthBz_m, idx, LX|LY|LZ) +
667 getWeightedData(FieldstrengthBz_m, idx, LX|LY|HZ) +
668 getWeightedData(FieldstrengthBz_m, idx, LX|HY|LZ) +
669 getWeightedData(FieldstrengthBz_m, idx, LX|HY|HZ) +
670 getWeightedData(FieldstrengthBz_m, idx, HX|LY|LZ) +
671 getWeightedData(FieldstrengthBz_m, idx, HX|LY|HZ) +
672 getWeightedData(FieldstrengthBz_m, idx, HX|HY|LZ) +
673 getWeightedData(FieldstrengthBz_m, idx, HX|HY|HZ));
674
675 return B;
676}
677
678double _FM3DMagnetoStaticExtended::getWeightedData(double *data, const IndexTriplet &idx, unsigned short corner) const {
679 unsigned short switchX = ((corner & HX) >> 2), switchY = ((corner & HY) >> 1), switchZ = (corner & HZ);
680 double factorX = 0.5 + (1 - 2 * switchX) * (0.5 - idx.weight(0));
681 double factorY = 0.5 + (1 - 2 * switchY) * (0.5 - idx.weight(1));
682 double factorZ = 0.5 + (1 - 2 * switchZ) * (0.5 - idx.weight(2));
683
684 unsigned long i = idx.i + switchX, j = idx.j + switchY, k = idx.k + switchZ;
685
686 return factorX * factorY * factorZ * data[getIndex(i, j, k)];
687}
688
689bool _FM3DMagnetoStaticExtended::getFieldstrength(const Vector_t &R, Vector_t &/*E*/, Vector_t &B) const {
690 if (isInside(R)) {
691 Vector_t suppB = interpolateTrilinearly(R);
692 suppB(0) *= copysign(1, R(1));
693 suppB(2) *= copysign(1, R(1));
694
695 B += suppB;
696 }
697
698 return false;
699}
700
701bool _FM3DMagnetoStaticExtended::getFieldDerivative(const Vector_t &/*R*/, Vector_t &/*E*/, Vector_t &/*B*/, const DiffDirection &/*dir*/) const {
702 return false;
703}
704
705void _FM3DMagnetoStaticExtended::getFieldDimensions(double &zBegin, double &zEnd) const {
706 zBegin = zbegin_m;
707 zEnd = zend_m;
708}
709
710void _FM3DMagnetoStaticExtended::getFieldDimensions(double &xIni, double &xFinal, double &yIni, double &yFinal, double &zIni, double &zFinal) const {
711 xIni = xbegin_m;
712 xFinal = xend_m;
713 yIni = -yend_m;
714 yFinal = yend_m;
715 zIni = zbegin_m;
716 zFinal = zend_m;
717}
718
719void _FM3DMagnetoStaticExtended::swap() {
720}
721
722void _FM3DMagnetoStaticExtended::getInfo(Inform *msg) {
723 (*msg) << Filename_m << " (3D magnetostatic, extended); zini= " << zbegin_m << " m; zfinal= " << zend_m << " m;" << endl;
724}
725
726double _FM3DMagnetoStaticExtended::getFrequency() const {
727 return 0.0;
728}
729
730void _FM3DMagnetoStaticExtended::setFrequency(double /*freq*/)
731{ ;}
gmsg
Inform * gmsg
Definition: Main.cpp:70
FM3DMagnetoStaticExtended
std::shared_ptr< _FM3DMagnetoStaticExtended > FM3DMagnetoStaticExtended
Definition: Definitions.h:72
T2DMagnetoStatic
@ T2DMagnetoStatic
Definition: Fieldmap.h:29
DiffDirection
DiffDirection
Definition: Fieldmap.h:55
pow
Tps< T > pow(const Tps< T > &x, int y)
Integer power.
Definition: TpsMath.h:76
max
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Definition: ReductionLoc.h:84
X
#define X(arg)
Definition: fftpack.cpp:112
sum
T::PETE_Expr_t::PETE_Return_t sum(const PETE_Expr< T > &expr)
Definition: PETE.h:1111
copysign
PETE_TBTree< FnCopysign, PETE_Scalar< Vektor< T1, Dim > >, typename T2::PETE_Expr_t > copysign(const Vektor< T1, Dim > &l, const PETE_Expr< T2 > &r)
Definition: IpplExpressions.h:244
abs
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
Definition: IpplExpressions.h:62
endl
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
level3
Inform & level3(Inform &inf)
Definition: Inform.cpp:47
INFOMSG
#define INFOMSG(msg)
Definition: IpplInfo.h:348
Physics::e
constexpr double e
The value of.
Definition: Physics.h:39
Units::cm2m
constexpr double cm2m
Definition: Units.h:35
Util::toUpper
std::string toUpper(const std::string &str)
Definition: Util.cpp:147
_Fieldmap::interpreteEOF
bool interpreteEOF(std::ifstream &in)
Definition: Fieldmap.cpp:555
_Fieldmap::normalize_m
bool normalize_m
Definition: Fieldmap.h:121
_Fieldmap::typeset_msg
static std::string typeset_msg(const std::string &msg, const std::string &title)
Definition: Fieldmap.cpp:649
_Fieldmap::disableFieldmapWarning
void disableFieldmapWarning()
Definition: Fieldmap.cpp:610
_Fieldmap::getLine
void getLine(std::ifstream &in, std::string &buffer)
Definition: Fieldmap.h:122
_Fieldmap::lines_read_m
int lines_read_m
Definition: Fieldmap.h:119
_Fieldmap::Filename_m
std::string Filename_m
Definition: Fieldmap.h:118
_Fieldmap::Type
MapType Type
Definition: Fieldmap.h:115
_FM3DMagnetoStaticExtended::isInside
virtual bool isInside(const Vector_t &r) const
Definition: FM3DMagnetoStaticExtended.h:84
_FM3DMagnetoStaticExtended::hy_m
double hy_m
length between points in grid, y-direction, m
Definition: FM3DMagnetoStaticExtended.h:74
_FM3DMagnetoStaticExtended::hz_m
double hz_m
length between points in grid, z-direction, m
Definition: FM3DMagnetoStaticExtended.h:75
_FM3DMagnetoStaticExtended::setFrequency
virtual void setFrequency(double freq)
Definition: FM3DMagnetoStaticExtended.cpp:730
_FM3DMagnetoStaticExtended::getFieldstrength
virtual bool getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &B) const
Definition: FM3DMagnetoStaticExtended.cpp:689
_FM3DMagnetoStaticExtended::num_gridpx_m
unsigned int num_gridpx_m
Read in number of points after 0(not counted here) in grid, x-direction.
Definition: FM3DMagnetoStaticExtended.h:76
_FM3DMagnetoStaticExtended::getWeightedData
double getWeightedData(double *data, const IndexTriplet &idx, unsigned short corner) const
Definition: FM3DMagnetoStaticExtended.cpp:678
_FM3DMagnetoStaticExtended::FieldstrengthBz_m
double * FieldstrengthBz_m
3D array with Bz, read in like Bx
Definition: FM3DMagnetoStaticExtended.h:64
_FM3DMagnetoStaticExtended::getIndex
IndexTriplet getIndex(const Vector_t &X) const
Definition: FM3DMagnetoStaticExtended.h:99
_FM3DMagnetoStaticExtended::smoothData
void smoothData(double *data, unsigned j)
Definition: FM3DMagnetoStaticExtended.cpp:505
_FM3DMagnetoStaticExtended::saveField
void saveField(const std::string &fname, unsigned int j) const
Definition: FM3DMagnetoStaticExtended.cpp:610
_FM3DMagnetoStaticExtended::getFieldDimensions
virtual void getFieldDimensions(double &zBegin, double &zEnd) const
Definition: FM3DMagnetoStaticExtended.cpp:705
_FM3DMagnetoStaticExtended::FieldstrengthBx_m
double * FieldstrengthBx_m
3D array with Bx, read in first along z0 - r0 to rN then z1 - r0 to rN until zN - r0 to rN
Definition: FM3DMagnetoStaticExtended.h:62
_FM3DMagnetoStaticExtended::hx_m
double hx_m
length between points in grid, x-direction, m
Definition: FM3DMagnetoStaticExtended.h:73
_FM3DMagnetoStaticExtended::num_gridpz_m
unsigned int num_gridpz_m
Read in number of points after 0(not counted here) in grid, z-direction.
Definition: FM3DMagnetoStaticExtended.h:78
_FM3DMagnetoStaticExtended::FieldstrengthBy_m
double * FieldstrengthBy_m
3D array with By, read in like Bx
Definition: FM3DMagnetoStaticExtended.h:63
_FM3DMagnetoStaticExtended::getFieldDerivative
virtual bool getFieldDerivative(const Vector_t &R, Vector_t &E, Vector_t &B, const DiffDirection &dir) const
Definition: FM3DMagnetoStaticExtended.cpp:701
_FM3DMagnetoStaticExtended::interpolateTrilinearly
Vector_t interpolateTrilinearly(const Vector_t &X) const
Definition: FM3DMagnetoStaticExtended.cpp:644
_FM3DMagnetoStaticExtended::_FM3DMagnetoStaticExtended
_FM3DMagnetoStaticExtended(const std::string &filename)
Definition: FM3DMagnetoStaticExtended.cpp:14
_FM3DMagnetoStaticExtended::num_gridpy_m
unsigned int num_gridpy_m
Read in number of points after 0(not counted here) in grid, y-direction.
Definition: FM3DMagnetoStaticExtended.h:77
_FM3DMagnetoStaticExtended::create
static FM3DMagnetoStaticExtended create(const std::string &filename)
Definition: FM3DMagnetoStaticExtended.cpp:94
Inform
Definition: Inform.h:42