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