Garfield++ v1r0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
ComponentTcad3d.cc
Go to the documentation of this file.
1#include <iostream>
2#include <fstream>
3#include <sstream>
4#include <cmath>
5#include <string>
6#include <algorithm>
7
8#include "ComponentTcad3d.hh"
9#include "GarfieldConstants.hh"
10
11namespace Garfield {
12
13ComponentTcad3d::ComponentTcad3d()
14 : ComponentBase(),
15 m_nRegions(0),
16 m_nVertices(0),
17 m_nElements(0),
18 m_lastElement(0) {
19
20 m_className = "ComponentTcad3d";
21
22 m_regions.clear();
23 m_vertices.clear();
24 m_elements.clear();
25
26 for (int i = nMaxVertices; i--;) m_w[i] = 0.;
27}
28
29void ComponentTcad3d::ElectricField(const double xin, const double yin,
30 const double zin, double& ex, double& ey,
31 double& ez, double& p, Medium*& m,
32 int& status) {
33
34 m = 0;
35 // Make sure the field map has been loaded.
36 if (!ready) {
37 std::cerr << m_className << "::ElectricField:\n";
38 std::cerr << " Field map is not available for interpolation.\n";
39 status = -10;
40 return;
41 }
42
43 // Initialise the electric field and potential.
44 ex = ey = ez = p = 0.;
45
46 double x = xin, y = yin, z = zin;
47 // In case of periodicity, reduce to the cell volume.
48 bool xMirrored = false;
49 const double cellsx = m_xMaxBoundingBox - m_xMinBoundingBox;
50 if (xPeriodic) {
51 x = m_xMinBoundingBox + fmod(x - m_xMinBoundingBox, cellsx);
52 if (x < m_xMinBoundingBox) x += cellsx;
53 } else if (xMirrorPeriodic) {
54 double xNew = m_xMinBoundingBox + fmod(x - m_xMinBoundingBox, cellsx);
55 if (xNew < m_xMinBoundingBox) xNew += cellsx;
56 int nx = int(floor(0.5 + (xNew - x) / cellsx));
57 if (nx != 2 * (nx / 2)) {
58 xNew = m_xMinBoundingBox + m_xMaxBoundingBox - xNew;
59 xMirrored = true;
60 }
61 x = xNew;
62 }
63 bool yMirrored = false;
64 const double cellsy = m_yMaxBoundingBox - m_yMinBoundingBox;
65 if (yPeriodic) {
66 y = m_yMinBoundingBox + fmod(y - m_yMinBoundingBox, cellsy);
67 if (y < m_yMinBoundingBox) y += cellsy;
68 } else if (yMirrorPeriodic) {
69 double yNew = m_yMinBoundingBox + fmod(y - m_yMinBoundingBox, cellsy);
70 if (yNew < m_yMinBoundingBox) yNew += cellsy;
71 int ny = int(floor(0.5 + (yNew - y) / cellsy));
72 if (ny != 2 * (ny / 2)) {
73 yNew = m_yMinBoundingBox + m_yMaxBoundingBox - yNew;
74 yMirrored = true;
75 }
76 y = yNew;
77 }
78 bool zMirrored = false;
79 const double cellsz = m_zMaxBoundingBox - m_zMinBoundingBox;
80 if (zPeriodic) {
81 z = m_zMinBoundingBox + fmod(z - m_zMinBoundingBox, cellsz);
82 if (z < m_zMinBoundingBox) z += cellsz;
83 } else if (zMirrorPeriodic) {
84 double zNew = m_zMinBoundingBox + fmod(z - m_zMinBoundingBox, cellsz);
85 if (zNew < m_zMinBoundingBox) zNew += cellsz;
86 int nz = int(floor(0.5 + (zNew - z) / cellsz));
87 if (nz != 2 * (nz / 2)) {
88 zNew = m_zMinBoundingBox + m_zMaxBoundingBox - zNew;
89 zMirrored = true;
90 }
91 z = zNew;
92 }
93
94 // Check if the point is inside the bounding box.
95 if (x < m_xMinBoundingBox || x > m_xMaxBoundingBox || y < m_yMinBoundingBox ||
96 y > m_yMaxBoundingBox || z < m_zMinBoundingBox || z > m_zMaxBoundingBox) {
97 if (debug) {
98 std::cerr << m_className << "::ElectricField:\n";
99 std::cerr << " Point (" << x << ", " << y << ", " << z
100 << ") is outside the bounding box.\n";
101 }
102 status = -11;
103 return;
104 }
105
106 // Assume this will work.
107 status = 0;
108 // Check if the point is still located in the previously found element.
109 int i = m_lastElement;
110 switch (m_elements[i].type) {
111 case 2:
112 if (CheckTriangle(x, y, z, i)) {
113 ex = m_w[0] * m_vertices[m_elements[i].vertex[0]].ex +
114 m_w[1] * m_vertices[m_elements[i].vertex[1]].ex +
115 m_w[2] * m_vertices[m_elements[i].vertex[2]].ex;
116 ey = m_w[0] * m_vertices[m_elements[i].vertex[0]].ey +
117 m_w[1] * m_vertices[m_elements[i].vertex[1]].ey +
118 m_w[2] * m_vertices[m_elements[i].vertex[2]].ey;
119 ez = m_w[0] * m_vertices[m_elements[i].vertex[0]].ez +
120 m_w[1] * m_vertices[m_elements[i].vertex[1]].ez +
121 m_w[2] * m_vertices[m_elements[i].vertex[2]].ez;
122 p = m_w[0] * m_vertices[m_elements[i].vertex[0]].p +
123 m_w[1] * m_vertices[m_elements[i].vertex[1]].p +
124 m_w[2] * m_vertices[m_elements[i].vertex[2]].p;
125 if (xMirrored) ex = -ex;
126 if (yMirrored) ey = -ey;
127 if (zMirrored) ez = -ez;
128 m = m_regions[m_elements[i].region].medium;
129 if (!m_regions[m_elements[i].region].drift || m == 0) status = -5;
130 return;
131 }
132 break;
133 case 5:
134 if (CheckTetrahedron(x, y, z, i)) {
135 ex = m_w[0] * m_vertices[m_elements[i].vertex[0]].ex +
136 m_w[1] * m_vertices[m_elements[i].vertex[1]].ex +
137 m_w[2] * m_vertices[m_elements[i].vertex[2]].ex +
138 m_w[3] * m_vertices[m_elements[i].vertex[3]].ex;
139 ey = m_w[0] * m_vertices[m_elements[i].vertex[0]].ey +
140 m_w[1] * m_vertices[m_elements[i].vertex[1]].ey +
141 m_w[2] * m_vertices[m_elements[i].vertex[2]].ey +
142 m_w[3] * m_vertices[m_elements[i].vertex[3]].ey;
143 ez = m_w[0] * m_vertices[m_elements[i].vertex[0]].ez +
144 m_w[1] * m_vertices[m_elements[i].vertex[1]].ez +
145 m_w[2] * m_vertices[m_elements[i].vertex[2]].ez +
146 m_w[3] * m_vertices[m_elements[i].vertex[3]].ez;
147 p = m_w[0] * m_vertices[m_elements[i].vertex[0]].p +
148 m_w[1] * m_vertices[m_elements[i].vertex[1]].p +
149 m_w[2] * m_vertices[m_elements[i].vertex[2]].p +
150 m_w[3] * m_vertices[m_elements[i].vertex[3]].p;
151 if (xMirrored) ex = -ex;
152 if (yMirrored) ey = -ey;
153 if (zMirrored) ez = -ez;
154 m = m_regions[m_elements[i].region].medium;
155 if (!m_regions[m_elements[i].region].drift || m == 0) status = -5;
156 return;
157 }
158 break;
159 default:
160 std::cerr << m_className << "::ElectricField:\n";
161 std::cerr << " Unknown element type (" << m_elements[i].type << ").\n";
162 status = -11;
163 return;
164 break;
165 }
166
167 // The point is not in the previous element.
168 // We have to loop over all elements.
169 for (i = m_nElements; i--;) {
170 switch (m_elements[i].type) {
171 case 2:
172 if (CheckTriangle(x, y, z, i)) {
173 ex = m_w[0] * m_vertices[m_elements[i].vertex[0]].ex +
174 m_w[1] * m_vertices[m_elements[i].vertex[1]].ex +
175 m_w[2] * m_vertices[m_elements[i].vertex[2]].ex;
176 ey = m_w[0] * m_vertices[m_elements[i].vertex[0]].ey +
177 m_w[1] * m_vertices[m_elements[i].vertex[1]].ey +
178 m_w[2] * m_vertices[m_elements[i].vertex[2]].ey;
179 ez = m_w[0] * m_vertices[m_elements[i].vertex[0]].ez +
180 m_w[1] * m_vertices[m_elements[i].vertex[1]].ez +
181 m_w[2] * m_vertices[m_elements[i].vertex[2]].ez;
182 p = m_w[0] * m_vertices[m_elements[i].vertex[0]].p +
183 m_w[1] * m_vertices[m_elements[i].vertex[1]].p +
184 m_w[2] * m_vertices[m_elements[i].vertex[2]].p;
185 if (xMirrored) ex = -ex;
186 if (yMirrored) ey = -ey;
187 if (zMirrored) ez = -ez;
188 m_lastElement = i;
189 m = m_regions[m_elements[i].region].medium;
190 if (!m_regions[m_elements[i].region].drift || m == 0) status = -5;
191 return;
192 }
193 break;
194 case 5:
195 if (CheckTetrahedron(x, y, z, i)) {
196 ex = m_w[0] * m_vertices[m_elements[i].vertex[0]].ex +
197 m_w[1] * m_vertices[m_elements[i].vertex[1]].ex +
198 m_w[2] * m_vertices[m_elements[i].vertex[2]].ex +
199 m_w[3] * m_vertices[m_elements[i].vertex[3]].ex;
200 ey = m_w[0] * m_vertices[m_elements[i].vertex[0]].ey +
201 m_w[1] * m_vertices[m_elements[i].vertex[1]].ey +
202 m_w[2] * m_vertices[m_elements[i].vertex[2]].ey +
203 m_w[3] * m_vertices[m_elements[i].vertex[3]].ey;
204 ez = m_w[0] * m_vertices[m_elements[i].vertex[0]].ez +
205 m_w[1] * m_vertices[m_elements[i].vertex[1]].ez +
206 m_w[2] * m_vertices[m_elements[i].vertex[2]].ez +
207 m_w[3] * m_vertices[m_elements[i].vertex[3]].ez;
208 p = m_w[0] * m_vertices[m_elements[i].vertex[0]].p +
209 m_w[1] * m_vertices[m_elements[i].vertex[1]].p +
210 m_w[2] * m_vertices[m_elements[i].vertex[2]].p +
211 m_w[3] * m_vertices[m_elements[i].vertex[3]].p;
212 if (xMirrored) ex = -ex;
213 if (yMirrored) ey = -ey;
214 if (zMirrored) ez = -ez;
215 m_lastElement = i;
216 m = m_regions[m_elements[i].region].medium;
217 if (!m_regions[m_elements[i].region].drift || m == 0) status = -5;
218 return;
219 }
220 break;
221 default:
222 std::cerr << m_className << "::ElectricField:\n";
223 std::cerr << " Invalid element type (" << m_elements[i].type
224 << ").\n";
225 status = -11;
226 return;
227 break;
228 }
229 }
230 // Point is outside the mesh.
231 if (debug) {
232 std::cerr << m_className << "::ElectricField:\n";
233 std::cerr << " Point (" << x << ", " << y << ", " << z
234 << ") is outside the mesh.\n";
235 }
236 status = -6;
237 return;
238}
239
240void ComponentTcad3d::ElectricField(const double x, const double y,
241 const double z, double& ex, double& ey,
242 double& ez, Medium*& m, int& status) {
243
244 double v = 0.;
245 ElectricField(x, y, z, ex, ey, ez, v, m, status);
246}
247
248Medium* ComponentTcad3d::GetMedium(const double& xin, const double& yin,
249 const double& zin) {
250
251 // Make sure the field map has been loaded.
252 if (!ready) {
253 std::cerr << m_className << "::GetMedium:\n";
254 std::cerr << " Field map not available for interpolation.\n";
255 return NULL;
256 }
257
258 double x = xin, y = yin, z = zin;
259 // In case of periodicity, reduce to the cell volume.
260 const double cellsx = m_xMaxBoundingBox - m_xMinBoundingBox;
261 if (xPeriodic) {
262 x = m_xMinBoundingBox + fmod(x - m_xMinBoundingBox, cellsx);
263 if (x < m_xMinBoundingBox) x += cellsx;
264 } else if (xMirrorPeriodic) {
265 double xNew = m_xMinBoundingBox + fmod(x - m_xMinBoundingBox, cellsx);
266 if (xNew < m_xMinBoundingBox) xNew += cellsx;
267 int nx = int(floor(0.5 + (xNew - x) / cellsx));
268 if (nx != 2 * (nx / 2)) {
269 xNew = m_xMinBoundingBox + m_xMaxBoundingBox - xNew;
270 }
271 x = xNew;
272 }
273 const double cellsy = m_yMaxBoundingBox - m_yMinBoundingBox;
274 if (yPeriodic) {
275 y = m_yMinBoundingBox + fmod(y - m_yMinBoundingBox, cellsy);
276 if (y < m_yMinBoundingBox) y += cellsy;
277 } else if (yMirrorPeriodic) {
278 double yNew = m_yMinBoundingBox + fmod(y - m_yMinBoundingBox, cellsy);
279 if (yNew < m_yMinBoundingBox) yNew += cellsy;
280 int ny = int(floor(0.5 + (yNew - y) / cellsy));
281 if (ny != 2 * (ny / 2)) {
282 yNew = m_yMinBoundingBox + m_yMaxBoundingBox - yNew;
283 }
284 y = yNew;
285 }
286 const double cellsz = m_zMaxBoundingBox - m_zMinBoundingBox;
287 if (zPeriodic) {
288 z = m_zMinBoundingBox + fmod(z - m_zMinBoundingBox, cellsz);
289 if (z < m_zMinBoundingBox) z += cellsz;
290 } else if (zMirrorPeriodic) {
291 double zNew = m_zMinBoundingBox + fmod(z - m_zMinBoundingBox, cellsz);
292 if (zNew < m_zMinBoundingBox) zNew += cellsz;
293 int nz = int(floor(0.5 + (zNew - z) / cellsz));
294 if (nz != 2 * (nz / 2)) {
295 zNew = m_zMinBoundingBox + m_zMaxBoundingBox - zNew;
296 }
297 z = zNew;
298 }
299
300 // Check if the point is inside the bounding box.
301 if (x < m_xMinBoundingBox || x > m_xMaxBoundingBox || y < m_yMinBoundingBox ||
302 y > m_yMaxBoundingBox || z < m_zMinBoundingBox || z > m_zMaxBoundingBox) {
303 return NULL;
304 }
305
306 // Check if the point is still located in the previous element.
307 int i = m_lastElement;
308 switch (m_elements[i].type) {
309 case 2:
310 if (CheckTriangle(x, y, z, i)) {
311 return m_regions[m_elements[i].region].medium;
312 }
313 break;
314 case 5:
315 if (CheckTetrahedron(x, y, z, i)) {
316 return m_regions[m_elements[i].region].medium;
317 }
318 break;
319 default:
320 std::cerr << m_className << "::GetMedium:\n";
321 std::cerr << " Invalid element type (" << m_elements[i].type << ").\n";
322 return NULL;
323 break;
324 }
325
326 // The point is not in the previous element.
327 // We have to loop over all elements.
328 for (i = m_nElements; i--;) {
329 switch (m_elements[i].type) {
330 case 2:
331 if (CheckTriangle(x, y, z, i)) {
332 m_lastElement = i;
333 return m_regions[m_elements[i].region].medium;
334 }
335 break;
336 case 5:
337 if (CheckTetrahedron(x, y, z, i)) {
338 m_lastElement = i;
339 return m_regions[m_elements[i].region].medium;
340 }
341 break;
342 default:
343 std::cerr << m_className << "::GetMedium:\n";
344 std::cerr << " Invalid element type (" << m_elements[i].type
345 << ").\n";
346 return NULL;
347 break;
348 }
349 }
350 // The point is outside the mesh.
351 return NULL;
352}
353
354bool ComponentTcad3d::Initialise(const std::string gridfilename,
355 const std::string datafilename) {
356
357 ready = false;
358 // Import mesh data from .grd file.
359 if (!LoadGrid(gridfilename)) {
360 std::cerr << m_className << "::Initialise:\n";
361 std::cerr << " Importing mesh data failed.\n";
362 return false;
363 }
364
365 // Import electric field and potential from .dat file.
366 if (!LoadData(datafilename)) {
367 std::cerr << m_className << "::Initialise:\n";
368 std::cerr << " Importing electric field and potential failed.\n";
369 return false;
370 }
371
372 // Find min./max. coordinates and potentials.
373 m_xMaxBoundingBox = m_vertices[m_elements[0].vertex[0]].x;
374 m_yMaxBoundingBox = m_vertices[m_elements[0].vertex[0]].y;
375 m_zMaxBoundingBox = m_vertices[m_elements[0].vertex[0]].z;
376 m_xMinBoundingBox = m_xMaxBoundingBox;
377 m_yMinBoundingBox = m_yMaxBoundingBox;
378 m_zMinBoundingBox = m_zMaxBoundingBox;
379 m_pMax = m_pMin = m_vertices[m_elements[0].vertex[0]].p;
380 for (int i = m_nElements; i--;) {
381 for (int j = 0; j <= m_elements[i].type; ++j) {
382 if (m_vertices[m_elements[i].vertex[j]].x < m_xMinBoundingBox) {
383 m_xMinBoundingBox = m_vertices[m_elements[i].vertex[j]].x;
384 } else if (m_vertices[m_elements[i].vertex[j]].x > m_xMaxBoundingBox) {
385 m_xMaxBoundingBox = m_vertices[m_elements[i].vertex[j]].x;
386 }
387 if (m_vertices[m_elements[i].vertex[j]].y < m_yMinBoundingBox) {
388 m_yMinBoundingBox = m_vertices[m_elements[i].vertex[j]].y;
389 } else if (m_vertices[m_elements[i].vertex[j]].y > m_yMaxBoundingBox) {
390 m_yMaxBoundingBox = m_vertices[m_elements[i].vertex[j]].y;
391 }
392 if (m_vertices[m_elements[i].vertex[j]].z < m_zMinBoundingBox) {
393 m_zMinBoundingBox = m_vertices[m_elements[i].vertex[j]].z;
394 } else if (m_vertices[m_elements[i].vertex[j]].z > m_zMaxBoundingBox) {
395 m_zMaxBoundingBox = m_vertices[m_elements[i].vertex[j]].z;
396 }
397 if (m_vertices[m_elements[i].vertex[j]].p < m_pMin) {
398 m_pMin = m_vertices[m_elements[i].vertex[j]].p;
399 } else if (m_vertices[m_elements[i].vertex[j]].p > m_pMax) {
400 m_pMax = m_vertices[m_elements[i].vertex[j]].p;
401 }
402 }
403 }
404
405 std::cout << m_className << "::Initialise:\n";
406 std::cout << " Bounding box:\n";
407 std::cout << " " << m_xMinBoundingBox << " < x [cm] < "
408 << m_xMaxBoundingBox << "\n";
409 std::cout << " " << m_yMinBoundingBox << " < y [cm] < "
410 << m_yMaxBoundingBox << "\n";
411 std::cout << " " << m_zMinBoundingBox << " < z [cm] < "
412 << m_zMaxBoundingBox << "\n";
413 std::cout << " Voltage range:\n";
414 std::cout << " " << m_pMin << " < V < " << m_pMax << "\n";
415
416 bool ok = true;
417
418 // Count the number of elements belonging to a region.
419 std::vector<int> nElementsRegion;
420 nElementsRegion.resize(m_nRegions);
421 for (int i = m_nRegions; i--;) nElementsRegion[i] = 0;
422
423 // Count the different element shapes.
424 int nTriangles = 0;
425 int nTetrahedra = 0;
426 int nOtherShapes = 0;
427
428 // Check if there are elements which are not part of any region.
429 int nLoose = 0;
430 std::vector<int> looseElements;
431 looseElements.clear();
432
433 // Check if there are degenerate elements.
434 int nDegenerate = 0;
435 std::vector<int> degenerateElements;
436 degenerateElements.clear();
437
438 for (int i = m_nElements; i--;) {
439 if (m_elements[i].type == 2) {
440 ++nTriangles;
441 if (m_elements[i].vertex[0] == m_elements[i].vertex[1] ||
442 m_elements[i].vertex[1] == m_elements[i].vertex[2] ||
443 m_elements[i].vertex[2] == m_elements[i].vertex[0]) {
444 degenerateElements.push_back(i);
445 ++nDegenerate;
446 }
447 } else if (m_elements[i].type == 5) {
448 if (m_elements[i].vertex[0] == m_elements[i].vertex[1] ||
449 m_elements[i].vertex[0] == m_elements[i].vertex[2] ||
450 m_elements[i].vertex[0] == m_elements[i].vertex[3] ||
451 m_elements[i].vertex[1] == m_elements[i].vertex[2] ||
452 m_elements[i].vertex[1] == m_elements[i].vertex[3] ||
453 m_elements[i].vertex[2] == m_elements[i].vertex[3]) {
454 degenerateElements.push_back(i);
455 ++nDegenerate;
456 }
457 ++nTetrahedra;
458 } else {
459 // Other shapes should not occur, since they were excluded in LoadGrid.
460 ++nOtherShapes;
461 }
462 if (m_elements[i].region >= 0 && m_elements[i].region < m_nRegions) {
463 ++nElementsRegion[m_elements[i].region];
464 } else {
465 looseElements.push_back(i);
466 ++nLoose;
467 }
468 }
469
470 if (nDegenerate > 0) {
471 std::cerr << m_className << "::Initialise:\n";
472 std::cerr << " The following elements are degenerate:\n";
473 for (int i = nDegenerate; i--;) {
474 std::cerr << " " << degenerateElements[i] << "\n";
475 }
476 ok = false;
477 }
478
479 if (nLoose > 0) {
480 std::cerr << m_className << "::Initialise:\n";
481 std::cerr << " The following elements are not part of any region:\n";
482 for (int i = nLoose; i--;) {
483 std::cerr << " " << looseElements[i] << "\n";
484 }
485 ok = false;
486 }
487
488 std::cout << m_className << "::Initialise:\n";
489 std::cout << " Number of regions: " << m_nRegions << "\n";
490 for (int i = 0; i < m_nRegions; ++i) {
491 std::cout << " " << i << ": " << m_regions[i].name << ", "
492 << nElementsRegion[i] << " elements\n";
493 }
494
495 std::cout << " Number of elements: " << m_nElements << "\n";
496 if (nTriangles > 0) {
497 std::cout << " " << nTriangles << " triangles\n";
498 }
499 if (nTetrahedra > 0) {
500 std::cout << " " << nTetrahedra << " tetrahedra\n";
501 }
502 if (nOtherShapes > 0) {
503 std::cout << " " << nOtherShapes << " elements of unknown type\n";
504 std::cout << " Program bug!\n";
505 ready = false;
506 Cleanup();
507 return false;
508 }
509 if (debug) {
510 // For each element, print the indices of the constituting vertices.
511 for (int i = 0; i < m_nElements; ++i) {
512 if (m_elements[i].type == 2) {
513 std::cout << " " << i << ": " << m_elements[i].vertex[0] << " "
514 << m_elements[i].vertex[1] << " " << m_elements[i].vertex[2]
515 << " (triangle, region " << m_elements[i].region << ")\n";
516 } else if (m_elements[i].type == 5) {
517 std::cout << " " << i << ": " << m_elements[i].vertex[0] << " "
518 << m_elements[i].vertex[1] << " " << m_elements[i].vertex[2]
519 << " " << m_elements[i].vertex[3] << " (tetrahedron, region "
520 << m_elements[i].region << ")\n";
521 }
522 }
523 }
524
525 std::cout << " Number of vertices: " << m_nVertices << "\n";
526 if (debug) {
527 for (int i = 0; i < m_nVertices; ++i) {
528 std::cout << " " << i << ": (x, y, z) = (" << m_vertices[i].x << ", "
529 << m_vertices[i].y << ", " << m_vertices[i].z
530 << "), V = " << m_vertices[i].p << "\n";
531 }
532 }
533
534 if (!ok) {
535 ready = false;
536 Cleanup();
537 return false;
538 }
539
540 ready = true;
541 UpdatePeriodicity();
542 return true;
543}
544
545bool ComponentTcad3d::GetBoundingBox(double& xmin, double& ymin, double& zmin,
546 double& xmax, double& ymax, double& zmax) {
547
548 if (!ready) return false;
549 xmin = m_xMinBoundingBox;
550 ymin = m_yMinBoundingBox;
551 zmin = m_zMinBoundingBox;
552 xmax = m_xMaxBoundingBox;
553 ymax = m_yMaxBoundingBox;
554 zmax = m_zMaxBoundingBox;
555 if (xPeriodic || xMirrorPeriodic) {
556 xmin = -INFINITY;
557 xmax = +INFINITY;
558 }
559 if (yPeriodic || yMirrorPeriodic) {
560 ymin = -INFINITY;
561 ymax = +INFINITY;
562 }
563 if (zPeriodic || zMirrorPeriodic) {
564 zmin = -INFINITY;
565 zmax = +INFINITY;
566 }
567 return true;
568}
569
570bool ComponentTcad3d::GetVoltageRange(double& vmin, double& vmax) {
571
572 if (!ready) return false;
573 vmin = m_pMin;
574 vmax = m_pMax;
575 return true;
576}
577
578void ComponentTcad3d::PrintRegions() {
579
580 // Do not proceed if not properly initialised.
581 if (!ready) {
582 std::cerr << m_className << "::PrintRegions:\n";
583 std::cerr << " Field map not yet initialised.\n";
584 return;
585 }
586
587 if (m_nRegions < 0) {
588 std::cerr << m_className << "::PrintRegions:\n";
589 std::cerr << " No regions are currently defined.\n";
590 return;
591 }
592
593 std::cout << m_className << "::PrintRegions:\n";
594 std::cout << " Currently " << m_nRegions << " regions are defined.\n";
595 std::cout << " Index Name Medium\n";
596 for (int i = 0; i < m_nRegions; ++i) {
597 std::cout << " " << i << " " << m_regions[i].name;
598 if (m_regions[i].medium == 0) {
599 std::cout << " none ";
600 } else {
601 std::cout << " " << m_regions[i].medium->GetName();
602 }
603 if (m_regions[i].drift) {
604 std::cout << " (active region)\n";
605 } else {
606 std::cout << "\n";
607 }
608 }
609}
610
611void ComponentTcad3d::GetRegion(const int i, std::string& name, bool& active) {
612
613 if (i < 0 || i >= m_nRegions) {
614 std::cerr << m_className << "::GetRegion:\n";
615 std::cerr << " Region " << i << " does not exist.\n";
616 return;
617 }
618 name = m_regions[i].name;
619 active = m_regions[i].drift;
620}
621
622void ComponentTcad3d::SetDriftRegion(const int i) {
623
624 if (i < 0 || i >= m_nRegions) {
625 std::cerr << m_className << "::SetDriftRegion:\n";
626 std::cerr << " Region " << i << " does not exist.\n";
627 return;
628 }
629 m_regions[i].drift = true;
630}
631
632void ComponentTcad3d::UnsetDriftRegion(const int i) {
633
634 if (i < 0 || i >= m_nRegions) {
635 std::cerr << m_className << "::UnsetDriftRegion:\n";
636 std::cerr << " Region " << i << " does not exist.\n";
637 return;
638 }
639 m_regions[i].drift = false;
640}
641
642void ComponentTcad3d::SetMedium(const int i, Medium* medium) {
643
644 if (i < 0 || i >= m_nRegions) {
645 std::cerr << m_className << "::SetMedium:\n";
646 std::cerr << " Region " << i << " does not exist.\n";
647 return;
648 }
649
650 if (medium == 0) {
651 std::cerr << m_className << "::SetMedium:\n";
652 std::cerr << " Medium pointer is null.\n";
653 return;
654 }
655
656 m_regions[i].medium = medium;
657}
658
659bool ComponentTcad3d::GetMedium(const int i, Medium*& m) const {
660
661 if (i < 0 || i >= m_nRegions) {
662 std::cerr << m_className << "::GetMedium:\n";
663 std::cerr << " Region " << i << " does not exist.\n";
664 return false;
665 }
666
667 m = m_regions[i].medium;
668 if (m == 0) return false;
669 return true;
670}
671
672bool ComponentTcad3d::GetElement(const int i, double& vol, double& dmin,
673 double& dmax, int& type) {
674
675 if (i < 0 || i >= m_nElements) {
676 std::cerr << m_className << "::GetElement:\n";
677 std::cerr << " Element index (" << i << ") out of range.\n";
678 return false;
679 }
680
681 type = m_elements[i].type;
682 if (m_elements[i].type == 2) {
683 // Triangle
684 const double vx = (m_vertices[m_elements[i].vertex[1]].y -
685 m_vertices[m_elements[i].vertex[0]].y) *
686 (m_vertices[m_elements[i].vertex[2]].z -
687 m_vertices[m_elements[i].vertex[0]].z) -
688 (m_vertices[m_elements[i].vertex[1]].z -
689 m_vertices[m_elements[i].vertex[0]].z) *
690 (m_vertices[m_elements[i].vertex[2]].y -
691 m_vertices[m_elements[i].vertex[0]].y);
692 const double vy = (m_vertices[m_elements[i].vertex[1]].z -
693 m_vertices[m_elements[i].vertex[0]].z) *
694 (m_vertices[m_elements[i].vertex[2]].x -
695 m_vertices[m_elements[i].vertex[0]].x) -
696 (m_vertices[m_elements[i].vertex[1]].x -
697 m_vertices[m_elements[i].vertex[0]].x) *
698 (m_vertices[m_elements[i].vertex[2]].z -
699 m_vertices[m_elements[i].vertex[0]].z);
700 const double vz = (m_vertices[m_elements[i].vertex[1]].x -
701 m_vertices[m_elements[i].vertex[0]].x) *
702 (m_vertices[m_elements[i].vertex[2]].y -
703 m_vertices[m_elements[i].vertex[0]].y) -
704 (m_vertices[m_elements[i].vertex[1]].y -
705 m_vertices[m_elements[i].vertex[0]].y) *
706 (m_vertices[m_elements[i].vertex[2]].x -
707 m_vertices[m_elements[i].vertex[0]].x);
708 vol = sqrt(vx * vx + vy * vy + vz * vz);
709 const double a = sqrt(pow(m_vertices[m_elements[i].vertex[1]].x -
710 m_vertices[m_elements[i].vertex[0]].x,
711 2) +
712 pow(m_vertices[m_elements[i].vertex[1]].y -
713 m_vertices[m_elements[i].vertex[0]].y,
714 2) +
715 pow(m_vertices[m_elements[i].vertex[1]].z -
716 m_vertices[m_elements[i].vertex[0]].z,
717 2));
718 const double b = sqrt(pow(m_vertices[m_elements[i].vertex[2]].x -
719 m_vertices[m_elements[i].vertex[0]].x,
720 2) +
721 pow(m_vertices[m_elements[i].vertex[2]].y -
722 m_vertices[m_elements[i].vertex[0]].y,
723 2) +
724 pow(m_vertices[m_elements[i].vertex[2]].z -
725 m_vertices[m_elements[i].vertex[0]].z,
726 2));
727 const double c = sqrt(pow(m_vertices[m_elements[i].vertex[1]].x -
728 m_vertices[m_elements[i].vertex[2]].x,
729 2) +
730 pow(m_vertices[m_elements[i].vertex[1]].y -
731 m_vertices[m_elements[i].vertex[2]].y,
732 2) +
733 pow(m_vertices[m_elements[i].vertex[1]].z -
734 m_vertices[m_elements[i].vertex[2]].z,
735 2));
736 dmin = dmax = a;
737 if (b < dmin) dmin = b;
738 if (c < dmin) dmin = c;
739 if (b > dmax) dmax = b;
740 if (c > dmax) dmax = c;
741 } else if (m_elements[i].type == 5) {
742 // Tetrahedron
743 vol = fabs((m_vertices[m_elements[i].vertex[3]].x -
744 m_vertices[m_elements[i].vertex[0]].x) *
745 ((m_vertices[m_elements[i].vertex[1]].y -
746 m_vertices[m_elements[i].vertex[0]].y) *
747 (m_vertices[m_elements[i].vertex[2]].z -
748 m_vertices[m_elements[i].vertex[0]].z) -
749 (m_vertices[m_elements[i].vertex[2]].y -
750 m_vertices[m_elements[i].vertex[0]].y) *
751 (m_vertices[m_elements[i].vertex[1]].z -
752 m_vertices[m_elements[i].vertex[0]].z)) +
753 (m_vertices[m_elements[i].vertex[3]].y -
754 m_vertices[m_elements[i].vertex[0]].y) *
755 ((m_vertices[m_elements[i].vertex[1]].z -
756 m_vertices[m_elements[i].vertex[0]].z) *
757 (m_vertices[m_elements[i].vertex[2]].x -
758 m_vertices[m_elements[i].vertex[0]].x) -
759 (m_vertices[m_elements[i].vertex[2]].z -
760 m_vertices[m_elements[i].vertex[0]].z) *
761 (m_vertices[m_elements[i].vertex[1]].x -
762 m_vertices[m_elements[i].vertex[0]].x)) +
763 (m_vertices[m_elements[i].vertex[3]].z -
764 m_vertices[m_elements[i].vertex[0]].z) *
765 ((m_vertices[m_elements[i].vertex[1]].x -
766 m_vertices[m_elements[i].vertex[0]].x) *
767 (m_vertices[m_elements[i].vertex[2]].y -
768 m_vertices[m_elements[i].vertex[0]].y) -
769 (m_vertices[m_elements[i].vertex[3]].x -
770 m_vertices[m_elements[i].vertex[0]].x) *
771 (m_vertices[m_elements[i].vertex[1]].y -
772 m_vertices[m_elements[i].vertex[0]].y))) /
773 6.;
774 // Loop over all pairs of m_vertices.
775 for (int j = 0; j < nMaxVertices - 1; ++j) {
776 for (int k = j + 1; k < nMaxVertices; ++k) {
777 // Compute distance.
778 const double dist = sqrt(pow(m_vertices[m_elements[i].vertex[j]].x -
779 m_vertices[m_elements[i].vertex[k]].x,
780 2) +
781 pow(m_vertices[m_elements[i].vertex[j]].y -
782 m_vertices[m_elements[i].vertex[k]].y,
783 2) +
784 pow(m_vertices[m_elements[i].vertex[j]].z -
785 m_vertices[m_elements[i].vertex[k]].z,
786 2));
787 if (k == 1) {
788 dmin = dist;
789 dmax = dist;
790 } else {
791 if (dist < dmin) dmin = dist;
792 if (dist > dmax) dmax = dist;
793 }
794 }
795 }
796 } else {
797 std::cerr << m_className << "::GetElement:\n";
798 std::cerr << " Unexpected element type (" << type << ").\n";
799 return false;
800 }
801 return true;
802}
803
804bool ComponentTcad3d::GetElement(const int i, double& vol, double& dmin,
805 double& dmax, int& type, int& node1,
806 int& node2, int& node3, int& node4, int& node5,
807 int& node6, int& node7, int& reg) {
808
809 if (!GetElement(i, vol, dmin, dmax, type)) return false;
810 node1 = m_elements[i].vertex[0];
811 node2 = m_elements[i].vertex[1];
812 node3 = m_elements[i].vertex[2];
813 node4 = m_elements[i].vertex[3];
814 node5 = m_elements[i].vertex[4];
815 node6 = m_elements[i].vertex[5];
816 node7 = m_elements[i].vertex[6];
817 reg = m_elements[i].region;
818 return true;
819}
820
821bool ComponentTcad3d::GetNode(const int i, double& x, double& y, double& z,
822 double& v, double& ex, double& ey, double& ez) {
823
824 if (i < 0 || i >= m_nVertices) {
825 std::cerr << m_className << "::GetNode:\n";
826 std::cerr << " Node index (" << i << ") out of range.\n";
827 return false;
828 }
829
830 x = m_vertices[i].x;
831 y = m_vertices[i].y;
832 z = m_vertices[i].z;
833 v = m_vertices[i].p;
834 ex = m_vertices[i].ex;
835 ey = m_vertices[i].ey;
836 ez = m_vertices[i].ez;
837 return true;
838}
839
840bool ComponentTcad3d::LoadData(const std::string datafilename) {
841
842 std::ifstream datafile;
843 datafile.open(datafilename.c_str(), std::ios::in);
844 if (!datafile) {
845 std::cerr << m_className << "::LoadData:\n";
846 std::cerr << " Could not open file " << datafilename << ".\n";
847 return false;
848 }
849
850 std::string line;
851 std::istringstream data;
852
853 std::vector<bool> isInRegion(m_nVertices);
854 std::vector<int> fillCount(m_nVertices, 0);
855 for (int i = m_nVertices; i--;) {
856 m_vertices[i].p = 0.;
857 m_vertices[i].ex = 0.;
858 m_vertices[i].ey = 0.;
859 m_vertices[i].ez = 0.;
860 m_vertices[i].isShared = false;
861 }
862
863 std::string::size_type pBra, pKet, pEq;
864
865 while (!datafile.fail()) {
866 // Read one line.
867 std::getline(datafile, line);
868 // Strip white space from beginning of line.
869 line.erase(line.begin(),
870 std::find_if(line.begin(), line.end(),
871 not1(std::ptr_fun<int, int>(isspace))));
872 // Find data section.
873 if (line.substr(0, 8) == "function") {
874 // Read type of data set.
875 pEq = line.find('=');
876 if (pEq == std::string::npos) {
877 // No "=" found.
878 std::cerr << m_className << "::LoadData:\n";
879 std::cerr << " Error reading file " << datafilename << ".\n";
880 std::cerr << " Line:\n";
881 std::cerr << " " << line << "\n";
882 datafile.close();
883 Cleanup();
884 return false;
885 }
886 line = line.substr(pEq + 1);
887 std::string dataset;
888 data.str(line);
889 data >> dataset;
890 data.clear();
891 if (dataset == "ElectrostaticPotential") {
892 std::getline(datafile, line);
893 std::getline(datafile, line);
894 std::getline(datafile, line);
895 std::getline(datafile, line);
896 // Get the region name (given in brackets).
897 pBra = line.find('[');
898 pKet = line.find(']');
899 if (pKet < pBra || pBra == std::string::npos ||
900 pKet == std::string::npos) {
901 std::cerr << m_className << "::LoadData:\n";
902 std::cerr << " Error reading file " << datafilename << "\n";
903 std::cerr << " Line:\n";
904 std::cerr << " " << line << "\n";
905 datafile.close();
906 Cleanup();
907 return false;
908 }
909 line = line.substr(pBra + 1, pKet - pBra - 1);
910 std::string name;
911 data.str(line);
912 data >> name;
913 data.clear();
914 // Check if the region name matches one from the mesh file.
915 int index = -1;
916 for (int j = 0; j < m_nRegions; ++j) {
917 if (name == m_regions[j].name) {
918 index = j;
919 break;
920 }
921 }
922 if (index == -1) {
923 std::cerr << m_className << "::LoadData:\n";
924 std::cerr << " Error reading file " << datafilename << "\n";
925 std::cerr << " Unknown region " << name << ".\n";
926 continue;
927 }
928 // Get the number of values.
929 std::getline(datafile, line);
930 pBra = line.find('(');
931 pKet = line.find(')');
932 if (pKet < pBra || pBra == std::string::npos ||
933 pKet == std::string::npos) {
934 std::cerr << m_className << "::LoadData:\n";
935 std::cerr << " Error reading file " << datafilename << "\n";
936 std::cerr << " Line:\n";
937 std::cerr << " " << line << "\n";
938 datafile.close();
939 Cleanup();
940 return false;
941 }
942 line = line.substr(pBra + 1, pKet - pBra - 1);
943 int nValues;
944 data.str(line);
945 data >> nValues;
946 data.clear();
947 // Mark the vertices belonging to this region.
948 for (int j = m_nVertices; j--;) isInRegion[j] = false;
949 for (int j = 0; j < m_nElements; ++j) {
950 if (m_elements[j].region != index) continue;
951 for (int k = 0; k <= m_elements[j].type; ++k) {
952 isInRegion[m_elements[j].vertex[k]] = true;
953 }
954 }
955 int ivertex = 0;
956 double val;
957 for (int j = 0; j < nValues; ++j) {
958 // Read the next value.
959 datafile >> val;
960 // Find the next vertex belonging to the region.
961 while (ivertex < m_nVertices) {
962 if (isInRegion[ivertex]) break;
963 ++ivertex;
964 }
965 // Check if there is a mismatch between the number of m_vertices
966 // and the number of potential values.
967 if (ivertex >= m_nVertices) {
968 std::cerr << m_className << "::LoadData:\n";
969 std::cerr << " Error reading file " << datafilename << "\n";
970 std::cerr << " Dataset has more values than "
971 << "there are vertices in region " << name << "\n";
972 datafile.close();
973 Cleanup();
974 return false;
975 }
976 m_vertices[ivertex].p = val;
977 ++fillCount[ivertex];
978 ++ivertex;
979 }
980 } else if (dataset == "ElectricField") {
981 // Same procedure as for the potential.
982 std::getline(datafile, line);
983 std::getline(datafile, line);
984 std::getline(datafile, line);
985 std::getline(datafile, line);
986 pBra = line.find('[');
987 pKet = line.find(']');
988 if (pKet < pBra || pBra == std::string::npos ||
989 pKet == std::string::npos) {
990 std::cerr << m_className << "::LoadData:\n";
991 std::cerr << " Error reading file " << datafilename << ".\n";
992 std::cerr << " Line:\n";
993 std::cerr << " " << line << "\n";
994 datafile.close();
995 Cleanup();
996 return false;
997 }
998 line = line.substr(pBra + 1, pKet - pBra - 1);
999 std::string name;
1000 data.str(line);
1001 data >> name;
1002 data.clear();
1003 int index = -1;
1004 for (int j = 0; j < m_nRegions; ++j) {
1005 if (name == m_regions[j].name) {
1006 index = j;
1007 break;
1008 }
1009 }
1010 if (index == -1) {
1011 std::cerr << m_className << "::LoadData:\n";
1012 std::cerr << " Error reading file " << datafilename << "\n";
1013 std::cerr << " Unknown region " << name << ".\n";
1014 continue;
1015 }
1016 std::getline(datafile, line);
1017 pBra = line.find('(');
1018 pKet = line.find(')');
1019 if (pKet < pBra || pBra == std::string::npos ||
1020 pKet == std::string::npos) {
1021 std::cerr << m_className << "::LoadData\n";
1022 std::cerr << " Error reading file " << datafilename << "\n";
1023 std::cerr << " Line:\n";
1024 std::cerr << " " << line << "\n";
1025 datafile.close();
1026 Cleanup();
1027 return false;
1028 }
1029 line = line.substr(pBra + 1, pKet - pBra - 1);
1030 int nValues;
1031 data.str(line);
1032 data >> nValues;
1033 data.clear();
1034 // In case of the electric field, there are three values per vertex.
1035 nValues = nValues / 3;
1036 for (int j = m_nVertices; j--;) isInRegion[j] = false;
1037 for (int j = 0; j < m_nElements; ++j) {
1038 if (m_elements[j].region != index) continue;
1039 for (int k = 0; k <= m_elements[j].type; ++k) {
1040 isInRegion[m_elements[j].vertex[k]] = true;
1041 }
1042 }
1043 int ivertex = 0;
1044 double val1, val2, val3;
1045 for (int j = 0; j < nValues; ++j) {
1046 datafile >> val1 >> val2 >> val3;
1047 while (ivertex < m_nVertices) {
1048 if (isInRegion[ivertex]) break;
1049 ++ivertex;
1050 }
1051 if (ivertex >= m_nVertices) {
1052 std::cerr << m_className << "::LoadData\n"
1053 << " Error reading file " << datafilename << "\n"
1054 << " Dataset has more values than"
1055 << " there are vertices in region " << name << ".\n";
1056 datafile.close();
1057 Cleanup();
1058 return false;
1059 }
1060 m_vertices[ivertex].ex = val1;
1061 m_vertices[ivertex].ey = val2;
1062 m_vertices[ivertex].ez = val3;
1063 ++ivertex;
1064 }
1065 }
1066 }
1067 }
1068 if (datafile.fail() && !datafile.eof()) {
1069 std::cerr << m_className << "::LoadData\n";
1070 std::cerr << " Error reading file " << datafilename << "\n";
1071 datafile.close();
1072 Cleanup();
1073 return false;
1074 }
1075
1076 for (int i = m_nVertices; i--;) {
1077 if (fillCount[i] > 1) m_vertices[i].isShared = true;
1078 }
1079
1080 datafile.close();
1081
1082 return true;
1083}
1084
1085bool ComponentTcad3d::LoadGrid(const std::string gridfilename) {
1086
1087 // Open the file containing the mesh description.
1088 std::ifstream gridfile;
1089 gridfile.open(gridfilename.c_str(), std::ios::in);
1090 if (!gridfile) {
1091 std::cerr << m_className << "::LoadGrid:\n";
1092 std::cerr << " Could not open file " << gridfilename << ".\n";
1093 return false;
1094 }
1095
1096 std::string line;
1097 std::istringstream data;
1098
1099 // Delete existing mesh information.
1100 Cleanup();
1101 std::string::size_type pBra, pKet, pEq;
1102 // Count line numbers.
1103 int iLine = 0;
1104
1105 // Get the number of regions.
1106 while (!gridfile.fail()) {
1107 // Read one line.
1108 std::getline(gridfile, line);
1109 ++iLine;
1110 // Strip white space from the beginning of the line.
1111 line.erase(line.begin(), find_if(line.begin(), line.end(),
1112 not1(std::ptr_fun<int, int>(isspace))));
1113 // Find entry 'nb_regions'.
1114 if (line.substr(0, 10) == "nb_regions") {
1115 pEq = line.find('=');
1116 if (pEq == std::string::npos) {
1117 // No "=" sign found.
1118 std::cerr << m_className << "::LoadGrid:\n";
1119 std::cerr << " Could not read number of regions.\n";
1120 Cleanup();
1121 gridfile.close();
1122 return false;
1123 }
1124 line = line.substr(pEq + 1);
1125 data.str(line);
1126 data >> m_nRegions;
1127 data.clear();
1128 break;
1129 }
1130 if (gridfile.fail()) break;
1131 }
1132 if (gridfile.eof()) {
1133 // Reached end of file.
1134 std::cerr << m_className << "::LoadGrid:\n";
1135 std::cerr << " Could not find entry 'nb_regions' in file\n";
1136 std::cerr << " " << gridfilename << ".\n";
1137 Cleanup();
1138 gridfile.close();
1139 return false;
1140 } else if (gridfile.fail()) {
1141 // Error reading from the file.
1142 std::cerr << m_className << "::LoadGrid:\n";
1143 std::cerr << " Error reading file " << gridfilename << " (line " << iLine
1144 << ").\n";
1145 Cleanup();
1146 gridfile.close();
1147 return false;
1148 }
1149 m_regions.resize(m_nRegions);
1150 for (int j = m_nRegions; j--;) {
1151 m_regions[j].name = "";
1152 m_regions[j].drift = false;
1153 m_regions[j].medium = 0;
1154 }
1155
1156 if (debug) {
1157 std::cout << m_className << "::LoadGrid:\n";
1158 std::cout << " Found " << m_nRegions << " regions.\n";
1159 }
1160
1161 // Get the region names.
1162 while (!gridfile.fail()) {
1163 std::getline(gridfile, line);
1164 ++iLine;
1165 line.erase(line.begin(), find_if(line.begin(), line.end(),
1166 not1(std::ptr_fun<int, int>(isspace))));
1167 // Find entry 'regions'.
1168 if (line.substr(0, 7) == "regions") {
1169 // Get region names (given in brackets).
1170 pBra = line.find('[');
1171 pKet = line.find(']');
1172 if (pKet < pBra || pBra == std::string::npos ||
1173 pKet == std::string::npos) {
1174 // No closed brackets [].
1175 std::cerr << m_className << "::LoadGrid:\n";
1176 std::cerr << " Could not read region names.\n";
1177 Cleanup();
1178 gridfile.close();
1179 return false;
1180 }
1181 line = line.substr(pBra + 1, pKet - pBra - 1);
1182 data.str(line);
1183 for (int j = 0; j < m_nRegions; ++j) {
1184 data >> m_regions[j].name;
1185 data.clear();
1186 // Assume by default that all regions are active.
1187 m_regions[j].drift = true;
1188 m_regions[j].medium = 0;
1189 }
1190 break;
1191 }
1192 }
1193 if (gridfile.eof()) {
1194 // Reached end of file.
1195 std::cerr << m_className << "::LoadGrid:\n";
1196 std::cerr << " Could not find entry 'regions' in file\n";
1197 std::cerr << " " << gridfilename << ".\n";
1198 Cleanup();
1199 gridfile.close();
1200 return false;
1201 } else if (gridfile.fail()) {
1202 // Error reading from the file.
1203 std::cerr << m_className << "::LoadGrid:\n";
1204 std::cerr << " Error reading file " << gridfilename << " (line " << iLine
1205 << ").\n";
1206 Cleanup();
1207 gridfile.close();
1208 return false;
1209 }
1210
1211 // Get the vertices.
1212 while (!gridfile.fail()) {
1213 std::getline(gridfile, line);
1214 ++iLine;
1215 line.erase(line.begin(), find_if(line.begin(), line.end(),
1216 not1(std::ptr_fun<int, int>(isspace))));
1217 // Find section 'Vertices'.
1218 if (line.substr(0, 8) == "Vertices") {
1219 // Get number of vertices (given in brackets).
1220 pBra = line.find('(');
1221 pKet = line.find(')');
1222 if (pKet < pBra || pBra == std::string::npos ||
1223 pKet == std::string::npos) {
1224 // No closed brackets [].
1225 std::cerr << m_className << "::LoadGrid:\n";
1226 std::cerr << " Could not read number of vertices.\n";
1227 Cleanup();
1228 gridfile.close();
1229 return false;
1230 }
1231 line = line.substr(pBra + 1, pKet - pBra - 1);
1232 data.str(line);
1233 data >> m_nVertices;
1234 data.clear();
1235 m_vertices.resize(m_nVertices);
1236 // Get the coordinates of this vertex.
1237 for (int j = 0; j < m_nVertices; ++j) {
1238 gridfile >> m_vertices[j].x >> m_vertices[j].y >> m_vertices[j].z;
1239 // Change units from micron to cm.
1240 m_vertices[j].x *= 1.e-4;
1241 m_vertices[j].y *= 1.e-4;
1242 m_vertices[j].z *= 1.e-4;
1243 }
1244 iLine += m_nVertices - 1;
1245 break;
1246 }
1247 }
1248 if (gridfile.eof()) {
1249 std::cerr << m_className << "::LoadGrid:\n";
1250 std::cerr << " Could not find section 'Vertices' in file\n";
1251 std::cerr << " " << gridfilename << ".\n";
1252 Cleanup();
1253 gridfile.close();
1254 return false;
1255 } else if (gridfile.fail()) {
1256 std::cerr << m_className << "::LoadGrid:\n";
1257 std::cerr << " Error reading file " << gridfilename << " (line " << iLine
1258 << ").\n";
1259 Cleanup();
1260 gridfile.close();
1261 return false;
1262 }
1263
1264 // Get the "edges" (lines connecting two vertices).
1265 int nEdges = 0;
1266 // Temporary arrays for storing edge points.
1267 std::vector<int> edgeP1;
1268 std::vector<int> edgeP2;
1269 while (!gridfile.fail()) {
1270 std::getline(gridfile, line);
1271 ++iLine;
1272 line.erase(line.begin(), find_if(line.begin(), line.end(),
1273 not1(std::ptr_fun<int, int>(isspace))));
1274 // Find section 'Edges'.
1275 if (line.substr(0, 5) == "Edges") {
1276 // Get the number of edges (given in brackets).
1277 pBra = line.find('(');
1278 pKet = line.find(')');
1279 if (pKet < pBra || pBra == std::string::npos ||
1280 pKet == std::string::npos) {
1281 // No closed brackets ()
1282 std::cerr << m_className << "::LoadGrid:\n";
1283 std::cerr << " Could not read number of edges.\n";
1284 Cleanup();
1285 gridfile.close();
1286 return false;
1287 }
1288 line = line.substr(pBra + 1, pKet - pBra - 1);
1289 data.str(line);
1290 data >> nEdges;
1291 data.clear();
1292 edgeP1.resize(nEdges);
1293 edgeP2.resize(nEdges);
1294 // Get the indices of the two endpoints.
1295 for (int j = 0; j < nEdges; ++j) {
1296 gridfile >> edgeP1[j] >> edgeP2[j];
1297 }
1298 iLine += nEdges - 1;
1299 break;
1300 }
1301 }
1302 if (gridfile.eof()) {
1303 std::cerr << m_className << "::LoadGrid:\n";
1304 std::cerr << " Could not find section 'Edges' in file\n";
1305 std::cerr << " " << gridfilename << ".\n";
1306 Cleanup();
1307 gridfile.close();
1308 return false;
1309 } else if (gridfile.fail()) {
1310 std::cerr << m_className << "::LoadGrid:\n";
1311 std::cerr << " Error reading file " << gridfilename << " (line " << iLine
1312 << ").\n";
1313 Cleanup();
1314 gridfile.close();
1315 return false;
1316 }
1317
1318 for (int i = nEdges; i--;) {
1319 // Make sure the indices of the edge endpoints are not out of range.
1320 if (edgeP1[i] < 0 || edgeP1[i] >= m_nVertices || edgeP2[i] < 0 ||
1321 edgeP2[i] >= m_nVertices) {
1322 std::cerr << m_className << "::LoadGrid:\n";
1323 std::cerr << " Vertex index of edge " << i << " out of range.\n";
1324 Cleanup();
1325 gridfile.close();
1326 return false;
1327 }
1328 // Make sure the edge is non-degenerate.
1329 if (edgeP1[i] == edgeP2[i]) {
1330 std::cerr << m_className << "::LoadGrid:\n";
1331 std::cerr << " Edge " << i << " is degenerate.\n";
1332 Cleanup();
1333 gridfile.close();
1334 return false;
1335 }
1336 }
1337
1338 // Get the "faces".
1339 int nFaces = 0;
1340 std::vector<face> faces;
1341 faces.clear();
1342
1343 while (!gridfile.fail()) {
1344 std::getline(gridfile, line);
1345 ++iLine;
1346 line.erase(line.begin(), find_if(line.begin(), line.end(),
1347 not1(std::ptr_fun<int, int>(isspace))));
1348 // Find section 'Faces'.
1349 if (line.substr(0, 5) == "Faces") {
1350 // Get the number of faces (given in brackets).
1351 pBra = line.find('(');
1352 pKet = line.find(')');
1353 if (pKet < pBra || pBra == std::string::npos ||
1354 pKet == std::string::npos) {
1355 // No closed brackets ()
1356 std::cerr << m_className << "::LoadGrid:\n";
1357 std::cerr << " Could not read number of faces.\n";
1358 Cleanup();
1359 gridfile.close();
1360 return false;
1361 }
1362 line = line.substr(pBra + 1, pKet - pBra - 1);
1363 data.str(line);
1364 data >> nFaces;
1365 data.clear();
1366 faces.resize(nFaces);
1367 // Get the indices of the edges constituting this face.
1368 for (int j = 0; j < nFaces; ++j) {
1369 gridfile >> faces[j].type;
1370 if (faces[j].type != 3 && faces[j].type != 4) {
1371 std::cerr << m_className << "::LoadGrid:\n";
1372 std::cerr << " Face with index " << j
1373 << " has invalid number of edges (" << faces[j].type
1374 << ").\n";
1375 Cleanup();
1376 gridfile.close();
1377 return false;
1378 }
1379 for (int k = 0; k < faces[j].type; ++k) {
1380 gridfile >> faces[j].edge[k];
1381 }
1382 }
1383 iLine += nFaces - 1;
1384 break;
1385 }
1386 }
1387 if (gridfile.eof()) {
1388 std::cerr << m_className << "::LoadGrid:\n";
1389 std::cerr << " Could not find section 'Faces' in file\n";
1390 std::cerr << " " << gridfilename << ".\n";
1391 Cleanup();
1392 gridfile.close();
1393 return false;
1394 } else if (gridfile.fail()) {
1395 std::cerr << m_className << "::LoadGrid:\n";
1396 std::cerr << " Error reading file " << gridfilename << " (line " << iLine
1397 << ").\n";
1398 Cleanup();
1399 gridfile.close();
1400 return false;
1401 }
1402
1403 // Get the elements.
1404 int edge0, edge1, edge2;
1405 int face0, face1, face2, face3;
1406 int type;
1407 while (!gridfile.fail()) {
1408 std::getline(gridfile, line);
1409 ++iLine;
1410 line.erase(line.begin(), find_if(line.begin(), line.end(),
1411 not1(std::ptr_fun<int, int>(isspace))));
1412 // Find section 'Elements'.
1413 if (line.substr(0, 8) == "Elements") {
1414 // Get number of elements (given in brackets).
1415 pBra = line.find('(');
1416 pKet = line.find(')');
1417 if (pKet < pBra || pBra == std::string::npos ||
1418 pKet == std::string::npos) {
1419 // No closed brackets ().
1420 std::cerr << m_className << "::LoadGrid:\n";
1421 std::cerr << " Could not read number of elements.\n";
1422 Cleanup();
1423 gridfile.close();
1424 return false;
1425 }
1426 line = line.substr(pBra + 1, pKet - pBra - 1);
1427 data.str(line);
1428 data >> m_nElements;
1429 data.clear();
1430 // Resize array of elements.
1431 m_elements.resize(m_nElements);
1432 // Get type and constituting edges of each element.
1433 for (int j = 0; j < m_nElements; ++j) {
1434 ++iLine;
1435 gridfile >> type;
1436 if (type == 2) {
1437 // Triangle
1438 gridfile >> edge0 >> edge1 >> edge2;
1439 // Get the vertices.
1440 // Negative edge index means that the sequence of the two points
1441 // is supposed to be inverted.
1442 // The actual index is then given by "-index - 1".
1443 // For our purposes, the orientation does not matter.
1444 // Make sure the indices are not out of range.
1445 if (edge0 >= nEdges || -edge0 - 1 >= nEdges || edge1 >= nEdges ||
1446 -edge1 - 1 >= nEdges || edge2 >= nEdges || -edge2 - 1 >= nEdges) {
1447 std::cerr << m_className << "::LoadGrid:\n";
1448 std::cerr << " Error reading file " << gridfilename << " (line "
1449 << iLine << ").\n";
1450 std::cerr << " Edge index out of range.\n";
1451 Cleanup();
1452 gridfile.close();
1453 return false;
1454 }
1455 if (edge0 < 0) edge0 = -edge0 - 1;
1456 if (edge1 < 0) edge1 = -edge1 - 1;
1457 m_elements[j].vertex[0] = edgeP1[edge0];
1458 m_elements[j].vertex[1] = edgeP2[edge0];
1459 if (edgeP1[edge1] != m_elements[j].vertex[0] &&
1460 edgeP1[edge1] != m_elements[j].vertex[1]) {
1461 m_elements[j].vertex[2] = edgeP1[edge1];
1462 } else {
1463 m_elements[j].vertex[2] = edgeP2[edge1];
1464 }
1465 } else if (type == 5) {
1466 // Tetrahedron
1467 // Get the faces.
1468 // Negative face index means that the sequence of the edges
1469 // is supposed to be inverted.
1470 // For our purposes, the orientation does not matter.
1471 gridfile >> face0 >> face1 >> face2 >> face3;
1472 // Make sure the face indices are not out of range.
1473 if (face0 >= nFaces || -face0 - 1 >= nFaces || face1 >= nFaces ||
1474 -face1 - 1 >= nFaces || face2 >= nFaces || -face2 - 1 >= nFaces ||
1475 face3 >= nFaces || -face3 - 1 >= nFaces) {
1476 std::cerr << m_className << "::LoadGrid:\n";
1477 std::cerr << " Error reading file " << gridfilename << " (line "
1478 << iLine << ").\n";
1479 std::cerr << " Face index out of range.\n";
1480 Cleanup();
1481 gridfile.close();
1482 return false;
1483 }
1484 if (face0 < 0) face0 = -face0 - 1;
1485 if (face1 < 0) face1 = -face1 - 1;
1486 // Get the edges of the first face.
1487 edge0 = faces[face0].edge[0];
1488 edge1 = faces[face0].edge[1];
1489 edge2 = faces[face0].edge[2];
1490 if (edge0 < 0) edge0 = -edge0 - 1;
1491 if (edge1 < 0) edge1 = -edge1 - 1;
1492 if (edge2 < 0) edge2 = -edge2 - 1;
1493 // Make sure the edge indices are not out of range.
1494 if (edge0 >= nEdges || edge1 >= nEdges || edge2 >= nEdges) {
1495 std::cerr << m_className << "::LoadGrid:\n";
1496 std::cerr << " Error reading file " << gridfilename << "\n";
1497 std::cerr << " Edge index in element " << j
1498 << " out of range.\n";
1499 Cleanup();
1500 gridfile.close();
1501 return false;
1502 }
1503 // Get the first three vertices.
1504 m_elements[j].vertex[0] = edgeP1[edge0];
1505 m_elements[j].vertex[1] = edgeP2[edge0];
1506 if (edgeP1[edge1] != m_elements[j].vertex[0] &&
1507 edgeP1[edge1] != m_elements[j].vertex[1]) {
1508 m_elements[j].vertex[2] = edgeP1[edge1];
1509 } else {
1510 m_elements[j].vertex[2] = edgeP2[edge1];
1511 }
1512 // Get the fourth vertex from face 1.
1513 edge0 = faces[face1].edge[0];
1514 edge1 = faces[face1].edge[1];
1515 edge2 = faces[face1].edge[2];
1516 if (edge0 < 0) edge0 = -edge0 - 1;
1517 if (edge1 < 0) edge1 = -edge1 - 1;
1518 if (edge2 < 0) edge2 = -edge2 - 1;
1519 if (edgeP1[edge0] != m_elements[j].vertex[0] &&
1520 edgeP1[edge0] != m_elements[j].vertex[1] &&
1521 edgeP1[edge0] != m_elements[j].vertex[2]) {
1522 m_elements[j].vertex[3] = edgeP1[edge0];
1523 } else if (edgeP2[edge0] != m_elements[j].vertex[0] &&
1524 edgeP2[edge0] != m_elements[j].vertex[1] &&
1525 edgeP2[edge0] != m_elements[j].vertex[2]) {
1526 m_elements[j].vertex[3] = edgeP2[edge0];
1527 } else if (edgeP1[edge1] != m_elements[j].vertex[0] &&
1528 edgeP1[edge1] != m_elements[j].vertex[1] &&
1529 edgeP1[edge1] != m_elements[j].vertex[2]) {
1530 m_elements[j].vertex[3] = edgeP1[edge1];
1531 } else if (edgeP2[edge1] != m_elements[j].vertex[0] &&
1532 edgeP2[edge1] != m_elements[j].vertex[1] &&
1533 edgeP2[edge1] != m_elements[j].vertex[2]) {
1534 m_elements[j].vertex[3] = edgeP2[edge1];
1535 } else {
1536 std::cerr << m_className << "::LoadGrid:\n";
1537 std::cerr << " Error reading file " << gridfilename << "\n";
1538 std::cerr << " Face 1 of element " << j << " is degenerate.\n";
1539 Cleanup();
1540 gridfile.close();
1541 return false;
1542 }
1543 } else {
1544 // Other element types are not allowed.
1545 std::cerr << m_className << "::LoadGrid:\n"
1546 << " Error reading file " << gridfilename << " (line "
1547 << iLine << ").\n";
1548 if (type == 0 || type == 1) {
1549 std::cerr << " Invalid element type (" << type
1550 << ") for 3d mesh.\n";
1551 } else {
1552 std::cerr << " Element type " << type << " is not supported.\n";
1553 std::cerr << " Remesh with option -t to create only"
1554 << " triangles and tetrahedra.\n";
1555 }
1556 Cleanup();
1557 gridfile.close();
1558 return false;
1559 }
1560 m_elements[j].type = type;
1561 m_elements[j].region = -1;
1562 }
1563 break;
1564 }
1565 }
1566 if (gridfile.eof()) {
1567 std::cerr << m_className << "::LoadGrid:\n";
1568 std::cerr << " Could not find section 'Elements' in file\n";
1569 std::cerr << " " << gridfilename << ".\n";
1570 Cleanup();
1571 gridfile.close();
1572 return false;
1573 } else if (gridfile.fail()) {
1574 std::cerr << m_className << "::LoadGrid:\n";
1575 std::cerr << " Error reading file " << gridfilename << " (line " << iLine
1576 << ").\n";
1577 Cleanup();
1578 gridfile.close();
1579 return false;
1580 }
1581
1582 // Assign regions to elements.
1583 std::string name;
1584 while (!gridfile.fail()) {
1585 std::getline(gridfile, line);
1586 line.erase(line.begin(), find_if(line.begin(), line.end(),
1587 not1(std::ptr_fun<int, int>(isspace))));
1588 // Find section 'Region'.
1589 if (line.substr(0, 6) == "Region") {
1590 // Get region name (given in brackets).
1591 pBra = line.find('(');
1592 pKet = line.find(')');
1593 if (pKet < pBra || pBra == std::string::npos ||
1594 pKet == std::string::npos) {
1595 std::cerr << m_className << "::LoadGrid:\n";
1596 std::cerr << " Could not read region name.\n";
1597 Cleanup();
1598 gridfile.close();
1599 return false;
1600 }
1601 line = line.substr(pBra + 1, pKet - pBra - 1);
1602 data.str(line);
1603 data >> name;
1604 data.clear();
1605 int index = -1;
1606 for (int j = 0; j < m_nRegions; ++j) {
1607 if (name == m_regions[j].name) {
1608 index = j;
1609 break;
1610 }
1611 }
1612 if (index == -1) {
1613 // Specified region name is not in the list.
1614 std::cerr << m_className << "::LoadGrid:\n";
1615 std::cerr << " Error reading file " << gridfilename << ".\n";
1616 std::cerr << " Unknown region " << name << ".\n";
1617 continue;
1618 }
1619 std::getline(gridfile, line);
1620 std::getline(gridfile, line);
1621 pBra = line.find('(');
1622 pKet = line.find(')');
1623 if (pKet < pBra || pBra == std::string::npos ||
1624 pKet == std::string::npos) {
1625 // No closed brackets ().
1626 std::cerr << m_className << "::LoadGrid:\n";
1627 std::cerr << " Error reading file " << gridfilename << ".\n";
1628 std::cerr << " Could not read number of elements in region " << name
1629 << ".\n";
1630 Cleanup();
1631 gridfile.close();
1632 return false;
1633 }
1634 line = line.substr(pBra + 1, pKet - pBra - 1);
1635 int nElementsRegion;
1636 int iElement;
1637 data.str(line);
1638 data >> nElementsRegion;
1639 data.clear();
1640 for (int j = 0; j < nElementsRegion; ++j) {
1641 gridfile >> iElement;
1642 m_elements[iElement].region = index;
1643 }
1644 }
1645 }
1646
1647 gridfile.close();
1648 if (gridfile.fail() && !gridfile.eof()) {
1649 std::cerr << m_className << "::LoadGrid:\n";
1650 std::cerr << " Error reading file " << gridfilename << ".\n";
1651 Cleanup();
1652 return false;
1653 }
1654
1655 return true;
1656}
1657
1658void ComponentTcad3d::Cleanup() {
1659
1660 // Vertices
1661 m_vertices.clear();
1662 m_nVertices = 0;
1663
1664 // Elements
1665 m_elements.clear();
1666 m_nElements = 0;
1667
1668 // Regions
1669 m_regions.clear();
1670 m_nRegions = 0;
1671}
1672
1673bool ComponentTcad3d::CheckTetrahedron(const double x, const double y,
1674 const double z, const int i) {
1675
1676 m_w[0] = m_w[1] = m_w[2] = m_w[3] = 0.;
1677
1678 const double x10 = m_vertices[m_elements[i].vertex[1]].x -
1679 m_vertices[m_elements[i].vertex[0]].x;
1680 const double y10 = m_vertices[m_elements[i].vertex[1]].y -
1681 m_vertices[m_elements[i].vertex[0]].y;
1682 const double z10 = m_vertices[m_elements[i].vertex[1]].z -
1683 m_vertices[m_elements[i].vertex[0]].z;
1684
1685 const double x20 = m_vertices[m_elements[i].vertex[2]].x -
1686 m_vertices[m_elements[i].vertex[0]].x;
1687 const double y20 = m_vertices[m_elements[i].vertex[2]].y -
1688 m_vertices[m_elements[i].vertex[0]].y;
1689 const double z20 = m_vertices[m_elements[i].vertex[2]].z -
1690 m_vertices[m_elements[i].vertex[0]].z;
1691
1692 const double x30 = m_vertices[m_elements[i].vertex[3]].x -
1693 m_vertices[m_elements[i].vertex[0]].x;
1694 const double y30 = m_vertices[m_elements[i].vertex[3]].y -
1695 m_vertices[m_elements[i].vertex[0]].y;
1696 const double z30 = m_vertices[m_elements[i].vertex[3]].z -
1697 m_vertices[m_elements[i].vertex[0]].z;
1698
1699 const double x21 = m_vertices[m_elements[i].vertex[2]].x -
1700 m_vertices[m_elements[i].vertex[1]].x;
1701 const double y21 = m_vertices[m_elements[i].vertex[2]].y -
1702 m_vertices[m_elements[i].vertex[1]].y;
1703 const double z21 = m_vertices[m_elements[i].vertex[2]].z -
1704 m_vertices[m_elements[i].vertex[1]].z;
1705
1706 const double x31 = m_vertices[m_elements[i].vertex[3]].x -
1707 m_vertices[m_elements[i].vertex[1]].x;
1708 const double y31 = m_vertices[m_elements[i].vertex[3]].y -
1709 m_vertices[m_elements[i].vertex[1]].y;
1710 const double z31 = m_vertices[m_elements[i].vertex[3]].z -
1711 m_vertices[m_elements[i].vertex[1]].z;
1712
1713 const double x32 = m_vertices[m_elements[i].vertex[3]].x -
1714 m_vertices[m_elements[i].vertex[2]].x;
1715 const double y32 = m_vertices[m_elements[i].vertex[3]].y -
1716 m_vertices[m_elements[i].vertex[2]].y;
1717 const double z32 = m_vertices[m_elements[i].vertex[3]].z -
1718 m_vertices[m_elements[i].vertex[2]].z;
1719
1720 m_w[0] =
1721 (x - m_vertices[m_elements[i].vertex[1]].x) * (y21 * z31 - y31 * z21) +
1722 (y - m_vertices[m_elements[i].vertex[1]].y) * (z21 * x31 - z31 * x21) +
1723 (z - m_vertices[m_elements[i].vertex[1]].z) * (x21 * y31 - x31 * y21);
1724
1725 m_w[0] /= x10 * (y31 * z21 - y21 * z31) + y10 * (z31 * x21 - z21 * x31) +
1726 z10 * (x31 * y21 - x21 * y31);
1727
1728 if (m_w[0] < 0.) return false;
1729
1730 m_w[1] =
1731 (x - m_vertices[m_elements[i].vertex[2]].x) * (-y20 * z32 + y32 * z20) +
1732 (y - m_vertices[m_elements[i].vertex[2]].y) * (-z20 * x32 + z32 * x20) +
1733 (z - m_vertices[m_elements[i].vertex[2]].z) * (-x20 * y32 + x32 * y20);
1734
1735 m_w[1] /= x21 * (y20 * z32 - y32 * z20) + y21 * (z20 * x32 - z32 * x20) +
1736 z21 * (x20 * y32 - x32 * y20);
1737
1738 if (m_w[1] < 0.) return false;
1739
1740 m_w[2] =
1741 (x - m_vertices[m_elements[i].vertex[3]].x) * (y30 * z31 - y31 * z30) +
1742 (y - m_vertices[m_elements[i].vertex[3]].y) * (z30 * x31 - z31 * x30) +
1743 (z - m_vertices[m_elements[i].vertex[3]].z) * (x30 * y31 - x31 * y30);
1744
1745 m_w[2] /= x32 * (y31 * z30 - y30 * z31) + y32 * (z31 * x30 - z30 * x31) +
1746 z32 * (x31 * y30 - x30 * y31);
1747
1748 if (m_w[2] < 0.) return false;
1749
1750 m_w[3] =
1751 (x - m_vertices[m_elements[i].vertex[0]].x) * (y20 * z10 - y10 * z20) +
1752 (y - m_vertices[m_elements[i].vertex[0]].y) * (z20 * x10 - z10 * x20) +
1753 (z - m_vertices[m_elements[i].vertex[0]].z) * (x20 * y10 - x10 * y20);
1754
1755 m_w[3] /= x30 * (y20 * z10 - y10 * z20) + y30 * (z20 * x10 - z10 * x20) +
1756 z30 * (x20 * y10 - x10 * y20);
1757
1758 if (m_w[3] < 0.) return false;
1759
1760 if (debug) {
1761 // Reconstruct the point from the local coordinates.
1762 const double xr = m_w[0] * m_vertices[m_elements[i].vertex[0]].x +
1763 m_w[1] * m_vertices[m_elements[i].vertex[1]].x +
1764 m_w[2] * m_vertices[m_elements[i].vertex[2]].x +
1765 m_w[3] * m_vertices[m_elements[i].vertex[3]].x;
1766 const double yr = m_w[0] * m_vertices[m_elements[i].vertex[0]].y +
1767 m_w[1] * m_vertices[m_elements[i].vertex[1]].y +
1768 m_w[2] * m_vertices[m_elements[i].vertex[2]].y +
1769 m_w[3] * m_vertices[m_elements[i].vertex[3]].y;
1770 const double zr = m_w[0] * m_vertices[m_elements[i].vertex[0]].z +
1771 m_w[1] * m_vertices[m_elements[i].vertex[1]].z +
1772 m_w[2] * m_vertices[m_elements[i].vertex[2]].z +
1773 m_w[3] * m_vertices[m_elements[i].vertex[3]].z;
1774 std::cout << m_className << "::CheckTetrahedron:\n";
1775 std::cout << " Original coordinates: (" << x << ", " << y << ", "
1776 << z << ")\n";
1777 std::cout << " Local coordinates: (" << m_w[0] << ", " << m_w[1]
1778 << ", " << m_w[2] << ", " << m_w[3] << ")\n";
1779 std::cerr << " Reconstructed coordinates: (" << xr << ", " << yr << ", "
1780 << zr << ")\n";
1781 std::cerr << " Checksum: " << m_w[0] + m_w[1] + m_w[2] + m_w[3] - 1.
1782 << "\n";
1783 }
1784
1785 return true;
1786}
1787
1788bool ComponentTcad3d::CheckTriangle(const double x, const double y,
1789 const double z, const int i) {
1790
1791 const double v1x = m_vertices[m_elements[i].vertex[1]].x -
1792 m_vertices[m_elements[i].vertex[0]].x;
1793 const double v2x = m_vertices[m_elements[i].vertex[2]].x -
1794 m_vertices[m_elements[i].vertex[0]].x;
1795 const double v1y = m_vertices[m_elements[i].vertex[1]].y -
1796 m_vertices[m_elements[i].vertex[0]].y;
1797 const double v2y = m_vertices[m_elements[i].vertex[2]].y -
1798 m_vertices[m_elements[i].vertex[0]].y;
1799 const double v1z = m_vertices[m_elements[i].vertex[1]].z -
1800 m_vertices[m_elements[i].vertex[0]].z;
1801 const double v2z = m_vertices[m_elements[i].vertex[2]].z -
1802 m_vertices[m_elements[i].vertex[0]].z;
1803
1804 // Check whether the point lies in the plane of the triangle.
1805 // Compute the coefficients of the plane equation.
1806 const double a = v1y * v2z - v2y * v1z;
1807 const double b = v1z * v2x - v2z * v1x;
1808 const double c = v1x * v2y - v2x * v1y;
1809 const double d = a * m_vertices[m_elements[i].vertex[0]].x +
1810 b * m_vertices[m_elements[i].vertex[0]].y +
1811 c * m_vertices[m_elements[i].vertex[0]].z;
1812 // Check if the point satisfies the plane equation.
1813 if (a * x + b * y + c * z != d) return false;
1814
1815 // Map (x, y) onto local variables (b, c) such that
1816 // P = A + b * (B - A) + c * (C - A)
1817 // A point P is inside the triangle ABC if b, c > 0 and b + c < 1;
1818 // b, c are also weighting factors for points B, C
1819
1820 m_w[1] = ((x - m_vertices[m_elements[i].vertex[0]].x) * v2y -
1821 (y - m_vertices[m_elements[i].vertex[0]].y) * v2x) /
1822 (v1x * v2y - v1y * v2x);
1823 if (m_w[1] < 0. || m_w[1] > 1.) return false;
1824
1825 m_w[2] = ((m_vertices[m_elements[i].vertex[0]].x - x) * v1y -
1826 (m_vertices[m_elements[i].vertex[0]].y - y) * v1x) /
1827 (v1x * v2y - v1y * v2x);
1828 if (m_w[2] < 0. || m_w[1] + m_w[2] > 1.) return false;
1829
1830 // Weighting factor for point A
1831 m_w[0] = 1. - m_w[1] - m_w[2];
1832
1833 return true;
1834}
1835
1836void ComponentTcad3d::Reset() {
1837
1838 Cleanup();
1839 ready = false;
1840}
1841
1842void ComponentTcad3d::UpdatePeriodicity() {
1843
1844 if (!ready) {
1845 std::cerr << m_className << "::UpdatePeriodicity:\n";
1846 std::cerr << " Field map not available.\n";
1847 return;
1848 }
1849
1850 // Check for conflicts.
1851 if (xPeriodic && xMirrorPeriodic) {
1852 std::cerr << m_className << "::UpdatePeriodicity:\n";
1853 std::cerr << " Both simple and mirror periodicity\n";
1854 std::cerr << " along x requested; reset.\n";
1855 xPeriodic = xMirrorPeriodic = false;
1856 }
1857
1858 if (yPeriodic && yMirrorPeriodic) {
1859 std::cerr << m_className << "::UpdatePeriodicity:\n";
1860 std::cerr << " Both simple and mirror periodicity\n";
1861 std::cerr << " along y requested; reset.\n";
1862 yPeriodic = yMirrorPeriodic = false;
1863 }
1864
1865 if (zPeriodic && zMirrorPeriodic) {
1866 std::cerr << m_className << "::UpdatePeriodicity:\n";
1867 std::cerr << " Both simple and mirror periodicity\n";
1868 std::cerr << " along z requested; reset.\n";
1869 zPeriodic = zMirrorPeriodic = false;
1870 }
1871
1872 if (xAxiallyPeriodic || yAxiallyPeriodic || zAxiallyPeriodic) {
1873 std::cerr << m_className << "::UpdatePeriodicity:\n";
1874 std::cerr << " Axial symmetry is not supported; reset.\n";
1875 xAxiallyPeriodic = yAxiallyPeriodic = zAxiallyPeriodic = false;
1876 }
1877
1878 if (xRotationSymmetry || yRotationSymmetry || zRotationSymmetry) {
1879 std::cerr << m_className << "::UpdatePeriodicity:\n";
1880 std::cerr << " Rotation symmetry is not supported; reset.\n";
1881 xRotationSymmetry = yRotationSymmetry = zRotationSymmetry = false;
1882 }
1883}
1884}
pow
DoubleAc pow(const DoubleAc &f, double p)
Definition: DoubleAc.cpp:336
sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:313
fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:616
Garfield::ComponentTcad3d::GetBoundingBox
bool GetBoundingBox(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
Definition: ComponentTcad3d.cc:545
Garfield::ComponentTcad3d::Initialise
bool Initialise(const std::string gridfilename, const std::string datafilename)
Definition: ComponentTcad3d.cc:354
Garfield::ComponentTcad3d::UnsetDriftRegion
void UnsetDriftRegion(const int ireg)
Definition: ComponentTcad3d.cc:632
Garfield::ComponentTcad3d::SetMedium
void SetMedium(const int ireg, Medium *m)
Definition: ComponentTcad3d.cc:642
Garfield::ComponentTcad3d::GetRegion
void GetRegion(const int ireg, std::string &name, bool &active)
Definition: ComponentTcad3d.cc:611
Garfield::ComponentTcad3d::GetMedium
Medium * GetMedium(const double &x, const double &y, const double &z)
Definition: ComponentTcad3d.cc:248
Garfield::ComponentTcad3d::GetElement
bool GetElement(const int i, double &vol, double &dmin, double &dmax, int &type)
Definition: ComponentTcad3d.cc:672
Garfield::ComponentTcad3d::GetVoltageRange
bool GetVoltageRange(double &vmin, double &vmax)
Definition: ComponentTcad3d.cc:570
Garfield::ComponentTcad3d::GetNode
bool GetNode(const int i, double &x, double &y, double &z, double &v, double &ex, double &ey, double &ez)
Definition: ComponentTcad3d.cc:821
Garfield::ComponentTcad3d::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: ComponentTcad3d.cc:29
Garfield::ComponentTcad3d::SetDriftRegion
void SetDriftRegion(const int ireg)
Definition: ComponentTcad3d.cc:622