Garfield++ 5.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
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.
 
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.
 
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 CrossedWire (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 InTrapRadius (const double q0, const double x0, const double y0, const double z0, double &xw, double &yx, double &rw) override
 
std::array< double, 3 > ElectricField (const double x, const double y, const double z)
 Calculate the drift field [V/cm] at (x, y, z).
 
- 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.
 
std::array< double, 3 > ElectricField (const double x, const double y, const double z)
 Calculate the drift field [V/cm] at (x, y, z).
 
virtual double ElectricPotential (const double x, const double y, const double z)
 Calculate the (drift) electrostatic potential [V] at (x, y, z).
 
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?
 
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 CrossedPlane (const double x0, const double y0, const double z0, const double x1, const double y1, const double z1, double &xc, double &yc, double &zc)
 
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 HasMagneticField () const
 Does the component have a non-zero magnetic field?
 
virtual bool HasTownsendMap () const
 Does the component have maps of the Townsend coefficient?
 
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 ElectronTownsend (const double, const double, const double, double &alpha)
 Get the electron Townsend coefficient.
 
virtual bool HoleTownsend (const double, const double, const double, double &alpha)
 Get the hole Townsend 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)
 
virtual double StepSizeHint ()
 

Additional Inherited Members

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

665 {
666 if (a < Small || b < Small) {
667 std::cerr << m_className << "::AddChargeDistribution:\n"
668 << " Lengths must be > 0.\n";
669 return;
670 }
671 const double a2 = a * a;
672 const double b2 = b * b;
673 const double v0 = -2 * (Pi * b2 + 2 * atan(b / a) * (a2 - b2));
674 SpaceCharge box;
675 box.x = x;
676 box.y = y;
677 box.a = a;
678 box.b = b;
679 box.q = rho;
680 box.v0 = v0;
681 m_spaceCharge.push_back(std::move(box));
682}
Garfield::Component::m_className
std::string m_className
Class name.
Definition Component.hh:359

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

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

525 {
526 const double dx = x1 - x0;
527 const double dy = y1 - y0;
528 if (dx * dx + dy * dy < Small) {
529 std::cerr << m_className << "::AddSegment: Length must be > 0.\n";
530 return false;
531 }
532
533 Segment segment;
534 segment.x0 = {x0, y0};
535 segment.x1 = {x1, y1};
536 segment.bc = std::make_pair(Voltage, v);
537 segment.region1 = -1;
538 segment.region2 = -1;
539 segment.ndiv = ndiv;
540 m_segments.push_back(std::move(segment));
541
542 if (m_debug) {
543 std::cout << m_className << "::AddSegment:\n ("
544 << x0 << ", " << y0 << ") - (" << x1 << ", " << y1 << ")\n"
545 << " Potential: " << v << " V\n";
546 }
547
548 m_ready = false;
549 return true;
550}
Garfield::Component::m_ready
bool m_ready
Ready for use?
Definition Component.hh:368

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

553 {
554 if (d < Small) {
555 std::cerr << m_className << "::AddWire: Diameter must be > 0.\n";
556 return false;
557 }
558 if (ntrap <= 0) {
559 std::cerr << m_className << "::AddWire: Nbr. of trap radii must be > 0.\n";
560 return false;
561 }
562
563 Wire wire;
564 wire.x = x;
565 wire.y = y;
566 wire.r = 0.5 * d;
567 wire.v = v;
568 wire.q = 0.;
569 wire.ntrap = ntrap;
570 m_wires.push_back(std::move(wire));
571
572 if (m_debug) {
573 std::cout << m_className << "::AddWire:\n"
574 << " Centre: (" << x << ", " << y << ")\n"
575 << " Diameter: " << d << " cm\n"
576 << " Potential: " << v << " V\n";
577 }
578
579 m_ready = false;
580 return true;
581}

◆ CrossedWire()

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

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

◆ ElectricField() [1/3]

std::array< double, 3 > Garfield::Component::ElectricField ( const double x,
const double y,
const double z )

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

Definition at line 55 of file Component.cc.

43 {
44 double ex = 0., ey = 0., ez = 0.;
45 Medium* medium = nullptr;
46 int status = 0;
47 ElectricField(x, y, z, ex, ey, ez, medium, status);
48 std::array<double, 3> efield = {ex, ey, ez};
49 return efield;
50}
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:220

◆ ElectricField() [2/3]

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() [3/3]

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 60 of file ComponentNeBem2d.hh.

60{ m_autoSize = on; }

◆ EnableRandomCollocation()

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

Definition at line 61 of file ComponentNeBem2d.hh.

61 {
62 m_randomCollocation = on;
63 }

◆ 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:362

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

760 {
761
762 if (i >= m_elements.size()) return false;
763 const auto& element = m_elements[i];
764 ToGlobal(-element.a, 0., element.cphi, element.sphi, x0, y0);
765 ToGlobal( element.a, 0., element.cphi, element.sphi, x1, y1);
766 x0 += element.x;
767 y0 += element.y;
768 x1 += element.x;
769 y1 += element.y;
770 q = element.q;
771 return true;
772}

◆ 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::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 83 of file ComponentNeBem2d.hh.

83{ return m_elements.size(); }

◆ GetNumberOfRegions()

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

Return the number of regions.

Definition at line 67 of file ComponentNeBem2d.hh.

67{ return m_regions.size(); }

◆ GetNumberOfSegments()

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

Return the number of conducting straight-line segments.

Definition at line 73 of file ComponentNeBem2d.hh.

73{ return m_segments.size(); }

◆ GetNumberOfWires()

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

Return the number of wires.

Definition at line 78 of file ComponentNeBem2d.hh.

78{ 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 715 of file ComponentNeBem2d.cc.

719 {
720 if (i >= m_regions.size()) return false;
721 if (!m_ready) {
722 if (!Initialise()) return false;
723 }
724 const auto& region = m_regions[i];
725 xv = region.xv;
726 yv = region.yv;
727 medium = region.medium;
728 bctype = region.bc.first == Voltage ? 1 : 4;
729 v = region.bc.second;
730 return true;
731}
Garfield::ComponentNeBem2d::Initialise
bool Initialise()
Discretise the geometry and compute the solution.
Definition ComponentNeBem2d.cc:774

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

734 {
735
736 if (i >= m_segments.size()) return false;
737 const auto& seg = m_segments[i];
738 x0 = seg.x0[0];
739 y0 = seg.x0[1];
740 x1 = seg.x1[0];
741 y1 = seg.x1[1];
742 v = seg.bc.second;
743 return true;
744}

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

747 {
748
749 if (i >= m_wires.size()) return false;
750 const auto& wire = m_wires[i];
751 x = wire.x;
752 y = wire.y;
753 d = 2 * wire.r;
754 v = wire.v;
755 q = wire.q;
756 return true;
757}

◆ Initialise()

bool Garfield::ComponentNeBem2d::Initialise ( )

Discretise the geometry and compute the solution.

Definition at line 774 of file ComponentNeBem2d.cc.

774 {
775
776 m_ready = false;
777 m_elements.clear();
778
779 double vmin = 0., vmax = 0.;
780 if (!GetVoltageRange(vmin, vmax)) {
781 std::cerr << m_className << "::Initialise:\n"
782 << " Could not determine the voltage range.\n";
783 return false;
784 }
785 if (fabs(vmin - vmax) < 1.e-6 * (vmin + vmax)) {
786 std::cerr << m_className << "::Initialise:\n"
787 << " All potentials are the same.\n";
788 return false;
789 }
790
791 if (m_debug) std::cout << m_className << "::Initialise:\n";
792 // Loop over the regions.
793 const unsigned int nRegions = m_regions.size();
794 if (m_debug) std::cout << " " << nRegions << " regions.\n";
795 std::vector<std::vector<unsigned int> > motherRegions(nRegions);
796 for (unsigned int i = 0; i < nRegions; ++i) {
797 auto& region = m_regions[i];
798 // Check if the region is fully enclosed by other ones.
799 for (unsigned int j = 0; j < nRegions; ++j) {
800 if (i == j) continue;
801 const auto& other = m_regions[j];
802 if (!Enclosed(region.xv, region.yv, other.xv, other.yv)) continue;
803 motherRegions[i].push_back(j);
804 if (m_debug) {
805 std::cout << " Region " << i << " is enclosed by region "
806 << j << ".\n";
807 }
808 }
809 region.depth = motherRegions[i].size();
810 }
811
812 std::vector<Segment> segments;
813 for (unsigned int i = 0; i < nRegions; ++i) {
814 const auto& region = m_regions[i];
815 int outerRegion = -1;
816 for (const unsigned int k : motherRegions[i]) {
817 if (outerRegion < 0) {
818 outerRegion = k;
819 } else if (m_regions[outerRegion].depth < m_regions[k].depth) {
820 outerRegion = k;
821 }
822 }
823 // Add the segments bounding this region.
824 const unsigned int n = region.xv.size();
825 for (unsigned int j = 0; j < n; ++j) {
826 const unsigned int k = j < n - 1 ? j + 1 : 0;
827 Segment seg;
828 seg.x0 = {region.xv[j], region.yv[j]};
829 seg.x1 = {region.xv[k], region.yv[k]};
830 seg.region1 = i;
831 seg.region2 = outerRegion;
832 seg.bc = region.bc;
833 seg.ndiv = region.ndiv;
834 segments.push_back(std::move(seg));
835 }
836 }
837 // Add the segments specified by the user.
838 segments.insert(segments.end(), m_segments.begin(), m_segments.end());
839 const unsigned int nSegments = segments.size();
840 if (m_debug) std::cout << " " << nSegments << " segments.\n";
841 std::vector<bool> done(nSegments, false);
842 // Look for overlaps.
843 for (unsigned int i = 0; i < nSegments; ++i) {
844 if (done[i]) continue;
845 if (m_debug) {
846 std::cout << " Segment " << i << ". ("
847 << segments[i].x0[0] << ", " << segments[i].x0[1] << ") - ("
848 << segments[i].x1[0] << ", " << segments[i].x1[1] << ")\n";
849 }
850 const double x0 = segments[i].x0[0];
851 const double x1 = segments[i].x1[0];
852 const double y0 = segments[i].x0[1];
853 const double y1 = segments[i].x1[1];
854 // Pick up all collinear segments.
855 std::vector<Segment> newSegments;
856 for (unsigned int j = i + 1; j < nSegments; ++j) {
857 const double u0 = segments[j].x0[0];
858 const double u1 = segments[j].x1[0];
859 const double v0 = segments[j].x0[1];
860 const double v1 = segments[j].x1[1];
861 const double epsx = std::max(1.e-10 * std::max({fabs(x0), fabs(x1),
862 fabs(u0), fabs(u1)}),
863 1.e-10);
864 const double epsy = std::max(1.e-10 * std::max({fabs(y0), fabs(y1),
865 fabs(v0), fabs(v1)}),
866 1.e-10);
867 const double a00 = y1 - y0;
868 const double a01 = v1 - v0;
869 const double a10 = x0 - x1;
870 const double a11 = u0 - u1;
871 const double det = a00 * a11 - a10 * a01;
872 const double tol = epsx * epsy;
873 // Skip non-parallel segments.
874 if (std::abs(det) > tol) continue;
875
876 if (std::abs(x0 * (y1 - v0) + x1 * (v0 - y0) + u0 * (y0 - y1)) > tol) {
877 continue;
878 }
879 newSegments.push_back(segments[j]);
880 done[j] = true;
881 }
882 newSegments.push_back(segments[i]);
883 if (newSegments.size() > 1) {
884 if (m_debug) {
885 std::cout << " Determining overlaps of " << newSegments.size()
886 << " collinear segments.\n";
887 }
888 EliminateOverlaps(newSegments);
889 if (m_debug) {
890 std::cout << " " << newSegments.size()
891 << " segments after splitting/merging.\n";
892 }
893 }
894 for (const auto& segment : newSegments) {
895 double lambda = 0.;
896 if (segment.bc.first == Dielectric) {
897 // Dielectric-dielectric interface.
898 const int reg1 = segment.region1;
899 const int reg2 = segment.region2;
900 double eps1 = 1.;
901 if (reg1 < 0) {
902 if (m_medium) eps1 = m_medium->GetDielectricConstant();
903 } else if (m_regions[reg1].medium) {
904 eps1 = m_regions[reg1].medium->GetDielectricConstant();
905 }
906 double eps2 = 1.;
907 if (reg2 < 0) {
908 if (m_medium) eps2 = m_medium->GetDielectricConstant();
909 } else if (m_regions[reg2].medium) {
910 eps2 = m_regions[reg2].medium->GetDielectricConstant();
911 }
912 if (fabs(eps1 - eps2) < 1.e-6 * (1. + fabs(eps1) + fabs(eps2))) {
913 if (m_debug) std::cout << " Same epsilon. Skip.\n";
914 continue;
915 }
916 // Compute lambda.
917 lambda = (eps2 - eps1) / (eps1 + eps2);
918 if (m_debug) std::cout << " Lambda = " << lambda << "\n";
919 }
920 const int ndiv = segment.ndiv <= 0 ? m_nDivisions : segment.ndiv;
921 Discretise(segment, m_elements, lambda, ndiv);
922 }
923 }
924 std::vector<std::vector<double> > influenceMatrix;
925 std::vector<std::vector<double> > inverseMatrix;
926
927 bool converged = false;
928 unsigned int nIter = 0;
929 while (!converged) {
930 ++nIter;
931 if (m_autoSize) {
932 std::cout << m_className << "::Initialise: Iteration " << nIter << "\n";
933 }
934 const unsigned int nElements = m_elements.size();
935 const unsigned int nEntries = nElements + m_wires.size() + 1;
936 if (m_debug) {
937 std::cout << " " << nElements << " elements.\n"
938 << " Matrix has " << nEntries << " rows/columns.\n";
939 }
940 // Compute the influence matrix.
941 influenceMatrix.assign(nEntries, std::vector<double>(nEntries, 0.));
942 if (!ComputeInfluenceMatrix(influenceMatrix)) {
943 std::cerr << m_className << "::Initialise:\n"
944 << " Error computing the influence matrix.\n";
945 return false;
946 }
947
948 // Invert the influence matrix.
949 inverseMatrix.assign(nEntries, std::vector<double>(nEntries, 0.));
950 if (!InvertMatrix(influenceMatrix, inverseMatrix)) {
951 std::cerr << m_className << "::Initialise: Matrix inversion failed.\n";
952 return false;
953 }
954 if (m_debug) std::cout << " Matrix inversion ok.\n";
955
956 // Compute the right hand side vector (boundary conditions).
957 std::vector<double> boundaryConditions(nEntries, 0.);
958 for (unsigned int i = 0; i < nElements; ++i) {
959 if (m_elements[i].bc.first == Voltage) {
960 boundaryConditions[i] = m_elements[i].bc.second;
961 for (const auto& box : m_spaceCharge) {
962 const double x = m_elements[i].x - box.x;
963 const double y = m_elements[i].y - box.y;
964 const double vs = BoxPotential(box.a, box.b, x, y, box.v0) * box.q;
965 boundaryConditions[i] -= vs;
966 }
967 } else {
968 for (const auto& box : m_spaceCharge) {
969 const double x = m_elements[i].x - box.x;
970 const double y = m_elements[i].y - box.y;
971 double fx = 0., fy = 0.;
972 BoxField(box.a, box.b, x, y, fx, fy);
973 // Rotate to the local frame of the target element.
974 ToLocal(fx, fy, m_elements[i].cphi, m_elements[i].sphi, fx, fy);
975 boundaryConditions[i] -= box.q * fy;
976 }
977 }
978 }
979 const unsigned int nWires = m_wires.size();
980 for (unsigned int i = 0; i < nWires; ++i) {
981 boundaryConditions[nElements + i] = m_wires[i].v;
982 for (const auto& box : m_spaceCharge) {
983 const double x = m_wires[i].x - box.x;
984 const double y = m_wires[i].y - box.y;
985 const double vs = BoxPotential(box.a, box.b, x, y, box.v0) * box.q;
986 boundaryConditions[nElements + i] -= vs;
987 }
988 }
989
990 double qsum = 0.;
991 for (const auto& box : m_spaceCharge) {
992 qsum += 4 * box.q * box.a * box.b;
993 }
994 boundaryConditions.back() = -qsum;
995
996 // Solve for the charge distribution.
997 if (!Solve(inverseMatrix, boundaryConditions)) {
998 std::cerr << m_className << "::Initialise: Solution failed.\n";
999 return false;
1000 }
1001 if (m_debug) std::cout << " Solution ok.\n";
1002 const double tol = 1.e-6 * fabs(vmax - vmin);
1003 std::vector<bool> ok(nElements, true);
1004 converged = CheckConvergence(tol, ok);
1005 if (!m_autoSize) break;
1006 if (nIter >= m_nMaxIterations) break;
1007 for (unsigned int j = 0; j < nElements; ++j) {
1008 if (!ok[j]) {
1009 SplitElement(m_elements[j], m_elements);
1010 if (m_debug) std::cout << " Splitting element " << j << ".\n";
1011 }
1012 }
1013 }
1014 // Sort the regions by depth (innermost first).
1015 std::sort(m_regions.begin(), m_regions.end(),
1016 [](const Region& lhs, const Region& rhs) {
1017 return (lhs.depth > rhs.depth);
1018 });
1019 m_ready = true;
1020 return true;
1021}
Garfield::ComponentNeBem2d::GetVoltageRange
bool GetVoltageRange(double &vmin, double &vmax) override
Calculate the voltage range [V].
Definition ComponentNeBem2d.cc:316
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition DoubleAc.h:615
Solve
int Solve(void)
Definition neBEM.c:2784
InvertMatrix
int InvertMatrix(void)
Definition neBEM.c:1309

Referenced by GetRegion().

◆ InTrapRadius()

bool Garfield::ComponentNeBem2d::InTrapRadius ( 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 492 of file ComponentNeBem2d.cc.

494 {
495
496 for (const auto& wire : m_wires) {
497 // Skip wires with the wrong charge.
498 if (q * wire.q > 0.) continue;
499 const double dx = wire.x - x;
500 const double dy = wire.y - y;
501 const double rt = wire.r * wire.ntrap;
502 if (dx * dx + dy * dy < rt * rt) {
503 xw = wire.x;
504 yw = wire.y;
505 rw = wire.r;
506 return true;
507 }
508 }
509 return false;
510}

◆ SetMaxNumberOfIterations()

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

Definition at line 706 of file ComponentNeBem2d.cc.

706 {
707 if (niter == 0) {
708 std::cerr << m_className << "::SetMaxNumberOfIterations:\n"
709 << " Number of iterations must be greater than zero.\n";
710 return;
711 }
712 m_nMaxIterations = niter;
713}

◆ SetMedium()

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

Set the "background" medium.

Definition at line 18 of file ComponentNeBem2d.hh.

18{ m_medium = medium; }

◆ SetNumberOfCollocationPoints()

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

Definition at line 695 of file ComponentNeBem2d.cc.

695 {
696 if (ncoll == 0) {
697 std::cerr << m_className << "::SetNumberOfCollocationPoints:\n"
698 << " Number of points must be greater than zero.\n";
699 return;
700 }
701
702 m_nCollocationPoints = ncoll;
703 m_ready = false;
704}

◆ SetNumberOfDivisions()

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

Set the default number of elements per segment.

Definition at line 684 of file ComponentNeBem2d.cc.

684 {
685 if (ndiv == 0) {
686 std::cerr << m_className << "::SetNumberOfDivisions:\n"
687 << " Number of divisions must be greater than zero.\n";
688 return;
689 }
690
691 m_nDivisions = ndiv;
692 m_ready = false;
693}

◆ SetRangeZ()

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

Set the extent of the drift region along z.

Definition at line 512 of file ComponentNeBem2d.cc.

512 {
513
514 if (fabs(zmax - zmin) <= 0.) {
515 std::cerr << m_className << "::SetRangeZ: Zero range is not permitted.\n";
516 return;
517 }
518 m_zmin = std::min(zmin, zmax);
519 m_zmax = std::max(zmin, zmax);
520 m_useRangeZ = true;
521}
The documentation for this class was generated from the following files: