Garfield++ 3.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Classes | Public Member Functions | List of all members
Garfield::ComponentNeBem2d Class Reference

Two-dimensional implementation of the nearly exact Boundary Element Method. More...

#include <ComponentNeBem2d.hh>

+ Inheritance diagram for Garfield::ComponentNeBem2d:

Public Member Functions

 ComponentNeBem2d ()
 Constructor.
 
 ~ComponentNeBem2d ()
 Destructor.
 
Medium * GetMedium (const double x, const double y, const double z) override
 Get the medium at a given location (x, y, z).
 
void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status) override
 
void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status) override
 Calculate the drift field [V/cm] and potential [V] at (x, y, z).
 
bool GetVoltageRange (double &vmin, double &vmax) override
 Calculate the voltage range [V].
 
bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax) override
 Get the bounding box coordinates.
 
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 IsInTrapRadius (const double q0, const double x0, const double y0, const double z0, double &xw, double &yx, double &rw) override
 
void SetMedium (Medium *medium)
 Set the "background" medium.
 
bool AddSegment (const double x0, const double y0, const double x1, const double y1, const double v, const int ndiv=-1)
 
bool AddWire (const double x, const double y, const double d, const double v, const int ntrap=5)
 
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)
 
void SetRangeZ (const double zmin, const double zmax)
 Set the extent of the drift region along z.
 
bool Initialise ()
 Discretise the geometry and compute the solution.
 
void SetNumberOfDivisions (const unsigned int ndiv)
 Set the default number of elements per segment.
 
void SetNumberOfCollocationPoints (const unsigned int ncoll)
 
void EnableAutoResizing (const bool on=true)
 
void EnableRandomCollocation (const bool on=true)
 
void SetMaxNumberOfIterations (const unsigned int niter)
 
unsigned int GetNumberOfRegions () const
 Return the number of regions.
 
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.
 
unsigned int GetNumberOfSegments () const
 Return the number of conducting straight-line segments.
 
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.
 
unsigned int GetNumberOfWires () const
 Return the number of wires.
 
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.
 
unsigned int GetNumberOfElements () const
 Return the number of boundary elements.
 
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.
 
- Public Member Functions inherited from Garfield::ComponentBase
 ComponentBase ()
 Constructor.
 
virtual ~ComponentBase ()
 Destructor.
 
virtual void SetGeometry (GeometryBase *geo)
 Define the geometry.
 
virtual void Clear ()
 Reset.
 
virtual Medium * GetMedium (const double x, const double y, const double z)
 Get the medium at a given location (x, y, z).
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status)=0
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status)=0
 Calculate the drift field [V/cm] and potential [V] at (x, y, z).
 
virtual bool GetVoltageRange (double &vmin, double &vmax)=0
 Calculate the voltage range [V].
 
virtual void WeightingField (const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label)
 
virtual double WeightingPotential (const double x, const double y, const double z, const std::string &label)
 
virtual void DelayedWeightingField (const double x, const double y, const double z, const double t, double &wx, double &wy, double &wz, const std::string &label)
 
virtual void MagneticField (const double x, const double y, const double z, double &bx, double &by, double &bz, int &status)
 
void SetMagneticField (const double bx, const double by, const double bz)
 Set a constant magnetic field.
 
virtual bool IsReady ()
 Ready for use?
 
virtual bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
 Get the bounding box coordinates.
 
double IntegrateFluxCircle (const double xc, const double yc, const double r, const unsigned int nI=50)
 
double IntegrateFluxSphere (const double xc, const double yc, const double zc, const double r, const unsigned int nI=20)
 
double IntegrateFlux (const double x0, const double y0, const double z0, const double dx1, const double dy1, const double dz1, const double dx2, const double dy2, const double dz2, const unsigned int nU=20, const unsigned int nV=20)
 
virtual 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)
 
virtual bool IsInTrapRadius (const double q0, const double x0, const double y0, const double z0, double &xw, double &yw, double &rw)
 
void EnablePeriodicityX (const bool on=true)
 Enable simple periodicity in the $x$ direction.
 
void DisablePeriodicityX ()
 
void EnablePeriodicityY (const bool on=true)
 Enable simple periodicity in the $y$ direction.
 
void DisablePeriodicityY ()
 
void EnablePeriodicityZ (const bool on=true)
 Enable simple periodicity in the $z$ direction.
 
void DisablePeriodicityZ ()
 
void EnableMirrorPeriodicityX (const bool on=true)
 Enable mirror periodicity in the $x$ direction.
 
void DisableMirrorPeriodicityX ()
 
void EnableMirrorPeriodicityY (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void DisableMirrorPeriodicityY ()
 
void EnableMirrorPeriodicityZ (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void DisableMirrorPeriodicityZ ()
 
void EnableAxialPeriodicityX (const bool on=true)
 Enable axial periodicity in the $x$ direction.
 
void DisableAxialPeriodicityX ()
 
void EnableAxialPeriodicityY (const bool on=true)
 Enable axial periodicity in the $y$ direction.
 
void DisableAxialPeriodicityY ()
 
void EnableAxialPeriodicityZ (const bool on=true)
 Enable axial periodicity in the $z$ direction.
 
void DisableAxialPeriodicityZ ()
 
void EnableRotationSymmetryX (const bool on=true)
 Enable rotation symmetry around the $x$ axis.
 
void DisableRotationSymmetryX ()
 
void EnableRotationSymmetryY (const bool on=true)
 Enable rotation symmetry around the $y$ axis.
 
void DisableRotationSymmetryY ()
 
void EnableRotationSymmetryZ (const bool on=true)
 Enable rotation symmetry around the $z$ axis.
 
void DisableRotationSymmetryZ ()
 
void EnableDebugging ()
 Switch on debugging messages.
 
void DisableDebugging ()
 Switch off debugging messages.
 
void ActivateTraps ()
 Request trapping to be taken care of by the component (for TCAD).
 
void DeactivateTraps ()
 
bool IsTrapActive ()
 
void ActivateVelocityMap ()
 Request velocity to be taken care of by the component (for TCAD).
 
void DectivateVelocityMap ()
 
bool IsVelocityActive ()
 
virtual bool ElectronAttachment (const double, const double, const double, double &eta)
 Get the electron attachment coefficient.
 
virtual bool HoleAttachment (const double, const double, const double, double &eta)
 Get the hole attachment coefficient.
 
virtual void ElectronVelocity (const double, const double, const double, double &vx, double &vy, double &vz, Medium *&, int &status)
 Get the electron drift velocity.
 
virtual void HoleVelocity (const double, const double, const double, double &vx, double &vy, double &vz, Medium *&, int &status)
 Get the hole drift velocity.
 
virtual bool GetElectronLifetime (const double, const double, const double, double &etau)
 
virtual bool GetHoleLifetime (const double, const double, const double, double &htau)
 

Additional Inherited Members

virtual void Reset ()=0
 Reset the component.
 
virtual void UpdatePeriodicity ()=0
 Verify periodicities.
 
- Protected Attributes inherited from Garfield::ComponentBase
std::string m_className = "ComponentBase"
 Class name.
 
GeometryBase * m_geometry = nullptr
 Pointer to the geometry.
 
bool m_ready = false
 Ready for use?
 
bool m_activeTraps = false
 Does the component have traps?
 
bool m_hasVelocityMap = false
 Does the component have velocity maps?
 
std::array< bool, 3 > m_periodic = {{false, false, false}}
 Simple periodicity in x, y, z.
 
std::array< bool, 3 > m_mirrorPeriodic = {{false, false, false}}
 Mirror periodicity in x, y, z.
 
std::array< bool, 3 > m_axiallyPeriodic = {{false, false, false}}
 Axial periodicity in x, y, z.
 
std::array< bool, 3 > m_rotationSymmetric = {{false, false, false}}
 Rotation symmetry around x-axis, y-axis, z-axis.
 
double m_bx0 = 0.
 
double m_by0 = 0.
 
double m_bz0 = 0.
 
bool m_debug = false
 Switch on/off debugging messages.
 

Detailed Description

Two-dimensional implementation of the nearly exact Boundary Element Method.

Definition at line 10 of file ComponentNeBem2d.hh.

Constructor & Destructor Documentation

◆ ComponentNeBem2d()

Garfield::ComponentNeBem2d::ComponentNeBem2d ( )

Constructor.

Definition at line 210 of file ComponentNeBem2d.cc.

210 : ComponentBase() {
211 m_className = "ComponentNeBem2d";
212}
Garfield::ComponentBase::ComponentBase
ComponentBase()
Constructor.
Definition: ComponentBase.cc:9
Garfield::ComponentBase::m_className
std::string m_className
Class name.
Definition: ComponentBase.hh:285

◆ ~ComponentNeBem2d()

Garfield::ComponentNeBem2d::~ComponentNeBem2d ( )
inline

Destructor.

Definition at line 15 of file ComponentNeBem2d.hh.

15{}

Member Function Documentation

◆ AddRegion()

bool Garfield::ComponentNeBem2d::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 
)

Add a region bounded by a polygon.

Parameters
xp,ypx/y-coordinates of the vertices of the polygon.
mediumpointer to the medium associated to the region.
bctype1: fixed voltage, 4: dielectric-dielectric interface.
vapplied potential.
ndivnumber of elements on each edge segment.

Definition at line 566 of file ComponentNeBem2d.cc.

569 {
570
571 if (xp.size() != yp.size()) {
572 std::cerr << m_className << "::AddRegion:\n"
573 << " Mismatch between number of x- and y-coordinates.\n";
574 return false;
575 }
576 if (xp.size() < 3) {
577 std::cerr << m_className << "::AddRegion: Too few points.\n";
578 return false;
579 }
580 if (bctype != 1 && bctype != 4) {
581 std::cerr << m_className << "::AddRegion: Invalid boundary condition.\n";
582 return false;
583 }
584
585 // Check if this is a valid polygon (no self-crossing).
586 const unsigned int np = xp.size();
587 if (np > 3) {
588 for (unsigned int i0 = 0; i0 < np; ++i0) {
589 const unsigned int i1 = i0 < np - 1 ? i0 + 1 : 0;
590 for (unsigned int j = 0; j < np - 3; ++j) {
591 const unsigned int j0 = i1 < np - 1 ? i1 + 1 : 0;
592 const unsigned int j1 = j0 < np - 1 ? j0 + 1 : 0;
593 double xc = 0., yc = 0.;
594 if (Crossing(xp[i0], yp[i0], xp[i1], yp[i1],
595 xp[j0], yp[j0], xp[j1], yp[j1], xc, yc)) {
596 std::cerr << m_className << "::AddRegion: Edges cross each other.\n";
597 return false;
598 }
599 }
600 }
601 }
602 std::vector<double> xv = xp;
603 std::vector<double> yv = yp;
604 const double xmin = *std::min_element(std::begin(xv), std::end(xv));
605 const double ymin = *std::min_element(std::begin(yv), std::end(yv));
606 const double xmax = *std::max_element(std::begin(xv), std::end(xv));
607 const double ymax = *std::max_element(std::begin(yv), std::end(yv));
608
609 const double epsx = 1.e-6 * std::max(std::abs(xmax), std::abs(xmin));
610 const double epsy = 1.e-6 * std::max(std::abs(ymax), std::abs(ymin));
611
612 const double f = Polygon::Area(xp, yp);
613 if (std::abs(f) < std::max(1.e-10, epsx * epsy)) {
614 std::cerr << m_className << "::AddRegion: Degenerate polygon.\n";
615 return false;
616 } else if (f > 0.) {
617 // Make sure all polygons have the same "handedness".
618 if (m_debug) {
619 std::cout << m_className << "::AddRegion: Reversing orientation.\n";
620 }
621 std::reverse(xv.begin(), xv.end());
622 std::reverse(yv.begin(), yv.end());
623 }
624 for (const auto& region : m_regions) {
625 if (Intersecting(xv, yv, region.xv, region.yv)) {
626 std::cerr << m_className << "::AddRegion:\n"
627 << " Polygon intersects an existing region.\n";
628 return false;
629 }
630 }
631 Region region;
632 region.xv = xv;
633 region.yv = yv;
634 region.medium = medium;
635 if (bctype == 1) {
636 region.bc = std::make_pair(Voltage, v);
637 } else if (bctype == 4) {
638 region.bc = std::make_pair(Dielectric, v);
639 }
640 region.depth = 0;
641 region.ndiv = ndiv;
642 m_regions.push_back(std::move(region));
643 return true;
644}
Garfield::ComponentBase::m_debug
bool m_debug
Switch on/off debugging messages.
Definition: ComponentBase.hh:310
Garfield::Polygon::Area
double Area(const std::vector< double > &xp, const std::vector< double > &yp)
Determine the (signed) area of a polygon.
Definition: Polygon.cc:222

◆ AddSegment()

bool Garfield::ComponentNeBem2d::AddSegment ( const double  x0,
const double  y0,
const double  x1,
const double  y1,
const double  v,
const int  ndiv = -1 
)

Add a conducting straight-line segment.

Parameters
x0,y0,x1,y1coordinates of start and end point.
vapplied potential.
ndivnumber of elements in which to split the segment.

Definition at line 506 of file ComponentNeBem2d.cc.

508 {
509 const double dx = x1 - x0;
510 const double dy = y1 - y0;
511 if (dx * dx + dy * dy < Small) {
512 std::cerr << m_className << "::AddSegment: Length must be > 0.\n";
513 return false;
514 }
515
516 Segment segment;
517 segment.x0 = {x0, y0};
518 segment.x1 = {x1, y1};
519 segment.bc = std::make_pair(Voltage, v);
520 segment.region1 = -1;
521 segment.region2 = -1;
522 segment.ndiv = ndiv;
523 m_segments.push_back(std::move(segment));
524
525 if (m_debug) {
526 std::cout << m_className << "::AddSegment:\n ("
527 << x0 << ", " << y0 << ") - (" << x1 << ", " << y1 << ")\n"
528 << " Potential: " << v << " V\n";
529 }
530
531 m_ready = false;
532 return true;
533}
Garfield::ComponentBase::m_ready
bool m_ready
Ready for use?
Definition: ComponentBase.hh:291

◆ AddWire()

bool Garfield::ComponentNeBem2d::AddWire ( const double  x,
const double  y,
const double  d,
const double  v,
const int  ntrap = 5 
)

Add a wire.

Parameters
x,ycentre of the wire.
dwire diameter.
vapplied potential.
ntrapmultiple of the wire radius within which a particle is considered to be trapped by the wire.

Definition at line 535 of file ComponentNeBem2d.cc.

536 {
537 if (d < Small) {
538 std::cerr << m_className << "::AddWire: Diameter must be > 0.\n";
539 return false;
540 }
541 if (ntrap <= 0) {
542 std::cerr << m_className << "::AddWire: Nbr. of trap radii must be > 0.\n";
543 return false;
544 }
545
546 Wire wire;
547 wire.x = x;
548 wire.y = y;
549 wire.r = 0.5 * d;
550 wire.v = v;
551 wire.q = 0.;
552 wire.ntrap = ntrap;
553 m_wires.push_back(std::move(wire));
554
555 if (m_debug) {
556 std::cout << m_className << "::AddWire:\n"
557 << " Centre: (" << x << ", " << y << ")\n"
558 << " Diameter: " << d << " cm\n"
559 << " Potential: " << v << " V\n";
560 }
561
562 m_ready = false;
563 return true;
564}

◆ ElectricField() [1/2]

void Garfield::ComponentNeBem2d::ElectricField ( const double  x,
const double  y,
const double  z,
double &  ex,
double &  ey,
double &  ez,
double &  v,
Medium *&  m,
int &  status 
)
overridevirtual

Calculate the drift field [V/cm] and potential [V] at (x, y, z).

Implements Garfield::ComponentBase.

Definition at line 214 of file ComponentNeBem2d.cc.

217 {
218 ex = ey = ez = v = 0.;
219 status = 0;
220 // Check if the requested point is inside the z-range.
221 if (m_useRangeZ && (z < m_zmin || z > m_zmax)) {
222 status = -6;
223 return;
224 }
225
226 // Check if the requested point is inside a medium.
227 m = GetMedium(x, y, z);
228 if (!m) {
229 status = -6;
230 return;
231 }
232 // Inside a conductor?
233 if (m->IsConductor()) {
234 status = -5;
235 // Find the potential.
236 for (const auto& region : m_regions) {
237 bool inside = false, edge = false;
238 Garfield::Polygon::Inside(region.xv, region.yv, x, y, inside, edge);
239 if (inside || edge) {
240 v = region.bc.second;
241 break;
242 }
243 }
244 return;
245 }
246
247 if (!m_ready) {
248 if (!Initialise()) {
249 std::cerr << m_className << "::ElectricField: Initialisation failed.\n";
250 status = -11;
251 return;
252 }
253 }
254
255 // See whether we are inside a wire.
256 const unsigned int nWires = m_wires.size();
257 for (unsigned int i = 0; i < nWires; ++i) {
258 const double dx = x - m_wires[i].x;
259 const double dy = y - m_wires[i].y;
260 if (dx * dx + dy * dy < m_wires[i].r * m_wires[i].r) {
261 v = m_wires[i].v;
262 status = i + 1;
263 return;
264 }
265 }
266
267 // Sum up the contributions from all straight-line elements.
268 for (const auto& element : m_elements) {
269 const double cphi = element.cphi;
270 const double sphi = element.sphi;
271 // Transform to local coordinates.
272 double xL = 0., yL = 0.;
273 ToLocal(x - element.x, y - element.y, cphi, sphi, xL, yL);
274 // Compute the potential.
275 v += LinePotential(element.a, xL, yL) * element.q;
276 // Compute the field in local coordinates.
277 double fx = 0., fy = 0.;
278 LineFlux(element.a, xL, yL, fx, fy);
279 // Rotate to the global frame.
280 ToGlobal(fx, fy, cphi, sphi, fx, fy);
281 ex += fx * element.q;
282 ey += fy * element.q;
283 }
284
285 // Add the contributions from the wires.
286 for (const auto& wire : m_wires) {
287 // Compute the potential.
288 v += WirePotential(wire.r, x - wire.x, y - wire.y) * wire.q;
289 // Compute the field.
290 double fx = 0., fy = 0.;
291 WireFlux(wire.x, x - wire.x, y - wire.y, fx, fy);
292 ex += fx * wire.q;
293 ey += fy * wire.q;
294 }
295
296}
Garfield::ComponentNeBem2d::Initialise
bool Initialise()
Discretise the geometry and compute the solution.
Definition: ComponentNeBem2d.cc:736
Garfield::ComponentNeBem2d::GetMedium
Medium * GetMedium(const double x, const double y, const double z) override
Get the medium at a given location (x, y, z).
Definition: ComponentNeBem2d.cc:341
Garfield::Polygon::Inside
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:135

◆ ElectricField() [2/2]

void Garfield::ComponentNeBem2d::ElectricField ( const double  x,
const double  y,
const double  z,
double &  ex,
double &  ey,
double &  ez,
Medium *&  m,
int &  status 
)
overridevirtual

Calculate the drift field at given point.

Parameters
x,y,zcoordinates [cm].
ex,ey,ezcomponents of the electric field [V/cm].
mpointer to the medium at this location.
statusstatus flag

Status flags: 

        0: Inside an active medium
      > 0: Inside a wire of type X
-4 ... -1: On the side of a plane where no wires are
       -5: Inside the mesh but not in an active medium
       -6: Outside the mesh
      -10: Unknown potential type (should not occur)
    other: Other cases (should not occur)

Implements Garfield::ComponentBase.

Definition at line 298 of file ComponentNeBem2d.cc.

300 {
301 double v = 0.;
302 ElectricField(x, y, z, ex, ey, ez, v, m, status);
303}
Garfield::ComponentNeBem2d::ElectricField
void ElectricField(const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status) override
Definition: ComponentNeBem2d.cc:298

Referenced by ElectricField().

◆ EnableAutoResizing()

void Garfield::ComponentNeBem2d::EnableAutoResizing ( const bool  on = true)
inline

Definition at line 77 of file ComponentNeBem2d.hh.

77{ m_autoSize = on; }

◆ EnableRandomCollocation()

void Garfield::ComponentNeBem2d::EnableRandomCollocation ( const bool  on = true)
inline

Definition at line 78 of file ComponentNeBem2d.hh.

78 {
79 m_randomCollocation = on;
80 }

◆ GetBoundingBox()

bool Garfield::ComponentNeBem2d::GetBoundingBox ( double &  xmin,
double &  ymin,
double &  zmin,
double &  xmax,
double &  ymax,
double &  zmax 
)
overridevirtual

Get the bounding box coordinates.

Reimplemented from Garfield::ComponentBase.

Definition at line 353 of file ComponentNeBem2d.cc.

355 {
356
357 if (m_geometry) {
358 return m_geometry->GetBoundingBox(xmin, ymin, zmin, xmax, ymax, zmax);
359 }
360 zmin = m_zmin;
361 zmax = m_zmax;
362 bool gotValue = false;
363 for (const auto& region : m_regions) {
364 const auto& xv = region.xv;
365 const auto& yv = region.yv;
366 if (!gotValue) {
367 xmin = *std::min_element(std::begin(xv), std::end(xv));
368 ymin = *std::min_element(std::begin(yv), std::end(yv));
369 xmax = *std::max_element(std::begin(xv), std::end(xv));
370 ymax = *std::max_element(std::begin(yv), std::end(yv));
371 gotValue = true;
372 } else {
373 xmin = std::min(*std::min_element(std::begin(xv), std::end(xv)), xmin);
374 ymin = std::min(*std::min_element(std::begin(yv), std::end(yv)), ymin);
375 xmax = std::max(*std::max_element(std::begin(xv), std::end(xv)), xmax);
376 ymax = std::max(*std::max_element(std::begin(yv), std::end(yv)), ymax);
377 }
378 }
379 for (const auto& seg : m_segments) {
380 if (!gotValue) {
381 xmin = std::min(seg.x0[0], seg.x1[0]);
382 xmax = std::max(seg.x0[0], seg.x1[0]);
383 ymin = std::min(seg.x0[1], seg.x1[1]);
384 ymax = std::max(seg.x0[1], seg.x1[1]);
385 gotValue = true;
386 } else {
387 xmin = std::min({xmin, seg.x0[0], seg.x1[0]});
388 xmax = std::max({xmax, seg.x0[0], seg.x1[0]});
389 ymin = std::min({ymin, seg.x0[1], seg.x1[1]});
390 ymax = std::max({ymax, seg.x0[1], seg.x1[1]});
391 }
392 }
393 for (const auto& wire : m_wires) {
394 if (!gotValue) {
395 xmin = xmax = wire.x;
396 ymin = ymax = wire.y;
397 } else {
398 xmin = std::min(xmin, wire.x);
399 xmax = std::max(xmax, wire.x);
400 ymin = std::min(ymin, wire.y);
401 ymax = std::max(ymax, wire.y);
402 }
403 }
404 return gotValue;
405}
Garfield::ComponentBase::m_geometry
GeometryBase * m_geometry
Pointer to the geometry.
Definition: ComponentBase.hh:288
Garfield::GeometryBase::GetBoundingBox
virtual bool GetBoundingBox(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)=0
Get the bounding box (envelope of the geometry).

◆ GetElement()

bool Garfield::ComponentNeBem2d::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.

Definition at line 721 of file ComponentNeBem2d.cc.

722 {
723
724 if (i >= m_elements.size()) return false;
725 const auto& element = m_elements[i];
726 ToGlobal(-element.a, 0., element.cphi, element.sphi, x0, y0);
727 ToGlobal( element.a, 0., element.cphi, element.sphi, x1, y1);
728 x0 += element.x;
729 y0 += element.y;
730 x1 += element.x;
731 y1 += element.y;
732 q = element.q;
733 return true;
734}

◆ GetMedium()

Medium * Garfield::ComponentNeBem2d::GetMedium ( const double  x,
const double  y,
const double  z 
)
overridevirtual

Get the medium at a given location (x, y, z).

Reimplemented from Garfield::ComponentBase.

Definition at line 341 of file ComponentNeBem2d.cc.

342 {
343
344 if (m_geometry) return m_geometry->GetMedium(x, y, z);
345 for (const auto& region : m_regions) {
346 bool inside = false, edge = false;
347 Garfield::Polygon::Inside(region.xv, region.yv, x, y, inside, edge);
348 if (inside || edge) return region.medium;
349 }
350 return m_medium;
351}
Garfield::GeometryBase::GetMedium
virtual Medium * GetMedium(const double x, const double y, const double z) const =0
Retrieve the medium at a given point.

Referenced by ElectricField().

◆ GetNumberOfElements()

unsigned int Garfield::ComponentNeBem2d::GetNumberOfElements ( ) const
inline

Return the number of boundary elements.

Definition at line 100 of file ComponentNeBem2d.hh.

100{ return m_elements.size(); }

◆ GetNumberOfRegions()

unsigned int Garfield::ComponentNeBem2d::GetNumberOfRegions ( ) const
inline

Return the number of regions.

Definition at line 84 of file ComponentNeBem2d.hh.

84{ return m_regions.size(); }

◆ GetNumberOfSegments()

unsigned int Garfield::ComponentNeBem2d::GetNumberOfSegments ( ) const
inline

Return the number of conducting straight-line segments.

Definition at line 90 of file ComponentNeBem2d.hh.

90{ return m_segments.size(); }

◆ GetNumberOfWires()

unsigned int Garfield::ComponentNeBem2d::GetNumberOfWires ( ) const
inline

Return the number of wires.

Definition at line 95 of file ComponentNeBem2d.hh.

95{ return m_wires.size(); }

◆ GetRegion()

bool Garfield::ComponentNeBem2d::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.

Definition at line 677 of file ComponentNeBem2d.cc.

681 {
682 if (i >= m_regions.size()) return false;
683 if (!m_ready) {
684 if (!Initialise()) return false;
685 }
686 const auto& region = m_regions[i];
687 xv = region.xv;
688 yv = region.yv;
689 medium = region.medium;
690 bctype = region.bc.first == Voltage ? 1 : 4;
691 v = region.bc.second;
692 return true;
693}

◆ GetSegment()

bool Garfield::ComponentNeBem2d::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.

Definition at line 695 of file ComponentNeBem2d.cc.

696 {
697
698 if (i >= m_segments.size()) return false;
699 const auto& seg = m_segments[i];
700 x0 = seg.x0[0];
701 y0 = seg.x0[1];
702 x1 = seg.x1[0];
703 y1 = seg.x1[1];
704 v = seg.bc.second;
705 return true;
706}

◆ GetVoltageRange()

bool Garfield::ComponentNeBem2d::GetVoltageRange ( double &  vmin,
double &  vmax 
)
overridevirtual

Calculate the voltage range [V].

Implements Garfield::ComponentBase.

Definition at line 305 of file ComponentNeBem2d.cc.

305 {
306 bool gotValue = false;
307 for (const auto& region : m_regions) {
308 if (region.bc.first != Voltage) continue;
309 if (!gotValue) {
310 vmin = vmax = region.bc.second;
311 gotValue = true;
312 } else {
313 vmin = std::min(vmin, region.bc.second);
314 vmax = std::max(vmax, region.bc.second);
315 }
316 }
317
318 for (const auto& segment : m_segments) {
319 if (segment.bc.first != Voltage) continue;
320 if (!gotValue) {
321 vmin = vmax = segment.bc.second;
322 gotValue = true;
323 } else {
324 vmin = std::min(vmin, segment.bc.second);
325 vmax = std::max(vmax, segment.bc.second);
326 }
327 }
328
329 for (const auto& wire : m_wires) {
330 if (!gotValue) {
331 vmin = vmax = wire.v;
332 gotValue = true;
333 } else {
334 vmin = std::min(vmin, wire.v);
335 vmax = std::max(vmax, wire.v);
336 }
337 }
338 return gotValue;
339}

Referenced by Initialise().

◆ GetWire()

bool Garfield::ComponentNeBem2d::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.

Definition at line 708 of file ComponentNeBem2d.cc.

709 {
710
711 if (i >= m_wires.size()) return false;
712 const auto& wire = m_wires[i];
713 x = wire.x;
714 y = wire.y;
715 d = 2 * wire.r;
716 v = wire.v;
717 q = wire.q;
718 return true;
719}

◆ Initialise()

bool Garfield::ComponentNeBem2d::Initialise ( )

Discretise the geometry and compute the solution.

Definition at line 736 of file ComponentNeBem2d.cc.

736 {
737
738 m_ready = false;
739 m_elements.clear();
740
741 double vmin = 0., vmax = 0.;
742 if (!GetVoltageRange(vmin, vmax)) {
743 std::cerr << m_className << "::Initialise:\n"
744 << " Could not determine the voltage range.\n";
745 return false;
746 }
747 if (fabs(vmin - vmax) < 1.e-6 * (vmin + vmax)) {
748 std::cerr << m_className << "::Initialise:\n"
749 << " All potentials are the same.\n";
750 return false;
751 }
752
753 if (m_debug) std::cout << m_className << "::Initialise:\n";
754 // Loop over the regions.
755 const unsigned int nRegions = m_regions.size();
756 if (m_debug) std::cout << " " << nRegions << " regions.\n";
757 std::vector<std::vector<unsigned int> > motherRegions(nRegions);
758 for (unsigned int i = 0; i < nRegions; ++i) {
759 auto& region = m_regions[i];
760 // Check if the region is fully enclosed by other ones.
761 for (unsigned int j = 0; j < nRegions; ++j) {
762 if (i == j) continue;
763 const auto& other = m_regions[j];
764 if (!Enclosed(region.xv, region.yv, other.xv, other.yv)) continue;
765 motherRegions[i].push_back(j);
766 if (m_debug) {
767 std::cout << " Region " << i << " is enclosed by region "
768 << j << ".\n";
769 }
770 }
771 region.depth = motherRegions[i].size();
772 }
773
774 std::vector<Segment> segments;
775 for (unsigned int i = 0; i < nRegions; ++i) {
776 const auto& region = m_regions[i];
777 int outerRegion = -1;
778 for (const unsigned int k : motherRegions[i]) {
779 if (outerRegion < 0) {
780 outerRegion = k;
781 } else if (m_regions[outerRegion].depth < m_regions[k].depth) {
782 outerRegion = k;
783 }
784 }
785 // Add the segments bounding this region.
786 const unsigned int n = region.xv.size();
787 for (unsigned int j = 0; j < n; ++j) {
788 const unsigned int k = j < n - 1 ? j + 1 : 0;
789 Segment seg;
790 seg.x0 = {region.xv[j], region.yv[j]};
791 seg.x1 = {region.xv[k], region.yv[k]};
792 seg.region1 = i;
793 seg.region2 = outerRegion;
794 seg.bc = region.bc;
795 seg.ndiv = region.ndiv;
796 segments.push_back(std::move(seg));
797 }
798 }
799 // Add the segments specified by the user.
800 segments.insert(segments.end(), m_segments.begin(), m_segments.end());
801 const unsigned int nSegments = segments.size();
802 if (m_debug) std::cout << " " << nSegments << " segments.\n";
803 std::vector<bool> done(nSegments, false);
804 // Look for overlaps.
805 for (unsigned int i = 0; i < nSegments; ++i) {
806 if (done[i]) continue;
807 if (m_debug) {
808 std::cout << " Segment " << i << ". ("
809 << segments[i].x0[0] << ", " << segments[i].x0[1] << ") - ("
810 << segments[i].x1[0] << ", " << segments[i].x1[1] << ")\n";
811 }
812 const double x0 = segments[i].x0[0];
813 const double x1 = segments[i].x1[0];
814 const double y0 = segments[i].x0[1];
815 const double y1 = segments[i].x1[1];
816 // Pick up all collinear segments.
817 std::vector<Segment> newSegments;
818 for (unsigned int j = i + 1; j < nSegments; ++j) {
819 const double u0 = segments[j].x0[0];
820 const double u1 = segments[j].x1[0];
821 const double v0 = segments[j].x0[1];
822 const double v1 = segments[j].x1[1];
823 const double epsx = std::max(1.e-10 * std::max({fabs(x0), fabs(x1),
824 fabs(u0), fabs(u1)}),
825 1.e-10);
826 const double epsy = std::max(1.e-10 * std::max({fabs(y0), fabs(y1),
827 fabs(v0), fabs(v1)}),
828 1.e-10);
829 const double a00 = y1 - y0;
830 const double a01 = v1 - v0;
831 const double a10 = x0 - x1;
832 const double a11 = u0 - u1;
833 const double det = a00 * a11 - a10 * a01;
834 const double tol = epsx * epsy;
835 // Skip non-parallel segments.
836 if (std::abs(det) > tol) continue;
837
838 if (std::abs(x0 * (y1 - v0) + x1 * (v0 - y0) + u0 * (y0 - y1)) > tol) {
839 continue;
840 }
841 newSegments.push_back(segments[j]);
842 done[j] = true;
843 }
844 newSegments.push_back(segments[i]);
845 if (newSegments.size() > 1) {
846 if (m_debug) {
847 std::cout << " Determining overlaps of " << newSegments.size()
848 << " collinear segments.\n";
849 }
850 EliminateOverlaps(newSegments);
851 if (m_debug) {
852 std::cout << " " << newSegments.size()
853 << " segments after splitting/merging.\n";
854 }
855 }
856 for (const auto& segment : newSegments) {
857 double lambda = 0.;
858 if (segment.bc.first == Dielectric) {
859 // Dielectric-dielectric interface.
860 const int reg1 = segment.region1;
861 const int reg2 = segment.region2;
862 double eps1 = 1.;
863 if (reg1 < 0) {
864 if (m_medium) eps1 = m_medium->GetDielectricConstant();
865 } else if (m_regions[reg1].medium) {
866 eps1 = m_regions[reg1].medium->GetDielectricConstant();
867 }
868 double eps2 = 1.;
869 if (reg2 < 0) {
870 if (m_medium) eps2 = m_medium->GetDielectricConstant();
871 } else if (m_regions[reg2].medium) {
872 eps2 = m_regions[reg2].medium->GetDielectricConstant();
873 }
874 if (fabs(eps1 - eps2) < 1.e-6 * (1. + fabs(eps1) + fabs(eps2))) {
875 if (m_debug) std::cout << " Same epsilon. Skip.\n";
876 continue;
877 }
878 // Compute lambda.
879 lambda = (eps2 - eps1) / (eps1 + eps2);
880 if (m_debug) std::cout << " Lambda = " << lambda << "\n";
881 }
882 const int ndiv = segment.ndiv <= 0 ? m_nDivisions : segment.ndiv;
883 Discretise(segment, m_elements, lambda, ndiv);
884 }
885 }
886 std::vector<std::vector<double> > influenceMatrix;
887 std::vector<std::vector<double> > inverseMatrix;
888
889 bool converged = false;
890 unsigned int nIter = 0;
891 while (!converged) {
892 ++nIter;
893 if (m_autoSize) {
894 std::cout << m_className << "::Initialise: Iteration " << nIter << "\n";
895 }
896 const unsigned int nElements = m_elements.size();
897 const unsigned int nEntries = nElements + m_wires.size() + 1;
898 if (m_debug) {
899 std::cout << " " << nElements << " elements.\n"
900 << " Matrix has " << nEntries << " rows/columns.\n";
901 }
902 // Compute the influence matrix.
903 influenceMatrix.assign(nEntries, std::vector<double>(nEntries, 0.));
904 if (!ComputeInfluenceMatrix(influenceMatrix)) {
905 std::cerr << m_className << "::Initialise:\n"
906 << " Error computing the influence matrix.\n";
907 return false;
908 }
909
910 // Invert the influence matrix.
911 inverseMatrix.assign(nEntries, std::vector<double>(nEntries, 0.));
912 if (!InvertMatrix(influenceMatrix, inverseMatrix)) {
913 std::cerr << m_className << "::Initialise: Matrix inversion failed.\n";
914 return false;
915 }
916 if (m_debug) std::cout << " Matrix inversion ok.\n";
917
918 // Compute the right hand side vector (boundary conditions).
919 std::vector<double> boundaryConditions(nEntries, 0.);
920 for (unsigned int i = 0; i < nElements; ++i) {
921 if (m_elements[i].bc.first != Voltage) continue;
922 boundaryConditions[i] = m_elements[i].bc.second;
923 }
924 const unsigned int nWires = m_wires.size();
925 for (unsigned int i = 0; i < nWires; ++i) {
926 boundaryConditions[nElements + i] = m_wires[i].v;
927 }
928
929 // Solve for the charge distribution.
930 if (!Solve(inverseMatrix, boundaryConditions)) {
931 std::cerr << m_className << "::Initialise: Solution failed.\n";
932 return false;
933 }
934 if (m_debug) std::cout << " Solution ok.\n";
935 const double tol = 1.e-6 * fabs(vmax - vmin);
936 std::vector<bool> ok(nElements, true);
937 converged = CheckConvergence(tol, ok);
938 if (!m_autoSize) break;
939 if (nIter >= m_nMaxIterations) break;
940 for (unsigned int j = 0; j < nElements; ++j) {
941 if (!ok[j]) {
942 SplitElement(m_elements[j], m_elements);
943 if (m_debug) std::cout << " Splitting element " << j << ".\n";
944 }
945 }
946 }
947 // Sort the regions by depth (innermost first).
948 std::sort(m_regions.begin(), m_regions.end(),
949 [](const Region& lhs, const Region& rhs) {
950 return (lhs.depth > rhs.depth);
951 });
952 m_ready = true;
953 return true;
954}
Garfield::ComponentNeBem2d::GetVoltageRange
bool GetVoltageRange(double &vmin, double &vmax) override
Calculate the voltage range [V].
Definition: ComponentNeBem2d.cc:305
Garfield::Medium::GetDielectricConstant
double GetDielectricConstant() const
Get the relative static dielectric constant.
Definition: Medium.hh:42
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615

Referenced by ElectricField(), and GetRegion().

◆ IsInTrapRadius()

bool Garfield::ComponentNeBem2d::IsInTrapRadius ( const double  q0,
const double  x0,
const double  y0,
const double  z0,
double &  xw,
double &  yw,
double &  rw 
)
overridevirtual

Determine whether a particle is inside the trap radius of a wire.

Parameters
q0charge of the particle [in elementary charges].
x0,y0,z0position [cm] of the particle.
xw,ywcoordinates of the wire (if applicable).
rwradius of the wire (if applicable).

Reimplemented from Garfield::ComponentBase.

Definition at line 475 of file ComponentNeBem2d.cc.

477 {
478
479 for (const auto& wire : m_wires) {
480 // Skip wires with the wrong charge.
481 if (q * wire.q > 0.) continue;
482 const double dx = wire.x - x;
483 const double dy = wire.y - y;
484 const double rt = wire.r * wire.ntrap;
485 if (dx * dx + dy * dy < rt * rt) {
486 xw = wire.x;
487 yw = wire.y;
488 rw = wire.r;
489 return true;
490 }
491 }
492 return false;
493}

◆ IsWireCrossed()

bool Garfield::ComponentNeBem2d::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 
)
overridevirtual

Determine whether the line between two points crosses a wire.

Parameters
x0,y0,z0first point [cm].
x1,y1,z1second point [cm]
xc,yc,zcpoint [cm] where the line crosses the wire or the coordinates of the wire centre.
centreflag whether to return the coordinates of the line-wire crossing point or of the wire centre.
rcradius [cm] of the wire.

Reimplemented from Garfield::ComponentBase.

Definition at line 407 of file ComponentNeBem2d.cc.

411 {
412 xc = x0;
413 yc = y0;
414 zc = z0;
415
416 if (m_wires.empty()) return false;
417
418 const double dx = x1 - x0;
419 const double dy = y1 - y0;
420 const double d2 = dx * dx + dy * dy;
421 // Make sure the step length is non-zero.
422 if (d2 < Small) return false;
423 const double invd2 = 1. / d2;
424
425 // Both coordinates are assumed to be located inside
426 // the drift area and inside a drift medium.
427 // This should have been checked before this call.
428
429 double dMin2 = 0.;
430 for (const auto& wire : m_wires) {
431 const double xw = wire.x;
432 const double yw = wire.y;
433 // Calculate the smallest distance between track and wire.
434 const double xIn0 = dx * (xw - x0) + dy * (yw - y0);
435 // Check if the minimum is located before (x0, y0).
436 if (xIn0 < 0.) continue;
437 const double xIn1 = -(dx * (xw - x1) + dy * (yw - y1));
438 // Check if the minimum is located behind (x1, y1).
439 if (xIn1 < 0.) continue;
440 // Minimum is located between (x0, y0) and (x1, y1).
441 const double xw0 = xw - x0;
442 const double xw1 = xw - x1;
443 const double yw0 = yw - y0;
444 const double yw1 = yw - y1;
445 const double dw02 = xw0 * xw0 + yw0 * yw0;
446 const double dw12 = xw1 * xw1 + yw1 * yw1;
447 if (xIn1 * xIn1 * dw02 > xIn0 * xIn0 * dw12) {
448 dMin2 = dw02 - xIn0 * xIn0 * invd2;
449 } else {
450 dMin2 = dw12 - xIn1 * xIn1 * invd2;
451 }
452 const double r2 = wire.r * wire.r;
453 if (dMin2 < r2) {
454 // Wire has been crossed.
455 if (centre) {
456 xc = xw;
457 yc = yw;
458 } else {
459 // Find the point of intersection.
460 const double p = -xIn0 * invd2;
461 const double q = (dw02 - r2) * invd2;
462 const double s = sqrt(p * p - q);
463 const double t = std::min(-p + s, -p - s);
464 xc = x0 + t * dx;
465 yc = y0 + t * dy;
466 zc = z0 + t * (z1 - z0);
467 }
468 rc = wire.r;
469 return true;
470 }
471 }
472 return false;
473}
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314

◆ SetMaxNumberOfIterations()

void Garfield::ComponentNeBem2d::SetMaxNumberOfIterations ( const unsigned int  niter)

Definition at line 668 of file ComponentNeBem2d.cc.

668 {
669 if (niter == 0) {
670 std::cerr << m_className << "::SetMaxNumberOfIterations:\n"
671 << " Number of iterations must be greater than zero.\n";
672 return;
673 }
674 m_nMaxIterations = niter;
675}

◆ SetMedium()

void Garfield::ComponentNeBem2d::SetMedium ( Medium *  medium)
inline

Set the "background" medium.

Definition at line 38 of file ComponentNeBem2d.hh.

38{ m_medium = medium; }

◆ SetNumberOfCollocationPoints()

void Garfield::ComponentNeBem2d::SetNumberOfCollocationPoints ( const unsigned int  ncoll)

Definition at line 657 of file ComponentNeBem2d.cc.

657 {
658 if (ncoll == 0) {
659 std::cerr << m_className << "::SetNumberOfCollocationPoints:\n"
660 << " Number of points must be greater than zero.\n";
661 return;
662 }
663
664 m_nCollocationPoints = ncoll;
665 m_ready = false;
666}

◆ SetNumberOfDivisions()

void Garfield::ComponentNeBem2d::SetNumberOfDivisions ( const unsigned int  ndiv)

Set the default number of elements per segment.

Definition at line 646 of file ComponentNeBem2d.cc.

646 {
647 if (ndiv == 0) {
648 std::cerr << m_className << "::SetNumberOfDivisions:\n"
649 << " Number of divisions must be greater than zero.\n";
650 return;
651 }
652
653 m_nDivisions = ndiv;
654 m_ready = false;
655}

◆ SetRangeZ()

void Garfield::ComponentNeBem2d::SetRangeZ ( const double  zmin,
const double  zmax 
)

Set the extent of the drift region along z.

Definition at line 495 of file ComponentNeBem2d.cc.

495 {
496
497 if (fabs(zmax - zmin) <= 0.) {
498 std::cerr << m_className << "::SetRangeZ: Zero range is not permitted.\n";
499 return;
500 }
501 m_zmin = std::min(zmin, zmax);
502 m_zmax = std::max(zmin, zmax);
503 m_useRangeZ = true;
504}
The documentation for this class was generated from the following files: