Garfield++ v2r0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
ComponentFieldMap.cc
Go to the documentation of this file.
1#include <stdio.h>
2#include <iostream>
3#include <fstream>
4
5#include <math.h>
6#include <string>
7
8#include "ComponentFieldMap.hh"
9#include "FundamentalConstants.hh"
10
11namespace Garfield {
12
13ComponentFieldMap::ComponentFieldMap() : ComponentBase(),
14 m_is3d(true),
15 nElements(-1),
16 nNodes(-1),
17 m_nMaterials(0),
18 nWeightingFields(0),
19 hasBoundingBox(false),
20 m_deleteBackground(true),
21 m_warning(false), m_nWarnings(0),
22 m_checkMultipleElement(false),
23 m_useTetrahedralTree(false),
24 m_isTreeInitialized(false),
25 m_tetTree(NULL),
26 m_cacheElemBoundingBoxes(false),
27 m_lastElement(-1) {
28
29 m_className = "ComponentFieldMap";
30}
31
32ComponentFieldMap::~ComponentFieldMap() {
33 if (m_tetTree) delete m_tetTree;
34}
35
36void ComponentFieldMap::PrintMaterials() {
37
38 // Do not proceed if not properly initialised.
39 if (!m_ready) PrintNotReady("PrintMaterials");
40
41 if (materials.empty()) {
42 std::cerr << m_className << "::PrintMaterials:\n"
43 << " No materials are currently defined.\n";
44 return;
45 }
46
47 std::cout << m_className << "::PrintMaterials:\n"
48 << " Currently " << m_nMaterials << " materials are defined.\n"
49 << " Index Permittivity Resistivity Notes\n";
50 for (unsigned int i = 0; i < m_nMaterials; ++i) {
51 printf(" %5d %12g %12g", i, materials[i].eps, materials[i].ohm);
52 if (materials[i].medium) {
53 std::string name = materials[i].medium->GetName();
54 std::cout << " " << name;
55 if (materials[i].medium->IsDriftable()) std::cout << ", drift medium";
56 if (materials[i].medium->IsIonisable()) std::cout << ", ionisable";
57 }
58 if (materials[i].driftmedium) {
59 std::cout << " (drift medium)\n";
60 } else {
61 std::cout << "\n";
62 }
63 }
64}
65
66void ComponentFieldMap::DriftMedium(const unsigned int imat) {
67
68 // Do not proceed if not properly initialised.
69 if (!m_ready) PrintNotReady("DriftMedium");
70
71 // Check value
72 if (imat >= m_nMaterials) {
73 std::cerr << m_className << "::DriftMedium:\n";
74 std::cerr << " Material index " << imat << " is out of range.\n";
75 return;
76 }
77
78 // Make drift medium
79 materials[imat].driftmedium = true;
80}
81
82void ComponentFieldMap::NotDriftMedium(const unsigned int imat) {
83
84 // Do not proceed if not properly initialised.
85 if (!m_ready) PrintNotReady("NotDriftMedium");
86
87 // Check value
88 if (imat >= m_nMaterials) {
89 std::cerr << m_className << "::NotDriftMedium:\n";
90 std::cerr << " Material index " << imat << " is out of range.\n";
91 return;
92 }
93
94 // Make drift medium
95 materials[imat].driftmedium = false;
96}
97
98double ComponentFieldMap::GetPermittivity(const unsigned int imat) const {
99
100 if (imat >= m_nMaterials) {
101 std::cerr << m_className << "::GetPermittivity:\n"
102 << " Material index " << imat << " is out of range.\n";
103 return -1.;
104 }
105
106 return materials[imat].eps;
107}
108
109double ComponentFieldMap::GetConductivity(const unsigned int imat) const {
110
111 if (imat >= m_nMaterials) {
112 std::cerr << m_className << "::GetConductivity:\n"
113 << " Material index " << imat << " is out of range.\n";
114 return -1.;
115 }
116
117 return materials[imat].ohm;
118}
119
120void ComponentFieldMap::SetMedium(const unsigned int imat, Medium* m) {
121
122 if (imat >= m_nMaterials) {
123 std::cerr << m_className << "::SetMedium:\n";
124 std::cerr << " Material index " << imat << " is out of range.\n";
125 return;
126 }
127
128 if (!m) {
129 std::cerr << m_className << "::SetMedium: Null pointer.\n";
130 return;
131 }
132
133 if (m_debug) {
134 std::cout << m_className << "::SetMedium:\n Associated material "
135 << imat << " with medium " << m->GetName() << ".\n";
136 }
137
138 materials[imat].medium = m;
139}
140
141Medium* ComponentFieldMap::GetMedium(const unsigned int imat) const {
142
143 if (imat >= m_nMaterials) {
144 std::cerr << m_className << "::GetMedium:\n"
145 << " Material index " << imat << " is out of range.\n";
146 return NULL;
147 }
148
149 return materials[imat].medium;
150}
151
152bool ComponentFieldMap::GetElement(const unsigned int i, double& vol,
153 double& dmin, double& dmax) {
154
155 if ((int)i >= nElements) {
156 std::cerr << m_className << "::GetElement:\n";
157 std::cerr << " Element index (" << i << ") out of range.\n";
158 return false;
159 }
160
161 vol = GetElementVolume(i);
162 GetAspectRatio(i, dmin, dmax);
163 return true;
164}
165
166int ComponentFieldMap::FindElement5(const double x, const double y,
167 double const z, double& t1, double& t2,
168 double& t3, double& t4, double jac[4][4],
169 double& det) {
170
171 // Check if bounding boxes of elements have been computed
172 if (!m_cacheElemBoundingBoxes) {
173 std::cout << m_className << "::FindElement5:\n"
174 << " Caching the bounding boxes of all elements...";
175 CalculateElementBoundingBoxes();
176 std::cout << " done.\n";
177 m_cacheElemBoundingBoxes = true;
178 }
179
180 // Tetra list in the block that contains the input 3D point.
181 std::vector<int> tetList;
182 if (m_useTetrahedralTree) {
183 if (!m_isTreeInitialized) {
184 if (!InitializeTetrahedralTree()) {
185 std::cerr << m_className << "::FindElement5:\n";
186 std::cerr << " Tetrahedral tree initialization failed.\n";
187 return -1;
188 }
189 }
190 tetList = m_tetTree->GetTetListInBlock(Vec3(x, y, z));
191 }
192 // Backup
193 double jacbak[4][4], detbak = 1.;
194 double t1bak = 0., t2bak = 0., t3bak = 0., t4bak = 0.;
195 int imapbak = -1;
196
197 // Initial values.
198 t1 = t2 = t3 = t4 = 0;
199
200 // Check previously used element
201 if (m_lastElement > -1 && !m_checkMultipleElement) {
202 if (elements[m_lastElement].degenerate) {
203 if (Coordinates3(x, y, z, t1, t2, t3, t4, jac, det, m_lastElement) == 0) {
204 if (t1 >= 0 && t1 <= +1 && t2 >= 0 && t2 <= +1 && t3 >= 0 && t3 <= +1) {
205 return m_lastElement;
206 }
207 }
208 } else {
209 if (Coordinates5(x, y, z, t1, t2, t3, t4, jac, det, m_lastElement) == 0) {
210 if (t1 >= -1 && t1 <= +1 && t2 >= -1 && t2 <= +1) return m_lastElement;
211 }
212 }
213 }
214
215 // Verify the count of volumes that contain the point.
216 int nfound = 0;
217 int imap = -1;
218
219 // Number of elements to scan.
220 // With tetra tree disabled, all elements are scanned.
221 const int numElemToSearch =
222 m_useTetrahedralTree ? tetList.size() : nElements;
223 for (int i = 0; i < numElemToSearch; ++i) {
224 const int idxToElemList = m_useTetrahedralTree ? tetList[i] : i;
225 const Element& element = elements[idxToElemList];
226 // Tolerance
227 const double f = 0.2;
228 const double tolx = f * (element.xmax - element.xmin);
229 if (x < element.xmin - tolx || x > element.xmax + tolx) continue;
230 const double toly = f * (element.ymax - element.ymin);
231 if (y < element.ymin - toly || y > element.ymax + toly) continue;
232 const double tolz = f * (element.zmax - element.zmin);
233 if (z < element.zmin - tolz || z > element.zmax + tolz) continue;
234
235 if (element.degenerate) {
236 // Degenerate element
237 if (Coordinates3(x, y, z, t1, t2, t3, t4, jac, det, idxToElemList) != 0) {
238 continue;
239 }
240 if (t1 < 0 || t1 > 1 || t2 < 0 || t2 > 1 || t3 < 0 || t3 > 1) continue;
241 ++nfound;
242 imap = idxToElemList;
243 m_lastElement = idxToElemList;
244 if (m_debug) {
245 std::cout << m_className << "::FindElement5:\n";
246 std::cout << " Found matching degenerate element " << idxToElemList
247 << ".\n";
248 }
249 if (!m_checkMultipleElement) return idxToElemList;
250 for (int j = 0; j < 4; ++j) {
251 for (int k = 0; k < 4; ++k) jacbak[j][k] = jac[j][k];
252 }
253 detbak = det;
254 t1bak = t1;
255 t2bak = t2;
256 t3bak = t3;
257 t4bak = t4;
258 imapbak = imap;
259 if (m_debug) {
260 PrintElement("FindElement5", x, y, z, t1, t2, t3, t4, imap, 6);
261 }
262 } else {
263 // Non-degenerate element
264 if (Coordinates5(x, y, z, t1, t2, t3, t4, jac, det, idxToElemList) != 0) {
265 continue;
266 }
267 if (t1 < -1 || t1 > 1 || t2 < -1 || t2 > 1) continue;
268 ++nfound;
269 imap = idxToElemList;
270 m_lastElement = idxToElemList;
271 if (m_debug) {
272 std::cout << m_className << "::FindElement5:\n";
273 std::cout << " Found matching non-degenerate element "
274 << idxToElemList << ".\n";
275 }
276 if (!m_checkMultipleElement) return idxToElemList;
277 for (int j = 0; j < 4; ++j) {
278 for (int k = 0; k < 4; ++k) jacbak[j][k] = jac[j][k];
279 }
280 detbak = det;
281 t1bak = t1;
282 t2bak = t2;
283 t3bak = t3;
284 t4bak = t4;
285 imapbak = imap;
286 if (m_debug) {
287 PrintElement("FindElement5", x, y, z, t1, t2, t3, t4, imap, 8);
288 }
289 }
290 }
291
292 // In checking mode, verify the tetrahedron/triangle count.
293 if (m_checkMultipleElement) {
294 if (nfound < 1) {
295 if (m_debug) {
296 std::cout << m_className << "::FindElement5:\n";
297 std::cout << " No element matching point (" << x << ", " << y
298 << ") found.\n";
299 }
300 m_lastElement = -1;
301 return -1;
302 }
303 if (nfound > 1) {
304 std::cout << m_className << "::FindElement5:\n";
305 std::cout << " Found " << nfound << " elements matching point (" << x
306 << ", " << y << ").\n";
307 }
308 if (nfound > 0) {
309 for (int j = 0; j < 4; ++j) {
310 for (int k = 0; k < 4; ++k) jac[j][k] = jacbak[j][k];
311 }
312 det = detbak;
313 t1 = t1bak;
314 t2 = t2bak;
315 t3 = t3bak;
316 t4 = t4bak;
317 imap = imapbak;
318 m_lastElement = imap;
319 return imap;
320 }
321 }
322
323 if (m_debug) {
324 std::cout << m_className << "::FindElement5:\n";
325 std::cout << " No element matching point (" << x << ", " << y
326 << ") found.\n";
327 }
328 return -1;
329}
330
331int ComponentFieldMap::FindElement13(const double x, const double y,
332 const double z, double& t1, double& t2,
333 double& t3, double& t4, double jac[4][4],
334 double& det) {
335 // Check if bounding boxes of elements have been computed
336 if (!m_cacheElemBoundingBoxes) {
337 std::cout << m_className << "::FindElement13:\n"
338 << " Caching the bounding boxes of all elements...";
339 CalculateElementBoundingBoxes();
340 std::cout << " done.\n";
341 m_cacheElemBoundingBoxes = true;
342 }
343
344 // Backup
345 double jacbak[4][4];
346 double detbak = 1.;
347 double t1bak = 0., t2bak = 0., t3bak = 0., t4bak = 0.;
348 int imapbak = -1;
349
350 // Initial values.
351 t1 = t2 = t3 = t4 = 0.;
352
353 // Check previously used element
354 if (m_lastElement > -1 && !m_checkMultipleElement) {
355 if (Coordinates13(x, y, z, t1, t2, t3, t4, jac, det, m_lastElement) == 0) {
356 if (t1 >= 0 && t1 <= +1 && t2 >= 0 && t2 <= +1 && t3 >= 0 && t3 <= +1 &&
357 t4 >= 0 && t4 <= +1) {
358 return m_lastElement;
359 }
360 }
361 }
362
363 // Tetra list in the block that contains the input 3D point.
364 std::vector<int> tetList;
365 if (m_useTetrahedralTree) {
366 if (!m_isTreeInitialized) {
367 if (!InitializeTetrahedralTree()) {
368 std::cerr << m_className << "::FindElement13:\n";
369 std::cerr << " Tetrahedral tree initialization failed.\n";
370 return -1;
371 }
372 }
373 tetList = m_tetTree->GetTetListInBlock(Vec3(x, y, z));
374 }
375 // Number of elements to scan.
376 // With tetra tree disabled, all elements are scanned.
377 const int numElemToSearch =
378 m_useTetrahedralTree ? tetList.size() : nElements;
379 // Verify the count of volumes that contain the point.
380 int nfound = 0;
381 int imap = -1;
382
383 // Scan all elements
384 for (int i = 0; i < numElemToSearch; i++) {
385 const int idxToElemList = m_useTetrahedralTree ? tetList[i] : i;
386 const Element& element = elements[idxToElemList];
387 // Tolerance
388 const double f = 0.2;
389 const double tolx = f * (element.xmax - element.xmin);
390 if (x < element.xmin - tolx || x > element.xmax + tolx) continue;
391 const double toly = f * (element.ymax - element.ymin);
392 if (y < element.ymin - toly || y > element.ymax + toly) continue;
393 const double tolz = f * (element.zmax - element.zmin);
394 if (z < element.zmin - tolz || z > element.zmax + tolz) continue;
395 if (Coordinates13(x, y, z, t1, t2, t3, t4, jac, det, idxToElemList) != 0) {
396 continue;
397 }
398 if (t1 < 0 || t1 > 1 || t2 < 0 || t2 > 1 || t3 < 0 || t3 > 1 || t4 < 0 ||
399 t4 > 1) {
400 continue;
401 }
402 ++nfound;
403 imap = idxToElemList;
404 m_lastElement = idxToElemList;
405 if (m_debug) {
406 std::cout << m_className << "::FindElement13:\n";
407 std::cout << " Found matching element " << i << ".\n";
408 }
409 if (!m_checkMultipleElement) return idxToElemList;
410 for (int j = 0; j < 4; ++j) {
411 for (int k = 0; k < 4; ++k) jacbak[j][k] = jac[j][k];
412 }
413 detbak = det;
414 t1bak = t1;
415 t2bak = t2;
416 t3bak = t3;
417 t4bak = t4;
418 imapbak = imap;
419 if (m_debug) {
420 PrintElement("FindElement13", x, y, z, t1, t2, t3, t4, imap, 10);
421 }
422 }
423
424 // In checking mode, verify the tetrahedron/triangle count.
425 if (m_checkMultipleElement) {
426 if (nfound < 1) {
427 if (m_debug) {
428 std::cout << m_className << "::FindElement13:\n";
429 std::cout << " No element matching point (" << x << ", " << y << ", "
430 << z << ") found.\n";
431 }
432 m_lastElement = -1;
433 return -1;
434 }
435 if (nfound > 1) {
436 std::cerr << m_className << "::FindElement13:\n";
437 std::cerr << " Found << " << nfound << " elements matching point ("
438 << x << ", " << y << ", " << z << ").\n";
439 }
440 if (nfound > 0) {
441 for (int j = 0; j < 4; ++j) {
442 for (int k = 0; k < 4; ++k) jac[j][k] = jacbak[j][k];
443 }
444 det = detbak;
445 t1 = t1bak;
446 t2 = t2bak;
447 t3 = t3bak;
448 t4 = t4bak;
449 imap = imapbak;
450 m_lastElement = imap;
451 return imap;
452 }
453 }
454
455 if (m_debug) {
456 std::cout << m_className << "::FindElement13:\n";
457 std::cout << " No element matching point (" << x << ", " << y << ", "
458 << z << ") found.\n";
459 }
460 return -1;
461}
462
463int ComponentFieldMap::FindElementCube(const double x, const double y,
464 const double z, double& t1, double& t2,
465 double& t3, TMatrixD*& jac,
466 std::vector<TMatrixD*>& dN) {
467
468 int imap = -1;
469 if (m_lastElement >= 0) {
470 const Element& element = elements[m_lastElement];
471 const Node& n3 = nodes[element.emap[3]];
472 if (x >= n3.x && y >= n3.y && z >= n3.z) {
473 const Node& n0 = nodes[element.emap[0]];
474 const Node& n2 = nodes[element.emap[2]];
475 const Node& n7 = nodes[element.emap[7]];
476 if (x < n0.x && y < n2.y && z < n7.z) {
477 imap = m_lastElement;
478 }
479 }
480 }
481
482 // Default element loop
483 if (imap == -1) {
484 for (int i = 0; i < nElements; ++i) {
485 const Element& element = elements[i];
486 const Node& n3 = nodes[element.emap[3]];
487 if (x < n3.x || y < n3.y || z < n3.z) continue;
488 const Node& n0 = nodes[element.emap[0]];
489 const Node& n2 = nodes[element.emap[2]];
490 const Node& n7 = nodes[element.emap[7]];
491 if (x < n0.x && y < n2.y && z < n7.z) {
492 imap = i;
493 break;
494 }
495 }
496 }
497
498 if (imap < 0) {
499 if (m_debug) {
500 std::cout << m_className << "::FindElementCube:\n";
501 std::cout << " Point (" << x << "," << y << "," << z
502 << ") not in the mesh, it is background or PEC.\n";
503 const Node& first0 = nodes[elements.front().emap[0]];
504 const Node& first2 = nodes[elements.front().emap[2]];
505 const Node& first3 = nodes[elements.front().emap[3]];
506 const Node& first7 = nodes[elements.front().emap[7]];
507 std::cout << " First node (" << first3.x << ","
508 << first3.y << "," << first3.z << ") in the mesh.\n";
509 std::cout << " dx= " << (first0.x - first3.x)
510 << ", dy= " << (first2.y - first3.y)
511 << ", dz= " << (first7.z - first3.z) << "\n";
512 const Node& last0 = nodes[elements.back().emap[0]];
513 const Node& last2 = nodes[elements.back().emap[2]];
514 const Node& last3 = nodes[elements.back().emap[3]];
515 const Node& last5 = nodes[elements.back().emap[5]];
516 const Node& last7 = nodes[elements.back().emap[7]];
517 std::cout << " Last node (" << last5.x << ","
518 << last5.y << "," << last5.z << ") in the mesh.\n";
519 std::cout << " dx= " << (last0.x - last3.x)
520 << ", dy= " << (last2.y - last3.y)
521 << ", dz= " << (last7.z - last3.z) << "\n";
522 }
523 return -1;
524 }
525 CoordinatesCube(x, y, z, t1, t2, t3, jac, dN, imap);
526 if (m_debug) {
527 PrintElement("FindElementCube", x, y, z, t1, t2, t3, 0., imap, 8);
528 }
529 return imap;
530}
531
532void ComponentFieldMap::Jacobian3(const unsigned int i, const double u,
533 const double v, const double w, double& det,
534 double jac[4][4]) const {
535
536 // Initial values
537 det = 0;
538 jac[0][0] = 0;
539 jac[0][1] = 0;
540 jac[1][0] = 0;
541 jac[1][1] = 0;
542
543 const Element& element = elements[i];
544 const Node& n0 = nodes[element.emap[0]];
545 const Node& n1 = nodes[element.emap[1]];
546 const Node& n2 = nodes[element.emap[2]];
547 const Node& n3 = nodes[element.emap[3]];
548 const Node& n4 = nodes[element.emap[4]];
549 const Node& n5 = nodes[element.emap[5]];
550
551 // Determinant of the quadratic triangular Jacobian
552 det = -(((-1 + 4 * v) * n1.x + n2.x - 4 * w * n2.x + 4 * u * n3.x -
553 4 * u * n4.x - 4 * v * n5.x + 4 * w * n5.x) *
554 (-n0.y + 4 * u * n0.y + 4 * v * n3.y + 4 * w * n4.y)) -
555 ((-1 + 4 * u) * n0.x + n1.x - 4 * v * n1.x - 4 * u * n3.x +
556 4 * v * n3.x + 4 * w * n4.x - 4 * w * n5.x) *
557 (-n2.y + 4 * w * n2.y + 4 * u * n4.y + 4 * v * n5.y) +
558 ((-1 + 4 * u) * n0.x + n2.x - 4 * w * n2.x + 4 * v * n3.x -
559 4 * u * n4.x + 4 * w * n4.x - 4 * v * n5.x) *
560 (-n1.y + 4 * v * n1.y + 4 * u * n3.y + 4 * w * n5.y);
561
562 // Terms of the quadratic triangular Jacobian
563 jac[0][0] = (-n1.x + 4 * v * n1.x + 4 * u * n3.x + 4 * w * n5.x) *
564 (-n2.y + 4 * w * n2.y + 4 * u * n4.y + 4 * v * n5.y) -
565 (-n2.x + 4 * w * n2.x + 4 * u * n4.x + 4 * v * n5.x) *
566 (-n1.y + 4 * v * n1.y + 4 * u * n3.y + 4 * w * n5.y);
567 jac[0][1] = (-1 + 4 * v) * n1.y + n2.y - 4 * w * n2.y + 4 * u * n3.y -
568 4 * u * n4.y - 4 * v * n5.y + 4 * w * n5.y;
569 jac[0][2] = n1.x - 4 * v * n1.x + (-1 + 4 * w) * n2.x - 4 * u * n3.x +
570 4 * u * n4.x + 4 * v * n5.x - 4 * w * n5.x;
571 jac[1][0] = (-n2.x + 4 * w * n2.x + 4 * u * n4.x + 4 * v * n5.x) *
572 (-n0.y + 4 * u * n0.y + 4 * v * n3.y + 4 * w * n4.y) -
573 (-n0.x + 4 * u * n0.x + 4 * v * n3.x + 4 * w * n4.x) *
574 (-n2.y + 4 * w * n2.y + 4 * u * n4.y + 4 * v * n5.y);
575 jac[1][1] = n0.y - 4 * u * n0.y - n2.y + 4 * w * n2.y - 4 * v * n3.y +
576 4 * u * n4.y - 4 * w * n4.y + 4 * v * n5.y;
577 jac[1][2] = (-1 + 4 * u) * n0.x + n2.x - 4 * w * n2.x + 4 * v * n3.x -
578 4 * u * n4.x + 4 * w * n4.x - 4 * v * n5.x;
579 jac[2][0] = -((-n1.x + 4 * v * n1.x + 4 * u * n3.x + 4 * w * n5.x) *
580 (-n0.y + 4 * u * n0.y + 4 * v * n3.y + 4 * w * n4.y)) +
581 (-n0.x + 4 * u * n0.x + 4 * v * n3.x + 4 * w * n4.x) *
582 (-n1.y + 4 * v * n1.y + 4 * u * n3.y + 4 * w * n5.y);
583 jac[2][1] = (-1 + 4 * u) * n0.y + n1.y - 4 * v * n1.y - 4 * u * n3.y +
584 4 * v * n3.y + 4 * w * n4.y - 4 * w * n5.y;
585 jac[2][2] = n0.x - 4 * u * n0.x - n1.x + 4 * v * n1.x + 4 * u * n3.x -
586 4 * v * n3.x - 4 * w * n4.x + 4 * w * n5.x;
587}
588
589void ComponentFieldMap::Jacobian5(const unsigned int i, const double u,
590 const double v, double& det,
591 double jac[4][4]) const {
592
593 // Initial values
594 det = 0;
595 jac[0][0] = 0;
596 jac[0][1] = 0;
597 jac[1][0] = 0;
598 jac[1][1] = 0;
599
600 const Element& element = elements[i];
601 const Node& n0 = nodes[element.emap[0]];
602 const Node& n1 = nodes[element.emap[1]];
603 const Node& n2 = nodes[element.emap[2]];
604 const Node& n3 = nodes[element.emap[3]];
605 const Node& n4 = nodes[element.emap[4]];
606 const Node& n5 = nodes[element.emap[5]];
607 const Node& n6 = nodes[element.emap[6]];
608 const Node& n7 = nodes[element.emap[7]];
609 // Determinant of the quadrilateral serendipity Jacobian
610 det =
611 (-2 * u * u * u * ((n2.x + n3.x - 2 * n6.x) * (n0.y + n1.y - 2 * n4.y) -
612 (n0.x + n1.x - 2 * n4.x) * (n2.y + n3.y - 2 * n6.y)) +
613 2 * v * v * v * (-((n0.x + n3.x - 2 * n7.x) * (n1.y + n2.y - 2 * n5.y)) +
614 (n1.x + n2.x - 2 * n5.x) * (n0.y + n3.y - 2 * n7.y)) +
615 2 * (-((n5.x - n7.x) * (n4.y - n6.y)) + (n4.x - n6.x) * (n5.y - n7.y)) +
616 v * (-(n6.x * n0.y) - 2 * n7.x * n0.y + n6.x * n1.y - 2 * n7.x * n1.y -
617 n6.x * n2.y - 2 * n7.x * n2.y + n4.x * (n0.y - n1.y + n2.y - n3.y) +
618 n6.x * n3.y - 2 * n7.x * n3.y - n0.x * n4.y + n1.x * n4.y -
619 n2.x * n4.y + n3.x * n4.y - 2 * n0.x * n5.y - 2 * n1.x * n5.y -
620 2 * n2.x * n5.y - 2 * n3.x * n5.y + 8 * n7.x * n5.y + n0.x * n6.y -
621 n1.x * n6.y + n2.x * n6.y - n3.x * n6.y +
622 2 * n5.x * (n0.y + n1.y + n2.y + n3.y - 4 * n7.y) +
623 2 * (n0.x + n1.x + n2.x + n3.x) * n7.y) +
624 v * v *
625 (-(n4.x * n0.y) + 2 * n5.x * n0.y + n6.x * n0.y + 2 * n7.x * n0.y +
626 n4.x * n1.y - 2 * n5.x * n1.y - n6.x * n1.y - 2 * n7.x * n1.y +
627 n4.x * n2.y + 2 * n5.x * n2.y - n6.x * n2.y + 2 * n7.x * n2.y -
628 n4.x * n3.y - 2 * n5.x * n3.y + n6.x * n3.y - 2 * n7.x * n3.y +
629 2 * n2.x * (n1.y + n3.y) - n2.x * n4.y + 2 * n5.x * n4.y -
630 2 * n7.x * n4.y - 2 * n2.x * n5.y - 2 * n4.x * n5.y +
631 2 * n6.x * n5.y + n2.x * n6.y - 2 * n5.x * n6.y + 2 * n7.x * n6.y +
632 n0.x * (2 * n1.y + 2 * n3.y + n4.y - 2 * n5.y - n6.y - 2 * n7.y) -
633 2 * (n2.x - n4.x + n6.x) * n7.y +
634 n3.x * (-2 * n0.y - 2 * n2.y + n4.y + 2 * n5.y - n6.y + 2 * n7.y) +
635 n1.x * (-2 * n0.y - 2 * n2.y - n4.y + 2 * n5.y + n6.y + 2 * n7.y)) +
636 u * (n5.x * n0.y - 2 * n6.x * n0.y - n7.x * n0.y - n5.x * n1.y -
637 2 * n6.x * n1.y + n7.x * n1.y + n5.x * n2.y - 2 * n6.x * n2.y -
638 n7.x * n2.y - n5.x * n3.y - 2 * n6.x * n3.y + n7.x * n3.y -
639 2 * n1.x * n4.y - 2 * n2.x * n4.y - 2 * n3.x * n4.y +
640 8 * n6.x * n4.y + n1.x * n5.y - n2.x * n5.y + n3.x * n5.y +
641 2 * n4.x * (n0.y + n1.y + n2.y + n3.y - 4 * n6.y) +
642 2 * n1.x * n6.y + 2 * n2.x * n6.y + 2 * n3.x * n6.y -
643 (n1.x - n2.x + n3.x) * n7.y +
644 n0.x * (-2 * n4.y - n5.y + 2 * n6.y + n7.y) +
645 v * v * (4 * n4.x * n0.y - 3 * n5.x * n0.y - 4 * n6.x * n0.y -
646 5 * n7.x * n0.y + 4 * n4.x * n1.y - 5 * n5.x * n1.y -
647 4 * n6.x * n1.y - 3 * n7.x * n1.y + 4 * n4.x * n2.y +
648 5 * n5.x * n2.y - 4 * n6.x * n2.y + 3 * n7.x * n2.y +
649 4 * n4.x * n3.y + 3 * n5.x * n3.y - 4 * n6.x * n3.y +
650 5 * n7.x * n3.y + 8 * n5.x * n4.y + 8 * n7.x * n4.y -
651 8 * n4.x * n5.y + 8 * n6.x * n5.y - 8 * n5.x * n6.y -
652 8 * n7.x * n6.y + n3.x * (5 * n0.y + 3 * n1.y - 4 * n4.y -
653 3 * n5.y + 4 * n6.y - 5 * n7.y) +
654 n2.x * (3 * n0.y + 5 * n1.y - 4 * n4.y - 5 * n5.y +
655 4 * n6.y - 3 * n7.y) -
656 8 * n4.x * n7.y + 8 * n6.x * n7.y +
657 n1.x * (-5 * n2.y - 3 * n3.y - 4 * n4.y + 5 * n5.y +
658 4 * n6.y + 3 * n7.y) +
659 n0.x * (-3 * n2.y - 5 * n3.y - 4 * n4.y + 3 * n5.y +
660 4 * n6.y + 5 * n7.y)) -
661 2 * v * (n6.x * n0.y - 3 * n7.x * n0.y + n6.x * n1.y - n7.x * n1.y +
662 3 * n6.x * n2.y - n7.x * n2.y + 3 * n6.x * n3.y -
663 3 * n7.x * n3.y - 3 * n0.x * n4.y - 3 * n1.x * n4.y -
664 n2.x * n4.y - n3.x * n4.y + 4 * n7.x * n4.y + n0.x * n5.y +
665 3 * n1.x * n5.y + 3 * n2.x * n5.y + n3.x * n5.y -
666 4 * n6.x * n5.y - n0.x * n6.y - n1.x * n6.y -
667 3 * n2.x * n6.y - 3 * n3.x * n6.y + 4 * n7.x * n6.y -
668 n5.x * (n0.y + 3 * n1.y + 3 * n2.y + n3.y -
669 4 * (n4.y + n6.y)) +
670 (3 * n0.x + n1.x + n2.x + 3 * n3.x - 4 * n6.x) * n7.y +
671 n4.x * (3 * n0.y + 3 * n1.y + n2.y + n3.y -
672 4 * (n5.y + n7.y)))) +
673 u * u *
674 (2 * n3.x * n0.y - 2 * n4.x * n0.y - n5.x * n0.y - 2 * n6.x * n0.y +
675 n7.x * n0.y - 2 * n0.x * n1.y + 2 * n4.x * n1.y - n5.x * n1.y +
676 2 * n6.x * n1.y + n7.x * n1.y + 2 * n3.x * n2.y - 2 * n4.x * n2.y +
677 n5.x * n2.y - 2 * n6.x * n2.y - n7.x * n2.y + 2 * n4.x * n3.y +
678 n5.x * n3.y + 2 * n6.x * n3.y - n7.x * n3.y - 2 * n3.x * n4.y +
679 2 * n5.x * n4.y - 2 * n7.x * n4.y - n3.x * n5.y - 2 * n4.x * n5.y +
680 2 * n6.x * n5.y - 2 * n3.x * n6.y - 2 * n5.x * n6.y +
681 2 * n7.x * n6.y +
682 n0.x * (-2 * n3.y + 2 * n4.y + n5.y + 2 * n6.y - n7.y) +
683 (n3.x + 2 * n4.x - 2 * n6.x) * n7.y +
684 n2.x * (-2 * n1.y - 2 * n3.y + 2 * n4.y - n5.y + 2 * n6.y + n7.y) -
685 3 * v * v *
686 (n5.x * n0.y - n6.x * n0.y - n7.x * n0.y + n5.x * n1.y +
687 n6.x * n1.y - n7.x * n1.y - n5.x * n2.y + n6.x * n2.y +
688 n7.x * n2.y - n5.x * n3.y - n6.x * n3.y + n7.x * n3.y -
689 2 * n5.x * n4.y + 2 * n7.x * n4.y - 2 * n6.x * n5.y +
690 2 * n5.x * n6.y - 2 * n7.x * n6.y +
691 n4.x * (n0.y - n1.y - n2.y + n3.y + 2 * n5.y - 2 * n7.y) +
692 n3.x * (n0.y - n2.y - n4.y + n5.y + n6.y - n7.y) +
693 2 * n6.x * n7.y +
694 (n0.x - n2.x) * (n1.y - n3.y - n4.y - n5.y + n6.y + n7.y)) +
695 v * (4 * n5.x * n0.y + 3 * n6.x * n0.y - 4 * n7.x * n0.y +
696 4 * n5.x * n1.y - 3 * n6.x * n1.y - 4 * n7.x * n1.y +
697 4 * n5.x * n2.y - 5 * n6.x * n2.y - 4 * n7.x * n2.y +
698 4 * n5.x * n3.y + 5 * n6.x * n3.y - 4 * n7.x * n3.y -
699 8 * n5.x * n4.y + 8 * n7.x * n4.y + 8 * n6.x * n5.y -
700 8 * n5.x * n6.y + 8 * n7.x * n6.y +
701 n4.x * (5 * n0.y - 5 * n1.y - 3 * n2.y + 3 * n3.y + 8 * n5.y -
702 8 * n7.y) -
703 8 * n6.x * n7.y + n3.x * (3 * n1.y + 5 * n2.y - 3 * n4.y -
704 4 * n5.y - 5 * n6.y + 4 * n7.y) +
705 n0.x * (5 * n1.y + 3 * n2.y - 5 * n4.y - 4 * n5.y - 3 * n6.y +
706 4 * n7.y) +
707 n2.x * (-3 * n0.y - 5 * n3.y + 3 * n4.y - 4 * n5.y + 5 * n6.y +
708 4 * n7.y)) +
709 n1.x * ((-1 + v) * (-2 + 3 * v) * n0.y + 2 * n2.y - 2 * n4.y +
710 n5.y - 2 * n6.y - n7.y +
711 v * (-3 * n3.y + 5 * n4.y - 4 * n5.y + 3 * n6.y + 4 * n7.y -
712 3 * v * (n2.y + n4.y - n5.y - n6.y + n7.y))))) /
713 8;
714 // Jacobian terms
715 jac[0][0] =
716 (u * u * (-n0.y - n1.y + n2.y + n3.y + 2 * n4.y - 2 * n6.y) +
717 2 * (-n4.y + n6.y +
718 v * (n0.y + n1.y + n2.y + n3.y - 2 * n5.y - 2 * n7.y)) +
719 u * (n0.y - 2 * v * n0.y - n1.y + 2 * v * n1.y + n2.y + 2 * v * n2.y -
720 n3.y - 2 * v * n3.y - 4 * v * n5.y + 4 * v * n7.y)) /
721 4;
722 jac[0][1] =
723 (u * u * (n0.x + n1.x - n2.x - n3.x - 2 * n4.x + 2 * n6.x) -
724 2 * (-n4.x + n6.x +
725 v * (n0.x + n1.x + n2.x + n3.x - 2 * n5.x - 2 * n7.x)) +
726 u * ((-1 + 2 * v) * n0.x + n1.x - 2 * v * n1.x - n2.x - 2 * v * n2.x +
727 n3.x + 2 * v * n3.x + 4 * v * n5.x - 4 * v * n7.x)) /
728 4;
729 jac[1][0] =
730 (v * (-n0.y + n1.y - n2.y + n3.y) - 2 * n5.y +
731 2 * u * ((-1 + v) * n0.y + (-1 + v) * n1.y - n2.y - v * n2.y - n3.y -
732 v * n3.y + 2 * n4.y - 2 * v * n4.y + 2 * n6.y + 2 * v * n6.y) +
733 v * v * (n0.y - n1.y - n2.y + n3.y + 2 * n5.y - 2 * n7.y) + 2 * n7.y) /
734 4;
735 jac[1][1] =
736 (v * (n0.x - n1.x + n2.x - n3.x) +
737 2 * u * (n0.x - v * n0.x + n1.x - v * n1.x + n2.x + v * n2.x + n3.x +
738 v * n3.x - 2 * n4.x + 2 * v * n4.x - 2 * n6.x - 2 * v * n6.x) +
739 2 * (n5.x - n7.x) +
740 v * v * (-n0.x + n1.x + n2.x - n3.x - 2 * n5.x + 2 * n7.x)) /
741 4;
742}
743
744void ComponentFieldMap::Jacobian13(const unsigned int i, const double t,
745 const double u, const double v,
746 const double w, double& det,
747 double jac[4][4]) const {
748
749 // Initial values
750 det = 0;
751 for (int j = 0; j < 4; ++j) {
752 for (int k = 0; k < 4; ++k) jac[j][k] = 0;
753 }
754
755 const Element& element = elements[i];
756 const Node& n0 = nodes[element.emap[0]];
757 const Node& n1 = nodes[element.emap[1]];
758 const Node& n2 = nodes[element.emap[2]];
759 const Node& n3 = nodes[element.emap[3]];
760 const Node& n4 = nodes[element.emap[4]];
761 const Node& n5 = nodes[element.emap[5]];
762 const Node& n6 = nodes[element.emap[6]];
763 const Node& n7 = nodes[element.emap[7]];
764 const Node& n8 = nodes[element.emap[8]];
765 const Node& n9 = nodes[element.emap[9]];
766 // Determinant of the quadrilateral serendipity Jacobian
767 det =
768 -(((-4 * v * n9.x - n1.x + 4 * u * n1.x + n3.x - 4 * w * n3.x +
769 4 * t * n4.x - 4 * t * n6.x + 4 * v * n7.x - 4 * u * n8.x +
770 4 * w * n8.x) *
771 (4 * w * n9.y - n2.y + 4 * v * n2.y + 4 * t * n5.y +
772 4 * u * n7.y) +
773 (n1.x - 4 * u * n1.x - n2.x + 4 * v * n2.x - 4 * t * n4.x +
774 4 * t * n5.x + 4 * u * n7.x - 4 * v * n7.x + 4 * w * (n9.x - n8.x)) *
775 (4 * v * n9.y - n3.y + 4 * w * n3.y + 4 * t * n6.y +
776 4 * u * n8.y) +
777 (-4 * w * n9.x + 4 * v * (n9.x - n2.x) + n2.x - n3.x + 4 * w * n3.x -
778 4 * t * n5.x + 4 * t * n6.x - 4 * u * n7.x + 4 * u * n8.x) *
779 ((-1 + 4 * u) * n1.y + 4 * (t * n4.y + v * n7.y + w * n8.y))) *
780 ((-1 + 4 * t) * n0.z + 4 * (u * n4.z + v * n5.z + w * n6.z))) -
781 ((n1.x - 4 * u * n1.x - n3.x + 4 * w * n3.x - 4 * t * n4.x +
782 4 * t * n6.x + 4 * v * (n9.x - n7.x) + 4 * u * n8.x - 4 * w * n8.x) *
783 ((-1 + 4 * t) * n0.y + 4 * (u * n4.y + v * n5.y + w * n6.y)) -
784 ((-1 + 4 * t) * n0.x + n1.x - 4 * u * n1.x +
785 4 * (-(t * n4.x) + u * n4.x + v * n5.x + w * n6.x - v * n7.x -
786 w * n8.x)) *
787 (4 * v * n9.y - n3.y + 4 * w * n3.y + 4 * t * n6.y + 4 * u * n8.y) +
788 ((-1 + 4 * t) * n0.x - 4 * v * n9.x + n3.x - 4 * w * n3.x +
789 4 * u * n4.x + 4 * v * n5.x - 4 * t * n6.x + 4 * w * n6.x -
790 4 * u * n8.x) *
791 ((-1 + 4 * u) * n1.y + 4 * (t * n4.y + v * n7.y + w * n8.y))) *
792 (4 * w * n9.z - n2.z + 4 * v * n2.z + 4 * t * n5.z + 4 * u * n7.z) +
793 ((n1.x - 4 * u * n1.x - n2.x + 4 * v * n2.x - 4 * t * n4.x +
794 4 * t * n5.x + 4 * u * n7.x - 4 * v * n7.x + 4 * w * (n9.x - n8.x)) *
795 ((-1 + 4 * t) * n0.y + 4 * (u * n4.y + v * n5.y + w * n6.y)) -
796 ((-1 + 4 * t) * n0.x + n1.x - 4 * u * n1.x +
797 4 * (-(t * n4.x) + u * n4.x + v * n5.x + w * n6.x - v * n7.x -
798 w * n8.x)) *
799 (4 * w * n9.y - n2.y + 4 * v * n2.y + 4 * t * n5.y + 4 * u * n7.y) +
800 ((-1 + 4 * t) * n0.x - 4 * w * n9.x + n2.x - 4 * v * n2.x +
801 4 * u * n4.x - 4 * t * n5.x + 4 * v * n5.x + 4 * w * n6.x -
802 4 * u * n7.x) *
803 ((-1 + 4 * u) * n1.y + 4 * (t * n4.y + v * n7.y + w * n8.y))) *
804 (4 * v * n9.z - n3.z + 4 * w * n3.z + 4 * t * n6.z + 4 * u * n8.z) +
805 ((-4 * w * n9.x + 4 * v * (n9.x - n2.x) + n2.x - n3.x + 4 * w * n3.x -
806 4 * t * n5.x + 4 * t * n6.x - 4 * u * n7.x + 4 * u * n8.x) *
807 ((-1 + 4 * t) * n0.y + 4 * (u * n4.y + v * n5.y + w * n6.y)) +
808 ((-1 + 4 * t) * n0.x - 4 * v * n9.x + n3.x - 4 * w * n3.x +
809 4 * u * n4.x + 4 * v * n5.x - 4 * t * n6.x + 4 * w * n6.x -
810 4 * u * n8.x) *
811 (4 * w * n9.y - n2.y + 4 * v * n2.y + 4 * t * n5.y + 4 * u * n7.y) -
812 ((-1 + 4 * t) * n0.x - 4 * w * n9.x + n2.x - 4 * v * n2.x +
813 4 * u * n4.x - 4 * t * n5.x + 4 * v * n5.x + 4 * w * n6.x -
814 4 * u * n7.x) *
815 (4 * v * n9.y - n3.y + 4 * w * n3.y + 4 * t * n6.y + 4 * u * n8.y)) *
816 ((-1 + 4 * u) * n1.z + 4 * (t * n4.z + v * n7.z + w * n8.z));
817
818 jac[0][0] =
819 -((((-1 + 4 * u) * n1.x + 4 * (t * n4.x + v * n7.x + w * n8.x)) *
820 (4 * v * n9.y - n3.y + 4 * w * n3.y + 4 * t * n6.y +
821 4 * u * n8.y) -
822 (4 * v * n9.x - n3.x + 4 * w * n3.x + 4 * t * n6.x + 4 * u * n8.x) *
823 ((-1 + 4 * u) * n1.y + 4 * (t * n4.y + v * n7.y + w * n8.y))) *
824 (4 * w * n9.z - n2.z + 4 * v * n2.z + 4 * t * n5.z + 4 * u * n7.z)) +
825 (((-1 + 4 * u) * n1.x + 4 * (t * n4.x + v * n7.x + w * n8.x)) *
826 (4 * w * n9.y - n2.y + 4 * v * n2.y + 4 * t * n5.y + 4 * u * n7.y) -
827 (4 * w * n9.x - n2.x + 4 * v * n2.x + 4 * t * n5.x + 4 * u * n7.x) *
828 ((-1 + 4 * u) * n1.y + 4 * (t * n4.y + v * n7.y + w * n8.y))) *
829 (4 * v * n9.z - n3.z + 4 * w * n3.z + 4 * t * n6.z + 4 * u * n8.z) +
830 (-((4 * v * n9.x - n3.x + 4 * w * n3.x + 4 * t * n6.x + 4 * u * n8.x) *
831 (4 * w * n9.y - n2.y + 4 * v * n2.y + 4 * t * n5.y + 4 * u * n7.y)) +
832 (4 * w * n9.x - n2.x + 4 * v * n2.x + 4 * t * n5.x + 4 * u * n7.x) *
833 (4 * v * n9.y - n3.y + 4 * w * n3.y + 4 * t * n6.y + 4 * u * n8.y)) *
834 ((-1 + 4 * u) * n1.z + 4 * (t * n4.z + v * n7.z + w * n8.z));
835
836 jac[0][1] =
837 (n1.y - 4 * u * n1.y - n3.y + 4 * w * n3.y - 4 * t * n4.y + 4 * t * n6.y +
838 4 * v * (n9.y - n7.y) + 4 * u * n8.y - 4 * w * n8.y) *
839 (4 * w * n9.z - n2.z + 4 * v * n2.z + 4 * t * n5.z + 4 * u * n7.z) +
840 (-4 * w * n9.y - n1.y + 4 * u * n1.y + n2.y - 4 * v * n2.y +
841 4 * t * n4.y - 4 * t * n5.y - 4 * u * n7.y + 4 * v * n7.y +
842 4 * w * n8.y) *
843 (4 * v * n9.z - n3.z + 4 * w * n3.z + 4 * t * n6.z + 4 * u * n8.z) +
844 (-4 * v * n9.y + 4 * w * n9.y - n2.y + 4 * v * n2.y + n3.y -
845 4 * w * n3.y + 4 * t * n5.y - 4 * t * n6.y + 4 * u * n7.y -
846 4 * u * n8.y) *
847 ((-1 + 4 * u) * n1.z + 4 * (t * n4.z + v * n7.z + w * n8.z));
848
849 jac[0][2] =
850 (-4 * v * n9.x - n1.x + 4 * u * n1.x + n3.x - 4 * w * n3.x +
851 4 * t * n4.x - 4 * t * n6.x + 4 * v * n7.x - 4 * u * n8.x +
852 4 * w * n8.x) *
853 (4 * w * n9.z - n2.z + 4 * v * n2.z + 4 * t * n5.z + 4 * u * n7.z) +
854 (n1.x - 4 * u * n1.x - n2.x + 4 * v * n2.x - 4 * t * n4.x + 4 * t * n5.x +
855 4 * u * n7.x - 4 * v * n7.x + 4 * w * (n9.x - n8.x)) *
856 (4 * v * n9.z - n3.z + 4 * w * n3.z + 4 * t * n6.z + 4 * u * n8.z) +
857 (-4 * w * n9.x + 4 * v * (n9.x - n2.x) + n2.x - n3.x + 4 * w * n3.x -
858 4 * t * n5.x + 4 * t * n6.x - 4 * u * n7.x + 4 * u * n8.x) *
859 ((-1 + 4 * u) * n1.z + 4 * (t * n4.z + v * n7.z + w * n8.z));
860
861 jac[0][3] =
862 (n1.x - 4 * u * n1.x - n3.x + 4 * w * n3.x - 4 * t * n4.x + 4 * t * n6.x +
863 4 * v * (n9.x - n7.x) + 4 * u * n8.x - 4 * w * n8.x) *
864 (4 * w * n9.y - n2.y + 4 * v * n2.y + 4 * t * n5.y + 4 * u * n7.y) +
865 (-4 * w * n9.x - n1.x + 4 * u * n1.x + n2.x - 4 * v * n2.x +
866 4 * t * n4.x - 4 * t * n5.x - 4 * u * n7.x + 4 * v * n7.x +
867 4 * w * n8.x) *
868 (4 * v * n9.y - n3.y + 4 * w * n3.y + 4 * t * n6.y + 4 * u * n8.y) +
869 (-4 * v * n9.x + 4 * w * n9.x - n2.x + 4 * v * n2.x + n3.x -
870 4 * w * n3.x + 4 * t * n5.x - 4 * t * n6.x + 4 * u * n7.x -
871 4 * u * n8.x) *
872 ((-1 + 4 * u) * n1.y + 4 * (t * n4.y + v * n7.y + w * n8.y));
873
874 jac[1][0] =
875 -((-((4 * v * n9.x - n3.x + 4 * w * n3.x + 4 * t * n6.x + 4 * u * n8.x) *
876 (4 * w * n9.y - n2.y + 4 * v * n2.y + 4 * t * n5.y + 4 * u * n7.y)) +
877 (4 * w * n9.x - n2.x + 4 * v * n2.x + 4 * t * n5.x + 4 * u * n7.x) *
878 (4 * v * n9.y - n3.y + 4 * w * n3.y + 4 * t * n6.y +
879 4 * u * n8.y)) *
880 ((-1 + 4 * t) * n0.z + 4 * (u * n4.z + v * n5.z + w * n6.z))) +
881 (-((4 * v * n9.x - n3.x + 4 * w * n3.x + 4 * t * n6.x + 4 * u * n8.x) *
882 ((-1 + 4 * t) * n0.y + 4 * (u * n4.y + v * n5.y + w * n6.y))) +
883 ((-1 + 4 * t) * n0.x + 4 * (u * n4.x + v * n5.x + w * n6.x)) *
884 (4 * v * n9.y - n3.y + 4 * w * n3.y + 4 * t * n6.y + 4 * u * n8.y)) *
885 (4 * w * n9.z - n2.z + 4 * v * n2.z + 4 * t * n5.z + 4 * u * n7.z) -
886 (-((4 * w * n9.x - n2.x + 4 * v * n2.x + 4 * t * n5.x + 4 * u * n7.x) *
887 ((-1 + 4 * t) * n0.y + 4 * (u * n4.y + v * n5.y + w * n6.y))) +
888 ((-1 + 4 * t) * n0.x + 4 * (u * n4.x + v * n5.x + w * n6.x)) *
889 (4 * w * n9.y - n2.y + 4 * v * n2.y + 4 * t * n5.y + 4 * u * n7.y)) *
890 (4 * v * n9.z - n3.z + 4 * w * n3.z + 4 * t * n6.z + 4 * u * n8.z);
891
892 jac[1][1] =
893 (-4 * w * n9.y + 4 * v * (n9.y - n2.y) + n2.y - n3.y + 4 * w * n3.y -
894 4 * t * n5.y + 4 * t * n6.y - 4 * u * n7.y + 4 * u * n8.y) *
895 ((-1 + 4 * t) * n0.z + 4 * (u * n4.z + v * n5.z + w * n6.z)) +
896 ((-1 + 4 * t) * n0.y - 4 * v * n9.y + n3.y - 4 * w * n3.y + 4 * u * n4.y +
897 4 * v * n5.y - 4 * t * n6.y + 4 * w * n6.y - 4 * u * n8.y) *
898 (4 * w * n9.z - n2.z + 4 * v * n2.z + 4 * t * n5.z + 4 * u * n7.z) -
899 ((-1 + 4 * t) * n0.y - 4 * w * n9.y + n2.y - 4 * v * n2.y + 4 * u * n4.y -
900 4 * t * n5.y + 4 * v * n5.y + 4 * w * n6.y - 4 * u * n7.y) *
901 (4 * v * n9.z - n3.z + 4 * w * n3.z + 4 * t * n6.z + 4 * u * n8.z);
902
903 jac[1][2] =
904 (-4 * v * n9.x + 4 * w * n9.x - n2.x + 4 * v * n2.x + n3.x -
905 4 * w * n3.x + 4 * t * n5.x - 4 * t * n6.x + 4 * u * n7.x -
906 4 * u * n8.x) *
907 ((-1 + 4 * t) * n0.z + 4 * (u * n4.z + v * n5.z + w * n6.z)) -
908 ((-1 + 4 * t) * n0.x - 4 * v * n9.x + n3.x - 4 * w * n3.x + 4 * u * n4.x +
909 4 * v * n5.x - 4 * t * n6.x + 4 * w * n6.x - 4 * u * n8.x) *
910 (4 * w * n9.z - n2.z + 4 * v * n2.z + 4 * t * n5.z + 4 * u * n7.z) +
911 ((-1 + 4 * t) * n0.x - 4 * w * n9.x + n2.x - 4 * v * n2.x + 4 * u * n4.x -
912 4 * t * n5.x + 4 * v * n5.x + 4 * w * n6.x - 4 * u * n7.x) *
913 (4 * v * n9.z - n3.z + 4 * w * n3.z + 4 * t * n6.z + 4 * u * n8.z);
914
915 jac[1][3] =
916 (-4 * w * n9.x + 4 * v * (n9.x - n2.x) + n2.x - n3.x + 4 * w * n3.x -
917 4 * t * n5.x + 4 * t * n6.x - 4 * u * n7.x + 4 * u * n8.x) *
918 ((-1 + 4 * t) * n0.y + 4 * (u * n4.y + v * n5.y + w * n6.y)) +
919 ((-1 + 4 * t) * n0.x - 4 * v * n9.x + n3.x - 4 * w * n3.x + 4 * u * n4.x +
920 4 * v * n5.x - 4 * t * n6.x + 4 * w * n6.x - 4 * u * n8.x) *
921 (4 * w * n9.y - n2.y + 4 * v * n2.y + 4 * t * n5.y + 4 * u * n7.y) -
922 ((-1 + 4 * t) * n0.x - 4 * w * n9.x + n2.x - 4 * v * n2.x + 4 * u * n4.x -
923 4 * t * n5.x + 4 * v * n5.x + 4 * w * n6.x - 4 * u * n7.x) *
924 (4 * v * n9.y - n3.y + 4 * w * n3.y + 4 * t * n6.y + 4 * u * n8.y);
925
926 jac[2][0] =
927 (((-1 + 4 * u) * n1.x + 4 * (t * n4.x + v * n7.x + w * n8.x)) *
928 (4 * v * n9.y - n3.y + 4 * w * n3.y + 4 * t * n6.y + 4 * u * n8.y) -
929 (4 * v * n9.x - n3.x + 4 * w * n3.x + 4 * t * n6.x + 4 * u * n8.x) *
930 ((-1 + 4 * u) * n1.y + 4 * (t * n4.y + v * n7.y + w * n8.y))) *
931 ((-1 + 4 * t) * n0.z + 4 * (u * n4.z + v * n5.z + w * n6.z)) +
932 (-(((-1 + 4 * u) * n1.x + 4 * (t * n4.x + v * n7.x + w * n8.x)) *
933 ((-1 + 4 * t) * n0.y + 4 * (u * n4.y + v * n5.y + w * n6.y))) +
934 ((-1 + 4 * t) * n0.x + 4 * (u * n4.x + v * n5.x + w * n6.x)) *
935 ((-1 + 4 * u) * n1.y + 4 * (t * n4.y + v * n7.y + w * n8.y))) *
936 (4 * v * n9.z - n3.z + 4 * w * n3.z + 4 * t * n6.z + 4 * u * n8.z) -
937 (-((4 * v * n9.x - n3.x + 4 * w * n3.x + 4 * t * n6.x + 4 * u * n8.x) *
938 ((-1 + 4 * t) * n0.y + 4 * (u * n4.y + v * n5.y + w * n6.y))) +
939 ((-1 + 4 * t) * n0.x + 4 * (u * n4.x + v * n5.x + w * n6.x)) *
940 (4 * v * n9.y - n3.y + 4 * w * n3.y + 4 * t * n6.y + 4 * u * n8.y)) *
941 ((-1 + 4 * u) * n1.z + 4 * (t * n4.z + v * n7.z + w * n8.z));
942
943 jac[2][1] =
944 (-4 * v * n9.y - n1.y + 4 * u * n1.y + n3.y - 4 * w * n3.y +
945 4 * t * n4.y - 4 * t * n6.y + 4 * v * n7.y - 4 * u * n8.y +
946 4 * w * n8.y) *
947 ((-1 + 4 * t) * n0.z + 4 * (u * n4.z + v * n5.z + w * n6.z)) +
948 ((-1 + 4 * t) * n0.y + n1.y - 4 * u * n1.y +
949 4 * (-(t * n4.y) + u * n4.y + v * n5.y + w * n6.y - v * n7.y -
950 w * n8.y)) *
951 (4 * v * n9.z - n3.z + 4 * w * n3.z + 4 * t * n6.z + 4 * u * n8.z) -
952 ((-1 + 4 * t) * n0.y - 4 * v * n9.y + n3.y - 4 * w * n3.y + 4 * u * n4.y +
953 4 * v * n5.y - 4 * t * n6.y + 4 * w * n6.y - 4 * u * n8.y) *
954 ((-1 + 4 * u) * n1.z + 4 * (t * n4.z + v * n7.z + w * n8.z));
955
956 jac[2][2] =
957 (n1.x - 4 * u * n1.x - n3.x + 4 * w * n3.x - 4 * t * n4.x + 4 * t * n6.x +
958 4 * v * (n9.x - n7.x) + 4 * u * n8.x - 4 * w * n8.x) *
959 ((-1 + 4 * t) * n0.z + 4 * (u * n4.z + v * n5.z + w * n6.z)) -
960 ((-1 + 4 * t) * n0.x + n1.x - 4 * u * n1.x +
961 4 * (-(t * n4.x) + u * n4.x + v * n5.x + w * n6.x - v * n7.x -
962 w * n8.x)) *
963 (4 * v * n9.z - n3.z + 4 * w * n3.z + 4 * t * n6.z + 4 * u * n8.z) +
964 ((-1 + 4 * t) * n0.x - 4 * v * n9.x + n3.x - 4 * w * n3.x + 4 * u * n4.x +
965 4 * v * n5.x - 4 * t * n6.x + 4 * w * n6.x - 4 * u * n8.x) *
966 ((-1 + 4 * u) * n1.z + 4 * (t * n4.z + v * n7.z + w * n8.z));
967
968 jac[2][3] =
969 (-4 * v * n9.x - n1.x + 4 * u * n1.x + n3.x - 4 * w * n3.x +
970 4 * t * n4.x - 4 * t * n6.x + 4 * v * n7.x - 4 * u * n8.x +
971 4 * w * n8.x) *
972 ((-1 + 4 * t) * n0.y + 4 * (u * n4.y + v * n5.y + w * n6.y)) +
973 ((-1 + 4 * t) * n0.x + n1.x - 4 * u * n1.x +
974 4 * (-(t * n4.x) + u * n4.x + v * n5.x + w * n6.x - v * n7.x -
975 w * n8.x)) *
976 (4 * v * n9.y - n3.y + 4 * w * n3.y + 4 * t * n6.y + 4 * u * n8.y) -
977 ((-1 + 4 * t) * n0.x - 4 * v * n9.x + n3.x - 4 * w * n3.x + 4 * u * n4.x +
978 4 * v * n5.x - 4 * t * n6.x + 4 * w * n6.x - 4 * u * n8.x) *
979 ((-1 + 4 * u) * n1.y + 4 * (t * n4.y + v * n7.y + w * n8.y));
980
981 jac[3][0] =
982 -((((-1 + 4 * u) * n1.x + 4 * (t * n4.x + v * n7.x + w * n8.x)) *
983 (4 * w * n9.y - n2.y + 4 * v * n2.y + 4 * t * n5.y +
984 4 * u * n7.y) -
985 (4 * w * n9.x - n2.x + 4 * v * n2.x + 4 * t * n5.x + 4 * u * n7.x) *
986 ((-1 + 4 * u) * n1.y + 4 * (t * n4.y + v * n7.y + w * n8.y))) *
987 ((-1 + 4 * t) * n0.z + 4 * (u * n4.z + v * n5.z + w * n6.z))) -
988 (-(((-1 + 4 * u) * n1.x + 4 * (t * n4.x + v * n7.x + w * n8.x)) *
989 ((-1 + 4 * t) * n0.y + 4 * (u * n4.y + v * n5.y + w * n6.y))) +
990 ((-1 + 4 * t) * n0.x + 4 * (u * n4.x + v * n5.x + w * n6.x)) *
991 ((-1 + 4 * u) * n1.y + 4 * (t * n4.y + v * n7.y + w * n8.y))) *
992 (4 * w * n9.z - n2.z + 4 * v * n2.z + 4 * t * n5.z + 4 * u * n7.z) +
993 (-((4 * w * n9.x - n2.x + 4 * v * n2.x + 4 * t * n5.x + 4 * u * n7.x) *
994 ((-1 + 4 * t) * n0.y + 4 * (u * n4.y + v * n5.y + w * n6.y))) +
995 ((-1 + 4 * t) * n0.x + 4 * (u * n4.x + v * n5.x + w * n6.x)) *
996 (4 * w * n9.y - n2.y + 4 * v * n2.y + 4 * t * n5.y + 4 * u * n7.y)) *
997 ((-1 + 4 * u) * n1.z + 4 * (t * n4.z + v * n7.z + w * n8.z));
998
999 jac[3][1] =
1000 (n1.y - 4 * u * n1.y - n2.y + 4 * v * n2.y - 4 * t * n4.y + 4 * t * n5.y +
1001 4 * u * n7.y - 4 * v * n7.y + 4 * w * (n9.y - n8.y)) *
1002 ((-1 + 4 * t) * n0.z + 4 * (u * n4.z + v * n5.z + w * n6.z)) -
1003 ((-1 + 4 * t) * n0.y + n1.y - 4 * u * n1.y +
1004 4 * (-(t * n4.y) + u * n4.y + v * n5.y + w * n6.y - v * n7.y -
1005 w * n8.y)) *
1006 (4 * w * n9.z - n2.z + 4 * v * n2.z + 4 * t * n5.z + 4 * u * n7.z) +
1007 ((-1 + 4 * t) * n0.y - 4 * w * n9.y + n2.y - 4 * v * n2.y + 4 * u * n4.y -
1008 4 * t * n5.y + 4 * v * n5.y + 4 * w * n6.y - 4 * u * n7.y) *
1009 ((-1 + 4 * u) * n1.z + 4 * (t * n4.z + v * n7.z + w * n8.z));
1010
1011 jac[3][2] =
1012 (-4 * w * n9.x - n1.x + 4 * u * n1.x + n2.x - 4 * v * n2.x +
1013 4 * t * n4.x - 4 * t * n5.x - 4 * u * n7.x + 4 * v * n7.x +
1014 4 * w * n8.x) *
1015 ((-1 + 4 * t) * n0.z + 4 * (u * n4.z + v * n5.z + w * n6.z)) +
1016 ((-1 + 4 * t) * n0.x + n1.x - 4 * u * n1.x +
1017 4 * (-(t * n4.x) + u * n4.x + v * n5.x + w * n6.x - v * n7.x -
1018 w * n8.x)) *
1019 (4 * w * n9.z - n2.z + 4 * v * n2.z + 4 * t * n5.z + 4 * u * n7.z) -
1020 ((-1 + 4 * t) * n0.x - 4 * w * n9.x + n2.x - 4 * v * n2.x + 4 * u * n4.x -
1021 4 * t * n5.x + 4 * v * n5.x + 4 * w * n6.x - 4 * u * n7.x) *
1022 ((-1 + 4 * u) * n1.z + 4 * (t * n4.z + v * n7.z + w * n8.z));
1023
1024 jac[3][3] =
1025 (n1.x - 4 * u * n1.x - n2.x + 4 * v * n2.x - 4 * t * n4.x + 4 * t * n5.x +
1026 4 * u * n7.x - 4 * v * n7.x + 4 * w * (n9.x - n8.x)) *
1027 ((-1 + 4 * t) * n0.y + 4 * (u * n4.y + v * n5.y + w * n6.y)) -
1028 ((-1 + 4 * t) * n0.x + n1.x - 4 * u * n1.x +
1029 4 * (-(t * n4.x) + u * n4.x + v * n5.x + w * n6.x - v * n7.x -
1030 w * n8.x)) *
1031 (4 * w * n9.y - n2.y + 4 * v * n2.y + 4 * t * n5.y + 4 * u * n7.y) +
1032 ((-1 + 4 * t) * n0.x - 4 * w * n9.x + n2.x - 4 * v * n2.x + 4 * u * n4.x -
1033 4 * t * n5.x + 4 * v * n5.x + 4 * w * n6.x - 4 * u * n7.x) *
1034 ((-1 + 4 * u) * n1.y + 4 * (t * n4.y + v * n7.y + w * n8.y));
1035}
1036
1037void ComponentFieldMap::JacobianCube(const unsigned int i, const double t1,
1038 const double t2, const double t3,
1039 TMatrixD*& jac,
1040 std::vector<TMatrixD*>& dN) const {
1041 if (!jac) {
1042 std::cerr << m_className << "::JacobianCube:\n";
1043 std::cerr << " Pointer to Jacobian matrix is empty!\n";
1044 return;
1045 }
1046 dN.clear();
1047
1048 // Here the partial derivatives of the 8 shaping functions are calculated
1049 double N1[3] = {-1 * (1 - t2) * (1 - t3), (1 - t1) * -1 * (1 - t3),
1050 (1 - t1) * (1 - t2) * -1};
1051 double N2[3] = {+1 * (1 - t2) * (1 - t3), (1 + t1) * -1 * (1 - t3),
1052 (1 + t1) * (1 - t2) * -1};
1053 double N3[3] = {+1 * (1 + t2) * (1 - t3), (1 + t1) * +1 * (1 - t3),
1054 (1 + t1) * (1 + t2) * -1};
1055 double N4[3] = {-1 * (1 + t2) * (1 - t3), (1 - t1) * +1 * (1 - t3),
1056 (1 - t1) * (1 + t2) * -1};
1057 double N5[3] = {-1 * (1 - t2) * (1 + t3), (1 - t1) * -1 * (1 + t3),
1058 (1 - t1) * (1 - t2) * +1};
1059 double N6[3] = {+1 * (1 - t2) * (1 + t3), (1 + t1) * -1 * (1 + t3),
1060 (1 + t1) * (1 - t2) * +1};
1061 double N7[3] = {+1 * (1 + t2) * (1 + t3), (1 + t1) * +1 * (1 + t3),
1062 (1 + t1) * (1 + t2) * +1};
1063 double N8[3] = {-1 * (1 + t2) * (1 + t3), (1 - t1) * +1 * (1 + t3),
1064 (1 - t1) * (1 + t2) * +1};
1065 // Partial derivatives are stored in dN
1066 TMatrixD* m_N1 = new TMatrixD(3, 1, N1);
1067 *m_N1 = (1. / 8. * (*m_N1));
1068 dN.push_back(m_N1);
1069 TMatrixD* m_N2 = new TMatrixD(3, 1, N2);
1070 *m_N2 = (1. / 8. * (*m_N2));
1071 dN.push_back(m_N2);
1072 TMatrixD* m_N3 = new TMatrixD(3, 1, N3);
1073 *m_N3 = (1. / 8. * (*m_N3));
1074 dN.push_back(m_N3);
1075 TMatrixD* m_N4 = new TMatrixD(3, 1, N4);
1076 *m_N4 = (1. / 8. * (*m_N4));
1077 dN.push_back(m_N4);
1078 TMatrixD* m_N5 = new TMatrixD(3, 1, N5);
1079 *m_N5 = (1. / 8. * (*m_N5));
1080 dN.push_back(m_N5);
1081 TMatrixD* m_N6 = new TMatrixD(3, 1, N6);
1082 *m_N6 = (1. / 8. * (*m_N6));
1083 dN.push_back(m_N6);
1084 TMatrixD* m_N7 = new TMatrixD(3, 1, N7);
1085 *m_N7 = (1. / 8. * (*m_N7));
1086 dN.push_back(m_N7);
1087 TMatrixD* m_N8 = new TMatrixD(3, 1, N8);
1088 *m_N8 = (1. / 8. * (*m_N8));
1089 dN.push_back(m_N8);
1090 // Calculation of the jacobian using dN
1091 const Element& element = elements[i];
1092 for (int j = 0; j < 8; ++j) {
1093 const Node& node = nodes[element.emap[j]];
1094 (*jac)(0, 0) += node.x * ((*dN.at(j))(0, 0));
1095 (*jac)(0, 1) += node.y * ((*dN.at(j))(0, 0));
1096 (*jac)(0, 2) += node.z * ((*dN.at(j))(0, 0));
1097 (*jac)(1, 0) += node.x * ((*dN.at(j))(1, 0));
1098 (*jac)(1, 1) += node.y * ((*dN.at(j))(1, 0));
1099 (*jac)(1, 2) += node.z * ((*dN.at(j))(1, 0));
1100 (*jac)(2, 0) += node.x * ((*dN.at(j))(2, 0));
1101 (*jac)(2, 1) += node.y * ((*dN.at(j))(2, 0));
1102 (*jac)(2, 2) += node.z * ((*dN.at(j))(2, 0));
1103 }
1104
1105 // compute determinant
1106 if (m_debug) {
1107 std::cout << m_className << "::JacobianCube:" << std::endl;
1108 std::cout << " Det.: " << jac->Determinant() << std::endl;
1109 std::cout << " Jacobian matrix.: " << std::endl;
1110 jac->Print("%11.10g");
1111 std::cout << " Hexahedral coordinates (t, u, v) = (" << t1 << "," << t2
1112 << "," << t3 << ")" << std::endl;
1113 std::cout << " Node xyzV" << std::endl;
1114 for (int j = 0; j < 8; ++j) {
1115 const Node& node = nodes[element.emap[j]];
1116 std::cout << " " << element.emap[j] << " " << node.x
1117 << " " << node.y << " " << node.z << " "
1118 << node.v << std::endl;
1119 }
1120 }
1121}
1122
1123int ComponentFieldMap::Coordinates3(const double x, const double y,
1124 const double z, double& t1, double& t2,
1125 double& t3, double& t4, double jac[4][4],
1126 double& det,
1127 const unsigned int imap) const {
1128
1129 if (m_debug) {
1130 std::cout << m_className << "::Coordinates3:\n";
1131 std::cout << " Point (" << x << ", " << y << ", " << z << ")\n";
1132 }
1133
1134 // Failure flag
1135 int ifail = 1;
1136
1137 // Provisional values
1138 t1 = t2 = t3 = t4 = 0;
1139
1140 const Element& element = elements[imap];
1141 // Make a first order approximation, using the linear triangle.
1142 const Node& n0 = nodes[element.emap[0]];
1143 const Node& n1 = nodes[element.emap[1]];
1144 const Node& n2 = nodes[element.emap[2]];
1145 const double tt1 = (x - n1.x) * (n2.y - n1.y) - (y - n1.y) * (n2.x - n1.x);
1146 const double tt2 = (x - n2.x) * (n0.y - n2.y) - (y - n2.y) * (n0.x - n2.x);
1147 const double tt3 = (x - n0.x) * (n1.y - n0.y) - (y - n0.y) * (n1.x - n0.x);
1148 const double f1 = (n0.x - n1.x) * (n2.y - n1.y) - (n2.x - n1.x) * (n0.y - n1.y);
1149 const double f2 = (n1.x - n2.x) * (n0.y - n2.y) - (n0.x - n2.x) * (n1.y - n2.y);
1150 const double f3 = (n2.x - n0.x) * (n1.y - n0.y) - (n1.x - n0.x) * (n2.y - n0.y);
1151 if (f1 == 0 || f2 == 0 || f3 == 0) {
1152 std::cerr << m_className << "::Coordinates3:\n";
1153 std::cerr << " Calculation of linear coordinates failed; abandoned.\n";
1154 return ifail;
1155 } else {
1156 t1 = tt1 / f1;
1157 t2 = tt2 / f2;
1158 t3 = tt3 / f3;
1159 }
1160 const Node& n3 = nodes[element.emap[3]];
1161 const Node& n4 = nodes[element.emap[4]];
1162 const Node& n5 = nodes[element.emap[5]];
1163
1164 // Start iterative refinement.
1165 double td1 = t1, td2 = t2, td3 = t3;
1166 bool converged = false;
1167 for (int iter = 0; iter < 10; iter++) {
1168 if (m_debug) {
1169 std::cout << m_className << "::Coordinates3:\n";
1170 std::cout << " Iteration " << iter << ": (u, v, w) = (" << td1
1171 << ", " << td2 << ", " << td3 << "), sum = " << td1 + td2 + td3
1172 << "\n";
1173 }
1174 // Re-compute the (x,y,z) position for this coordinate.
1175 const double xr = n0.x * td1 * (2 * td1 - 1) + n1.x * td2 * (2 * td2 - 1) +
1176 n2.x * td3 * (2 * td3 - 1) + n3.x * 4 * td1 * td2 +
1177 n4.x * 4 * td1 * td3 + n5.x * 4 * td2 * td3;
1178 const double yr = n0.y * td1 * (2 * td1 - 1) + n1.y * td2 * (2 * td2 - 1) +
1179 n2.y * td3 * (2 * td3 - 1) + n3.y * 4 * td1 * td2 +
1180 n4.y * 4 * td1 * td3 + n5.y * 4 * td2 * td3;
1181 const double sr = td1 + td2 + td3;
1182 // Compute the Jacobian.
1183 Jacobian3(imap, td1, td2, td3, det, jac);
1184 // Compute the difference vector.
1185 const double diff[3] = {1 - sr, x - xr, y - yr};
1186 // Update the estimate.
1187 const double invdet = 1. / det;
1188 double corr[3] = {0., 0., 0.};
1189 for (int l = 0; l < 3; l++) {
1190 for (int k = 0; k < 3; k++) {
1191 corr[l] += jac[l][k] * diff[k] * invdet;
1192 }
1193 }
1194 // Debugging
1195 if (m_debug) {
1196 std::cout << m_className << "::Coordinates3:\n";
1197 std::cout << " Difference vector: (1, x, y) = ("
1198 << diff[0] << ", " << diff[1] << ", " << diff[2] << ").\n";
1199 std::cout << " Correction vector: (u, v, w) = ("
1200 << corr[0] << ", " << corr[1] << ", " << corr[2] << ").\n";
1201 }
1202 // Update the vector.
1203 td1 += corr[0];
1204 td2 += corr[1];
1205 td3 += corr[2];
1206 // Check for convergence.
1207 if (fabs(corr[0]) < 1.0e-5 && fabs(corr[1]) < 1.0e-5 &&
1208 fabs(corr[2]) < 1.0e-5) {
1209 if (m_debug) {
1210 std::cout << m_className << "::Coordinates3: Convergence reached.";
1211 }
1212 converged = true;
1213 break;
1214 }
1215 }
1216 // No convergence reached
1217 if (!converged) {
1218 double xmin = n0.x;
1219 double xmax = n0.x;
1220 if (n1.x < xmin) xmin = n1.x;
1221 if (n1.x > xmax) xmax = n1.x;
1222 if (n2.x < xmin) xmin = n2.x;
1223 if (n2.x > xmax) xmax = n2.x;
1224 double ymin = n0.y;
1225 double ymax = n0.y;
1226 if (n1.y < ymin) ymin = n1.y;
1227 if (n1.y > ymax) ymax = n1.y;
1228 if (n2.y < ymin) ymin = n2.y;
1229 if (n2.y > ymax) ymax = n2.y;
1230
1231 if (x >= xmin && x <= xmax && y >= ymin && y <= ymax) {
1232 std::cout << m_className << "::Coordinates3:\n";
1233 std::cout << " No convergence achieved "
1234 << "when refining internal isoparametric coordinates\n";
1235 std::cout << " in element " << imap << " at position (" << x << ", "
1236 << y << ").\n";
1237 t1 = t2 = t3 = t4 = 0;
1238 return ifail;
1239 }
1240 }
1241
1242 // Convergence reached.
1243 t1 = td1;
1244 t2 = td2;
1245 t3 = td3;
1246 t4 = 0;
1247 if (m_debug) {
1248 std::cout << m_className << "::Coordinates3:\n";
1249 std::cout << " Convergence reached at (t1, t2, t3) = (" << t1 << ", "
1250 << t2 << ", " << t3 << ").\n";
1251 }
1252
1253 // For debugging purposes, show position
1254 if (m_debug) {
1255 double xr = n0.x * td1 * (2 * td1 - 1) + n1.x * td2 * (2 * td2 - 1) +
1256 n2.x * td3 * (2 * td3 - 1) + n3.x * 4 * td1 * td2 +
1257 n4.x * 4 * td1 * td3 + n5.x * 4 * td2 * td3;
1258 double yr = n0.y * td1 * (2 * td1 - 1) + n1.y * td2 * (2 * td2 - 1) +
1259 n2.y * td3 * (2 * td3 - 1) + n3.y * 4 * td1 * td2 +
1260 n4.y * 4 * td1 * td3 + n5.y * 4 * td2 * td3;
1261 double sr = td1 + td2 + td3;
1262 std::cout << m_className << "::Coordinates3:\n";
1263 std::cout << " Position requested: (" << x << ", " << y << ")\n";
1264 std::cout << " Reconstructed: (" << xr << ", " << yr << ")\n";
1265 std::cout << " Difference: (" << x - xr << ", " << y - yr
1266 << ")\n";
1267 std::cout << " Checksum - 1: " << sr - 1 << "\n";
1268 }
1269
1270 // Success
1271 ifail = 0;
1272 return ifail;
1273}
1274
1275int ComponentFieldMap::Coordinates4(const double x, const double y,
1276 const double z, double& t1, double& t2,
1277 double& t3, double& t4, double jac[4][4],
1278 double& det,
1279 const unsigned int imap) const {
1280
1281 // Debugging
1282 if (m_debug) {
1283 std::cout << m_className << "::Coordinates4:\n";
1284 std::cout << " Point (" << x << ", " << y << ", " << z << ")\n";
1285 }
1286
1287 // Failure flag
1288 int ifail = 1;
1289
1290 // Provisional values
1291 t1 = t2 = t3 = t4 = 0.;
1292
1293 const Element& element = elements[imap];
1294 const Node& n0 = nodes[element.emap[0]];
1295 const Node& n1 = nodes[element.emap[1]];
1296 const Node& n2 = nodes[element.emap[2]];
1297 const Node& n3 = nodes[element.emap[3]];
1298 // Compute determinant.
1299 det = -(-((n0.x - n3.x) * (n1.y - n2.y)) + (n1.x - n2.x) * (n0.y - n3.y)) *
1300 (2 * x * (-n0.y + n1.y + n2.y - n3.y) -
1301 (n0.x + n3.x) * (n1.y + n2.y - 2 * y) +
1302 n1.x * (n0.y + n3.y - 2 * y) + n2.x * (n0.y + n3.y - 2 * y)) +
1303 pow(-(n0.x * n1.y) + n3.x * n2.y - n2.x * n3.y +
1304 x * (-n0.y + n1.y - n2.y + n3.y) + n1.x * (n0.y - y) +
1305 (n0.x + n2.x - n3.x) * y,
1306 2);
1307
1308 // Check that the determinant is non-negative
1309 // (this can happen if the point is out of range).
1310 if (det < 0) {
1311 if (m_debug) {
1312 std::cerr << m_className << "::Coordinates4:\n";
1313 std::cerr << " No solution found for isoparametric coordinates\n";
1314 std::cerr << " in element " << imap << " because the determinant "
1315 << det << " is < 0.\n";
1316 }
1317 return ifail;
1318 }
1319
1320 // Vector products for evaluation of T1.
1321 double prod = ((n2.x - n3.x) * (n0.y - n1.y) - (n0.x - n1.x) * (n2.y - n3.y));
1322 if (prod * prod >
1323 1.0e-12 *
1324 ((n0.x - n1.x) * (n0.x - n1.x) + (n0.y - n1.y) * (n0.y - n1.y)) *
1325 ((n2.x - n3.x) * (n2.x - n3.x) + (n2.y - n3.y) * (n2.y - n3.y))) {
1326 t1 = (-(n3.x * n0.y) + x * n0.y + n2.x * n1.y - x * n1.y - n1.x * n2.y +
1327 x * n2.y + n0.x * n3.y - x * n3.y - n0.x * y + n1.x * y - n2.x * y +
1328 n3.x * y + sqrt(det)) /
1329 prod;
1330 } else {
1331 double xp = n0.y - n1.y;
1332 double yp = n1.x - n0.x;
1333 double dn = sqrt(xp * xp + yp * yp);
1334 if (dn <= 0) {
1335 std::cerr << m_className << "::Coordinates4:\n";
1336 std::cerr << " Element " << imap
1337 << " appears to be degenerate in the 1 - 2 axis.\n";
1338 return ifail;
1339 }
1340 xp = xp / dn;
1341 yp = yp / dn;
1342 double dpoint = xp * (x - n0.x) + yp * (y - n0.y);
1343 double dbox = xp * (n3.x - n0.x) + yp * (n3.y - n0.y);
1344 if (dbox == 0) {
1345 std::cerr << m_className << "::Coordinates4:\n";
1346 std::cerr << " Element " << imap
1347 << " appears to be degenerate in the 1 - 3 axis.\n";
1348 return ifail;
1349 }
1350 double t = -1 + 2 * dpoint / dbox;
1351 double xt1 = n0.x + 0.5 * (t + 1) * (n3.x - n0.x);
1352 double yt1 = n0.y + 0.5 * (t + 1) * (n3.y - n0.y);
1353 double xt2 = n1.x + 0.5 * (t + 1) * (n2.x - n1.x);
1354 double yt2 = n1.y + 0.5 * (t + 1) * (n2.y - n1.y);
1355 dn = (xt1 - xt2) * (xt1 - xt2) + (yt1 - yt2) * (yt1 - yt2);
1356 if (dn <= 0) {
1357 std::cout << m_className << "::Coordinates4:\n";
1358 std::cout << " Coordinate requested at convergence point of element "
1359 << imap << ".\n";
1360 return ifail;
1361 }
1362 t1 = -1 + 2 * ((x - xt1) * (xt2 - xt1) + (y - yt1) * (yt2 - yt1)) / dn;
1363 }
1364
1365 // Vector products for evaluation of T2.
1366 prod = ((n0.x - n3.x) * (n1.y - n2.y) - (n1.x - n2.x) * (n0.y - n3.y));
1367 if (prod * prod >
1368 1.0e-12 *
1369 ((n0.x - n3.x) * (n0.x - n3.x) + (n0.y - n3.y) * (n0.y - n3.y)) *
1370 ((n1.x - n2.x) * (n1.x - n2.x) + (n1.y - n2.y) * (n1.y - n2.y))) {
1371 t2 = (-(n1.x * n0.y) + x * n0.y + n0.x * n1.y - x * n1.y - n3.x * n2.y +
1372 x * n2.y + n2.x * n3.y - x * n3.y - n0.x * y + n1.x * y - n2.x * y +
1373 n3.x * y - sqrt(det)) /
1374 prod;
1375 } else {
1376 double xp = n0.y - n3.y;
1377 double yp = n3.x - n0.x;
1378 double dn = sqrt(xp * xp + yp * yp);
1379 if (dn <= 0) {
1380 std::cerr << m_className << "Coordinates4:\n";
1381 std::cerr << " Element " << imap
1382 << " appears to be degenerate in the 1 - 4 axis.\n";
1383 return ifail;
1384 }
1385 xp = xp / dn;
1386 yp = yp / dn;
1387 double dpoint = xp * (x - n0.x) + yp * (y - n0.y);
1388 double dbox = xp * (n1.x - n0.x) + yp * (n1.y - n0.y);
1389 if (dbox == 0) {
1390 std::cerr << m_className << "::Coordinates4:\n";
1391 std::cerr << " Element " << imap
1392 << " appears to be degenerate in the 1 - 2 axis.\n";
1393 return ifail;
1394 }
1395 double t = -1 + 2 * dpoint / dbox;
1396 double xt1 = n0.x + 0.5 * (t + 1) * (n1.x - n0.x);
1397 double yt1 = n0.y + 0.5 * (t + 1) * (n1.y - n0.y);
1398 double xt2 = n3.x + 0.5 * (t + 1) * (n2.x - n3.x);
1399 double yt2 = n3.y + 0.5 * (t + 1) * (n2.y - n3.y);
1400 dn = (xt1 - xt2) * (xt1 - xt2) + (yt1 - yt2) * (yt1 - yt2);
1401 if (dn <= 0) {
1402 std::cout << m_className << "::Coordinates4:\n";
1403 std::cout << " Coordinate requested at convergence point of element "
1404 << imap << ".\n";
1405 return ifail;
1406 }
1407 t2 = -1 + 2 * ((x - xt1) * (xt2 - xt1) + (y - yt1) * (yt2 - yt1)) / dn;
1408 }
1409 if (m_debug) {
1410 std::cout << m_className << "::Coordinates4:\n";
1411 std::cout << " Isoparametric (u, v): (" << t1 << ", " << t2 << ").\n";
1412 }
1413
1414 // Re-compute the (x,y,z) position for this coordinate.
1415 if (m_debug) {
1416 double xr = n0.x * (1 - t1) * (1 - t2) * 0.25 +
1417 n1.x * (1 + t1) * (1 - t2) * 0.25 +
1418 n2.x * (1 + t1) * (1 + t2) * 0.25 +
1419 n3.x * (1 - t1) * (1 + t2) * 0.25;
1420 double yr = n0.y * (1 - t1) * (1 - t2) * 0.25 +
1421 n1.y * (1 + t1) * (1 - t2) * 0.25 +
1422 n2.y * (1 + t1) * (1 + t2) * 0.25 +
1423 n3.y * (1 - t1) * (1 + t2) * 0.25;
1424 std::cout << m_className << "::Coordinates4: \n";
1425 std::cout << " Position requested: (" << x << ", " << y << ")\n";
1426 std::cout << " Reconstructed: (" << xr << ", " << yr << ")\n";
1427 std::cout << " Difference: (" << x - xr << ", " << y - yr
1428 << ")\n";
1429 }
1430
1431 // This should have worked if we get this far.
1432 ifail = 0;
1433 return ifail;
1434 // Variable jac is not used.
1435 // The following lines are just for quieting the compiler.
1436 jac[0][0] = jac[0][1] = jac[0][2] = jac[0][3] = 0.;
1437 jac[1][0] = jac[1][1] = jac[1][2] = jac[1][3] = 0.;
1438 jac[2][0] = jac[2][1] = jac[2][2] = jac[2][3] = 0.;
1439 jac[3][0] = jac[3][1] = jac[3][2] = jac[3][3] = 0.;
1440}
1441
1442int ComponentFieldMap::Coordinates5(const double x, const double y,
1443 const double z, double& t1, double& t2,
1444 double& t3, double& t4, double jac[4][4],
1445 double& det,
1446 const unsigned int imap) const {
1447
1448 // Debugging
1449 if (m_debug) {
1450 std::cout << m_className << "::Coordinates5:\n";
1451 std::cout << " Point (" << x << ", " << y << ", " << z << ")\n";
1452 }
1453
1454 // Failure flag
1455 int ifail = 1;
1456
1457 // Provisional values
1458 t1 = t2 = t3 = t4 = 0;
1459
1460 const Element& element = elements[imap];
1461 // Degenerate elements should have been treated as triangles.
1462 if (element.degenerate) {
1463 std::cerr << m_className << "::Coordinates5:\n";
1464 std::cerr << " Received degenerate element " << imap << ".\n";
1465 return ifail;
1466 }
1467
1468 // Set tolerance parameter.
1469 double f = 0.5;
1470
1471 // Make a first order approximation.
1472 if (Coordinates4(x, y, z, t1, t2, t3, t4, jac, det, imap) > 0) {
1473 if (m_debug) {
1474 std::cout << m_className << "::Coordinates5:\n";
1475 std::cout << " Failure to obtain linear estimate of isoparametric "
1476 "coordinates\n";
1477 std::cout << " in element " << imap << ".\n";
1478 }
1479 return ifail;
1480 }
1481
1482 // Check whether the point is far outside.
1483 if (t1 < -(1 + f) || t1 > (1 + f) || t2 < -(1 + f) || t2 > (1 + f)) {
1484 if (m_debug) {
1485 std::cout << m_className << "::Coordinates5:\n";
1486 std::cout << " Point far outside, (t1,t2) = (" << t1 << ", " << t2
1487 << ").\n";
1488 }
1489 return ifail;
1490 }
1491
1492 // Start iteration
1493 double td1 = t1, td2 = t2;
1494 const Node& n0 = nodes[element.emap[0]];
1495 const Node& n1 = nodes[element.emap[1]];
1496 const Node& n2 = nodes[element.emap[2]];
1497 const Node& n3 = nodes[element.emap[3]];
1498 const Node& n4 = nodes[element.emap[4]];
1499 const Node& n5 = nodes[element.emap[5]];
1500 const Node& n6 = nodes[element.emap[6]];
1501 const Node& n7 = nodes[element.emap[7]];
1502 bool converged = false;
1503 for (int iter = 0; iter < 10; iter++) {
1504 if (m_debug) {
1505 std::cout << m_className << "::Coordinates5:\n";
1506 std::cout << " Iteration " << iter << ": (t1, t2) = (" << td1
1507 << ", " << td2 << ").\n";
1508 }
1509 // Re-compute the (x,y,z) position for this coordinate.
1510 double xr = n0.x * (-(1 - td1) * (1 - td2) * (1 + td1 + td2)) * 0.25 +
1511 n1.x * (-(1 + td1) * (1 - td2) * (1 - td1 + td2)) * 0.25 +
1512 n2.x * (-(1 + td1) * (1 + td2) * (1 - td1 - td2)) * 0.25 +
1513 n3.x * (-(1 - td1) * (1 + td2) * (1 + td1 - td2)) * 0.25 +
1514 n4.x * (1 - td1) * (1 + td1) * (1 - td2) * 0.5 +
1515 n5.x * (1 + td1) * (1 + td2) * (1 - td2) * 0.5 +
1516 n6.x * (1 - td1) * (1 + td1) * (1 + td2) * 0.5 +
1517 n7.x * (1 - td1) * (1 + td2) * (1 - td2) * 0.5;
1518 double yr = n0.y * (-(1 - td1) * (1 - td2) * (1 + td1 + td2)) * 0.25 +
1519 n1.y * (-(1 + td1) * (1 - td2) * (1 - td1 + td2)) * 0.25 +
1520 n2.y * (-(1 + td1) * (1 + td2) * (1 - td1 - td2)) * 0.25 +
1521 n3.y * (-(1 - td1) * (1 + td2) * (1 + td1 - td2)) * 0.25 +
1522 n4.y * (1 - td1) * (1 + td1) * (1 - td2) * 0.5 +
1523 n5.y * (1 + td1) * (1 + td2) * (1 - td2) * 0.5 +
1524 n6.y * (1 - td1) * (1 + td1) * (1 + td2) * 0.5 +
1525 n7.y * (1 - td1) * (1 + td2) * (1 - td2) * 0.5;
1526 // Compute the Jacobian.
1527 Jacobian5(imap, td1, td2, det, jac);
1528 // Compute the difference vector.
1529 double diff[2] = {x - xr, y - yr};
1530 // Update the estimate.
1531 double corr[2] = {0., 0.};
1532 for (int l = 0; l < 2; ++l) {
1533 for (int k = 0; k < 2; ++k) {
1534 corr[l] += jac[l][k] * diff[k] / det;
1535 }
1536 }
1537 // Debugging
1538 if (m_debug) {
1539 std::cout << m_className << "::Coordinates5:\n";
1540 std::cout << " Difference vector: (x, y) = (" << diff[0] << ", "
1541 << diff[1] << ").\n";
1542 std::cout << " Correction vector: (t1, t2) = (" << corr[0] << ", "
1543 << corr[1] << ").\n";
1544 }
1545 // Update the vector.
1546 td1 += corr[0];
1547 td2 += corr[1];
1548 // Check for convergence.
1549 if (fabs(corr[0]) < 1.0e-5 && fabs(corr[1]) < 1.0e-5) {
1550 if (m_debug) {
1551 std::cout << m_className << "::Coordinates5:\n";
1552 std::cout << " Convergence reached.\n";
1553 }
1554 converged = true;
1555 break;
1556 }
1557 }
1558 // No convergence reached.
1559 if (!converged) {
1560 double xmin = n0.x;
1561 double xmax = n0.x;
1562 if (n1.x < xmin) xmin = n1.x;
1563 if (n1.x > xmax) xmax = n1.x;
1564 if (n2.x < xmin) xmin = n2.x;
1565 if (n2.x > xmax) xmax = n2.x;
1566 if (n3.x < xmin) xmin = n3.x;
1567 if (n3.x > xmax) xmax = n3.x;
1568 if (n4.x < xmin) xmin = n4.x;
1569 if (n4.x > xmax) xmax = n4.x;
1570 if (n5.x < xmin) xmin = n5.x;
1571 if (n5.x > xmax) xmax = n5.x;
1572 if (n6.x < xmin) xmin = n6.x;
1573 if (n6.x > xmax) xmax = n6.x;
1574 if (n7.x < xmin) xmin = n7.x;
1575 if (n7.x > xmax) xmax = n7.x;
1576 double ymin = n0.y;
1577 double ymax = n0.y;
1578 if (n1.y < ymin) ymin = n1.y;
1579 if (n1.y > ymax) ymax = n1.y;
1580 if (n2.y < ymin) ymin = n2.y;
1581 if (n2.y > ymax) ymax = n2.y;
1582 if (n3.y < ymin) ymin = n3.y;
1583 if (n3.y > ymax) ymax = n3.y;
1584 if (n4.y < ymin) ymin = n4.y;
1585 if (n4.y > ymax) ymax = n4.y;
1586 if (n5.y < ymin) ymin = n5.y;
1587 if (n5.y > ymax) ymax = n5.y;
1588 if (n6.y < ymin) ymin = n6.y;
1589 if (n6.y > ymax) ymax = n6.y;
1590 if (n7.y < ymin) ymin = n7.y;
1591 if (n7.y > ymax) ymax = n7.y;
1592
1593 if (x >= xmin && x <= xmax && y >= ymin && y <= ymax) {
1594 std::cout << m_className << "::Coordinates5:\n";
1595 std::cout << " No convergence achieved "
1596 << "when refining internal isoparametric coordinates\n";
1597 std::cout << " in element " << imap << " at position (" << x << ", "
1598 << y << ").\n";
1599 t1 = t2 = 0;
1600 return ifail;
1601 }
1602 }
1603
1604 // Convergence reached.
1605 t1 = td1;
1606 t2 = td2;
1607 t3 = 0;
1608 t4 = 0;
1609 if (m_debug) {
1610 std::cout << m_className << "::Coordinates5:\n";
1611 std::cout << " Convergence reached at (t1, t2) = (" << t1 << ", " << t2
1612 << ").\n";
1613 }
1614
1615 // For debugging purposes, show position.
1616 if (m_debug) {
1617 double xr = n0.x * (-(1 - t1) * (1 - t2) * (1 + t1 + t2)) * 0.25 +
1618 n1.x * (-(1 + t1) * (1 - t2) * (1 - t1 + t2)) * 0.25 +
1619 n2.x * (-(1 + t1) * (1 + t2) * (1 - t1 - t2)) * 0.25 +
1620 n3.x * (-(1 - t1) * (1 + t2) * (1 + t1 - t2)) * 0.25 +
1621 n4.x * (1 - t1) * (1 + t1) * (1 - t2) * 0.5 +
1622 n5.x * (1 + t1) * (1 + t2) * (1 - t2) * 0.5 +
1623 n6.x * (1 - t1) * (1 + t1) * (1 + t2) * 0.5 +
1624 n7.x * (1 - t1) * (1 + t2) * (1 - t2) * 0.5;
1625 double yr = n0.y * (-(1 - t1) * (1 - t2) * (1 + t1 + t2)) * 0.25 +
1626 n1.y * (-(1 + t1) * (1 - t2) * (1 - t1 + t2)) * 0.25 +
1627 n2.y * (-(1 + t1) * (1 + t2) * (1 - t1 - t2)) * 0.25 +
1628 n3.y * (-(1 - t1) * (1 + t2) * (1 + t1 - t2)) * 0.25 +
1629 n4.y * (1 - t1) * (1 + t1) * (1 - t2) * 0.5 +
1630 n5.y * (1 + t1) * (1 + t2) * (1 - t2) * 0.5 +
1631 n6.y * (1 - t1) * (1 + t1) * (1 + t2) * 0.5 +
1632 n7.y * (1 - t1) * (1 + t2) * (1 - t2) * 0.5;
1633 std::cout << m_className << "::Coordinates5:\n";
1634 std::cout << " Position requested: (" << x << ", " << y << ")\n";
1635 std::cout << " Reconstructed: (" << xr << ", " << yr << ")\n";
1636 std::cout << " Difference: (" << x - xr << ", " << y - yr
1637 << ")\n";
1638 }
1639
1640 // Success
1641 ifail = 0;
1642 return ifail;
1643}
1644
1645int ComponentFieldMap::Coordinates12(const double x, const double y,
1646 const double z, double& t1, double& t2,
1647 double& t3, double& t4,
1648 const unsigned int imap) const {
1649
1650 if (m_debug) {
1651 std::cout << m_className << "::Coordinates12:\n";
1652 std::cout << " Point (" << x << ", " << y << ", " << z << ") for element "
1653 << imap << "\n";
1654 }
1655
1656 // Failure flag
1657 int ifail = 1;
1658 const Element& element = elements[imap];
1659 const Node& n0 = nodes[element.emap[0]];
1660 const Node& n1 = nodes[element.emap[1]];
1661 const Node& n2 = nodes[element.emap[2]];
1662 const Node& n3 = nodes[element.emap[3]];
1663 // Compute tetrahedral coordinates.
1664 const double f1x = (n2.y - n1.y) * (n3.z - n1.z) - (n3.y - n1.y) * (n2.z - n1.z);
1665 const double f1y = (n2.z - n1.z) * (n3.x - n1.x) - (n3.z - n1.z) * (n2.x - n1.x);
1666 const double f1z = (n2.x - n1.x) * (n3.y - n1.y) - (n3.x - n1.x) * (n2.y - n1.y);
1667 t1 = (x - n1.x) * f1x + (y - n1.y) * f1y + (z - n1.z) * f1z;
1668 t1 = t1 / ((n0.x - n1.x) * f1x + (n0.y - n1.y) * f1y + (n0.z - n1.z) * f1z);
1669 const double f2x = (n0.y - n2.y) * (n3.z - n2.z) - (n3.y - n2.y) * (n0.z - n2.z);
1670 const double f2y = (n0.z - n2.z) * (n3.x - n2.x) - (n3.z - n2.z) * (n0.x - n2.x);
1671 const double f2z = (n0.x - n2.x) * (n3.y - n2.y) - (n3.x - n2.x) * (n0.y - n2.y);
1672 t2 = (x - n2.x) * f2x + (y - n2.y) * f2y + (z - n2.z) * f2z;
1673 t2 = t2 / ((n1.x - n2.x) * f2x + (n1.y - n2.y) * f2y + (n1.z - n2.z) * f2z);
1674 const double f3x = (n0.y - n3.y) * (n1.z - n3.z) - (n1.y - n3.y) * (n0.z - n3.z);
1675 const double f3y = (n0.z - n3.z) * (n1.x - n3.x) - (n1.z - n3.z) * (n0.x - n3.x);
1676 const double f3z = (n0.x - n3.x) * (n1.y - n3.y) - (n1.x - n3.x) * (n0.y - n3.y);
1677 t3 = (x - n3.x) * f3x + (y - n3.y) * f3y + (z - n3.z) * f3z;
1678 t3 = t3 / ((n2.x - n3.x) * f3x + (n2.y - n3.y) * f3y + (n2.z - n3.z) * f3z);
1679 const double f4x = (n2.y - n0.y) * (n1.z - n0.z) - (n1.y - n0.y) * (n2.z - n0.z);
1680 const double f4y = (n2.z - n0.z) * (n1.x - n0.x) - (n1.z - n0.z) * (n2.x - n0.x);
1681 const double f4z = (n2.x - n0.x) * (n1.y - n0.y) - (n1.x - n0.x) * (n2.y - n0.y);
1682 t4 = (x - n0.x) * f4x + (y - n0.y) * f4y + (z - n0.z) * f4z;
1683 t4 = t4 / ((n3.x - n0.x) * f4x + (n3.y - n0.y) * f4y + (n3.z - n0.z) * f4z);
1684
1685 // Result
1686 if (m_debug) {
1687 std::cout << m_className << "::Coordinates12:\n";
1688 std::cout << " Tetrahedral coordinates (t, u, v, w) = (" << t1 << ", "
1689 << t2 << ", " << t3 << ", " << t4
1690 << ") sum = " << t1 + t2 + t3 + t4 << ".\n";
1691 }
1692 // Re-compute the (x,y,z) position for this coordinate.
1693 if (m_debug) {
1694 const double xr = n0.x * t1 + n1.x * t2 + n2.x * t3 + n3.x * t4;
1695 const double yr = n0.y * t1 + n1.y * t2 + n2.y * t3 + n3.y * t4;
1696 const double zr = n0.z * t1 + n1.z * t2 + n2.z * t3 + n3.z * t4;
1697 const double sr = t1 + t2 + t3 + t4;
1698 std::cout << m_className << "::Coordinates12:\n";
1699 std::cout << " Position requested: (" << x << ", " << y << ", " << z
1700 << ")\n";
1701 std::cout << " Reconstructed: (" << xr << ", " << yr << ", "
1702 << zr << ")\n";
1703 std::cout << " Difference: (" << x - xr << ", " << y - yr
1704 << ", " << z - zr << ")\n";
1705 std::cout << " Checksum - 1: " << sr - 1 << "\n";
1706 }
1707
1708 // This should always work.
1709 ifail = 0;
1710 return ifail;
1711}
1712
1713int ComponentFieldMap::Coordinates13(const double x, const double y,
1714 const double z, double& t1, double& t2,
1715 double& t3, double& t4, double jac[4][4],
1716 double& det,
1717 const unsigned int imap) const {
1718
1719 if (m_debug) {
1720 std::cout << m_className << "::Coordinates13:\n";
1721 std::cout << " Point (" << x << ", " << y << ", " << z << ")\n";
1722 }
1723
1724 // Failure flag
1725 int ifail = 1;
1726
1727 // Provisional values
1728 t1 = t2 = t3 = t4 = 0.;
1729
1730 // Make a first order approximation.
1731 if (Coordinates12(x, y, z, t1, t2, t3, t4, imap) > 0) {
1732 if (m_debug) {
1733 std::cout << m_className << "::Coordinates13:\n";
1734 std::cout << " Failure to obtain linear estimate of isoparametric "
1735 "coordinates\n";
1736 std::cout << " in element " << imap << ".\n";
1737 }
1738 return ifail;
1739 }
1740
1741 // Set tolerance parameter.
1742 const double f = 0.5;
1743 if (t1 < -f || t2 < -f || t3 < -f || t4 < -f || t1 > 1 + f || t2 > 1 + f ||
1744 t3 > 1 + f || t4 > 1 + f) {
1745 if (m_debug) {
1746 std::cout << m_className << "::Coordinates13:\n";
1747 std::cout << " Linear isoparametric coordinates more than\n";
1748 std::cout << " f (" << f << ") out of range in element " << imap
1749 << ".\n";
1750 }
1751 ifail = 0;
1752 return ifail;
1753 }
1754
1755 // Start iteration.
1756 double td1 = t1, td2 = t2, td3 = t3, td4 = t4;
1757 if (m_debug) {
1758 std::cout << m_className << "::Coordinates13:\n";
1759 std::cout << " Iteration starts at (t1,t2,t3,t4) = (" << td1 << ", "
1760 << td2 << ", " << td3 << ", " << td4 << ").\n";
1761 }
1762 const Element& element = elements[imap];
1763 const Node& n0 = nodes[element.emap[0]];
1764 const Node& n1 = nodes[element.emap[1]];
1765 const Node& n2 = nodes[element.emap[2]];
1766 const Node& n3 = nodes[element.emap[3]];
1767 const Node& n4 = nodes[element.emap[4]];
1768 const Node& n5 = nodes[element.emap[5]];
1769 const Node& n6 = nodes[element.emap[6]];
1770 const Node& n7 = nodes[element.emap[7]];
1771 const Node& n8 = nodes[element.emap[8]];
1772 const Node& n9 = nodes[element.emap[9]];
1773
1774 // Loop
1775 bool converged = false;
1776 double diff[4], corr[4];
1777 for (int iter = 0; iter < 10; iter++) {
1778 if (m_debug) {
1779 std::cout << m_className << "::Coordinates13:\n";
1780 std::cout << " Iteration " << iter << ": (t1,t2,t3,t4) = (" << td1
1781 << ", " << td2 << ", " << td3 << ", " << td4 << ").\n";
1782 }
1783 // Re-compute the (x,y,z) position for this coordinate.
1784 const double xr = n0.x * td1 * (2 * td1 - 1) + n1.x * td2 * (2 * td2 - 1) +
1785 n2.x * td3 * (2 * td3 - 1) + n3.x * td4 * (2 * td4 - 1) +
1786 n4.x * 4 * td1 * td2 + n5.x * 4 * td1 * td3 +
1787 n6.x * 4 * td1 * td4 + n7.x * 4 * td2 * td3 +
1788 n8.x * 4 * td2 * td4 + n9.x * 4 * td3 * td4;
1789 const double yr = n0.y * td1 * (2 * td1 - 1) + n1.y * td2 * (2 * td2 - 1) +
1790 n2.y * td3 * (2 * td3 - 1) + n3.y * td4 * (2 * td4 - 1) +
1791 n4.y * 4 * td1 * td2 + n5.y * 4 * td1 * td3 +
1792 n6.y * 4 * td1 * td4 + n7.y * 4 * td2 * td3 +
1793 n8.y * 4 * td2 * td4 + n9.y * 4 * td3 * td4;
1794 const double zr = n0.z * td1 * (2 * td1 - 1) + n1.z * td2 * (2 * td2 - 1) +
1795 n2.z * td3 * (2 * td3 - 1) + n3.z * td4 * (2 * td4 - 1) +
1796 n4.z * 4 * td1 * td2 + n5.z * 4 * td1 * td3 +
1797 n6.z * 4 * td1 * td4 + n7.z * 4 * td2 * td3 +
1798 n8.z * 4 * td2 * td4 + n9.z * 4 * td3 * td4;
1799 const double sr = td1 + td2 + td3 + td4;
1800
1801 // Compute the Jacobian.
1802 Jacobian13(imap, td1, td2, td3, td4, det, jac);
1803 // Compute the difference vector.
1804 diff[0] = 1 - sr;
1805 diff[1] = x - xr;
1806 diff[2] = y - yr;
1807 diff[3] = z - zr;
1808
1809 // Update the estimate.
1810 const double invdet = 1. / det;
1811 for (int l = 0; l < 4; ++l) {
1812 corr[l] = 0;
1813 for (int k = 0; k < 4; ++k) {
1814 corr[l] += jac[l][k] * diff[k] * invdet;
1815 }
1816 }
1817
1818 // Debugging
1819 if (m_debug) {
1820 std::cout << m_className << "::Coordinates13:\n";
1821 std::cout << " Difference vector: (1, x, y, z) = (" << diff[0]
1822 << ", " << diff[1] << ", " << diff[2] << ", " << diff[3]
1823 << ").\n";
1824 std::cout << " Correction vector: (t1,t2,t3,t4) = (" << corr[0]
1825 << ", " << corr[1] << ", " << corr[2] << ", " << corr[3]
1826 << ").\n";
1827 }
1828
1829 // Update the vector.
1830 td1 += corr[0];
1831 td2 += corr[1];
1832 td3 += corr[2];
1833 td4 += corr[3];
1834
1835 // Check for convergence.
1836 if (fabs(corr[0]) < 1.0e-5 && fabs(corr[1]) < 1.0e-5 &&
1837 fabs(corr[2]) < 1.0e-5 && fabs(corr[3]) < 1.0e-5) {
1838 if (m_debug) {
1839 std::cout << m_className << "::Coordinates13:\n";
1840 std::cout << " Convergence reached.\n";
1841 }
1842 converged = true;
1843 break;
1844 }
1845 }
1846
1847 // No convergence reached.
1848 if (!converged) {
1849 double xmin = n0.x;
1850 double xmax = n0.x;
1851 if (n1.x < xmin) xmin = n1.x;
1852 if (n1.x > xmax) xmax = n1.x;
1853 if (n2.x < xmin) xmin = n2.x;
1854 if (n2.x > xmax) xmax = n2.x;
1855 if (n3.x < xmin) xmin = n3.x;
1856 if (n3.x > xmax) xmax = n3.x;
1857 double ymin = n0.y;
1858 double ymax = n0.y;
1859 if (n1.y < ymin) ymin = n1.y;
1860 if (n1.y > ymax) ymax = n1.y;
1861 if (n2.y < ymin) ymin = n2.y;
1862 if (n2.y > ymax) ymax = n2.y;
1863 if (n3.y < ymin) ymin = n3.y;
1864 if (n3.y > ymax) ymax = n3.y;
1865 double zmin = n0.z;
1866 double zmax = n0.z;
1867 if (n1.z < zmin) zmin = n1.z;
1868 if (n1.z > zmax) zmax = n1.z;
1869 if (n2.z < zmin) zmin = n2.z;
1870 if (n2.z > zmax) zmax = n2.z;
1871 if (n3.z < zmin) zmin = n3.z;
1872 if (n3.z > zmax) zmax = n3.z;
1873
1874 if (x >= xmin && x <= xmax && y >= ymin && y <= ymax && z >= zmin &&
1875 z <= zmax) {
1876 std::cout << m_className << "::Coordinates13:\n";
1877 std::cout << " No convergence achieved "
1878 << "when refining internal isoparametric coordinates\n";
1879 std::cout << " in element " << imap << " at position (" << x << ", "
1880 << y << ", " << z << ").\n";
1881 t1 = t2 = t3 = t4 = -1;
1882 return ifail;
1883 }
1884 }
1885
1886 // Convergence reached.
1887 t1 = td1;
1888 t2 = td2;
1889 t3 = td3;
1890 t4 = td4;
1891 if (m_debug) {
1892 std::cout << m_className << "::Coordinates13:\n";
1893 std::cout << " Convergence reached at (t1, t2, t3, t4) = (" << t1 << ", "
1894 << t2 << ", " << t3 << ", " << t4 << ").\n";
1895 }
1896
1897 // For debugging purposes, show position.
1898 if (m_debug) {
1899 // Re-compute the (x,y,z) position for this coordinate.
1900 double xr = n0.x * td1 * (2 * td1 - 1) + n1.x * td2 * (2 * td2 - 1) +
1901 n2.x * td3 * (2 * td3 - 1) + n3.x * td4 * (2 * td4 - 1) +
1902 n4.x * 4 * td1 * td2 + n5.x * 4 * td1 * td3 +
1903 n6.x * 4 * td1 * td4 + n7.x * 4 * td2 * td3 +
1904 n8.x * 4 * td2 * td4 + n9.x * 4 * td3 * td4;
1905 double yr = n0.y * td1 * (2 * td1 - 1) + n1.y * td2 * (2 * td2 - 1) +
1906 n2.y * td3 * (2 * td3 - 1) + n3.y * td4 * (2 * td4 - 1) +
1907 n4.y * 4 * td1 * td2 + n5.y * 4 * td1 * td3 +
1908 n6.y * 4 * td1 * td4 + n7.y * 4 * td2 * td3 +
1909 n8.y * 4 * td2 * td4 + n9.y * 4 * td3 * td4;
1910 double zr = n0.z * td1 * (2 * td1 - 1) + n1.z * td2 * (2 * td2 - 1) +
1911 n2.z * td3 * (2 * td3 - 1) + n3.z * td4 * (2 * td4 - 1) +
1912 n4.z * 4 * td1 * td2 + n5.z * 4 * td1 * td3 +
1913 n6.z * 4 * td1 * td4 + n7.z * 4 * td2 * td3 +
1914 n8.z * 4 * td2 * td4 + n9.z * 4 * td3 * td4;
1915 double sr = td1 + td2 + td3 + td4;
1916 std::cout << m_className << "::Coordinates13:\n";
1917 std::cout << " Position requested: (" << x << ", " << y << ", " << z
1918 << ")\n";
1919 std::cout << " Reconstructed: (" << xr << ", " << yr << ", "
1920 << zr << ")\n";
1921 std::cout << " Difference: (" << x - xr << ", " << y - yr
1922 << ", " << z - zr << ")\n";
1923 std::cout << " Checksum - 1: " << sr - 1 << "\n";
1924 }
1925
1926 // Success
1927 ifail = 0;
1928 return ifail;
1929}
1930
1931int ComponentFieldMap::CoordinatesCube(const double x, const double y,
1932 const double z, double& t1, double& t2,
1933 double& t3, TMatrixD*& jac,
1934 std::vector<TMatrixD*>& dN,
1935 const unsigned int imap) const {
1936
1937 /*
1938 global coordinates 7__ _ _ 6 t3 t2
1939 / /| ^ /|
1940 ^ z / / | | /
1941 | 4_______5 | | /
1942 | | | | | /
1943 | | 3 | 2 |/ t1
1944 -------> | | / ------->
1945 / y | |/ local coordinates
1946 / 0--------1
1947 /
1948 v x
1949 */
1950
1951 // Failure flag
1952 int ifail = 1;
1953
1954 const Element& element = elements[imap];
1955 const Node& n0 = nodes[element.emap[0]];
1956 const Node& n2 = nodes[element.emap[2]];
1957 const Node& n3 = nodes[element.emap[3]];
1958 const Node& n7 = nodes[element.emap[7]];
1959
1960 // Compute hexahedral coordinates (t1->[-1,1],t2->[-1,1],t3->[-1,1]) and
1961 // t1 (zeta) is in y-direction
1962 // t2 (eta) is in opposite x-direction
1963 // t3 (mu) is in z-direction
1964 // Nodes are set in that way, that node [0] has always lowest x,y,z!
1965 t2 = (2. * (x - n3.x) / (n0.x - n3.x) - 1) * -1.;
1966 t1 = 2. * (y - n3.y) / (n2.y - n3.y) - 1;
1967 t3 = 2. * (z - n3.z) / (n7.z - n3.z) - 1;
1968 // Re-compute the (x,y,z) position for this coordinate.
1969 if (m_debug) {
1970 double n[8];
1971 n[0] = 1. / 8 * (1 - t1) * (1 - t2) * (1 - t3);
1972 n[1] = 1. / 8 * (1 + t1) * (1 - t2) * (1 - t3);
1973 n[2] = 1. / 8 * (1 + t1) * (1 + t2) * (1 - t3);
1974 n[3] = 1. / 8 * (1 - t1) * (1 + t2) * (1 - t3);
1975 n[4] = 1. / 8 * (1 - t1) * (1 - t2) * (1 + t3);
1976 n[5] = 1. / 8 * (1 + t1) * (1 - t2) * (1 + t3);
1977 n[6] = 1. / 8 * (1 + t1) * (1 + t2) * (1 + t3);
1978 n[7] = 1. / 8 * (1 - t1) * (1 + t2) * (1 + t3);
1979
1980 double xr = 0;
1981 double yr = 0;
1982 double zr = 0;
1983
1984 for (int i = 0; i < 8; i++) {
1985 const Node& node = nodes[element.emap[i]];
1986 xr += node.x * n[i];
1987 yr += node.y * n[i];
1988 zr += node.z * n[i];
1989 }
1990 double sr = n[0] + n[1] + n[2] + n[3] + n[4] + n[5] + n[6] + n[7];
1991 std::cout << m_className << "::CoordinatesCube:\n";
1992 std::cout << " Position requested: (" << x << "," << y << "," << z
1993 << ")\n";
1994 std::cout << " Position reconstructed: (" << xr << "," << yr << "," << zr
1995 << ")\n";
1996 std::cout << " Difference: (" << (x - xr) << "," << (y - yr)
1997 << "," << (z - zr) << ")\n";
1998 std::cout << " Hexahedral coordinates (t, u, v) = (" << t1 << "," << t2
1999 << "," << t3 << ")\n";
2000 std::cout << " Checksum - 1: " << (sr - 1) << "\n";
2001 }
2002 if (jac != 0) JacobianCube(imap, t1, t2, t3, jac, dN);
2003 // This should always work.
2004 ifail = 0;
2005 return ifail;
2006}
2007
2008void ComponentFieldMap::UpdatePeriodicityCommon() {
2009
2010 // Check the required data is available.
2011 if (!m_ready) {
2012 std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
2013 std::cerr << " No valid field map available.\n";
2014 return;
2015 }
2016
2017 // No regular and mirror periodicity at the same time.
2018 if (m_xPeriodic && m_xMirrorPeriodic) {
2019 std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
2020 std::cerr << " Both simple and mirror periodicity\n";
2021 std::cerr << " along x requested; reset.\n";
2022 m_xPeriodic = false;
2023 m_xMirrorPeriodic = false;
2024 m_warning = true;
2025 }
2026 if (m_yPeriodic && m_yMirrorPeriodic) {
2027 std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
2028 std::cerr << " Both simple and mirror periodicity\n";
2029 std::cerr << " along y requested; reset.\n";
2030 m_yPeriodic = false;
2031 m_yMirrorPeriodic = false;
2032 m_warning = true;
2033 }
2034 if (m_zPeriodic && m_zMirrorPeriodic) {
2035 std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
2036 std::cerr << " Both simple and mirror periodicity\n";
2037 std::cerr << " along z requested; reset.\n";
2038 m_zPeriodic = false;
2039 m_zMirrorPeriodic = false;
2040 m_warning = true;
2041 }
2042
2043 // In case of axial periodicity,
2044 // the range must be an integral part of two pi.
2045 if (m_xAxiallyPeriodic) {
2046 if (mapxamin >= mapxamax) {
2047 mapnxa = 0;
2048 } else {
2049 mapnxa = TwoPi / (mapxamax - mapxamin);
2050 }
2051 if (fabs(mapnxa - int(0.5 + mapnxa)) > 0.001 || mapnxa < 1.5) {
2052 std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
2053 std::cerr << " X-axial symmetry has been requested but the map\n";
2054 std::cerr << " does not cover an integral fraction of 2 pi; reset.\n";
2055 m_xAxiallyPeriodic = false;
2056 m_warning = true;
2057 }
2058 }
2059
2060 if (m_yAxiallyPeriodic) {
2061 if (mapyamin >= mapyamax) {
2062 mapnya = 0;
2063 } else {
2064 mapnya = TwoPi / (mapyamax - mapyamin);
2065 }
2066 if (fabs(mapnya - int(0.5 + mapnya)) > 0.001 || mapnya < 1.5) {
2067 std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
2068 std::cerr << " Y-axial symmetry has been requested but the map\n";
2069 std::cerr << " does not cover an integral fraction of 2 pi; reset.\n";
2070 m_yAxiallyPeriodic = false;
2071 m_warning = true;
2072 }
2073 }
2074
2075 if (m_zAxiallyPeriodic) {
2076 if (mapzamin >= mapzamax) {
2077 mapnza = 0;
2078 } else {
2079 mapnza = TwoPi / (mapzamax - mapzamin);
2080 }
2081 if (fabs(mapnza - int(0.5 + mapnza)) > 0.001 || mapnza < 1.5) {
2082 std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
2083 std::cerr << " Z-axial symmetry has been requested but the map\n";
2084 std::cerr << " does not cover an integral fraction of 2 pi; reset.\n";
2085 m_zAxiallyPeriodic = false;
2086 m_warning = true;
2087 }
2088 }
2089
2090 // Not more than 1 rotational symmetry
2091 if ((m_xRotationSymmetry && m_yRotationSymmetry) ||
2092 (m_xRotationSymmetry && m_zRotationSymmetry) ||
2093 (m_yRotationSymmetry && m_zRotationSymmetry)) {
2094 std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
2095 std::cerr << " Only 1 rotational symmetry allowed; reset.\n";
2096 m_xRotationSymmetry = false;
2097 m_yRotationSymmetry = false;
2098 m_zRotationSymmetry = false;
2099 m_warning = true;
2100 }
2101
2102 // No rotational symmetry as well as axial periodicity
2103 if ((m_xRotationSymmetry || m_yRotationSymmetry || m_zRotationSymmetry) &&
2104 (m_xAxiallyPeriodic || m_yAxiallyPeriodic || m_zAxiallyPeriodic)) {
2105 std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
2106 std::cerr << " Not allowed to combine rotational symmetry\n";
2107 std::cerr << " and axial periodicity; reset.\n";
2108 m_xAxiallyPeriodic = false;
2109 m_yAxiallyPeriodic = false;
2110 m_zAxiallyPeriodic = false;
2111 m_xRotationSymmetry = false;
2112 m_yRotationSymmetry = false;
2113 m_zRotationSymmetry = false;
2114 m_warning = true;
2115 }
2116
2117 // In case of rotational symmetry, the x-range should not straddle 0.
2118 if (m_xRotationSymmetry || m_yRotationSymmetry || m_zRotationSymmetry) {
2119 if (mapxmin * mapxmax < 0) {
2120 std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
2121 std::cerr << " Rotational symmetry requested, \n";
2122 std::cerr << " but x-range straddles 0; reset.\n";
2123 m_xRotationSymmetry = false;
2124 m_yRotationSymmetry = false;
2125 m_zRotationSymmetry = false;
2126 m_warning = true;
2127 }
2128 }
2129
2130 // Recompute the cell ranges.
2131 xMinBoundingBox = mapxmin;
2132 xMaxBoundingBox = mapxmax;
2133 yMinBoundingBox = mapymin;
2134 yMaxBoundingBox = mapymax;
2135 zMinBoundingBox = mapzmin;
2136 zMaxBoundingBox = mapzmax;
2137 cellsx = fabs(mapxmax - mapxmin);
2138 cellsy = fabs(mapymax - mapymin);
2139 cellsz = fabs(mapzmax - mapzmin);
2140 if (m_xRotationSymmetry) {
2141 xMinBoundingBox = mapymin;
2142 xMaxBoundingBox = mapymax;
2143 yMinBoundingBox = -std::max(fabs(mapxmin), fabs(mapxmax));
2144 yMaxBoundingBox = +std::max(fabs(mapxmin), fabs(mapxmax));
2145 zMinBoundingBox = -std::max(fabs(mapxmin), fabs(mapxmax));
2146 zMaxBoundingBox = +std::max(fabs(mapxmin), fabs(mapxmax));
2147 } else if (m_yRotationSymmetry) {
2148 xMinBoundingBox = -std::max(fabs(mapxmin), fabs(mapxmax));
2149 xMaxBoundingBox = +std::max(fabs(mapxmin), fabs(mapxmax));
2150 yMinBoundingBox = mapymin;
2151 yMaxBoundingBox = mapymax;
2152 zMinBoundingBox = -std::max(fabs(mapxmin), fabs(mapxmax));
2153 zMaxBoundingBox = +std::max(fabs(mapxmin), fabs(mapxmax));
2154 } else if (m_zRotationSymmetry) {
2155 xMinBoundingBox = -std::max(fabs(mapxmin), fabs(mapxmax));
2156 xMaxBoundingBox = +std::max(fabs(mapxmin), fabs(mapxmax));
2157 yMinBoundingBox = -std::max(fabs(mapxmin), fabs(mapxmax));
2158 yMaxBoundingBox = +std::max(fabs(mapxmin), fabs(mapxmax));
2159 zMinBoundingBox = mapymin;
2160 zMaxBoundingBox = mapymax;
2161 }
2162
2163 if (m_xAxiallyPeriodic) {
2164 yMinBoundingBox = -std::max(std::max(fabs(mapymin), fabs(mapymax)),
2165 std::max(fabs(mapzmin), fabs(mapzmax)));
2166 yMaxBoundingBox = +std::max(std::max(fabs(mapymin), fabs(mapymax)),
2167 std::max(fabs(mapzmin), fabs(mapzmax)));
2168 zMinBoundingBox = -std::max(std::max(fabs(mapymin), fabs(mapymax)),
2169 std::max(fabs(mapzmin), fabs(mapzmax)));
2170 zMaxBoundingBox = +std::max(std::max(fabs(mapymin), fabs(mapymax)),
2171 std::max(fabs(mapzmin), fabs(mapzmax)));
2172 } else if (m_yAxiallyPeriodic) {
2173 xMinBoundingBox = -std::max(std::max(fabs(mapxmin), fabs(mapxmax)),
2174 std::max(fabs(mapzmin), fabs(mapzmax)));
2175 xMaxBoundingBox = +std::max(std::max(fabs(mapxmin), fabs(mapxmax)),
2176 std::max(fabs(mapzmin), fabs(mapzmax)));
2177 zMinBoundingBox = -std::max(std::max(fabs(mapxmin), fabs(mapxmax)),
2178 std::max(fabs(mapzmin), fabs(mapzmax)));
2179 zMaxBoundingBox = +std::max(std::max(fabs(mapxmin), fabs(mapxmax)),
2180 std::max(fabs(mapzmin), fabs(mapzmax)));
2181 } else if (m_zAxiallyPeriodic) {
2182 xMinBoundingBox = -std::max(std::max(fabs(mapxmin), fabs(mapxmax)),
2183 std::max(fabs(mapymin), fabs(mapymax)));
2184 xMaxBoundingBox = +std::max(std::max(fabs(mapxmin), fabs(mapxmax)),
2185 std::max(fabs(mapymin), fabs(mapymax)));
2186 yMinBoundingBox = -std::max(std::max(fabs(mapxmin), fabs(mapxmax)),
2187 std::max(fabs(mapymin), fabs(mapymax)));
2188 yMaxBoundingBox = +std::max(std::max(fabs(mapxmin), fabs(mapxmax)),
2189 std::max(fabs(mapymin), fabs(mapymax)));
2190 }
2191
2192 if (m_xPeriodic || m_xMirrorPeriodic) {
2193 xMinBoundingBox = -INFINITY;
2194 xMaxBoundingBox = +INFINITY;
2195 }
2196 if (m_yPeriodic || m_yMirrorPeriodic) {
2197 yMinBoundingBox = -INFINITY;
2198 yMaxBoundingBox = +INFINITY;
2199 }
2200 if (m_zPeriodic || m_zMirrorPeriodic) {
2201 zMinBoundingBox = -INFINITY;
2202 zMaxBoundingBox = +INFINITY;
2203 }
2204
2205 // Display the range if requested.
2206 if (m_debug) PrintRange();
2207}
2208
2209void ComponentFieldMap::UpdatePeriodicity2d() {
2210
2211 // Check the required data is available.
2212 if (!m_ready) {
2213 std::cerr << m_className << "::UpdatePeriodicity2d:\n";
2214 std::cerr << " No valid field map available.\n";
2215 return;
2216 }
2217
2218 // No z-periodicity in 2d
2219 if (m_zPeriodic || m_zMirrorPeriodic) {
2220 std::cerr << m_className << "::UpdatePeriodicity2d:\n";
2221 std::cerr << " Simple or mirror periodicity along z\n";
2222 std::cerr << " requested for a 2d map; reset.\n";
2223 m_zPeriodic = false;
2224 m_zMirrorPeriodic = false;
2225 m_warning = true;
2226 }
2227
2228 // Only z-axial periodicity in 2d maps
2229 if (m_xAxiallyPeriodic || m_yAxiallyPeriodic) {
2230 std::cerr << m_className << "::UpdatePeriodicity2d:\n";
2231 std::cerr << " Axial symmetry has been requested \n";
2232 std::cerr << " around x or y for a 2D map; reset.\n";
2233 m_xAxiallyPeriodic = false;
2234 m_yAxiallyPeriodic = false;
2235 m_warning = true;
2236 }
2237}
2238
2239void ComponentFieldMap::SetRange() {
2240
2241 // Initial values
2242 mapxmin = mapymin = mapzmin = 0.;
2243 mapxmax = mapymax = mapzmax = 0.;
2244 mapxamin = mapyamin = mapzamin = 0.;
2245 mapxamax = mapyamax = mapzamax = 0.;
2246 mapvmin = mapvmax = 0.;
2247 setangx = setangy = setangz = false;
2248
2249 // Make sure the required data is available.
2250 if (!m_ready || nNodes < 1) {
2251 std::cerr << m_className << "::SetRange:\n";
2252 std::cerr << " Field map not yet set.\n";
2253 return;
2254 }
2255 if (nNodes < 1) {
2256 std::cerr << m_className << "::SetRange:\n";
2257 std::cerr << " Number of nodes < 1.\n";
2258 return;
2259 }
2260
2261 // Loop over the nodes.
2262 mapxmin = mapxmax = nodes[0].x;
2263 mapymin = mapymax = nodes[0].y;
2264 mapzmin = mapzmax = nodes[0].z;
2265 mapvmin = mapvmax = nodes[0].v;
2266
2267 double ang;
2268 for (int i = 1; i < nNodes; i++) {
2269 if (mapxmin > nodes[i].x) mapxmin = nodes[i].x;
2270 if (mapxmax < nodes[i].x) mapxmax = nodes[i].x;
2271 if (mapymin > nodes[i].y) mapymin = nodes[i].y;
2272 if (mapymax < nodes[i].y) mapymax = nodes[i].y;
2273 if (mapzmin > nodes[i].z) mapzmin = nodes[i].z;
2274 if (mapzmax < nodes[i].z) mapzmax = nodes[i].z;
2275 if (mapvmin > nodes[i].v) mapvmin = nodes[i].v;
2276 if (mapvmax < nodes[i].v) mapvmax = nodes[i].v;
2277
2278 if (nodes[i].y != 0 || nodes[i].z != 0) {
2279 ang = atan2(nodes[i].z, nodes[i].y);
2280 if (setangx) {
2281 if (ang < mapxamin) mapxamin = ang;
2282 if (ang > mapxamax) mapxamax = ang;
2283 } else {
2284 mapxamin = mapxamax = ang;
2285 setangx = true;
2286 }
2287 }
2288
2289 if (nodes[i].z != 0 || nodes[i].x != 0) {
2290 ang = atan2(nodes[i].x, nodes[i].z);
2291 if (setangy) {
2292 if (ang < mapyamin) mapyamin = ang;
2293 if (ang > mapyamax) mapyamax = ang;
2294 } else {
2295 mapyamin = mapyamax = ang;
2296 setangy = true;
2297 }
2298 }
2299
2300 if (nodes[i].x != 0 || nodes[i].y != 0) {
2301 ang = atan2(nodes[i].y, nodes[i].x);
2302 if (setangz) {
2303 if (ang < mapzamin) mapzamin = ang;
2304 if (ang > mapzamax) mapzamax = ang;
2305 } else {
2306 mapzamin = mapzamax = ang;
2307 setangz = true;
2308 }
2309 }
2310 }
2311
2312 // Fix the angular ranges.
2313 if (mapxamax - mapxamin > Pi) {
2314 double aux = mapxamin;
2315 mapxamin = mapxamax;
2316 mapxamax = aux + TwoPi;
2317 }
2318
2319 if (mapyamax - mapyamin > Pi) {
2320 double aux = mapyamin;
2321 mapyamin = mapyamax;
2322 mapyamax = aux + TwoPi;
2323 }
2324
2325 if (mapzamax - mapzamin > Pi) {
2326 double aux = mapzamin;
2327 mapzamin = mapzamax;
2328 mapzamax = aux + TwoPi;
2329 }
2330
2331 // Set the periodicity length (maybe not needed).
2332 mapsx = fabs(mapxmax - mapxmin);
2333 mapsy = fabs(mapymax - mapymin);
2334 mapsz = fabs(mapzmax - mapzmin);
2335
2336 // Set provisional cell dimensions.
2337 xMinBoundingBox = mapxmin;
2338 xMaxBoundingBox = mapxmax;
2339 yMinBoundingBox = mapymin;
2340 yMaxBoundingBox = mapymax;
2341 if (m_is3d) {
2342 zMinBoundingBox = mapzmin;
2343 zMaxBoundingBox = mapzmax;
2344 } else {
2345 mapzmin = zMinBoundingBox;
2346 mapzmax = zMaxBoundingBox;
2347 }
2348 hasBoundingBox = true;
2349
2350 // Display the range if requested.
2351 if (m_debug) PrintRange();
2352}
2353
2354void ComponentFieldMap::PrintRange() {
2355
2356 std::cout << m_className << "::PrintRange:\n";
2357 std::cout << " Dimensions of the elementary block\n";
2358 printf(" %15g < x < %-15g cm,\n", mapxmin, mapxmax);
2359 printf(" %15g < y < %-15g cm,\n", mapymin, mapymax);
2360 printf(" %15g < z < %-15g cm,\n", mapzmin, mapzmax);
2361 printf(" %15g < V < %-15g V.\n", mapvmin, mapvmax);
2362
2363 std::cout << " Periodicities\n";
2364
2365 std::cout << " x:";
2366 if (m_xPeriodic) {
2367 std::cout << " simple with length " << cellsx << " cm";
2368 }
2369 if (m_xMirrorPeriodic) {
2370 std::cout << " mirror with length " << cellsx << " cm";
2371 }
2372 if (m_xAxiallyPeriodic) {
2373 std::cout << " axial " << int(0.5 + mapnxa) << "-fold repetition";
2374 }
2375 if (m_xRotationSymmetry) std::cout << " rotational symmetry";
2376 if (!(m_xPeriodic || m_xMirrorPeriodic || m_xAxiallyPeriodic ||
2377 m_xRotationSymmetry))
2378 std::cout << " none";
2379 std::cout << "\n";
2380
2381 std::cout << " y:";
2382 if (m_yPeriodic) {
2383 std::cout << " simple with length " << cellsy << " cm";
2384 }
2385 if (m_yMirrorPeriodic) {
2386 std::cout << " mirror with length " << cellsy << " cm";
2387 }
2388 if (m_yAxiallyPeriodic) {
2389 std::cout << " axial " << int(0.5 + mapnya) << "-fold repetition";
2390 }
2391 if (m_yRotationSymmetry) {
2392 std::cout << " rotational symmetry";
2393 }
2394 if (!(m_yPeriodic || m_yMirrorPeriodic || m_yAxiallyPeriodic ||
2395 m_yRotationSymmetry))
2396 std::cout << " none";
2397 std::cout << "\n";
2398
2399 std::cout << " z:";
2400 if (m_zPeriodic) {
2401 std::cout << " simple with length " << cellsz << " cm";
2402 }
2403 if (m_zMirrorPeriodic) {
2404 std::cout << " mirror with length " << cellsz << " cm";
2405 }
2406 if (m_zAxiallyPeriodic) {
2407 std::cout << " axial " << int(0.5 + mapnza) << "-fold repetition";
2408 }
2409 if (m_zRotationSymmetry) {
2410 std::cout << " rotational symmetry";
2411 }
2412 if (!(m_zPeriodic || m_zMirrorPeriodic || m_zAxiallyPeriodic ||
2413 m_zRotationSymmetry))
2414 std::cout << " none";
2415 std::cout << "\n";
2416}
2417
2418bool ComponentFieldMap::IsInBoundingBox(const double x, const double y,
2419 const double z) const {
2420
2421 return (x >= xMinBoundingBox && x <= xMaxBoundingBox &&
2422 y >= yMinBoundingBox && y <= yMaxBoundingBox &&
2423 z >= zMinBoundingBox && z <= zMaxBoundingBox);
2424}
2425
2426bool ComponentFieldMap::GetBoundingBox(double& xmin, double& ymin, double& zmin,
2427 double& xmax, double& ymax,
2428 double& zmax) {
2429
2430 if (!m_ready) return false;
2431
2432 xmin = xMinBoundingBox;
2433 xmax = xMaxBoundingBox;
2434 ymin = yMinBoundingBox;
2435 ymax = yMaxBoundingBox;
2436 zmin = zMinBoundingBox;
2437 zmax = zMaxBoundingBox;
2438 return true;
2439}
2440
2441void ComponentFieldMap::MapCoordinates(double& xpos, double& ypos, double& zpos,
2442 bool& xmirrored, bool& ymirrored,
2443 bool& zmirrored, double& rcoordinate,
2444 double& rotation) const {
2445
2446 // Initial values
2447 rotation = 0;
2448
2449 // If chamber is periodic, reduce to the cell volume.
2450 xmirrored = false;
2451 double auxr, auxphi;
2452 if (m_xPeriodic) {
2453 xpos = mapxmin + fmod(xpos - mapxmin, mapxmax - mapxmin);
2454 if (xpos < mapxmin) xpos += mapxmax - mapxmin;
2455 } else if (m_xMirrorPeriodic) {
2456 double xnew = mapxmin + fmod(xpos - mapxmin, mapxmax - mapxmin);
2457 if (xnew < mapxmin) xnew += mapxmax - mapxmin;
2458 int nx = int(floor(0.5 + (xnew - xpos) / (mapxmax - mapxmin)));
2459 if (nx != 2 * (nx / 2)) {
2460 xnew = mapxmin + mapxmax - xnew;
2461 xmirrored = true;
2462 }
2463 xpos = xnew;
2464 }
2465 if (m_xAxiallyPeriodic && (zpos != 0 || ypos != 0)) {
2466 auxr = sqrt(zpos * zpos + ypos * ypos);
2467 auxphi = atan2(zpos, ypos);
2468 rotation = (mapxamax - mapxamin) *
2469 floor(0.5 + (auxphi - 0.5 * (mapxamin + mapxamax)) /
2470 (mapxamax - mapxamin));
2471 if (auxphi - rotation < mapxamin)
2472 rotation = rotation - (mapxamax - mapxamin);
2473 if (auxphi - rotation > mapxamax)
2474 rotation = rotation + (mapxamax - mapxamin);
2475 auxphi = auxphi - rotation;
2476 ypos = auxr * cos(auxphi);
2477 zpos = auxr * sin(auxphi);
2478 }
2479
2480 ymirrored = false;
2481 if (m_yPeriodic) {
2482 ypos = mapymin + fmod(ypos - mapymin, mapymax - mapymin);
2483 if (ypos < mapymin) ypos += mapymax - mapymin;
2484 } else if (m_yMirrorPeriodic) {
2485 double ynew = mapymin + fmod(ypos - mapymin, mapymax - mapymin);
2486 if (ynew < mapymin) ynew += mapymax - mapymin;
2487 int ny = int(floor(0.5 + (ynew - ypos) / (mapymax - mapymin)));
2488 if (ny != 2 * (ny / 2)) {
2489 ynew = mapymin + mapymax - ynew;
2490 ymirrored = true;
2491 }
2492 ypos = ynew;
2493 }
2494 if (m_yAxiallyPeriodic && (xpos != 0 || zpos != 0)) {
2495 auxr = sqrt(xpos * xpos + zpos * zpos);
2496 auxphi = atan2(xpos, zpos);
2497 rotation = (mapyamax - mapyamin) *
2498 floor(0.5 + (auxphi - 0.5 * (mapyamin + mapyamax)) /
2499 (mapyamax - mapyamin));
2500 if (auxphi - rotation < mapyamin)
2501 rotation = rotation - (mapyamax - mapyamin);
2502 if (auxphi - rotation > mapyamax)
2503 rotation = rotation + (mapyamax - mapyamin);
2504 auxphi = auxphi - rotation;
2505 zpos = auxr * cos(auxphi);
2506 xpos = auxr * sin(auxphi);
2507 }
2508
2509 zmirrored = false;
2510 if (m_zPeriodic) {
2511 zpos = mapzmin + fmod(zpos - mapzmin, mapzmax - mapzmin);
2512 if (zpos < mapzmin) zpos += mapzmax - mapzmin;
2513 } else if (m_zMirrorPeriodic) {
2514 double znew = mapzmin + fmod(zpos - mapzmin, mapzmax - mapzmin);
2515 if (znew < mapzmin) znew += mapzmax - mapzmin;
2516 int nz = int(floor(0.5 + (znew - zpos) / (mapzmax - mapzmin)));
2517 if (nz != 2 * (nz / 2)) {
2518 znew = mapzmin + mapzmax - znew;
2519 zmirrored = true;
2520 }
2521 zpos = znew;
2522 }
2523 if (m_zAxiallyPeriodic && (ypos != 0 || xpos != 0)) {
2524 auxr = sqrt(ypos * ypos + xpos * xpos);
2525 auxphi = atan2(ypos, xpos);
2526 rotation = (mapzamax - mapzamin) *
2527 floor(0.5 + (auxphi - 0.5 * (mapzamin + mapzamax)) /
2528 (mapzamax - mapzamin));
2529 if (auxphi - rotation < mapzamin)
2530 rotation = rotation - (mapzamax - mapzamin);
2531 if (auxphi - rotation > mapzamax)
2532 rotation = rotation + (mapzamax - mapzamin);
2533 auxphi = auxphi - rotation;
2534 xpos = auxr * cos(auxphi);
2535 ypos = auxr * sin(auxphi);
2536 }
2537
2538 // If we have a rotationally symmetric field map, store coordinates.
2539 rcoordinate = 0;
2540 double zcoordinate = 0;
2541 if (m_xRotationSymmetry) {
2542 rcoordinate = sqrt(ypos * ypos + zpos * zpos);
2543 zcoordinate = xpos;
2544 } else if (m_yRotationSymmetry) {
2545 rcoordinate = sqrt(xpos * xpos + zpos * zpos);
2546 zcoordinate = ypos;
2547 } else if (m_zRotationSymmetry) {
2548 rcoordinate = sqrt(xpos * xpos + ypos * ypos);
2549 zcoordinate = zpos;
2550 }
2551
2552 if (m_xRotationSymmetry || m_yRotationSymmetry || m_zRotationSymmetry) {
2553 xpos = rcoordinate;
2554 ypos = zcoordinate;
2555 zpos = 0;
2556 }
2557}
2558
2559void ComponentFieldMap::UnmapFields(double& ex, double& ey, double& ez,
2560 double& xpos, double& ypos, double& zpos,
2561 bool& xmirrored, bool& ymirrored,
2562 bool& zmirrored, double& rcoordinate,
2563 double& rotation) const {
2564
2565 // Apply mirror imaging.
2566 if (xmirrored) ex = -ex;
2567 if (ymirrored) ey = -ey;
2568 if (zmirrored) ez = -ez;
2569
2570 // Rotate the field.
2571 double er, theta;
2572 if (m_xAxiallyPeriodic) {
2573 er = sqrt(ey * ey + ez * ez);
2574 theta = atan2(ez, ey);
2575 theta += rotation;
2576 ey = er * cos(theta);
2577 ez = er * sin(theta);
2578 }
2579 if (m_yAxiallyPeriodic) {
2580 er = sqrt(ez * ez + ex * ex);
2581 theta = atan2(ex, ez);
2582 theta += rotation;
2583 ez = er * cos(theta);
2584 ex = er * sin(theta);
2585 }
2586 if (m_zAxiallyPeriodic) {
2587 er = sqrt(ex * ex + ey * ey);
2588 theta = atan2(ey, ex);
2589 theta += rotation;
2590 ex = er * cos(theta);
2591 ey = er * sin(theta);
2592 }
2593
2594 // Take care of symmetry.
2595 double eaxis;
2596 er = ex;
2597 eaxis = ey;
2598
2599 // Rotational symmetry
2600 if (m_xRotationSymmetry) {
2601 if (rcoordinate <= 0) {
2602 ex = eaxis;
2603 ey = 0;
2604 ez = 0;
2605 } else {
2606 ex = eaxis;
2607 ey = er * ypos / rcoordinate;
2608 ez = er * zpos / rcoordinate;
2609 }
2610 }
2611 if (m_yRotationSymmetry) {
2612 if (rcoordinate <= 0) {
2613 ex = 0;
2614 ey = eaxis;
2615 ez = 0;
2616 } else {
2617 ex = er * xpos / rcoordinate;
2618 ey = eaxis;
2619 ez = er * zpos / rcoordinate;
2620 }
2621 }
2622 if (m_zRotationSymmetry) {
2623 if (rcoordinate <= 0) {
2624 ex = 0;
2625 ey = 0;
2626 ez = eaxis;
2627 } else {
2628 ex = er * xpos / rcoordinate;
2629 ey = er * ypos / rcoordinate;
2630 ez = eaxis;
2631 }
2632 }
2633}
2634
2635int ComponentFieldMap::ReadInteger(char* token, int def, bool& error) {
2636
2637 if (!token) {
2638 error = true;
2639 return def;
2640 }
2641
2642 return atoi(token);
2643}
2644
2645double ComponentFieldMap::ReadDouble(char* token, double def, bool& error) {
2646
2647 if (!token) {
2648 error = true;
2649 return def;
2650 }
2651 return atof(token);
2652}
2653
2654void ComponentFieldMap::CalculateElementBoundingBoxes(void) {
2655
2656 // Do not proceed if not properly initialised.
2657 if (!m_ready) {
2658 PrintNotReady("CalculateElementBoundingBoxes");
2659 return;
2660 }
2661
2662 // Calculate the bounding boxes of all elements
2663 for (int i = 0; i < nElements; ++i) {
2664 Element& elem = elements[i];
2665 const Node& n0 = nodes[elem.emap[0]];
2666 const Node& n1 = nodes[elem.emap[1]];
2667 const Node& n2 = nodes[elem.emap[2]];
2668 const Node& n3 = nodes[elem.emap[3]];
2669 elem.xmin = std::min(std::min(n0.x, n1.x), std::min(n2.x, n3.x));
2670 elem.xmax = std::max(std::max(n0.x, n1.x), std::max(n2.x, n3.x));
2671 elem.ymin = std::min(std::min(n0.y, n1.y), std::min(n2.y, n3.y));
2672 elem.ymax = std::max(std::max(n0.y, n1.y), std::max(n2.y, n3.y));
2673 elem.zmin = std::min(std::min(n0.z, n1.z), std::min(n2.z, n3.z));
2674 elem.zmax = std::max(std::max(n0.z, n1.z), std::max(n2.z, n3.z));
2675 }
2676}
2677
2678bool ComponentFieldMap::InitializeTetrahedralTree() {
2679
2680 // Do not proceed if not properly initialised.
2681 if (!m_ready) {
2682 PrintNotReady("InitializeTetrahedralTree");
2683 return false;
2684 }
2685
2686 std::cout << m_className << "::InitializeTetrahedralTree:\n"
2687 << " About to initialize the tetrahedral tree.\n";
2688
2689 // check if the caching has not been done before
2690 if (!m_cacheElemBoundingBoxes) CalculateElementBoundingBoxes();
2691
2692 // Determine the bounding box
2693 double xmin = 0., ymin = 0., zmin = 0., xmax = 0., ymax = 0., zmax = 0.;
2694 for (unsigned int i = 0; i < nodes.size(); i++) {
2695 const Node& n = nodes[i];
2696 if (n.x <= xmin) xmin = n.x;
2697 if (n.x > xmax) xmax = n.x;
2698 if (n.y <= ymin) ymin = n.y;
2699 if (n.y > ymax) ymax = n.y;
2700 if (n.z <= zmin) zmin = n.z;
2701 if (n.z > zmax) zmax = n.z;
2702 }
2703
2704 std::cout << " Bounding box:\n"
2705 << std::scientific << "\tx: " << xmin << " -> " << xmax << "\n"
2706 << std::scientific << "\ty: " << ymin << " -> " << ymax << "\n"
2707 << std::scientific << "\tz: " << zmin << " -> " << zmax << "\n";
2708
2709 const double hx = 0.5 * (xmax - xmin);
2710 const double hy = 0.5 * (ymax - ymin);
2711 const double hz = 0.5 * (zmax - zmin);
2712 m_tetTree = new TetrahedralTree(Vec3(xmin + hx, ymin + hy, zmin + hz),
2713 Vec3(hx, hy, hz));
2714
2715 std::cout << " Tree instantiated.\n";
2716
2717 // insert all mesh nodes in the tree
2718 for (unsigned int i = 0; i < nodes.size(); i++) {
2719 const Node& n = nodes[i];
2720 m_tetTree->InsertMeshNode(Vec3(n.x, n.y, n.z), i);
2721 }
2722
2723 std::cout << m_className << "::InitializeTetrahedralTree:\n"
2724 << " Tetrahedral tree nodes initialized successfully.\n";
2725
2726 // insert all mesh elements (tetrahedrons) in the tree
2727 for (unsigned int i = 0; i < elements.size(); i++) {
2728 const Element& e = elements[i];
2729 const double bb[6] = {e.xmin, e.ymin, e.zmin, e.xmax, e.ymax, e.zmax};
2730 m_tetTree->InsertTetrahedron(bb, i);
2731 }
2732
2733 std::cerr << m_className << "::InitializeTetrahedralTree:\n";
2734 std::cerr << " Tetrahedral tree initialized successfully.\n";
2735
2736 m_isTreeInitialized = true;
2737 return true;
2738}
2739
2740void ComponentFieldMap::PrintElement(const std::string& header, const double x,
2741 const double y, const double z,
2742 const double t1, const double t2,
2743 const double t3, const double t4,
2744 const unsigned int i, const unsigned int n,
2745 const int iw) const {
2746
2747 const Element& element = elements[i];
2748 std::cout << m_className << "::" << header << ":\n"
2749 << " Global = (" << x << ", " << y << ", " << z << ")\n"
2750 << " Local = (" << t1 << ", " << t2 << ", " << t3 << ", " << t4
2751 << ")\n"
2752 << " Element = " << i << " (degenerate: " << element.degenerate
2753 << ")\n "
2754 << " Node x y z V\n";
2755 for (unsigned int ii = 0; ii < n; ++ii) {
2756 const Node& node = nodes[element.emap[ii]];
2757 const double v = iw < 0 ? node.v : node.w[iw];
2758 printf(" %-5d %12g %12g %12g %12g\n", element.emap[i], node.x, node.y,
2759 node.z, v);
2760 }
2761}
2762}
Garfield::ComponentBase
Abstract base class for components.
Definition: ComponentBase.hh:13
Garfield::ComponentBase::m_zMirrorPeriodic
bool m_zMirrorPeriodic
Mirror periodicity in z.
Definition: ComponentBase.hh:257
Garfield::ComponentBase::m_yAxiallyPeriodic
bool m_yAxiallyPeriodic
Axial periodicity in y.
Definition: ComponentBase.hh:261
Garfield::ComponentBase::m_zRotationSymmetry
bool m_zRotationSymmetry
Rotation symmetry around z-axis.
Definition: ComponentBase.hh:269
Garfield::ComponentBase::m_yRotationSymmetry
bool m_yRotationSymmetry
Rotation symmetry around y-axis.
Definition: ComponentBase.hh:267
Garfield::ComponentBase::m_ready
bool m_ready
Ready for use?
Definition: ComponentBase.hh:239
Garfield::ComponentBase::m_zAxiallyPeriodic
bool m_zAxiallyPeriodic
Axial periodicity in z.
Definition: ComponentBase.hh:263
Garfield::ComponentBase::m_xRotationSymmetry
bool m_xRotationSymmetry
Rotation symmetry around x-axis.
Definition: ComponentBase.hh:265
Garfield::ComponentBase::m_yPeriodic
bool m_yPeriodic
Simple periodicity in y.
Definition: ComponentBase.hh:249
Garfield::ComponentBase::m_yMirrorPeriodic
bool m_yMirrorPeriodic
Mirror periodicity in y.
Definition: ComponentBase.hh:255
Garfield::ComponentBase::m_xPeriodic
bool m_xPeriodic
Simple periodicity in x.
Definition: ComponentBase.hh:247
Garfield::ComponentBase::m_zPeriodic
bool m_zPeriodic
Simple periodicity in z.
Definition: ComponentBase.hh:251
Garfield::ComponentBase::m_className
std::string m_className
Class name.
Definition: ComponentBase.hh:233
Garfield::ComponentBase::m_xAxiallyPeriodic
bool m_xAxiallyPeriodic
Axial periodicity in x.
Definition: ComponentBase.hh:259
Garfield::ComponentBase::m_debug
bool m_debug
Switch on/off debugging messages.
Definition: ComponentBase.hh:274
Garfield::ComponentBase::m_xMirrorPeriodic
bool m_xMirrorPeriodic
Mirror periodicity in x.
Definition: ComponentBase.hh:253
Garfield::ComponentFieldMap::DriftMedium
void DriftMedium(const unsigned int imat)
Flag a field map material as a drift medium.
Definition: ComponentFieldMap.cc:66
Garfield::ComponentFieldMap::PrintRange
void PrintRange()
Show x, y, z, V and angular ranges.
Definition: ComponentFieldMap.cc:2354
Garfield::ComponentFieldMap::GetElementVolume
virtual double GetElementVolume(const unsigned int i)=0
Garfield::ComponentFieldMap::SetMedium
void SetMedium(const unsigned int imat, Medium *medium)
Associate a field map material with a Medium class.
Definition: ComponentFieldMap.cc:120
Garfield::ComponentFieldMap::PrintMaterials
void PrintMaterials()
List all currently defined materials.
Definition: ComponentFieldMap.cc:36
Garfield::ComponentFieldMap::ReadDouble
double ReadDouble(char *token, double def, bool &error)
Definition: ComponentFieldMap.cc:2645
Garfield::ComponentFieldMap::GetBoundingBox
virtual bool GetBoundingBox(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
Get the bounding box coordinates.
Definition: ComponentFieldMap.cc:2426
Garfield::ComponentFieldMap::MapCoordinates
void MapCoordinates(double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
Move (xpos, ypos, zpos) to field map coordinates.
Definition: ComponentFieldMap.cc:2441
Garfield::ComponentFieldMap::GetAspectRatio
virtual void GetAspectRatio(const unsigned int i, double &dmin, double &dmax)=0
Garfield::ComponentFieldMap::ReadInteger
int ReadInteger(char *token, int def, bool &error)
Definition: ComponentFieldMap.cc:2635
Garfield::ComponentFieldMap::FindElementCube
int FindElementCube(const double x, const double y, const double z, double &t1, double &t2, double &t3, TMatrixD *&jac, std::vector< TMatrixD * > &dN)
Find the element for a point in a cube.
Definition: ComponentFieldMap.cc:463
Garfield::ComponentFieldMap::PrintElement
void PrintElement(const std::string &header, const double x, const double y, const double z, const double t1, const double t2, const double t3, const double t4, const unsigned int i, const unsigned int n, const int iw=-1) const
Definition: ComponentFieldMap.cc:2740
Garfield::ComponentFieldMap::~ComponentFieldMap
virtual ~ComponentFieldMap()
Destructor.
Definition: ComponentFieldMap.cc:32
Garfield::ComponentFieldMap::GetConductivity
double GetConductivity(const unsigned int imat) const
Return the conductivity of a field map material.
Definition: ComponentFieldMap.cc:109
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)
Find the element for a point in curved quadratic quadrilaterals.
Definition: ComponentFieldMap.cc:166
Garfield::ComponentFieldMap::GetPermittivity
double GetPermittivity(const unsigned int imat) const
Return the permittivity of a field map material.
Definition: ComponentFieldMap.cc:98
Garfield::ComponentFieldMap::materials
std::vector< Material > materials
Definition: ComponentFieldMap.hh:133
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)
Find the element for a point in curved quadratic tetrahedra.
Definition: ComponentFieldMap.cc:331
Garfield::ComponentFieldMap::GetMedium
Medium * GetMedium(const unsigned int i) const
Return the Medium associated to a field map material.
Definition: ComponentFieldMap.cc:141
Garfield::ComponentFieldMap::NotDriftMedium
void NotDriftMedium(const unsigned int imat)
Flag a field map materials as a non-drift medium.
Definition: ComponentFieldMap.cc:82
Garfield::ComponentFieldMap::GetElement
bool GetElement(const unsigned int i, double &vol, double &dmin, double &dmax)
Return the volume and aspect ratio of a mesh element.
Definition: ComponentFieldMap.cc:152
Garfield::ComponentFieldMap::SetRange
virtual void SetRange()
Calculate x, y, z, V and angular ranges.
Definition: ComponentFieldMap.cc:2239
Garfield::ComponentFieldMap::elements
std::vector< Element > elements
Definition: ComponentFieldMap.hh:108
Garfield::ComponentFieldMap::IsInBoundingBox
virtual bool IsInBoundingBox(const double x, const double y, const double z) const
Definition: ComponentFieldMap.cc:2418
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) const
Move (ex, ey, ez) to global coordinates.
Definition: ComponentFieldMap.cc:2559
Garfield::ComponentFieldMap::PrintNotReady
void PrintNotReady(const std::string &header) const
Definition: ComponentFieldMap.hh:207
Garfield::Medium
Abstract base class for media.
Definition: Medium.hh:11
Garfield::Medium::GetName
const std::string & GetName() const
Definition: Medium.hh:22
Garfield::TetrahedralTree::GetTetListInBlock
std::vector< int > GetTetListInBlock(const Vec3 &point)
Definition: TetrahedralTree.cc:129
Garfield::TetrahedralTree::InsertTetrahedron
void InsertTetrahedron(const double elemBoundingBox[6], const int elemIndex)
Definition: TetrahedralTree.cc:108
Garfield::TetrahedralTree::InsertMeshNode
void InsertMeshNode(Vec3 point, const int nodeIndex)
Definition: TetrahedralTree.cc:68
Heed::pow
DoubleAc pow(const DoubleAc &f, double p)
Definition: DoubleAc.cpp:337
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314