OPAL (Object Oriented Parallel Accelerator Library)  2021.1.99
OPAL
FM3DMagnetoStatic.cpp
Go to the documentation of this file.
2 #include "Fields/Fieldmap.hpp"
4 #include "Utilities/Util.h"
5 #include "Physics/Physics.h"
6 
7 #include <fstream>
8 #include <ios>
9 
10 
11 FM3DMagnetoStatic::FM3DMagnetoStatic(std::string aFilename):
12  Fieldmap(aFilename),
13  FieldstrengthBz_m(NULL),
14  FieldstrengthBx_m(NULL),
15  FieldstrengthBy_m(NULL) {
16 
17  std::string tmpString;
18  double tmpDouble;
19 
21  std::ifstream file(Filename_m.c_str());
22 
23  if(file.good()) {
24  bool parsing_passed = true;
25  try {
26  interpretLine<std::string>(file, tmpString);
27  } catch (GeneralClassicException &e) {
28  parsing_passed = interpretLine<std::string, std::string>(file, tmpString, tmpString);
29 
30  tmpString = Util::toUpper(tmpString);
31  if (tmpString != "TRUE" &&
32  tmpString != "FALSE")
33  throw GeneralClassicException("FM3DMagnetoStatic::FM3DMagnetoStatic",
34  "The second string on the first line of 3D field "
35  "maps has to be either TRUE or FALSE");
36 
37  normalize_m = (tmpString == "TRUE");
38  }
39 
40  parsing_passed = parsing_passed &&
41  interpretLine<double, double, unsigned int>(file, xbegin_m, xend_m, num_gridpx_m);
42  parsing_passed = parsing_passed &&
43  interpretLine<double, double, unsigned int>(file, ybegin_m, yend_m, num_gridpy_m);
44  parsing_passed = parsing_passed &&
45  interpretLine<double, double, unsigned int>(file, zbegin_m, zend_m, num_gridpz_m);
46 
47  for(unsigned long i = 0; (i < (num_gridpz_m + 1) * (num_gridpx_m + 1) * (num_gridpy_m + 1)) && parsing_passed; ++ i) {
48  parsing_passed = parsing_passed &&
49  interpretLine<double>(file,
50  tmpDouble,
51  tmpDouble,
52  tmpDouble);
53  }
54 
55  parsing_passed = parsing_passed &&
56  interpreteEOF(file);
57 
58  file.close();
59 
60  if(!parsing_passed) {
61  throw GeneralClassicException("FM3DMagnetoStatic::FM3DMagnetoStatic",
62  "An error occured when reading the fieldmap '" + Filename_m + "'");
63  } else {
64  xbegin_m /= 100.0;
65  xend_m /= 100.0;
66  ybegin_m /= 100.0;
67  yend_m /= 100.0;
68  zbegin_m /= 100.0;
69  zend_m /= 100.0;
70 
74 
75  ++ num_gridpx_m;
76  ++ num_gridpy_m;
77  ++ num_gridpz_m;
78 
79  }
80  } else {
81  throw GeneralClassicException("FM3DMagnetoStatic::FM3DMagnetoStatic",
82  "An error occured when reading the fieldmap '" + Filename_m + "'");
83  }
84 }
85 
86 
88  freeMap();
89 }
90 
92  if(FieldstrengthBz_m == NULL) {
93 
94  std::ifstream in(Filename_m.c_str());
95  std::string tmpString;
96  const size_t totalSize = num_gridpx_m * num_gridpy_m * num_gridpz_m;
97 
98  getLine(in, tmpString);
99  getLine(in, tmpString);
100  getLine(in, tmpString);
101  getLine(in, tmpString);
102 
103  FieldstrengthBz_m = new double[totalSize];
104  FieldstrengthBx_m = new double[totalSize];
105  FieldstrengthBy_m = new double[totalSize];
106 
107  long ii = 0;
108  for(unsigned int i = 0; i < num_gridpx_m; ++ i) {
109  for(unsigned int j = 0; j < num_gridpy_m; ++ j) {
110  for(unsigned int k = 0; k < num_gridpz_m; ++ k) {
111  interpretLine<double>(in,
112  FieldstrengthBx_m[ii],
113  FieldstrengthBy_m[ii],
114  FieldstrengthBz_m[ii]);
115  ++ ii;
116  }
117  }
118  }
119  in.close();
120 
121  if (normalize_m) {
122  double Bymax = 0.0;
123  // find maximum field
124  unsigned int centerX = static_cast<unsigned int>(std::round(-xbegin_m / hx_m));
125  unsigned int centerY = static_cast<unsigned int>(std::round(-ybegin_m / hy_m));
126  for(unsigned int k = 0; k < num_gridpz_m; ++ k) {
127  double By = FieldstrengthBy_m[getIndex(centerX, centerY, k)];
128  if(std::abs(By) > Bymax) {
129  Bymax = std::abs(By);
130  }
131  }
132 
133  // normalize field
134  for(size_t i = 0; i < totalSize; ++ i) {
135  FieldstrengthBx_m[i] /= Bymax;
136  FieldstrengthBy_m[i] /= Bymax;
137  FieldstrengthBz_m[i] /= Bymax;
138  }
139  }
140 
141  INFOMSG(level3 << typeset_msg("read in fieldmap '" + Filename_m + "'", "info") << "\n"
142  << endl);
143  }
144 }
145 
147  if(FieldstrengthBz_m != NULL) {
148  delete[] FieldstrengthBz_m;
149  delete[] FieldstrengthBx_m;
150  delete[] FieldstrengthBy_m;
151 
152  FieldstrengthBz_m = NULL;
153  FieldstrengthBx_m = NULL;
154  FieldstrengthBy_m = NULL;
155 
156  INFOMSG(level3 << typeset_msg("freed fieldmap '" + Filename_m + "'", "info") << "\n"
157  << endl);
158  }
159 }
160 
162  if (!isInside(R)) {
163  return true;
164  }
165 
167 
168  return false;
169 }
170 
172  IndexTriplet idx = getIndex(X);
173  Vector_t B(0.0);
174 
175  B(0) = (getWeightedData(FieldstrengthBx_m, idx, LX|LY|LZ) +
183 
184  B(1) = (getWeightedData(FieldstrengthBy_m, idx, LX|LY|LZ) +
192 
193  B(2) = (getWeightedData(FieldstrengthBz_m, idx, LX|LY|LZ) +
201 
202  return B;
203 }
204 
205 double FM3DMagnetoStatic::getWeightedData(double *data, const IndexTriplet &idx, unsigned short corner) const {
206  unsigned short switchX = ((corner & HX) >> 2), switchY = ((corner & HY) >> 1), switchZ = (corner & HZ);
207  double factorX = 0.5 + (1 - 2 * switchX) * (0.5 - idx.weight(0));
208  double factorY = 0.5 + (1 - 2 * switchY) * (0.5 - idx.weight(1));
209  double factorZ = 0.5 + (1 - 2 * switchZ) * (0.5 - idx.weight(2));
210 
211  unsigned long i = idx.i + switchX, j = idx.j + switchY, k = idx.k + switchZ;
212 
213  return factorX * factorY * factorZ * data[getIndex(i, j, k)];
214 }
215 
216 bool FM3DMagnetoStatic::getFieldDerivative(const Vector_t &/*R*/, Vector_t &/*E*/, Vector_t &/*B*/, const DiffDirection &/*dir*/) const {
217  return false;
218 }
219 
220 void FM3DMagnetoStatic::getFieldDimensions(double &zBegin, double &zEnd) const {
221  zBegin = zbegin_m;
222  zEnd = zend_m;
223 }
224 void FM3DMagnetoStatic::getFieldDimensions(double &/*xIni*/, double &/*xFinal*/, double &/*yIni*/, double &/*yFinal*/, double &/*zIni*/, double &/*zFinal*/) const {}
225 
227 
229  (*msg) << Filename_m << " (3D magnetostatic) "
230  << " xini= " << xbegin_m << " xfinal= " << xend_m
231  << " yini= " << ybegin_m << " yfinal= " << yend_m
232  << " zini= " << zbegin_m << " zfinal= " << zend_m << " (m) " << endl;
233 }
@ T3DMagnetoStatic
Definition: Fieldmap.h:32
DiffDirection
Definition: Fieldmap.h:54
#define X(arg)
Definition: fftpack.cpp:112
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
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
bool normalize_m
Definition: Fieldmap.h:120
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
unsigned int num_gridpz_m
virtual void getFieldDimensions(double &zBegin, double &zEnd) const
IndexTriplet getIndex(const Vector_t &X) const
virtual bool isInside(const Vector_t &r) const
virtual bool getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &B) const
virtual void freeMap()
double getWeightedData(double *data, const IndexTriplet &idx, unsigned short corner) const
virtual bool getFieldDerivative(const Vector_t &R, Vector_t &E, Vector_t &B, const DiffDirection &dir) const
unsigned int num_gridpx_m
virtual void readMap()
Vector_t interpolateTrilinearly(const Vector_t &X) const
virtual void getInfo(Inform *msg)
FM3DMagnetoStatic(std::string aFilename)
unsigned int num_gridpy_m
Definition: Inform.h:42