Garfield++ v1r0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
ComponentAnsys123.cc
Go to the documentation of this file.
1#include <stdio.h>
2#include <string.h>
3#include <iostream>
4#include <fstream>
5#include <stdlib.h>
6#include <math.h>
7
8#include "ComponentAnsys123.hh"
9
10namespace Garfield {
11
12ComponentAnsys123::ComponentAnsys123() : ComponentFieldMap() {
13
14 m_className = "ComponentAnsys123";
15 ready = false;
16}
17
18bool ComponentAnsys123::Initialise(std::string elist, std::string nlist,
19 std::string mplist, std::string prnsol,
20 std::string unit) {
21
22 ready = false;
23 // Keep track of the success.
24 bool ok = true;
25
26 // Buffer for reading
27 const int size = 100;
28 char line[size];
29
30 // Open the material list.
31 std::ifstream fmplist;
32 fmplist.open(mplist.c_str(), std::ios::in);
33 if (fmplist.fail()) {
34 std::cerr << m_className << "::Initialise:\n";
35 std::cerr << " Could not open material file " << mplist
36 << " for reading.\n",
37 std::cerr << " The file perhaps does not exist.\n";
38 return false;
39 }
40
41 // Read the material list.
42 nMaterials = 0;
43 int il = 0, icurrmat = -1;
44 bool readerror = false;
45 while (fmplist.getline(line, size, '\n')) {
46 il++;
47 // Skip page feed
48 if (strcmp(line, "1") == 0) {
49 fmplist.getline(line, size, '\n');
50 il++;
51 fmplist.getline(line, size, '\n');
52 il++;
53 fmplist.getline(line, size, '\n');
54 il++;
55 fmplist.getline(line, size, '\n');
56 il++;
57 fmplist.getline(line, size, '\n');
58 il++;
59 continue;
60 }
61 // Split the line in tokens
62 char* token = NULL;
63 token = strtok(line, " ");
64 // Skip blank lines and headers
65 if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
66 strcmp(token, "TEMPERATURE") == 0 || strcmp(token, "PROPERTY=") == 0 ||
67 int(token[0]) == 10 || int(token[0]) == 13)
68 continue;
69 // Read number of materials,
70 // ensure it does not exceed the maximum and initialise the list
71 if (strcmp(token, "LIST") == 0) {
72 token = strtok(NULL, " ");
73 token = strtok(NULL, " ");
74 token = strtok(NULL, " ");
75 token = strtok(NULL, " ");
76 nMaterials = ReadInteger(token, -1, readerror);
77 if (readerror) {
78 std::cerr << m_className << "::Initialise:\n";
79 std::cerr << " Error reading file " << mplist << " (line " << il
80 << ").\n";
81 fmplist.close();
82 ok = false;
83 return false;
84 }
85 materials.resize(nMaterials);
86 for (int i = 0; i < nMaterials; ++i) {
87 materials[i].ohm = -1;
88 materials[i].eps = -1;
89 materials[i].medium = NULL;
90 }
91 if (debug) {
92 std::cout << m_className << "::Initialise:\n";
93 std::cout << " Number of materials: " << nMaterials << "\n";
94 }
95 } else if (strcmp(token, "MATERIAL") == 0) {
96 // Version 12 format: read material number
97 token = strtok(NULL, " ");
98 token = strtok(NULL, " ");
99 icurrmat = ReadInteger(token, -1, readerror);
100 if (readerror) {
101 std::cerr << m_className << "::Initialise:\n";
102 std::cerr << " Error reading file " << mplist << " (line " << il
103 << ").\n";
104 fmplist.close();
105 ok = false;
106 return false;
107 }
108 } else if (strcmp(token, "TEMP") == 0) {
109 // Version 12 format: read property tag and value
110 token = strtok(NULL, " ");
111 int itype = 0;
112 if (strncmp(token, "PERX", 4) == 0) {
113 itype = 1;
114 } else if (strncmp(token, "RSVX", 4) == 0) {
115 itype = 2;
116 } else {
117 std::cerr << m_className << "::Initialise:\n";
118 std::cerr << " Found unknown material property flag " << token
119 << "\n";
120 std::cerr << " on material properties file " << mplist << " (line "
121 << il << ").\n";
122 ok = false;
123 }
124 fmplist.getline(line, size, '\n');
125 il++;
126 token = NULL;
127 token = strtok(line, " ");
128 if (icurrmat < 1 || icurrmat > nMaterials) {
129 std::cerr << m_className << "::Initialise:\n";
130 std::cerr << " Found out-of-range current material index "
131 << icurrmat << "\n";
132 std::cerr << " in material properties file " << mplist << ".\n";
133 ok = false;
134 readerror = false;
135 } else if (itype == 1) {
136 materials[icurrmat - 1].eps = ReadDouble(token, -1, readerror);
137 } else if (itype == 2) {
138 materials[icurrmat - 1].ohm = ReadDouble(token, -1, readerror);
139 }
140 if (readerror) {
141 std::cerr << m_className << "::Initialise:\n";
142 std::cerr << " Error reading file " << mplist << " line " << il
143 << ").\n";
144 fmplist.close();
145 ok = false;
146 return false;
147 }
148 } else if (strcmp(token, "PROPERTY") == 0) {
149 // Version 11 format
150 token = strtok(NULL, " ");
151 token = strtok(NULL, " ");
152 int itype = 0;
153 if (strcmp(token, "PERX") == 0) {
154 itype = 1;
155 } else if (strcmp(token, "RSVX") == 0) {
156 itype = 2;
157 } else {
158 std::cerr << m_className << "::Initialise:\n";
159 std::cerr << " Found unknown material property flag " << token
160 << "\n";
161 std::cerr << " on material properties file " << mplist << " (line "
162 << il << ").\n";
163 ok = false;
164 }
165 token = strtok(NULL, " ");
166 token = strtok(NULL, " ");
167 int imat = ReadInteger(token, -1, readerror);
168 if (readerror) {
169 std::cerr << m_className << "::Initialise:\n";
170 std::cerr << " Error reading file " << mplist << " (line " << il
171 << ").\n";
172 fmplist.close();
173 ok = false;
174 return false;
175 } else if (imat < 1 || imat > nMaterials) {
176 std::cerr << m_className << "::Initialise:\n";
177 std::cerr << " Found out-of-range material index " << imat << "\n";
178 std::cerr << " in material properties file " << mplist << ".\n";
179 ok = false;
180 } else {
181 fmplist.getline(line, size, '\n');
182 il++;
183 fmplist.getline(line, size, '\n');
184 il++;
185 token = NULL;
186 token = strtok(line, " ");
187 token = strtok(NULL, " ");
188 if (itype == 1) {
189 materials[imat - 1].eps = ReadDouble(token, -1, readerror);
190 } else if (itype == 2) {
191 materials[imat - 1].ohm = ReadDouble(token, -1, readerror);
192 }
193 if (readerror) {
194 std::cerr << m_className << "::Initialise:\n";
195 std::cerr << " Error reading file " << mplist << " (line " << il
196 << ").\n";
197 fmplist.close();
198 ok = false;
199 return false;
200 }
201 }
202 }
203 }
204
205 // Close the file
206 fmplist.close();
207
208 // Find the lowest epsilon, check for eps = 0, set default drift media
209 double epsmin = -1;
210 int iepsmin = -1;
211 for (int imat = 0; imat < nMaterials; ++imat) {
212 if (materials[imat].eps < 0) continue;
213 if (materials[imat].eps == 0) {
214 std::cerr << m_className << "::Initialise:\n";
215 std::cerr << " Material " << imat
216 << " has been assigned a permittivity\n";
217 std::cerr << " equal to zero in " << mplist << ".\n";
218 ok = false;
219 } else if (iepsmin < 0 || epsmin > materials[imat].eps) {
220 epsmin = materials[imat].eps;
221 iepsmin = imat;
222 }
223 }
224
225 if (iepsmin < 0) {
226 std::cerr << m_className << "::Initialise:\n";
227 std::cerr << " No material with positive permittivity found \n";
228 std::cerr << " in material list " << mplist << ".\n";
229 ok = false;
230 } else {
231 for (int imat = 0; imat < nMaterials; ++imat) {
232 if (imat == iepsmin) {
233 materials[imat].driftmedium = true;
234 } else {
235 materials[imat].driftmedium = false;
236 }
237 }
238 }
239
240 // Tell how many lines read
241 std::cout << m_className << "::Initialise:\n";
242 std::cout << " Read properties of " << nMaterials
243 << " materials from file " << mplist << ".\n";
244 if (debug) PrintMaterials();
245
246 // Open the element list
247 std::ifstream felist;
248 felist.open(elist.c_str(), std::ios::in);
249 if (felist.fail()) {
250 std::cerr << m_className << "::Initialise:\n";
251 std::cerr << " Could not open element file " << elist
252 << " for reading.\n";
253 std::cerr << " The file perhaps does not exist.\n";
254 return false;
255 }
256
257 // Read the element list
258 elements.clear();
259 nElements = 0;
260 element newElement;
261 int nbackground = 0;
262 il = 0;
263 int highestnode = 0;
264 while (felist.getline(line, size, '\n')) {
265 il++;
266 // Skip page feed
267 if (strcmp(line, "1") == 0) {
268 felist.getline(line, size, '\n');
269 il++;
270 felist.getline(line, size, '\n');
271 il++;
272 felist.getline(line, size, '\n');
273 il++;
274 felist.getline(line, size, '\n');
275 il++;
276 felist.getline(line, size, '\n');
277 il++;
278 continue;
279 }
280 // Split the line in tokens
281 char* token = NULL;
282 // Split into tokens
283 token = strtok(line, " ");
284 // Skip blank lines and headers
285 if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
286 int(token[0]) == 10 || int(token[0]) == 13 ||
287 strcmp(token, "LIST") == 0 || strcmp(token, "ELEM") == 0) {
288 continue;
289 }
290 // Read the element
291 int ielem = ReadInteger(token, -1, readerror);
292 token = strtok(NULL, " ");
293 int imat = ReadInteger(token, -1, readerror);
294 token = strtok(NULL, " ");
295 token = strtok(NULL, " ");
296 token = strtok(NULL, " ");
297 token = strtok(NULL, " ");
298 token = strtok(NULL, " ");
299 int in0 = ReadInteger(token, -1, readerror);
300 token = strtok(NULL, " ");
301 int in1 = ReadInteger(token, -1, readerror);
302 token = strtok(NULL, " ");
303 int in2 = ReadInteger(token, -1, readerror);
304 token = strtok(NULL, " ");
305 int in3 = ReadInteger(token, -1, readerror);
306 token = strtok(NULL, " ");
307 int in4 = ReadInteger(token, -1, readerror);
308 token = strtok(NULL, " ");
309 int in5 = ReadInteger(token, -1, readerror);
310 token = strtok(NULL, " ");
311 int in6 = ReadInteger(token, -1, readerror);
312 token = strtok(NULL, " ");
313 int in7 = ReadInteger(token, -1, readerror);
314 if (!felist.getline(line, size, '\n')) {
315 std::cerr << m_className << "::Initialise:\n";
316 std::cerr << " Error reading element " << ielem << ".\n";
317 ok = false;
318 break;
319 }
320 token = NULL;
321 token = strtok(line, " ");
322 int in8 = ReadInteger(token, -1, readerror);
323 token = strtok(NULL, " ");
324 int in9 = ReadInteger(token, -1, readerror);
325
326 // Check synchronisation
327 if (readerror) {
328 std::cerr << m_className << "::Initialise:\n";
329 std::cerr << " Error reading file " << elist << " (line " << il
330 << ").\n";
331 felist.close();
332 ok = false;
333 return false;
334 } else if (ielem - 1 != nElements + nbackground) {
335 std::cerr << m_className << "::Initialise:\n";
336 std::cerr << " Synchronisation lost on file " << elist << " (line "
337 << il << ").\n";
338 std::cerr << " Element: " << ielem << " (expected " << nElements
339 << "), material: " << imat << ",\n";
340 std::cerr << " nodes: (" << in0 << ", " << in1 << ", " << in2 << ", "
341 << in3 << ", " << in4 << ", " << in5 << ", " << in6 << ", "
342 << in7 << ", " << in8 << ", " << in9 << ")\n";
343 ok = false;
344 }
345
346 // Check the material number and ensure that epsilon is non-negative
347 if (imat < 1 || imat > nMaterials) {
348 std::cerr << m_className << "::Initialise:\n";
349 std::cerr << " Out-of-range material number on file " << elist
350 << " (line " << il << ").\n";
351 std::cerr << " Element: " << ielem << ", material: " << imat << ",\n";
352 std::cerr << " nodes: (" << in0 << ", " << in1 << ", " << in2 << ", "
353 << in3 << ", " << in4 << ", " << in5 << ", " << in6 << ", "
354 << in7 << ", " << in8 << ", " << in9 << ")\n";
355 ok = false;
356 }
357 if (materials[imat - 1].eps < 0) {
358 std::cerr << m_className << "::Initialise:\n";
359 std::cerr << " Element " << ielem << " in element list " << elist
360 << "\n";
361 std::cerr << " uses material " << imat
362 << " which has not been assigned\n";
363 std::cerr << " a positive permittivity in material list " << mplist
364 << ".\n";
365 ok = false;
366 }
367
368 // Check the node numbers
369 if (in0 < 1 || in1 < 1 || in2 < 1 || in3 < 1 || in4 < 1 || in5 < 1 ||
370 in6 < 1 || in7 < 1 || in8 < 1 || in9 < 1) {
371 std::cerr << m_className << "::Initialise:\n";
372 std::cerr << " Found a node number < 1 on file " << elist << " (line "
373 << il << ").\n";
374 std::cerr << " Element: " << ielem << ", material: " << imat << ",\n";
375 std::cerr << " nodes: (" << in0 << ", " << in1 << ", " << in2 << ", "
376 << in3 << ", " << in4 << ", " << in5 << ", " << in6 << ", "
377 << in7 << ", " << in8 << ", " << in9 << ")\n";
378 ok = false;
379 }
380 if (in0 > highestnode) highestnode = in0;
381 if (in1 > highestnode) highestnode = in1;
382 if (in2 > highestnode) highestnode = in2;
383 if (in3 > highestnode) highestnode = in3;
384 if (in4 > highestnode) highestnode = in4;
385 if (in5 > highestnode) highestnode = in5;
386 if (in6 > highestnode) highestnode = in6;
387 if (in7 > highestnode) highestnode = in7;
388 if (in8 > highestnode) highestnode = in8;
389 if (in9 > highestnode) highestnode = in9;
390
391 // Skip quadrilaterals which are background.
392 if (deleteBackground && materials[imat - 1].ohm == 0) {
393 nbackground++;
394 continue;
395 }
396
397 // These elements must not be degenerate.
398 if (in0 == in1 || in0 == in2 || in0 == in3 || in0 == in4 || in0 == in5 ||
399 in0 == in6 || in0 == in7 || in0 == in8 || in0 == in9 || in1 == in2 ||
400 in1 == in3 || in1 == in4 || in1 == in5 || in1 == in6 || in1 == in7 ||
401 in1 == in8 || in1 == in9 || in2 == in3 || in2 == in4 || in2 == in5 ||
402 in2 == in6 || in2 == in7 || in2 == in8 || in2 == in9 || in3 == in4 ||
403 in3 == in5 || in3 == in6 || in3 == in7 || in3 == in8 || in3 == in9 ||
404 in4 == in5 || in4 == in6 || in4 == in7 || in4 == in8 || in4 == in9 ||
405 in5 == in6 || in5 == in7 || in5 == in8 || in5 == in9 || in6 == in7 ||
406 in6 == in8 || in6 == in9 || in7 == in8 || in7 == in9 || in8 == in9) {
407 std::cerr << m_className << "::Initialise:\n";
408 std::cerr << " Element " << ielem << " of file " << elist
409 << " is degenerate,\n";
410 std::cerr << " no such elements allowed in this type of map.\n";
411 ok = false;
412 }
413
414 newElement.degenerate = false;
415
416 // Store the material reference
417 newElement.matmap = imat - 1;
418
419 // Node references
420 newElement.emap[0] = in0 - 1;
421 newElement.emap[1] = in1 - 1;
422 newElement.emap[2] = in2 - 1;
423 newElement.emap[3] = in3 - 1;
424 newElement.emap[4] = in4 - 1;
425 newElement.emap[7] = in5 - 1;
426 newElement.emap[5] = in6 - 1;
427 newElement.emap[6] = in7 - 1;
428 newElement.emap[8] = in8 - 1;
429 newElement.emap[9] = in9 - 1;
430 elements.push_back(newElement);
431 nElements++;
432 }
433 // Close the file
434 felist.close();
435
436 // Tell how many lines read.
437 std::cout << m_className << "::Initialise:\n";
438 std::cout << " Read " << nElements << " elements from file " << elist
439 << ",\n";
440 std::cout << " highest node number: " << highestnode << ",\n";
441 std::cout << " background elements skipped: " << nbackground << "\n";
442 // Check the value of the unit
443 double funit;
444 if (strcmp(unit.c_str(), "mum") == 0 || strcmp(unit.c_str(), "micron") == 0 ||
445 strcmp(unit.c_str(), "micrometer") == 0) {
446 funit = 0.0001;
447 } else if (strcmp(unit.c_str(), "mm") == 0 ||
448 strcmp(unit.c_str(), "millimeter") == 0) {
449 funit = 0.1;
450 } else if (strcmp(unit.c_str(), "cm") == 0 ||
451 strcmp(unit.c_str(), "centimeter") == 0) {
452 funit = 1.0;
453 } else if (strcmp(unit.c_str(), "m") == 0 ||
454 strcmp(unit.c_str(), "meter") == 0) {
455 funit = 100.0;
456 } else {
457 std::cerr << m_className << "::Initialise:\n";
458 std::cerr << " Unknown length unit " << unit << ".\n";
459 ok = false;
460 funit = 1.0;
461 }
462 if (debug) {
463 std::cout << m_className << ":Initialise:\n";
464 std::cout << " Unit scaling factor = " << funit << ".\n";
465 }
466
467 // Open the node list
468 std::ifstream fnlist;
469 fnlist.open(nlist.c_str(), std::ios::in);
470 if (fnlist.fail()) {
471 std::cerr << m_className << "::Initialise:\n";
472 std::cerr << " Could not open nodes file " << nlist << " for reading.\n";
473 std::cerr << " The file perhaps does not exist.\n";
474 return false;
475 }
476
477 // Read the node list
478 nodes.clear();
479 nNodes = 0;
480 node newNode;
481 newNode.w.clear();
482 il = 0;
483 while (fnlist.getline(line, size, '\n')) {
484 il++;
485 // Skip page feed
486 if (strcmp(line, "1") == 0) {
487 fnlist.getline(line, size, '\n');
488 il++;
489 fnlist.getline(line, size, '\n');
490 il++;
491 fnlist.getline(line, size, '\n');
492 il++;
493 fnlist.getline(line, size, '\n');
494 il++;
495 fnlist.getline(line, size, '\n');
496 il++;
497 continue;
498 }
499 // Split the line in tokens
500 char* token = NULL;
501 token = strtok(line, " ");
502 // Skip blank lines and headers
503 if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
504 int(token[0]) == 10 || int(token[0]) == 13 ||
505 strcmp(token, "LIST") == 0 || strcmp(token, "NODE") == 0)
506 continue;
507 // Read the element
508 int inode = ReadInteger(token, -1, readerror);
509 token = strtok(NULL, " ");
510 double xnode = ReadDouble(token, -1, readerror);
511 token = strtok(NULL, " ");
512 double ynode = ReadDouble(token, -1, readerror);
513 token = strtok(NULL, " ");
514 double znode = ReadDouble(token, -1, readerror);
515 // Check syntax
516 if (readerror) {
517 std::cerr << m_className << "::Initialise:\n";
518 std::cerr << " Error reading file " << nlist << " (line " << il
519 << ").\n";
520 fnlist.close();
521 ok = false;
522 return false;
523 }
524 // Check synchronisation
525 if (inode - 1 != nNodes) {
526 std::cerr << m_className << "::Initialise:\n";
527 std::cerr << " Synchronisation lost on file " << nlist << " (line "
528 << il << ").\n";
529 std::cerr << " Node: " << inode << " (expected " << nNodes
530 << "), (x,y,z) = (" << xnode << ", " << ynode << ", " << znode
531 << ")\n";
532 ok = false;
533 }
534 // Store the point coordinates
535 newNode.x = xnode * funit;
536 newNode.y = ynode * funit;
537 newNode.z = znode * funit;
538 nodes.push_back(newNode);
539 nNodes++;
540 }
541 // Close the file
542 fnlist.close();
543 // Tell how many lines read
544 std::cout << m_className << "::Initialise:\n";
545 std::cout << " Read " << nNodes << " nodes from file " << nlist << ".\n";
546 // Check number of nodes
547 if (nNodes != highestnode) {
548 std::cerr << m_className << "::Initialise:\n";
549 std::cerr << " Number of nodes read (" << nNodes << ") on " << nlist
550 << "\n";
551 std::cerr << " does not match element list (" << highestnode << ").\n";
552 ok = false;
553 }
554
555 // Open the voltage list
556 std::ifstream fprnsol;
557 fprnsol.open(prnsol.c_str(), std::ios::in);
558 if (fprnsol.fail()) {
559 std::cerr << m_className << "::Initialise:\n";
560 std::cerr << " Could not open potential file " << prnsol
561 << " for reading.\n";
562 std::cerr << " The file perhaps does not exist.\n";
563 return false;
564 }
565
566 // Read the voltage list
567 il = 0;
568 int nread = 0;
569 while (fprnsol.getline(line, size, '\n')) {
570 il++;
571 // Skip page feed
572 if (strcmp(line, "1") == 0) {
573 fprnsol.getline(line, size, '\n');
574 il++;
575 fprnsol.getline(line, size, '\n');
576 il++;
577 fprnsol.getline(line, size, '\n');
578 il++;
579 fprnsol.getline(line, size, '\n');
580 il++;
581 fprnsol.getline(line, size, '\n');
582 il++;
583 continue;
584 }
585 // Split the line in tokens
586 char* token = NULL;
587 token = strtok(line, " ");
588 // Skip blank lines and headers
589 if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
590 int(token[0]) == 10 || int(token[0]) == 13 ||
591 strcmp(token, "PRINT") == 0 || strcmp(token, "*****") == 0 ||
592 strcmp(token, "LOAD") == 0 || strcmp(token, "TIME=") == 0 ||
593 strcmp(token, "MAXIMUM") == 0 || strcmp(token, "VALUE") == 0 ||
594 strcmp(token, "NODE") == 0)
595 continue;
596 // Read the element
597 int inode = ReadInteger(token, -1, readerror);
598 token = strtok(NULL, " ");
599 double volt = ReadDouble(token, -1, readerror);
600 // Check syntax
601 if (readerror) {
602 std::cerr << m_className << "::Initialise:\n";
603 std::cerr << " Error reading file " << prnsol << " (line << " << il
604 << ").\n";
605 fprnsol.close();
606 ok = false;
607 return false;
608 }
609 // Check node number and store if OK
610 if (inode < 1 || inode > highestnode) {
611 std::cerr << m_className << "::Initialise:\n";
612 std::cerr << " Node number " << inode << " out of range\n";
613 std::cerr << " on potential file " << prnsol << " (line " << il
614 << ").\n";
615 ok = false;
616 } else {
617 nodes[inode - 1].v = volt;
618 nread++;
619 }
620 }
621 // Close the file
622 fprnsol.close();
623 // Tell how many lines read
624 std::cout << m_className << "::Initialise:\n";
625 std::cout << " Read " << nread << " potentials from file " << prnsol
626 << ".\n";
627 // Check number of nodes
628 if (nread != nNodes) {
629 std::cerr << m_className << "::Initialise:\n";
630 std::cerr << " Number of nodes read (" << nread << ") on potential file "
631 << prnsol << "\n";
632 std::cerr << " does not match the node list (" << nNodes << ").\n";
633 ok = false;
634 }
635
636 // Set the ready flag
637 if (ok) {
638 ready = true;
639 } else {
640 std::cerr << m_className << "::Initialise:\n";
641 std::cerr
642 << " Field map could not be read and can not be interpolated.\n";
643 return false;
644 }
645
646 // Remove weighting fields (if any).
647 wfields.clear();
648 wfieldsOk.clear();
649 nWeightingFields = 0;
650
651 // Establish the ranges
652 SetRange();
653 UpdatePeriodicity();
654 return true;
655}
656
657bool ComponentAnsys123::SetWeightingField(std::string prnsol,
658 std::string label) {
659
660 if (!ready) {
661 std::cerr << m_className << "::SetWeightingField:\n";
662 std::cerr << " No valid field map is present.\n";
663 std::cerr << " Weighting field cannot be added.\n";
664 return false;
665 }
666
667 // Open the voltage list.
668 std::ifstream fprnsol;
669 fprnsol.open(prnsol.c_str(), std::ios::in);
670 if (fprnsol.fail()) {
671 std::cerr << m_className << "::SetWeightingField:\n";
672 std::cerr << " Could not open potential file " << prnsol
673 << " for reading.\n";
674 std::cerr << " The file perhaps does not exist.\n";
675 return false;
676 }
677
678 // Check if a weighting field with the same label already exists.
679 int iw = nWeightingFields;
680 for (int i = nWeightingFields; i--;) {
681 if (wfields[i] == label) {
682 iw = i;
683 break;
684 }
685 }
686 if (iw == nWeightingFields) {
687 ++nWeightingFields;
688 wfields.resize(nWeightingFields);
689 wfieldsOk.resize(nWeightingFields);
690 for (int j = nNodes; j--;) {
691 nodes[j].w.resize(nWeightingFields);
692 }
693 } else {
694 std::cout << m_className << "::SetWeightingField:\n";
695 std::cout << " Replacing existing weighting field " << label << ".\n";
696 }
697 wfields[iw] = label;
698 wfieldsOk[iw] = false;
699
700 // Buffer for reading
701 const int size = 100;
702 char line[size];
703
704 bool ok = true;
705 // Read the voltage list.
706 int il = 0;
707 int nread = 0;
708 bool readerror = false;
709
710 while (fprnsol.getline(line, size, '\n')) {
711 il++;
712 // Skip page feed
713 if (strcmp(line, "1") == 0) {
714 fprnsol.getline(line, size, '\n');
715 il++;
716 fprnsol.getline(line, size, '\n');
717 il++;
718 fprnsol.getline(line, size, '\n');
719 il++;
720 fprnsol.getline(line, size, '\n');
721 il++;
722 fprnsol.getline(line, size, '\n');
723 il++;
724 continue;
725 }
726 // Split the line in tokens.
727 char* token = NULL;
728 token = strtok(line, " ");
729 // Skip blank lines and headers.
730 if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
731 int(token[0]) == 10 || int(token[0]) == 13 ||
732 strcmp(token, "PRINT") == 0 || strcmp(token, "*****") == 0 ||
733 strcmp(token, "LOAD") == 0 || strcmp(token, "TIME=") == 0 ||
734 strcmp(token, "MAXIMUM") == 0 || strcmp(token, "VALUE") == 0 ||
735 strcmp(token, "NODE") == 0)
736 continue;
737 // Read the element.
738 int inode = ReadInteger(token, -1, readerror);
739 token = strtok(NULL, " ");
740 double volt = ReadDouble(token, -1, readerror);
741 // Check the syntax.
742 if (readerror) {
743 std::cerr << m_className << "::SetWeightingField:\n";
744 std::cerr << " Error reading file " << prnsol.c_str() << " (line "
745 << il << ").\n";
746 fprnsol.close();
747 return false;
748 }
749 // Check node number and store if OK.
750 if (inode < 1 || inode > nNodes) {
751 std::cerr << m_className << "::SetWeightingField:\n";
752 std::cerr << " Node number " << inode << " out of range\n";
753 std::cerr << " on potential file " << prnsol.c_str() << " (line " << il
754 << ").\n";
755 ok = false;
756 } else {
757 nodes[inode - 1].w[iw] = volt;
758 nread++;
759 }
760 }
761 // Close the file.
762 fprnsol.close();
763
764 std::cout << m_className << "::SetWeightingField:\n";
765 std::cout << " Read " << nread << " potentials from file "
766 << prnsol.c_str() << ".\n";
767 // Check the number of nodes.
768 if (nread != nNodes) {
769 std::cerr << m_className << "::SetWeightingField:\n";
770 std::cerr << " Number of nodes read (" << nread << ") "
771 << " on potential file " << prnsol.c_str() << "\n";
772 std::cerr << " does not match the node list (" << nNodes << ").\n";
773 ok = false;
774 }
775
776 // Set the ready flag.
777 wfieldsOk[iw] = ok;
778 if (!ok) {
779 std::cerr << m_className << "::SetWeightingField:\n";
780 std::cerr << " Field map could not be read "
781 << "and cannot be interpolated.\n";
782 return false;
783 }
784 return true;
785}
786
787void ComponentAnsys123::ElectricField(const double x, const double y,
788 const double z, double& ex, double& ey,
789 double& ez, Medium*& m, int& status) {
790
791 double v = 0.;
792 ElectricField(x, y, z, ex, ey, ez, v, m, status);
793}
794
795void ComponentAnsys123::ElectricField(const double xin, const double yin,
796 const double zin, double& ex, double& ey,
797 double& ez, double& volt, Medium*& m,
798 int& status) {
799
800 // Copy the coordinates
801 double x = xin, y = yin, z = zin;
802
803 // Map the coordinates onto field map coordinates
804 bool xmirrored, ymirrored, zmirrored;
805 double rcoordinate, rotation;
806 MapCoordinates(x, y, z, xmirrored, ymirrored, zmirrored, rcoordinate,
807 rotation);
808
809 // Initial values
810 ex = ey = ez = volt = 0.;
811 status = 0;
812 m = 0;
813
814 // Do not proceed if not properly initialised.
815 if (!ready) {
816 status = -10;
817 std::cerr << m_className << "::ElectricField:\n";
818 std::cerr << " Field map not available for interpolation.\n";
819 return;
820 }
821
822 if (warning) {
823 std::cerr << m_className << "::ElectricField:\n";
824 std::cerr << " Warnings have been issued for this field map.\n";
825 }
826
827 // Find the element that contains this point
828 double t1, t2, t3, t4, jac[4][4], det;
829 int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
830 if (imap < 0) {
831 if (debug) {
832 std::cerr << m_className << "::ElectricField:\n";
833 std::cerr << " Point (" << x << ", " << y << ", " << z
834 << ") not in the mesh.\n";
835 }
836 status = -6;
837 return;
838 }
839
840 if (debug) {
841 std::cout << m_className << "::ElectricField:\n";
842 std::cout << " Global: (" << x << ", " << y << ", " << z << "),\n";
843 std::cout << " Local: (" << t1 << ", " << t2 << ", " << t3 << ", " << t4
844 << ") in element " << imap << ".\n";
845 std::cout
846 << " Node x y z V\n";
847 for (int i = 0; i < 10; i++) {
848 printf(" %-5d %12g %12g %12g %12g\n", elements[imap].emap[i],
849 nodes[elements[imap].emap[i]].x, nodes[elements[imap].emap[i]].y,
850 nodes[elements[imap].emap[i]].z, nodes[elements[imap].emap[i]].v);
851 }
852 }
853
854 // Tetrahedral field
855 volt = nodes[elements[imap].emap[0]].v * t1 * (2 * t1 - 1) +
856 nodes[elements[imap].emap[1]].v * t2 * (2 * t2 - 1) +
857 nodes[elements[imap].emap[2]].v * t3 * (2 * t3 - 1) +
858 nodes[elements[imap].emap[3]].v * t4 * (2 * t4 - 1) +
859 4 * nodes[elements[imap].emap[4]].v * t1 * t2 +
860 4 * nodes[elements[imap].emap[5]].v * t1 * t3 +
861 4 * nodes[elements[imap].emap[6]].v * t1 * t4 +
862 4 * nodes[elements[imap].emap[7]].v * t2 * t3 +
863 4 * nodes[elements[imap].emap[8]].v * t2 * t4 +
864 4 * nodes[elements[imap].emap[9]].v * t3 * t4;
865
866 ex = -(nodes[elements[imap].emap[0]].v * (4 * t1 - 1) * jac[0][1] +
867 nodes[elements[imap].emap[1]].v * (4 * t2 - 1) * jac[1][1] +
868 nodes[elements[imap].emap[2]].v * (4 * t3 - 1) * jac[2][1] +
869 nodes[elements[imap].emap[3]].v * (4 * t4 - 1) * jac[3][1] +
870 nodes[elements[imap].emap[4]].v *
871 (4 * t2 * jac[0][1] + 4 * t1 * jac[1][1]) +
872 nodes[elements[imap].emap[5]].v *
873 (4 * t3 * jac[0][1] + 4 * t1 * jac[2][1]) +
874 nodes[elements[imap].emap[6]].v *
875 (4 * t4 * jac[0][1] + 4 * t1 * jac[3][1]) +
876 nodes[elements[imap].emap[7]].v *
877 (4 * t3 * jac[1][1] + 4 * t2 * jac[2][1]) +
878 nodes[elements[imap].emap[8]].v *
879 (4 * t4 * jac[1][1] + 4 * t2 * jac[3][1]) +
880 nodes[elements[imap].emap[9]].v *
881 (4 * t4 * jac[2][1] + 4 * t3 * jac[3][1])) /
882 det;
883
884 ey = -(nodes[elements[imap].emap[0]].v * (4 * t1 - 1) * jac[0][2] +
885 nodes[elements[imap].emap[1]].v * (4 * t2 - 1) * jac[1][2] +
886 nodes[elements[imap].emap[2]].v * (4 * t3 - 1) * jac[2][2] +
887 nodes[elements[imap].emap[3]].v * (4 * t4 - 1) * jac[3][2] +
888 nodes[elements[imap].emap[4]].v *
889 (4 * t2 * jac[0][2] + 4 * t1 * jac[1][2]) +
890 nodes[elements[imap].emap[5]].v *
891 (4 * t3 * jac[0][2] + 4 * t1 * jac[2][2]) +
892 nodes[elements[imap].emap[6]].v *
893 (4 * t4 * jac[0][2] + 4 * t1 * jac[3][2]) +
894 nodes[elements[imap].emap[7]].v *
895 (4 * t3 * jac[1][2] + 4 * t2 * jac[2][2]) +
896 nodes[elements[imap].emap[8]].v *
897 (4 * t4 * jac[1][2] + 4 * t2 * jac[3][2]) +
898 nodes[elements[imap].emap[9]].v *
899 (4 * t4 * jac[2][2] + 4 * t3 * jac[3][2])) /
900 det;
901
902 ez = -(nodes[elements[imap].emap[0]].v * (4 * t1 - 1) * jac[0][3] +
903 nodes[elements[imap].emap[1]].v * (4 * t2 - 1) * jac[1][3] +
904 nodes[elements[imap].emap[2]].v * (4 * t3 - 1) * jac[2][3] +
905 nodes[elements[imap].emap[3]].v * (4 * t4 - 1) * jac[3][3] +
906 nodes[elements[imap].emap[4]].v *
907 (4 * t2 * jac[0][3] + 4 * t1 * jac[1][3]) +
908 nodes[elements[imap].emap[5]].v *
909 (4 * t3 * jac[0][3] + 4 * t1 * jac[2][3]) +
910 nodes[elements[imap].emap[6]].v *
911 (4 * t4 * jac[0][3] + 4 * t1 * jac[3][3]) +
912 nodes[elements[imap].emap[7]].v *
913 (4 * t3 * jac[1][3] + 4 * t2 * jac[2][3]) +
914 nodes[elements[imap].emap[8]].v *
915 (4 * t4 * jac[1][3] + 4 * t2 * jac[3][3]) +
916 nodes[elements[imap].emap[9]].v *
917 (4 * t4 * jac[2][3] + 4 * t3 * jac[3][3])) /
918 det;
919
920 // Transform field to global coordinates
921 UnmapFields(ex, ey, ez, x, y, z, xmirrored, ymirrored, zmirrored, rcoordinate,
922 rotation);
923
924 // Drift medium?
925 if (debug) {
926 std::cout << m_className << "::ElectricField:\n";
927 std::cout << " Material " << elements[imap].matmap << ", drift flag "
928 << materials[elements[imap].matmap].driftmedium << ".\n";
929 }
930 m = materials[elements[imap].matmap].medium;
931 status = -5;
932 if (materials[elements[imap].matmap].driftmedium) {
933 if (m != 0) {
934 if (m->IsDriftable()) status = 0;
935 }
936 }
937}
938
939void ComponentAnsys123::WeightingField(const double xin, const double yin,
940 const double zin, double& wx, double& wy,
941 double& wz, const std::string label) {
942
943 // Initial values
944 wx = wy = wz = 0;
945
946 // Do not proceed if not properly initialised.
947 if (!ready) return;
948
949 // Look for the label.
950 int iw = 0;
951 bool found = false;
952 for (int i = nWeightingFields; i--;) {
953 if (wfields[i] == label) {
954 iw = i;
955 found = true;
956 break;
957 }
958 }
959
960 // Do not proceed if the requested weighting field does not exist.
961 if (!found) return;
962 // Check if the weighting field is properly initialised.
963 if (!wfieldsOk[iw]) return;
964
965 // Copy the coordinates.
966 double x = xin, y = yin, z = zin;
967
968 // Map the coordinates onto field map coordinates
969 bool xmirrored, ymirrored, zmirrored;
970 double rcoordinate, rotation;
971 MapCoordinates(x, y, z, xmirrored, ymirrored, zmirrored, rcoordinate,
972 rotation);
973
974 if (warning) {
975 std::cerr << m_className << "::WeightingField:\n";
976 std::cerr << " Warnings have been issued for this field map.\n";
977 }
978
979 // Find the element that contains this point.
980 double t1, t2, t3, t4, jac[4][4], det;
981 int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
982 // Check if the point is in the mesh.
983 if (imap < 0) return;
984
985 if (debug) {
986 std::cout << m_className << "::WeightingField:\n";
987 std::cout << " Global: (" << x << ", " << y << ", " << z << ")\n";
988 std::cout << " Local: (" << t1 << ", " << t2 << ", " << t3 << ", " << t4
989 << ") in element " << imap << "\n",
990 std::cout << " Node x y z "
991 " V\n";
992 for (int i = 0; i < 10; i++) {
993 printf(" %-5d %12g %12g %12g %12g\n", elements[imap].emap[i],
994 nodes[elements[imap].emap[i]].x, nodes[elements[imap].emap[i]].y,
995 nodes[elements[imap].emap[i]].z,
996 nodes[elements[imap].emap[i]].w[iw]);
997 }
998 }
999
1000 // Tetrahedral field
1001 wx = -(nodes[elements[imap].emap[0]].w[iw] * (4 * t1 - 1) * jac[0][1] +
1002 nodes[elements[imap].emap[1]].w[iw] * (4 * t2 - 1) * jac[1][1] +
1003 nodes[elements[imap].emap[2]].w[iw] * (4 * t3 - 1) * jac[2][1] +
1004 nodes[elements[imap].emap[3]].w[iw] * (4 * t4 - 1) * jac[3][1] +
1005 nodes[elements[imap].emap[4]].w[iw] *
1006 (4 * t2 * jac[0][1] + 4 * t1 * jac[1][1]) +
1007 nodes[elements[imap].emap[5]].w[iw] *
1008 (4 * t3 * jac[0][1] + 4 * t1 * jac[2][1]) +
1009 nodes[elements[imap].emap[6]].w[iw] *
1010 (4 * t4 * jac[0][1] + 4 * t1 * jac[3][1]) +
1011 nodes[elements[imap].emap[7]].w[iw] *
1012 (4 * t3 * jac[1][1] + 4 * t2 * jac[2][1]) +
1013 nodes[elements[imap].emap[8]].w[iw] *
1014 (4 * t4 * jac[1][1] + 4 * t2 * jac[3][1]) +
1015 nodes[elements[imap].emap[9]].w[iw] *
1016 (4 * t4 * jac[2][1] + 4 * t3 * jac[3][1])) /
1017 det;
1018
1019 wy = -(nodes[elements[imap].emap[0]].w[iw] * (4 * t1 - 1) * jac[0][2] +
1020 nodes[elements[imap].emap[1]].w[iw] * (4 * t2 - 1) * jac[1][2] +
1021 nodes[elements[imap].emap[2]].w[iw] * (4 * t3 - 1) * jac[2][2] +
1022 nodes[elements[imap].emap[3]].w[iw] * (4 * t4 - 1) * jac[3][2] +
1023 nodes[elements[imap].emap[4]].w[iw] *
1024 (4 * t2 * jac[0][2] + 4 * t1 * jac[1][2]) +
1025 nodes[elements[imap].emap[5]].w[iw] *
1026 (4 * t3 * jac[0][2] + 4 * t1 * jac[2][2]) +
1027 nodes[elements[imap].emap[6]].w[iw] *
1028 (4 * t4 * jac[0][2] + 4 * t1 * jac[3][2]) +
1029 nodes[elements[imap].emap[7]].w[iw] *
1030 (4 * t3 * jac[1][2] + 4 * t2 * jac[2][2]) +
1031 nodes[elements[imap].emap[8]].w[iw] *
1032 (4 * t4 * jac[1][2] + 4 * t2 * jac[3][2]) +
1033 nodes[elements[imap].emap[9]].w[iw] *
1034 (4 * t4 * jac[2][2] + 4 * t3 * jac[3][2])) /
1035 det;
1036
1037 wz = -(nodes[elements[imap].emap[0]].w[iw] * (4 * t1 - 1) * jac[0][3] +
1038 nodes[elements[imap].emap[1]].w[iw] * (4 * t2 - 1) * jac[1][3] +
1039 nodes[elements[imap].emap[2]].w[iw] * (4 * t3 - 1) * jac[2][3] +
1040 nodes[elements[imap].emap[3]].w[iw] * (4 * t4 - 1) * jac[3][3] +
1041 nodes[elements[imap].emap[4]].w[iw] *
1042 (4 * t2 * jac[0][3] + 4 * t1 * jac[1][3]) +
1043 nodes[elements[imap].emap[5]].w[iw] *
1044 (4 * t3 * jac[0][3] + 4 * t1 * jac[2][3]) +
1045 nodes[elements[imap].emap[6]].w[iw] *
1046 (4 * t4 * jac[0][3] + 4 * t1 * jac[3][3]) +
1047 nodes[elements[imap].emap[7]].w[iw] *
1048 (4 * t3 * jac[1][3] + 4 * t2 * jac[2][3]) +
1049 nodes[elements[imap].emap[8]].w[iw] *
1050 (4 * t4 * jac[1][3] + 4 * t2 * jac[3][3]) +
1051 nodes[elements[imap].emap[9]].w[iw] *
1052 (4 * t4 * jac[2][3] + 4 * t3 * jac[3][3])) /
1053 det;
1054
1055 // Transform field to global coordinates
1056 UnmapFields(wx, wy, wz, x, y, z, xmirrored, ymirrored, zmirrored, rcoordinate,
1057 rotation);
1058}
1059
1060double ComponentAnsys123::WeightingPotential(const double xin, const double yin,
1061 const double zin,
1062 const std::string label) {
1063
1064 // Do not proceed if not properly initialised.
1065 if (!ready) return 0.;
1066
1067 // Look for the label.
1068 int iw = 0;
1069 bool found = false;
1070 for (int i = nWeightingFields; i--;) {
1071 if (wfields[i] == label) {
1072 iw = i;
1073 found = true;
1074 break;
1075 }
1076 }
1077
1078 // Do not proceed if the requested weighting field does not exist.
1079 if (!found) return 0.;
1080 // Check if the weighting field is properly initialised.
1081 if (!wfieldsOk[iw]) return 0.;
1082
1083 // Copy the coordinates.
1084 double x = xin, y = yin, z = zin;
1085
1086 // Map the coordinates onto field map coordinates.
1087 bool xmirrored, ymirrored, zmirrored;
1088 double rcoordinate, rotation;
1089 MapCoordinates(x, y, z, xmirrored, ymirrored, zmirrored, rcoordinate,
1090 rotation);
1091
1092 if (warning) {
1093 std::cerr << m_className << "::WeightingPotential:\n";
1094 std::cerr << " Warnings have been issued for this field map.\n";
1095 }
1096
1097 // Find the element that contains this point.
1098 double t1, t2, t3, t4, jac[4][4], det;
1099 int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
1100 if (imap < 0) return 0.;
1101
1102 if (debug) {
1103 std::cerr << m_className << "::WeightingPotential:\n";
1104 std::cout << " Global: (" << x << ", " << y << ", "
1105 << ", " << z << ")\n";
1106 std::cout << " Local: (" << t1 << ", " << t2 << ", " << t3 << ", " << t4
1107 << ") in element " << imap << "\n";
1108 std::cout
1109 << " Node x y z V\n";
1110 for (int i = 0; i < 10; i++) {
1111 printf(" %-5d %12g %12g %12g %12g\n", elements[imap].emap[i],
1112 nodes[elements[imap].emap[i]].x, nodes[elements[imap].emap[i]].y,
1113 nodes[elements[imap].emap[i]].z, nodes[elements[imap].emap[i]].v);
1114 }
1115 }
1116
1117 // Tetrahedral field
1118 return nodes[elements[imap].emap[0]].w[iw] * t1 * (2 * t1 - 1) +
1119 nodes[elements[imap].emap[1]].w[iw] * t2 * (2 * t2 - 1) +
1120 nodes[elements[imap].emap[2]].w[iw] * t3 * (2 * t3 - 1) +
1121 nodes[elements[imap].emap[3]].w[iw] * t4 * (2 * t4 - 1) +
1122 4 * nodes[elements[imap].emap[4]].w[iw] * t1 * t2 +
1123 4 * nodes[elements[imap].emap[5]].w[iw] * t1 * t3 +
1124 4 * nodes[elements[imap].emap[6]].w[iw] * t1 * t4 +
1125 4 * nodes[elements[imap].emap[7]].w[iw] * t2 * t3 +
1126 4 * nodes[elements[imap].emap[8]].w[iw] * t2 * t4 +
1127 4 * nodes[elements[imap].emap[9]].w[iw] * t3 * t4;
1128}
1129
1130Medium* ComponentAnsys123::GetMedium(const double& xin, const double& yin,
1131 const double& zin) {
1132
1133 // Copy the coordinates.
1134 double x = xin, y = yin, z = zin;
1135
1136 // Map the coordinates onto field map coordinates.
1137 bool xmirrored, ymirrored, zmirrored;
1138 double rcoordinate, rotation;
1139 MapCoordinates(x, y, z, xmirrored, ymirrored, zmirrored, rcoordinate,
1140 rotation);
1141
1142 // Do not proceed if not properly initialised.
1143 if (!ready) {
1144 std::cerr << m_className << "::GetMedium:\n";
1145 std::cerr << " Field map not available for interpolation.\n";
1146 return NULL;
1147 }
1148 if (warning) {
1149 std::cerr << m_className << "::GetMedium:\n";
1150 std::cerr << " Warnings have been issued for this field map.\n";
1151 }
1152
1153 // Find the element that contains this point.
1154 double t1, t2, t3, t4, jac[4][4], det;
1155 int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
1156 if (imap < 0) {
1157 if (debug) {
1158 std::cerr << m_className << "::GetMedium:\n";
1159 std::cerr << " Point (" << x << ", " << y << ", " << z
1160 << ") not in the mesh.\n";
1161 }
1162 return NULL;
1163 }
1164 if (elements[imap].matmap < 0 || elements[imap].matmap >= nMaterials) {
1165 if (debug) {
1166 std::cerr << m_className << "::GetMedium:\n";
1167 std::cerr << " Point (" << x << ", " << y << ", " << z << ")"
1168 << " has out of range material number " << imap << ".\n";
1169 }
1170 return NULL;
1171 }
1172
1173 if (debug) {
1174 std::cout << m_className << "::GetMedium:\n";
1175 std::cout << " Global: (" << x << ", " << y << ", " << z << ")\n";
1176 std::cout << " Local: (" << t1 << ", " << t2 << ", " << t3 << ", " << t4
1177 << ") in element " << imap << "\n";
1178 std::cout
1179 << " Node x y z V\n";
1180 for (int ii = 0; ii < 10; ++ii) {
1181 printf(" %-5d %12g %12g %12g %12g\n", elements[imap].emap[ii],
1182 nodes[elements[imap].emap[ii]].x, nodes[elements[imap].emap[ii]].y,
1183 nodes[elements[imap].emap[ii]].z,
1184 nodes[elements[imap].emap[ii]].v);
1185 }
1186 }
1187
1188 // Assign a medium.
1189 return materials[elements[imap].matmap].medium;
1190}
1191
1192double ComponentAnsys123::GetElementVolume(const int i) {
1193
1194 if (i < 0 || i >= nElements) return 0.;
1195
1196 const double vol =
1197 fabs((nodes[elements[i].emap[3]].x - nodes[elements[i].emap[0]].x) *
1198 ((nodes[elements[i].emap[1]].y - nodes[elements[i].emap[0]].y) *
1199 (nodes[elements[i].emap[2]].z -
1200 nodes[elements[i].emap[0]].z) -
1201 (nodes[elements[i].emap[2]].y - nodes[elements[i].emap[0]].y) *
1202 (nodes[elements[i].emap[1]].z -
1203 nodes[elements[i].emap[0]].z)) +
1204 (nodes[elements[i].emap[3]].y - nodes[elements[i].emap[0]].y) *
1205 ((nodes[elements[i].emap[1]].z - nodes[elements[i].emap[0]].z) *
1206 (nodes[elements[i].emap[2]].x -
1207 nodes[elements[i].emap[0]].x) -
1208 (nodes[elements[i].emap[2]].z - nodes[elements[i].emap[0]].z) *
1209 (nodes[elements[i].emap[1]].x -
1210 nodes[elements[i].emap[0]].x)) +
1211 (nodes[elements[i].emap[3]].z - nodes[elements[i].emap[0]].z) *
1212 ((nodes[elements[i].emap[1]].x - nodes[elements[i].emap[0]].x) *
1213 (nodes[elements[i].emap[2]].y -
1214 nodes[elements[i].emap[0]].y) -
1215 (nodes[elements[i].emap[3]].x - nodes[elements[i].emap[0]].x) *
1216 (nodes[elements[i].emap[1]].y -
1217 nodes[elements[i].emap[0]].y))) /
1218 6.;
1219 return vol;
1220}
1221
1222void ComponentAnsys123::GetAspectRatio(const int i, double& dmin,
1223 double& dmax) {
1224
1225 if (i < 0 || i >= nElements) {
1226 dmin = dmax = 0.;
1227 return;
1228 }
1229
1230 const int np = 4;
1231 // Loop over all pairs of vertices.
1232 for (int j = 0; j < np - 1; ++j) {
1233 for (int k = j + 1; k < np; ++k) {
1234 // Compute distance.
1235 const double dist = sqrt(
1236 pow(nodes[elements[i].emap[j]].x - nodes[elements[i].emap[k]].x, 2) +
1237 pow(nodes[elements[i].emap[j]].y - nodes[elements[i].emap[k]].y, 2) +
1238 pow(nodes[elements[i].emap[j]].z - nodes[elements[i].emap[k]].z, 2));
1239 if (k == 1) {
1240 dmin = dist;
1241 dmax = dist;
1242 } else {
1243 if (dist < dmin) dmin = dist;
1244 if (dist > dmax) dmax = dist;
1245 }
1246 }
1247 }
1248}
1249}
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::ComponentAnsys123::Initialise
bool Initialise(std::string elist="ELIST.lis", std::string nlist="NLIST.lis", std::string mplist="MPLIST.lis", std::string prnsol="PRNSOL.lis", std::string unit="cm")
Definition: ComponentAnsys123.cc:18
Garfield::ComponentAnsys123::GetElementVolume
double GetElementVolume(const int i)
Definition: ComponentAnsys123.cc:1192
Garfield::ComponentAnsys123::GetAspectRatio
void GetAspectRatio(const int i, double &dmin, double &dmax)
Definition: ComponentAnsys123.cc:1222
Garfield::ComponentAnsys123::GetMedium
Medium * GetMedium(const double &x, const double &y, const double &z)
Definition: ComponentAnsys123.cc:1130
Garfield::ComponentAnsys123::WeightingPotential
double WeightingPotential(const double x, const double y, const double z, const std::string label)
Definition: ComponentAnsys123.cc:1060
Garfield::ComponentAnsys123::ElectricField
void ElectricField(const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status)
Definition: ComponentAnsys123.cc:787
Garfield::ComponentAnsys123::SetWeightingField
bool SetWeightingField(std::string prnsol, std::string label)
Definition: ComponentAnsys123.cc:657
Garfield::ComponentAnsys123::WeightingField
void WeightingField(const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string label)
Definition: ComponentAnsys123.cc:939
Garfield::ComponentFieldMap::ReadDouble
double ReadDouble(char *token, double def, bool &error)
Definition: ComponentFieldMap.cc:4240
Garfield::ComponentFieldMap::MapCoordinates
void MapCoordinates(double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
Definition: ComponentFieldMap.cc:4036
Garfield::ComponentFieldMap::materials
std::vector< material > materials
Definition: ComponentFieldMap.hh:124
Garfield::ComponentFieldMap::UnmapFields
void UnmapFields(double &ex, double &ey, double &ez, double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation)
Definition: ComponentFieldMap.cc:4154
Garfield::ComponentFieldMap::ReadInteger
int ReadInteger(char *token, int def, bool &error)
Definition: ComponentFieldMap.cc:4230
Garfield::ComponentFieldMap::elements
std::vector< element > elements
Definition: ComponentFieldMap.hh:96
Garfield::ComponentFieldMap::FindElement13
int FindElement13(const double x, const double y, const double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det)
Definition: ComponentFieldMap.cc:356
Garfield::ComponentFieldMap::wfields
std::vector< std::string > wfields
Definition: ComponentFieldMap.hh:127
Garfield::Medium::IsDriftable
bool IsDriftable() const
Definition: Medium.hh:57