OPAL (Object Oriented Parallel Accelerator Library)  2024.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 
12 extern 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 
94 FM3DMagnetoStaticExtended _FM3DMagnetoStaticExtended::create(const std::string& filename)
95 {
96  return FM3DMagnetoStaticExtended(new _FM3DMagnetoStaticExtended(filename));
97 }
98 
99 void _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 
163 void _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 
230 void _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 
288 void _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 
505 void _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 
610 void _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 
632 void _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 
644 Vector_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 
678 double _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 
689 bool _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 
701 bool _FM3DMagnetoStaticExtended::getFieldDerivative(const Vector_t &/*R*/, Vector_t &/*E*/, Vector_t &/*B*/, const DiffDirection &/*dir*/) const {
702  return false;
703 }
704 
705 void _FM3DMagnetoStaticExtended::getFieldDimensions(double &zBegin, double &zEnd) const {
706  zBegin = zbegin_m;
707  zEnd = zend_m;
708 }
709 
710 void _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 
719 void _FM3DMagnetoStaticExtended::swap() {
720 }
721 
722 void _FM3DMagnetoStaticExtended::getInfo(Inform *msg) {
723  (*msg) << Filename_m << " (3D magnetostatic, extended); zini= " << zbegin_m << " m; zfinal= " << zend_m << " m;" << endl;
724 }
725 
726 double _FM3DMagnetoStaticExtended::getFrequency() const {
727  return 0.0;
728 }
729 
730 void _FM3DMagnetoStaticExtended::setFrequency(double /*freq*/)
731 { ;}
_Fieldmap::Type
MapType Type
Definition: Fieldmap.h:115
_FM3DMagnetoStaticExtended::saveField
void saveField(const std::string &fname, unsigned int j) const
Definition: FM3DMagnetoStaticExtended.cpp:610
DiffDirection
DiffDirection
Definition: Fieldmap.h:55
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
Units::cm2m
constexpr double cm2m
Definition: Units.h:35
_FM3DMagnetoStaticExtended::interpolateTrilinearly
Vector_t interpolateTrilinearly(const Vector_t &X) const
Definition: FM3DMagnetoStaticExtended.cpp:644
_Fieldmap::lines_read_m
int lines_read_m
Definition: Fieldmap.h:119
_Fieldmap::typeset_msg
static std::string typeset_msg(const std::string &msg, const std::string &title)
Definition: Fieldmap.cpp:649
_FM3DMagnetoStaticExtended::getWeightedData
double getWeightedData(double *data, const IndexTriplet &idx, unsigned short corner) const
Definition: FM3DMagnetoStaticExtended.cpp:678
Util::toUpper
std::string toUpper(const std::string &str)
Definition: Util.cpp:147
INFOMSG
#define INFOMSG(msg)
Definition: IpplInfo.h:348
endl
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
_Fieldmap::interpreteEOF
bool interpreteEOF(std::ifstream &in)
Definition: Fieldmap.cpp:555
_FM3DMagnetoStaticExtended::getIndex
IndexTriplet getIndex(const Vector_t &X) const
Definition: FM3DMagnetoStaticExtended.h:99
_FM3DMagnetoStaticExtended::getFieldDerivative
virtual bool getFieldDerivative(const Vector_t &R, Vector_t &E, Vector_t &B, const DiffDirection &dir) const
Definition: FM3DMagnetoStaticExtended.cpp:701
max
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Definition: ReductionLoc.h:84
_FM3DMagnetoStaticExtended::getFieldDimensions
virtual void getFieldDimensions(double &zBegin, double &zEnd) const
Definition: FM3DMagnetoStaticExtended.cpp:705
sum
T::PETE_Expr_t::PETE_Return_t sum(const PETE_Expr< T > &expr)
Definition: PETE.h:1111
_Fieldmap::normalize_m
bool normalize_m
Definition: Fieldmap.h:121
Inform
Definition: Inform.h:42
_FM3DMagnetoStaticExtended::getFieldstrength
virtual bool getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &B) const
Definition: FM3DMagnetoStaticExtended.cpp:689
FM3DMagnetoStaticExtended
std::shared_ptr< _FM3DMagnetoStaticExtended > FM3DMagnetoStaticExtended
Definition: Definitions.h:72
X
#define X(arg)
Definition: fftpack.cpp:112
_Fieldmap::disableFieldmapWarning
void disableFieldmapWarning()
Definition: Fieldmap.cpp:610
_FM3DMagnetoStaticExtended::isInside
virtual bool isInside(const Vector_t &r) const
Definition: FM3DMagnetoStaticExtended.h:84
_Fieldmap::getLine
void getLine(std::ifstream &in, std::string &buffer)
Definition: Fieldmap.h:122
level3
Inform & level3(Inform &inf)
Definition: Inform.cpp:47
Physics::e
constexpr double e
The value of .
Definition: Physics.h:39
_Fieldmap::Filename_m
std::string Filename_m
Definition: Fieldmap.h:118
_FM3DMagnetoStaticExtended::setFrequency
virtual void setFrequency(double freq)
Definition: FM3DMagnetoStaticExtended.cpp:730
_FM3DMagnetoStaticExtended::smoothData
void smoothData(double *data, unsigned j)
Definition: FM3DMagnetoStaticExtended.cpp:505
_FM3DMagnetoStaticExtended::_FM3DMagnetoStaticExtended
_FM3DMagnetoStaticExtended(const std::string &filename)
Definition: FM3DMagnetoStaticExtended.cpp:14
pow
Tps< T > pow(const Tps< T > &x, int y)
Integer power.
Definition: TpsMath.h:76
_FM3DMagnetoStaticExtended::create
static FM3DMagnetoStaticExtended create(const std::string &filename)
Definition: FM3DMagnetoStaticExtended.cpp:94
gmsg
Inform * gmsg
Definition: Main.cpp:70
Generated on Tue Jul 2 2024 15:05:21 for OPAL (Object Oriented Parallel Accelerator Library) by   doxygen 1.8.5