Garfield++ v2r0
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 <iostream>
2#include <fstream>
3#include <sstream>
4#include <string>
5#include <algorithm>
6#include <cmath>
7
8#include "ComponentVoxel.hh"
9
10namespace Garfield {
11
12ComponentVoxel::ComponentVoxel()
13 : ComponentBase(),
14 m_nX(0),
15 m_nY(0),
16 m_nZ(0),
17 m_xMin(0.),
18 m_yMin(0.),
19 m_zMin(0.),
20 m_xMax(0.),
21 m_yMax(0.),
22 m_zMax(0.),
23 m_hasMesh(false),
24 m_hasPotential(false),
25 m_hasEfield(false),
26 m_hasBfield(false),
27 m_pMin(0.),
28 m_pMax(0.) {
29
30 m_className = "ComponentVoxel";
31}
32
33void ComponentVoxel::ElectricField(const double x, const double y,
34 const double z, double& ex, double& ey,
35 double& ez, double& p, Medium*& m,
36 int& status) {
37
38 m = NULL;
39 // Make sure the field map has been loaded.
40 if (!m_ready) {
41 std::cerr << m_className << "::ElectricField:\n"
42 << " Field map is not available for interpolation.\n";
43 status = -10;
44 return;
45 }
46
47 // Get the mesh element.
48 unsigned int i = 0, j = 0, k = 0;
49 bool xMirrored = false, yMirrored = false, zMirrored = false;
50 if (!GetElement(x, y, z, i, j, k, xMirrored, yMirrored, zMirrored)) {
51 status = -11;
52 return;
53 }
54 status = 0;
55 // Get the electric field and potential.
56 const Element& element = m_efields[i][j][k];
57 ex = element.fx;
58 ey = element.fy;
59 ez = element.fz;
60 if (xMirrored) ex = -ex;
61 if (yMirrored) ey = -ey;
62 if (zMirrored) ez = -ez;
63 p = element.v;
64 // Get the medium.
65 const int region = m_regions[i][j][k];
66 if (region < 0 || region > (int)m_media.size()) {
67 m = NULL;
68 status = -5;
69 return;
70 }
71 m = m_media[region];
72 if (!m) status = -5;
73}
74
75void ComponentVoxel::ElectricField(const double x, const double y,
76 const double z, double& ex, double& ey,
77 double& ez, Medium*& m, int& status) {
78
79 double v = 0.;
80 ElectricField(x, y, z, ex, ey, ez, v, m, status);
81}
82
83void ComponentVoxel::WeightingField(const double x, const double y,
84 const double z, double& wx, double& wy,
85 double& wz, const std::string& /*label*/) {
86 int status = 0;
87 Medium* med = NULL;
88 double v = 0.;
89 double x1 = x - m_wField_xOffset;
90 double y1 = y - m_wField_yOffset;
91 double z1 = z - m_wField_zOffset;
92 ElectricField(x1, y1, z1, wx, wy, wz, v, med, status);
93}
94
95void ComponentVoxel::SetWeightingFieldOffset(const double x, const double y,
96 const double z) {
97 m_wField_xOffset = x;
98 m_wField_yOffset = y;
99 m_wField_zOffset = z;
100}
101
102void ComponentVoxel::MagneticField(const double x, const double y,
103 const double z,
104 double& bx, double& by, double& bz,
105 int& status) {
106
107 if (!m_hasBfield) {
108 return ComponentBase::MagneticField(x, y, z, bx, by, bz, status);
109 }
110
111 // Get the mesh element.
112 unsigned int i = 0, j = 0, k = 0;
113 bool xMirrored = false, yMirrored = false, zMirrored = false;
114 if (!GetElement(x, y, z, i, j, k, xMirrored, yMirrored, zMirrored)) {
115 status = -11;
116 return;
117 }
118 status = 0;
119 // Get the field.
120 const Element& element = m_bfields[i][j][k];
121 bx = element.fx;
122 by = element.fy;
123 bz = element.fz;
124 if (xMirrored) bx = -bx;
125 if (yMirrored) by = -by;
126 if (zMirrored) bz = -bz;
127}
128
129Medium* ComponentVoxel::GetMedium(const double x, const double y,
130 const double z) {
131
132 // Make sure the field map has been loaded.
133 if (!m_ready) {
134 std::cerr << m_className << "::GetMedium:\n"
135 << " Field map is not available for interpolation.\n";
136 return NULL;
137 }
138
139 unsigned int i, j, k;
140 bool xMirrored, yMirrored, zMirrored;
141 if (!GetElement(x, y, z, i, j, k, xMirrored, yMirrored, zMirrored)) {
142 return NULL;
143 }
144 const int region = m_regions[i][j][k];
145 if (region < 0 || region > (int)m_media.size()) return NULL;
146 return m_media[region];
147}
148
149void ComponentVoxel::SetMesh(const unsigned int nx, const unsigned int ny,
150 const unsigned int nz,
151 const double xmin, const double xmax,
152 const double ymin, const double ymax,
153 const double zmin, const double zmax) {
154
155 Reset();
156 if (nx == 0 || ny == 0 || nz == 0) {
157 std::cerr << m_className << "::SetMesh:\n"
158 << " Number of mesh elements must be positive.\n";
159 return;
160 }
161 if (xmin >= xmax) {
162 std::cerr << m_className << "::SetMesh:\n Invalid x range.\n";
163 return;
164 } else if (ymin >= ymax) {
165 std::cerr << m_className << "::SetMesh:\n Invalid y range.\n";
166 return;
167 } else if (zmin >= zmax) {
168 std::cerr << m_className << "::SetMesh:\n Invalid z range.\n";
169 return;
170 }
171 m_nX = nx;
172 m_nY = ny;
173 m_nZ = nz;
174 m_xMin = xmin;
175 m_yMin = ymin;
176 m_zMin = zmin;
177 m_xMax = xmax;
178 m_yMax = ymax;
179 m_zMax = zmax;
180 m_hasMesh = true;
181}
182
183bool ComponentVoxel::LoadElectricField(const std::string& filename,
184 const std::string& format,
185 const bool withPotential,
186 const bool withRegion,
187 const double scaleX,
188 const double scaleE,
189 const double scaleP) {
190
191 m_ready = false;
192 m_hasPotential = m_hasEfield = false;
193 if (!m_hasMesh) {
194 std::cerr << m_className << "::LoadElectricField:\n"
195 << " Mesh is not set. Call SetMesh first.\n";
196 return false;
197 }
198
199 // Set up the grid.
200 m_efields.resize(m_nX);
201 m_regions.resize(m_nX);
202 for (unsigned int i = 0; i < m_nX; ++i) {
203 m_efields[i].resize(m_nY);
204 m_regions[i].resize(m_nY);
205 for (unsigned int j = 0; j < m_nY; ++j) {
206 m_efields[i][j].resize(m_nZ);
207 m_regions[i][j].resize(m_nZ);
208 for (unsigned int k = 0; k < m_nZ; ++k) {
209 m_efields[i][j][k].fx = 0.;
210 m_efields[i][j][k].fy = 0.;
211 m_efields[i][j][k].fz = 0.;
212 m_efields[i][j][k].v = 0.;
213 m_regions[i][j][k] = 0;
214 }
215 }
216 }
217
218 m_pMin = m_pMax = 0.;
219 if (withPotential) {
220 m_pMin = 1.;
221 m_pMax = -1.;
222 }
223 return LoadData(filename, format, withPotential, withRegion,
224 scaleX, scaleE, scaleP, 'e');
225}
226
227bool ComponentVoxel::LoadMagneticField(const std::string& filename,
228 const std::string& format,
229 const double scaleX,
230 const double scaleB) {
231
232 m_hasBfield = false;
233 if (!m_hasMesh) {
234 std::cerr << m_className << "::LoadMagneticField:\n"
235 << " Mesh is not set. Call SetMesh first.\n";
236 return false;
237 }
238
239 // Set up the grid.
240 m_bfields.resize(m_nX);
241 for (unsigned int i = 0; i < m_nX; ++i) {
242 m_bfields[i].resize(m_nY);
243 for (unsigned int j = 0; j < m_nY; ++j) {
244 m_bfields[i][j].resize(m_nZ);
245 for (unsigned int k = 0; k < m_nZ; ++k) {
246 m_bfields[i][j][k].fx = 0.;
247 m_bfields[i][j][k].fy = 0.;
248 m_bfields[i][j][k].fz = 0.;
249 m_bfields[i][j][k].v = 0.;
250 }
251 }
252 }
253
254 return LoadData(filename, format, false, false, scaleX, scaleB, 1., 'b');
255}
256
257bool ComponentVoxel::LoadData(const std::string& filename, std::string format,
258 const bool withPotential, const bool withRegion,
259 const double scaleX, const double scaleF,
260 const double scaleP, const char field) {
261
262 if (!m_hasMesh) {
263 std::cerr << m_className << "::LoadData:\n Mesh has not been set.\n";
264 return false;
265 }
266
267 unsigned int nValues = 0;
268 // Keep track of which elements have been read.
269 std::vector<std::vector<std::vector<bool> > > isSet(m_nX,
270 std::vector<std::vector<bool> >(m_nY, std::vector<bool>(m_nZ, false)));
271
272 std::ifstream infile;
273 infile.open(filename.c_str(), std::ios::in);
274 if (!infile) {
275 std::cerr << m_className << "::LoadData:\n"
276 << " Could not open file " << filename << ".\n";
277 return false;
278 }
279
280 std::transform(format.begin(), format.end(), format.begin(), toupper);
281 unsigned int fmt = 0;
282 if (format == "XY") {
283 fmt = 1;
284 } else if (format == "XYZ") {
285 fmt = 2;
286 } else if (format == "IJ") {
287 fmt = 3;
288 } else if (format == "IJK") {
289 fmt = 4;
290 } else if (format == "YXZ") {
291 fmt = 5;
292 } else {
293 std::cerr << m_className << "::LoadData:\n"
294 << " Unkown format (" << format << ").\n";
295 return false;
296 }
297 std::string line;
298 unsigned int nLines = 0;
299 bool bad = false;
300 while (!infile.fail()) {
301 // Read one line.
302 std::getline(infile, line);
303 ++nLines;
304 // Strip white space from beginning of line.
305 line.erase(line.begin(),
306 std::find_if(line.begin(), line.end(),
307 not1(std::ptr_fun<int, int>(isspace))));
308 // Skip empty lines.
309 if (line.empty()) continue;
310 // Skip comments.
311 if (line[0] == '#') continue;
312 if (line[0] == '/' && line[1] == '/') continue;
313 unsigned int i = 0;
314 unsigned int j = 0;
315 unsigned int k = 0;
316 double fx = 0.;
317 double fy = 0.;
318 double fz = 0.;
319 double v = 0.;
320 int region = 0;
321 std::istringstream data;
322 data.str(line);
323 if (fmt == 1) {
324 // "XY"
325 double x, y;
326 data >> x >> y;
327 if (data.fail()) {
328 std::cerr << m_className << "::LoadData:\n"
329 << " Error reading line " << nLines << ".\n"
330 << " Cannot retrieve element coordinates.\n";
331 bad = true;
332 break;
333 }
334 x *= scaleX;
335 y *= scaleX;
336 const double z = 0.5 * (m_zMin + m_zMax);
337 bool xMirrored, yMirrored, zMirrored;
338 if (!GetElement(x, y, z, i, j, k, xMirrored, yMirrored, zMirrored)) {
339 std::cerr << m_className << "::LoadData:\n"
340 << " Error reading line " << nLines << ".\n"
341 << " Point is outside mesh.\n";
342 bad = true;
343 break;
344 }
345 } else if (fmt == 2) {
346 // "XYZ"
347 double x, y, z;
348 data >> x >> y >> z;
349 if (data.fail()) {
350 std::cerr << m_className << "::LoadData:\n"
351 << " Error reading line " << nLines << ".\n"
352 << " Cannot retrieve element coordinates.\n";
353 bad = true;
354 break;
355 }
356 x *= scaleX;
357 y *= scaleX;
358 z *= scaleX;
359 bool xMirrored, yMirrored, zMirrored;
360 if (!GetElement(x, y, z, i, j, k, xMirrored, yMirrored, zMirrored)) {
361 std::cerr << m_className << "::LoadData:\n"
362 << " Error reading line " << nLines << ".\n"
363 << " Point is outside mesh.\n";
364 bad = true;
365 break;
366 }
367 } else if (fmt == 3) {
368 // "IJ"
369 k = 0;
370 data >> i >> j;
371 if (data.fail()) {
372 std::cerr << m_className << "::LoadData:\n"
373 << " Error reading line " << nLines << ".\n"
374 << " Cannot retrieve element index.\n";
375 bad = true;
376 break;
377 }
378 } else if (fmt == 4) {
379 // "IJK"
380 data >> i >> j >> k;
381 if (data.fail()) {
382 std::cerr << m_className << "::LoadData:\n"
383 << " Error reading line " << nLines << ".\n"
384 << " Cannot retrieve element index.\n";
385 bad = true;
386 break;
387 }
388 } else if (fmt == 5) {
389 // "YXZ"
390 double x, y, z, temp;
391 data >> y >> x >> temp;
392 z = temp;
393 if (data.fail()) {
394 std::cerr << m_className << "::LoadData:\n"
395 << " Error reading line " << nLines << ".\n"
396 << " Cannot retrieve element coordinates.\n";
397 bad = true;
398 break;
399 }
400 x *= scaleX;
401 y *= scaleX;
402 z *= scaleX;
403 bool xMirrored, yMirrored, zMirrored;
404 if (!GetElement(x, y, z, i, j, k, xMirrored, yMirrored, zMirrored)) {
405 std::cerr << m_className << "::LoadData:\n"
406 << " Error reading line " << nLines << ".\n"
407 << " Point is outside mesh.\n";
408 bad = true;
409 break;
410 }
411 }
412 // Check the indices.
413 if (i >= m_nX || j >= m_nY || k >= m_nZ) {
414 std::cerr << m_className << "::LoadData:\n"
415 << " Error reading line " << nLines << ".\n"
416 << " Index (" << i << ", " << j << ", " << k
417 << ") out of range.\n";
418 continue;
419 }
420 if (isSet[i][j][k]) {
421 std::cerr << m_className << "::LoadData:\n"
422 << " Error reading line " << nLines << ".\n"
423 << " Mesh element (" << i << ", " << j << ", " << k
424 << ") has already been set.\n";
425 continue;
426 }
427 // Get the field values.
428 if (fmt == 1 || fmt == 3) {
429 // Two-dimensional field-map
430 fz = 0.;
431 data >> fx >> fy;
432 } else if (fmt == 5) {
433 double temp;
434 data >> fy >> fx >> temp;
435 fz = temp;
436 } else {
437 data >> fx >> fy >> fz;
438 }
439 if (data.fail()) {
440 std::cerr << m_className << "::LoadData:\n"
441 << " Error reading line " << nLines << ".\n"
442 << " Cannot read field values.\n";
443 bad = true;
444 break;
445 }
446 fx *= scaleF;
447 fy *= scaleF;
448 fz *= scaleF;
449 if (withPotential) {
450 data >> v;
451 if (data.fail()) {
452 std::cerr << m_className << "::LoadData:\n"
453 << " Error reading line " << nLines << ".\n"
454 << " Cannot read potential.\n";
455 bad = true;
456 break;
457 }
458 v *= scaleP;
459 if (m_pMin > m_pMax) {
460 // First value.
461 m_pMin = v;
462 m_pMax = v;
463 } else {
464 if (v < m_pMin) m_pMin = v;
465 if (v > m_pMax) m_pMax = v;
466 }
467 }
468 if (withRegion) {
469 data >> region;
470 if (data.fail()) {
471 std::cerr << m_className << "::LoadData:\n"
472 << " Error reading line " << nLines << ".\n"
473 << " Cannot read region.\n";
474 bad = true;
475 break;
476 }
477 }
478 if (fmt == 1 || fmt == 3) {
479 // Two-dimensional field-map
480 for (unsigned int kk = 0; kk < m_nZ; ++kk) {
481 if (field == 'e') {
482 m_efields[i][j][kk].fx = fx;
483 m_efields[i][j][kk].fy = fy;
484 m_efields[i][j][kk].fz = fz;
485 m_efields[i][j][kk].v = v;
486 m_regions[i][j][kk] = region;
487 } else if (field == 'b') {
488 m_bfields[i][j][kk].fx = fx;
489 m_bfields[i][j][kk].fy = fy;
490 m_bfields[i][j][kk].fz = fz;
491 }
492 isSet[i][j][kk] = true;
493 }
494 } else {
495 if (field == 'e') {
496 m_efields[i][j][k].fx = fx;
497 m_efields[i][j][k].fy = fy;
498 m_efields[i][j][k].fz = fz;
499 m_efields[i][j][k].v = v;
500 m_regions[i][j][k] = region;
501 } else if (field == 'b') {
502 m_bfields[i][j][k].fx = fx;
503 m_bfields[i][j][k].fy = fy;
504 m_bfields[i][j][k].fz = fz;
505 }
506 isSet[i][j][k] = true;
507 }
508 ++nValues;
509 }
510 if (bad) return false;
511 std::cout << m_className << "::LoadData:\n"
512 << " Read " << nValues << " values from " << filename << ".\n";
513 unsigned int nExpected = m_nX * m_nY;
514 if (fmt == 2 || fmt == 4 || fmt == 5) nExpected *= m_nZ;
515 if (nExpected != nValues) {
516 std::cerr << m_className << "::LoadData:\n"
517 << " Expected " << nExpected << " values.\n";
518 }
519 if (field == 'e') {
520 m_hasEfield = true;
521 m_ready = true;
522 if (withPotential) m_hasPotential = true;
523 } else if (field == 'b') {
524 m_hasBfield = true;
525 }
526 return true;
527}
528
529bool ComponentVoxel::GetBoundingBox(double& xmin, double& ymin, double& zmin,
530 double& xmax, double& ymax, double& zmax) {
531
532 if (!m_ready) return false;
533 if (m_xPeriodic || m_xMirrorPeriodic) {
534 xmin = -INFINITY;
535 xmax = +INFINITY;
536 } else {
537 xmin = m_xMin;
538 xmax = m_xMax;
539 }
540
541 if (m_yPeriodic || m_yMirrorPeriodic) {
542 ymin = -INFINITY;
543 ymax = +INFINITY;
544 } else {
545 ymin = m_yMin;
546 ymax = m_yMax;
547 }
548
549 if (m_zPeriodic || m_zMirrorPeriodic) {
550 zmin = -INFINITY;
551 zmax = +INFINITY;
552 } else {
553 zmin = m_zMin;
554 zmax = m_zMax;
555 }
556 return true;
557}
558
559bool ComponentVoxel::GetVoltageRange(double& vmin, double& vmax) {
560
561 if (!m_ready) return false;
562 vmin = m_pMin;
563 vmax = m_pMax;
564 return true;
565}
566
567bool ComponentVoxel::GetElectricFieldRange(double& exmin, double& exmax,
568 double& eymin, double& eymax,
569 double& ezmin, double& ezmax) {
570
571 if (!m_ready) {
572 std::cerr << m_className << "::GetElectricFieldRange:\n";
573 std::cerr << " Field map not available.\n";
574 return false;
575 }
576
577 exmin = exmax = m_efields[0][0][0].fx;
578 eymin = eymax = m_efields[0][0][0].fy;
579 ezmin = ezmax = m_efields[0][0][0].fz;
580 for (unsigned int i = 0; i < m_nX; ++i) {
581 for (unsigned int j = 0; j < m_nY; ++j) {
582 for (unsigned int k = 0; k < m_nZ; ++k) {
583 const Element& element = m_efields[i][j][k];
584 if (element.fx < exmin) exmin = element.fx;
585 if (element.fx > exmax) exmax = element.fx;
586 if (element.fy < eymin) eymin = element.fy;
587 if (element.fy > eymax) eymax = element.fy;
588 if (element.fz < ezmin) ezmin = element.fz;
589 if (element.fz > ezmax) ezmax = element.fz;
590 }
591 }
592 }
593 return true;
594}
595
596void ComponentVoxel::PrintRegions() const {
597
598 // Do not proceed if not properly initialised.
599 if (!m_ready) {
600 std::cerr << m_className << "::PrintRegions:\n";
601 std::cerr << " Field map not yet initialised.\n";
602 return;
603 }
604
605 if (m_media.empty()) {
606 std::cerr << m_className << "::PrintRegions:\n No regions defined.\n";
607 return;
608 }
609
610 std::cout << m_className << "::PrintRegions:\n";
611 std::cout << " Index Medium\n";
612 const unsigned int nMedia = m_media.size();
613 for (unsigned int i = 0; i < nMedia; ++i) {
614 const std::string name = m_media[i] ? m_media[i]->GetName() : "none";
615 std::cout << " " << i << " " << name << "\n";
616 }
617}
618
619void ComponentVoxel::SetMedium(const unsigned int i, Medium* m) {
620
621 if (!m) {
622 std::cerr << m_className << "::SetMedium:\n Null pointer.\n";
623 if (m_media.empty()) return;
624 }
625 if (i >= m_media.size()) m_media.resize(i + 1, NULL);
626 m_media[i] = m;
627}
628
629Medium* ComponentVoxel::GetMedium(const unsigned int i) const {
630
631 if (i > m_media.size()) {
632 std::cerr << m_className << "::GetMedium:\n Index out of range.\n";
633 return NULL;
634 }
635 return m_media[i];
636}
637
638bool ComponentVoxel::GetElement(const double xi, const double yi,
639 const double zi, unsigned int& i,
640 unsigned int& j, unsigned int& k,
641 bool& xMirrored, bool& yMirrored,
642 bool& zMirrored) const {
643
644 if (!m_hasMesh) {
645 std::cerr << m_className << "::GetElement:\n Mesh is not set.\n";
646 return false;
647 }
648
649 // Reduce the point to the basic cell (in case of periodicity) and
650 // check if it is inside the mesh.
651 const double x = Reduce(xi, m_xMin, m_xMax, m_xPeriodic,
652 m_xMirrorPeriodic, xMirrored);
653 if (x < m_xMin || x > m_xMax) return false;
654 const double y = Reduce(yi, m_yMin, m_yMax, m_yPeriodic,
655 m_yMirrorPeriodic, yMirrored);
656 if (y < m_yMin || y > m_yMax) return false;
657 const double z = Reduce(zi, m_zMin, m_zMax, m_zPeriodic,
658 m_zMirrorPeriodic, zMirrored);
659 if (z < m_zMin || z > m_zMax) return false;
660
661 // Get the indices.
662 const double dx = (m_xMax - m_xMin) / m_nX;
663 const double dy = (m_yMax - m_yMin) / m_nY;
664 const double dz = (m_zMax - m_zMin) / m_nZ;
665 i = (unsigned int)((x - m_xMin) / dx);
666 j = (unsigned int)((y - m_yMin) / dy);
667 k = (unsigned int)((z - m_zMin) / dz);
668 if (i >= m_nX) i = m_nX - 1;
669 if (j >= m_nY) j = m_nY - 1;
670 if (k >= m_nZ) k = m_nZ - 1;
671 return true;
672}
673
674bool ComponentVoxel::GetElement(const unsigned int i, const unsigned int j,
675 const unsigned int k, double& v, double& ex,
676 double& ey, double& ez) const {
677
678 v = ex = ey = ez = 0.;
679 if (!m_ready) {
680 if (!m_hasMesh) {
681 std::cerr << m_className << "::GetElement:\n Mesh not set.\n";
682 return false;
683 }
684 std::cerr << m_className << "::GetElement:\n Field map not set.\n";
685 return false;
686 }
687 if (i >= m_nX || j >= m_nY || k >= m_nZ) {
688 std::cerr << m_className << "::GetElement:\n Index out of range.\n";
689 return false;
690 }
691 const Element& element = m_efields[i][j][k];
692 v = element.v;
693 ex = element.fx;
694 ey = element.fy;
695 ez = element.fz;
696 return true;
697}
698
699void ComponentVoxel::Reset() {
700
701 m_efields.clear();
702 m_bfields.clear();
703 m_regions.clear();
704 m_nX = m_nY = m_nZ = 0;
705 m_xMin = m_yMin = m_zMin = 0.;
706 m_xMax = m_yMax = m_zMax = 0.;
707 m_pMin = m_pMax = 0.;
708 m_media.clear();
709
710 m_hasMesh = false;
711 m_hasPotential = false;
712 m_hasEfield = false;
713 m_hasBfield = false;
714 m_ready = false;
715}
716
717void ComponentVoxel::UpdatePeriodicity() {
718
719 if (!m_ready) {
720 std::cerr << m_className << "::UpdatePeriodicity:\n";
721 std::cerr << " Field map not available.\n";
722 return;
723 }
724
725 // Check for conflicts.
726 if (m_xPeriodic && m_xMirrorPeriodic) {
727 std::cerr << m_className << "::UpdatePeriodicity:\n";
728 std::cerr << " Both simple and mirror periodicity\n";
729 std::cerr << " along x requested; reset.\n";
730 m_xPeriodic = m_xMirrorPeriodic = false;
731 }
732
733 if (m_yPeriodic && m_yMirrorPeriodic) {
734 std::cerr << m_className << "::UpdatePeriodicity:\n";
735 std::cerr << " Both simple and mirror periodicity\n";
736 std::cerr << " along y requested; reset.\n";
737 m_yPeriodic = m_yMirrorPeriodic = false;
738 }
739
740 if (m_zPeriodic && m_zMirrorPeriodic) {
741 std::cerr << m_className << "::UpdatePeriodicity:\n";
742 std::cerr << " Both simple and mirror periodicity\n";
743 std::cerr << " along z requested; reset.\n";
744 m_zPeriodic = m_zMirrorPeriodic = false;
745 }
746
747 if (m_xAxiallyPeriodic || m_yAxiallyPeriodic || m_zAxiallyPeriodic) {
748 std::cerr << m_className << "::UpdatePeriodicity:\n";
749 std::cerr << " Axial symmetry is not supported; reset.\n";
750 m_xAxiallyPeriodic = m_yAxiallyPeriodic = m_zAxiallyPeriodic = false;
751 }
752
753 if (m_xRotationSymmetry || m_yRotationSymmetry || m_zRotationSymmetry) {
754 std::cerr << m_className << "::UpdatePeriodicity:\n";
755 std::cerr << " Rotation symmetry is not supported; reset.\n";
756 m_xRotationSymmetry = m_yRotationSymmetry = m_zRotationSymmetry = false;
757 }
758}
759
760double ComponentVoxel::Reduce(const double xin,
761 const double xmin, const double xmax,
762 const bool simplePeriodic,
763 const bool mirrorPeriodic, bool& mirrored) const {
764
765 // In case of periodicity, reduce the coordinate to the basic cell.
766 double x = xin;
767 const double lx = xmax - xmin;
768 if (simplePeriodic) {
769 x = xmin + fmod(x - xmin, lx);
770 if (x < xmin) x += lx;
771 } else if (mirrorPeriodic) {
772 double xNew = xmin + fmod(x - xmin, lx);
773 if (xNew < xmin) xNew += lx;
774 const int nx = int(floor(0.5 + (xNew - x) / lx));
775 if (nx != 2 * (nx / 2)) {
776 xNew = xmin + xmax - xNew;
777 mirrored = true;
778 }
779 x = xNew;
780 }
781 return x;
782}
783
784}
Garfield::ComponentBase
Abstract base class for components.
Definition: ComponentBase.hh:13
Garfield::ComponentBase::m_zMirrorPeriodic
bool m_zMirrorPeriodic
Mirror periodicity in z.
Definition: ComponentBase.hh:257
Garfield::ComponentBase::m_yAxiallyPeriodic
bool m_yAxiallyPeriodic
Axial periodicity in y.
Definition: ComponentBase.hh:261
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:78
Garfield::ComponentBase::m_zRotationSymmetry
bool m_zRotationSymmetry
Rotation symmetry around z-axis.
Definition: ComponentBase.hh:269
Garfield::ComponentBase::m_yRotationSymmetry
bool m_yRotationSymmetry
Rotation symmetry around y-axis.
Definition: ComponentBase.hh:267
Garfield::ComponentBase::m_ready
bool m_ready
Ready for use?
Definition: ComponentBase.hh:239
Garfield::ComponentBase::m_zAxiallyPeriodic
bool m_zAxiallyPeriodic
Axial periodicity in z.
Definition: ComponentBase.hh:263
Garfield::ComponentBase::m_xRotationSymmetry
bool m_xRotationSymmetry
Rotation symmetry around x-axis.
Definition: ComponentBase.hh:265
Garfield::ComponentBase::m_yPeriodic
bool m_yPeriodic
Simple periodicity in y.
Definition: ComponentBase.hh:249
Garfield::ComponentBase::m_yMirrorPeriodic
bool m_yMirrorPeriodic
Mirror periodicity in y.
Definition: ComponentBase.hh:255
Garfield::ComponentBase::m_xPeriodic
bool m_xPeriodic
Simple periodicity in x.
Definition: ComponentBase.hh:247
Garfield::ComponentBase::m_zPeriodic
bool m_zPeriodic
Simple periodicity in z.
Definition: ComponentBase.hh:251
Garfield::ComponentBase::m_className
std::string m_className
Class name.
Definition: ComponentBase.hh:233
Garfield::ComponentBase::m_xAxiallyPeriodic
bool m_xAxiallyPeriodic
Axial periodicity in x.
Definition: ComponentBase.hh:259
Garfield::ComponentBase::m_xMirrorPeriodic
bool m_xMirrorPeriodic
Mirror periodicity in x.
Definition: ComponentBase.hh:253
Garfield::ComponentVoxel::ComponentVoxel
ComponentVoxel()
Constructor.
Definition: ComponentVoxel.cc:12
Garfield::ComponentVoxel::WeightingField
void WeightingField(const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label)
Definition: ComponentVoxel.cc:83
Garfield::ComponentVoxel::SetWeightingFieldOffset
void SetWeightingFieldOffset(const double x, const double y, const double z)
Definition: ComponentVoxel.cc:95
Garfield::ComponentVoxel::GetBoundingBox
bool GetBoundingBox(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
Get the bounding box coordinates.
Definition: ComponentVoxel.cc:529
Garfield::ComponentVoxel::PrintRegions
void PrintRegions() const
Print all regions.
Definition: ComponentVoxel.cc:596
Garfield::ComponentVoxel::GetElectricFieldRange
bool GetElectricFieldRange(double &exmin, double &exmax, double &eymin, double &eymax, double &ezmin, double &ezmax)
Definition: ComponentVoxel.cc:567
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:638
Garfield::ComponentVoxel::GetVoltageRange
bool GetVoltageRange(double &vmin, double &vmax)
Calculate the voltage range [V].
Definition: ComponentVoxel.cc:559
Garfield::ComponentVoxel::SetMedium
void SetMedium(const unsigned int i, Medium *m)
Set the medium in region i.
Definition: ComponentVoxel.cc:619
Garfield::ComponentVoxel::GetMedium
Medium * GetMedium(const double x, const double y, const double z)
Get the medium at a given location (x, y, z).
Definition: ComponentVoxel.cc:129
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)
Definition: ComponentVoxel.cc:33
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:149
Garfield::ComponentVoxel::MagneticField
void MagneticField(const double x, const double y, const double z, double &bx, double &by, double &bz, int &status)
Definition: ComponentVoxel.cc:102
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:183
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:227
Garfield::Medium
Abstract base class for media.
Definition: Medium.hh:11