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