Garfield++ 3.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
ComponentNeBem3dMap.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/ComponentNeBem3dMap.hh"
9#include "Garfield/Utilities.hh"
10
11namespace Garfield {
12
13ComponentNeBem3dMap::ComponentNeBem3dMap() : ComponentBase() {
14 m_className = "ComponentNeBem3dMap";
15}
16
17void ComponentNeBem3dMap::ElectricField(const double x, const double y,
18 const double z, double& ex, double& ey,
19 double& ez, double& p, Medium*& m,
20 int& status) {
21 if (m_debug) std::cout << m_className << ": In ElectricField\n";
22
23 m = nullptr;
24 // Make sure the field map has been loaded.
25 if (!m_ready) {
26 std::cerr << m_className << "::ElectricField:\n"
27 << " Field map is not available for interpolation.\n";
28 status = -10;
29 return;
30 }
31
32 // Get the mesh element - values of the lowest corner, i,j,k are returned
33 // For trilinear interpolation, we need to use till the i+1, j+1, k+1 corner
34 unsigned int i = 0, j = 0, k = 0;
35 bool xMirrored = false, yMirrored = false, zMirrored = false;
36 if (!GetElement(x, y, z, i, j, k, xMirrored, yMirrored, zMirrored)) {
37 status = -11;
38 return;
39 }
40 status = 0;
41 if (m_debug) {
42 std::cout << "x, y, z: " << x << ", " << y << ", " << z << "\n"
43 << "i, j, k: " << i << ", " << j << ", " << k << std::endl;
44 }
45 { // adjustment, if needed. int adj not needed outside this block
46 int adj = 0;
47 if (i >= m_nX - 1) {
48 i = m_nX - 1; // CHECK! arbitrary adjustment to avoid index overflow
49 adj = 1;
50 }
51 if (j >= m_nY - 1) {
52 j = m_nY - 1; // This should not happen
53 adj = 1;
54 }
55 if (k >= m_nZ - 1) {
56 k = m_nZ - 1;
57 adj = 1;
58 }
59 if (adj) {
60 std::cout << "x, y, z: " << x << ", " << y << ", " << z << "\n"
61 << "adjusted indices:\n "
62 << "i, j, k: " << i << ", " << j << ", " << k << std::endl;
63 }
64 } // adjustment block
65
66 // Get the electric field and potential using values of six corners.
67 Element& element = m_efields[i][j][k];
68 double ex000 = element.fx;
69 double ey000 = element.fy;
70 double ez000 = element.fz;
71 if (xMirrored) ex = -ex;
72 if (yMirrored) ey = -ey;
73 if (zMirrored) ez = -ez;
74 double p000 = element.v;
75 // Get the medium.
76 int region = m_regions[i][j][k];
77 if (region < 0 || region > (int)m_media.size()) {
78 m = nullptr;
79 status = -5;
80 return;
81 }
82 Medium* m000 = m_media[region];
83 if (!m000) status = -5;
84
85 element = m_efields[i + 1][j][k];
86 double ex100 = element.fx;
87 double ey100 = element.fy;
88 double ez100 = element.fz;
89 if (xMirrored) ex = -ex;
90 if (yMirrored) ey = -ey;
91 if (zMirrored) ez = -ez;
92 double p100 = element.v;
93 // Get the medium.
94 region = m_regions[i + 1][j][k];
95 if (region < 0 || region > (int)m_media.size()) {
96 m = nullptr;
97 status = -5;
98 return;
99 }
100 Medium* m100 = m_media[region];
101 if (!m100) status = -5;
102
103 element = m_efields[i][j + 1][k];
104 double ex010 = element.fx;
105 double ey010 = element.fy;
106 double ez010 = element.fz;
107 if (xMirrored) ex = -ex;
108 if (yMirrored) ey = -ey;
109 if (zMirrored) ez = -ez;
110 double p010 = element.v;
111 // Get the medium.
112 region = m_regions[i][j + 1][k];
113 if (region < 0 || region > (int)m_media.size()) {
114 m = nullptr;
115 status = -5;
116 return;
117 }
118 Medium* m010 = m_media[region];
119 if (!m010) status = -5;
120
121 element = m_efields[i][j][k + 1];
122 double ex001 = element.fx;
123 double ey001 = element.fy;
124 double ez001 = element.fz;
125 if (xMirrored) ex = -ex;
126 if (yMirrored) ey = -ey;
127 if (zMirrored) ez = -ez;
128 double p001 = element.v;
129 // Get the medium.
130 region = m_regions[i][j][k + 1];
131 if (region < 0 || region > (int)m_media.size()) {
132 m = nullptr;
133 status = -5;
134 return;
135 }
136 Medium* m001 = m_media[region];
137 if (!m001) status = -5;
138
139 element = m_efields[i + 1][j + 1][k];
140 double ex110 = element.fx;
141 double ey110 = element.fy;
142 double ez110 = element.fz;
143 if (xMirrored) ex = -ex;
144 if (yMirrored) ey = -ey;
145 if (zMirrored) ez = -ez;
146 double p110 = element.v;
147 // Get the medium.
148 region = m_regions[i + 1][j + 1][k];
149 if (region < 0 || region > (int)m_media.size()) {
150 m = nullptr;
151 status = -5;
152 return;
153 }
154 Medium* m110 = m_media[region];
155 if (!m110) status = -5;
156
157 element = m_efields[i + 1][j][k + 1];
158 double ex101 = element.fx;
159 double ey101 = element.fy;
160 double ez101 = element.fz;
161 if (xMirrored) ex = -ex;
162 if (yMirrored) ey = -ey;
163 if (zMirrored) ez = -ez;
164 double p101 = element.v;
165 // Get the medium.
166 region = m_regions[i + 1][j][k + 1];
167 if (region < 0 || region > (int)m_media.size()) {
168 m = nullptr;
169 status = -5;
170 return;
171 }
172 Medium* m101 = m_media[region];
173 if (!m101) status = -5;
174
175 element = m_efields[i][j + 1][k + 1];
176 double ex011 = element.fx;
177 double ey011 = element.fy;
178 double ez011 = element.fz;
179 if (xMirrored) ex = -ex;
180 if (yMirrored) ey = -ey;
181 if (zMirrored) ez = -ez;
182 double p011 = element.v;
183 // Get the medium.
184 region = m_regions[i][j + 1][k + 1];
185 if (region < 0 || region > (int)m_media.size()) {
186 m = nullptr;
187 status = -5;
188 return;
189 }
190 Medium* m011 = m_media[region];
191 if (!m011) status = -5;
192
193 element = m_efields[i + 1][j + 1][k + 1];
194 double ex111 = element.fx;
195 double ey111 = element.fy;
196 double ez111 = element.fz;
197 if (xMirrored) ex = -ex;
198 if (yMirrored) ey = -ey;
199 if (zMirrored) ez = -ez;
200 double p111 = element.v;
201 // Get the medium.
202 region = m_regions[i + 1][j + 1][k + 1];
203 if (region < 0 || region > (int)m_media.size()) {
204 m = nullptr;
205 status = -5;
206 return;
207 }
208 Medium* m111 = m_media[region];
209 if (!m111) status = -5;
210
211 double delx = (m_xMax - m_xMin) / double(m_nX - 1);
212 double x0 = m_xMin + double(i) * delx;
213 double x1 = m_xMin + double(i + 1) * delx;
214 double dely = (m_yMax - m_yMin) / double(m_nY - 1);
215 double y0 = m_yMin + double(j) * dely;
216 double y1 = m_yMin + double(j + 1) * dely;
217 double delz = (m_zMax - m_zMin) / double(m_nZ - 1);
218 double z0 = m_zMin + double(k) * delz;
219 double z1 = m_zMin + double(k + 1) * delz;
220 double dx0 = (x - x0);
221 double dx1 = (x1 - x);
222 double dy0 = (y - y0);
223 double dy1 = (y1 - y);
224 double dz0 = (z - z0);
225 double dz1 = (z1 - z);
226 double xd = (x - x0) / (x1 - x0);
227 if (xd < 0.0) xd = 0.0;
228 if (xd > 1.0) xd = 1.0;
229 double yd = (y - y0) / (y1 - y0);
230 if (yd < 0.0) yd = 0.0;
231 if (yd > 1.0) yd = 1.0;
232 double zd = (z - z0) / (z1 - z0);
233 if (zd < 0.0) zd = 0.0;
234 if (zd > 1.0) zd = 1.0;
235 ex = TriLinInt(xd, yd, zd, ex000, ex100, ex010, ex001, ex110, ex101, ex011,
236 ex111);
237 ey = TriLinInt(xd, yd, zd, ey000, ey100, ey010, ey001, ey110, ey101, ey011,
238 ey111);
239 ez = TriLinInt(xd, yd, zd, ez000, ez100, ez010, ez001, ez110, ez101, ez011,
240 ez111);
241 p = TriLinInt(xd, yd, zd, p000, p100, p010, p001, p110, p101, p011, p111);
242 /*
243 m = either this material or that!
244 Material value is that of the nearest node
245 If there are equidistant nodes, m is the lower value to give drift a chance.
246 */
247 double d000 = sqrt(dx0 * dx0 + dy0 * dy0 + dz0 * dz0);
248 double d100 = sqrt(dx1 * dx1 + dy0 * dy0 + dz0 * dz0);
249 double d010 = sqrt(dx0 * dx0 + dy1 * dy1 + dz0 * dz0);
250 double d001 = sqrt(dx0 * dx0 + dy0 * dy0 + dz1 * dz1);
251 double d110 = sqrt(dx1 * dx1 + dy1 * dy1 + dz0 * dz0);
252 double d101 = sqrt(dx1 * dx1 + dy0 * dy0 + dz1 * dz1);
253 double d011 = sqrt(dx0 * dx0 + dy1 * dy1 + dz1 * dz1);
254 double d111 = sqrt(dx1 * dx1 + dy1 * dy1 + dz1 * dz1);
255
256 // The lower value notion is not implemented yet
257 // At present the algo works benefiting the last match for a condition.
258 // int mindx = 0;
259 double mindst = d000;
260 m = m000;
261 if (d100 <= mindst) {
262 // mindx = 1;
263 mindst = d100;
264 m = m100;
265 } else if (d010 <= mindst) {
266 // mindx = 2;
267 mindst = d010;
268 m = m010;
269 } else if (d001 <= mindst) {
270 // mindx = 3;
271 mindst = d001;
272 m = m001;
273 } else if (d110 <= mindst) {
274 // mindx = 4;
275 mindst = d110;
276 m = m110;
277 } else if (d101 <= mindst) {
278 // mindx = 5;
279 mindst = d101;
280 m = m101;
281 } else if (d011 <= mindst) {
282 // mindx = 6;
283 mindst = d011;
284 m = m011;
285 } else if (d111 <= mindst) {
286 // mindx = 7;
287 mindst = d111;
288 m = m111;
289 }
290
291 if (m_debug) {
292 std::cout << "x, y, z: " << x << ", " << y << ", " << z << "\n"
293 << "i, j, k: " << i << ", " << j << ", " << k << "\n"
294 << "000=> ex, ey, ez, p, m: " << ex000 << ", " << ey000 << ", "
295 << ez000 << ", " << p000 << ", " << m000 << "\n"
296 << "100=> ex, ey, ez, p, m: " << ex100 << ", " << ey100 << ", "
297 << ez100 << ", " << p100 << ", " << m100 << std::endl
298 << "010=> ex, ey, ez, p, m: " << ex010 << ", " << ey010 << ", "
299 << ez010 << ", " << p010 << ", " << m010 << std::endl
300 << "001=> ex, ey, ez, p, m: " << ex001 << ", " << ey001 << ", "
301 << ez001 << ", " << p001 << ", " << m001 << std::endl
302 << "110=> ex, ey, ez, p, m: " << ex110 << ", " << ey110 << ", "
303 << ez110 << ", " << p110 << ", " << m110 << std::endl
304 << "101=> ex, ey, ez, p, m: " << ex101 << ", " << ey101 << ", "
305 << ez101 << ", " << p101 << ", " << m101 << std::endl
306 << "011=> ex, ey, ez, p, m: " << ex011 << ", " << ey011 << ", "
307 << ez011 << ", " << p011 << ", " << m011 << std::endl
308 << "111=> ex, ey, ez, p, m: " << ex111 << ", " << ey111 << ", "
309 << ez111 << ", " << p111 << ", " << m111 << std::endl
310 << "delx, x, x0, x1, dx0, dx1, xd: " << delx << ", " << x << ", "
311 << x0 << ", " << x1 << ", " << dx0 << ", " << dx1 << ", " << xd
312 << std::endl
313 << "dely, y, y0, y1, dy0, dy1, yd: " << dely << ", " << y << ", "
314 << y0 << ", " << y1 << ", " << dy0 << ", " << dy1 << ", " << yd
315 << std::endl
316 << "delz, z, z0, z1, dz0, dz1, zd: " << delz << ", " << z << ", "
317 << z0 << ", " << z1 << ", " << dz0 << ", " << dz1 << ", " << zd
318 << std::endl
319 << "Values after LinInt=> ex, ey, ez, p, m: " << ex << ", " << ey
320 << ", " << ez << ", " << p << ", " << m << std::endl;
321 }
322 /*
323 ex = ex000; // manual override interpolation
324 ey = ey000;
325 ez = ez000;
326 p = p000;
327 m = m000; // override the nearest node material assignment
328 */
329}
330
331void ComponentNeBem3dMap::ElectricField(const double x, const double y,
332 const double z, double& ex, double& ey,
333 double& ez, Medium*& m, int& status) {
334 double v = 0.;
335 ElectricField(x, y, z, ex, ey, ez, v, m, status);
336}
337
338void ComponentNeBem3dMap::WeightingField(const double x, const double y,
339 const double z, double& wx, double& wy,
340 double& wz,
341 const std::string& /*label*/) {
342 int status = 0;
343 Medium* med = nullptr;
344 double v = 0.;
345 const double x1 = x - m_wField_xOffset;
346 const double y1 = y - m_wField_yOffset;
347 const double z1 = z - m_wField_zOffset;
348 ElectricField(x1, y1, z1, wx, wy, wz, v, med, status);
349}
350
351double ComponentNeBem3dMap::WeightingPotential(const double x, const double y,
352 const double z,
353 const std::string& /*label*/) {
354 int status = 0;
355 Medium* med = nullptr;
356 double v = 0.;
357 const double x1 = x - m_wField_xOffset;
358 const double y1 = y - m_wField_yOffset;
359 const double z1 = z - m_wField_zOffset;
360 double wx = 0., wy = 0., wz = 0.;
361 ElectricField(x1, y1, z1, wx, wy, wz, v, med, status);
362 return v;
363}
364
365void ComponentNeBem3dMap::SetWeightingFieldOffset(const double x,
366 const double y,
367 const double z) {
368 m_wField_xOffset = x;
369 m_wField_yOffset = y;
370 m_wField_zOffset = z;
371}
372
373void ComponentNeBem3dMap::MagneticField(const double x, const double y,
374 const double z, double& bx, double& by,
375 double& bz, int& status) {
376 if (!m_hasBfield) {
377 return ComponentBase::MagneticField(x, y, z, bx, by, bz, status);
378 }
379
380 // Get the mesh element.
381 unsigned int i = 0, j = 0, k = 0;
382 bool xMirrored = false, yMirrored = false, zMirrored = false;
383 if (!GetElement(x, y, z, i, j, k, xMirrored, yMirrored, zMirrored)) {
384 status = -11;
385 return;
386 }
387 status = 0;
388 // Get the field.
389 const Element& element = m_bfields[i][j][k];
390 bx = element.fx;
391 by = element.fy;
392 bz = element.fz;
393 if (xMirrored) bx = -bx;
394 if (yMirrored) by = -by;
395 if (zMirrored) bz = -bz;
396}
397
398Medium* ComponentNeBem3dMap::GetMedium(const double x, const double y,
399 const double z) {
400 // Make sure the field map has been loaded.
401 if (!m_ready) {
402 std::cerr << m_className << "::GetMedium:\n"
403 << " Field map is not available for interpolation.\n";
404 return nullptr;
405 }
406
407 unsigned int i, j, k;
408 bool xMirrored, yMirrored, zMirrored;
409 if (!GetElement(x, y, z, i, j, k, xMirrored, yMirrored, zMirrored)) {
410 return nullptr;
411 }
412 const int region = m_regions[i][j][k];
413 if (region < 0 || region > (int)m_media.size()) return nullptr;
414 return m_media[region];
415}
416
417/// Read map information from file.
418bool ComponentNeBem3dMap::LoadMapInfo(
419 const std::string& MapInfoFile, std::string& MapVersion, int& OptMap,
420 int& OptStaggerMap, unsigned int& NbOfXCells, unsigned int& NbOfYCells,
421 unsigned int& NbOfZCells, double& Xmin, double& Xmax, double& Ymin,
422 double& Ymax, double& Zmin, double& Zmax, double& XStagger,
423 double& YStagger, double& ZStagger, std::string& MapDataFile) {
424 std::ifstream infile;
425 infile.open(MapInfoFile.c_str(), std::ios::in);
426 if (!infile) {
427 std::cerr << m_className << "::LoadMapInfo:\n"
428 << " Could not open file " << MapInfoFile << ".\n";
429 return false;
430 }
431
432 std::string line;
433 unsigned int nLines = 0;
434
435 // read version string
436 if (!infile.fail()) {
437 // Read one line.
438 std::getline(infile, line);
439 ++nLines;
440
441 std::istringstream data;
442 data.str(line);
443 data >> MapVersion;
444
445 if (data.fail()) {
446 std::cerr << m_className << "::LoadMapInfo:\n"
447 << " Error reading line " << nLines << ".\n"
448 << " Cannot retrieve MapVersion.\n";
449 return false;
450 }
451 } // version string
452
453 // read OptMap
454 if (!infile.fail()) {
455 // Read one line.
456 std::getline(infile, line);
457 ++nLines;
458
459 std::istringstream data;
460 data.str(line);
461 data >> OptMap;
462
463 if (data.fail()) {
464 std::cerr << m_className << "::LoadMapInfo:\n"
465 << " Error reading line " << nLines << ".\n"
466 << " Cannot retrieve OptMap.\n";
467 return false;
468 }
469 } // OptMap
470
471 // read OptStaggerMap
472 if (!infile.fail()) {
473 // Read one line.
474 std::getline(infile, line);
475 ++nLines;
476
477 std::istringstream data;
478 data.str(line);
479 data >> OptStaggerMap;
480
481 if (data.fail()) {
482 std::cerr << m_className << "::LoadMapInfo:\n"
483 << " Error reading line " << nLines << ".\n"
484 << " Cannot retrieve OptStaggerMap.\n";
485 return false;
486 }
487 } // OptStaggerMap
488
489 // read NbOfXCells
490 if (!infile.fail()) {
491 // Read one line.
492 std::getline(infile, line);
493 ++nLines;
494
495 std::istringstream data;
496 data.str(line);
497 data >> NbOfXCells;
498
499 if (data.fail()) {
500 std::cerr << m_className << "::LoadMapInfo:\n"
501 << " Error reading line " << nLines << ".\n"
502 << " Cannot retrieve NbOfXCells.\n";
503 return false;
504 }
505 } // NbOfXCells
506
507 // read NbOfYCells
508 if (!infile.fail()) {
509 // Read one line.
510 std::getline(infile, line);
511 ++nLines;
512
513 std::istringstream data;
514 data.str(line);
515 data >> NbOfYCells;
516
517 if (data.fail()) {
518 std::cerr << m_className << "::LoadMapInfo:\n"
519 << " Error reading line " << nLines << ".\n"
520 << " Cannot retrieve NbOfYCells.\n";
521 return false;
522 }
523 } // NbOfYCells
524
525 // read NbOfZCells
526 if (!infile.fail()) {
527 // Read one line.
528 std::getline(infile, line);
529 ++nLines;
530
531 std::istringstream data;
532 data.str(line);
533 data >> NbOfZCells;
534
535 if (data.fail()) {
536 std::cerr << m_className << "::LoadMapInfo:\n"
537 << " Error reading line " << nLines << ".\n"
538 << " Cannot retrieve NbOfZCells.\n";
539 return false;
540 }
541 } // NbOfZCells
542
543 // read Xmin, Xmax
544 if (!infile.fail()) {
545 // Read one line.
546 std::getline(infile, line);
547 ++nLines;
548
549 std::istringstream data;
550 data.str(line);
551 data >> Xmin >> Xmax;
552
553 if (data.fail()) {
554 std::cerr << m_className << "::LoadMapInfo:\n"
555 << " Error reading line " << nLines << ".\n"
556 << " Cannot retrieve Xmin, Xmax.\n";
557 return false;
558 }
559 } // Xmin, Xmax
560
561 // read Ymin, Ymax
562 if (!infile.fail()) {
563 // Read one line.
564 std::getline(infile, line);
565 ++nLines;
566
567 std::istringstream data;
568 data.str(line);
569 data >> Ymin >> Ymax;
570
571 if (data.fail()) {
572 std::cerr << m_className << "::LoadMapInfo:\n"
573 << " Error reading line " << nLines << ".\n"
574 << " Cannot retrieve Ymin, Ymax.\n";
575 return false;
576 }
577 } // Ymin, Ymax
578
579 // read Zmin, Zmax
580 if (!infile.fail()) {
581 // Read one line.
582 std::getline(infile, line);
583 ++nLines;
584
585 std::istringstream data;
586 data.str(line);
587 data >> Zmin >> Zmax;
588
589 if (data.fail()) {
590 std::cerr << m_className << "::LoadMapInfo:\n"
591 << " Error reading line " << nLines << ".\n"
592 << " Cannot retrieve Zmin, Zmax.\n";
593 return false;
594 }
595 } // Zmin, Zmax
596
597 // read XStagger
598 if (!infile.fail()) {
599 // Read one line.
600 std::getline(infile, line);
601 ++nLines;
602
603 std::istringstream data;
604 data.str(line);
605 data >> XStagger;
606
607 if (data.fail()) {
608 std::cerr << m_className << "::LoadMapInfo:\n"
609 << " Error reading line " << nLines << ".\n"
610 << " Cannot retrieve XStagger.\n";
611 return false;
612 }
613 } // XStagger
614
615 // read YStagger
616 if (!infile.fail()) {
617 // Read one line.
618 std::getline(infile, line);
619 ++nLines;
620
621 std::istringstream data;
622 data.str(line);
623 data >> YStagger;
624
625 if (data.fail()) {
626 std::cerr << m_className << "::LoadMapInfo:\n"
627 << " Error reading line " << nLines << ".\n"
628 << " Cannot retrieve YStagger.\n";
629 return false;
630 }
631 } // YStagger
632
633 // read ZStagger
634 if (!infile.fail()) {
635 // Read one line.
636 std::getline(infile, line);
637 ++nLines;
638
639 std::istringstream data;
640 data.str(line);
641 data >> ZStagger;
642
643 if (data.fail()) {
644 std::cerr << m_className << "::LoadMapInfo:\n"
645 << " Error reading line " << nLines << ".\n"
646 << " Cannot retrieve ZStagger.\n";
647 return false;
648 }
649 } // ZStagger
650
651 // read MapDataFile
652 if (!infile.fail()) {
653 // Read one line.
654 std::getline(infile, line);
655 ++nLines;
656
657 std::istringstream data;
658 data.str(line);
659 data >> MapDataFile;
660
661 if (data.fail()) {
662 std::cerr << m_className << "::LoadMapInfo:\n"
663 << " Error reading line " << nLines << ".\n"
664 << " Cannot retrieve MapDataFile.\n";
665 return false;
666 }
667 } // MapDataFile
668
669 return true;
670} // end of LoadMapInfo
671
672void ComponentNeBem3dMap::SetMesh(const unsigned int nx, const unsigned int ny,
673 const unsigned int nz, const double xmin,
674 const double xmax, const double ymin,
675 const double ymax, const double zmin,
676 const double zmax) {
677 Reset();
678 if (nx == 0 || ny == 0 || nz == 0) {
679 std::cerr << m_className << "::SetMesh:\n"
680 << " Number of mesh elements must be positive.\n";
681 return;
682 }
683 if (xmin >= xmax) {
684 std::cerr << m_className << "::SetMesh: Invalid x range.\n";
685 return;
686 } else if (ymin >= ymax) {
687 std::cerr << m_className << "::SetMesh: Invalid y range.\n";
688 return;
689 } else if (zmin >= zmax) {
690 std::cerr << m_className << "::SetMesh: Invalid z range.\n";
691 return;
692 }
693 m_nX = nx; // Note that these are the number of nodes
694 m_nY = ny; // Number of elements is one less than the number of nodes
695 m_nZ = nz;
696 m_xMin = xmin;
697 m_yMin = ymin;
698 m_zMin = zmin;
699 m_xMax = xmax;
700 m_yMax = ymax;
701 m_zMax = zmax;
702 m_hasMesh = true;
703}
704
705bool ComponentNeBem3dMap::LoadElectricField(
706 const std::string& filename, const std::string& format,
707 const bool withPotential, const bool withRegion, const double scaleX,
708 const double scaleE, const double scaleP) {
709 m_ready = false;
710 m_hasPotential = m_hasEfield = false;
711 if (!m_hasMesh) {
712 std::cerr << m_className << "::LoadElectricField:\n"
713 << " Mesh is not set. Call SetMesh first.\n";
714 return false;
715 }
716
717 // Set up the grid.
718 m_efields.resize(m_nX);
719 m_regions.resize(m_nX);
720 for (unsigned int i = 0; i < m_nX; ++i) {
721 m_efields[i].resize(m_nY);
722 m_regions[i].resize(m_nY);
723 for (unsigned int j = 0; j < m_nY; ++j) {
724 m_efields[i][j].resize(m_nZ);
725 m_regions[i][j].resize(m_nZ);
726 for (unsigned int k = 0; k < m_nZ; ++k) {
727 m_efields[i][j][k].fx = 0.;
728 m_efields[i][j][k].fy = 0.;
729 m_efields[i][j][k].fz = 0.;
730 m_efields[i][j][k].v = 0.;
731 m_regions[i][j][k] = 0;
732 }
733 }
734 }
735
736 m_pMin = m_pMax = 0.;
737 if (withPotential) {
738 m_pMin = 1.;
739 m_pMax = -1.;
740 }
741 return LoadData(filename, format, withPotential, withRegion, scaleX, scaleE,
742 scaleP, 'e');
743}
744
745bool ComponentNeBem3dMap::LoadMagneticField(const std::string& filename,
746 const std::string& format,
747 const double scaleX,
748 const double scaleB) {
749 m_hasBfield = false;
750 if (!m_hasMesh) {
751 std::cerr << m_className << "::LoadMagneticField:\n"
752 << " Mesh is not set. Call SetMesh first.\n";
753 return false;
754 }
755
756 // Set up the grid.
757 m_bfields.resize(m_nX);
758 for (unsigned int i = 0; i < m_nX; ++i) {
759 m_bfields[i].resize(m_nY);
760 for (unsigned int j = 0; j < m_nY; ++j) {
761 m_bfields[i][j].resize(m_nZ);
762 for (unsigned int k = 0; k < m_nZ; ++k) {
763 m_bfields[i][j][k].fx = 0.;
764 m_bfields[i][j][k].fy = 0.;
765 m_bfields[i][j][k].fz = 0.;
766 m_bfields[i][j][k].v = 0.;
767 }
768 }
769 }
770
771 return LoadData(filename, format, false, false, scaleX, scaleB, 1., 'b');
772}
773
774bool ComponentNeBem3dMap::LoadData(const std::string& filename,
775 std::string format, const bool withPotential,
776 const bool withRegion, const double scaleX,
777 const double scaleF, const double scaleP,
778 const char field) {
779 if (!m_hasMesh) {
780 std::cerr << m_className << "::LoadData: Mesh has not been set.\n";
781 return false;
782 }
783
784 unsigned int nValues = 0;
785 // Keep track of which elements have been read.
786 std::vector<std::vector<std::vector<bool> > > isSet(
787 m_nX,
788 std::vector<std::vector<bool> >(m_nY, std::vector<bool>(m_nZ, false)));
789
790 std::ifstream infile;
791 infile.open(filename.c_str(), std::ios::in);
792 if (!infile) {
793 std::cerr << m_className << "::LoadData:\n"
794 << " Could not open file " << filename << ".\n";
795 return false;
796 }
797
798 std::transform(format.begin(), format.end(), format.begin(), toupper);
799 unsigned int fmt = 0;
800 if (format == "XY") {
801 fmt = 1;
802 } else if (format == "XYZ") {
803 fmt = 2;
804 } else if (format == "IJ") {
805 fmt = 3;
806 } else if (format == "IJK") {
807 fmt = 4;
808 } else if (format == "YXZ") {
809 fmt = 5;
810 } else {
811 std::cerr << m_className << "::LoadData:\n"
812 << " Unkown format (" << format << ").\n";
813 return false;
814 }
815 std::string line;
816 unsigned int nLines = 0;
817 bool bad = false;
818 while (!infile.fail()) {
819 // Read one line.
820 std::getline(infile, line);
821 ++nLines;
822 // Strip white space from beginning of line.
823 ltrim(line);
824 // Skip empty lines.
825 if (line.empty()) continue;
826 // Skip comments.
827 if (line[0] == '#') continue;
828 if (line[0] == '/' && line[1] == '/') continue;
829 unsigned int i = 0;
830 unsigned int j = 0;
831 unsigned int k = 0;
832 double fx = 0.;
833 double fy = 0.;
834 double fz = 0.;
835 double v = 0.;
836 int region = 0;
837 std::istringstream data;
838 data.str(line);
839 if (fmt == 1) {
840 // "XY"
841 double x, y;
842 data >> x >> y;
843 if (data.fail()) {
844 std::cerr << m_className << "::LoadData:\n"
845 << " Error reading line " << nLines << ".\n"
846 << " Cannot retrieve element coordinates.\n";
847 bad = true;
848 break;
849 }
850 x *= scaleX;
851 y *= scaleX;
852 const double z = 0.5 * (m_zMin + m_zMax);
853 bool xMirrored, yMirrored, zMirrored;
854 if (!GetElement(x, y, z, i, j, k, xMirrored, yMirrored, zMirrored)) {
855 std::cerr << m_className << "::LoadData:\n"
856 << " Error reading line " << nLines << ".\n"
857 << " Point is outside mesh.\n";
858 bad = true;
859 break;
860 }
861 } else if (fmt == 2) {
862 // "XYZ"
863 double x, y, z;
864 data >> x >> y >> z;
865 if (data.fail()) {
866 std::cerr << m_className << "::LoadData:\n"
867 << " Error reading line " << nLines << ".\n"
868 << " Cannot retrieve element coordinates.\n";
869 bad = true;
870 break;
871 }
872 x *= scaleX;
873 y *= scaleX;
874 z *= scaleX;
875 bool xMirrored, yMirrored, zMirrored;
876 if (!GetElement(x, y, z, i, j, k, xMirrored, yMirrored, zMirrored)) {
877 std::cerr << m_className << "::LoadData:\n"
878 << " Error reading line " << nLines << ".\n"
879 << " Point is outside mesh.\n";
880 bad = true;
881 break;
882 }
883 } else if (fmt == 3) {
884 // "IJ"
885 k = 0;
886 data >> i >> j;
887 if (data.fail()) {
888 std::cerr << m_className << "::LoadData:\n"
889 << " Error reading line " << nLines << ".\n"
890 << " Cannot retrieve element index.\n";
891 bad = true;
892 break;
893 }
894 } else if (fmt == 4) {
895 // "IJK"
896 data >> i >> j >> k;
897 if (data.fail()) {
898 std::cerr << m_className << "::LoadData:\n"
899 << " Error reading line " << nLines << ".\n"
900 << " Cannot retrieve element index.\n";
901 bad = true;
902 break;
903 }
904 } else if (fmt == 5) {
905 // "YXZ"
906 double x, y, z, temp;
907 data >> y >> x >> temp;
908 z = temp;
909 if (data.fail()) {
910 std::cerr << m_className << "::LoadData:\n"
911 << " Error reading line " << nLines << ".\n"
912 << " Cannot retrieve element coordinates.\n";
913 bad = true;
914 break;
915 }
916 x *= scaleX;
917 y *= scaleX;
918 z *= scaleX;
919 bool xMirrored, yMirrored, zMirrored;
920 if (!GetElement(x, y, z, i, j, k, xMirrored, yMirrored, zMirrored)) {
921 std::cerr << m_className << "::LoadData:\n"
922 << " Error reading line " << nLines << ".\n"
923 << " Point is outside mesh.\n";
924 bad = true;
925 break;
926 }
927 }
928 // Check the indices.
929 if (i >= m_nX || j >= m_nY || k >= m_nZ) {
930 std::cerr << m_className << "::LoadData:\n"
931 << " Error reading line " << nLines << ".\n"
932 << " Index (" << i << ", " << j << ", " << k
933 << ") out of range.\n";
934 continue;
935 }
936 if (isSet[i][j][k]) {
937 std::cerr << m_className << "::LoadData:\n"
938 << " Error reading line " << nLines << ".\n"
939 << " Mesh element (" << i << ", " << j << ", " << k
940 << ") has already been set.\n";
941 continue;
942 }
943 // Get the field values.
944 if (fmt == 1 || fmt == 3) {
945 // Two-dimensional field-map
946 fz = 0.;
947 data >> fx >> fy;
948 } else if (fmt == 5) {
949 double temp;
950 data >> fy >> fx >> temp;
951 fz = temp;
952 } else {
953 data >> fx >> fy >> fz;
954 }
955 if (data.fail()) {
956 std::cerr << m_className << "::LoadData:\n"
957 << " Error reading line " << nLines << ".\n"
958 << " Cannot read field values.\n";
959 bad = true;
960 break;
961 }
962 fx *= scaleF;
963 fy *= scaleF;
964 fz *= scaleF;
965 if (withPotential) {
966 data >> v;
967 if (data.fail()) {
968 std::cerr << m_className << "::LoadData:\n"
969 << " Error reading line " << nLines << ".\n"
970 << " Cannot read potential.\n";
971 bad = true;
972 break;
973 }
974 v *= scaleP;
975 if (m_pMin > m_pMax) {
976 // First value.
977 m_pMin = v;
978 m_pMax = v;
979 } else {
980 if (v < m_pMin) m_pMin = v;
981 if (v > m_pMax) m_pMax = v;
982 }
983 }
984 if (withRegion) {
985 data >> region;
986 if (data.fail()) {
987 std::cerr << m_className << "::LoadData:\n"
988 << " Error reading line " << nLines << ".\n"
989 << " Cannot read region.\n";
990 bad = true;
991 break;
992 }
993 }
994 if (fmt == 1 || fmt == 3) {
995 // Two-dimensional field-map
996 for (unsigned int kk = 0; kk < m_nZ; ++kk) {
997 if (field == 'e') {
998 m_efields[i][j][kk].fx = fx;
999 m_efields[i][j][kk].fy = fy;
1000 m_efields[i][j][kk].fz = fz;
1001 m_efields[i][j][kk].v = v;
1002 m_regions[i][j][kk] = region;
1003 } else if (field == 'b') {
1004 m_bfields[i][j][kk].fx = fx;
1005 m_bfields[i][j][kk].fy = fy;
1006 m_bfields[i][j][kk].fz = fz;
1007 }
1008 isSet[i][j][kk] = true;
1009 }
1010 } else {
1011 if (field == 'e') {
1012 m_efields[i][j][k].fx = fx;
1013 m_efields[i][j][k].fy = fy;
1014 m_efields[i][j][k].fz = fz;
1015 m_efields[i][j][k].v = v;
1016 m_regions[i][j][k] = region;
1017 } else if (field == 'b') {
1018 m_bfields[i][j][k].fx = fx;
1019 m_bfields[i][j][k].fy = fy;
1020 m_bfields[i][j][k].fz = fz;
1021 }
1022 isSet[i][j][k] = true;
1023 }
1024 ++nValues;
1025 }
1026 if (bad) return false;
1027 std::cout << m_className << "::LoadData:\n"
1028 << " Read " << nValues << " values from " << filename << ".\n";
1029 unsigned int nExpected = m_nX * m_nY;
1030 if (fmt == 2 || fmt == 4 || fmt == 5) nExpected *= m_nZ;
1031 if (nExpected != nValues) {
1032 std::cerr << m_className << "::LoadData:\n"
1033 << " Expected " << nExpected << " values.\n";
1034 }
1035 if (field == 'e') {
1036 m_hasEfield = true;
1037 m_ready = true;
1038 if (withPotential) m_hasPotential = true;
1039 } else if (field == 'b') {
1040 m_hasBfield = true;
1041 }
1042 return true;
1043}
1044
1045bool ComponentNeBem3dMap::GetBoundingBox(double& xmin, double& ymin,
1046 double& zmin, double& xmax,
1047 double& ymax, double& zmax) {
1048 if (!m_ready) return false;
1049 if (m_periodic[0] || m_mirrorPeriodic[0]) {
1050 xmin = -INFINITY;
1051 xmax = +INFINITY;
1052 } else {
1053 xmin = m_xMin;
1054 xmax = m_xMax;
1055 }
1056
1057 if (m_periodic[1] || m_mirrorPeriodic[1]) {
1058 ymin = -INFINITY;
1059 ymax = +INFINITY;
1060 } else {
1061 ymin = m_yMin;
1062 ymax = m_yMax;
1063 }
1064
1065 if (m_periodic[2] || m_mirrorPeriodic[2]) {
1066 zmin = -INFINITY;
1067 zmax = +INFINITY;
1068 } else {
1069 zmin = m_zMin;
1070 zmax = m_zMax;
1071 }
1072 return true;
1073}
1074
1075bool ComponentNeBem3dMap::GetVoltageRange(double& vmin, double& vmax) {
1076 if (!m_ready) return false;
1077 vmin = m_pMin;
1078 vmax = m_pMax;
1079 return true;
1080}
1081
1082bool ComponentNeBem3dMap::GetElectricFieldRange(double& exmin, double& exmax,
1083 double& eymin, double& eymax,
1084 double& ezmin, double& ezmax) {
1085 if (!m_ready) {
1086 std::cerr << m_className << "::GetElectricFieldRange:\n";
1087 std::cerr << " Field map not available.\n";
1088 return false;
1089 }
1090
1091 exmin = exmax = m_efields[0][0][0].fx;
1092 eymin = eymax = m_efields[0][0][0].fy;
1093 ezmin = ezmax = m_efields[0][0][0].fz;
1094 for (unsigned int i = 0; i < m_nX; ++i) {
1095 for (unsigned int j = 0; j < m_nY; ++j) {
1096 for (unsigned int k = 0; k < m_nZ; ++k) {
1097 const Element& element = m_efields[i][j][k];
1098 if (element.fx < exmin) exmin = element.fx;
1099 if (element.fx > exmax) exmax = element.fx;
1100 if (element.fy < eymin) eymin = element.fy;
1101 if (element.fy > eymax) eymax = element.fy;
1102 if (element.fz < ezmin) ezmin = element.fz;
1103 if (element.fz > ezmax) ezmax = element.fz;
1104 }
1105 }
1106 }
1107 return true;
1108}
1109
1110void ComponentNeBem3dMap::PrintRegions() const {
1111 // Do not proceed if not properly initialised.
1112 if (!m_ready) {
1113 std::cerr << m_className << "::PrintRegions:\n"
1114 << " Field map not yet initialised.\n";
1115 return;
1116 }
1117
1118 if (m_media.empty()) {
1119 std::cerr << m_className << "::PrintRegions: No regions defined.\n";
1120 return;
1121 }
1122
1123 std::cout << m_className << "::PrintRegions:\n";
1124 std::cout << " Index Medium\n";
1125 const unsigned int nMedia = m_media.size();
1126 for (unsigned int i = 0; i < nMedia; ++i) {
1127 const std::string name = m_media[i] ? m_media[i]->GetName() : "none";
1128 std::cout << " " << i << " " << name << "\n";
1129 }
1130}
1131
1132void ComponentNeBem3dMap::SetMedium(const unsigned int i, Medium* m) {
1133 if (!m) {
1134 std::cerr << m_className << "::SetMedium: Null pointer.\n";
1135 if (m_media.empty()) return;
1136 }
1137 if (i >= m_media.size()) m_media.resize(i + 1, nullptr);
1138 m_media[i] = m;
1139}
1140
1141Medium* ComponentNeBem3dMap::GetMedium(const unsigned int i) const {
1142 if (i > m_media.size()) {
1143 std::cerr << m_className << "::GetMedium: Index out of range.\n";
1144 return nullptr;
1145 }
1146 return m_media[i];
1147}
1148
1149bool ComponentNeBem3dMap::GetElement(const double xi, const double yi,
1150 const double zi, unsigned int& i,
1151 unsigned int& j, unsigned int& k,
1152 bool& xMirrored, bool& yMirrored,
1153 bool& zMirrored) const {
1154 if (!m_hasMesh) {
1155 std::cerr << m_className << "::GetElement: Mesh is not set.\n";
1156 return false;
1157 }
1158
1159 // Reduce the point to the basic cell (in case of periodicity) and
1160 // check if it is inside the mesh.
1161 const double x =
1162 Reduce(xi, m_xMin, m_xMax, m_periodic[0], m_mirrorPeriodic[0], xMirrored);
1163 if (x < m_xMin || x > m_xMax) return false;
1164 const double y =
1165 Reduce(yi, m_yMin, m_yMax, m_periodic[1], m_mirrorPeriodic[1], yMirrored);
1166 if (y < m_yMin || y > m_yMax) return false;
1167 const double z =
1168 Reduce(zi, m_zMin, m_zMax, m_periodic[2], m_mirrorPeriodic[2], zMirrored);
1169 if (z < m_zMin || z > m_zMax) return false;
1170
1171 // Get the indices.
1172 const double dx =
1173 (m_xMax - m_xMin) /
1174 (m_nX - 1); // m_nX is number of nodes, (m_nX-1) number of elements
1175 const double dy = (m_yMax - m_yMin) / (m_nY - 1);
1176 const double dz = (m_zMax - m_zMin) / (m_nZ - 1);
1177 i = (unsigned int)((x - m_xMin) / dx);
1178 j = (unsigned int)((y - m_yMin) / dy);
1179 k = (unsigned int)((z - m_zMin) / dz);
1180 if (i >= m_nX) i = m_nX - 1;
1181 if (j >= m_nY) j = m_nY - 1;
1182 if (k >= m_nZ) k = m_nZ - 1;
1183 if (m_debug) {
1184 std::cout << m_className << ":In GetElement\n"
1185 << "x, y, z: " << x << ", " << y << ", " << z << std::endl
1186 << "m_xMax, m_yMax, m_zMax: " << m_xMax << ", " << m_yMax << ", "
1187 << m_zMax << std::endl
1188 << "m_xMin, m_yMin, m_zMin: " << m_xMin << ", " << m_yMin << ", "
1189 << m_zMin << std::endl
1190 << "m_nX, m_nY, m_nZ: " << m_nX << ", " << m_nY << ", " << m_nZ
1191 << std::endl
1192 << "dx, dy, dz: " << dx << ", " << dy << ", " << dz << std::endl
1193 << "x-m_xMin, y-m_yMin, z-m_zMin: " << x - m_xMin << ", "
1194 << y - m_yMin << ", " << z - m_zMin << std::endl
1195 << "i, j, k: " << i << ", " << j << ", " << k << std::endl
1196 << "End GetElement" << std::endl;
1197 }
1198 return true;
1199}
1200
1201bool ComponentNeBem3dMap::GetElement(const unsigned int i, const unsigned int j,
1202 const unsigned int k, double& v,
1203 double& ex, double& ey, double& ez) const {
1204 v = ex = ey = ez = 0.;
1205 if (!m_ready) {
1206 if (!m_hasMesh) {
1207 std::cerr << m_className << "::GetElement: Mesh not set.\n";
1208 return false;
1209 }
1210 std::cerr << m_className << "::GetElement: Field map not set.\n";
1211 return false;
1212 }
1213 if (i >= m_nX || j >= m_nY || k >= m_nZ) {
1214 std::cerr << m_className << "::GetElement: Index out of range.\n";
1215 return false;
1216 }
1217 const Element& element = m_efields[i][j][k];
1218 v = element.v;
1219 ex = element.fx;
1220 ey = element.fy;
1221 ez = element.fz;
1222 return true;
1223}
1224
1225void ComponentNeBem3dMap::Reset() {
1226 m_efields.clear();
1227 m_bfields.clear();
1228 m_regions.clear();
1229 m_nX = m_nY = m_nZ = 0;
1230 m_xMin = m_yMin = m_zMin = 0.;
1231 m_xMax = m_yMax = m_zMax = 0.;
1232 m_pMin = m_pMax = 0.;
1233 m_media.clear();
1234
1235 m_hasMesh = false;
1236 m_hasPotential = false;
1237 m_hasEfield = false;
1238 m_hasBfield = false;
1239 m_ready = false;
1240}
1241
1242void ComponentNeBem3dMap::UpdatePeriodicity() {
1243 if (!m_ready) {
1244 std::cerr << m_className << "::UpdatePeriodicity:\n"
1245 << " Field map not available.\n";
1246 return;
1247 }
1248
1249 // Check for conflicts.
1250 for (unsigned int i = 0; i < 3; ++i) {
1251 if (m_periodic[i] && m_mirrorPeriodic[i]) {
1252 std::cerr << m_className << "::UpdatePeriodicity:\n"
1253 << " Both simple and mirror periodicity requested. Reset.\n";
1254 m_periodic[i] = m_mirrorPeriodic[i] = false;
1255 }
1256 }
1257
1258 if (m_axiallyPeriodic[0] || m_axiallyPeriodic[1] || m_axiallyPeriodic[2]) {
1259 std::cerr << m_className << "::UpdatePeriodicity:\n"
1260 << " Axial symmetry is not supported. Reset.\n";
1261 m_axiallyPeriodic.fill(false);
1262 }
1263
1264 if (m_rotationSymmetric[0] || m_rotationSymmetric[1] ||
1265 m_rotationSymmetric[2]) {
1266 std::cerr << m_className << "::UpdatePeriodicity:\n"
1267 << " Rotation symmetry is not supported. Reset.\n";
1268 m_rotationSymmetric.fill(false);
1269 }
1270}
1271
1272double ComponentNeBem3dMap::Reduce(const double xin, const double xmin,
1273 const double xmax, const bool simplePeriodic,
1274 const bool mirrorPeriodic,
1275 bool& mirrored) const {
1276 // In case of periodicity, reduce the coordinate to the basic cell.
1277 double x = xin;
1278 const double lx = xmax - xmin;
1279 if (simplePeriodic) {
1280 x = xmin + fmod(x - xmin, lx);
1281 if (x < xmin) x += lx;
1282 } else if (mirrorPeriodic) {
1283 double xNew = xmin + fmod(x - xmin, lx);
1284 if (xNew < xmin) xNew += lx;
1285 const int nx = int(floor(0.5 + (xNew - x) / lx));
1286 if (nx != 2 * (nx / 2)) {
1287 xNew = xmin + xmax - xNew;
1288 mirrored = true;
1289 }
1290 x = xNew;
1291 }
1292 return x;
1293}
1294
1295double ComponentNeBem3dMap::TriLinInt(const double xd, const double yd,
1296 const double zd, const double c000,
1297 const double c100, const double c010,
1298 const double c001, const double c110,
1299 const double c101, const double c011,
1300 const double c111) {
1301 double c00 = c000 * (1.0 - xd) + c100 * xd;
1302 double c10 = c010 * (1.0 - xd) + c110 * xd;
1303 double c01 = c001 * (1.0 - xd) + c101 * xd;
1304 double c11 = c011 * (1.0 - xd) + c111 * xd;
1305 double c0 = c00 * (1.0 - yd) + c10 * yd;
1306 double c1 = c01 * (1.0 - yd) + c11 * yd;
1307 return (c0 * (1.0 - zd) + c1 * zd);
1308}
1309}
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::ComponentNeBem3dMap::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: ComponentNeBem3dMap.cc:1149
Garfield::ComponentNeBem3dMap::SetMedium
void SetMedium(const unsigned int i, Medium *m)
Set the medium in region i.
Definition: ComponentNeBem3dMap.cc:1132
Garfield::ComponentNeBem3dMap::MagneticField
void MagneticField(const double x, const double y, const double z, double &bx, double &by, double &bz, int &status) override
Definition: ComponentNeBem3dMap.cc:373
Garfield::ComponentNeBem3dMap::PrintRegions
void PrintRegions() const
Print all regions.
Definition: ComponentNeBem3dMap.cc:1110
Garfield::ComponentNeBem3dMap::LoadMapInfo
bool LoadMapInfo(const std::string &MapInfoFile, std::string &MapVersion, int &OptMap, int &OptStaggerMap, unsigned int &NbOfXCells, unsigned int &NbOfYCells, unsigned int &NbOfZCells, double &Xmin, double &Xmax, double &Ymin, double &Ymax, double &Zmin, double &Zmax, double &XStagger, double &YStagger, double &ZStagger, std::string &MapDataFile)
Map related information.
Definition: ComponentNeBem3dMap.cc:418
Garfield::ComponentNeBem3dMap::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: ComponentNeBem3dMap.cc:17
Garfield::ComponentNeBem3dMap::GetBoundingBox
bool GetBoundingBox(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax) override
Get the bounding box coordinates.
Definition: ComponentNeBem3dMap.cc:1045
Garfield::ComponentNeBem3dMap::WeightingPotential
double WeightingPotential(const double x, const double y, const double z, const std::string &label) override
Definition: ComponentNeBem3dMap.cc:351
Garfield::ComponentNeBem3dMap::WeightingField
void WeightingField(const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label) override
Definition: ComponentNeBem3dMap.cc:338
Garfield::ComponentNeBem3dMap::GetMedium
Medium * GetMedium(const double x, const double y, const double z) override
Get the medium at a given location (x, y, z).
Definition: ComponentNeBem3dMap.cc:398
Garfield::ComponentNeBem3dMap::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: ComponentNeBem3dMap.cc:745
Garfield::ComponentNeBem3dMap::GetElectricFieldRange
bool GetElectricFieldRange(double &exmin, double &exmax, double &eymin, double &eymax, double &ezmin, double &ezmax)
Definition: ComponentNeBem3dMap.cc:1082
Garfield::ComponentNeBem3dMap::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: ComponentNeBem3dMap.cc:705
Garfield::ComponentNeBem3dMap::GetVoltageRange
bool GetVoltageRange(double &vmin, double &vmax) override
Calculate the voltage range [V].
Definition: ComponentNeBem3dMap.cc:1075
Garfield::ComponentNeBem3dMap::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: ComponentNeBem3dMap.cc:672
Garfield::ComponentNeBem3dMap::SetWeightingFieldOffset
void SetWeightingFieldOffset(const double x, const double y, const double z)
Definition: ComponentNeBem3dMap.cc:365
Garfield::Medium
Abstract base class for media.
Definition: Medium.hh:13
Garfield::ltrim
void ltrim(std::string &line)
Definition: Utilities.hh:10