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 | Protected Member Functions | List of all members
Garfield::ComponentNeBem3d Class Reference

Interface to neBEM. More...

#include <ComponentNeBem3d.hh>

+ Inheritance diagram for Garfield::ComponentNeBem3d:

Public Member Functions

 ComponentNeBem3d ()
 Constructor.
 
 ~ComponentNeBem3d ()
 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].
 
void WeightingField (const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label) override
 
double WeightingPotential (const double x, const double y, const double z, const std::string &label) override
 
void AddPlaneX (const double x, const double voltage)
 Add a plane at constant x.
 
void AddPlaneY (const double y, const double voltage)
 Add a plane at constant y.
 
void AddPlaneZ (const double z, const double voltage)
 Add a plane at constant z.
 
unsigned int GetNumberOfPlanesX () const
 Get the number of equipotential planes at constant x.
 
unsigned int GetNumberOfPlanesY () const
 Get the number of equipotential planes at constant y.
 
unsigned int GetNumberOfPlanesZ () const
 Get the number of equipotential planes at constant z.
 
bool GetPlaneX (const unsigned int i, double &x, double &v) const
 Retrieve the parameters of a plane at constant x.
 
bool GetPlaneY (const unsigned int i, double &y, double &v) const
 Retrieve the parameters of a plane at constant y.
 
bool GetPlaneZ (const unsigned int i, double &z, double &v) const
 Retrieve the parameters of a plane at constant z.
 
unsigned int GetNumberOfPrimitives () const
 
bool GetPrimitive (const unsigned int i, double &a, double &b, double &c, std::vector< double > &xv, std::vector< double > &yv, std::vector< double > &zv, int &interface, double &v, double &q, double &lambda) const
 
bool GetPrimitive (const unsigned int i, double &a, double &b, double &c, std::vector< double > &xv, std::vector< double > &yv, std::vector< double > &zv, int &vol1, int &vol2) const
 
bool GetVolume (const unsigned int vol, int &shape, int &material, double &eps, double &potential, double &charge, int &bc)
 
int GetVolume (const double x, const double y, const double z)
 
unsigned int GetNumberOfElements () const
 
bool GetElement (const unsigned int i, std::vector< double > &xv, std::vector< double > &yv, std::vector< double > &zv, int &interface, double &bc, double &lambda) const
 
bool Initialise ()
 
void SetTargetElementSize (const double length)
 
void SetMinMaxNumberOfElements (const unsigned int nmin, const unsigned int nmax)
 
void UseLUInversion ()
 Invert the influence matrix using lower-upper (LU) decomposition.
 
void UseSVDInversion ()
 Invert the influence matrix using singular value decomposition.
 
void SetPeriodicCopies (const unsigned int nx, const unsigned int ny, const unsigned int nz)
 
void GetPeriodicCopies (unsigned int &nx, unsigned int &ny, unsigned int &nz) const
 Retrieve the number of periodic copies used by neBEM.
 
void SetPeriodicityX (const double s)
 Set the periodic length [cm] in the x-direction.
 
void SetPeriodicityY (const double s)
 Set the periodic length [cm] in the y-direction.
 
void SetPeriodicityZ (const double s)
 Set the periodic length [cm] in the z-direction.
 
void SetMirrorPeriodicityX (const double s)
 Set the periodic length [cm] in the x-direction.
 
void SetMirrorPeriodicityY (const double s)
 Set the periodic length [cm] in the y-direction.
 
void SetMirrorPeriodicityZ (const double s)
 Set the periodic length [cm] in the z-direction.
 
bool GetPeriodicityX (double &s) const
 Get the periodic length in the x-direction.
 
bool GetPeriodicityY (double &s) const
 Get the periodic length in the y-direction.
 
bool GetPeriodicityZ (double &s) const
 Get the periodic length in the z-direction.
 
void SetNumberOfThreads (const unsigned int n)
 Set the number of threads to be used by neBEM.
 
- 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)
 

Protected Member Functions

void Reset () override
 Reset the component.
 
void UpdatePeriodicity () override
 Verify periodicities.
 
virtual void Reset ()=0
 Reset the component.
 
virtual void UpdatePeriodicity ()=0
 Verify periodicities.
 

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

Interface to neBEM.

Definition at line 12 of file ComponentNeBem3d.hh.

Constructor & Destructor Documentation

◆ ComponentNeBem3d()

Garfield::ComponentNeBem3d::ComponentNeBem3d ( )

Constructor.

Definition at line 384 of file ComponentNeBem3d.cc.

384: Component("NeBem3d") {}
Garfield::Component::Component
Component()=delete
Default constructor.

◆ ~ComponentNeBem3d()

Garfield::ComponentNeBem3d::~ComponentNeBem3d ( )
inline

Destructor.

Definition at line 17 of file ComponentNeBem3d.hh.

17{}

Member Function Documentation

◆ AddPlaneX()

void Garfield::ComponentNeBem3d::AddPlaneX ( const double  x,
const double  voltage 
)

Add a plane at constant x.

Definition at line 480 of file ComponentNeBem3d.cc.

480 {
481 if (m_ynplan[0] && m_ynplan[1]) {
482 std::cerr << m_className << "::AddPlaneX:\n"
483 << " Cannot have more than two planes at constant x.\n";
484 return;
485 }
486
487 if (m_ynplan[0]) {
488 m_ynplan[1] = true;
489 if (x < m_coplan[0]) {
490 m_coplan[1] = m_coplan[0];
491 m_vtplan[1] = m_vtplan[0];
492 m_coplan[0] = x;
493 m_vtplan[0] = v;
494 } else {
495 m_coplan[1] = x;
496 m_vtplan[1] = v;
497 }
498 } else {
499 m_ynplan[0] = true;
500 m_coplan[0] = x;
501 m_vtplan[0] = v;
502 }
503 m_ready = false;
504}
Garfield::Component::m_className
std::string m_className
Class name.
Definition: Component.hh:329
Garfield::Component::m_ready
bool m_ready
Ready for use?
Definition: Component.hh:338

◆ AddPlaneY()

void Garfield::ComponentNeBem3d::AddPlaneY ( const double  y,
const double  voltage 
)

Add a plane at constant y.

Definition at line 506 of file ComponentNeBem3d.cc.

506 {
507 if (m_ynplan[2] && m_ynplan[3]) {
508 std::cerr << m_className << "::AddPlaneY:\n"
509 << " Cannot have more than two planes at constant y.\n";
510 return;
511 }
512
513 if (m_ynplan[2]) {
514 m_ynplan[3] = true;
515 if (y < m_coplan[2]) {
516 m_coplan[3] = m_coplan[2];
517 m_vtplan[3] = m_vtplan[2];
518 m_coplan[2] = y;
519 m_vtplan[2] = v;
520 } else {
521 m_coplan[3] = y;
522 m_vtplan[3] = v;
523 }
524 } else {
525 m_ynplan[2] = true;
526 m_coplan[2] = y;
527 m_vtplan[2] = v;
528 }
529 m_ready = false;
530}

◆ AddPlaneZ()

void Garfield::ComponentNeBem3d::AddPlaneZ ( const double  z,
const double  voltage 
)

Add a plane at constant z.

Definition at line 532 of file ComponentNeBem3d.cc.

532 {
533 if (m_ynplan[4] && m_ynplan[5]) {
534 std::cerr << m_className << "::AddPlaneZ:\n"
535 << " Cannot have more than two planes at constant z.\n";
536 return;
537 }
538
539 if (m_ynplan[4]) {
540 m_ynplan[5] = true;
541 if (z < m_coplan[4]) {
542 m_coplan[5] = m_coplan[4];
543 m_vtplan[5] = m_vtplan[4];
544 m_coplan[4] = z;
545 m_vtplan[4] = v;
546 } else {
547 m_coplan[5] = z;
548 m_vtplan[5] = v;
549 }
550 } else {
551 m_ynplan[4] = true;
552 m_coplan[4] = z;
553 m_vtplan[4] = v;
554 }
555 m_ready = false;
556}

◆ ElectricField() [1/2]

void Garfield::ComponentNeBem3d::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 392 of file ComponentNeBem3d.cc.

395 {
396 ex = ey = ez = v = 0.;
397 status = 0;
398 // Check if the requested point is inside a medium
399 m = GetMedium(x, y, z);
400 if (!m) {
401 status = -6;
402 } else if (!m->IsDriftable()) {
403 status = -5;
404 }
405
406 if (!m_ready) {
407 if (!Initialise()) {
408 std::cerr << m_className << "::ElectricField: Initialisation failed.\n";
409 status = -11;
410 return;
411 }
412 m_ready = true;
413 }
414
415 // Construct a point.
416 neBEM::Point3D point;
417 point.X = 0.01 * x;
418 point.Y = 0.01 * y;
419 point.Z = 0.01 * z;
420
421 // Compute the field.
422 neBEM::Vector3D field;
423 if (neBEM::neBEMField(&point, &v, &field) != 0) {
424 status = -10;
425 return;
426 }
427 ex = 0.01 * field.X;
428 ey = 0.01 * field.Y;
429 ez = 0.01 * field.Z;
430}
Garfield::ComponentNeBem3d::GetMedium
Medium * GetMedium(const double x, const double y, const double z) override
Get the medium at a given location (x, y, z).
Definition: ComponentNeBem3d.cc:386

◆ ElectricField() [2/2]

void Garfield::ComponentNeBem3d::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 432 of file ComponentNeBem3d.cc.

434 {
435 double v = 0.;
436 ElectricField(x, y, z, ex, ey, ez, v, m, status);
437}
Garfield::ComponentNeBem3d::ElectricField
void ElectricField(const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status) override
Definition: ComponentNeBem3d.cc:432

Referenced by ElectricField().

◆ GetElement()

bool Garfield::ComponentNeBem3d::GetElement ( const unsigned int  i,
std::vector< double > &  xv,
std::vector< double > &  yv,
std::vector< double > &  zv,
int &  interface,
double &  bc,
double &  lambda 
) const

Definition at line 3214 of file ComponentNeBem3d.cc.

3218 {
3219
3220 if (i >= m_elements.size()) {
3221 std::cerr << m_className << "::GetElement: Index out of range.\n";
3222 return false;
3223 }
3224 const auto& element = m_elements[i];
3225 xv = element.xv;
3226 yv = element.yv;
3227 zv = element.zv;
3228 interface = element.interface;
3229 bc = element.bc;
3230 lambda = element.lambda;
3231 return true;
3232}

◆ GetMedium()

Medium * Garfield::ComponentNeBem3d::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 386 of file ComponentNeBem3d.cc.

387 {
388 if (!m_geometry) return nullptr;
389 return m_geometry->GetMedium(x, y, z, true);
390}
Garfield::Component::m_geometry
Geometry * m_geometry
Pointer to the geometry.
Definition: Component.hh:332
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.

Referenced by ElectricField().

◆ GetNumberOfElements()

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

Definition at line 64 of file ComponentNeBem3d.hh.

64{ return m_elements.size(); }

◆ GetNumberOfPlanesX()

unsigned int Garfield::ComponentNeBem3d::GetNumberOfPlanesX ( ) const

Get the number of equipotential planes at constant x.

Definition at line 558 of file ComponentNeBem3d.cc.

558 {
559 if (m_ynplan[0] && m_ynplan[1]) {
560 return 2;
561 } else if (m_ynplan[0] || m_ynplan[1]) {
562 return 1;
563 }
564 return 0;
565}

Referenced by neBEM::neBEMGetBoundingPlanes().

◆ GetNumberOfPlanesY()

unsigned int Garfield::ComponentNeBem3d::GetNumberOfPlanesY ( ) const

Get the number of equipotential planes at constant y.

Definition at line 567 of file ComponentNeBem3d.cc.

567 {
568 if (m_ynplan[2] && m_ynplan[3]) {
569 return 2;
570 } else if (m_ynplan[2] || m_ynplan[3]) {
571 return 1;
572 }
573 return 0;
574}

Referenced by neBEM::neBEMGetBoundingPlanes().

◆ GetNumberOfPlanesZ()

unsigned int Garfield::ComponentNeBem3d::GetNumberOfPlanesZ ( ) const

Get the number of equipotential planes at constant z.

Definition at line 576 of file ComponentNeBem3d.cc.

576 {
577 if (m_ynplan[4] && m_ynplan[5]) {
578 return 2;
579 } else if (m_ynplan[4] || m_ynplan[5]) {
580 return 1;
581 }
582 return 0;
583}

Referenced by neBEM::neBEMGetBoundingPlanes().

◆ GetNumberOfPrimitives()

unsigned int Garfield::ComponentNeBem3d::GetNumberOfPrimitives ( ) const
inline

Definition at line 52 of file ComponentNeBem3d.hh.

52{ return m_primitives.size(); }

Referenced by neBEM::neBEMGetNbPrimitives().

◆ GetPeriodicCopies()

void Garfield::ComponentNeBem3d::GetPeriodicCopies ( unsigned int &  nx,
unsigned int &  ny,
unsigned int &  nz 
) const
inline

Retrieve the number of periodic copies used by neBEM.

Definition at line 93 of file ComponentNeBem3d.hh.

94 {
95 nx = m_nCopiesX;
96 ny = m_nCopiesY;
97 nz = m_nCopiesZ;
98 }

Referenced by neBEM::neBEMGetPeriodicities().

◆ GetPeriodicityX()

bool Garfield::ComponentNeBem3d::GetPeriodicityX ( double &  s) const

Get the periodic length in the x-direction.

Definition at line 720 of file ComponentNeBem3d.cc.

720 {
721 if (!m_periodic[0] && !m_mirrorPeriodic[0]) {
722 s = 0.;
723 return false;
724 }
725 s = m_periodicLength[0];
726 return true;
727}
Garfield::Component::m_mirrorPeriodic
std::array< bool, 3 > m_mirrorPeriodic
Mirror periodicity in x, y, z.
Definition: Component.hh:346
Garfield::Component::m_periodic
std::array< bool, 3 > m_periodic
Simple periodicity in x, y, z.
Definition: Component.hh:344

Referenced by neBEM::neBEMGetPeriodicities().

◆ GetPeriodicityY()

bool Garfield::ComponentNeBem3d::GetPeriodicityY ( double &  s) const

Get the periodic length in the y-direction.

Definition at line 729 of file ComponentNeBem3d.cc.

729 {
730 if (!m_periodic[1] && !m_mirrorPeriodic[1]) {
731 s = 0.;
732 return false;
733 }
734 s = m_periodicLength[1];
735 return true;
736}

Referenced by neBEM::neBEMGetPeriodicities().

◆ GetPeriodicityZ()

bool Garfield::ComponentNeBem3d::GetPeriodicityZ ( double &  s) const

Get the periodic length in the z-direction.

Definition at line 738 of file ComponentNeBem3d.cc.

738 {
739 if (!m_periodic[2] && !m_mirrorPeriodic[2]) {
740 s = 0.;
741 return false;
742 }
743 s = m_periodicLength[2];
744 return true;
745}

Referenced by neBEM::neBEMGetPeriodicities().

◆ GetPlaneX()

bool Garfield::ComponentNeBem3d::GetPlaneX ( const unsigned int  i,
double &  x,
double &  v 
) const

Retrieve the parameters of a plane at constant x.

Definition at line 585 of file ComponentNeBem3d.cc.

586 {
587 if (i >= 2 || (i == 1 && !m_ynplan[1])) {
588 std::cerr << m_className << "::GetPlaneX: Index out of range.\n";
589 return false;
590 }
591 x = m_coplan[i];
592 v = m_vtplan[i];
593 return true;
594}

Referenced by neBEM::neBEMGetBoundingPlanes().

◆ GetPlaneY()

bool Garfield::ComponentNeBem3d::GetPlaneY ( const unsigned int  i,
double &  y,
double &  v 
) const

Retrieve the parameters of a plane at constant y.

Definition at line 596 of file ComponentNeBem3d.cc.

597 {
598 if (i >= 2 || (i == 1 && !m_ynplan[3])) {
599 std::cerr << m_className << "::GetPlaneY: Index out of range.\n";
600 return false;
601 }
602 y = m_coplan[i + 2];
603 v = m_vtplan[i + 2];
604 return true;
605}

Referenced by neBEM::neBEMGetBoundingPlanes().

◆ GetPlaneZ()

bool Garfield::ComponentNeBem3d::GetPlaneZ ( const unsigned int  i,
double &  z,
double &  v 
) const

Retrieve the parameters of a plane at constant z.

Definition at line 607 of file ComponentNeBem3d.cc.

608 {
609 if (i >= 2 || (i == 1 && !m_ynplan[5])) {
610 std::cerr << m_className << "::GetPlaneZ: Index out of range.\n";
611 return false;
612 }
613 z = m_coplan[i + 4];
614 v = m_vtplan[i + 4];
615 return true;
616}

Referenced by neBEM::neBEMGetBoundingPlanes().

◆ GetPrimitive() [1/2]

bool Garfield::ComponentNeBem3d::GetPrimitive ( const unsigned int  i,
double &  a,
double &  b,
double &  c,
std::vector< double > &  xv,
std::vector< double > &  yv,
std::vector< double > &  zv,
int &  interface,
double &  v,
double &  q,
double &  lambda 
) const

Definition at line 3119 of file ComponentNeBem3d.cc.

3124 {
3125 if (i >= m_primitives.size()) {
3126 std::cerr << m_className << "::GetPrimitive: Index out of range.\n";
3127 return false;
3128 }
3129 const auto& primitive = m_primitives[i];
3130 a = primitive.a;
3131 b = primitive.b;
3132 c = primitive.c;
3133 xv = primitive.xv;
3134 yv = primitive.yv;
3135 zv = primitive.zv;
3136 interface = primitive.interface;
3137 v = primitive.v;
3138 q = primitive.q;
3139 lambda = primitive.lambda;
3140 return true;
3141}

◆ GetPrimitive() [2/2]

bool Garfield::ComponentNeBem3d::GetPrimitive ( const unsigned int  i,
double &  a,
double &  b,
double &  c,
std::vector< double > &  xv,
std::vector< double > &  yv,
std::vector< double > &  zv,
int &  vol1,
int &  vol2 
) const

Definition at line 3143 of file ComponentNeBem3d.cc.

3146 {
3147 if (i >= m_primitives.size()) {
3148 std::cerr << m_className << "::GetPrimitive: Index out of range.\n";
3149 return false;
3150 }
3151 const auto& primitive = m_primitives[i];
3152 a = primitive.a;
3153 b = primitive.b;
3154 c = primitive.c;
3155 xv = primitive.xv;
3156 yv = primitive.yv;
3157 zv = primitive.zv;
3158 vol1 = primitive.vol1;
3159 vol2 = primitive.vol2;
3160 return true;
3161}

◆ GetVoltageRange()

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

Calculate the voltage range [V].

Implements Garfield::Component.

Definition at line 439 of file ComponentNeBem3d.cc.

439 {
440 // Voltage and other bc have to come from the solids
441 vmin = vmax = 0;
442 return true;
443}

◆ GetVolume() [1/2]

int Garfield::ComponentNeBem3d::GetVolume ( const double  x,
const double  y,
const double  z 
)

Definition at line 3201 of file ComponentNeBem3d.cc.

3202 {
3203 if (!m_geometry) return -1;
3204 const unsigned int nSolids = m_geometry->GetNumberOfSolids();
3205 for (unsigned int i = 0; i < nSolids; ++i) {
3206 Medium* medium = nullptr;
3207 const auto solid = m_geometry->GetSolid(i, medium);
3208 if (!solid) continue;
3209 if (solid->IsInside(x, y, z)) return solid->GetId();
3210 }
3211 return -1;
3212}
Garfield::Geometry::GetSolid
virtual Solid * GetSolid(const size_t) const
Get a solid from the list.
Definition: Geometry.hh:29
Garfield::Geometry::GetNumberOfSolids
virtual size_t GetNumberOfSolids() const
Return the number of solids in the geometry.
Definition: Geometry.hh:27
Garfield::Solid::GetId
unsigned int GetId() const
Get the ID of the solid.
Definition: Solid.hh:137

◆ GetVolume() [2/2]

bool Garfield::ComponentNeBem3d::GetVolume ( const unsigned int  vol,
int &  shape,
int &  material,
double &  eps,
double &  potential,
double &  charge,
int &  bc 
)

Definition at line 3163 of file ComponentNeBem3d.cc.

3166 {
3167
3168 if (!m_geometry) return false;
3169 const unsigned int nSolids = m_geometry->GetNumberOfSolids();
3170 for (unsigned int i = 0; i < nSolids; ++i) {
3171 Medium* medium = nullptr;
3172 const auto solid = m_geometry->GetSolid(i, medium);
3173 if (!solid) continue;
3174 if (solid->GetId() != vol) continue;
3175 if (solid->IsTube() || solid->IsWire()) {
3176 shape = 1;
3177 } else if (solid->IsHole()) {
3178 shape = 2;
3179 } else if (solid->IsBox()) {
3180 shape = 3;
3181 } else if (solid->IsSphere()) {
3182 shape = 4;
3183 } else if (solid->IsRidge()) {
3184 shape = 5;
3185 } else if (solid->IsExtrusion()) {
3186 shape = 6;
3187 } else {
3188 std::cerr << m_className << "::GetVolume: Unknown solid shape.\n";
3189 return false;
3190 }
3191 material = medium ? medium->GetId() : 11;
3192 epsilon = medium ? medium->GetDielectricConstant() : 1.;
3193 potential = solid->GetBoundaryPotential();
3194 charge = solid->GetBoundaryChargeDensity();
3195 bc = solid->GetBoundaryConditionType();
3196 return true;
3197 }
3198 return false;
3199}

Referenced by neBEM::neBEMVolumePoint().

◆ Initialise()

bool Garfield::ComponentNeBem3d::Initialise ( )

Retrieve surface panels, remove contacts and cut polygons to rectangles and right-angle triangles.

Definition at line 747 of file ComponentNeBem3d.cc.

747 {
748 // Reset the lists.
749 m_primitives.clear();
750 m_elements.clear();
751
752 if (!m_geometry) {
753 std::cerr << m_className << "::Initialise: Geometry not set.\n";
754 return false;
755 }
756 // Be sure we won't have intersections with the bounding box.
757 // TODO! Loop over the solids and call PLACYE, PLACHE, PLABXE, PLASPE, ...
758 // Loop over the solids.
759 std::map<int, Solid*> solids;
760 std::map<int, Solid::BoundaryCondition> bc;
761 std::map<int, double> volt;
762 std::map<int, double> eps;
763 std::map<int, double> charge;
764 const unsigned int nSolids = m_geometry->GetNumberOfSolids();
765 std::vector<Panel> panelsIn;
766 for (unsigned int i = 0; i < nSolids; ++i) {
767 Medium* medium = nullptr;
768 const auto solid = m_geometry->GetSolid(i, medium);
769 if (!solid) continue;
770 // Get the panels.
771 solid->SolidPanels(panelsIn);
772 // Get the boundary condition.
773 const auto id = solid->GetId();
774 solids[id] = solid;
775 bc[id] = solid->GetBoundaryConditionType();
776 volt[id] = solid->GetBoundaryPotential();
777 if (bc[id] == Solid::Unknown) {
778 std::cout << m_className << "::Initialise:\n"
779 << " Boundary conditions for solid " << id << " not set.\n";
780 if (medium && medium->IsConductor()) {
781 std::cout << " Assuming the panels to be grounded.\n";
782 bc[id] = Solid::Voltage;
783 volt[id] = 0.;
784 } else {
785 std::cout << " Assuming dielectric-dielectric interfaces.\n";
786 bc[id] = Solid::Dielectric;
787 }
788 }
789 charge[id] = solid->GetBoundaryChargeDensity();
790 if (!medium) {
791 eps[id] = 1.;
792 } else {
793 eps[id] = medium->GetDielectricConstant();
794 }
795 }
796 // Apply cuts.
797 // CALL CELSCT('APPLY')
798 // Reduce to basic periodic copy.
799 ShiftPanels(panelsIn);
800
801 // Find contact panels and split into primitives.
802
803 // *---------------------------------------------------------------------
804 // * PLABEM - Prepares panels for BEM applications: removes the contacts
805 // * and cuts polygons to rectangles and right-angle triangles.
806 // *---------------------------------------------------------------------
807
808 // Establish tolerances.
809 const double epsang = 1.e-6; // BEMEPA
810 const double epsxyz = 1.e-6; // BEMEPD
811 // CALL EPSSET('SET',EPSXYZ,EPSXYZ,EPSXYZ)
812
813 const unsigned int nIn = panelsIn.size();
814 if (m_debug) {
815 std::cout << m_className << "::Initialise: Retrieved "
816 << nIn << " panels from the solids.\n";
817 }
818 // Keep track of which panels have been processed.
819 std::vector<bool> mark(nIn, false);
820 // Count the number of interface panels that have been discarded.
821 unsigned int nTrivial = 0;
822 unsigned int nConflicting = 0;
823 unsigned int nNotImplemented = 0;
824 // Pick up panels which coincide potentially.
825 for (unsigned int i = 0; i < nIn; ++i) {
826 // Skip panels already done.
827 if (mark[i]) continue;
828 // Fetch panel parameters.
829 const double a1 = panelsIn[i].a;
830 const double b1 = panelsIn[i].b;
831 const double c1 = panelsIn[i].c;
832 const auto& xp1 = panelsIn[i].xv;
833 const auto& yp1 = panelsIn[i].yv;
834 const auto& zp1 = panelsIn[i].zv;
835 const unsigned int np1 = xp1.size();
836 // Establish its norm and offset.
837 const double d1 = a1 * xp1[0] + b1 * yp1[0] + c1 * zp1[0];
838 if (m_debug) {
839 std::cout << " Panel " << i << "\n Norm vector: "
840 << a1 << ", " << b1 << ", " << c1 << ", " << d1 << ".\n";
841 }
842 // Rotation matrix.
843 std::array<std::array<double, 3>, 3> rot;
844 if (fabs(c1) <= fabs(a1) && fabs(c1) <= fabs(b1)) {
845 // Rotation: removing C
846 rot[0][0] = b1 / sqrt(a1 * a1 + b1 * b1);
847 rot[0][1] = -a1 / sqrt(a1 * a1 + b1 * b1);
848 rot[0][2] = 0.;
849 } else if (fabs(b1) <= fabs(a1) && fabs(b1) <= fabs(c1)) {
850 // Rotation: removing B
851 rot[0][0] = c1 / sqrt(a1 * a1 + c1 * c1);
852 rot[0][1] = 0.;
853 rot[0][2] = -a1 / sqrt(a1 * a1 + c1 * c1);
854 } else {
855 // Rotation: removing A
856 rot[0][0] = 0.;
857 rot[0][1] = c1 / sqrt(b1 * b1 + c1 * c1);
858 rot[0][2] = -b1 / sqrt(b1 * b1 + c1 * c1);
859 }
860 rot[2][0] = a1;
861 rot[2][1] = b1;
862 rot[2][2] = c1;
863 rot[1][0] = rot[2][1] * rot[0][2] - rot[2][2] * rot[0][1];
864 rot[1][1] = rot[2][2] * rot[0][0] - rot[2][0] * rot[0][2];
865 rot[1][2] = rot[2][0] * rot[0][1] - rot[2][1] * rot[0][0];
866 // Rotate to the x, y plane.
867 std::vector<double> xp(np1, 0.);
868 std::vector<double> yp(np1, 0.);
869 std::vector<double> zp(np1, 0.);
870 double zm = 0.;
871 for (unsigned int k = 0; k < np1; ++k) {
872 xp[k] = rot[0][0] * xp1[k] + rot[0][1] * yp1[k] + rot[0][2] * zp1[k];
873 yp[k] = rot[1][0] * xp1[k] + rot[1][1] * yp1[k] + rot[1][2] * zp1[k];
874 zp[k] = rot[2][0] * xp1[k] + rot[2][1] * yp1[k] + rot[2][2] * zp1[k];
875 zm += zp[k];
876 }
877 zm /= np1;
878 // Store it.
879 std::vector<Panel> newPanels;
880 std::vector<int> vol1;
881 std::vector<int> vol2;
882 Panel panel1 = panelsIn[i];
883 panel1.xv = xp;
884 panel1.yv = yp;
885 panel1.zv = zp;
886 vol1.push_back(panel1.volume);
887 vol2.push_back(-1);
888 newPanels.push_back(std::move(panel1));
889 // Pick up all matching planes.
890 for (unsigned int j = i + 1; j < nIn; ++j) {
891 if (mark[j]) continue;
892 const double a2 = panelsIn[j].a;
893 const double b2 = panelsIn[j].b;
894 const double c2 = panelsIn[j].c;
895 const auto& xp2 = panelsIn[j].xv;
896 const auto& yp2 = panelsIn[j].yv;
897 const auto& zp2 = panelsIn[j].zv;
898 const unsigned int np2 = xp2.size();
899 // See whether this matches the first.
900 const double d2 = a2 * xp2[0] + b2 * yp2[0] + c2 * zp2[0];
901 // Inner product.
902 const double dot = a1 * a2 + b1 * b2 + c1 * c2;
903 // Offset between the two planes.
904 const double offset = d1 - d2 * dot;
905 if (fabs(fabs(dot) - 1.) > epsang || fabs(offset) > epsxyz) continue;
906 // Found a match.
907 mark[j] = true;
908 if (m_debug) std::cout << " Match with panel " << j << ".\n";
909 // Rotate this plane too.
910 xp.assign(np2, 0.);
911 yp.assign(np2, 0.);
912 zp.assign(np2, 0.);
913 zm = 0.;
914 for (unsigned int k = 0; k < np2; ++k) {
915 xp[k] = rot[0][0] * xp2[k] + rot[0][1] * yp2[k] + rot[0][2] * zp2[k];
916 yp[k] = rot[1][0] * xp2[k] + rot[1][1] * yp2[k] + rot[1][2] * zp2[k];
917 zp[k] = rot[2][0] * xp2[k] + rot[2][1] * yp2[k] + rot[2][2] * zp2[k];
918 zm += zp[k];
919 }
920 zm /= np2;
921 // Store it.
922 Panel panel2 = panelsIn[j];
923 panel2.xv = xp;
924 panel2.yv = yp;
925 panel2.zv = zp;
926 vol1.push_back(panel2.volume);
927 vol2.push_back(-1);
928 newPanels.push_back(std::move(panel2));
929 }
930 std::vector<bool> obsolete(newPanels.size(), false);
931 // Cut them as long as needed till no contacts remain.
932 unsigned int jmin = 0;
933 bool change = true;
934 while (change) {
935 change = false;
936 const unsigned int n = newPanels.size();
937 for (unsigned int j = 0; j < n; ++j) {
938 if (obsolete[j] || j < jmin) continue;
939 if (vol1[j] >= 0 && vol2[j] >= 0) continue;
940 const auto& panelj = newPanels[j];
941 for (unsigned int k = j + 1; k < n; ++k) {
942 if (obsolete[k]) continue;
943 if (vol1[k] >= 0 && vol2[k] >= 0) continue;
944 const auto& panelk = newPanels[k];
945 if (m_debug) std::cout << " Cutting " << j << ", " << k << ".\n";
946 // Separate contact and non-contact areas.
947 std::vector<Panel> panelsOut;
948 std::vector<int> itypo;
949 EliminateOverlaps(panelj, panelk, panelsOut, itypo);
950 const unsigned int nOut = panelsOut.size();
951 if (nOut == 2) {
952 // TODO: retrieve epsx, epsy from overlap finding?
953 const double epsx = epsxyz;
954 const double epsy = epsxyz;
955 // If there are just 2 panels, see whether there is a new one.
956 const bool equal1 = Equal(panelj, panelsOut[0], epsx, epsy);
957 const bool equal2 = Equal(panelj, panelsOut[1], epsx, epsy);
958 const bool equal3 = Equal(panelk, panelsOut[0], epsx, epsy);
959 const bool equal4 = Equal(panelk, panelsOut[1], epsx, epsy);
960 if ((equal1 || equal3) && (equal2 || equal4)) {
961 if (m_debug) {
962 std::cout << " Original and new panels are identical.\n";
963 }
964 } else {
965 change = true;
966 }
967 } else {
968 change = true;
969 }
970 if (m_debug) std::cout << " Change flag: " << change << ".\n";
971 // If there is no change, keep the two panels and proceed.
972 if (!change) continue;
973 // Flag the existing panels as inactive.
974 obsolete[j] = true;
975 obsolete[k] = true;
976
977 // Add the new panels.
978 for (unsigned int l = 0; l < nOut; ++l) {
979 if (itypo[l] == 1) {
980 vol1.push_back(std::max(vol1[j], vol2[j]));
981 vol2.push_back(-1);
982 } else if (itypo[l] == 2) {
983 vol1.push_back(std::max(vol1[k], vol2[k]));
984 vol2.push_back(-1);
985 } else {
986 vol1.push_back(std::max(vol1[j], vol2[j]));
987 vol2.push_back(std::max(vol1[k], vol2[k]));
988 }
989 newPanels.push_back(std::move(panelsOut[l]));
990 obsolete.push_back(false);
991 }
992 jmin = j + 1;
993 // Restart the loops.
994 break;
995 }
996 if (change) break;
997 }
998 }
999 // And rotate the panels back in place.
1000 const unsigned int nNew = newPanels.size();
1001 for (unsigned int j = 0; j < nNew; ++j) {
1002 if (obsolete[j]) continue;
1003 // Examine the boundary conditions.
1004 int interfaceType = 0;
1005 double potential = 0.;
1006 double lambda = 0.;
1007 double chargeDensity = 0.;
1008 if (m_debug) {
1009 std::cout << " Volume 1: " << vol1[j]
1010 << ". Volume 2: " << vol2[j] << ".\n";
1011 }
1012 if (vol1[j] < 0 && vol2[j] < 0) {
1013 // Shouldn't happen.
1014 continue;
1015 } else if (vol1[j] < 0 || vol2[j] < 0) {
1016 // Interface between a solid and vacuum/background.
1017 const auto vol = vol1[j] < 0 ? vol2[j] : vol1[j];
1018 interfaceType = InterfaceType(bc[vol]);
1019 if (bc[vol] == Solid::Dielectric ||
1020 bc[vol] == Solid::DielectricCharge) {
1021 if (fabs(eps[vol] - 1.) < 1.e-6) {
1022 // Same epsilon on both sides. Skip.
1023 interfaceType = 0;
1024 } else {
1025 lambda = (eps[vol] - 1.) / (eps[vol] + 1.);
1026 }
1027 } else if (bc[vol] == Solid::Voltage) {
1028 potential = volt[vol];
1029 }
1030 if (bc[vol] == Solid::Charge || bc[vol] == Solid::DielectricCharge) {
1031 chargeDensity = charge[vol];
1032 }
1033 } else {
1034 const auto bc1 = bc[vol1[j]];
1035 const auto bc2 = bc[vol2[j]];
1036 if (bc1 == Solid::Voltage || bc1 == Solid::Charge ||
1037 bc1 == Solid::Float) {
1038 interfaceType = InterfaceType(bc1);
1039 // First volume is a conductor. Other volume must be a dielectric.
1040 if (bc2 == Solid::Dielectric || bc2 == Solid::DielectricCharge) {
1041 if (bc1 == Solid::Voltage) {
1042 potential = volt[vol1[j]];
1043 } else if (bc1 == Solid::Charge) {
1044 chargeDensity = charge[vol1[j]];
1045 }
1046 } else {
1047 interfaceType = -1;
1048 }
1049 if (bc1 == Solid::Voltage && bc2 == Solid::Voltage) {
1050 const double v1 = volt[vol1[j]];
1051 const double v2 = volt[vol2[j]];
1052 if (fabs(v1 - v2) < 1.e-6 * (1. + fabs(v1) + fabs(v2))) {
1053 interfaceType = 0;
1054 }
1055 }
1056 } else if (bc1 == Solid::Dielectric ||
1057 bc1 == Solid::DielectricCharge) {
1058 interfaceType = InterfaceType(bc2);
1059 // First volume is a dielectric.
1060 if (bc2 == Solid::Voltage) {
1061 potential = volt[vol2[j]];
1062 } else if (bc2 == Solid::Charge) {
1063 chargeDensity = charge[vol2[j]];
1064 } else if (bc2 == Solid::Dielectric ||
1065 bc2 == Solid::DielectricCharge) {
1066 const double eps1 = eps[vol1[j]];
1067 const double eps2 = eps[vol2[j]];
1068 if (fabs(eps1 - eps2) < 1.e-6 * (1. + fabs(eps1) + fabs(eps2))) {
1069 // Same epsilon. Skip.
1070 interfaceType = 0;
1071 } else {
1072 lambda = (eps1 - eps2) / (eps1 + eps2);
1073 }
1074 }
1075 }
1076 }
1077 if (m_debug) {
1078 if (interfaceType < 0) {
1079 std::cout << " Conflicting boundary conditions. Skip.\n";
1080 } else if (interfaceType < 1) {
1081 std::cout << " Trivial interface. Skip.\n";
1082 } else if (interfaceType > 5) {
1083 std::cout << " Interface type " << interfaceType
1084 << " is not implemented. Skip.\n";
1085 }
1086 }
1087 if (interfaceType < 0) {
1088 ++nConflicting;
1089 } else if (interfaceType < 1) {
1090 ++nTrivial;
1091 } else if (interfaceType > 5) {
1092 ++nNotImplemented;
1093 }
1094 if (interfaceType < 1 || interfaceType > 5) continue;
1095
1096 std::vector<Panel> panelsOut;
1097 // Reduce to rectangles and right-angle triangles.
1098 if (m_debug) std::cout << " Creating primitives.\n";
1099 MakePrimitives(newPanels[j], panelsOut);
1100 // Loop over the rectangles and triangles.
1101 for (auto& panel : panelsOut) {
1102 const auto& up = panel.xv;
1103 const auto& vp = panel.yv;
1104 const auto& wp = panel.zv;
1105 const unsigned int np = up.size();
1106 // Rotate.
1107 xp.assign(np, 0.);
1108 yp.assign(np, 0.);
1109 zp.assign(np, 0.);
1110 for (unsigned int k = 0; k < np; ++k) {
1111 xp[k] = rot[0][0] * up[k] + rot[1][0] * vp[k] + rot[2][0] * wp[k];
1112 yp[k] = rot[0][1] * up[k] + rot[1][1] * vp[k] + rot[2][1] * wp[k];
1113 zp[k] = rot[0][2] * up[k] + rot[1][2] * vp[k] + rot[2][2] * wp[k];
1114 }
1115 Primitive primitive;
1116 primitive.a = panel.a;
1117 primitive.b = panel.b;
1118 primitive.c = panel.c;
1119 primitive.xv = xp;
1120 primitive.yv = yp;
1121 primitive.zv = zp;
1122 primitive.v = potential;
1123 primitive.q = chargeDensity;
1124 primitive.lambda = lambda;
1125 primitive.interface = interfaceType;
1126 // Set the requested discretization level (target element size).
1127 primitive.elementSize = -1.;
1128 if (solids.find(vol1[j]) != solids.end()) {
1129 const auto solid = solids[vol1[j]];
1130 if (solid) {
1131 primitive.elementSize = solid->GetDiscretisationLevel(panel);
1132 }
1133 }
1134 primitive.vol1 = vol1[j];
1135 primitive.vol2 = vol2[j];
1136 m_primitives.push_back(std::move(primitive));
1137 }
1138 }
1139 }
1140
1141 // Add the wires.
1142 for (unsigned int i = 0; i < nSolids; ++i) {
1143 const auto solid = m_geometry->GetSolid(i);
1144 if (!solid) continue;
1145 if (!solid->IsWire()) continue;
1146 double x0 = 0., y0 = 0., z0 = 0.;
1147 solid->GetCentre(x0, y0, z0);
1148 double dx = 0., dy = 0., dz = 1.;
1149 solid->GetDirection(dx, dy, dz);
1150 const double dnorm = sqrt(dx * dx + dy * dy + dz * dz);
1151 if (dnorm < Small) {
1152 std::cerr << m_className << "::Initialise:\n"
1153 << " Wire has zero norm direction vector; skipped.\n";
1154 continue;
1155 }
1156 dx /= dnorm;
1157 dy /= dnorm;
1158 dz /= dnorm;
1159 const double h = solid->GetHalfLengthZ();
1160 Primitive primitive;
1161 primitive.a = solid->GetRadius();
1162 primitive.b = 0.;
1163 primitive.c = 0.;
1164 primitive.xv = {x0 - h * dx, x0 + h * dx};
1165 primitive.yv = {y0 - h * dy, y0 + h * dy};
1166 primitive.zv = {z0 - h * dz, z0 + h * dz};
1167 primitive.v = solid->GetBoundaryPotential();
1168 primitive.q = solid->GetBoundaryChargeDensity();
1169 primitive.lambda = 0.;
1170 primitive.interface = InterfaceType(solid->GetBoundaryConditionType());
1171 // Set the requested discretization level (target element size).
1172 Panel panel;
1173 primitive.elementSize = solid->GetDiscretisationLevel(panel);
1174 primitive.vol1 = solid->GetId();
1175 primitive.vol2 = -1;
1176 m_primitives.push_back(std::move(primitive));
1177 }
1178 // Print a warning if we have discarded some panels during the process.
1179 if (nTrivial > 0 || nConflicting > 0 || nNotImplemented > 0) {
1180 std::cerr << m_className << "::Initialise:\n";
1181 if (nConflicting > 0) {
1182 std::cerr << " Skipped " << nConflicting
1183 << " panels with conflicting boundary conditions.\n";
1184 }
1185 if (nNotImplemented > 0) {
1186 std::cerr << " Skipped " << nNotImplemented
1187 << " panels with not yet available boundary conditions.\n";
1188 }
1189 if (nTrivial > 0) {
1190 std::cerr << " Skipped " << nTrivial
1191 << " panels with trivial boundary conditions.\n";
1192 }
1193 }
1194 if (m_debug) {
1195 std::cout << m_className << "::Initialise:\n"
1196 << " Created " << m_primitives.size() << " primitives.\n";
1197 }
1198
1199 // Discretize the primitives.
1200 for (const auto& primitive : m_primitives) {
1201 const auto nVertices = primitive.xv.size();
1202 if (nVertices < 2 || nVertices > 4) continue;
1203 std::vector<Element> elements;
1204 // Get the target element size.
1205 double targetSize = primitive.elementSize;
1206 if (targetSize < MinDist) targetSize = m_targetElementSize;
1207 if (nVertices == 2) {
1208 DiscretizeWire(primitive, targetSize, elements);
1209 } else if (nVertices == 3) {
1210 DiscretizeTriangle(primitive, targetSize, elements);
1211 } else if (nVertices == 4) {
1212 DiscretizeRectangle(primitive, targetSize, elements);
1213 }
1214 for (auto& element: elements) {
1215 element.interface = primitive.interface;
1216 element.lambda = primitive.lambda;
1217 element.bc = primitive.v;
1218 element.assigned = primitive.q;
1219 element.solution = 0.;
1220 }
1221 m_elements.insert(m_elements.end(),
1222 std::make_move_iterator(elements.begin()),
1223 std::make_move_iterator(elements.end()));
1224 }
1225 neBEM::NbThreads = m_nThreads;
1226 if (neBEM::neBEMInitialize() != 0) {
1227 std::cerr << m_className << "::Initialise:\n"
1228 << " Initialising neBEM failed.\n";
1229 return false;
1230 }
1231 gComponentNeBem3d = this;
1232 // Set the user options.
1233 neBEM::MinNbElementsOnLength = m_minNbElementsOnLength;
1234 neBEM::MaxNbElementsOnLength = m_maxNbElementsOnLength;
1235 neBEM::ElementLengthRqstd = m_targetElementSize * 0.01;
1236
1237 // New model / reuse existing model flag.
1238 // TODO!
1239 neBEM::NewModel = 1;
1240 neBEM::NewMesh = 1;
1241 neBEM::NewBC = 1;
1242 neBEM::NewPP = 1;
1243 // Pass debug level.
1244 neBEM::DebugLevel = m_debug ? 101 : 0;
1245 // Store inverted matrix or not.
1246 // TODO!
1247 neBEM::OptStoreInvMatrix = 0;
1248 neBEM::OptFormattedFile = 0;
1249 // Matrix inversion method (LU or SVD).
1250 if (m_inversion == Inversion::LU) {
1251 neBEM::OptInvMatProc = 0;
1252 } else {
1253 neBEM::OptInvMatProc = 1;
1254 }
1255 // Delete existing weighting fields (if any).
1256 if (neBEM::WtFieldChDen != NULL) {
1257 neBEM::neBEMDeleteAllWeightingFields();
1258 }
1259 // Transfer the geometry.
1260 if (neBEM::neBEMReadGeometry() != 0) {
1261 std::cerr << m_className << "::Initialise:\n"
1262 << " Transferring the geometry to neBEM failed.\n";
1263 return false;
1264 }
1265
1266 // Discretization.
1267 int** elementNbs = neBEM::imatrix(1, neBEM::NbPrimitives, 1, 2);
1268 for (int i = 1; i <= neBEM::NbPrimitives; ++i) {
1269 const int vol1 = neBEM::VolRef1[i];
1270 double size1 = -1.;
1271 if (solids.find(vol1) != solids.end()) {
1272 const auto solid = solids[vol1];
1273 if (solid) {
1274 Panel panel;
1275 panel.a = neBEM::XNorm[i];
1276 panel.b = neBEM::YNorm[i];
1277 panel.c = neBEM::ZNorm[i];
1278 const int nv = neBEM::NbVertices[i];
1279 panel.xv.resize(nv);
1280 panel.yv.resize(nv);
1281 panel.zv.resize(nv);
1282 for (int j = 0; j < nv; ++j) {
1283 panel.xv[j] = neBEM::XVertex[i][j];
1284 panel.yv[j] = neBEM::YVertex[i][j];
1285 panel.zv[j] = neBEM::ZVertex[i][j];
1286 }
1287 size1 = solid->GetDiscretisationLevel(panel);
1288 }
1289 }
1290 int vol2 = neBEM::VolRef2[i];
1291 double size2 = -1.;
1292 if (solids.find(vol2) != solids.end()) {
1293 const auto solid = solids[vol2];
1294 if (solid) {
1295 Panel panel;
1296 panel.a = -neBEM::XNorm[i];
1297 panel.b = -neBEM::YNorm[i];
1298 panel.c = -neBEM::ZNorm[i];
1299 const int nv = neBEM::NbVertices[i];
1300 panel.xv.resize(nv);
1301 panel.yv.resize(nv);
1302 panel.zv.resize(nv);
1303 for (int j = 0; j < nv; ++j) {
1304 panel.xv[j] = neBEM::XVertex[i][j];
1305 panel.yv[j] = neBEM::YVertex[i][j];
1306 panel.zv[j] = neBEM::ZVertex[i][j];
1307 }
1308 size2 = solid->GetDiscretisationLevel(panel);
1309 }
1310 }
1311 double size = m_targetElementSize;
1312 if (size1 > 0. && size2 > 0.) {
1313 size = std::min(size1, size2);
1314 } else if (size1 > 0.) {
1315 size = size1;
1316 } else if (size2 > 0.) {
1317 size = size2;
1318 }
1319 // Convert from cm to m.
1320 size *= 0.01;
1321
1322 // Work out the element dimensions.
1323 const double dx1 = neBEM::XVertex[i][0] - neBEM::XVertex[i][1];
1324 const double dy1 = neBEM::YVertex[i][0] - neBEM::YVertex[i][1];
1325 const double dz1 = neBEM::ZVertex[i][0] - neBEM::ZVertex[i][1];
1326 const double l1 = dx1 * dx1 + dy1 * dy1 + dz1 * dz1;
1327 int nb1 = (int)(sqrt(l1) / size);
1328 // Truncate to desired range.
1329 if (nb1 < neBEM::MinNbElementsOnLength) {
1330 nb1 = neBEM::MinNbElementsOnLength;
1331 } else if (nb1 > neBEM::MaxNbElementsOnLength) {
1332 nb1 = neBEM::MaxNbElementsOnLength;
1333 }
1334 int nb2 = 0;
1335 if (neBEM::NbVertices[i] > 2) {
1336 const double dx2 = neBEM::XVertex[i][2] - neBEM::XVertex[i][1];
1337 const double dy2 = neBEM::YVertex[i][2] - neBEM::YVertex[i][1];
1338 const double dz2 = neBEM::ZVertex[i][2] - neBEM::ZVertex[i][1];
1339 const double l2 = dx2 * dx2 + dy2 * dy2 + dz2 * dz2;
1340 nb2 = (int)(sqrt(l2) / size);
1341 if (nb2 < neBEM::MinNbElementsOnLength) {
1342 nb2 = neBEM::MinNbElementsOnLength;
1343 } else if (nb2 > neBEM::MaxNbElementsOnLength) {
1344 nb2 = neBEM::MaxNbElementsOnLength;
1345 }
1346 }
1347 elementNbs[i][1] = nb1;
1348 elementNbs[i][2] = nb2;
1349 }
1350
1351 if (neBEM::neBEMDiscretize(elementNbs) != 0) {
1352 std::cerr << m_className << "::Initialise: Discretization failed.\n";
1353 neBEM::free_imatrix(elementNbs, 1, neBEM::NbPrimitives, 1, 2);
1354 return false;
1355 }
1356 neBEM::free_imatrix(elementNbs, 1, neBEM::NbPrimitives, 1, 2);
1357 if (neBEM::neBEMBoundaryConditions() != 0) {
1358 std::cerr << m_className << "::Initialise:\n"
1359 << " Setting the boundary conditions failed.\n";
1360 return false;
1361 }
1362 if (neBEM::neBEMSolve() != 0) {
1363 std::cerr << m_className << "::Initialise: Solution failed.\n";
1364 return false;
1365 }
1366 // Now the weighting fields.
1367 std::set<std::string> labels;
1368 for (unsigned int i = 0; i < nSolids; ++i) {
1369 const auto solid = m_geometry->GetSolid(i);
1370 if (!solid) continue;
1371 const std::string label = solid->GetLabel();
1372 if (!label.empty()) labels.insert(label);
1373 }
1374 for (const auto& label : labels) {
1375 std::vector<int> primitives;
1376 for (unsigned int i = 0; i < nSolids; ++i) {
1377 const auto solid = m_geometry->GetSolid(i);
1378 if (!solid) continue;
1379 if (solid->GetLabel() != label) continue;
1380 const int id = solid->GetId();
1381 // Add the primitives associated to this solid to the list.
1382 for (int j = 1; j <= neBEM::NbPrimitives; ++j) {
1383 if (neBEM::VolRef1[j] == id || neBEM::VolRef2[j] == id) {
1384 primitives.push_back(j);
1385 }
1386 }
1387 }
1388 // Request the weighting field for this list of primitives.
1389 const int np = primitives.size();
1390 const int id = neBEM::neBEMPrepareWeightingField(np, primitives.data());
1391 if (id < 0) {
1392 std::cerr << m_className << "::Initialise:\n"
1393 << " Weighting field calculation for readout group \""
1394 << label << "\" failed.\n";
1395 continue;
1396 } else {
1397 std::cout << m_className << "::Initialise:\n"
1398 << " Prepared weighting field for readout group \""
1399 << label << "\".\n";
1400 m_wfields[label] = id;
1401 }
1402 }
1403 // TODO! Not sure if we should call this here.
1404 // neBEM::neBEMEnd();
1405 m_ready = true;
1406 return true;
1407}
Garfield::Component::m_debug
bool m_debug
Switch on/off debugging messages.
Definition: Component.hh:341
Garfield::Solid::GetLabel
std::string GetLabel() const
Return the label.
Definition: Solid.hh:68
Garfield::Solid::GetCentre
bool GetCentre(double &x, double &y, double &z) const
Retrieve the centre point of the solid.
Definition: Solid.hh:71
Garfield::Solid::DielectricCharge
@ DielectricCharge
Definition: Solid.hh:157
Garfield::Solid::SolidPanels
virtual bool SolidPanels(std::vector< Panel > &panels)=0
Retrieve the surface panels of the solid.
Garfield::gComponentNeBem3d
ComponentNeBem3d * gComponentNeBem3d
Definition: ComponentNeBem3d.cc:382
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314
InterfaceType
neBEMGLOBAL int * InterfaceType
Definition: neBEM.h:68

Referenced by ElectricField().

◆ Reset()

void Garfield::ComponentNeBem3d::Reset ( )
overrideprotectedvirtual

Reset the component.

Implements Garfield::Component.

Definition at line 3234 of file ComponentNeBem3d.cc.

3234 {
3235 m_primitives.clear();
3236 m_elements.clear();
3237 m_ynplan.fill(false);
3238 m_coplan.fill(0.);
3239 m_vtplan.fill(0.);
3240 m_ready = false;
3241}

◆ SetMinMaxNumberOfElements()

void Garfield::ComponentNeBem3d::SetMinMaxNumberOfElements ( const unsigned int  nmin,
const unsigned int  nmax 
)

Set the smallest and largest allowed number of elements along the lenght of a primitive.

Definition at line 628 of file ComponentNeBem3d.cc.

629 {
630
631 if (nmin == 0 || nmax == 0) {
632 std::cerr << m_className << "::SetMinMaxNumberOfElements:\n"
633 << " Values must be non-zero.\n";
634 return;
635 }
636 m_minNbElementsOnLength = std::min(nmin, nmax);
637 m_maxNbElementsOnLength = std::max(nmin, nmax);
638}

◆ SetMirrorPeriodicityX()

void Garfield::ComponentNeBem3d::SetMirrorPeriodicityX ( const double  s)

Set the periodic length [cm] in the x-direction.

Definition at line 684 of file ComponentNeBem3d.cc.

684 {
685 if (s < Small) {
686 std::cerr << m_className << "::SetMirrorPeriodicityX:\n"
687 << " Periodic length must be greater than zero.\n";
688 return;
689 }
690 m_periodicLength[0] = s;
691 m_periodic[0] = false;
692 m_mirrorPeriodic[0] = true;
693 UpdatePeriodicity();
694}
Garfield::ComponentNeBem3d::UpdatePeriodicity
void UpdatePeriodicity() override
Verify periodicities.
Definition: ComponentNeBem3d.cc:3243

◆ SetMirrorPeriodicityY()

void Garfield::ComponentNeBem3d::SetMirrorPeriodicityY ( const double  s)

Set the periodic length [cm] in the y-direction.

Definition at line 696 of file ComponentNeBem3d.cc.

696 {
697 if (s < Small) {
698 std::cerr << m_className << "::SetMirrorPeriodicityY:\n"
699 << " Periodic length must be greater than zero.\n";
700 return;
701 }
702 m_periodicLength[1] = s;
703 m_periodic[1] = false;
704 m_mirrorPeriodic[1] = true;
705 UpdatePeriodicity();
706}

◆ SetMirrorPeriodicityZ()

void Garfield::ComponentNeBem3d::SetMirrorPeriodicityZ ( const double  s)

Set the periodic length [cm] in the z-direction.

Definition at line 708 of file ComponentNeBem3d.cc.

708 {
709 if (s < Small) {
710 std::cerr << m_className << "::SetMirrorPeriodicityZ:\n"
711 << " Periodic length must be greater than zero.\n";
712 return;
713 }
714 m_periodicLength[2] = s;
715 m_periodic[2] = false;
716 m_mirrorPeriodic[2] = true;
717 UpdatePeriodicity();
718}

◆ SetNumberOfThreads()

void Garfield::ComponentNeBem3d::SetNumberOfThreads ( const unsigned int  n)
inline

Set the number of threads to be used by neBEM.

Definition at line 119 of file ComponentNeBem3d.hh.

119 {
120 m_nThreads = n > 0 ? n : 1;
121 }

◆ SetPeriodicCopies()

void Garfield::ComponentNeBem3d::SetPeriodicCopies ( const unsigned int  nx,
const unsigned int  ny,
const unsigned int  nz 
)

Set the parameters $n_x, n_y, n_z$ defining the number of periodic copies that neBEM will use when dealing with periodic configurations. neBEM will use $2 \times n + 1$ copies (default: $n = 5$).

Definition at line 640 of file ComponentNeBem3d.cc.

642 {
643 m_nCopiesX = nx;
644 m_nCopiesY = ny;
645 m_nCopiesZ = nz;
646}

◆ SetPeriodicityX()

void Garfield::ComponentNeBem3d::SetPeriodicityX ( const double  s)

Set the periodic length [cm] in the x-direction.

Definition at line 648 of file ComponentNeBem3d.cc.

648 {
649 if (s < Small) {
650 std::cerr << m_className << "::SetPeriodicityX:\n"
651 << " Periodic length must be greater than zero.\n";
652 return;
653 }
654 m_periodicLength[0] = s;
655 m_periodic[0] = true;
656 m_mirrorPeriodic[0] = false;
657 UpdatePeriodicity();
658}

◆ SetPeriodicityY()

void Garfield::ComponentNeBem3d::SetPeriodicityY ( const double  s)

Set the periodic length [cm] in the y-direction.

Definition at line 660 of file ComponentNeBem3d.cc.

660 {
661 if (s < Small) {
662 std::cerr << m_className << "::SetPeriodicityY:\n"
663 << " Periodic length must be greater than zero.\n";
664 return;
665 }
666 m_periodicLength[1] = s;
667 m_periodic[1] = true;
668 m_mirrorPeriodic[1] = false;
669 UpdatePeriodicity();
670}

◆ SetPeriodicityZ()

void Garfield::ComponentNeBem3d::SetPeriodicityZ ( const double  s)

Set the periodic length [cm] in the z-direction.

Definition at line 672 of file ComponentNeBem3d.cc.

672 {
673 if (s < Small) {
674 std::cerr << m_className << "::SetPeriodicityZ:\n"
675 << " Periodic length must be greater than zero.\n";
676 return;
677 }
678 m_periodicLength[2] = s;
679 m_periodic[2] = true;
680 m_mirrorPeriodic[2] = false;
681 UpdatePeriodicity();
682}

◆ SetTargetElementSize()

void Garfield::ComponentNeBem3d::SetTargetElementSize ( const double  length)

Set the default value of the target linear size of the elements produced by neBEM's discretisation process.

Definition at line 618 of file ComponentNeBem3d.cc.

618 {
619
620 if (length < MinDist) {
621 std::cerr << m_className << "::SetTargetElementSize: Value must be > "
622 << MinDist << ".\n";
623 return;
624 }
625 m_targetElementSize = length;
626}

◆ UpdatePeriodicity()

void Garfield::ComponentNeBem3d::UpdatePeriodicity ( )
overrideprotectedvirtual

Verify periodicities.

Implements Garfield::Component.

Definition at line 3243 of file ComponentNeBem3d.cc.

3243 {
3244
3245 for (unsigned int i = 0; i < 3; ++i) {
3246 // Cannot have regular and mirror periodicity at the same time.
3247 if (m_periodic[i] && m_mirrorPeriodic[i]) {
3248 std::cerr << m_className << "::UpdatePeriodicity:\n"
3249 << " Both simple and mirror periodicity requested. Reset.\n";
3250 m_periodic[i] = false;
3251 m_mirrorPeriodic[i] = false;
3252 continue;
3253 }
3254 if ((m_periodic[i] || m_mirrorPeriodic[i]) && m_periodicLength[i] < Small) {
3255 std::cerr << m_className << "::UpdatePeriodicity:\n"
3256 << " Periodic length is not set. Reset.\n";
3257 m_periodic[i] = m_mirrorPeriodic[i] = false;
3258 }
3259 }
3260
3261 if (m_axiallyPeriodic[0] || m_axiallyPeriodic[1] || m_axiallyPeriodic[2]) {
3262 std::cerr << m_className << "::UpdatePeriodicity:\n"
3263 << " Axial periodicity is not available.\n";
3264 m_axiallyPeriodic.fill(false);
3265 }
3266
3267 if (m_rotationSymmetric[0] || m_rotationSymmetric[1] ||
3268 m_rotationSymmetric[2]) {
3269 std::cerr << m_className << "::UpdatePeriodicity:\n"
3270 << " Rotation symmetry is not available.\n";
3271 m_rotationSymmetric.fill(false);
3272 }
3273
3274}
Garfield::Component::m_rotationSymmetric
std::array< bool, 3 > m_rotationSymmetric
Rotation symmetry around x-axis, y-axis, z-axis.
Definition: Component.hh:350
Garfield::Component::m_axiallyPeriodic
std::array< bool, 3 > m_axiallyPeriodic
Axial periodicity in x, y, z.
Definition: Component.hh:348

Referenced by SetMirrorPeriodicityX(), SetMirrorPeriodicityY(), SetMirrorPeriodicityZ(), SetPeriodicityX(), SetPeriodicityY(), and SetPeriodicityZ().

◆ UseLUInversion()

void Garfield::ComponentNeBem3d::UseLUInversion ( )
inline

Invert the influence matrix using lower-upper (LU) decomposition.

Definition at line 83 of file ComponentNeBem3d.hh.

83{ m_inversion = Inversion::LU; }

◆ UseSVDInversion()

void Garfield::ComponentNeBem3d::UseSVDInversion ( )
inline

Invert the influence matrix using singular value decomposition.

Definition at line 85 of file ComponentNeBem3d.hh.

85{ m_inversion = Inversion::SVD; }

◆ WeightingField()

void Garfield::ComponentNeBem3d::WeightingField ( const double  x,
const double  y,
const double  z,
double &  wx,
double &  wy,
double &  wz,
const std::string &  label 
)
overridevirtual

Calculate the weighting field at a given point and for a given electrode.

Parameters
x,y,zcoordinates [cm].
wx,wy,wzcomponents of the weighting field [1/cm].
labelname of the electrode

Reimplemented from Garfield::Component.

Definition at line 445 of file ComponentNeBem3d.cc.

448 {
449 wx = wy = wz = 0.;
450 if (m_wfields.count(label) == 0) return;
451 const int id = m_wfields[label];
452 neBEM::Point3D point;
453 point.X = 0.01 * x;
454 point.Y = 0.01 * y;
455 point.Z = 0.01 * z;
456 neBEM::Vector3D field;
457 if (neBEM::neBEMWeightingField(&point, &field, id) == DBL_MAX) {
458 std::cerr << m_className << "::WeightingField: Evaluation failed.\n";
459 return;
460 }
461 wx = 0.01 * field.X;
462 wy = 0.01 * field.Y;
463 wz = 0.01 * field.Z;
464}

◆ WeightingPotential()

double Garfield::ComponentNeBem3d::WeightingPotential ( const double  x,
const double  y,
const double  z,
const std::string &  label 
)
overridevirtual

Calculate the weighting potential at a given point.

Parameters
x,y,zcoordinates [cm].
labelname of the electrode.
Returns
weighting potential [dimensionless].

Reimplemented from Garfield::Component.

Definition at line 466 of file ComponentNeBem3d.cc.

468 {
469 if (m_wfields.count(label) == 0) return 0.;
470 const int id = m_wfields[label];
471 neBEM::Point3D point;
472 point.X = 0.01 * x;
473 point.Y = 0.01 * y;
474 point.Z = 0.01 * z;
475 neBEM::Vector3D field;
476 const double v = neBEM::neBEMWeightingField(&point, &field, id);
477 return v == DBL_MAX ? 0. : v;
478}
The documentation for this class was generated from the following files: