Garfield++ 4.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
ComponentNeBem2d.cc
Go to the documentation of this file.
1#include <cmath>
2#include <fstream>
3#include <iostream>
4#include <cstdio>
5#include <algorithm>
6#include <numeric>
7#include <complex>
8
12#include "Garfield/Polygon.hh"
13#include "Garfield/Random.hh"
14
15namespace {
16
17double Mag2(const std::array<double, 2>& x) {
18 return x[0] * x[0] + x[1] * x[1];
19}
20
21/// Determine whether a point (u, v) lies on a straight line
22/// (x1, y1) to (x2, y2).
23bool OnLine(const double x1, const double y1, const double x2, const double y2,
24 const double u, const double v) {
25 // Set tolerances.
26 double epsx = 1.e-10 * std::max({fabs(x1), fabs(x2), fabs(u)});
27 double epsy = 1.e-10 * std::max({fabs(y1), fabs(y2), fabs(v)});
28 epsx = std::max(1.e-10, epsx);
29 epsy = std::max(1.e-10, epsy);
30
31 if ((fabs(x1 - u) <= epsx && fabs(y1 - v) <= epsy) ||
32 (fabs(x2 - u) <= epsx && fabs(y2 - v) <= epsy)) {
33 // Point to be examined coincides with start or end.
34 return true;
35 } else if (fabs(x1 - x2) <= epsx && fabs(y1 - y2) <= epsy) {
36 // The line (x1, y1) to (x2, y2) is in fact a point.
37 return false;
38 }
39 double xc = 0., yc = 0.;
40 if (fabs(u - x1) + fabs(v - y1) < fabs(u - x2) + fabs(v - y2)) {
41 // (u, v) is nearer to (x1, y1).
42 const double dx = (x2 - x1);
43 const double dy = (y2 - y1);
44 const double xl = ((u - x1) * dx + (v - y1) * dy) / (dx * dx + dy * dy);
45 if (xl < 0.) {
46 xc = x1;
47 yc = y1;
48 } else if (xl > 1.) {
49 xc = x2;
50 yc = y2;
51 } else {
52 xc = x1 + xl * dx;
53 yc = y1 + xl * dy;
54 }
55 } else {
56 // (u, v) is nearer to (x2, y2).
57 const double dx = (x1 - x2);
58 const double dy = (y1 - y2);
59 const double xl = ((u - x2) * dx + (v - y2) * dy) / (dx * dx + dy * dy);
60 if (xl < 0.) {
61 xc = x2;
62 yc = y2;
63 } else if (xl > 1.) {
64 xc = x1;
65 yc = y1;
66 } else {
67 xc = x2 + xl * dx;
68 yc = y2 + xl * dy;
69 }
70 }
71 // See whether the point is on the line.
72 if (fabs(u - xc) < epsx && fabs(v - yc) < epsy) {
73 return true;
74 }
75 return false;
76}
77
78/// Determine whether the 2 straight lines (x1, y1) to (x2, y2)
79/// and (u1, v1) to (u2, v2) cross at an intermediate point for both lines.
80bool Crossing(const double x1, const double y1, const double x2,
81 const double y2, const double u1, const double v1,
82 const double u2, const double v2, double& xc, double& yc) {
83 // Initial values.
84 xc = 0.;
85 yc = 0.;
86 /// Matrix to compute the crossing point.
87 std::array<std::array<double, 2>, 2> a;
88 a[0][0] = y2 - y1;
89 a[0][1] = v2 - v1;
90 a[1][0] = x1 - x2;
91 a[1][1] = u1 - u2;
92 const double det = a[0][0] * a[1][1] - a[1][0] * a[0][1];
93 // Set tolerances.
94 const double epsx = std::max(1.e-10, 1.e-10 * std::max({fabs(x1), fabs(x2),
95 fabs(u1), fabs(u2)}));
96 const double epsy = std::max(1.e-10, 1.e-10 * std::max({fabs(y1), fabs(y2),
97 fabs(v1), fabs(v2)}));
98 // Check if the lines are parallel.
99 if (fabs(det) < epsx * epsy) return false;
100
101 // Crossing, non-trivial lines: solve crossing equations.
102 const double aux = a[1][1];
103 a[1][1] = a[0][0] / det;
104 a[0][0] = aux / det;
105 a[1][0] = -a[1][0] / det;
106 a[0][1] = -a[0][1] / det;
107 // Compute crossing point.
108 xc = a[0][0] * (x1 * y2 - x2 * y1) + a[1][0] * (u1 * v2 - u2 * v1);
109 yc = a[0][1] * (x1 * y2 - x2 * y1) + a[1][1] * (u1 * v2 - u2 * v1);
110 // See whether the crossing point is on both lines.
111 if (OnLine(x1, y1, x2, y2, xc, yc) && OnLine(u1, v1, u2, v2, xc, yc)) {
112 // Intersecting lines.
113 return true;
114 }
115 // Crossing point not on both lines.
116 return false;
117}
118
119/// Determine whether two polygons are overlapping.
120bool Intersecting(const std::vector<double>& xp1,
121 const std::vector<double>& yp1,
122 const std::vector<double>& xp2,
123 const std::vector<double>& yp2) {
124
125 const double xmin1 = *std::min_element(std::begin(xp1), std::end(xp1));
126 const double ymin1 = *std::min_element(std::begin(yp1), std::end(yp1));
127 const double xmax1 = *std::max_element(std::begin(xp1), std::end(xp1));
128 const double ymax1 = *std::max_element(std::begin(yp1), std::end(yp1));
129
130 const double xmin2 = *std::min_element(std::begin(xp2), std::end(xp2));
131 const double ymin2 = *std::min_element(std::begin(yp2), std::end(yp2));
132 const double xmax2 = *std::max_element(std::begin(xp2), std::end(xp2));
133 const double ymax2 = *std::max_element(std::begin(yp2), std::end(yp2));
134
135 const double epsx = 1.e-6 * std::max({std::abs(xmax1), std::abs(xmin1),
136 std::abs(xmax2), std::abs(xmin2)});
137 const double epsy = 1.e-6 * std::max({std::abs(ymax1), std::abs(ymin1),
138 std::abs(ymax2), std::abs(ymin2)});
139 // Check if the bounding boxes overlap.
140 if (xmax1 + epsx < xmin2 || xmax2 + epsx < xmin1) return false;
141 if (ymax1 + epsy < ymin2 || ymax2 + epsy < ymin1) return false;
142
143 const unsigned int n1 = xp1.size();
144 const unsigned int n2 = xp2.size();
145 for (unsigned int i = 0; i < n1; ++i) {
146 const double x0 = xp1[i];
147 const double y0 = yp1[i];
148 const unsigned int ii = i < n1 - 1 ? i + 1 : 0;
149 const double x1 = xp1[ii];
150 const double y1 = yp1[ii];
151 for (unsigned int j = 0; j < n2; ++j) {
152 const unsigned int jj = j < n2 - 1 ? j + 1 : 0;
153 const double u0 = xp2[j];
154 const double v0 = yp2[j];
155 const double u1 = xp2[jj];
156 const double v1 = yp2[jj];
157 double xc = 0., yc = 0.;
158 if (!Crossing(x0, y0, x1, y1, u0, v0, u1, v1, xc, yc)) continue;
159 if ((OnLine(x0, y0, x1, y1, u0, v0) || OnLine(x0, y0, x1, y1, u1, v1)) &&
160 (OnLine(u0, v0, u1, v1, x0, y0) || OnLine(u0, v0, u1, v1, x1, y1))) {
161 continue;
162 }
163 return true;
164 }
165 }
166 return false;
167}
168
169/// Determine whether polygon 1 is fully enclosed by polygon 2.
170bool Enclosed(const std::vector<double>& xp1,
171 const std::vector<double>& yp1,
172 const std::vector<double>& xp2,
173 const std::vector<double>& yp2) {
174
175 const double xmin1 = *std::min_element(std::begin(xp1), std::end(xp1));
176 const double ymin1 = *std::min_element(std::begin(yp1), std::end(yp1));
177 const double xmax1 = *std::max_element(std::begin(xp1), std::end(xp1));
178 const double ymax1 = *std::max_element(std::begin(yp1), std::end(yp1));
179
180 const double xmin2 = *std::min_element(std::begin(xp2), std::end(xp2));
181 const double ymin2 = *std::min_element(std::begin(yp2), std::end(yp2));
182 const double xmax2 = *std::max_element(std::begin(xp2), std::end(xp2));
183 const double ymax2 = *std::max_element(std::begin(yp2), std::end(yp2));
184
185 const double epsx = 1.e-6 * std::max({std::abs(xmax1), std::abs(xmin1),
186 std::abs(xmax2), std::abs(xmin2)});
187 const double epsy = 1.e-6 * std::max({std::abs(ymax1), std::abs(ymin1),
188 std::abs(ymax2), std::abs(ymin2)});
189 // Check the bounding boxes.
190 if (xmax1 + epsx < xmin2 || xmax2 + epsx < xmin1) return false;
191 if (ymax1 + epsy < ymin2 || ymax2 + epsy < ymin1) return false;
192 if (xmin1 + epsx < xmin2 || xmax1 > xmax2 + epsx) return false;
193 if (ymin1 + epsy < ymin2 || ymax1 > ymax2 + epsy) return false;
194
195 const unsigned int n1 = xp1.size();
196 for (unsigned int i = 0; i < n1; ++i) {
197 bool inside = false, edge = false;
198 Garfield::Polygon::Inside(xp2, yp2, xp1[i], yp1[i], inside, edge);
199 if (!inside) return false;
200 }
201 return true;
202}
203
204}
205
206namespace Garfield {
207
208const double ComponentNeBem2d::InvEpsilon0 = 1. / VacuumPermittivity;
209const double ComponentNeBem2d::InvTwoPiEpsilon0 = 1. / TwoPiEpsilon0;
210
212
213void ComponentNeBem2d::ElectricField(const double x, const double y,
214 const double z, double& ex, double& ey,
215 double& ez, double& v, Medium*& m,
216 int& status) {
217 status = Field(x, y, z, ex, ey, ez, v, m, true);
218}
219
220void ComponentNeBem2d::ElectricField(const double x, const double y,
221 const double z, double& ex, double& ey,
222 double& ez, Medium*& m, int& status) {
223 double v = 0.;
224 status = Field(x, y, z, ex, ey, ez, v, m, false);
225}
226
227int ComponentNeBem2d::Field(const double x, const double y, const double z, double& ex, double& ey, double& ez, double& v,
228 Medium*& m, const bool opt) {
229
230 ex = ey = ez = 0.;
231 // Check if the requested point is inside the z-range.
232 if (m_useRangeZ && (z < m_zmin || z > m_zmax)) return -6;
233
234 // Check if the requested point is inside a medium.
235 m = GetMedium(x, y, z);
236 if (!m) return -6;
237
238 // Inside a conductor?
239 if (m->IsConductor()) {
240 if (!opt) return -5;
241 // Find the potential.
242 for (const auto& region : m_regions) {
243 bool inside = false, edge = false;
244 Garfield::Polygon::Inside(region.xv, region.yv, x, y, inside, edge);
245 if (inside || edge) {
246 v = region.bc.second;
247 break;
248 }
249 }
250 return -5;
251 }
252
253 if (!m_ready) {
254 if (!Initialise()) {
255 std::cerr << m_className << "::ElectricField: Initialisation failed.\n";
256 return -11;
257 }
258 }
259
260 // See whether we are inside a wire.
261 const unsigned int nWires = m_wires.size();
262 for (unsigned int i = 0; i < nWires; ++i) {
263 const double dx = x - m_wires[i].x;
264 const double dy = y - m_wires[i].y;
265 if (dx * dx + dy * dy < m_wires[i].r * m_wires[i].r) {
266 v = m_wires[i].v;
267 return i + 1;
268 }
269 }
270
271 // Sum up the contributions from all straight-line elements.
272 for (const auto& element : m_elements) {
273 const double cphi = element.cphi;
274 const double sphi = element.sphi;
275 // Transform to local coordinates.
276 double xL = 0., yL = 0.;
277 ToLocal(x - element.x, y - element.y, cphi, sphi, xL, yL);
278 if (opt) {
279 // Compute the potential.
280 v += LinePotential(element.a, xL, yL) * element.q;
281 }
282 // Compute the field in local coordinates.
283 double fx = 0., fy = 0.;
284 LineField(element.a, xL, yL, fx, fy);
285 // Rotate to the global frame.
286 ToGlobal(fx, fy, cphi, sphi, fx, fy);
287 ex += fx * element.q;
288 ey += fy * element.q;
289 }
290
291 // Add the contributions from the wires.
292 for (const auto& wire : m_wires) {
293 if (opt) {
294 // Compute the potential.
295 v += WirePotential(wire.r, x - wire.x, y - wire.y) * wire.q;
296 }
297 // Compute the field.
298 double fx = 0., fy = 0.;
299 WireField(wire.x, x - wire.x, y - wire.y, fx, fy);
300 ex += fx * wire.q;
301 ey += fy * wire.q;
302 }
303
304 for (const auto& box : m_spaceCharge) {
305 if (opt) {
306 v += BoxPotential(box.a, box.b, x - box.x, y - box.y, box.v0) * box.q;
307 }
308 double fx = 0., fy = 0.;
309 BoxField(box.a, box.b, x - box.x, y - box.y, fx, fy);
310 ex += fx * box.q;
311 ey += fy * box.q;
312 }
313 return 0;
314}
315
316bool ComponentNeBem2d::GetVoltageRange(double& vmin, double& vmax) {
317 bool gotValue = false;
318 for (const auto& region : m_regions) {
319 if (region.bc.first != Voltage) continue;
320 if (!gotValue) {
321 vmin = vmax = region.bc.second;
322 gotValue = true;
323 } else {
324 vmin = std::min(vmin, region.bc.second);
325 vmax = std::max(vmax, region.bc.second);
326 }
327 }
328
329 for (const auto& segment : m_segments) {
330 if (segment.bc.first != Voltage) continue;
331 if (!gotValue) {
332 vmin = vmax = segment.bc.second;
333 gotValue = true;
334 } else {
335 vmin = std::min(vmin, segment.bc.second);
336 vmax = std::max(vmax, segment.bc.second);
337 }
338 }
339
340 for (const auto& wire : m_wires) {
341 if (!gotValue) {
342 vmin = vmax = wire.v;
343 gotValue = true;
344 } else {
345 vmin = std::min(vmin, wire.v);
346 vmax = std::max(vmax, wire.v);
347 }
348 }
349 return gotValue;
350}
351
352Medium* ComponentNeBem2d::GetMedium(const double x, const double y,
353 const double z) {
354
355 if (m_geometry) return m_geometry->GetMedium(x, y, z);
356 for (const auto& region : m_regions) {
357 bool inside = false, edge = false;
358 Garfield::Polygon::Inside(region.xv, region.yv, x, y, inside, edge);
359 if (inside || edge) return region.medium;
360 }
361 return m_medium;
362}
363
365 double& xmin, double& ymin, double& zmin,
366 double& xmax, double& ymax, double& zmax) {
367
368 if (m_geometry) {
369 return m_geometry->GetBoundingBox(xmin, ymin, zmin, xmax, ymax, zmax);
370 }
371 return GetElementaryCell(xmin, ymin, zmin, xmax, ymax, zmax);
372}
373
375 double& xmin, double& ymin, double& zmin,
376 double& xmax, double& ymax, double& zmax) {
377
378 zmin = m_zmin;
379 zmax = m_zmax;
380 bool gotValue = false;
381 for (const auto& region : m_regions) {
382 const auto& xv = region.xv;
383 const auto& yv = region.yv;
384 if (!gotValue) {
385 xmin = *std::min_element(std::begin(xv), std::end(xv));
386 ymin = *std::min_element(std::begin(yv), std::end(yv));
387 xmax = *std::max_element(std::begin(xv), std::end(xv));
388 ymax = *std::max_element(std::begin(yv), std::end(yv));
389 gotValue = true;
390 } else {
391 xmin = std::min(*std::min_element(std::begin(xv), std::end(xv)), xmin);
392 ymin = std::min(*std::min_element(std::begin(yv), std::end(yv)), ymin);
393 xmax = std::max(*std::max_element(std::begin(xv), std::end(xv)), xmax);
394 ymax = std::max(*std::max_element(std::begin(yv), std::end(yv)), ymax);
395 }
396 }
397 for (const auto& seg : m_segments) {
398 if (!gotValue) {
399 xmin = std::min(seg.x0[0], seg.x1[0]);
400 xmax = std::max(seg.x0[0], seg.x1[0]);
401 ymin = std::min(seg.x0[1], seg.x1[1]);
402 ymax = std::max(seg.x0[1], seg.x1[1]);
403 gotValue = true;
404 } else {
405 xmin = std::min({xmin, seg.x0[0], seg.x1[0]});
406 xmax = std::max({xmax, seg.x0[0], seg.x1[0]});
407 ymin = std::min({ymin, seg.x0[1], seg.x1[1]});
408 ymax = std::max({ymax, seg.x0[1], seg.x1[1]});
409 }
410 }
411 for (const auto& wire : m_wires) {
412 if (!gotValue) {
413 xmin = xmax = wire.x;
414 ymin = ymax = wire.y;
415 } else {
416 xmin = std::min(xmin, wire.x);
417 xmax = std::max(xmax, wire.x);
418 ymin = std::min(ymin, wire.y);
419 ymax = std::max(ymax, wire.y);
420 }
421 }
422 return gotValue;
423}
424
425bool ComponentNeBem2d::IsWireCrossed(const double x0, const double y0,
426 const double z0, const double x1,
427 const double y1, const double z1,
428 double& xc, double& yc, double& zc,
429 const bool centre, double& rc) {
430 xc = x0;
431 yc = y0;
432 zc = z0;
433
434 if (m_wires.empty()) return false;
435
436 const double dx = x1 - x0;
437 const double dy = y1 - y0;
438 const double d2 = dx * dx + dy * dy;
439 // Make sure the step length is non-zero.
440 if (d2 < Small) return false;
441 const double invd2 = 1. / d2;
442
443 // Both coordinates are assumed to be located inside
444 // the drift area and inside a drift medium.
445 // This should have been checked before this call.
446
447 double dMin2 = 0.;
448 for (const auto& wire : m_wires) {
449 const double xw = wire.x;
450 const double yw = wire.y;
451 // Calculate the smallest distance between track and wire.
452 const double xIn0 = dx * (xw - x0) + dy * (yw - y0);
453 // Check if the minimum is located before (x0, y0).
454 if (xIn0 < 0.) continue;
455 const double xIn1 = -(dx * (xw - x1) + dy * (yw - y1));
456 // Check if the minimum is located behind (x1, y1).
457 if (xIn1 < 0.) continue;
458 // Minimum is located between (x0, y0) and (x1, y1).
459 const double xw0 = xw - x0;
460 const double xw1 = xw - x1;
461 const double yw0 = yw - y0;
462 const double yw1 = yw - y1;
463 const double dw02 = xw0 * xw0 + yw0 * yw0;
464 const double dw12 = xw1 * xw1 + yw1 * yw1;
465 if (xIn1 * xIn1 * dw02 > xIn0 * xIn0 * dw12) {
466 dMin2 = dw02 - xIn0 * xIn0 * invd2;
467 } else {
468 dMin2 = dw12 - xIn1 * xIn1 * invd2;
469 }
470 const double r2 = wire.r * wire.r;
471 if (dMin2 < r2) {
472 // Wire has been crossed.
473 if (centre) {
474 xc = xw;
475 yc = yw;
476 } else {
477 // Find the point of intersection.
478 const double p = -xIn0 * invd2;
479 const double q = (dw02 - r2) * invd2;
480 const double s = sqrt(p * p - q);
481 const double t = std::min(-p + s, -p - s);
482 xc = x0 + t * dx;
483 yc = y0 + t * dy;
484 zc = z0 + t * (z1 - z0);
485 }
486 rc = wire.r;
487 return true;
488 }
489 }
490 return false;
491}
492
493bool ComponentNeBem2d::IsInTrapRadius(const double q, const double x,
494 const double y, const double /*z*/,
495 double& xw, double& yw, double& rw) {
496
497 for (const auto& wire : m_wires) {
498 // Skip wires with the wrong charge.
499 if (q * wire.q > 0.) continue;
500 const double dx = wire.x - x;
501 const double dy = wire.y - y;
502 const double rt = wire.r * wire.ntrap;
503 if (dx * dx + dy * dy < rt * rt) {
504 xw = wire.x;
505 yw = wire.y;
506 rw = wire.r;
507 return true;
508 }
509 }
510 return false;
511}
512
513void ComponentNeBem2d::SetRangeZ(const double zmin, const double zmax) {
514
515 if (fabs(zmax - zmin) <= 0.) {
516 std::cerr << m_className << "::SetRangeZ: Zero range is not permitted.\n";
517 return;
518 }
519 m_zmin = std::min(zmin, zmax);
520 m_zmax = std::max(zmin, zmax);
521 m_useRangeZ = true;
522}
523
524bool ComponentNeBem2d::AddSegment(const double x0, const double y0,
525 const double x1, const double y1,
526 const double v, const int ndiv) {
527 const double dx = x1 - x0;
528 const double dy = y1 - y0;
529 if (dx * dx + dy * dy < Small) {
530 std::cerr << m_className << "::AddSegment: Length must be > 0.\n";
531 return false;
532 }
533
534 Segment segment;
535 segment.x0 = {x0, y0};
536 segment.x1 = {x1, y1};
537 segment.bc = std::make_pair(Voltage, v);
538 segment.region1 = -1;
539 segment.region2 = -1;
540 segment.ndiv = ndiv;
541 m_segments.push_back(std::move(segment));
542
543 if (m_debug) {
544 std::cout << m_className << "::AddSegment:\n ("
545 << x0 << ", " << y0 << ") - (" << x1 << ", " << y1 << ")\n"
546 << " Potential: " << v << " V\n";
547 }
548
549 m_ready = false;
550 return true;
551}
552
553bool ComponentNeBem2d::AddWire(const double x, const double y, const double d,
554 const double v, const int ntrap) {
555 if (d < Small) {
556 std::cerr << m_className << "::AddWire: Diameter must be > 0.\n";
557 return false;
558 }
559 if (ntrap <= 0) {
560 std::cerr << m_className << "::AddWire: Nbr. of trap radii must be > 0.\n";
561 return false;
562 }
563
564 Wire wire;
565 wire.x = x;
566 wire.y = y;
567 wire.r = 0.5 * d;
568 wire.v = v;
569 wire.q = 0.;
570 wire.ntrap = ntrap;
571 m_wires.push_back(std::move(wire));
572
573 if (m_debug) {
574 std::cout << m_className << "::AddWire:\n"
575 << " Centre: (" << x << ", " << y << ")\n"
576 << " Diameter: " << d << " cm\n"
577 << " Potential: " << v << " V\n";
578 }
579
580 m_ready = false;
581 return true;
582}
583
584bool ComponentNeBem2d::AddRegion(const std::vector<double>& xp,
585 const std::vector<double>& yp,
586 Medium* medium, const unsigned int bctype,
587 const double v, const int ndiv) {
588
589 if (xp.size() != yp.size()) {
590 std::cerr << m_className << "::AddRegion:\n"
591 << " Mismatch between number of x- and y-coordinates.\n";
592 return false;
593 }
594 if (xp.size() < 3) {
595 std::cerr << m_className << "::AddRegion: Too few points.\n";
596 return false;
597 }
598 if (bctype != 1 && bctype != 4) {
599 std::cerr << m_className << "::AddRegion: Invalid boundary condition.\n";
600 return false;
601 }
602
603 // Check if this is a valid polygon (no self-crossing).
604 const unsigned int np = xp.size();
605 if (np > 3) {
606 for (unsigned int i0 = 0; i0 < np; ++i0) {
607 const unsigned int i1 = i0 < np - 1 ? i0 + 1 : 0;
608 for (unsigned int j = 0; j < np - 3; ++j) {
609 const unsigned int j0 = i1 < np - 1 ? i1 + 1 : 0;
610 const unsigned int j1 = j0 < np - 1 ? j0 + 1 : 0;
611 double xc = 0., yc = 0.;
612 if (Crossing(xp[i0], yp[i0], xp[i1], yp[i1],
613 xp[j0], yp[j0], xp[j1], yp[j1], xc, yc)) {
614 std::cerr << m_className << "::AddRegion: Edges cross each other.\n";
615 return false;
616 }
617 }
618 }
619 }
620 std::vector<double> xv = xp;
621 std::vector<double> yv = yp;
622 const double xmin = *std::min_element(std::begin(xv), std::end(xv));
623 const double ymin = *std::min_element(std::begin(yv), std::end(yv));
624 const double xmax = *std::max_element(std::begin(xv), std::end(xv));
625 const double ymax = *std::max_element(std::begin(yv), std::end(yv));
626
627 const double epsx = 1.e-6 * std::max(std::abs(xmax), std::abs(xmin));
628 const double epsy = 1.e-6 * std::max(std::abs(ymax), std::abs(ymin));
629
630 const double f = Polygon::Area(xp, yp);
631 if (std::abs(f) < std::max(1.e-10, epsx * epsy)) {
632 std::cerr << m_className << "::AddRegion: Degenerate polygon.\n";
633 return false;
634 } else if (f > 0.) {
635 // Make sure all polygons have the same "handedness".
636 if (m_debug) {
637 std::cout << m_className << "::AddRegion: Reversing orientation.\n";
638 }
639 std::reverse(xv.begin(), xv.end());
640 std::reverse(yv.begin(), yv.end());
641 }
642 for (const auto& region : m_regions) {
643 if (Intersecting(xv, yv, region.xv, region.yv)) {
644 std::cerr << m_className << "::AddRegion:\n"
645 << " Polygon intersects an existing region.\n";
646 return false;
647 }
648 }
649 Region region;
650 region.xv = xv;
651 region.yv = yv;
652 region.medium = medium;
653 if (bctype == 1) {
654 region.bc = std::make_pair(Voltage, v);
655 } else if (bctype == 4) {
656 region.bc = std::make_pair(Dielectric, v);
657 }
658 region.depth = 0;
659 region.ndiv = ndiv;
660 m_regions.push_back(std::move(region));
661 return true;
662}
663
664void ComponentNeBem2d::AddChargeDistribution(const double x, const double y,
665 const double a, const double b,
666 const double rho) {
667 if (a < Small || b < Small) {
668 std::cerr << m_className << "::AddChargeDistribution:\n"
669 << " Lengths must be > 0.\n";
670 return;
671 }
672 const double a2 = a * a;
673 const double b2 = b * b;
674 const double v0 = -2 * (Pi * b2 + 2 * atan(b / a) * (a2 - b2));
675 SpaceCharge box;
676 box.x = x;
677 box.y = y;
678 box.a = a;
679 box.b = b;
680 box.q = rho;
681 box.v0 = v0;
682 m_spaceCharge.push_back(std::move(box));
683}
684
685void ComponentNeBem2d::SetNumberOfDivisions(const unsigned int ndiv) {
686 if (ndiv == 0) {
687 std::cerr << m_className << "::SetNumberOfDivisions:\n"
688 << " Number of divisions must be greater than zero.\n";
689 return;
690 }
691
692 m_nDivisions = ndiv;
693 m_ready = false;
694}
695
697 if (ncoll == 0) {
698 std::cerr << m_className << "::SetNumberOfCollocationPoints:\n"
699 << " Number of points must be greater than zero.\n";
700 return;
701 }
702
703 m_nCollocationPoints = ncoll;
704 m_ready = false;
705}
706
707void ComponentNeBem2d::SetMaxNumberOfIterations(const unsigned int niter) {
708 if (niter == 0) {
709 std::cerr << m_className << "::SetMaxNumberOfIterations:\n"
710 << " Number of iterations must be greater than zero.\n";
711 return;
712 }
713 m_nMaxIterations = niter;
714}
715
716bool ComponentNeBem2d::GetRegion(const unsigned int i,
717 std::vector<double>& xv,
718 std::vector<double>& yv,
719 Medium*& medium, unsigned int& bctype,
720 double& v) {
721 if (i >= m_regions.size()) return false;
722 if (!m_ready) {
723 if (!Initialise()) return false;
724 }
725 const auto& region = m_regions[i];
726 xv = region.xv;
727 yv = region.yv;
728 medium = region.medium;
729 bctype = region.bc.first == Voltage ? 1 : 4;
730 v = region.bc.second;
731 return true;
732}
733
734bool ComponentNeBem2d::GetSegment(const unsigned int i,
735 double& x0, double& y0, double& x1, double& y1, double& v) const {
736
737 if (i >= m_segments.size()) return false;
738 const auto& seg = m_segments[i];
739 x0 = seg.x0[0];
740 y0 = seg.x0[1];
741 x1 = seg.x1[0];
742 y1 = seg.x1[1];
743 v = seg.bc.second;
744 return true;
745}
746
747bool ComponentNeBem2d::GetWire(const unsigned int i,
748 double& x, double& y, double& d, double& v, double& q) const {
749
750 if (i >= m_wires.size()) return false;
751 const auto& wire = m_wires[i];
752 x = wire.x;
753 y = wire.y;
754 d = 2 * wire.r;
755 v = wire.v;
756 q = wire.q;
757 return true;
758}
759
760bool ComponentNeBem2d::GetElement(const unsigned int i,
761 double& x0, double& y0, double& x1, double& y1, double& q) const {
762
763 if (i >= m_elements.size()) return false;
764 const auto& element = m_elements[i];
765 ToGlobal(-element.a, 0., element.cphi, element.sphi, x0, y0);
766 ToGlobal( element.a, 0., element.cphi, element.sphi, x1, y1);
767 x0 += element.x;
768 y0 += element.y;
769 x1 += element.x;
770 y1 += element.y;
771 q = element.q;
772 return true;
773}
774
776
777 m_ready = false;
778 m_elements.clear();
779
780 double vmin = 0., vmax = 0.;
781 if (!GetVoltageRange(vmin, vmax)) {
782 std::cerr << m_className << "::Initialise:\n"
783 << " Could not determine the voltage range.\n";
784 return false;
785 }
786 if (fabs(vmin - vmax) < 1.e-6 * (vmin + vmax)) {
787 std::cerr << m_className << "::Initialise:\n"
788 << " All potentials are the same.\n";
789 return false;
790 }
791
792 if (m_debug) std::cout << m_className << "::Initialise:\n";
793 // Loop over the regions.
794 const unsigned int nRegions = m_regions.size();
795 if (m_debug) std::cout << " " << nRegions << " regions.\n";
796 std::vector<std::vector<unsigned int> > motherRegions(nRegions);
797 for (unsigned int i = 0; i < nRegions; ++i) {
798 auto& region = m_regions[i];
799 // Check if the region is fully enclosed by other ones.
800 for (unsigned int j = 0; j < nRegions; ++j) {
801 if (i == j) continue;
802 const auto& other = m_regions[j];
803 if (!Enclosed(region.xv, region.yv, other.xv, other.yv)) continue;
804 motherRegions[i].push_back(j);
805 if (m_debug) {
806 std::cout << " Region " << i << " is enclosed by region "
807 << j << ".\n";
808 }
809 }
810 region.depth = motherRegions[i].size();
811 }
812
813 std::vector<Segment> segments;
814 for (unsigned int i = 0; i < nRegions; ++i) {
815 const auto& region = m_regions[i];
816 int outerRegion = -1;
817 for (const unsigned int k : motherRegions[i]) {
818 if (outerRegion < 0) {
819 outerRegion = k;
820 } else if (m_regions[outerRegion].depth < m_regions[k].depth) {
821 outerRegion = k;
822 }
823 }
824 // Add the segments bounding this region.
825 const unsigned int n = region.xv.size();
826 for (unsigned int j = 0; j < n; ++j) {
827 const unsigned int k = j < n - 1 ? j + 1 : 0;
828 Segment seg;
829 seg.x0 = {region.xv[j], region.yv[j]};
830 seg.x1 = {region.xv[k], region.yv[k]};
831 seg.region1 = i;
832 seg.region2 = outerRegion;
833 seg.bc = region.bc;
834 seg.ndiv = region.ndiv;
835 segments.push_back(std::move(seg));
836 }
837 }
838 // Add the segments specified by the user.
839 segments.insert(segments.end(), m_segments.begin(), m_segments.end());
840 const unsigned int nSegments = segments.size();
841 if (m_debug) std::cout << " " << nSegments << " segments.\n";
842 std::vector<bool> done(nSegments, false);
843 // Look for overlaps.
844 for (unsigned int i = 0; i < nSegments; ++i) {
845 if (done[i]) continue;
846 if (m_debug) {
847 std::cout << " Segment " << i << ". ("
848 << segments[i].x0[0] << ", " << segments[i].x0[1] << ") - ("
849 << segments[i].x1[0] << ", " << segments[i].x1[1] << ")\n";
850 }
851 const double x0 = segments[i].x0[0];
852 const double x1 = segments[i].x1[0];
853 const double y0 = segments[i].x0[1];
854 const double y1 = segments[i].x1[1];
855 // Pick up all collinear segments.
856 std::vector<Segment> newSegments;
857 for (unsigned int j = i + 1; j < nSegments; ++j) {
858 const double u0 = segments[j].x0[0];
859 const double u1 = segments[j].x1[0];
860 const double v0 = segments[j].x0[1];
861 const double v1 = segments[j].x1[1];
862 const double epsx = std::max(1.e-10 * std::max({fabs(x0), fabs(x1),
863 fabs(u0), fabs(u1)}),
864 1.e-10);
865 const double epsy = std::max(1.e-10 * std::max({fabs(y0), fabs(y1),
866 fabs(v0), fabs(v1)}),
867 1.e-10);
868 const double a00 = y1 - y0;
869 const double a01 = v1 - v0;
870 const double a10 = x0 - x1;
871 const double a11 = u0 - u1;
872 const double det = a00 * a11 - a10 * a01;
873 const double tol = epsx * epsy;
874 // Skip non-parallel segments.
875 if (std::abs(det) > tol) continue;
876
877 if (std::abs(x0 * (y1 - v0) + x1 * (v0 - y0) + u0 * (y0 - y1)) > tol) {
878 continue;
879 }
880 newSegments.push_back(segments[j]);
881 done[j] = true;
882 }
883 newSegments.push_back(segments[i]);
884 if (newSegments.size() > 1) {
885 if (m_debug) {
886 std::cout << " Determining overlaps of " << newSegments.size()
887 << " collinear segments.\n";
888 }
889 EliminateOverlaps(newSegments);
890 if (m_debug) {
891 std::cout << " " << newSegments.size()
892 << " segments after splitting/merging.\n";
893 }
894 }
895 for (const auto& segment : newSegments) {
896 double lambda = 0.;
897 if (segment.bc.first == Dielectric) {
898 // Dielectric-dielectric interface.
899 const int reg1 = segment.region1;
900 const int reg2 = segment.region2;
901 double eps1 = 1.;
902 if (reg1 < 0) {
903 if (m_medium) eps1 = m_medium->GetDielectricConstant();
904 } else if (m_regions[reg1].medium) {
905 eps1 = m_regions[reg1].medium->GetDielectricConstant();
906 }
907 double eps2 = 1.;
908 if (reg2 < 0) {
909 if (m_medium) eps2 = m_medium->GetDielectricConstant();
910 } else if (m_regions[reg2].medium) {
911 eps2 = m_regions[reg2].medium->GetDielectricConstant();
912 }
913 if (fabs(eps1 - eps2) < 1.e-6 * (1. + fabs(eps1) + fabs(eps2))) {
914 if (m_debug) std::cout << " Same epsilon. Skip.\n";
915 continue;
916 }
917 // Compute lambda.
918 lambda = (eps2 - eps1) / (eps1 + eps2);
919 if (m_debug) std::cout << " Lambda = " << lambda << "\n";
920 }
921 const int ndiv = segment.ndiv <= 0 ? m_nDivisions : segment.ndiv;
922 Discretise(segment, m_elements, lambda, ndiv);
923 }
924 }
925 std::vector<std::vector<double> > influenceMatrix;
926 std::vector<std::vector<double> > inverseMatrix;
927
928 bool converged = false;
929 unsigned int nIter = 0;
930 while (!converged) {
931 ++nIter;
932 if (m_autoSize) {
933 std::cout << m_className << "::Initialise: Iteration " << nIter << "\n";
934 }
935 const unsigned int nElements = m_elements.size();
936 const unsigned int nEntries = nElements + m_wires.size() + 1;
937 if (m_debug) {
938 std::cout << " " << nElements << " elements.\n"
939 << " Matrix has " << nEntries << " rows/columns.\n";
940 }
941 // Compute the influence matrix.
942 influenceMatrix.assign(nEntries, std::vector<double>(nEntries, 0.));
943 if (!ComputeInfluenceMatrix(influenceMatrix)) {
944 std::cerr << m_className << "::Initialise:\n"
945 << " Error computing the influence matrix.\n";
946 return false;
947 }
948
949 // Invert the influence matrix.
950 inverseMatrix.assign(nEntries, std::vector<double>(nEntries, 0.));
951 if (!InvertMatrix(influenceMatrix, inverseMatrix)) {
952 std::cerr << m_className << "::Initialise: Matrix inversion failed.\n";
953 return false;
954 }
955 if (m_debug) std::cout << " Matrix inversion ok.\n";
956
957 // Compute the right hand side vector (boundary conditions).
958 std::vector<double> boundaryConditions(nEntries, 0.);
959 for (unsigned int i = 0; i < nElements; ++i) {
960 if (m_elements[i].bc.first == Voltage) {
961 boundaryConditions[i] = m_elements[i].bc.second;
962 for (const auto& box : m_spaceCharge) {
963 const double x = m_elements[i].x - box.x;
964 const double y = m_elements[i].y - box.y;
965 const double vs = BoxPotential(box.a, box.b, x, y, box.v0) * box.q;
966 boundaryConditions[i] -= vs;
967 }
968 } else {
969 for (const auto& box : m_spaceCharge) {
970 const double x = m_elements[i].x - box.x;
971 const double y = m_elements[i].y - box.y;
972 double fx = 0., fy = 0.;
973 BoxField(box.a, box.b, x, y, fx, fy);
974 // Rotate to the local frame of the target element.
975 ToLocal(fx, fy, m_elements[i].cphi, m_elements[i].sphi, fx, fy);
976 boundaryConditions[i] -= box.q * fy;
977 }
978 }
979 }
980 const unsigned int nWires = m_wires.size();
981 for (unsigned int i = 0; i < nWires; ++i) {
982 boundaryConditions[nElements + i] = m_wires[i].v;
983 for (const auto& box : m_spaceCharge) {
984 const double x = m_wires[i].x - box.x;
985 const double y = m_wires[i].y - box.y;
986 const double vs = BoxPotential(box.a, box.b, x, y, box.v0) * box.q;
987 boundaryConditions[nElements + i] -= vs;
988 }
989 }
990
991 double qsum = 0.;
992 for (const auto& box : m_spaceCharge) {
993 qsum += 4 * box.q * box.a * box.b;
994 }
995 boundaryConditions.back() = -qsum;
996
997 // Solve for the charge distribution.
998 if (!Solve(inverseMatrix, boundaryConditions)) {
999 std::cerr << m_className << "::Initialise: Solution failed.\n";
1000 return false;
1001 }
1002 if (m_debug) std::cout << " Solution ok.\n";
1003 const double tol = 1.e-6 * fabs(vmax - vmin);
1004 std::vector<bool> ok(nElements, true);
1005 converged = CheckConvergence(tol, ok);
1006 if (!m_autoSize) break;
1007 if (nIter >= m_nMaxIterations) break;
1008 for (unsigned int j = 0; j < nElements; ++j) {
1009 if (!ok[j]) {
1010 SplitElement(m_elements[j], m_elements);
1011 if (m_debug) std::cout << " Splitting element " << j << ".\n";
1012 }
1013 }
1014 }
1015 // Sort the regions by depth (innermost first).
1016 std::sort(m_regions.begin(), m_regions.end(),
1017 [](const Region& lhs, const Region& rhs) {
1018 return (lhs.depth > rhs.depth);
1019 });
1020 m_ready = true;
1021 return true;
1022}
1023
1024void ComponentNeBem2d::EliminateOverlaps(std::vector<Segment>& segments) {
1025
1026 if (segments.empty()) return;
1027 const unsigned int nIn = segments.size();
1028 // Find the first/last point along the line.
1029 std::array<double, 2> x0 = segments[0].x0;
1030 std::array<double, 2> x1 = segments[0].x1;
1031 // Use x or y coordinate depending on the orientation of the line.
1032 const unsigned int ic = fabs(x1[1] - x0[1]) > fabs(x1[0] - x0[0]) ? 1 : 0;
1033 std::vector<bool> swapped(nIn, false);
1034 for (unsigned int i = 0; i < nIn; ++i) {
1035 const auto& seg = segments[i];
1036 std::array<double, 2> u0 = seg.x0;
1037 std::array<double, 2> u1 = seg.x1;
1038 if (u0[ic] > u1[ic]) {
1039 // Swap points.
1040 std::swap(u0, u1);
1041 swapped[i] = true;
1042 }
1043 if (u0[ic] < x0[ic]) x0 = u0;
1044 if (u1[ic] > x1[ic]) x1 = u1;
1045 }
1046 const std::array<double, 2> d = {x1[0] - x0[0], x1[1] - x0[1]};
1047
1048 // Make a list of all points and their linear coordinate.
1049 std::vector<std::pair<double, std::vector<unsigned int> > > points;
1050 for (unsigned int i = 0; i < nIn; ++i) {
1051 for (const auto& xl : {segments[i].x0, segments[i].x1}) {
1052 const std::array<double, 2> d0 = {xl[0] - x0[0], xl[1] - x0[1]};
1053 const std::array<double, 2> d1 = {x1[0] - xl[0], x1[1] - xl[1]};
1054 double lambda = 0.;
1055 if (Mag2(d0) < Mag2(d1)) {
1056 // Point nearer to x0.
1057 lambda = d0[ic] / d[ic];
1058 } else {
1059 // Point nearer to p1.
1060 lambda = 1. - d1[ic] / d[ic];
1061 }
1062 // Add the point to the list.
1063 bool found = false;
1064 for (auto& point : points) {
1065 if (fabs(point.first - lambda) < 1.e-6) {
1066 found = true;
1067 point.second.push_back(i);
1068 break;
1069 }
1070 }
1071 if (found) continue;
1072 points.push_back(std::make_pair(lambda,
1073 std::vector<unsigned int>({i})));
1074 }
1075 }
1076 // Sort the points by linear coordinate.
1077 std::sort(std::begin(points), std::end(points));
1078
1079 std::vector<Segment> newSegments;
1080 const unsigned int nPoints = points.size();
1081 std::array<double, 2> xl = {x0[0] + points[0].first * d[0],
1082 x0[1] + points[0].first * d[1]};
1083 std::vector<unsigned int> left = points[0].second;
1084 for (unsigned int i = 1; i < nPoints; ++i) {
1085 Segment seg = segments[left.front()];
1086 seg.x0 = xl;
1087 xl = {x0[0] + points[i].first * d[0], x0[1] + points[i].first * d[1]};
1088 seg.x1 = xl;
1089 if (swapped[left.front()]) std::swap(seg.x0, seg.x1);
1090 // Sort out the boundary conditions.
1091 if (left.size() > 1) {
1092 for (unsigned int j = 1; j < left.size(); ++j) {
1093 const auto& other = segments[left[j]];
1094 if (seg.bc.first == Dielectric) {
1095 if (other.bc.first == Dielectric) {
1096 // Dielectric-dielectric interface.
1097 if ((seg.x1[ic] - seg.x0[ic]) * (other.x1[ic] - other.x0[ic]) > 0) {
1098 // Same orientation.
1099 continue;
1100 }
1101 seg.region2 = other.region1;
1102 } else {
1103 seg.bc = other.bc;
1104 }
1105 } else if (seg.bc.first != other.bc.first) {
1106 std::cerr << m_className << "::EliminateOverlaps:\n"
1107 << " Warning: conflicting boundary conditions.\n";
1108 }
1109 }
1110 }
1111 newSegments.push_back(std::move(seg));
1112 for (unsigned int k : points[i].second) {
1113 const auto it = std::find(left.begin(), left.end(), k);
1114 if (it == left.end()) {
1115 left.push_back(k);
1116 } else {
1117 left.erase(it);
1118 }
1119 }
1120 }
1121 segments.swap(newSegments);
1122}
1123
1124bool ComponentNeBem2d::Discretise(const Segment& seg,
1125 std::vector<Element>& elements,
1126 const double lambda,
1127 const unsigned int ndiv) {
1128
1129 if (ndiv < 1) {
1130 std::cerr << m_className << "::Discretise: Number of elements < 1.\n";
1131 return false;
1132 }
1133 const double phi = atan2(seg.x1[1] - seg.x0[1], seg.x1[0] - seg.x0[0]);
1134 const double cphi = cos(phi);
1135 const double sphi = sin(phi);
1136 const double dx = (seg.x1[0] - seg.x0[0]) / ndiv;
1137 const double dy = (seg.x1[1] - seg.x0[1]) / ndiv;
1138 const double a = 0.5 * sqrt(dx * dx + dy * dy);
1139 double x = seg.x0[0] - 0.5 * dx;
1140 double y = seg.x0[1] - 0.5 * dy;
1141 for (unsigned int i = 0; i < ndiv; ++i) {
1142 x += dx;
1143 y += dy;
1144 Element element;
1145 element.cphi = cphi;
1146 element.sphi = sphi;
1147 element.x = x;
1148 element.y = y;
1149 element.a = a;
1150 element.bc = seg.bc;
1151 element.lambda = lambda;
1152 elements.push_back(std::move(element));
1153 }
1154 return true;
1155}
1156
1157bool ComponentNeBem2d::ComputeInfluenceMatrix(
1158 std::vector<std::vector<double> >& infmat) const {
1159
1160 const unsigned int nL = m_elements.size();
1161 const unsigned int nE = nL + m_wires.size();
1162 // Loop over the target elements (F).
1163 for (unsigned int iF = 0; iF < nE; ++iF) {
1164 const auto bcF = iF < nL ? m_elements[iF].bc.first : Voltage;
1165 const double cphiF = iF < nL ? m_elements[iF].cphi : 1.;
1166 const double sphiF = iF < nL ? m_elements[iF].sphi : 0.;
1167 // Collocation point.
1168 const double xF = iF < nL ? m_elements[iF].x : m_wires[iF - nL].x;
1169 const double yF = iF < nL ? m_elements[iF].y : m_wires[iF - nL].y;
1170
1171 // Loop over the source elements (S).
1172 for (unsigned int jS = 0; jS < nE; ++jS) {
1173 // Calculate the influence coefficient.
1174 double infCoeff = 0.;
1175 if (jS < nL) {
1176 // Straight line element.
1177 const auto& src = m_elements[jS];
1178 double xL = 0., yL = 0.;
1179 ToLocal(xF - src.x, yF - src.y, src.cphi, src.sphi, xL, yL);
1180 if (bcF == Voltage) {
1181 infCoeff = LinePotential(src.a, xL, yL);
1182 } else if (bcF == Dielectric) {
1183 // Dielectric-dielectric interface.
1184 // Normal component of the displacement vector is continuous.
1185 if (iF == jS) {
1186 // Self-influence.
1187 infCoeff = 1. / (2. * src.lambda * VacuumPermittivity);
1188 } else {
1189 // Compute flux at the collocation point.
1190 double fx = 0., fy = 0.;
1191 LineField(src.a, xL, yL, fx, fy);
1192 // Rotate to the global frame.
1193 ToGlobal(fx, fy, src.cphi, src.sphi, fx, fy);
1194 // Rotate to the local frame of the target element.
1195 ToLocal(fx, fy, cphiF, sphiF, fx, fy);
1196 infCoeff = fy;
1197 }
1198 }
1199 } else {
1200 // Wire.
1201 const auto& src = m_wires[jS - nL];
1202 if (bcF == Voltage) {
1203 infCoeff = WirePotential(src.r, xF - src.x, yF - src.y);
1204 } else if (bcF == Dielectric) {
1205 double fx = 0., fy = 0.;
1206 WireField(src.r, xF - src.x, yF - src.y, fx, fy);
1207 ToLocal(fx, fy, cphiF, sphiF, fx, fy);
1208 infCoeff = fy;
1209 }
1210 }
1211 infmat[iF][jS] = infCoeff;
1212 }
1213 }
1214
1215 // Add charge neutrality condition.
1216 for (unsigned int i = 0; i < nE; ++i) {
1217 if (i < nL) {
1218 infmat[nE][i] = m_elements[i].a;
1219 } else {
1220 infmat[nE][i] = m_wires[i - nL].r;
1221 }
1222 infmat[i][nE] = 0.;
1223 }
1224 infmat[nE][nE] = 0.;
1225
1226 return true;
1227}
1228
1229void ComponentNeBem2d::SplitElement(Element& oldElement,
1230 std::vector<Element>& elements) {
1231 oldElement.a *= 0.5;
1232
1233 Element newElement = oldElement;
1234 double dx = 0., dy = 0.;
1235 ToGlobal(newElement.a, 0., newElement.cphi, newElement.sphi, dx, dy);
1236 oldElement.x += dx;
1237 oldElement.y += dy;
1238 newElement.x -= dx;
1239 newElement.y -= dx;
1240
1241 elements.push_back(std::move(newElement));
1242}
1243
1244bool ComponentNeBem2d::InvertMatrix(
1245 std::vector<std::vector<double> >& influenceMatrix,
1246 std::vector<std::vector<double> >& inverseMatrix) const {
1247
1248 const unsigned int nEntries = influenceMatrix.size();
1249
1250 // Temporary arrays for LU decomposition/substitution
1251 std::vector<double> col(nEntries, 0.);
1252 std::vector<int> index(nEntries, 0);
1253
1254 // Decompose the influence matrix
1255 if (!LUDecomposition(influenceMatrix, index)) {
1256 std::cerr << m_className << "::InvertMatrix: LU decomposition failed.\n";
1257 return false;
1258 }
1259
1260 // Initialise the inverse influence matrix
1261 inverseMatrix.assign(nEntries, std::vector<double>(nEntries, 0.));
1262 // Invert the matrix.
1263 for (unsigned int j = 0; j < nEntries; ++j) {
1264 col.assign(nEntries, 0.);
1265 col[j] = 1.;
1266 LUSubstitution(influenceMatrix, index, col);
1267 for (unsigned int i = 0; i < nEntries; ++i) inverseMatrix[i][j] = col[i];
1268 }
1269
1270 // Clear the influence matrix.
1271 influenceMatrix.clear();
1272
1273 return true;
1274}
1275
1276bool ComponentNeBem2d::LUDecomposition(std::vector<std::vector<double> >& mat,
1277 std::vector<int>& index) const {
1278 // The influence matrix is replaced by the LU decomposition of a rowwise
1279 // permutation of itself. The implementation is based on:
1280 // W. H. Press,
1281 // Numerical recipes in C++: the Art of Scientific Computing (version 2.11)
1282
1283 const unsigned int n = m_elements.size() + m_wires.size();
1284 // v stores the implicit scaling of each row
1285 std::vector<double> v(n, 0.);
1286
1287 // Loop over rows to get the implicit scaling information.
1288 for (unsigned int i = 0; i < n; ++i) {
1289 double big = 0.;
1290 for (unsigned int j = 0; j < n; ++j) {
1291 big = std::max(big, fabs(mat[i][j]));
1292 }
1293 if (big == 0.) return false;
1294 // Save the scaling
1295 v[i] = 1. / big;
1296 }
1297
1298 // Loop over columns
1299 unsigned int imax = 0;
1300 for (unsigned int j = 0; j < n; ++j) {
1301 for (unsigned int i = 0; i < j; ++i) {
1302 double sum = mat[i][j];
1303 for (unsigned int k = 0; k < i; ++k) {
1304 sum -= mat[i][k] * mat[k][j];
1305 }
1306 mat[i][j] = sum;
1307 }
1308 // Initialise for the search for the largest pivot element
1309 double big = 0.;
1310 for (unsigned int i = j; i < n; ++i) {
1311 double sum = mat[i][j];
1312 for (unsigned int k = 0; k < j; ++k) {
1313 sum -= mat[i][k] * mat[k][j];
1314 }
1315 mat[i][j] = sum;
1316 // Is the figure of merit for the pivot better than the best so far?
1317 const double dum = v[i] * fabs(sum);
1318 if (dum >= big) {
1319 big = dum;
1320 imax = i;
1321 }
1322 }
1323 // Do we need to interchange rows?
1324 if (j != imax) {
1325 for (unsigned k = 0; k < n; ++k) {
1326 const double dum = mat[imax][k];
1327 mat[imax][k] = mat[j][k];
1328 mat[j][k] = dum;
1329 }
1330 // Interchange the scale factor
1331 v[imax] = v[j];
1332 }
1333 index[j] = imax;
1334 if (mat[j][j] == 0.) mat[j][j] = Small;
1335 if (j != n - 1) {
1336 // Divide by the pivot element
1337 const double dum = 1. / mat[j][j];
1338 for (unsigned int i = j + 1; i < n; ++i) {
1339 mat[i][j] *= dum;
1340 }
1341 }
1342 }
1343
1344 return true;
1345}
1346
1347void ComponentNeBem2d::LUSubstitution(
1348 const std::vector<std::vector<double> >& mat, const std::vector<int>& index,
1349 std::vector<double>& col) const {
1350
1351 const unsigned int n = m_elements.size() + m_wires.size();
1352 unsigned int ii = 0;
1353 // Forward substitution
1354 for (unsigned i = 0; i < n; ++i) {
1355 const unsigned int ip = index[i];
1356 double sum = col[ip];
1357 col[ip] = col[i];
1358 if (ii != 0) {
1359 for (unsigned j = ii - 1; j < i; ++j) {
1360 sum -= mat[i][j] * col[j];
1361 }
1362 } else if (sum != 0.) {
1363 ii = i + 1;
1364 }
1365 col[i] = sum;
1366 }
1367
1368 // Backsubstitution
1369 for (int i = n - 1; i >= 0; i--) {
1370 double sum = col[i];
1371 for (unsigned j = i + 1; j < n; ++j) {
1372 sum -= mat[i][j] * col[j];
1373 }
1374 col[i] = sum / mat[i][i];
1375 }
1376}
1377
1378bool ComponentNeBem2d::Solve(const std::vector<std::vector<double> >& invmat,
1379 const std::vector<double>& bc) {
1380 const unsigned int nEntries = bc.size();
1381 const unsigned int nElements = m_elements.size();
1382 for (unsigned int i = 0; i < nElements; ++i) {
1383 double solution = 0.;
1384 for (unsigned int j = 0; j < nEntries; ++j) {
1385 solution += invmat[i][j] * bc[j];
1386 }
1387 m_elements[i].q = solution;
1388 }
1389 const unsigned int nWires = m_wires.size();
1390 for (unsigned int i = 0; i < nWires; ++i) {
1391 double solution = 0.;
1392 for (unsigned int j = 0; j < nEntries; ++j) {
1393 solution += invmat[nElements + i][j] * bc[j];
1394 }
1395 m_wires[i].q = solution;
1396 }
1397
1398 if (m_debug) {
1399 std::cout << m_className << "::Solve:\n Element Solution\n";
1400 for (unsigned int i = 0; i < nElements; ++i) {
1401 std::printf(" %8u %15.5f\n", i, m_elements[i].q);
1402 }
1403 if (!m_wires.empty()) {
1404 std::cout << " Wire Solution\n";
1405 for (unsigned int i = 0; i < nWires; ++i) {
1406 std::printf(" %8u %15.5f\n", i, m_wires[i].q);
1407 }
1408 }
1409 }
1410 return true;
1411}
1412
1413bool ComponentNeBem2d::CheckConvergence(const double tol,
1414 std::vector<bool>& ok) {
1415
1416 // Potential and normal component of the electric field
1417 // evaluated at the collocation points.
1418 std::vector<double> v(m_nCollocationPoints, 0.);
1419 std::vector<double> n(m_nCollocationPoints, 0.);
1420
1421 if (m_debug) {
1422 std::cout << m_className << "::CheckConvergence:\n"
1423 << " element # type LHS RHS\n";
1424 }
1425 const double scale = 1. / m_nCollocationPoints;
1426 unsigned int i = 0;
1427 for (const auto& tgt : m_elements) {
1428 v.assign(m_nCollocationPoints, 0.);
1429 n.assign(m_nCollocationPoints, 0.);
1430 double dx = 0., dy = 0.;
1431 ToGlobal(2 * tgt.a, 0., tgt.cphi, tgt.sphi, dx, dy);
1432 const double x0 = tgt.x - 0.5 * dx;
1433 const double y0 = tgt.y - 0.5 * dy;
1434 // Loop over the collocation points.
1435 for (unsigned int k = 0; k < m_nCollocationPoints; ++k) {
1436 double xG = x0;
1437 double yG = y0;
1438 if (m_randomCollocation) {
1439 const double r = RndmUniformPos();
1440 xG += r * dx;
1441 yG += r * dy;
1442 } else {
1443 const double s = (k + 1.) / (m_nCollocationPoints + 1.);
1444 xG += s * dx;
1445 yG += s * dy;
1446 }
1447 // Sum up the contributions from all boundary elements.
1448 for (const auto& src : m_elements) {
1449 double xL = 0., yL = 0.;
1450 // Transform to local coordinate system.
1451 ToLocal(xG - src.x, yG - src.y, src.cphi, src.sphi, xL, yL);
1452 // Compute the potential.
1453 v[k] += LinePotential(src.a, xL, yL) * src.q;
1454 // Compute the field.
1455 double fx = 0., fy = 0.;
1456 LineField(src.a, xL, yL, fx, fy);
1457 // Rotate to the global frame.
1458 ToGlobal(fx, fy, src.cphi, src.sphi, fx, fy);
1459 // Rotate to the local frame of the test element.
1460 ToLocal(fx, fy, tgt.cphi, tgt.sphi, fx, fy);
1461 n[k] += fy * src.q;
1462 }
1463
1464 for (const auto& src : m_wires) {
1465 // Compute the potential.
1466 v[k] += WirePotential(src.r, xG - src.x, yG - src.y) * src.q;
1467 // Compute the field.
1468 double fx = 0., fy = 0.;
1469 WireField(src.r, xG - src.x, yG - src.y, fx, fy);
1470 // Rotate to the local frame of the test element.
1471 ToLocal(fx, fy, tgt.cphi, tgt.sphi, fx, fy);
1472 n[k] += fy * src.q;
1473 }
1474
1475 for (const auto& box : m_spaceCharge) {
1476 const double xL = xG - box.x;
1477 const double yL = yG - box.y;
1478 v[k] += BoxPotential(box.a, box.b, xL, yL, box.v0) * box.q;
1479 double fx = 0., fy = 0.;
1480 BoxField(box.a, box.b, xL, yL, fx, fy);
1481 ToLocal(fx, fy, tgt.cphi, tgt.sphi, fx, fy);
1482 n[k] += fy * box.q;
1483 }
1484 }
1485 const double v0 = scale * std::accumulate(v.begin(), v.end(), 0.);
1486 const double n0 = scale * std::accumulate(n.begin(), n.end(), 0.);
1487 double n1 = 0.;
1488 if (tgt.bc.first == Voltage) {
1489 const double dv = v0 - tgt.bc.second;
1490 if (fabs(dv) > tol) ok[i] = false;
1491 if (m_debug) {
1492 std::printf(" %8u cond. %15.5f %15.5f %15.5f\n",
1493 i, v0, tgt.bc.second, dv);
1494 }
1495 } else if (tgt.bc.first == Dielectric) {
1496 // Dielectric-dielectric interface
1497 // TODO.
1498 n1 = n0 + 0.5 * InvEpsilon0 * tgt.q / tgt.lambda;
1499 if (m_debug) std::printf(" %8u diel. %15.5f %15.5f\n", i, n0, n1);
1500 }
1501 ++i;
1502 }
1503
1504 for (const auto& tgt : m_wires) {
1505 v.assign(m_nCollocationPoints, 0.);
1506 double x0 = tgt.x;
1507 double y0 = tgt.y;
1508 // Loop over the collocation points.
1509 for (unsigned int k = 0; k < m_nCollocationPoints; ++k) {
1510 const double phi = TwoPi * RndmUniform();
1511 const double xG = x0 + tgt.r * cos(phi);
1512 const double yG = y0 + tgt.r * sin(phi);
1513 // Sum up the contributions from all boundary elements.
1514 for (const auto& src : m_elements) {
1515 // Transform to local coordinate system.
1516 double xL = 0., yL = 0.;
1517 ToLocal(xG - src.x, yG - src.y, src.cphi, src.sphi, xL, yL);
1518 // Compute the potential.
1519 v[k] += LinePotential(src.a, xL, yL) * src.q;
1520 }
1521 for (const auto& src : m_wires) {
1522 v[k] += WirePotential(src.r, xG - src.x, yG - src.y) * src.q;
1523 }
1524 for (const auto& box : m_spaceCharge) {
1525 const double xL = xG - box.x;
1526 const double yL = yG - box.y;
1527 v[k] += BoxPotential(box.a, box.b, xL, yL, box.v0) * box.q;
1528 }
1529 }
1530 const double v0 = scale * std::accumulate(v.begin(), v.end(), 0.);
1531 if (m_debug) {
1532 std::printf(" %8u wire %15.5f %15.5f\n", i, v0, tgt.v);
1533 }
1534 ++i;
1535 }
1536
1537 return true;
1538}
1539
1540double ComponentNeBem2d::LinePotential(const double a, const double x,
1541 const double y) const {
1542 double p = 0.;
1543 const double amx = a - x;
1544 const double apx = a + x;
1545 if (fabs(y) > Small) {
1546 const double y2 = y * y;
1547 p = 2. * a - y * (atan(amx / y) + atan(apx / y)) -
1548 0.5 * amx * log(amx * amx + y2) - 0.5 * apx * log(apx * apx + y2);
1549 } else if (fabs(x) != a) {
1550 p = 2. * a - 0.5 * amx * log(amx * amx) - 0.5 * apx * log(apx * apx);
1551 } else {
1552 p = 2. * a * (1. - log(2. * a));
1553 }
1554
1555 return InvTwoPiEpsilon0 * p;
1556}
1557
1558double ComponentNeBem2d::WirePotential(const double r0, const double x,
1559 const double y) const {
1560 const double r = sqrt(x * x + y * y);
1561 if (r >= r0) {
1562 return -log(r) * r0 * InvEpsilon0;
1563 }
1564
1565 // Inside the wire the potential is constant
1566 return -log(r0) * r0 * InvEpsilon0;
1567}
1568
1569void ComponentNeBem2d::LineField(const double a,
1570 const double x, const double y,
1571 double& ex, double& ey) const {
1572 const double amx = a - x;
1573 const double apx = a + x;
1574 if (fabs(y) > 0.) {
1575 const double y2 = y * y;
1576 ex = 0.5 * log((apx * apx + y2) / (amx * amx + y2));
1577 ey = atan(amx / y) + atan(apx / y);
1578 } else if (fabs(x) != a) {
1579 ex = 0.5 * log(apx * apx / (amx * amx));
1580 ey = 0.;
1581 } else {
1582 // Singularity at the end points of the line
1583 constexpr double eps2 = 1.e-24;
1584 ex = 0.25 * log(pow(apx * apx - eps2, 2) / pow(amx * amx - eps2, 2));
1585 ey = 0.;
1586 }
1587 ex *= InvTwoPiEpsilon0;
1588 ey *= InvTwoPiEpsilon0;
1589}
1590
1591void ComponentNeBem2d::WireField(const double r0,
1592 const double x, const double y,
1593 double& ex, double& ey) const {
1594 const double r02 = r0 * r0;
1595 const double r2 = x * x + y * y;
1596 if (r2 > r02) {
1597 ex = x * r0 / r2;
1598 ey = y * r0 / r2;
1599 } else if (r2 == r02) {
1600 ex = 0.5 * x / r0;
1601 ey = 0.5 * y / r0;
1602 } else {
1603 // Inside the wire the field is zero.
1604 ex = ey = 0.;
1605 return;
1606 }
1607 ex *= InvEpsilon0;
1608 ey *= InvEpsilon0;
1609}
1610
1611double ComponentNeBem2d::BoxPotential(const double a, const double b,
1612 const double x, const double y,
1613 const double v0) const {
1614
1615 const double invc2 = 1. / (a * a + b * b);
1616 double v1 = 0., v2 = 0.;
1617 if (fabs(x) > a || fabs(y) > b) {
1618 // Outside the rectangle.
1619 const std::array<double, 2> dx = {x - a, x + a};
1620 const std::array<double, 2> dy = {y - b, y + b};
1621 const std::array<double, 2> dx2 = {dx[0] * dx[0], dx[1] * dx[1]};
1622 const std::array<double, 2> dy2 = {dy[0] * dy[0], dy[1] * dy[1]};
1623 v1 = dx[0] * dy[0] * log((dx2[0] + dy2[0]) * invc2)
1624 - dx[1] * dy[0] * log((dx2[1] + dy2[0]) * invc2)
1625 + dx[1] * dy[1] * log((dx2[1] + dy2[1]) * invc2)
1626 - dx[0] * dy[1] * log((dx2[0] + dy2[1]) * invc2);
1627 std::array<double, 4> alpha = {atan2(dy[0], dx[0]), atan2(dy[0], dx[1]),
1628 atan2(dy[1], dx[1]), atan2(dy[1], dx[0])};
1629 if (x < 0.) {
1630 for (size_t i = 0; i < 4; ++i) if (alpha[i] < 0.) alpha[i] += TwoPi;
1631 }
1632 v2 = dx2[0] * (alpha[0] - alpha[3]) + dx2[1] * (alpha[2] - alpha[1]) +
1633 dy2[0] * (alpha[1] - alpha[0]) + dy2[1] * (alpha[3] - alpha[2]);
1634 } else {
1635 // Inside the rectangle.
1636 const std::array<double, 2> dx = {a - x, a + x};
1637 const std::array<double, 2> dy = {b - y, b + y};
1638 const std::array<double, 2> dx2 = {dx[0] * dx[0], dx[1] * dx[1]};
1639 const std::array<double, 2> dy2 = {dy[0] * dy[0], dy[1] * dy[1]};
1640 v1 = dx[0] * dy[0] * log((dx2[0] + dy2[0]) * invc2) +
1641 dy[0] * dx[1] * log((dx2[1] + dy2[0]) * invc2) +
1642 dx[1] * dy[1] * log((dx2[1] + dy2[1]) * invc2) +
1643 dy[1] * dx[0] * log((dx2[0] + dy2[1]) * invc2);
1644 const double beta1 = atan2(dy[0], dx[0]);
1645 const double beta2 = atan2(dx[1], dy[0]);
1646 const double beta3 = atan2(dy[1], dx[1]);
1647 const double beta4 = atan2(dx[0], dy[1]);
1648 v2 = dx2[0] * (beta1 - beta4) + dy2[0] * (beta2 - beta1) +
1649 dx2[1] * (beta3 - beta2) + dy2[1] * (beta4 - beta3);
1650 v2 += HalfPi * (dx2[0] + dy2[0] + dx2[1] + dy2[1]);
1651 }
1652 return -InvTwoPiEpsilon0 * 0.5 * (v1 + v2 - v0);
1653}
1654
1655void ComponentNeBem2d::BoxField(const double a, const double b,
1656 const double x, const double y,
1657 double& ex, double& ey) const {
1658 const std::array<double, 2> dx = {x - a, x + a};
1659 const std::array<double, 2> dy = {y - b, y + b};
1660 const std::array<double, 2> dx2 = {dx[0] * dx[0], dx[1] * dx[1]};
1661 const std::array<double, 2> dy2 = {dy[0] * dy[0], dy[1] * dy[1]};
1662 const double r1 = dx2[0] + dy2[0];
1663 const double r2 = dx2[1] + dy2[0];
1664 const double r3 = dx2[1] + dy2[1];
1665 const double r4 = dx2[0] + dy2[1];
1666 ex = 0.5 * (dy[0] * log(r1 / r2) + dy[1] * log(r3 / r4));
1667 ey = 0.5 * (dx[0] * log(r1 / r4) + dx[1] * log(r3 / r2));
1668 if (fabs(x) > a || fabs(y) > b) {
1669 std::array<double, 4> alpha = {atan2(dy[0], dx[0]), atan2(dy[0], dx[1]),
1670 atan2(dy[1], dx[1]), atan2(dy[1], dx[0])};
1671 if (x < 0.) {
1672 for (size_t i = 0; i < 4; ++i) if (alpha[i] < 0.) alpha[i] += TwoPi;
1673 }
1674 ex += dx[0] * (alpha[0] - alpha[3]) + dx[1] * (alpha[2] - alpha[1]);
1675 ey -= dy[0] * (alpha[0] - alpha[1]) + dy[1] * (alpha[2] - alpha[3]);
1676 } else {
1677 const double beta1 = atan2(-dy[0], -dx[0]);
1678 const double beta2 = atan2( dx[1], -dy[0]);
1679 const double beta3 = atan2( dy[1], dx[1]);
1680 const double beta4 = atan2(-dx[0], dy[1]);
1681 ex += dx[0] * (HalfPi + beta1 - beta4) + dx[1] * (HalfPi + beta3 - beta2);
1682 ey -= dy[0] * (beta1 - beta2 - HalfPi) + dy[1] * (beta3 - beta4 - HalfPi);
1683 }
1684 ex *= InvTwoPiEpsilon0;
1685 ey *= InvTwoPiEpsilon0;
1686}
1687
1688void ComponentNeBem2d::Reset() {
1689 m_regions.clear();
1690 m_segments.clear();
1691 m_wires.clear();
1692 m_elements.clear();
1693 m_spaceCharge.clear();
1694 m_ready = false;
1695}
1696
1697void ComponentNeBem2d::UpdatePeriodicity() {
1698 std::cerr << m_className << "::UpdatePeriodicity:\n"
1699 << " Periodicities are not supported.\n";
1700}
1701
1702void ComponentNeBem2d::ToLocal(const double xIn, const double yIn,
1703 const double cphi, const double sphi,
1704 double& xOut, double& yOut) const {
1705 xOut = +cphi * xIn + sphi * yIn;
1706 yOut = -sphi * xIn + cphi * yIn;
1707}
1708
1709void ComponentNeBem2d::ToGlobal(const double xIn, const double yIn,
1710 const double cphi, const double sphi,
1711 double& xOut, double& yOut) const {
1712 xOut = cphi * xIn - sphi * yIn;
1713 yOut = sphi * xIn + cphi * yIn;
1714}
1715
1716}
void SetNumberOfDivisions(const unsigned int ndiv)
Set the default number of elements per segment.
bool AddSegment(const double x0, const double y0, const double x1, const double y1, const double v, const int ndiv=-1)
void ElectricField(const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status) override
bool GetVoltageRange(double &vmin, double &vmax) override
Calculate the voltage range [V].
bool GetSegment(const unsigned int i, double &x0, double &y0, double &x1, double &x2, double &v) const
Return the coordinates and voltage of a given straight-line segment.
bool GetBoundingBox(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax) override
Get the bounding box coordinates.
void SetRangeZ(const double zmin, const double zmax)
Set the extent of the drift region along z.
void SetNumberOfCollocationPoints(const unsigned int ncoll)
bool IsWireCrossed(const double x0, const double y0, const double z0, const double x1, const double y1, const double z1, double &xc, double &yc, double &zc, const bool centre, double &rc) override
bool Initialise()
Discretise the geometry and compute the solution.
bool GetRegion(const unsigned int i, std::vector< double > &xv, std::vector< double > &yv, Medium *&medium, unsigned int &bctype, double &v)
Return the properties of a given region.
Medium * GetMedium(const double x, const double y, const double z) override
Get the medium at a given location (x, y, z).
void AddChargeDistribution(const double x, const double y, const double a, const double b, const double rho)
bool GetWire(const unsigned int i, double &x, double &y, double &d, double &v, double &q) const
Return the coordinates, diameter, potential and charge of a given wire.
bool AddWire(const double x, const double y, const double d, const double v, const int ntrap=5)
void SetMaxNumberOfIterations(const unsigned int niter)
bool IsInTrapRadius(const double q0, const double x0, const double y0, const double z0, double &xw, double &yx, double &rw) override
bool AddRegion(const std::vector< double > &xp, const std::vector< double > &yp, Medium *medium, const unsigned int bctype=4, const double v=0., const int ndiv=-1)
bool GetElementaryCell(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax) override
Get the coordinates of the elementary cell.
bool GetElement(const unsigned int i, double &x0, double &y0, double &x1, double &y1, double &q) const
Return the coordinates and charge of a given boundary element.
Abstract base class for components.
Definition: Component.hh:13
bool m_debug
Switch on/off debugging messages.
Definition: Component.hh:341
std::string m_className
Class name.
Definition: Component.hh:329
Geometry * m_geometry
Pointer to the geometry.
Definition: Component.hh:332
bool m_ready
Ready for use?
Definition: Component.hh:338
virtual bool GetBoundingBox(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)=0
Get the bounding box (envelope of the geometry).
virtual Medium * GetMedium(const double x, const double y, const double z, const bool tesselated=false) const =0
Retrieve the medium at a given point.
Abstract base class for media.
Definition: Medium.hh:13
double GetDielectricConstant() const
Get the relative static dielectric constant.
Definition: Medium.hh:42
virtual bool IsConductor() const
Is this medium a conductor?
Definition: Medium.hh:29
double Area(const std::vector< double > &xp, const std::vector< double > &yp)
Determine the (signed) area of a polygon.
Definition: Polygon.cc:274
void Inside(const std::vector< double > &xpl, const std::vector< double > &ypl, const double x, const double y, bool &inside, bool &edge)
Definition: Polygon.cc:187
double RndmUniform()
Draw a random number uniformly distributed in the range [0, 1).
Definition: Random.hh:14
double RndmUniformPos()
Draw a random number uniformly distributed in the range (0, 1).
Definition: Random.hh:17
DoubleAc cos(const DoubleAc &f)
Definition: DoubleAc.cpp:432
DoubleAc pow(const DoubleAc &f, double p)
Definition: DoubleAc.cpp:337
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615
DoubleAc sin(const DoubleAc &f)
Definition: DoubleAc.cpp:384
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314
int Solve(void)
Definition: neBEM.c:2441
int InvertMatrix(void)
Definition: neBEM.c:1236
Definition: neBEM.h:155