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