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