Garfield++ 5.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
ComponentVoxel.cc
Go to the documentation of this file.
1#include <algorithm>
2#include <cmath>
3#include <fstream>
4#include <iostream>
5#include <sstream>
6#include <string>
7
8#include "Garfield/ComponentVoxel.hh"
9#include "Garfield/Utilities.hh"
10#include "Garfield/GarfieldConstants.hh"
11
12namespace Garfield {
13
14ComponentVoxel::ComponentVoxel() : Component("Voxel") {}
15
16void ComponentVoxel::ElectricField(
17 const double x, const double y, const double z,
18 double& ex, double& ey, double& ez, double& p, Medium*& m, int& status) {
19 m = nullptr;
20 status = 0;
21
22 // Make sure the field map has been loaded.
23 if (!m_ready) {
24 std::cerr << m_className << "::ElectricField:\n"
25 << " Field map is not available for interpolation.\n";
26 status = -10;
27 return;
28 }
29
30 status = 0;
31 int region = -1;
32 if (!GetField(x, y, z, m_efields, ex, ey, ez, p, region)) {
33 status = -6;
34 return;
35 }
36
37 if (region < 0 || region > (int)m_media.size()) {
38 m = nullptr;
39 status = -5;
40 return;
41 }
42 m = m_media[region];
43 if (!m) status = -5;
44}
45
46void ComponentVoxel::ElectricField(
47 const double x, const double y, const double z,
48 double& ex, double& ey, double& ez, Medium*& m, int& status) {
49 double v = 0.;
50 ElectricField(x, y, z, ex, ey, ez, v, m, status);
51}
52
53void ComponentVoxel::WeightingField(
54 const double x, const double y, const double z,
55 double& wx, double& wy, double& wz, const std::string& /*label*/) {
56
57 wx = wy = wz = 0.;
58 if (!m_hasWfield) return;
59 const double xx = x - m_wField_xOffset;
60 const double yy = y - m_wField_yOffset;
61 const double zz = z - m_wField_zOffset;
62 double wp = 0.;
63 int region = 0;
64 GetField(xx, yy, zz, m_wfields, wx, wy, wz, wp, region);
65}
66
67double ComponentVoxel::WeightingPotential(
68 const double x, const double y, const double z,
69 const std::string& /*label*/) {
70
71 if (!m_hasWfield) return 0.;
72 const double xx = x - m_wField_xOffset;
73 const double yy = y - m_wField_yOffset;
74 const double zz = z - m_wField_zOffset;
75 double wx = 0., wy = 0., wz = 0.;
76 double wp = 0.;
77 int region = 0;
78 if (!GetField(xx, yy, zz, m_wfields, wx, wy, wz, wp, region)) return 0.;
79 return wp;
80}
81
82void ComponentVoxel::DelayedWeightingField(
83 const double x, const double y, const double z, const double t,
84 double& wx, double& wy, double& wz, const std::string& /*label*/) {
85
86 wx = wy = wz = 0.;
87 if (m_wdtimes.empty()) return;
88 // Assume no weighting field for times outside the range of available maps.
89 if (t < m_wdtimes.front() || t > m_wdtimes.back()) return;
90
91 const double xx = x - m_wField_xOffset;
92 const double yy = y - m_wField_yOffset;
93 const double zz = z - m_wField_zOffset;
94
95 const auto it1 = std::upper_bound(m_wdtimes.cbegin(), m_wdtimes.cend(), t);
96 const auto it0 = std::prev(it1);
97
98 const double dt = t - *it0;
99 double wp = 0.;
100 int region = 0;
101 const unsigned int i0 = it0 - m_wdtimes.cbegin();
102 double wx0 = 0., wy0 = 0., wz0 = 0.;
103 if (!GetField(xx, yy, zz, m_wdfields[i0], wx0, wy0, wz0, wp, region)) {
104 return;
105 }
106 if (dt < Small || it1 == m_wdtimes.cend()) {
107 wx = wx0;
108 wy = wy0;
109 wz = wz0;
110 return;
111 }
112 const unsigned int i1 = it1 - m_wdtimes.cbegin();
113 double wx1 = 0., wy1 = 0., wz1 = 0.;
114 if (!GetField(xx, yy, zz, m_wdfields[i1], wx1, wy1, wz1, wp, region)) {
115 return;
116 }
117 const double f1 = dt / (*it1 - *it0);
118 const double f0 = 1. - f1;
119 wx = f0 * wx0 + f1 * wx1;
120 wy = f0 * wy0 + f1 * wy1;
121 wz = f0 * wz0 + f1 * wz1;
122}
123
124double ComponentVoxel::DelayedWeightingPotential(
125 const double x, const double y, const double z, const double t,
126 const std::string& /*label*/) {
127
128 if (m_wdtimes.empty()) return 0.;
129 // Outside the range of the available maps?
130 if (t < m_wdtimes.front() || t > m_wdtimes.back()) return 0.;
131
132 const double xx = x - m_wField_xOffset;
133 const double yy = y - m_wField_yOffset;
134 const double zz = z - m_wField_zOffset;
135
136 const auto it1 = std::upper_bound(m_wdtimes.cbegin(), m_wdtimes.cend(), t);
137 const auto it0 = std::prev(it1);
138
139 const double dt = t - *it0;
140 int region = 0;
141 const unsigned int i0 = it0 - m_wdtimes.cbegin();
142 double wp0 = 0., wx0 = 0., wy0 = 0., wz0 = 0.;
143 if (!GetField(xx, yy, zz, m_wdfields[i0], wx0, wy0, wz0, wp0, region)) {
144 return 0.;
145 }
146 if (dt < Small || it1 == m_wdtimes.cend()) return 0.;
147
148 const unsigned int i1 = it1 - m_wdtimes.cbegin();
149 double wp1 = 0., wx1 = 0., wy1 = 0., wz1 = 0.;
150 if (!GetField(xx, yy, zz, m_wdfields[i1], wx1, wy1, wz1, wp1, region)) {
151 return 0.;
152 }
153 const double f1 = dt / (*it1 - *it0);
154 const double f0 = 1. - f1;
155 return f0 * wp0 + f1 * wp1;
156}
157
158void ComponentVoxel::SetWeightingFieldOffset(const double x, const double y,
159 const double z) {
160 m_wField_xOffset = x;
161 m_wField_yOffset = y;
162 m_wField_zOffset = z;
163}
164
165void ComponentVoxel::MagneticField(const double x, const double y,
166 const double z, double& bx, double& by,
167 double& bz, int& status) {
168 status = 0;
169 if (!m_hasBfield) {
170 return Component::MagneticField(x, y, z, bx, by, bz, status);
171 }
172
173 int region = -1;
174 double p = 0.;
175 if (!GetField(x, y, z, m_bfields, bx, by, bz, p, region)) {
176 status = -6;
177 }
178}
179
180bool ComponentVoxel::HasMagneticField() const {
181 return m_hasBfield ? true : Component::HasMagneticField();
182}
183
184Medium* ComponentVoxel::GetMedium(const double x, const double y,
185 const double z) {
186 // Make sure the field map has been loaded.
187 if (!m_ready) {
188 std::cerr << m_className << "::GetMedium:\n"
189 << " Field map is not available for interpolation.\n";
190 return nullptr;
191 }
192
193 unsigned int i, j, k;
194 bool xMirrored, yMirrored, zMirrored;
195 if (!GetElement(x, y, z, i, j, k, xMirrored, yMirrored, zMirrored)) {
196 return nullptr;
197 }
198 const int region = m_regions[i][j][k];
199 if (region < 0 || region > (int)m_media.size()) return nullptr;
200 return m_media[region];
201}
202
203void ComponentVoxel::SetMesh(const unsigned int nx, const unsigned int ny,
204 const unsigned int nz, const double xmin,
205 const double xmax, const double ymin,
206 const double ymax, const double zmin,
207 const double zmax) {
208 Reset();
209 if (nx == 0 || ny == 0 || nz == 0) {
210 std::cerr << m_className << "::SetMesh:\n"
211 << " Number of mesh elements must be positive.\n";
212 return;
213 }
214 if (xmin >= xmax) {
215 std::cerr << m_className << "::SetMesh: Invalid x range.\n";
216 return;
217 } else if (ymin >= ymax) {
218 std::cerr << m_className << "::SetMesh: Invalid y range.\n";
219 return;
220 } else if (zmin >= zmax) {
221 std::cerr << m_className << "::SetMesh: Invalid z range.\n";
222 return;
223 }
224 m_nX = nx;
225 m_nY = ny;
226 m_nZ = nz;
227 m_xMin = xmin;
228 m_yMin = ymin;
229 m_zMin = zmin;
230 m_xMax = xmax;
231 m_yMax = ymax;
232 m_zMax = zmax;
233 m_dx = (m_xMax - m_xMin) / m_nX;
234 m_dy = (m_yMax - m_yMin) / m_nY;
235 m_dz = (m_zMax - m_zMin) / m_nZ;
236 m_hasMesh = true;
237}
238
239bool ComponentVoxel::LoadElectricField(const std::string& fname,
240 const std::string& fmt,
241 const bool withP, const bool withR,
242 const double scaleX, const double scaleE,
243 const double scaleP) {
244 m_ready = false;
245 m_efields.clear();
246 m_hasPotential = m_hasEfield = false;
247 if (!m_hasMesh) {
248 std::cerr << m_className << "::LoadElectricField:\n"
249 << " Mesh is not set. Call SetMesh first.\n";
250 return false;
251 }
252
253 // Set up the grid.
254 Initialise(m_efields);
255 InitialiseRegions();
256
257 m_pMin = m_pMax = 0.;
258 if (withP) {
259 m_pMin = 1.;
260 m_pMax = -1.;
261 }
262 if (!LoadData(fname, fmt, withP, withR, scaleX, scaleE, scaleP, m_efields)) {
263 return false;
264 }
265 m_hasEfield = true;
266 m_ready = true;
267 if (withP) m_hasPotential = true;
268 return true;
269}
270
271bool ComponentVoxel::LoadWeightingField(const std::string& fname,
272 const std::string& fmt,
273 const bool withP,
274 const double scaleX,
275 const double scaleE,
276 const double scaleP) {
277 m_hasWfield = false;
278 if (!m_hasMesh) {
279 std::cerr << m_className << "::LoadWeightingField:\n"
280 << " Mesh is not set. Call SetMesh first.\n";
281 return false;
282 }
283
284 // Set up the grid.
285 Initialise(m_wfields);
286 if (m_regions.empty()) InitialiseRegions();
287
288 // Read the file.
289 if (!LoadData(fname, fmt, withP, false, scaleX, scaleE, scaleP, m_wfields)) {
290 return false;
291 }
292 m_hasWfield = true;
293 return true;
294}
295bool ComponentVoxel::LoadWeightingField(const std::string& fname,
296 const std::string& fmt,
297 const double t, const bool withP,
298 const double scaleX, const double scaleE,
299 const double scaleP) {
300
301 if (!m_hasMesh) {
302 std::cerr << m_className << "::LoadWeightingField:\n"
303 << " Mesh is not set. Call SetMesh first.\n";
304 return false;
305 }
306
307 std::vector<std::vector<std::vector<Element> > > wfield;
308 Initialise(wfield);
309 if (m_regions.empty()) InitialiseRegions();
310
311 // Read the file.
312 if (!LoadData(fname, fmt, withP, false, scaleX, scaleE, scaleP, wfield)) {
313 return false;
314 }
315 if (m_wdtimes.empty() || t > m_wdtimes.back()) {
316 m_wdtimes.push_back(t);
317 m_wdfields.push_back(std::move(wfield));
318 } else {
319 const auto it = std::upper_bound(m_wdtimes.begin(), m_wdtimes.end(), t);
320 const auto n = std::distance(m_wdtimes.begin(), it);
321 m_wdtimes.insert(it, t);
322 m_wdfields.insert(m_wdfields.begin() + n, std::move(wfield));
323 }
324 return true;
325}
326
327bool ComponentVoxel::LoadMagneticField(const std::string& fname,
328 const std::string& fmt,
329 const double scaleX,
330 const double scaleB) {
331 m_hasBfield = false;
332 if (!m_hasMesh) {
333 std::cerr << m_className << "::LoadMagneticField:\n"
334 << " Mesh is not set. Call SetMesh first.\n";
335 return false;
336 }
337
338 // Set up the grid.
339 Initialise(m_bfields);
340 InitialiseRegions();
341
342 // Read the file.
343 if (!LoadData(fname, fmt, false, false, scaleX, scaleB, 1., m_bfields)) {
344 return false;
345 }
346 m_hasBfield = true;
347 return true;
348}
349
350bool ComponentVoxel::LoadData(const std::string& filename, std::string format,
351 const bool withPotential, const bool withRegion,
352 const double scaleX, const double scaleF, const double scaleP,
353 std::vector<std::vector<std::vector<Element> > >& fields) {
354
355 if (!m_hasMesh) {
356 std::cerr << m_className << "::LoadData: Mesh has not been set.\n";
357 return false;
358 }
359
360 unsigned int nValues = 0;
361 // Keep track of which elements have been read.
362 std::vector<std::vector<std::vector<bool> > > isSet(
363 m_nX,
364 std::vector<std::vector<bool> >(m_nY, std::vector<bool>(m_nZ, false)));
365
366 std::ifstream infile(filename);
367 if (!infile) {
368 std::cerr << m_className << "::LoadData:\n"
369 << " Could not open file " << filename << ".\n";
370 return false;
371 }
372
373 std::transform(format.begin(), format.end(), format.begin(), toupper);
374 unsigned int fmt = 0;
375 if (format == "XY") {
376 fmt = 1;
377 } else if (format == "XYZ") {
378 fmt = 2;
379 } else if (format == "IJ") {
380 fmt = 3;
381 } else if (format == "IJK") {
382 fmt = 4;
383 } else if (format == "YXZ") {
384 fmt = 5;
385 } else {
386 std::cerr << m_className << "::LoadData:\n"
387 << " Unknown format (" << format << ").\n";
388 return false;
389 }
390 std::string line;
391 unsigned int nLines = 0;
392 bool bad = false;
393 // Read the file line by line.
394 while (std::getline(infile, line)) {
395 ++nLines;
396 // Strip white space from beginning of line.
397 ltrim(line);
398 // Skip empty lines.
399 if (line.empty()) continue;
400 // Skip comments.
401 if (line[0] == '#') continue;
402 if (line[0] == '/' && line[1] == '/') continue;
403 unsigned int i = 0;
404 unsigned int j = 0;
405 unsigned int k = 0;
406 double fx = 0.;
407 double fy = 0.;
408 double fz = 0.;
409 double v = 0.;
410 int region = 0;
411 std::istringstream data(line);
412 if (fmt == 1) {
413 // "XY"
414 double x, y;
415 data >> x >> y;
416 if (data.fail()) {
417 std::cerr << m_className << "::LoadData:\n"
418 << " Error reading line " << nLines << ".\n"
419 << " Cannot retrieve element coordinates.\n";
420 bad = true;
421 break;
422 }
423 x *= scaleX;
424 y *= scaleX;
425 const double z = 0.5 * (m_zMin + m_zMax);
426 bool xMirrored, yMirrored, zMirrored;
427 if (!GetElement(x, y, z, i, j, k, xMirrored, yMirrored, zMirrored)) {
428 std::cerr << m_className << "::LoadData:\n"
429 << " Error reading line " << nLines << ".\n"
430 << " Point is outside mesh.\n";
431 bad = true;
432 break;
433 }
434 } else if (fmt == 2) {
435 // "XYZ"
436 double x, y, z;
437 data >> x >> y >> z;
438 if (data.fail()) {
439 std::cerr << m_className << "::LoadData:\n"
440 << " Error reading line " << nLines << ".\n"
441 << " Cannot retrieve element coordinates.\n";
442 bad = true;
443 break;
444 }
445 x *= scaleX;
446 y *= scaleX;
447 z *= scaleX;
448 bool xMirrored, yMirrored, zMirrored;
449 if (!GetElement(x, y, z, i, j, k, xMirrored, yMirrored, zMirrored)) {
450 std::cerr << m_className << "::LoadData:\n"
451 << " Error reading line " << nLines << ".\n"
452 << " Point is outside mesh.\n";
453 bad = true;
454 break;
455 }
456 } else if (fmt == 3) {
457 // "IJ"
458 k = 0;
459 data >> i >> j;
460 if (data.fail()) {
461 std::cerr << m_className << "::LoadData:\n"
462 << " Error reading line " << nLines << ".\n"
463 << " Cannot retrieve element index.\n";
464 bad = true;
465 break;
466 }
467 } else if (fmt == 4) {
468 // "IJK"
469 data >> i >> j >> k;
470 if (data.fail()) {
471 std::cerr << m_className << "::LoadData:\n"
472 << " Error reading line " << nLines << ".\n"
473 << " Cannot retrieve element index.\n";
474 bad = true;
475 break;
476 }
477 } else if (fmt == 5) {
478 // "YXZ"
479 double x, y, z, temp;
480 data >> y >> x >> temp;
481 z = temp;
482 if (data.fail()) {
483 std::cerr << m_className << "::LoadData:\n"
484 << " Error reading line " << nLines << ".\n"
485 << " Cannot retrieve element coordinates.\n";
486 bad = true;
487 break;
488 }
489 x *= scaleX;
490 y *= scaleX;
491 z *= scaleX;
492 bool xMirrored, yMirrored, zMirrored;
493 if (!GetElement(x, y, z, i, j, k, xMirrored, yMirrored, zMirrored)) {
494 std::cerr << m_className << "::LoadData:\n"
495 << " Error reading line " << nLines << ".\n"
496 << " Point is outside mesh.\n";
497 bad = true;
498 break;
499 }
500 }
501 // Check the indices.
502 if (i >= m_nX || j >= m_nY || k >= m_nZ) {
503 std::cerr << m_className << "::LoadData:\n"
504 << " Error reading line " << nLines << ".\n"
505 << " Index (" << i << ", " << j << ", " << k
506 << ") out of range.\n";
507 continue;
508 }
509 if (isSet[i][j][k]) {
510 std::cerr << m_className << "::LoadData:\n"
511 << " Error reading line " << nLines << ".\n"
512 << " Mesh element (" << i << ", " << j << ", " << k
513 << ") has already been set.\n";
514 continue;
515 }
516 // Get the field values.
517 if (fmt == 1 || fmt == 3) {
518 // Two-dimensional field-map
519 fz = 0.;
520 data >> fx >> fy;
521 } else if (fmt == 5) {
522 double temp;
523 data >> fy >> fx >> temp;
524 fz = temp;
525 } else {
526 data >> fx >> fy >> fz;
527 }
528 if (data.fail()) {
529 std::cerr << m_className << "::LoadData:\n"
530 << " Error reading line " << nLines << ".\n"
531 << " Cannot read field values.\n";
532 bad = true;
533 break;
534 }
535 fx *= scaleF;
536 fy *= scaleF;
537 fz *= scaleF;
538 if (withPotential) {
539 data >> v;
540 if (data.fail()) {
541 std::cerr << m_className << "::LoadData:\n"
542 << " Error reading line " << nLines << ".\n"
543 << " Cannot read potential.\n";
544 bad = true;
545 break;
546 }
547 v *= scaleP;
548 if (m_pMin > m_pMax) {
549 // First value.
550 m_pMin = v;
551 m_pMax = v;
552 } else {
553 if (v < m_pMin) m_pMin = v;
554 if (v > m_pMax) m_pMax = v;
555 }
556 }
557 if (withRegion) {
558 data >> region;
559 if (data.fail()) {
560 std::cerr << m_className << "::LoadData:\n"
561 << " Error reading line " << nLines << ".\n"
562 << " Cannot read region.\n";
563 bad = true;
564 break;
565 }
566 }
567 if (fmt == 1 || fmt == 3) {
568 // Two-dimensional field-map
569 for (unsigned int kk = 0; kk < m_nZ; ++kk) {
570 fields[i][j][kk].fx = fx;
571 fields[i][j][kk].fy = fy;
572 fields[i][j][kk].fz = fz;
573 fields[i][j][kk].v = v;
574 if (withRegion) m_regions[i][j][kk] = region;
575 isSet[i][j][kk] = true;
576 }
577 } else {
578 fields[i][j][k].fx = fx;
579 fields[i][j][k].fy = fy;
580 fields[i][j][k].fz = fz;
581 fields[i][j][k].v = v;
582 if (withRegion) m_regions[i][j][k] = region;
583 isSet[i][j][k] = true;
584 }
585 ++nValues;
586 }
587 if (bad) return false;
588 std::cout << m_className << "::LoadData:\n"
589 << " Read " << nValues << " values from " << filename << ".\n";
590 unsigned int nExpected = m_nX * m_nY;
591 if (fmt == 2 || fmt == 4 || fmt == 5) nExpected *= m_nZ;
592 if (nExpected != nValues) {
593 std::cerr << m_className << "::LoadData:\n"
594 << " Expected " << nExpected << " values.\n";
595 }
596 return true;
597}
598
599bool ComponentVoxel::GetBoundingBox(double& xmin, double& ymin, double& zmin,
600 double& xmax, double& ymax, double& zmax) {
601 if (!m_ready) return false;
602 if (m_periodic[0] || m_mirrorPeriodic[0]) {
603 xmin = -INFINITY;
604 xmax = +INFINITY;
605 } else {
606 xmin = m_xMin;
607 xmax = m_xMax;
608 }
609
610 if (m_periodic[1] || m_mirrorPeriodic[1]) {
611 ymin = -INFINITY;
612 ymax = +INFINITY;
613 } else {
614 ymin = m_yMin;
615 ymax = m_yMax;
616 }
617
618 if (m_periodic[2] || m_mirrorPeriodic[2]) {
619 zmin = -INFINITY;
620 zmax = +INFINITY;
621 } else {
622 zmin = m_zMin;
623 zmax = m_zMax;
624 }
625 return true;
626}
627
628bool ComponentVoxel::GetElementaryCell(
629 double& xmin, double& ymin, double& zmin,
630 double& xmax, double& ymax, double& zmax) {
631 if (!m_ready) return false;
632 xmin = m_xMin;
633 xmax = m_xMax;
634 ymin = m_yMin;
635 ymax = m_yMax;
636 zmin = m_zMin;
637 zmax = m_zMax;
638 return true;
639}
640
641bool ComponentVoxel::GetVoltageRange(double& vmin, double& vmax) {
642 if (!m_ready) return false;
643 vmin = m_pMin;
644 vmax = m_pMax;
645 return true;
646}
647
648bool ComponentVoxel::GetElectricFieldRange(double& exmin, double& exmax,
649 double& eymin, double& eymax,
650 double& ezmin, double& ezmax) {
651 if (!m_ready) {
652 std::cerr << m_className << "::GetElectricFieldRange:\n"
653 << " Field map is not ready for interpolation.\n";
654 return false;
655 }
656
657 exmin = exmax = m_efields[0][0][0].fx;
658 eymin = eymax = m_efields[0][0][0].fy;
659 ezmin = ezmax = m_efields[0][0][0].fz;
660 for (unsigned int i = 0; i < m_nX; ++i) {
661 for (unsigned int j = 0; j < m_nY; ++j) {
662 for (unsigned int k = 0; k < m_nZ; ++k) {
663 const Element& element = m_efields[i][j][k];
664 if (element.fx < exmin) exmin = element.fx;
665 if (element.fx > exmax) exmax = element.fx;
666 if (element.fy < eymin) eymin = element.fy;
667 if (element.fy > eymax) eymax = element.fy;
668 if (element.fz < ezmin) ezmin = element.fz;
669 if (element.fz > ezmax) ezmax = element.fz;
670 }
671 }
672 }
673 return true;
674}
675
676void ComponentVoxel::PrintRegions() const {
677 // Do not proceed if not properly initialised.
678 if (!m_ready) {
679 std::cerr << m_className << "::PrintRegions:\n"
680 << " Field map not yet initialised.\n";
681 return;
682 }
683
684 if (m_media.empty()) {
685 std::cerr << m_className << "::PrintRegions: No regions defined.\n";
686 return;
687 }
688
689 std::cout << m_className << "::PrintRegions:\n";
690 std::cout << " Index Medium\n";
691 const unsigned int nMedia = m_media.size();
692 for (unsigned int i = 0; i < nMedia; ++i) {
693 const std::string name = m_media[i] ? m_media[i]->GetName() : "none";
694 std::cout << " " << i << " " << name << "\n";
695 }
696}
697
698void ComponentVoxel::SetMedium(const unsigned int i, Medium* m) {
699 if (!m) {
700 std::cerr << m_className << "::SetMedium: Null pointer.\n";
701 if (m_media.empty()) return;
702 }
703 if (i >= m_media.size()) m_media.resize(i + 1, nullptr);
704 m_media[i] = m;
705}
706
707Medium* ComponentVoxel::GetMedium(const unsigned int i) const {
708 if (i >= m_media.size()) {
709 std::cerr << m_className << "::GetMedium: Index out of range.\n";
710 return nullptr;
711 }
712 return m_media[i];
713}
714
715bool ComponentVoxel::GetField(
716 const double xi, const double yi, const double zi,
717 const std::vector<std::vector<std::vector<Element> > >& field, double& fx,
718 double& fy, double& fz, double& p, int& region) {
719 if (!m_hasMesh) {
720 std::cerr << m_className << "::GetField: Mesh is not set.\n";
721 return false;
722 }
723
724 // Reduce the point to the basic cell (in case of periodicity) and
725 // check if it is inside the mesh.
726 bool xMirrored = false;
727 const double x =
728 Reduce(xi, m_xMin, m_xMax, m_periodic[0], m_mirrorPeriodic[0], xMirrored);
729 if (x < m_xMin || x > m_xMax) return false;
730 bool yMirrored = false;
731 const double y =
732 Reduce(yi, m_yMin, m_yMax, m_periodic[1], m_mirrorPeriodic[1], yMirrored);
733 if (y < m_yMin || y > m_yMax) return false;
734 bool zMirrored = false;
735 const double z =
736 Reduce(zi, m_zMin, m_zMax, m_periodic[2], m_mirrorPeriodic[2], zMirrored);
737 if (z < m_zMin || z > m_zMax) return false;
738
739 // Get the indices.
740 const double sx = (x - m_xMin) / m_dx;
741 const double sy = (y - m_yMin) / m_dy;
742 const double sz = (z - m_zMin) / m_dz;
743 unsigned int i = static_cast<unsigned int>(sx);
744 unsigned int j = static_cast<unsigned int>(sy);
745 unsigned int k = static_cast<unsigned int>(sz);
746 if (i >= m_nX) i = m_nX - 1;
747 if (j >= m_nY) j = m_nY - 1;
748 if (k >= m_nZ) k = m_nZ - 1;
749 region = m_regions[i][j][k];
750
751 // Get the field and potential.
752 if (m_interpolate) {
753 // Get the "nodes" (voxel centres) surrounding the point.
754 const double tx = sx - 0.5;
755 const double ty = sy - 0.5;
756 const double tz = sz - 0.5;
757 int i0 = static_cast<int>(std::floor(tx));
758 int j0 = static_cast<int>(std::floor(ty));
759 int k0 = static_cast<int>(std::floor(tz));
760 double vx = tx - i0;
761 double vy = ty - j0;
762 double vz = tz - k0;
763 unsigned int i1 = i0 + 1;
764 unsigned int j1 = j0 + 1;
765 unsigned int k1 = k0 + 1;
766 const bool perx = m_periodic[0] || m_mirrorPeriodic[0];
767 const bool pery = m_periodic[1] || m_mirrorPeriodic[1];
768 const bool perz = m_periodic[2] || m_mirrorPeriodic[2];
769 if (i0 < 0) {
770 if (perx) {
771 i0 = m_nX - 1;
772 } else {
773 i0 = 0;
774 vx = 0.;
775 }
776 }
777 if (j0 < 0) {
778 if (pery) {
779 j0 = m_nY - 1;
780 } else {
781 j0 = 0;
782 vy = 0.;
783 }
784 }
785 if (k0 < 0) {
786 if (perz) {
787 k0 = m_nZ - 1;
788 } else {
789 k0 = 0;
790 vz = 0.;
791 }
792 }
793 if (i1 >= m_nX) i1 = perx ? 0 : m_nX - 1;
794 if (j1 >= m_nY) j1 = pery ? 0 : m_nY - 1;
795 if (k1 >= m_nZ) k1 = perz ? 0 : m_nZ - 1;
796 const Element& n000 = field[i0][j0][k0];
797 const Element& n100 = field[i1][j0][k0];
798 const Element& n010 = field[i0][j1][k0];
799 const Element& n110 = field[i1][j1][k0];
800 const Element& n001 = field[i0][j0][k1];
801 const Element& n101 = field[i1][j0][k1];
802 const Element& n011 = field[i0][j1][k1];
803 const Element& n111 = field[i1][j1][k1];
804
805 const double ux = 1. - vx;
806 const double uy = 1. - vy;
807 const double uz = 1. - vz;
808 if (m_debug) {
809 std::cout << m_className << "::GetField:\n Determining field at ("
810 << xi << ", " << yi << ", " << zi << ").\n"
811 << " X: " << i0 << " (" << ux << ") - "
812 << i1 << " (" << vx << ").\n"
813 << " Y: " << j0 << " (" << uy << ") - "
814 << j1 << " (" << vy << ").\n"
815 << " Z: " << k0 << " (" << uz << ") - "
816 << k1 << " (" << vz << ").\n";
817 }
818 fx = ((n000.fx * ux + n100.fx * vx) * uy +
819 (n010.fx * ux + n110.fx * vx) * vy) *
820 uz +
821 ((n001.fx * ux + n101.fx * vx) * uy +
822 (n011.fx * ux + n111.fx * vx) * vy) *
823 vz;
824 fy = ((n000.fy * ux + n100.fy * vx) * uy +
825 (n010.fy * ux + n110.fy * vx) * vy) *
826 uz +
827 ((n001.fy * ux + n101.fy * vx) * uy +
828 (n011.fy * ux + n111.fy * vx) * vy) *
829 vz;
830 fz = ((n000.fz * ux + n100.fz * vx) * uy +
831 (n010.fz * ux + n110.fz * vx) * vy) *
832 uz +
833 ((n001.fz * ux + n101.fz * vx) * uy +
834 (n011.fz * ux + n111.fz * vx) * vy) *
835 vz;
836 p = ((n000.v * ux + n100.v * vx) * uy + (n010.v * ux + n110.v * vx) * vy) *
837 uz +
838 ((n001.v * ux + n101.v * vx) * uy + (n011.v * ux + n111.v * vx) * vy) *
839 vz;
840 } else {
841 const Element& element = field[i][j][k];
842 fx = element.fx;
843 fy = element.fy;
844 fz = element.fz;
845 p = element.v;
846 }
847 if (xMirrored) fx = -fx;
848 if (yMirrored) fy = -fy;
849 if (zMirrored) fz = -fz;
850 return true;
851}
852
853bool ComponentVoxel::GetElement(const double xi, const double yi,
854 const double zi, unsigned int& i,
855 unsigned int& j, unsigned int& k,
856 bool& xMirrored, bool& yMirrored,
857 bool& zMirrored) const {
858 if (!m_hasMesh) {
859 std::cerr << m_className << "::GetElement: Mesh is not set.\n";
860 return false;
861 }
862
863 // Reduce the point to the basic cell (in case of periodicity) and
864 // check if it is inside the mesh.
865 const double x =
866 Reduce(xi, m_xMin, m_xMax, m_periodic[0], m_mirrorPeriodic[0], xMirrored);
867 if (x < m_xMin || x > m_xMax) return false;
868 const double y =
869 Reduce(yi, m_yMin, m_yMax, m_periodic[1], m_mirrorPeriodic[1], yMirrored);
870 if (y < m_yMin || y > m_yMax) return false;
871 const double z =
872 Reduce(zi, m_zMin, m_zMax, m_periodic[2], m_mirrorPeriodic[2], zMirrored);
873 if (z < m_zMin || z > m_zMax) return false;
874
875 // Get the indices.
876 i = (unsigned int)((x - m_xMin) / m_dx);
877 j = (unsigned int)((y - m_yMin) / m_dy);
878 k = (unsigned int)((z - m_zMin) / m_dz);
879 if (i >= m_nX) i = m_nX - 1;
880 if (j >= m_nY) j = m_nY - 1;
881 if (k >= m_nZ) k = m_nZ - 1;
882 return true;
883}
884
885bool ComponentVoxel::GetElement(const unsigned int i, const unsigned int j,
886 const unsigned int k, double& v, double& ex,
887 double& ey, double& ez) const {
888 v = ex = ey = ez = 0.;
889 if (!m_ready) {
890 if (!m_hasMesh) {
891 std::cerr << m_className << "::GetElement: Mesh not set.\n";
892 return false;
893 }
894 std::cerr << m_className << "::GetElement: Field map not set.\n";
895 return false;
896 }
897 if (i >= m_nX || j >= m_nY || k >= m_nZ) {
898 std::cerr << m_className << "::GetElement: Index out of range.\n";
899 return false;
900 }
901 const Element& element = m_efields[i][j][k];
902 v = element.v;
903 ex = element.fx;
904 ey = element.fy;
905 ez = element.fz;
906 return true;
907}
908
909void ComponentVoxel::Reset() {
910 m_regions.clear();
911 m_efields.clear();
912 m_bfields.clear();
913 m_wfields.clear();
914
915 m_wdfields.clear();
916 m_wdtimes.clear();
917
918 m_nX = m_nY = m_nZ = 0;
919 m_xMin = m_yMin = m_zMin = 0.;
920 m_xMax = m_yMax = m_zMax = 0.;
921 m_pMin = m_pMax = 0.;
922 m_media.clear();
923
924 m_hasMesh = false;
925 m_hasPotential = false;
926 m_hasEfield = false;
927 m_hasBfield = false;
928 m_hasWfield = false;
929 m_ready = false;
930
931 m_wField_xOffset = 0.;
932 m_wField_yOffset = 0.;
933 m_wField_zOffset = 0.;
934}
935
936void ComponentVoxel::UpdatePeriodicity() {
937 if (!m_ready) {
938 std::cerr << m_className << "::UpdatePeriodicity:\n"
939 << " Field map not available.\n";
940 return;
941 }
942
943 // Check for conflicts.
944 for (unsigned int i = 0; i < 3; ++i) {
945 if (m_periodic[i] && m_mirrorPeriodic[i]) {
946 std::cerr << m_className << "::UpdatePeriodicity:\n"
947 << " Both simple and mirror periodicity requested. Reset.\n";
948 m_periodic[i] = m_mirrorPeriodic[i] = false;
949 }
950 }
951
952 if (m_axiallyPeriodic[0] || m_axiallyPeriodic[1] || m_axiallyPeriodic[2]) {
953 std::cerr << m_className << "::UpdatePeriodicity:\n"
954 << " Axial symmetry is not supported. Reset.\n";
955 m_axiallyPeriodic.fill(false);
956 }
957
958 if (m_rotationSymmetric[0] || m_rotationSymmetric[1] ||
959 m_rotationSymmetric[2]) {
960 std::cerr << m_className << "::UpdatePeriodicity:\n"
961 << " Rotation symmetry is not supported. Reset.\n";
962 m_rotationSymmetric.fill(false);
963 }
964}
965
966double ComponentVoxel::Reduce(const double xin, const double xmin,
967 const double xmax, const bool simplePeriodic,
968 const bool mirrorPeriodic, bool& mirrored) const {
969 // In case of periodicity, reduce the coordinate to the basic cell.
970 double x = xin;
971 const double lx = xmax - xmin;
972 if (simplePeriodic) {
973 x = xmin + fmod(x - xmin, lx);
974 if (x < xmin) x += lx;
975 } else if (mirrorPeriodic) {
976 double xNew = xmin + fmod(x - xmin, lx);
977 if (xNew < xmin) xNew += lx;
978 const int nx = int(floor(0.5 + (xNew - x) / lx));
979 if (nx != 2 * (nx / 2)) {
980 xNew = xmin + xmax - xNew;
981 mirrored = true;
982 }
983 x = xNew;
984 }
985 return x;
986}
987
988void ComponentVoxel::Initialise(
989 std::vector<std::vector<std::vector<Element> > >& fields) {
990
991 fields.resize(m_nX);
992 for (unsigned int i = 0; i < m_nX; ++i) {
993 fields[i].resize(m_nY);
994 for (unsigned int j = 0; j < m_nY; ++j) {
995 fields[i][j].resize(m_nZ);
996 for (unsigned int k = 0; k < m_nZ; ++k) {
997 fields[i][j][k].fx = 0.;
998 fields[i][j][k].fy = 0.;
999 fields[i][j][k].fz = 0.;
1000 fields[i][j][k].v = 0.;
1001 }
1002 }
1003 }
1004}
1005
1006void ComponentVoxel::InitialiseRegions() {
1007 if (!m_hasMesh) return;
1008 m_regions.resize(m_nX);
1009 for (unsigned int i = 0; i < m_nX; ++i) {
1010 m_regions[i].resize(m_nY);
1011 for (unsigned int j = 0; j < m_nY; ++j) {
1012 m_regions[i][j].assign(m_nZ, 0);
1013 }
1014 }
1015}
1016}
Garfield::ComponentVoxel::HasMagneticField
bool HasMagneticField() const override
Does the component have a non-zero magnetic field?
Definition ComponentVoxel.cc:180
Garfield::ComponentVoxel::MagneticField
void MagneticField(const double x, const double y, const double z, double &bx, double &by, double &bz, int &status) override
Definition ComponentVoxel.cc:165
Garfield::ComponentVoxel::SetWeightingFieldOffset
void SetWeightingFieldOffset(const double x, const double y, const double z)
Definition ComponentVoxel.cc:158
Garfield::ComponentVoxel::PrintRegions
void PrintRegions() const
Print all regions.
Definition ComponentVoxel.cc:676
Garfield::ComponentVoxel::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 ComponentVoxel.cc:271
Garfield::ComponentVoxel::GetBoundingBox
bool GetBoundingBox(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax) override
Get the bounding box coordinates.
Definition ComponentVoxel.cc:599
Garfield::ComponentVoxel::GetElectricFieldRange
bool GetElectricFieldRange(double &exmin, double &exmax, double &eymin, double &eymax, double &ezmin, double &ezmax)
Definition ComponentVoxel.cc:648
Garfield::ComponentVoxel::GetMedium
Medium * GetMedium(const unsigned int i) const
Get the medium in region i.
Definition ComponentVoxel.cc:707
Garfield::ComponentVoxel::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 ComponentVoxel.cc:82
Garfield::ComponentVoxel::GetElement
bool GetElement(const double xi, const double yi, const double zi, unsigned int &i, unsigned int &j, unsigned int &k, bool &xMirrored, bool &yMirrored, bool &zMirrored) const
Return the indices of the element at a given point.
Definition ComponentVoxel.cc:853
Garfield::ComponentVoxel::GetVoltageRange
bool GetVoltageRange(double &vmin, double &vmax) override
Calculate the voltage range [V].
Definition ComponentVoxel.cc:641
Garfield::ComponentVoxel::GetElementaryCell
bool GetElementaryCell(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax) override
Get the coordinates of the elementary cell.
Definition ComponentVoxel.cc:628
Garfield::ComponentVoxel::SetMedium
void SetMedium(const unsigned int i, Medium *m)
Set the medium in region i.
Definition ComponentVoxel.cc:698
Garfield::ComponentVoxel::WeightingField
void WeightingField(const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label) override
Definition ComponentVoxel.cc:53
Garfield::ComponentVoxel::DelayedWeightingPotential
double DelayedWeightingPotential(const double x, const double y, const double z, const double t, const std::string &label) override
Definition ComponentVoxel.cc:124
Garfield::ComponentVoxel::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 ComponentVoxel.cc:16
Garfield::ComponentVoxel::SetMesh
void 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 ComponentVoxel.cc:203
Garfield::ComponentVoxel::WeightingPotential
double WeightingPotential(const double x, const double y, const double z, const std::string &label) override
Definition ComponentVoxel.cc:67
Garfield::ComponentVoxel::LoadElectricField
bool LoadElectricField(const std::string &filename, const std::string &format, const bool withPotential, const bool withRegion, const double scaleX=1., const double scaleE=1., const double scaleP=1.)
Definition ComponentVoxel.cc:239
Garfield::ComponentVoxel::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 ComponentVoxel.cc:327
Garfield::Component::m_rotationSymmetric
std::array< bool, 3 > m_rotationSymmetric
Rotation symmetry around x-axis, y-axis, z-axis.
Definition Component.hh:380
Garfield::Component::m_mirrorPeriodic
std::array< bool, 3 > m_mirrorPeriodic
Mirror periodicity in x, y, z.
Definition Component.hh:376
Garfield::Component::m_debug
bool m_debug
Switch on/off debugging messages.
Definition Component.hh:371
Garfield::Component::Component
Component()=delete
Default constructor.
Garfield::Component::m_periodic
std::array< bool, 3 > m_periodic
Simple periodicity in x, y, z.
Definition Component.hh:374
Garfield::Component::HasMagneticField
virtual bool HasMagneticField() const
Does the component have a non-zero magnetic field?
Definition Component.cc:155
Garfield::Component::m_className
std::string m_className
Class name.
Definition Component.hh:359
Garfield::Component::MagneticField
virtual void MagneticField(const double x, const double y, const double z, double &bx, double &by, double &bz, int &status)
Definition Component.cc:102
Garfield::Component::m_ready
bool m_ready
Ready for use?
Definition Component.hh:368
Garfield::Component::m_axiallyPeriodic
std::array< bool, 3 > m_axiallyPeriodic
Axial periodicity in x, y, z.
Definition Component.hh:378
Garfield::Medium
Abstract base class for media.
Definition Medium.hh:16
Garfield::ltrim
void ltrim(std::string &line)
Definition Utilities.hh:12