Garfield++ 4.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 GetElementaryCell (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax) override
 Get the coordinates of the elementary cell.
 
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 AddChargeDistribution (const double x, const double y, const double a, const double b, const double rho)
 
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::Component
 Component ()=delete
 Default constructor.
 
 Component (const std::string &name)
 Constructor.
 
virtual ~Component ()
 Destructor.
 
virtual void SetGeometry (Geometry *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 double DelayedWeightingPotential (const double x, const double y, const double z, const double t, 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.
 
virtual bool GetElementaryCell (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
 Get the coordinates of the elementary cell.
 
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 IntegrateFluxParallelogram (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)
 
double IntegrateWeightingFluxParallelogram (const std::string &label, 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)
 
double IntegrateFluxLine (const double x0, const double y0, const double z0, const double x1, const double y1, const double z1, const double xp, const double yp, const double zp, const unsigned int nI, const int isign=0)
 
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 EnablePeriodicityY (const bool on=true)
 Enable simple periodicity in the $y$ direction.
 
void EnablePeriodicityZ (const bool on=true)
 Enable simple periodicity in the $z$ direction.
 
void IsPeriodic (bool &perx, bool &pery, bool &perz)
 Return periodicity flags.
 
void EnableMirrorPeriodicityX (const bool on=true)
 Enable mirror periodicity in the $x$ direction.
 
void EnableMirrorPeriodicityY (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void EnableMirrorPeriodicityZ (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void IsMirrorPeriodic (bool &perx, bool &pery, bool &perz)
 Return mirror periodicity flags.
 
void EnableAxialPeriodicityX (const bool on=true)
 Enable axial periodicity in the $x$ direction.
 
void EnableAxialPeriodicityY (const bool on=true)
 Enable axial periodicity in the $y$ direction.
 
void EnableAxialPeriodicityZ (const bool on=true)
 Enable axial periodicity in the $z$ direction.
 
void IsAxiallyPeriodic (bool &perx, bool &pery, bool &perz)
 Return axial periodicity flags.
 
void EnableRotationSymmetryX (const bool on=true)
 Enable rotation symmetry around the $x$ axis.
 
void EnableRotationSymmetryY (const bool on=true)
 Enable rotation symmetry around the $y$ axis.
 
void EnableRotationSymmetryZ (const bool on=true)
 Enable rotation symmetry around the $z$ axis.
 
void IsRotationSymmetric (bool &rotx, bool &roty, bool &rotz)
 Return rotation symmetry flags.
 
void EnableDebugging ()
 Switch on debugging messages.
 
void DisableDebugging ()
 Switch off debugging messages.
 
virtual bool HasAttachmentMap () const
 Does the component have attachment maps?
 
virtual bool HasVelocityMap () const
 Does the component have velocity maps?
 
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 bool ElectronVelocity (const double, const double, const double, double &vx, double &vy, double &vz)
 Get the electron drift velocity.
 
virtual bool HoleVelocity (const double, const double, const double, double &vx, double &vy, double &vz)
 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::Component
std::string m_className = "Component"
 Class name.
 
Geometry * m_geometry = nullptr
 Pointer to the geometry.
 
std::array< double, 3 > m_b0 = {{0., 0., 0.}}
 Constant magnetic field.
 
bool m_ready = false
 Ready for use?
 
bool m_debug = false
 Switch on/off debugging messages.
 
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.
 

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 211 of file ComponentNeBem2d.cc.

211: Component("NeBem2d") {}
Garfield::Component::Component
Component()=delete
Default constructor.

◆ ~ComponentNeBem2d()

Garfield::ComponentNeBem2d::~ComponentNeBem2d ( )
inline

Destructor.

Definition at line 15 of file ComponentNeBem2d.hh.

15{}

Member Function Documentation

◆ AddChargeDistribution()

void Garfield::ComponentNeBem2d::AddChargeDistribution ( const double  x,
const double  y,
const double  a,
const double  b,
const double  rho 
)

Definition at line 664 of file ComponentNeBem2d.cc.

666 {
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}
Garfield::Component::m_className
std::string m_className
Class name.
Definition: Component.hh:329

◆ 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 584 of file ComponentNeBem2d.cc.

587 {
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}
Garfield::Component::m_debug
bool m_debug
Switch on/off debugging messages.
Definition: Component.hh:341
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:274

◆ 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 524 of file ComponentNeBem2d.cc.

526 {
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}
Garfield::Component::m_ready
bool m_ready
Ready for use?
Definition: Component.hh:338

◆ 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 553 of file ComponentNeBem2d.cc.

554 {
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}

◆ 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::Component.

Definition at line 213 of file ComponentNeBem2d.cc.

216 {
217 status = Field(x, y, z, ex, ey, ez, v, m, true);
218}

◆ 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::Component.

Definition at line 220 of file ComponentNeBem2d.cc.

222 {
223 double v = 0.;
224 status = Field(x, y, z, ex, ey, ez, v, m, false);
225}

◆ EnableAutoResizing()

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

Definition at line 82 of file ComponentNeBem2d.hh.

82{ m_autoSize = on; }

◆ EnableRandomCollocation()

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

Definition at line 83 of file ComponentNeBem2d.hh.

83 {
84 m_randomCollocation = on;
85 }

◆ 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::Component.

Definition at line 364 of file ComponentNeBem2d.cc.

366 {
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}
Garfield::ComponentNeBem2d::GetElementaryCell
bool GetElementaryCell(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax) override
Get the coordinates of the elementary cell.
Definition: ComponentNeBem2d.cc:374
Garfield::Component::m_geometry
Geometry * m_geometry
Pointer to the geometry.
Definition: Component.hh:332
Garfield::Geometry::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 760 of file ComponentNeBem2d.cc.

761 {
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}

◆ GetElementaryCell()

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

Get the coordinates of the elementary cell.

Reimplemented from Garfield::Component.

Definition at line 374 of file ComponentNeBem2d.cc.

376 {
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}

Referenced by GetBoundingBox().

◆ 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::Component.

Definition at line 352 of file ComponentNeBem2d.cc.

353 {
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}
Garfield::Geometry::GetMedium
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.
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:187

◆ GetNumberOfElements()

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

Return the number of boundary elements.

Definition at line 105 of file ComponentNeBem2d.hh.

105{ return m_elements.size(); }

◆ GetNumberOfRegions()

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

Return the number of regions.

Definition at line 89 of file ComponentNeBem2d.hh.

89{ return m_regions.size(); }

◆ GetNumberOfSegments()

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

Return the number of conducting straight-line segments.

Definition at line 95 of file ComponentNeBem2d.hh.

95{ return m_segments.size(); }

◆ GetNumberOfWires()

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

Return the number of wires.

Definition at line 100 of file ComponentNeBem2d.hh.

100{ 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 716 of file ComponentNeBem2d.cc.

720 {
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}
Garfield::ComponentNeBem2d::Initialise
bool Initialise()
Discretise the geometry and compute the solution.
Definition: ComponentNeBem2d.cc:775

◆ 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 734 of file ComponentNeBem2d.cc.

735 {
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}

◆ GetVoltageRange()

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

Calculate the voltage range [V].

Implements Garfield::Component.

Definition at line 316 of file ComponentNeBem2d.cc.

316 {
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}

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 747 of file ComponentNeBem2d.cc.

748 {
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}

◆ Initialise()

bool Garfield::ComponentNeBem2d::Initialise ( )

Discretise the geometry and compute the solution.

Definition at line 775 of file ComponentNeBem2d.cc.

775 {
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}
Garfield::ComponentNeBem2d::GetVoltageRange
bool GetVoltageRange(double &vmin, double &vmax) override
Calculate the voltage range [V].
Definition: ComponentNeBem2d.cc:316
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
Solve
int Solve(void)
Definition: neBEM.c:2441
InvertMatrix
int InvertMatrix(void)
Definition: neBEM.c:1236

Referenced by 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::Component.

Definition at line 493 of file ComponentNeBem2d.cc.

495 {
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}

◆ 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::Component.

Definition at line 425 of file ComponentNeBem2d.cc.

429 {
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}
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314

◆ SetMaxNumberOfIterations()

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

Definition at line 707 of file ComponentNeBem2d.cc.

707 {
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}

◆ SetMedium()

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

Set the "background" medium.

Definition at line 40 of file ComponentNeBem2d.hh.

40{ m_medium = medium; }

◆ SetNumberOfCollocationPoints()

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

Definition at line 696 of file ComponentNeBem2d.cc.

696 {
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}

◆ SetNumberOfDivisions()

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

Set the default number of elements per segment.

Definition at line 685 of file ComponentNeBem2d.cc.

685 {
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}

◆ SetRangeZ()

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

Set the extent of the drift region along z.

Definition at line 513 of file ComponentNeBem2d.cc.

513 {
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}
The documentation for this class was generated from the following files: