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