OPAL (Object Oriented Parallel Accelerator Library)  2.2.0
OPAL
FM3DMagnetoStaticExtended.cpp
Go to the documentation of this file.
1 #include "Fields/FM3DMagnetoStaticExtended.h"
2 #include "Fields/Fieldmap.hpp"
3 #include "Utilities/GeneralClassicException.h"
4 #include "Utilities/Util.h"
5 
6 #include <fstream>
7 #include <ios>
8 #include <algorithm>
9 
10 extern Inform *gmsg;
11 
12 FM3DMagnetoStaticExtended::FM3DMagnetoStaticExtended(std::string aFilename):
13  Fieldmap(aFilename),
14  FieldstrengthBx_m(NULL),
15  FieldstrengthBy_m(NULL),
16  FieldstrengthBz_m(NULL) {
17  std::ifstream file;
18  std::string tmpString;
19  double tmpDouble;
20 
21  Type = T2DMagnetoStatic;
22 
23  // open field map, parse it and disable element on error
24  file.open(Filename_m.c_str());
25  if(file.good()) {
26  bool parsing_passed = true;
27  try {
28  parsing_passed = interpreteLine<std::string>(file, tmpString);
29  } catch (GeneralClassicException &e) {
30  parsing_passed = interpreteLine<std::string, std::string>(file, tmpString, tmpString);
31 
32  tmpString = Util::toUpper(tmpString);
33  if (tmpString != "TRUE" &&
34  tmpString != "FALSE")
35  throw GeneralClassicException("FM3DMagnetoStaticExtended::FM3DMagnetoStaticExtended",
36  "The second string on the first line of 3D field "
37  "maps has to be either TRUE or FALSE");
38 
39  normalize_m = (tmpString == "TRUE");
40  }
41  parsing_passed = (parsing_passed &&
42  interpreteLine<double, double, unsigned int>(file, xbegin_m, xend_m, num_gridpx_m));
43  parsing_passed = (parsing_passed &&
44  interpreteLine<double, double, unsigned int>(file, ybegin_m, yend_m, num_gridpy_m));
45  parsing_passed = (parsing_passed &&
46  interpreteLine<double, double, unsigned int>(file, zbegin_m, zend_m, num_gridpz_m));
47 
48  for(unsigned long i = 0; (i < (num_gridpz_m + 1) * (num_gridpx_m + 1)) && parsing_passed; ++ i) {
49  parsing_passed = parsing_passed && interpreteLine<double>(file, tmpDouble);
50  }
51 
52  parsing_passed = parsing_passed &&
53  interpreteEOF(file);
54 
55  file.close();
56  lines_read_m = 0;
57 
58  if(!parsing_passed) {
59  disableFieldmapWarning();
60  zend_m = zbegin_m - 1e-3;
61  throw GeneralClassicException("FM3DMagnetoStaticExtended::FM3DMagnetoStaticExtended(std::string)",
62  "Format of fieldmap '" + Filename_m + "' didn't pass basic test");
63  } else {
64  // conversion from cm to m
65  xbegin_m /= 100.;
66  xend_m /= 100.;
67  yend_m = std::max(std::abs(ybegin_m), yend_m) / 100.;
68  ybegin_m = 0.0;
69  zbegin_m /= 100.;
70  zend_m /= 100.;
71  length_m = zend_m - zbegin_m;
72 
73  hx_m = (xend_m - xbegin_m) / num_gridpx_m;
74  hy_m = (yend_m - ybegin_m) / num_gridpy_m;
75  hz_m = (zend_m - zbegin_m) / num_gridpz_m;
76 
77  // num spacings -> num grid points
78  num_gridpx_m++;
79  num_gridpy_m++;
80  num_gridpz_m++;
81  }
82  } else {
83  throw GeneralClassicException("FM3DMagnetoStaticExtended::FM3DMagnetoStaticExtended(std::string)",
84  "Couldn't read fieldmap '" + Filename_m + "'");
85  }
86 }
87 
88 FM3DMagnetoStaticExtended::~FM3DMagnetoStaticExtended() {
89  freeMap();
90 }
91 
92 void FM3DMagnetoStaticExtended::readMap() {
93  if(FieldstrengthBz_m == NULL) {
94  // declare variables and allocate memory
95  std::ifstream in;
96  std::string tmpString;
97  const size_t totalSize = num_gridpx_m * num_gridpy_m * num_gridpz_m;
98  FieldstrengthBx_m = new double[totalSize];
99  FieldstrengthBy_m = new double[totalSize];
100  FieldstrengthBz_m = new double[totalSize];
101 
102  // read in and parse field map
103  in.open(Filename_m.c_str());
104  getLine(in, tmpString);
105  getLine(in, tmpString);
106  getLine(in, tmpString);
107  getLine(in, tmpString);
108 
109  for(unsigned int i = 0; i < num_gridpx_m; i++) {
110  for(unsigned int k = 0; k < num_gridpz_m; k++) {
111  unsigned long index = getIndex(i,0,k);
112  interpreteLine<double>(in, FieldstrengthBy_m[index]);
113  }
114  }
115  in.close();
116 
117  std::fill(FieldstrengthBx_m, FieldstrengthBx_m + totalSize, 0.0);
118  std::fill(FieldstrengthBz_m, FieldstrengthBz_m + totalSize, 0.0);
119 
120  if (normalize_m) {
121  double Bymax = 0.0;
122 
123  // find maximum field
124  unsigned int centerX = static_cast<unsigned int>(std::floor(-xbegin_m / hx_m + 0.5));
125  for(unsigned int k = 0; k < num_gridpz_m; ++ k) {
126  double By = FieldstrengthBy_m[getIndex(centerX, 0, k)];
127  if(std::abs(By) > std::abs(Bymax)) {
128  Bymax = By;
129  }
130  }
131 
132  // normalize field
133  for(unsigned int i = 0; i < num_gridpx_m; i ++) {
134  for(unsigned int k = 0; k < num_gridpz_m; k ++) {
135  unsigned long index = getIndex(i,0,k);
136  FieldstrengthBy_m[index] /= Bymax;
137  }
138  }
139  }
140 
141  smoothData(FieldstrengthBy_m, 0);
142  for (unsigned int j = 1; j < num_gridpy_m; j ++) {
143  integrateBx(j);
144  integrateBz(j);
145  integrateBy(j);
146 
147  smoothData(FieldstrengthBx_m, j);
148  smoothData(FieldstrengthBz_m, j);
149  smoothData(FieldstrengthBy_m, j);
150  }
151 
152  // saveField("data/" + Filename_m + ".ext", num_gridpy_m - 1);
153 
154  INFOMSG(level3 << typeset_msg("read in fieldmap '" + Filename_m + "'", "info") << endl);
155  }
156 }
157 
158 void 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 
225 void 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 
283 void 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 
500 void 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 
605 void 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 
627 void FM3DMagnetoStaticExtended::freeMap() {
628  if(FieldstrengthBz_m != NULL) {
629  delete[] FieldstrengthBx_m;
630  delete[] FieldstrengthBy_m;
631  delete[] FieldstrengthBz_m;
632 
633  FieldstrengthBx_m = NULL;
634  FieldstrengthBy_m = NULL;
635  FieldstrengthBz_m = NULL;
636 
637  INFOMSG(level3 << typeset_msg("freed fieldmap '" + Filename_m + "'", "info") << endl);
638  }
639 }
640 
641 Vector_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 
675 double 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 
686 bool 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 
698 bool FM3DMagnetoStaticExtended::getFieldDerivative(const Vector_t &R, Vector_t &E, Vector_t &B, const DiffDirection &dir) const {
699  return false;
700 }
701 
702 void FM3DMagnetoStaticExtended::getFieldDimensions(double &zBegin, double &zEnd, double &xbegin, double &xend) const {
703  zBegin = zbegin_m;
704  zEnd = zend_m;
705  xbegin = xbegin_m;
706  xend = xend_m;
707 }
708 
709 void FM3DMagnetoStaticExtended::getFieldDimensions(double &xIni, double &xFinal, double &yIni, double &yFinal, double &zIni, double &zFinal) const {
710  xIni = xbegin_m;
711  xFinal = xend_m;
712  yIni = -yend_m;
713  yFinal = yend_m;
714  zIni = zbegin_m;
715  zFinal = zend_m;
716 }
717 
718 void FM3DMagnetoStaticExtended::swap() {
719 }
720 
721 void FM3DMagnetoStaticExtended::getInfo(Inform *msg) {
722  (*msg) << Filename_m << " (3D magnetostatic, extended); zini= " << zbegin_m << " m; zfinal= " << zend_m << " m;" << endl;
723 }
724 
725 double FM3DMagnetoStaticExtended::getFrequency() const {
726  return 0.0;
727 }
728 
729 void FM3DMagnetoStaticExtended::setFrequency(double freq)
730 { ;}
FM3DMagnetoStaticExtended::getFieldDerivative
virtual bool getFieldDerivative(const Vector_t &R, Vector_t &E, Vector_t &B, const DiffDirection &dir) const
Definition: FM3DMagnetoStaticExtended.cpp:698
FM3DMagnetoStaticExtended::getIndex
IndexTriplet getIndex(const Vector_t &X) const
Definition: FM3DMagnetoStaticExtended.h:97
abs
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
Definition: IpplExpressions.h:62
FM3DMagnetoStaticExtended::interpolateTrilinearly
Vector_t interpolateTrilinearly(const Vector_t &X) const
Definition: FM3DMagnetoStaticExtended.cpp:641
Physics::e
constexpr double e
The value of .
Definition: Physics.h:40
DiffDirection
DiffDirection
Definition: Fieldmap.h:54
max
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Definition: ReductionLoc.h:123
X
Definition: rbendmap.h:8
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
gmsg
Inform * gmsg
Definition: Main.cpp:21
Util::toUpper
std::string toUpper(const std::string &str)
Definition: Util.cpp:130
Fieldmap::disableFieldmapWarning
void disableFieldmapWarning()
Definition: Fieldmap.cpp:610
FM3DMagnetoStaticExtended::isInside
virtual bool isInside(const Vector_t &r) const
Definition: FM3DMagnetoStaticExtended.h:81
Fieldmap::typeset_msg
static std::string typeset_msg(const std::string &msg, const std::string &title)
Definition: Fieldmap.cpp:647
FM3DMagnetoStaticExtended::getFieldDimensions
virtual void getFieldDimensions(double &zBegin, double &zEnd, double &rBegin, double &rEnd) const
Definition: FM3DMagnetoStaticExtended.cpp:702
sum
T::PETE_Expr_t::PETE_Return_t sum(const PETE_Expr< T > &expr)
Definition: PETE.h:1213
Fieldmap::normalize_m
bool normalize_m
Definition: Fieldmap.h:113
Fieldmap::lines_read_m
int lines_read_m
Definition: Fieldmap.h:111
FM3DMagnetoStaticExtended::saveField
void saveField(const std::string &fname, unsigned int j) const
Definition: FM3DMagnetoStaticExtended.cpp:605
FM3DMagnetoStaticExtended::getInfo
virtual void getInfo(Inform *msg)
Definition: FM3DMagnetoStaticExtended.cpp:721
FM3DMagnetoStaticExtended::getFrequency
virtual double getFrequency() const
Definition: FM3DMagnetoStaticExtended.cpp:725
FM3DMagnetoStaticExtended::FM3DMagnetoStaticExtended
FM3DMagnetoStaticExtended(std::string aFilename)
Definition: FM3DMagnetoStaticExtended.cpp:12
Fieldmap::interpreteEOF
bool interpreteEOF(std::ifstream &in)
Definition: Fieldmap.cpp:547
FM3DMagnetoStaticExtended::getWeightedData
double getWeightedData(double *data, const IndexTriplet &idx, unsigned short corner) const
Definition: FM3DMagnetoStaticExtended.cpp:675
INFOMSG
#define INFOMSG(msg)
Definition: IpplInfo.h:397
FM3DMagnetoStaticExtended::getFieldstrength
virtual bool getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &B) const
Definition: FM3DMagnetoStaticExtended.cpp:686
pow
Tps< T > pow(const Tps< T > &x, int y)
Integer power.
Definition: TpsMath.h:76
Fieldmap::getLine
void getLine(std::ifstream &in, std::string &buffer)
Definition: Fieldmap.h:198
Fieldmap::Filename_m
std::string Filename_m
Definition: Fieldmap.h:110
Fieldmap::Type
MapType Type
Definition: Fieldmap.h:107
level3
Inform & level3(Inform &inf)
Definition: Inform.cpp:47
Inform
Definition: Inform.h:41
FM3DMagnetoStaticExtended::setFrequency
virtual void setFrequency(double freq)
Definition: FM3DMagnetoStaticExtended.cpp:729
FM3DMagnetoStaticExtended::smoothData
void smoothData(double *data, unsigned j)
Definition: FM3DMagnetoStaticExtended.cpp:500
floor
PETE_TUTree< FnFloor, typename T::PETE_Expr_t > floor(const PETE_Expr< T > &l)
Definition: PETE.h:816
endl
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
Generated on Sat Jan 4 2020 14:47:55 for OPAL (Object Oriented Parallel Accelerator Library) by   doxygen 1.8.5