Garfield++ 3.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
ComponentGrid.cc
Go to the documentation of this file.
1#include <limits>
2#include <cmath>
3#include <algorithm>
4#include <fstream>
5#include <iostream>
6#include <cstdio>
7#include <sstream>
8#include <set>
9#include <bitset>
10
11#include "Garfield/ComponentGrid.hh"
12#include "Garfield/Utilities.hh"
13#include "Garfield/GarfieldConstants.hh"
14
15namespace {
16
17unsigned int GetFormat(std::string format) {
18
19 std::transform(format.begin(), format.end(), format.begin(), toupper);
20 unsigned int fmt = 0;
21 if (format == "XY") {
22 fmt = 1;
23 } else if (format == "XYZ") {
24 fmt = 2;
25 } else if (format == "IJ") {
26 fmt = 3;
27 } else if (format == "IJK") {
28 fmt = 4;
29 } else if (format == "YXZ") {
30 fmt = 5;
31 }
32 return fmt;
33}
34
35void PrintError(const std::string& fcn, const unsigned int line,
36 const std::string& par) {
37
38 std::cerr << fcn << ": Error reading line " << line << ".\n"
39 << " Could not read " << par << ".\n";
40}
41
42void PrintNotReady(const std::string& fcn) {
43 std::cerr << fcn << ": Field map not available.\n";
44}
45
46void PrintProgress(const double f) {
47
48 if (f < 0.) return;
49 constexpr unsigned int width = 70;
50 const unsigned int n = static_cast<unsigned int>(std::floor(width * f));
51 std::string bar = "[";
52 if (n < 1) {
53 bar += std::string(width, ' ');
54 } else if (n >= width) {
55 bar += std::string(width, '=');
56 } else {
57 bar += std::string(n, '=') + ">" + std::string(width - n - 1, ' ');
58 }
59 bar += "]";
60 std::cout << bar << "\r" << std::flush;
61}
62
63}
64
65namespace Garfield {
66
67ComponentGrid::ComponentGrid() : ComponentBase() {
68 m_className = "ComponentGrid";
69}
70
71void ComponentGrid::ElectricField(const double x, const double y,
72 const double z, double& ex, double& ey,
73 double& ez, double& p, Medium*& m,
74 int& status) {
75 m = nullptr;
76 status = 0;
77
78 // Make sure the field map has been loaded.
79 if (!m_ready) {
80 PrintNotReady(m_className + "::ElectricField");
81 status = -10;
82 return;
83 }
84
85 status = 0;
86 bool active = true;
87 if (!GetField(x, y, z, m_efields, ex, ey, ez, p, active)) {
88 status = -11;
89 return;
90 }
91 if (!active) {
92 status = -5;
93 return;
94 }
95 m = m_medium;
96 if (!m) status = -5;
97}
98
99void ComponentGrid::ElectricField(const double x, const double y,
100 const double z, double& ex, double& ey,
101 double& ez, Medium*& m, int& status) {
102 double v = 0.;
103 ElectricField(x, y, z, ex, ey, ez, v, m, status);
104}
105
106void ComponentGrid::WeightingField(const double x, const double y,
107 const double z, double& wx, double& wy,
108 double& wz, const std::string& /*label*/) {
109
110 wx = wy = wz = 0.;
111 if (!m_hasWfield) return;
112 const double xx = x - m_wField_xOffset;
113 const double yy = y - m_wField_yOffset;
114 const double zz = z - m_wField_zOffset;
115 double wp = 0.;
116 bool active = true;
117 GetField(xx, yy, zz, m_wfields, wx, wy, wz, wp, active);
118}
119
120double ComponentGrid::WeightingPotential(const double x, const double y,
121 const double z,
122 const std::string& /*label*/) {
123 if (!m_hasWfield) return 0.;
124 const double xx = x - m_wField_xOffset;
125 const double yy = y - m_wField_yOffset;
126 const double zz = z - m_wField_zOffset;
127 double wx = 0., wy = 0., wz = 0.;
128 double wp = 0.;
129 bool active = true;
130 if (!GetField(xx, yy, zz, m_wfields, wx, wy, wz, wp, active)) return 0.;
131 return wp;
132}
133
134void ComponentGrid::DelayedWeightingField(const double x, const double y,
135 const double z, const double t,
136 double& wx, double& wy, double& wz,
137 const std::string& /*label*/) {
138
139 wx = wy = wz = 0.;
140 if (m_wdtimes.empty()) return;
141 // Assume no weighting field for times outside the range of available maps.
142 if (t < m_wdtimes.front() || t > m_wdtimes.back()) return;
143
144 const double xx = x - m_wField_xOffset;
145 const double yy = y - m_wField_yOffset;
146 const double zz = z - m_wField_zOffset;
147
148 const auto it1 = std::upper_bound(m_wdtimes.cbegin(), m_wdtimes.cend(), t);
149 const auto it0 = std::prev(it1);
150
151 const double dt = t - *it0;
152 double wp = 0.;
153 const unsigned int i0 = it0 - m_wdtimes.cbegin();
154 double wx0 = 0., wy0 = 0., wz0 = 0.;
155 bool active = true;
156 if (!GetField(xx, yy, zz, m_wdfields[i0], wx0, wy0, wz0, wp, active)) return;
157
158 if (dt < Small || it1 == m_wdtimes.cend()) {
159 wx = wx0;
160 wy = wy0;
161 wz = wz0;
162 return;
163 }
164 const unsigned int i1 = it1 - m_wdtimes.cbegin();
165 double wx1 = 0., wy1 = 0., wz1 = 0.;
166 if (!GetField(xx, yy, zz, m_wdfields[i1], wx1, wy1, wz1, wp, active)) return;
167
168 const double f1 = dt / (*it1 - *it0);
169 const double f0 = 1. - f1;
170 wx = f0 * wx0 + f1 * wx1;
171 wy = f0 * wy0 + f1 * wy1;
172 wz = f0 * wz0 + f1 * wz1;
173}
174
175void ComponentGrid::SetWeightingFieldOffset(const double x, const double y,
176 const double z) {
177 m_wField_xOffset = x;
178 m_wField_yOffset = y;
179 m_wField_zOffset = z;
180}
181
182void ComponentGrid::MagneticField(const double x, const double y,
183 const double z, double& bx, double& by,
184 double& bz, int& status) {
185 status = 0;
186 if (!m_hasBfield) {
187 return ComponentBase::MagneticField(x, y, z, bx, by, bz, status);
188 }
189
190 double p = 0.;
191 bool active = true;
192 if (!GetField(x, y, z, m_bfields, bx, by, bz, p, active)) {
193 status = -11;
194 }
195}
196
197Medium* ComponentGrid::GetMedium(const double x, const double y,
198 const double z) {
199 // Make sure the field map has been loaded.
200 if (!m_ready) {
201 PrintNotReady(m_className + "::GetMedium");
202 return nullptr;
203 }
204
205 if (!m_periodic[0] && !m_mirrorPeriodic[0] && (x < m_xMin || x > m_xMax)) {
206 return nullptr;
207 }
208 if (!m_periodic[1] && !m_mirrorPeriodic[1] && (y < m_yMin || x > m_yMax)) {
209 return nullptr;
210 }
211 if (!m_periodic[2] && !m_mirrorPeriodic[2] && (z < m_zMin || x > m_zMax)) {
212 return nullptr;
213 }
214 if (m_active.empty()) return m_medium;
215
216 bool mirrored = false;
217 const double xx =
218 Reduce(x, m_xMin, m_xMax, m_periodic[0], m_mirrorPeriodic[0], mirrored);
219 const double yy =
220 Reduce(y, m_yMin, m_yMax, m_periodic[1], m_mirrorPeriodic[1], mirrored);
221 const double zz =
222 Reduce(z, m_zMin, m_zMax, m_periodic[2], m_mirrorPeriodic[2], mirrored);
223 // Get the indices.
224 const double sx = (xx - m_xMin) / m_dx;
225 const double sy = (yy - m_yMin) / m_dy;
226 const double sz = (zz - m_zMin) / m_dz;
227 const unsigned int i0 = static_cast<unsigned int>(std::floor(sx));
228 const unsigned int j0 = static_cast<unsigned int>(std::floor(sy));
229 const unsigned int k0 = static_cast<unsigned int>(std::floor(sz));
230 const unsigned int i1 = std::min(i0 + 1, m_nX - 1);
231 const unsigned int j1 = std::min(j0 + 1, m_nY - 1);
232 const unsigned int k1 = std::min(k0 + 1, m_nZ - 1);
233 if (m_active[i0][j0][k0] && m_active[i0][j0][k1] &&
234 m_active[i0][j1][k0] && m_active[i0][j1][k1] &&
235 m_active[i1][j0][k0] && m_active[i1][j0][k1] &&
236 m_active[i1][j1][k0] && m_active[i1][j1][k1]) {
237 return m_medium;
238 }
239 return nullptr;
240}
241
242bool ComponentGrid::SetMesh(const unsigned int nx, const unsigned int ny,
243 const unsigned int nz, const double xmin,
244 const double xmax, const double ymin,
245 const double ymax, const double zmin,
246 const double zmax) {
247 Reset();
248 if (nx == 0 || ny == 0 || nz == 0) {
249 std::cerr << m_className << "::SetMesh:\n"
250 << " Number of mesh elements must be positive.\n";
251 return false;
252 }
253 if (xmin >= xmax) {
254 std::cerr << m_className << "::SetMesh: Invalid x range.\n";
255 return false;
256 } else if (ymin >= ymax) {
257 std::cerr << m_className << "::SetMesh: Invalid y range.\n";
258 return false;
259 } else if (zmin >= zmax) {
260 std::cerr << m_className << "::SetMesh: Invalid z range.\n";
261 return false;
262 }
263 m_nX = nx;
264 m_nY = ny;
265 m_nZ = nz;
266 m_xMin = xmin;
267 m_yMin = ymin;
268 m_zMin = zmin;
269 m_xMax = xmax;
270 m_yMax = ymax;
271 m_zMax = zmax;
272 m_dx = m_nX > 1 ? (m_xMax - m_xMin) / (m_nX - 1) : (m_xMax - m_xMin);
273 m_dy = m_nY > 1 ? (m_yMax - m_yMin) / (m_nY - 1) : (m_yMax - m_yMin);
274 m_dz = m_nZ > 1 ? (m_zMax - m_zMin) / (m_nZ - 1) : (m_zMax - m_zMin);
275 m_hasMesh = true;
276 return true;
277}
278
279bool ComponentGrid::GetMesh(
280 unsigned int& nx, unsigned int& ny, unsigned int& nz,
281 double& xmin, double& xmax, double& ymin, double& ymax,
282 double& zmin, double& zmax) const {
283
284 if (!m_hasMesh) return false;
285 nx = m_nX;
286 ny = m_nY;
287 nz = m_nZ;
288 xmin = m_xMin;
289 ymin = m_yMin;
290 zmin = m_zMin;
291 xmax = m_xMax;
292 ymax = m_yMax;
293 zmax = m_zMax;
294 return true;
295}
296
297bool ComponentGrid::LoadElectricField(const std::string& fname,
298 const std::string& fmt, const bool withP,
299 const bool withFlag, const double scaleX,
300 const double scaleE,
301 const double scaleP) {
302 m_ready = false;
303 m_hasPotential = m_hasEfield = false;
304 m_active.assign(m_nX, std::vector<std::vector<bool> >(m_nY, std::vector<bool>(m_nZ, true)));
305 // Read the file.
306 m_pMin = withP ? +1. : 0.;
307 m_pMax = withP ? -1. : 0.;
308 if (!LoadData(fname, fmt, withP, withFlag,
309 scaleX, scaleE, scaleP, m_efields)) {
310 return false;
311 }
312 m_hasEfield = true;
313 m_ready = true;
314 if (withP) m_hasPotential = true;
315 return true;
316}
317
318bool ComponentGrid::LoadWeightingField(const std::string& fname,
319 const std::string& fmt,
320 const bool withP,
321 const double scaleX,
322 const double scaleE,
323 const double scaleP) {
324 m_hasWfield = false;
325 // Read the file.
326 if (!LoadData(fname, fmt, withP, false, scaleX, scaleE, scaleP, m_wfields)) {
327 return false;
328 }
329 m_hasWfield = true;
330 return true;
331}
332bool ComponentGrid::LoadWeightingField(const std::string& fname,
333 const std::string& fmt,
334 const double t, const bool withP,
335 const double scaleX, const double scaleE,
336 const double scaleP) {
337
338 std::vector<std::vector<std::vector<Node> > > wfield;
339 // Read the file.
340 if (!LoadData(fname, fmt, withP, false, scaleX, scaleE, scaleP, wfield)) {
341 return false;
342 }
343 if (m_wdtimes.empty() || t > m_wdtimes.back()) {
344 m_wdtimes.push_back(t);
345 m_wdfields.push_back(std::move(wfield));
346 } else {
347 const auto it = std::upper_bound(m_wdtimes.cbegin(), m_wdtimes.cend(), t);
348 const auto n = std::distance(m_wdtimes.cbegin(), it);
349 m_wdtimes.insert(it, t);
350 m_wdfields.insert(m_wdfields.cbegin() + n, std::move(wfield));
351 }
352 return true;
353}
354
355bool ComponentGrid::LoadMagneticField(const std::string& fname,
356 const std::string& fmt,
357 const double scaleX,
358 const double scaleB) {
359 m_hasBfield = false;
360 // Read the file.
361 if (!LoadData(fname, fmt, false, false, scaleX, scaleB, 1., m_bfields)) {
362 return false;
363 }
364 m_hasBfield = true;
365 return true;
366}
367
368bool ComponentGrid::SaveElectricField(ComponentBase* cmp,
369 const std::string& filename, const std::string& format) {
370
371 if (!cmp) {
372 std::cerr << m_className << "::SaveElectricField: Null pointer.\n";
373 return false;
374 }
375 if (!m_hasMesh) {
376 std::cerr << m_className << "::SaveElectricField: Mesh not set.\n";
377 return false;
378 }
379 const unsigned int fmt = GetFormat(format);
380 if (fmt == 0) {
381 std::cerr << m_className << "::SaveElectricField:\n"
382 << " Unknown format (" << format << ").\n";
383 return false;
384 }
385 std::ofstream outfile;
386 outfile.open(filename.c_str(), std::ios::out);
387 if (!outfile) {
388 std::cerr << m_className << "::SaveElectricField:\n"
389 << " Could not open file " << filename << ".\n";
390 return false;
391 }
392 std::cout << m_className << "::SaveElectricField:\n"
393 << " Exporting field/potential to " << filename << ".\n"
394 << " Be patient...\n";
395 PrintProgress(0.);
396 outfile << "# XMIN = " << m_xMin << ", XMAX = " << m_xMax << ", NX = "
397 << m_nX << "\n";
398 outfile << "# YMIN = " << m_yMin << ", YMAX = " << m_yMax << ", NY = "
399 << m_nY << "\n";
400 outfile << "# ZMIN = " << m_zMin << ", ZMAX = " << m_zMax << ", NZ = "
401 << m_nZ << "\n";
402
403 const unsigned int nValues = m_nX * m_nY * m_nZ;
404 const unsigned int nPrint = std::pow(10, static_cast<unsigned int>(
405 std::max(std::floor(std::log10(nValues)) - 1, 1.)));
406 unsigned int nLines = 0;
407 Medium* medium = nullptr;
408 int status = 0;
409 for (unsigned int i = 0; i < m_nX; ++i) {
410 const double x = m_xMin + i * m_dx;
411 for (unsigned int j = 0; j < m_nY; ++j) {
412 const double y = m_yMin + j * m_dy;
413 for (unsigned int k = 0; k < m_nZ; ++k) {
414 const double z = m_zMin + k * m_dz;
415 if (fmt == 1) {
416 outfile << x << " " << y << " ";
417 } else if (fmt == 2) {
418 outfile << x << " " << y << " " << z << " ";
419 } else if (fmt == 3) {
420 outfile << i << " " << j << " ";
421 } else if (fmt == 4) {
422 outfile << i << " " << j << " " << k << " ";
423 } else if (fmt == 5) {
424 outfile << y << " " << x << " " << z << " ";
425 }
426 double ex = 0., ey = 0., ez = 0., v = 0.;
427 cmp->ElectricField(x, y, z, ex, ey, ez, v, medium, status);
428 outfile << ex << " " << ey << " " << ez << " " << v << "\n";
429 ++nLines;
430 if (nLines % nPrint == 0) PrintProgress(double(nLines) / nValues);
431 }
432 }
433 }
434 outfile.close();
435 std::cout << std::endl << m_className << "::SaveElectricField: Done.\n";
436 return true;
437}
438
439bool ComponentGrid::SaveWeightingField(ComponentBase* cmp,
440 const std::string& id, const std::string& filename,
441 const std::string& format) {
442
443 if (!cmp) {
444 std::cerr << m_className << "::SaveWeightingField: Null pointer.\n";
445 return false;
446 }
447 if (!m_hasMesh) {
448 std::cerr << m_className << "::SaveWeightingField: Mesh not set.\n";
449 return false;
450 }
451 const unsigned int fmt = GetFormat(format);
452 if (fmt == 0) {
453 std::cerr << m_className << "::SaveWeightingField:\n"
454 << " Unknown format (" << format << ").\n";
455 return false;
456 }
457 std::ofstream outfile;
458 outfile.open(filename.c_str(), std::ios::out);
459 if (!outfile) {
460 std::cerr << m_className << "::SaveWeightingField:\n"
461 << " Could not open file " << filename << ".\n";
462 return false;
463 }
464 std::cout << m_className << "::SaveWeightingField:\n"
465 << " Exporting field/potential to " << filename << ".\n"
466 << " Be patient...\n";
467 PrintProgress(0.);
468 outfile << "# XMIN = " << m_xMin << ", XMAX = " << m_xMax << ", NX = "
469 << m_nX << "\n";
470 outfile << "# YMIN = " << m_yMin << ", YMAX = " << m_yMax << ", NY = "
471 << m_nY << "\n";
472 outfile << "# ZMIN = " << m_zMin << ", ZMAX = " << m_zMax << ", NZ = "
473 << m_nZ << "\n";
474 const unsigned int nValues = m_nX * m_nY * m_nZ;
475 const unsigned int nPrint = std::pow(10, static_cast<unsigned int>(
476 std::max(std::floor(std::log10(nValues)) - 1, 1.)));
477 unsigned int nLines = 0;
478 for (unsigned int i = 0; i < m_nX; ++i) {
479 const double x = m_xMin + i * m_dx;
480 for (unsigned int j = 0; j < m_nY; ++j) {
481 const double y = m_yMin + j * m_dy;
482 for (unsigned int k = 0; k < m_nZ; ++k) {
483 const double z = m_zMin + k * m_dz;
484 if (fmt == 1) {
485 outfile << x << " " << y << " ";
486 } else if (fmt == 2) {
487 outfile << x << " " << y << " " << z << " ";
488 } else if (fmt == 3) {
489 outfile << i << " " << j << " ";
490 } else if (fmt == 4) {
491 outfile << i << " " << j << " " << k << " ";
492 } else if (fmt == 5) {
493 outfile << y << " " << x << " " << z << " ";
494 }
495 double wx = 0., wy = 0., wz = 0.;
496 cmp->WeightingField(x, y, z, wx, wy, wz, id);
497 const double v = cmp->WeightingPotential(x, y, z, id);
498 outfile << wx << " " << wy << " " << wz << " " << v << "\n";
499 ++nLines;
500 if (nLines % nPrint == 0) PrintProgress(double(nLines) / nValues);
501 }
502 }
503 }
504 outfile.close();
505 std::cout << std::endl << m_className << "::SaveWeightingField: Done.\n";
506 return true;
507}
508
509bool ComponentGrid::LoadMesh(const std::string& filename, std::string format,
510 const double scaleX) {
511
512 const unsigned int fmt = GetFormat(format);
513 if (fmt == 0) {
514 std::cerr << m_className << "::LoadMesh:\n"
515 << " Unknown format (" << format << ").\n";
516 return false;
517 }
518
519 // Keep track of which mesh parameters we have found.
520 std::bitset<9> found;
521 found.reset();
522 double xmin = 0., ymin = 0., zmin = 0.;
523 double xmax = 0., ymax = 0., zmax = 0.;
524 unsigned int nx = 0, ny = 0, nz = 0;
525
526 // Parse the comment lines in the file.
527 std::ifstream infile;
528 infile.open(filename.c_str(), std::ios::in);
529 if (!infile) {
530 std::cerr << m_className << "::LoadMesh:\n"
531 << " Could not open file " << filename << ".\n";
532 return false;
533 }
534 std::string line;
535 unsigned int nLines = 0;
536 while (!infile.fail()) {
537 // Read one line.
538 std::getline(infile, line);
539 ++nLines;
540 // Strip white space from the beginning of the line.
541 ltrim(line);
542 // Skip empty lines.
543 if (line.empty()) continue;
544 // Skip lines that are not comments.
545 if (line[0] != '#' && !(line[0] == '/' && line[1] == '/')) {
546 continue;
547 }
548 std::size_t pos0 = 0;
549 std::size_t pos1 = line.find("=", pos0);
550 while (pos1 != std::string::npos) {
551 std::string key = line.substr(pos0, pos1 - pos0);
552 std::transform(key.begin(), key.end(), key.begin(), toupper);
553 const std::size_t pos2 = line.find_first_of(",;", pos1 + 1);
554 std::istringstream val(line.substr(pos1 + 1, pos2 - pos1 - 1));
555 if (key.find("XMIN") != std::string::npos) {
556 val >> xmin;
557 found.set(0);
558 } else if (key.find("YMIN") != std::string::npos) {
559 val >> ymin;
560 found.set(1);
561 } else if (key.find("ZMIN") != std::string::npos) {
562 val >> zmin;
563 found.set(2);
564 } else if (key.find("XMAX") != std::string::npos) {
565 val >> xmax;
566 found.set(3);
567 } else if (key.find("YMAX") != std::string::npos) {
568 val >> ymax;
569 found.set(4);
570 } else if (key.find("ZMAX") != std::string::npos) {
571 val >> zmax;
572 found.set(5);
573 } else if (key.find("NX") != std::string::npos) {
574 val >> nx;
575 found.set(6);
576 } else if (key.find("NY") != std::string::npos) {
577 val >> ny;
578 found.set(7);
579 } else if (key.find("NZ") != std::string::npos) {
580 val >> nz;
581 found.set(8);
582 }
583 if (pos2 == std::string::npos) break;
584 pos0 = pos2 + 1;
585 pos1 = line.find("=", pos0);
586 }
587 }
588 infile.close();
589
590 if (fmt == 1 || fmt == 3) {
591 // Try to complement missing information on the z-range.
592 if (!found[8]) {
593 nz = 1;
594 found.set(8);
595 }
596 if (!found[2]) {
597 if (found[0] || found[1] || found[3] || found[4] || found[5]) {
598 zmin = -std::max({fabs(xmin), fabs(xmax),
599 fabs(ymin), fabs(ymax), fabs(zmax)});
600 } else {
601 zmin = -100.;
602 }
603 found.set(2);
604 }
605 if (!found[5]) {
606 zmax = std::max({fabs(xmin), fabs(xmax),
607 fabs(ymin), fabs(ymax), fabs(zmin)});
608 found.set(5);
609 }
610 }
611 if (found.all()) {
612 std::cout << m_className << "::LoadMesh:\n";
613 std::printf("%12.6f < x [cm] < %12.6f, %5u points\n", xmin, xmax, nx);
614 std::printf("%12.6f < y [cm] < %12.6f, %5u points\n", ymin, ymax, ny);
615 std::printf("%12.6f < z [cm] < %12.6f, %5u points\n", zmin, zmax, nz);
616 return SetMesh(nx, ny, nz, xmin, xmax, ymin, ymax, zmin, zmax);
617 }
618
619 if ((fmt == 3 || fmt == 4) && !(found[0] && found[3])) {
620 std::cerr << m_className << "::LoadMesh: x-limits not found.\n";
621 return false;
622 } else if ((fmt == 3 || fmt == 4) && !(found[1] && found[4])) {
623 std::cerr << m_className << "::LoadMesh: y-limits not found.\n";
624 return false;
625 } else if (fmt == 4 && !(found[2] && found[5])) {
626 std::cerr << m_className << "::LoadMesh: z-limits not found.\n";
627 return false;
628 }
629
630 unsigned int nValues = 0;
631 infile.open(filename.c_str(), std::ios::in);
632 if (!infile) {
633 std::cerr << m_className << "::LoadMesh:\n"
634 << " Could not open file " << filename << ".\n";
635 return false;
636 }
637
638 if (!found[0]) xmin = std::numeric_limits<double>::max();
639 if (!found[1]) ymin = std::numeric_limits<double>::max();
640 if (!found[2]) zmin = std::numeric_limits<double>::max();
641 if (!found[3]) xmax = std::numeric_limits<double>::min();
642 if (!found[4]) ymax = std::numeric_limits<double>::min();
643 if (!found[5]) zmax = std::numeric_limits<double>::min();
644 constexpr double tol = 1.e-10;
645 auto cmp = [](double x, double y) {
646 return x < y - tol * (std::fabs(x) + std::fabs(y));
647 };
648 std::set<double, decltype(cmp)> xLines(cmp);
649 std::set<double, decltype(cmp)> yLines(cmp);
650 std::set<double, decltype(cmp)> zLines(cmp);
651 nLines = 0;
652 bool bad = false;
653 while (!infile.fail()) {
654 // Read one line.
655 std::getline(infile, line);
656 ++nLines;
657 // Strip white space from the beginning of the line.
658 ltrim(line);
659 // Skip empty lines.
660 if (line.empty()) continue;
661 // Skip comments.
662 if (line[0] == '#') continue;
663 if (line[0] == '/' && line[1] == '/') continue;
664 std::istringstream data;
665 data.str(line);
666 if (fmt == 1) {
667 // "XY"
668 double x, y;
669 data >> x >> y;
670 if (data.fail()) {
671 PrintError(m_className + "::LoadMesh", nLines, "coordinates");
672 bad = true;
673 break;
674 }
675 x *= scaleX;
676 y *= scaleX;
677 if (!found[0]) xmin = std::min(x, xmin);
678 if (!found[1]) ymin = std::min(y, ymin);
679 if (!found[3]) xmax = std::max(x, xmin);
680 if (!found[4]) ymax = std::max(y, ymin);
681 xLines.insert(x);
682 yLines.insert(y);
683 } else if (fmt == 2) {
684 // "XYZ"
685 double x, y, z;
686 data >> x >> y >> z;
687 if (data.fail()) {
688 PrintError(m_className + "::LoadMesh", nLines, "coordinates");
689 bad = true;
690 break;
691 }
692 x *= scaleX;
693 y *= scaleX;
694 z *= scaleX;
695 if (!found[0]) xmin = std::min(x, xmin);
696 if (!found[1]) ymin = std::min(y, ymin);
697 if (!found[2]) zmin = std::min(z, zmin);
698 if (!found[3]) xmax = std::max(x, xmax);
699 if (!found[4]) ymax = std::max(y, ymax);
700 if (!found[5]) zmax = std::max(z, zmax);
701 xLines.insert(x);
702 yLines.insert(y);
703 zLines.insert(z);
704 } else if (fmt == 3) {
705 // "IJ"
706 unsigned int i = 0, j = 0;
707 data >> i >> j;
708 if (data.fail()) {
709 PrintError(m_className + "::LoadMesh", nLines, "indices");
710 bad = true;
711 break;
712 }
713 if (!found[6]) nx = std::max(nx, i);
714 if (!found[7]) ny = std::max(ny, j);
715 } else if (fmt == 4) {
716 // "IJK"
717 unsigned int i = 0, j = 0, k = 0;
718 data >> i >> j >> k;
719 if (data.fail()) {
720 PrintError(m_className + "::LoadMesh", nLines, "indices");
721 bad = true;
722 break;
723 }
724 if (!found[6]) nx = std::max(nx, i);
725 if (!found[7]) ny = std::max(ny, j);
726 if (!found[8]) nz = std::max(nz, k);
727 } else if (fmt == 5) {
728 // "YXZ"
729 double x, y, z;
730 data >> y >> x >> z;
731 if (data.fail()) {
732 PrintError(m_className + "::LoadMesh", nLines, "coordinates");
733 bad = true;
734 break;
735 }
736 x *= scaleX;
737 y *= scaleX;
738 z *= scaleX;
739 if (!found[0]) xmin = std::min(x, xmin);
740 if (!found[1]) ymin = std::min(y, ymin);
741 if (!found[2]) zmin = std::min(z, zmin);
742 if (!found[3]) xmax = std::max(x, xmax);
743 if (!found[4]) ymax = std::max(y, ymax);
744 if (!found[5]) zmax = std::max(z, zmax);
745 xLines.insert(x);
746 yLines.insert(y);
747 zLines.insert(z);
748 }
749 ++nValues;
750 }
751 infile.close();
752 if (bad) return false;
753
754 if (fmt == 1 || fmt == 2 || fmt == 5) {
755 if (!found[6]) nx = xLines.size();
756 if (!found[7]) ny = yLines.size();
757 if (!found[8]) nz = zLines.size();
758 }
759
760 std::cout << m_className << "::LoadMesh:\n";
761 std::printf("%12.6f < x [cm] < %12.6f, %5u points\n", xmin, xmax, nx);
762 std::printf("%12.6f < y [cm] < %12.6f, %5u points\n", ymin, ymax, ny);
763 std::printf("%12.6f < z [cm] < %12.6f, %5u points\n", zmin, zmax, nz);
764 unsigned int nExpected = nx * ny;
765 if (fmt == 2 || fmt == 4 || fmt == 5) nExpected *= nz;
766 if (nExpected != nValues) {
767 std::cerr << m_className << "::LoadMesh:\n"
768 << " Warning: Expected " << nExpected << " lines, read "
769 << nValues << ".\n";
770 }
771 return SetMesh(nx, ny, nz, xmin, xmax, ymin, ymax, zmin, zmax);
772}
773
774bool ComponentGrid::LoadData(const std::string& filename, std::string format,
775 const bool withPotential, const bool withFlag,
776 const double scaleX, const double scaleF, const double scaleP,
777 std::vector<std::vector<std::vector<Node> > >& fields) {
778
779 if (!m_hasMesh) {
780 if (!LoadMesh(filename, format, scaleX)) {
781 std::cerr << m_className << "::LoadData: Mesh not set.\n";
782 return false;
783 }
784 }
785
786 const unsigned int fmt = GetFormat(format);
787 if (fmt == 0) {
788 std::cerr << m_className << "::LoadData:\n"
789 << " Unknown format (" << format << ").\n";
790 return false;
791 }
792
793 // Set up the grid.
794 Initialise(fields);
795
796 unsigned int nValues = 0;
797 // Keep track of which elements have been read.
798 std::vector<std::vector<std::vector<bool> > > isSet(
799 m_nX,
800 std::vector<std::vector<bool> >(m_nY, std::vector<bool>(m_nZ, false)));
801
802 std::ifstream infile;
803 infile.open(filename.c_str(), std::ios::in);
804 if (!infile) {
805 std::cerr << m_className << "::LoadData:\n"
806 << " Could not open file " << filename << ".\n";
807 return false;
808 }
809
810 std::string line;
811 unsigned int nLines = 0;
812 bool bad = false;
813 while (!infile.fail()) {
814 // Read one line.
815 std::getline(infile, line);
816 ++nLines;
817 // Strip white space from beginning of line.
818 ltrim(line);
819 // Skip empty lines.
820 if (line.empty()) continue;
821 // Skip comments.
822 if (line[0] == '#') continue;
823 if (line[0] == '/' && line[1] == '/') continue;
824 unsigned int i = 0;
825 unsigned int j = 0;
826 unsigned int k = 0;
827 double fx = 0.;
828 double fy = 0.;
829 double fz = 0.;
830 double p = 0.;
831 std::istringstream data;
832 data.str(line);
833 if (fmt == 1) {
834 // "XY"
835 double x, y;
836 data >> x >> y;
837 if (data.fail()) {
838 PrintError(m_className + "::LoadData", nLines, "coordinates");
839 bad = true;
840 break;
841 }
842 x *= scaleX;
843 y *= scaleX;
844 const double u = std::round((x - m_xMin) / m_dx);
845 const double v = std::round((y - m_yMin) / m_dy);
846 i = u < 0. ? 0 : static_cast<unsigned int>(u);
847 j = v < 0. ? 0 : static_cast<unsigned int>(v);
848 if (i >= m_nX) i = m_nX - 1;
849 if (j >= m_nY) j = m_nY - 1;
850 } else if (fmt == 2) {
851 // "XYZ"
852 double x, y, z;
853 data >> x >> y >> z;
854 if (data.fail()) {
855 PrintError(m_className + "::LoadData", nLines, "coordinates");
856 bad = true;
857 break;
858 }
859 x *= scaleX;
860 y *= scaleX;
861 z *= scaleX;
862 const double u = std::round((x - m_xMin) / m_dx);
863 const double v = std::round((y - m_yMin) / m_dy);
864 const double w = std::round((z - m_zMin) / m_dz);
865 i = u < 0. ? 0 : static_cast<unsigned int>(u);
866 j = v < 0. ? 0 : static_cast<unsigned int>(v);
867 j = w < 0. ? 0 : static_cast<unsigned int>(w);
868 if (i >= m_nX) i = m_nX - 1;
869 if (j >= m_nY) j = m_nY - 1;
870 if (k >= m_nZ) k = m_nZ - 1;
871 } else if (fmt == 3) {
872 // "IJ"
873 data >> i >> j;
874 if (data.fail()) {
875 PrintError(m_className + "::LoadData", nLines, "indices");
876 bad = true;
877 break;
878 }
879 } else if (fmt == 4) {
880 // "IJK"
881 data >> i >> j >> k;
882 if (data.fail()) {
883 PrintError(m_className + "::LoadData", nLines, "indices");
884 bad = true;
885 break;
886 }
887 } else if (fmt == 5) {
888 // "YXZ"
889 double x, y, z;
890 data >> y >> x >> z;
891 if (data.fail()) {
892 PrintError(m_className + "::LoadData", nLines, "coordinates");
893 bad = true;
894 break;
895 }
896 x *= scaleX;
897 y *= scaleX;
898 z *= scaleX;
899 const double u = std::round((x - m_xMin) / m_dx);
900 const double v = std::round((y - m_yMin) / m_dy);
901 const double w = std::round((z - m_zMin) / m_dz);
902 i = u < 0. ? 0 : static_cast<unsigned int>(u);
903 j = v < 0. ? 0 : static_cast<unsigned int>(v);
904 j = w < 0. ? 0 : static_cast<unsigned int>(w);
905 if (i >= m_nX) i = m_nX - 1;
906 if (j >= m_nY) j = m_nY - 1;
907 if (k >= m_nZ) k = m_nZ - 1;
908 }
909 // Check the indices.
910 if (i >= m_nX || j >= m_nY || k >= m_nZ) {
911 std::cerr << m_className << "::LoadData:\n"
912 << " Error reading line " << nLines << ".\n"
913 << " Index (" << i << ", " << j << ", " << k
914 << ") out of range.\n";
915 continue;
916 }
917 if (isSet[i][j][k]) {
918 std::cerr << m_className << "::LoadData:\n"
919 << " Error reading line " << nLines << ".\n"
920 << " Mesh element (" << i << ", " << j << ", " << k
921 << ") has already been set.\n";
922 continue;
923 }
924 // Get the field values.
925 if (fmt == 1 || fmt == 3) {
926 // Two-dimensional field-map
927 fz = 0.;
928 data >> fx >> fy;
929 } else if (fmt == 5) {
930 data >> fy >> fx >> fz;
931 } else {
932 data >> fx >> fy >> fz;
933 }
934 if (data.fail()) {
935 PrintError(m_className + "::LoadData", nLines, "field components");
936 bad = true;
937 break;
938 }
939 fx *= scaleF;
940 fy *= scaleF;
941 fz *= scaleF;
942 if (withPotential) {
943 data >> p;
944 if (data.fail()) {
945 PrintError(m_className + "::LoadData", nLines, "potential");
946 bad = true;
947 break;
948 }
949 p *= scaleP;
950 if (m_pMin > m_pMax) {
951 // First value.
952 m_pMin = p;
953 m_pMax = p;
954 } else {
955 if (p < m_pMin) m_pMin = p;
956 if (p > m_pMax) m_pMax = p;
957 }
958 }
959 int flag = 0;
960 if (withFlag) {
961 data >> flag;
962 if (data.fail()) {
963 PrintError(m_className + "::LoadData", nLines, "region");
964 bad = true;
965 break;
966 }
967 }
968 const bool isActive = flag == 0 ? false : true;
969 if (fmt == 1 || fmt == 3) {
970 // Two-dimensional field-map
971 for (unsigned int kk = 0; kk < m_nZ; ++kk) {
972 fields[i][j][kk].fx = fx;
973 fields[i][j][kk].fy = fy;
974 fields[i][j][kk].fz = fz;
975 fields[i][j][kk].v = p;
976 if (withFlag) m_active[i][j][kk] = isActive;
977 isSet[i][j][kk] = true;
978 }
979 } else {
980 fields[i][j][k].fx = fx;
981 fields[i][j][k].fy = fy;
982 fields[i][j][k].fz = fz;
983 fields[i][j][k].v = p;
984 isSet[i][j][k] = true;
985 }
986 ++nValues;
987 }
988 infile.close();
989 if (bad) return false;
990 std::cout << m_className << "::LoadData:\n"
991 << " Read " << nValues << " values from " << filename << ".\n";
992 unsigned int nExpected = m_nX * m_nY;
993 if (fmt == 2 || fmt == 4 || fmt == 5) nExpected *= m_nZ;
994 if (nExpected != nValues) {
995 std::cerr << m_className << "::LoadData:\n"
996 << " Expected " << nExpected << " values.\n";
997 }
998 return true;
999}
1000
1001bool ComponentGrid::GetBoundingBox(double& xmin, double& ymin, double& zmin,
1002 double& xmax, double& ymax, double& zmax) {
1003 if (!m_ready) return false;
1004 if (m_periodic[0] || m_mirrorPeriodic[0]) {
1005 xmin = -INFINITY;
1006 xmax = +INFINITY;
1007 } else {
1008 xmin = m_xMin;
1009 xmax = m_xMax;
1010 }
1011
1012 if (m_periodic[1] || m_mirrorPeriodic[1]) {
1013 ymin = -INFINITY;
1014 ymax = +INFINITY;
1015 } else {
1016 ymin = m_yMin;
1017 ymax = m_yMax;
1018 }
1019
1020 if (m_periodic[2] || m_mirrorPeriodic[2]) {
1021 zmin = -INFINITY;
1022 zmax = +INFINITY;
1023 } else {
1024 zmin = m_zMin;
1025 zmax = m_zMax;
1026 }
1027 return true;
1028}
1029
1030bool ComponentGrid::GetVoltageRange(double& vmin, double& vmax) {
1031 if (!m_ready) return false;
1032 vmin = m_pMin;
1033 vmax = m_pMax;
1034 return true;
1035}
1036
1037bool ComponentGrid::GetElectricFieldRange(double& exmin, double& exmax,
1038 double& eymin, double& eymax,
1039 double& ezmin, double& ezmax) {
1040 if (!m_ready) {
1041 PrintNotReady(m_className + "::GetElectricFieldRange");
1042 return false;
1043 }
1044
1045 exmin = exmax = m_efields[0][0][0].fx;
1046 eymin = eymax = m_efields[0][0][0].fy;
1047 ezmin = ezmax = m_efields[0][0][0].fz;
1048 for (unsigned int i = 0; i < m_nX; ++i) {
1049 for (unsigned int j = 0; j < m_nY; ++j) {
1050 for (unsigned int k = 0; k < m_nZ; ++k) {
1051 const Node& node = m_efields[i][j][k];
1052 if (node.fx < exmin) exmin = node.fx;
1053 if (node.fx > exmax) exmax = node.fx;
1054 if (node.fy < eymin) eymin = node.fy;
1055 if (node.fy > eymax) eymax = node.fy;
1056 if (node.fz < ezmin) ezmin = node.fz;
1057 if (node.fz > ezmax) ezmax = node.fz;
1058 }
1059 }
1060 }
1061 return true;
1062}
1063
1064void ComponentGrid::SetMedium(Medium* m) {
1065 if (!m) {
1066 std::cerr << m_className << "::SetMedium: Null pointer.\n";
1067 }
1068 m_medium = m;
1069}
1070
1071bool ComponentGrid::GetField(
1072 const double xi, const double yi, const double zi,
1073 const std::vector<std::vector<std::vector<Node> > >& field, double& fx,
1074 double& fy, double& fz, double& p, bool& active) {
1075 if (!m_hasMesh) {
1076 std::cerr << m_className << "::GetField: Mesh is not set.\n";
1077 return false;
1078 }
1079
1080 // Reduce the point to the basic cell (in case of periodicity) and
1081 // check if it is inside the mesh.
1082 bool xMirrored = false;
1083 const double x =
1084 Reduce(xi, m_xMin, m_xMax, m_periodic[0], m_mirrorPeriodic[0], xMirrored);
1085 if (x < m_xMin || x > m_xMax) return false;
1086 bool yMirrored = false;
1087 const double y =
1088 Reduce(yi, m_yMin, m_yMax, m_periodic[1], m_mirrorPeriodic[1], yMirrored);
1089 if (y < m_yMin || y > m_yMax) return false;
1090 bool zMirrored = false;
1091 const double z =
1092 Reduce(zi, m_zMin, m_zMax, m_periodic[2], m_mirrorPeriodic[2], zMirrored);
1093 if (z < m_zMin || z > m_zMax) return false;
1094
1095 // Get the indices.
1096 const double sx = (x - m_xMin) / m_dx;
1097 const double sy = (y - m_yMin) / m_dy;
1098 const double sz = (z - m_zMin) / m_dz;
1099 const unsigned int i0 = static_cast<unsigned int>(std::floor(sx));
1100 const unsigned int j0 = static_cast<unsigned int>(std::floor(sy));
1101 const unsigned int k0 = static_cast<unsigned int>(std::floor(sz));
1102 const double ux = sx - i0;
1103 const double uy = sy - j0;
1104 const double uz = sz - k0;
1105 const unsigned int i1 = std::min(i0 + 1, m_nX - 1);
1106 const unsigned int j1 = std::min(j0 + 1, m_nY - 1);
1107 const unsigned int k1 = std::min(k0 + 1, m_nZ - 1);
1108 const double vx = 1. - ux;
1109 const double vy = 1. - uy;
1110 const double vz = 1. - uz;
1111 if (!m_active.empty()) {
1112 active = m_active[i0][j0][k0] && m_active[i0][j0][k1] &&
1113 m_active[i0][j1][k0] && m_active[i0][j1][k1] &&
1114 m_active[i1][j0][k0] && m_active[i1][j0][k1] &&
1115 m_active[i1][j1][k0] && m_active[i1][j1][k1];
1116 }
1117 const Node& n000 = field[i0][j0][k0];
1118 const Node& n100 = field[i1][j0][k0];
1119 const Node& n010 = field[i0][j1][k0];
1120 const Node& n110 = field[i1][j1][k0];
1121 const Node& n001 = field[i0][j0][k1];
1122 const Node& n101 = field[i1][j0][k1];
1123 const Node& n011 = field[i0][j1][k1];
1124 const Node& n111 = field[i1][j1][k1];
1125
1126 if (m_debug) {
1127 std::cout << m_className << "::GetField: Determining field at ("
1128 << xi << ", " << yi << ", " << zi << ").\n"
1129 << " X: " << i0 << " (" << ux << ") - "
1130 << i1 << " (" << vx << ").\n"
1131 << " Y: " << j0 << " (" << uy << ") - "
1132 << j1 << " (" << vy << ").\n"
1133 << " Z: " << k0 << " (" << uz << ") - "
1134 << k1 << " (" << vz << ").\n";
1135 }
1136 fx = ((n000.fx * vx + n100.fx * ux) * vy +
1137 (n010.fx * vx + n110.fx * ux) * uy) *
1138 vz +
1139 ((n001.fx * vx + n101.fx * ux) * vy +
1140 (n011.fx * vx + n111.fx * ux) * uy) *
1141 uz;
1142 fy = ((n000.fy * vx + n100.fy * ux) * vy +
1143 (n010.fy * vx + n110.fy * ux) * uy) *
1144 vz +
1145 ((n001.fy * vx + n101.fy * ux) * vy +
1146 (n011.fy * vx + n111.fy * ux) * uy) *
1147 uz;
1148 fz = ((n000.fz * vx + n100.fz * ux) * vy +
1149 (n010.fz * vx + n110.fz * ux) * uy) *
1150 vz +
1151 ((n001.fz * vx + n101.fz * ux) * vy +
1152 (n011.fz * vx + n111.fz * ux) * uy) *
1153 uz;
1154 p = ((n000.v * vx + n100.v * ux) * vy + (n010.v * vx + n110.v * ux) * uy) *
1155 vz +
1156 ((n001.v * vx + n101.v * ux) * vy + (n011.v * vx + n111.v * ux) * uy) *
1157 uz;
1158 if (xMirrored) fx = -fx;
1159 if (yMirrored) fy = -fy;
1160 if (zMirrored) fz = -fz;
1161 return true;
1162}
1163
1164bool ComponentGrid::GetElectricField(
1165 const unsigned int i, const unsigned int j, const unsigned int k,
1166 double& v, double& ex, double& ey, double& ez) const {
1167 v = ex = ey = ez = 0.;
1168 if (!m_ready) {
1169 if (!m_hasMesh) {
1170 std::cerr << m_className << "::GetElectricField: Mesh not set.\n";
1171 return false;
1172 }
1173 PrintNotReady(m_className + "::GetElectricField");
1174 return false;
1175 }
1176 if (i >= m_nX || j >= m_nY || k >= m_nZ) {
1177 std::cerr << m_className << "::GetElectricField: Index out of range.\n";
1178 return false;
1179 }
1180 const Node& node = m_efields[i][j][k];
1181 v = node.v;
1182 ex = node.fx;
1183 ey = node.fy;
1184 ez = node.fz;
1185 return true;
1186}
1187
1188void ComponentGrid::Reset() {
1189 m_efields.clear();
1190 m_bfields.clear();
1191 m_wfields.clear();
1192
1193 m_wdfields.clear();
1194 m_wdtimes.clear();
1195
1196 m_active.clear();
1197
1198 m_nX = m_nY = m_nZ = 0;
1199 m_xMin = m_yMin = m_zMin = 0.;
1200 m_xMax = m_yMax = m_zMax = 0.;
1201 m_pMin = m_pMax = 0.;
1202 m_medium = nullptr;
1203
1204 m_hasMesh = false;
1205 m_hasPotential = false;
1206 m_hasEfield = false;
1207 m_hasBfield = false;
1208 m_hasWfield = false;
1209 m_ready = false;
1210
1211 m_wField_xOffset = 0.;
1212 m_wField_yOffset = 0.;
1213 m_wField_zOffset = 0.;
1214}
1215
1216void ComponentGrid::UpdatePeriodicity() {
1217 if (!m_ready) {
1218 PrintNotReady(m_className + "::UpdatePeriodicity");
1219 return;
1220 }
1221
1222 // Check for conflicts.
1223 for (unsigned int i = 0; i < 3; ++i) {
1224 if (m_periodic[i] && m_mirrorPeriodic[i]) {
1225 std::cerr << m_className << "::UpdatePeriodicity:\n"
1226 << " Both simple and mirror periodicity requested. Reset.\n";
1227 m_periodic[i] = m_mirrorPeriodic[i] = false;
1228 }
1229 }
1230
1231 if (m_axiallyPeriodic[0] || m_axiallyPeriodic[1] || m_axiallyPeriodic[2]) {
1232 std::cerr << m_className << "::UpdatePeriodicity:\n"
1233 << " Axial symmetry is not supported. Reset.\n";
1234 m_axiallyPeriodic.fill(false);
1235 }
1236
1237 if (m_rotationSymmetric[0] || m_rotationSymmetric[1] ||
1238 m_rotationSymmetric[2]) {
1239 std::cerr << m_className << "::UpdatePeriodicity:\n"
1240 << " Rotation symmetry is not supported. Reset.\n";
1241 m_rotationSymmetric.fill(false);
1242 }
1243}
1244
1245double ComponentGrid::Reduce(const double xin, const double xmin,
1246 const double xmax, const bool simplePeriodic,
1247 const bool mirrorPeriodic, bool& mirrored) const {
1248 // In case of periodicity, reduce the coordinate to the basic cell.
1249 double x = xin;
1250 const double lx = xmax - xmin;
1251 if (simplePeriodic) {
1252 x = xmin + fmod(x - xmin, lx);
1253 if (x < xmin) x += lx;
1254 } else if (mirrorPeriodic) {
1255 double xNew = xmin + fmod(x - xmin, lx);
1256 if (xNew < xmin) xNew += lx;
1257 const int nx = int(floor(0.5 + (xNew - x) / lx));
1258 if (nx != 2 * (nx / 2)) {
1259 xNew = xmin + xmax - xNew;
1260 mirrored = true;
1261 }
1262 x = xNew;
1263 }
1264 return x;
1265}
1266
1267void ComponentGrid::Initialise(
1268 std::vector<std::vector<std::vector<Node> > >& fields) {
1269
1270 fields.resize(m_nX);
1271 for (unsigned int i = 0; i < m_nX; ++i) {
1272 fields[i].resize(m_nY);
1273 for (unsigned int j = 0; j < m_nY; ++j) {
1274 fields[i][j].resize(m_nZ);
1275 for (unsigned int k = 0; k < m_nZ; ++k) {
1276 fields[i][j][k].fx = 0.;
1277 fields[i][j][k].fy = 0.;
1278 fields[i][j][k].fz = 0.;
1279 fields[i][j][k].v = 0.;
1280 }
1281 }
1282 }
1283}
1284
1285}
Garfield::ComponentBase
Abstract base class for components.
Definition: ComponentBase.hh:13
Garfield::ComponentBase::m_periodic
std::array< bool, 3 > m_periodic
Simple periodicity in x, y, z.
Definition: ComponentBase.hh:299
Garfield::ComponentBase::ElectricField
virtual void ElectricField(const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status)=0
Garfield::ComponentBase::m_rotationSymmetric
std::array< bool, 3 > m_rotationSymmetric
Rotation symmetry around x-axis, y-axis, z-axis.
Definition: ComponentBase.hh:305
Garfield::ComponentBase::WeightingPotential
virtual double WeightingPotential(const double x, const double y, const double z, const std::string &label)
Definition: ComponentBase.cc:54
Garfield::ComponentBase::m_axiallyPeriodic
std::array< bool, 3 > m_axiallyPeriodic
Axial periodicity in x, y, z.
Definition: ComponentBase.hh:303
Garfield::ComponentBase::MagneticField
virtual void MagneticField(const double x, const double y, const double z, double &bx, double &by, double &bz, int &status)
Definition: ComponentBase.cc:64
Garfield::ComponentBase::m_ready
bool m_ready
Ready for use?
Definition: ComponentBase.hh:291
Garfield::ComponentBase::m_mirrorPeriodic
std::array< bool, 3 > m_mirrorPeriodic
Mirror periodicity in x, y, z.
Definition: ComponentBase.hh:301
Garfield::ComponentBase::WeightingField
virtual void WeightingField(const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label)
Definition: ComponentBase.cc:32
Garfield::ComponentBase::m_className
std::string m_className
Class name.
Definition: ComponentBase.hh:285
Garfield::ComponentBase::m_debug
bool m_debug
Switch on/off debugging messages.
Definition: ComponentBase.hh:310
Garfield::ComponentGrid::SetMesh
bool SetMesh(const unsigned int nx, const unsigned int ny, const unsigned int nz, const double xmin, const double xmax, const double ymin, const double ymax, const double zmin, const double zmax)
Definition: ComponentGrid.cc:242
Garfield::ComponentGrid::SaveElectricField
bool SaveElectricField(ComponentBase *cmp, const std::string &filename, const std::string &format)
Definition: ComponentGrid.cc:368
Garfield::ComponentGrid::GetBoundingBox
bool GetBoundingBox(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax) override
Get the bounding box coordinates.
Definition: ComponentGrid.cc:1001
Garfield::ComponentGrid::ComponentGrid
ComponentGrid()
Constructor.
Definition: ComponentGrid.cc:67
Garfield::ComponentGrid::LoadElectricField
bool LoadElectricField(const std::string &filename, const std::string &format, const bool withPotential, const bool withFlag, const double scaleX=1., const double scaleE=1., const double scaleP=1.)
Definition: ComponentGrid.cc:297
Garfield::ComponentGrid::DelayedWeightingField
void DelayedWeightingField(const double x, const double y, const double z, const double t, double &wx, double &wy, double &wz, const std::string &label) override
Definition: ComponentGrid.cc:134
Garfield::ComponentGrid::GetMesh
bool GetMesh(unsigned int &nx, unsigned int &ny, unsigned int &nz, double &xmin, double &xmax, double &ymin, double &ymax, double &zmin, double &zmax) const
Retrieve the parameters of the grid.
Definition: ComponentGrid.cc:279
Garfield::ComponentGrid::WeightingPotential
double WeightingPotential(const double x, const double y, const double z, const std::string &label) override
Definition: ComponentGrid.cc:120
Garfield::ComponentGrid::GetVoltageRange
bool GetVoltageRange(double &vmin, double &vmax) override
Calculate the voltage range [V].
Definition: ComponentGrid.cc:1030
Garfield::ComponentGrid::SaveWeightingField
bool SaveWeightingField(ComponentBase *cmp, const std::string &id, const std::string &filename, const std::string &format)
Definition: ComponentGrid.cc:439
Garfield::ComponentGrid::SetMedium
void SetMedium(Medium *m)
Set the medium.
Definition: ComponentGrid.cc:1064
Garfield::ComponentGrid::GetMedium
Medium * GetMedium() const
Get the medium.
Definition: ComponentGrid.hh:126
Garfield::ComponentGrid::LoadMagneticField
bool LoadMagneticField(const std::string &filename, const std::string &format, const double scaleX=1., const double scaleB=1.)
Import magnetic field values from a file.
Definition: ComponentGrid.cc:355
Garfield::ComponentGrid::GetElectricField
bool GetElectricField(const unsigned int i, const unsigned int j, const unsigned int k, double &v, double &ex, double &ey, double &ez) const
Return the field at a given node.
Definition: ComponentGrid.cc:1164
Garfield::ComponentGrid::LoadWeightingField
bool LoadWeightingField(const std::string &filename, const std::string &format, const bool withPotential, const double scaleX=1., const double scaleE=1., const double scaleP=1.)
Import (prompt) weighting field from file.
Definition: ComponentGrid.cc:318
Garfield::ComponentGrid::ElectricField
void ElectricField(const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status) override
Calculate the drift field [V/cm] and potential [V] at (x, y, z).
Definition: ComponentGrid.cc:71
Garfield::ComponentGrid::GetElectricFieldRange
bool GetElectricFieldRange(double &exmin, double &exmax, double &eymin, double &eymax, double &ezmin, double &ezmax)
Definition: ComponentGrid.cc:1037
Garfield::ComponentGrid::MagneticField
void MagneticField(const double x, const double y, const double z, double &bx, double &by, double &bz, int &status) override
Definition: ComponentGrid.cc:182
Garfield::ComponentGrid::SetWeightingFieldOffset
void SetWeightingFieldOffset(const double x, const double y, const double z)
Definition: ComponentGrid.cc:175
Garfield::ComponentGrid::WeightingField
void WeightingField(const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label) override
Definition: ComponentGrid.cc:106
Garfield::Medium
Abstract base class for media.
Definition: Medium.hh:13
Garfield::ltrim
void ltrim(std::string &line)
Definition: Utilities.hh:10
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615