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