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 | 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 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 SetNewModel (const unsigned int NewModel)
 
void SetNewMesh (const unsigned int NewMesh)
 
void SetNewBC (const unsigned int NewBC)
 
void SetNewPP (const unsigned int NewPP)
 
void SetModelOptions (const unsigned int NewModel, const unsigned int NewMesh, const unsigned int NewBC, const unsigned int NewPP)
 
void SetStoreInflMatrix (const unsigned int OptStoreInflMatrix)
 
void SetReadInflMatrix (const unsigned int OptReadInflMatrix)
 
void SetStoreInvMatrix (const unsigned int OptStoreInvMatrix)
 
void SetReadInvMatrix (const unsigned int OptReadInvMatrix)
 
void SetStorePrimitives (const unsigned int OptStorePrimitives)
 
void SetReadPrimitives (const unsigned int OptReadPrimitives)
 
void SetStoreElements (const unsigned int OptStoreElements)
 
void SetReadElements (const unsigned int OptReadElements)
 
void SetFormattedFile (const unsigned int OptFormattedFile)
 
void SetUnformattedFile (const unsigned int OptUnformattedFile)
 
void SetStoreReadOptions (const unsigned int OptStoreInflMatrix, const unsigned int OptReadInflMatrix, const unsigned int OptStoreInvMatrix, const unsigned int OptReadInvMatrix, const unsigned int OptStorePrimitives, const unsigned int OptReadPrimitives, const unsigned int OptStoreElements, const unsigned int OptReadElements, const unsigned int OptFormattedFile, const unsigned int OptUnformattedFile)
 
void SetReuseModel (void)
 
void SetSystemChargeZero (const unsigned int OptSystemChargeZero)
 
void SetValidateSolution (const unsigned int OptValidateSolution)
 
void SetForceValidation (const unsigned int OptForceValidation)
 
void SetRepeatLHMatrix (const unsigned int OptRepeatLHMatrix)
 
void SetComputeOptions (const unsigned int OptSystemChargeZero, const unsigned int OptValidateSolution, const unsigned int OptForceValidation, const unsigned int OptRepeatLHMatrix)
 
void SetFastVolOptions (const unsigned int OptFastVol, const unsigned int OptCreateFastPF, const unsigned int OptReadFastPF)
 
void SetFastVolVersion (const unsigned int VersionFV)
 
void SetFastVolBlocks (const unsigned int NbBlocksFV)
 
void SetWtFldFastVolOptions (const unsigned int IdWtField, const unsigned int OptWtFldFastVol, const unsigned int OptCreateWtFldFastPF, const unsigned int OptReadWtFldFastPF)
 
void SetWtFldFastVolVersion (const unsigned int IdWtField, const unsigned int VersionWtFldFV)
 
void SetWtFldFastVolBlocks (const unsigned int IdWtField, const unsigned int NbBlocksWtFldFV)
 
void SetKnownChargeOptions (const unsigned int OptKnownCharge)
 
void SetChargingUpOptions (const unsigned int OptChargingUp)
 
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.
 
void SetPrimAfter (const int n)
 
void SetWtFldPrimAfter (const int n)
 
void SetOptRmPrim (const unsigned int n)
 Set option related to removal of primitives.
 
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
 
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 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 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)
 
virtual bool InTrapRadius (const double q0, const double x0, const double y0, const double z0, double &xw, double &yw, double &rw)
 
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 ()
 

Protected Member Functions

void Reset () override
 Reset the component.
 
void UpdatePeriodicity () override
 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 382 of file ComponentNeBem3d.cc.

382 : Component("NeBem3d") {
383 InitValues();
384}
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 492 of file ComponentNeBem3d.cc.

492 {
493 if (m_ynplan[0] && m_ynplan[1]) {
494 std::cerr << m_className << "::AddPlaneX:\n"
495 << " Cannot have more than two planes at constant x.\n";
496 return;
497 }
498
499 if (m_ynplan[0]) {
500 m_ynplan[1] = true;
501 if (x < m_coplan[0]) {
502 m_coplan[1] = m_coplan[0];
503 m_vtplan[1] = m_vtplan[0];
504 m_coplan[0] = x;
505 m_vtplan[0] = v;
506 } else {
507 m_coplan[1] = x;
508 m_vtplan[1] = v;
509 }
510 } else {
511 m_ynplan[0] = true;
512 m_coplan[0] = x;
513 m_vtplan[0] = v;
514 }
515 m_ready = false;
516}
Garfield::Component::m_className
std::string m_className
Class name.
Definition Component.hh:359
Garfield::Component::m_ready
bool m_ready
Ready for use?
Definition Component.hh:368

◆ AddPlaneY()

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

Add a plane at constant y.

Definition at line 518 of file ComponentNeBem3d.cc.

518 {
519 if (m_ynplan[2] && m_ynplan[3]) {
520 std::cerr << m_className << "::AddPlaneY:\n"
521 << " Cannot have more than two planes at constant y.\n";
522 return;
523 }
524
525 if (m_ynplan[2]) {
526 m_ynplan[3] = true;
527 if (y < m_coplan[2]) {
528 m_coplan[3] = m_coplan[2];
529 m_vtplan[3] = m_vtplan[2];
530 m_coplan[2] = y;
531 m_vtplan[2] = v;
532 } else {
533 m_coplan[3] = y;
534 m_vtplan[3] = v;
535 }
536 } else {
537 m_ynplan[2] = true;
538 m_coplan[2] = y;
539 m_vtplan[2] = v;
540 }
541 m_ready = false;
542}

◆ AddPlaneZ()

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

Add a plane at constant z.

Definition at line 544 of file ComponentNeBem3d.cc.

544 {
545 if (m_ynplan[4] && m_ynplan[5]) {
546 std::cerr << m_className << "::AddPlaneZ:\n"
547 << " Cannot have more than two planes at constant z.\n";
548 return;
549 }
550
551 if (m_ynplan[4]) {
552 m_ynplan[5] = true;
553 if (z < m_coplan[4]) {
554 m_coplan[5] = m_coplan[4];
555 m_vtplan[5] = m_vtplan[4];
556 m_coplan[4] = z;
557 m_vtplan[4] = v;
558 } else {
559 m_coplan[5] = z;
560 m_vtplan[5] = v;
561 }
562 } else {
563 m_ynplan[4] = true;
564 m_coplan[4] = z;
565 m_vtplan[4] = v;
566 }
567 m_ready = false;
568}

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

◆ ElectricField() [2/3]

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 405 of file ComponentNeBem3d.cc.

408 {
409 ex = ey = ez = v = 0.;
410 status = 0;
411 // Check if the requested point is inside a medium
412 m = GetMedium(x, y, z);
413 if (!m) {
414 status = -6;
415 } else if (!m->IsDriftable()) {
416 status = -5;
417 }
418
419 if (!m_ready) {
420 if (!Initialise()) {
421 std::cerr << m_className << "::ElectricField: Initialisation failed.\n";
422 status = -11;
423 return;
424 }
425 m_ready = true;
426 }
427
428 // Construct a point.
429 neBEM::Point3D point;
430 point.X = 0.01 * x;
431 point.Y = 0.01 * y;
432 point.Z = 0.01 * z;
433
434 // Compute the field.
435 neBEM::Vector3D field;
436 if (neBEM::neBEMPF(&point, &v, &field) != 0) {
437 status = -10;
438 return;
439 }
440 ex = 0.01 * field.X;
441 ey = 0.01 * field.Y;
442 ez = 0.01 * field.Z;
443}
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() [3/3]

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 445 of file ComponentNeBem3d.cc.

447 {
448 double v = 0.;
449 ElectricField(x, y, z, ex, ey, ez, v, m, status);
450}

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 3438 of file ComponentNeBem3d.cc.

3441 {
3442 if (i >= m_elements.size()) {
3443 std::cerr << m_className << "::GetElement: Index out of range.\n";
3444 return false;
3445 }
3446 const auto& element = m_elements[i];
3447 xv = element.xv;
3448 yv = element.yv;
3449 zv = element.zv;
3450 interface = element.interface;
3451 bc = element.bc;
3452 lambda = element.lambda;
3453 return true;
3454}

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

Referenced by ElectricField().

◆ GetNumberOfElements()

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

Definition at line 52 of file ComponentNeBem3d.hh.

52{ return m_elements.size(); }

◆ GetNumberOfPlanesX()

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

Get the number of equipotential planes at constant x.

Definition at line 570 of file ComponentNeBem3d.cc.

570 {
571 if (m_ynplan[0] && m_ynplan[1]) {
572 return 2;
573 } else if (m_ynplan[0] || m_ynplan[1]) {
574 return 1;
575 }
576 return 0;
577}

◆ GetNumberOfPlanesY()

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

Get the number of equipotential planes at constant y.

Definition at line 579 of file ComponentNeBem3d.cc.

579 {
580 if (m_ynplan[2] && m_ynplan[3]) {
581 return 2;
582 } else if (m_ynplan[2] || m_ynplan[3]) {
583 return 1;
584 }
585 return 0;
586}

◆ GetNumberOfPlanesZ()

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

Get the number of equipotential planes at constant z.

Definition at line 588 of file ComponentNeBem3d.cc.

588 {
589 if (m_ynplan[4] && m_ynplan[5]) {
590 return 2;
591 } else if (m_ynplan[4] || m_ynplan[5]) {
592 return 1;
593 }
594 return 0;
595}

◆ GetNumberOfPrimitives()

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

Definition at line 40 of file ComponentNeBem3d.hh.

40{ return m_primitives.size(); }

◆ 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 157 of file ComponentNeBem3d.hh.

158 {
159 nx = m_nCopiesX;
160 ny = m_nCopiesY;
161 nz = m_nCopiesZ;
162 }

◆ GetPeriodicityX()

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

Get the periodic length in the x-direction.

Definition at line 906 of file ComponentNeBem3d.cc.

906 {
907 if (!m_periodic[0] && !m_mirrorPeriodic[0]) {
908 s = 0.;
909 return false;
910 }
911 s = m_periodicLength[0];
912 return true;
913}
Garfield::Component::m_mirrorPeriodic
std::array< bool, 3 > m_mirrorPeriodic
Mirror periodicity in x, y, z.
Definition Component.hh:376
Garfield::Component::m_periodic
std::array< bool, 3 > m_periodic
Simple periodicity in x, y, z.
Definition Component.hh:374

◆ GetPeriodicityY()

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

Get the periodic length in the y-direction.

Definition at line 915 of file ComponentNeBem3d.cc.

915 {
916 if (!m_periodic[1] && !m_mirrorPeriodic[1]) {
917 s = 0.;
918 return false;
919 }
920 s = m_periodicLength[1];
921 return true;
922}

◆ GetPeriodicityZ()

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

Get the periodic length in the z-direction.

Definition at line 924 of file ComponentNeBem3d.cc.

924 {
925 if (!m_periodic[2] && !m_mirrorPeriodic[2]) {
926 s = 0.;
927 return false;
928 }
929 s = m_periodicLength[2];
930 return true;
931}

◆ 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 597 of file ComponentNeBem3d.cc.

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

◆ 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 608 of file ComponentNeBem3d.cc.

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

◆ 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 619 of file ComponentNeBem3d.cc.

620 {
621 if (i >= 2 || (i == 1 && !m_ynplan[5])) {
622 std::cerr << m_className << "::GetPlaneZ: Index out of range.\n";
623 return false;
624 }
625 z = m_coplan[i + 4];
626 v = m_vtplan[i + 4];
627 return true;
628}

◆ 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 3344 of file ComponentNeBem3d.cc.

3349 {
3350 if (i >= m_primitives.size()) {
3351 std::cerr << m_className << "::GetPrimitive: Index out of range.\n";
3352 return false;
3353 }
3354 const auto& primitive = m_primitives[i];
3355 a = primitive.a;
3356 b = primitive.b;
3357 c = primitive.c;
3358 xv = primitive.xv;
3359 yv = primitive.yv;
3360 zv = primitive.zv;
3361 interface = primitive.interface;
3362 v = primitive.v;
3363 q = primitive.q;
3364 lambda = primitive.lambda;
3365 return true;
3366}

◆ 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 3368 of file ComponentNeBem3d.cc.

3372 {
3373 if (i >= m_primitives.size()) {
3374 std::cerr << m_className << "::GetPrimitive: Index out of range.\n";
3375 return false;
3376 }
3377 const auto& primitive = m_primitives[i];
3378 a = primitive.a;
3379 b = primitive.b;
3380 c = primitive.c;
3381 xv = primitive.xv;
3382 yv = primitive.yv;
3383 zv = primitive.zv;
3384 vol1 = primitive.vol1;
3385 vol2 = primitive.vol2;
3386 return true;
3387}

◆ GetVoltageRange()

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

Calculate the voltage range [V].

Implements Garfield::Component.

Definition at line 452 of file ComponentNeBem3d.cc.

452 {
453 // Voltage and other bc have to come from the solids
454 vmin = vmax = 0;
455 return true;
456}

◆ GetVolume() [1/2]

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

Definition at line 3425 of file ComponentNeBem3d.cc.

3426 {
3427 if (!m_geometry) return -1;
3428 const unsigned int nSolids = m_geometry->GetNumberOfSolids();
3429 for (unsigned int i = 0; i < nSolids; ++i) {
3430 Medium* medium = nullptr;
3431 const auto solid = m_geometry->GetSolid(i, medium);
3432 if (!solid) continue;
3433 if (solid->IsInside(x, y, z)) return solid->GetId();
3434 }
3435 return -1;
3436}

◆ 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 3389 of file ComponentNeBem3d.cc.

3391 {
3392 if (!m_geometry) return false;
3393 const unsigned int nSolids = m_geometry->GetNumberOfSolids();
3394 for (unsigned int i = 0; i < nSolids; ++i) {
3395 Medium* medium = nullptr;
3396 const auto solid = m_geometry->GetSolid(i, medium);
3397 if (!solid) continue;
3398 if (solid->GetId() != vol) continue;
3399 if (solid->IsTube() || solid->IsWire()) {
3400 shape = 1;
3401 } else if (solid->IsHole()) {
3402 shape = 2;
3403 } else if (solid->IsBox()) {
3404 shape = 3;
3405 } else if (solid->IsSphere()) {
3406 shape = 4;
3407 } else if (solid->IsRidge()) {
3408 shape = 5;
3409 } else if (solid->IsExtrusion()) {
3410 shape = 6;
3411 } else {
3412 std::cerr << m_className << "::GetVolume: Unknown solid shape.\n";
3413 return false;
3414 }
3415 material = medium ? medium->GetId() : 11;
3416 epsilon = medium ? medium->GetDielectricConstant() : 1.;
3417 potential = solid->GetBoundaryPotential();
3418 charge = solid->GetBoundaryChargeDensity();
3419 bc = solid->GetBoundaryConditionType();
3420 return true;
3421 }
3422 return false;
3423}

◆ Initialise()

bool Garfield::ComponentNeBem3d::Initialise ( )

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

Definition at line 933 of file ComponentNeBem3d.cc.

933 {
934 // Reset the lists.
935 m_primitives.clear();
936 m_elements.clear();
937
938 if (!m_geometry) {
939 std::cerr << m_className << "::Initialise: Geometry not set.\n";
940 return false;
941 }
942 // Be sure we won't have intersections with the bounding box.
943 // TODO! Loop over the solids and call PLACYE, PLACHE, PLABXE, PLASPE, ...
944 // Loop over the solids.
945 std::map<int, Solid*> solids;
946 std::map<int, Solid::BoundaryCondition> bc;
947 std::map<int, double> volt;
948 std::map<int, double> eps;
949 std::map<int, double> charge;
950 const unsigned int nSolids = m_geometry->GetNumberOfSolids();
951 std::vector<Panel> panelsIn;
952 for (unsigned int i = 0; i < nSolids; ++i) {
953 Medium* medium = nullptr;
954 const auto solid = m_geometry->GetSolid(i, medium);
955 if (!solid) continue;
956 // Get the panels.
957 solid->SolidPanels(panelsIn);
958 // Get the boundary condition.
959 const auto id = solid->GetId();
960 solids[id] = solid;
961 bc[id] = solid->GetBoundaryConditionType();
962 volt[id] = solid->GetBoundaryPotential();
963 if (bc[id] == Solid::Unknown) {
964 std::cout << m_className << "::Initialise:\n"
965 << " Boundary conditions for solid " << id << " not set.\n";
966 if (medium && medium->IsConductor()) {
967 std::cout << " Assuming the panels to be grounded.\n";
968 bc[id] = Solid::Voltage;
969 volt[id] = 0.;
970 } else {
971 std::cout << " Assuming dielectric-dielectric interfaces.\n";
972 bc[id] = Solid::Dielectric;
973 }
974 }
975 charge[id] = solid->GetBoundaryChargeDensity();
976 if (!medium) {
977 eps[id] = 1.;
978 } else {
979 eps[id] = medium->GetDielectricConstant();
980 }
981 }
982 // Apply cuts.
983 // CALL CELSCT('APPLY')
984 // Reduce to basic periodic copy.
985 ShiftPanels(panelsIn);
986
987 // Find contact panels and split into primitives.
988
989 // *---------------------------------------------------------------------
990 // * PLABEM - Prepares panels for BEM applications: removes the contacts
991 // * and cuts polygons to rectangles and right-angle triangles.
992 // *---------------------------------------------------------------------
993
994 // Establish tolerances.
995 const double epsang = 1.e-6; // BEMEPA
996 const double epsxyz = 1.e-6; // BEMEPD
997 // CALL EPSSET('SET',EPSXYZ,EPSXYZ,EPSXYZ)
998
999 const unsigned int nIn = panelsIn.size();
1000 if (m_debug) {
1001 std::cout << m_className << "::Initialise: Retrieved " << nIn
1002 << " panels from the solids.\n";
1003 }
1004 // Keep track of which panels have been processed.
1005 std::vector<bool> mark(nIn, false);
1006 // Count the number of interface panels that have been discarded.
1007 unsigned int nTrivial = 0;
1008 unsigned int nConflicting = 0;
1009 unsigned int nNotImplemented = 0;
1010 // Pick up panels which coincide potentially.
1011 for (unsigned int i = 0; i < nIn; ++i) {
1012 // Skip panels already done.
1013 if (mark[i]) continue;
1014 // Fetch panel parameters.
1015 const double a1 = panelsIn[i].a;
1016 const double b1 = panelsIn[i].b;
1017 const double c1 = panelsIn[i].c;
1018 const auto& xp1 = panelsIn[i].xv;
1019 const auto& yp1 = panelsIn[i].yv;
1020 const auto& zp1 = panelsIn[i].zv;
1021 const unsigned int np1 = xp1.size();
1022 // Establish its norm and offset.
1023 const double d1 = a1 * xp1[0] + b1 * yp1[0] + c1 * zp1[0];
1024 if (m_debug) {
1025 std::cout << " Panel " << i << "\n Norm vector: " << a1 << ", " << b1
1026 << ", " << c1 << ", " << d1 << ".\n";
1027 }
1028 // Rotation matrix.
1029 std::array<std::array<double, 3>, 3> rot;
1030 if (fabs(c1) <= fabs(a1) && fabs(c1) <= fabs(b1)) {
1031 // Rotation: removing C
1032 rot[0][0] = b1 / sqrt(a1 * a1 + b1 * b1);
1033 rot[0][1] = -a1 / sqrt(a1 * a1 + b1 * b1);
1034 rot[0][2] = 0.;
1035 } else if (fabs(b1) <= fabs(a1) && fabs(b1) <= fabs(c1)) {
1036 // Rotation: removing B
1037 rot[0][0] = c1 / sqrt(a1 * a1 + c1 * c1);
1038 rot[0][1] = 0.;
1039 rot[0][2] = -a1 / sqrt(a1 * a1 + c1 * c1);
1040 } else {
1041 // Rotation: removing A
1042 rot[0][0] = 0.;
1043 rot[0][1] = c1 / sqrt(b1 * b1 + c1 * c1);
1044 rot[0][2] = -b1 / sqrt(b1 * b1 + c1 * c1);
1045 }
1046 rot[2][0] = a1;
1047 rot[2][1] = b1;
1048 rot[2][2] = c1;
1049 rot[1][0] = rot[2][1] * rot[0][2] - rot[2][2] * rot[0][1];
1050 rot[1][1] = rot[2][2] * rot[0][0] - rot[2][0] * rot[0][2];
1051 rot[1][2] = rot[2][0] * rot[0][1] - rot[2][1] * rot[0][0];
1052 // Rotate to the x, y plane.
1053 std::vector<double> xp(np1, 0.);
1054 std::vector<double> yp(np1, 0.);
1055 std::vector<double> zp(np1, 0.);
1056 for (unsigned int k = 0; k < np1; ++k) {
1057 xp[k] = rot[0][0] * xp1[k] + rot[0][1] * yp1[k] + rot[0][2] * zp1[k];
1058 yp[k] = rot[1][0] * xp1[k] + rot[1][1] * yp1[k] + rot[1][2] * zp1[k];
1059 zp[k] = rot[2][0] * xp1[k] + rot[2][1] * yp1[k] + rot[2][2] * zp1[k];
1060 }
1061 // Store it.
1062 std::vector<Panel> newPanels;
1063 std::vector<int> vol1;
1064 std::vector<int> vol2;
1065 Panel panel1 = panelsIn[i];
1066 panel1.xv = xp;
1067 panel1.yv = yp;
1068 panel1.zv = zp;
1069 vol1.push_back(panel1.volume);
1070 vol2.push_back(-1);
1071 newPanels.push_back(std::move(panel1));
1072 // Pick up all matching planes.
1073 for (unsigned int j = i + 1; j < nIn; ++j) {
1074 if (mark[j]) continue;
1075 const double a2 = panelsIn[j].a;
1076 const double b2 = panelsIn[j].b;
1077 const double c2 = panelsIn[j].c;
1078 const auto& xp2 = panelsIn[j].xv;
1079 const auto& yp2 = panelsIn[j].yv;
1080 const auto& zp2 = panelsIn[j].zv;
1081 const unsigned int np2 = xp2.size();
1082 // See whether this matches the first.
1083 const double d2 = a2 * xp2[0] + b2 * yp2[0] + c2 * zp2[0];
1084 // Inner product.
1085 const double dot = a1 * a2 + b1 * b2 + c1 * c2;
1086 // Offset between the two planes.
1087 const double offset = d1 - d2 * dot;
1088 if (fabs(fabs(dot) - 1.) > epsang || fabs(offset) > epsxyz) continue;
1089 // Found a match.
1090 mark[j] = true;
1091 if (m_debug) std::cout << " Match with panel " << j << ".\n";
1092 // Rotate this plane too.
1093 xp.assign(np2, 0.);
1094 yp.assign(np2, 0.);
1095 zp.assign(np2, 0.);
1096 for (unsigned int k = 0; k < np2; ++k) {
1097 xp[k] = rot[0][0] * xp2[k] + rot[0][1] * yp2[k] + rot[0][2] * zp2[k];
1098 yp[k] = rot[1][0] * xp2[k] + rot[1][1] * yp2[k] + rot[1][2] * zp2[k];
1099 zp[k] = rot[2][0] * xp2[k] + rot[2][1] * yp2[k] + rot[2][2] * zp2[k];
1100 }
1101 // Store it.
1102 Panel panel2 = panelsIn[j];
1103 panel2.xv = xp;
1104 panel2.yv = yp;
1105 panel2.zv = zp;
1106 vol1.push_back(panel2.volume);
1107 vol2.push_back(-1);
1108 newPanels.push_back(std::move(panel2));
1109 }
1110 std::vector<bool> obsolete(newPanels.size(), false);
1111 // Cut them as long as needed till no contacts remain.
1112 unsigned int jmin = 0;
1113 bool change = true;
1114 while (change) {
1115 change = false;
1116 const unsigned int n = newPanels.size();
1117 for (unsigned int j = 0; j < n; ++j) {
1118 if (obsolete[j] || j < jmin) continue;
1119 if (vol1[j] >= 0 && vol2[j] >= 0) continue;
1120 const auto& panelj = newPanels[j];
1121 for (unsigned int k = j + 1; k < n; ++k) {
1122 if (obsolete[k]) continue;
1123 if (vol1[k] >= 0 && vol2[k] >= 0) continue;
1124 const auto& panelk = newPanels[k];
1125 if (m_debug) std::cout << " Cutting " << j << ", " << k << ".\n";
1126 // Separate contact and non-contact areas.
1127 std::vector<Panel> panelsOut;
1128 std::vector<int> itypo;
1129 EliminateOverlaps(panelj, panelk, panelsOut, itypo);
1130 const unsigned int nOut = panelsOut.size();
1131 if (nOut == 2) {
1132 // TODO: retrieve epsx, epsy from overlap finding?
1133 const double epsx = epsxyz;
1134 const double epsy = epsxyz;
1135 // If there are just 2 panels, see whether there is a new one.
1136 const bool equal1 = Equal(panelj, panelsOut[0], epsx, epsy);
1137 const bool equal2 = Equal(panelj, panelsOut[1], epsx, epsy);
1138 const bool equal3 = Equal(panelk, panelsOut[0], epsx, epsy);
1139 const bool equal4 = Equal(panelk, panelsOut[1], epsx, epsy);
1140 if ((equal1 || equal3) && (equal2 || equal4)) {
1141 if (m_debug) {
1142 std::cout << " Original and new panels are identical.\n";
1143 }
1144 } else {
1145 change = true;
1146 }
1147 } else {
1148 change = true;
1149 }
1150 if (m_debug) std::cout << " Change flag: " << change << ".\n";
1151 // If there is no change, keep the two panels and proceed.
1152 if (!change) continue;
1153 // Flag the existing panels as inactive.
1154 obsolete[j] = true;
1155 obsolete[k] = true;
1156
1157 // Add the new panels.
1158 for (unsigned int l = 0; l < nOut; ++l) {
1159 if (itypo[l] == 1) {
1160 vol1.push_back(std::max(vol1[j], vol2[j]));
1161 vol2.push_back(-1);
1162 } else if (itypo[l] == 2) {
1163 vol1.push_back(std::max(vol1[k], vol2[k]));
1164 vol2.push_back(-1);
1165 } else {
1166 vol1.push_back(std::max(vol1[j], vol2[j]));
1167 vol2.push_back(std::max(vol1[k], vol2[k]));
1168 }
1169 newPanels.push_back(std::move(panelsOut[l]));
1170 obsolete.push_back(false);
1171 }
1172 jmin = j + 1;
1173 // Restart the loops.
1174 break;
1175 }
1176 if (change) break;
1177 }
1178 }
1179 // And rotate the panels back in place.
1180 const unsigned int nNew = newPanels.size();
1181 for (unsigned int j = 0; j < nNew; ++j) {
1182 if (obsolete[j]) continue;
1183 // Examine the boundary conditions.
1184 int interfaceType = 0;
1185 double potential = 0.;
1186 double lambda = 0.;
1187 double chargeDensity = 0.;
1188 if (m_debug) {
1189 std::cout << " Volume 1: " << vol1[j] << ". Volume 2: " << vol2[j]
1190 << ".\n";
1191 }
1192 if (vol1[j] < 0 && vol2[j] < 0) {
1193 // Shouldn't happen.
1194 continue;
1195 } else if (vol1[j] < 0 || vol2[j] < 0) {
1196 // Interface between a solid and vacuum/background.
1197 const auto vol = vol1[j] < 0 ? vol2[j] : vol1[j];
1198 interfaceType = InterfaceType(bc[vol]);
1199 if (bc[vol] == Solid::Dielectric ||
1200 bc[vol] == Solid::DielectricCharge) {
1201 if (fabs(eps[vol] - 1.) < 1.e-6) {
1202 // Same epsilon on both sides. Skip.
1203 interfaceType = 0;
1204 } else {
1205 lambda = (eps[vol] - 1.) / (eps[vol] + 1.);
1206 }
1207 } else if (bc[vol] == Solid::Voltage) {
1208 potential = volt[vol];
1209 }
1210 if (bc[vol] == Solid::Charge || bc[vol] == Solid::DielectricCharge) {
1211 chargeDensity = charge[vol];
1212 }
1213 } else {
1214 const auto bc1 = bc[vol1[j]];
1215 const auto bc2 = bc[vol2[j]];
1216 if (bc1 == Solid::Voltage || bc1 == Solid::Charge ||
1217 bc1 == Solid::Float) {
1218 interfaceType = InterfaceType(bc1);
1219 // First volume is a conductor. Other volume must be a dielectric.
1220 if (bc2 == Solid::Dielectric || bc2 == Solid::DielectricCharge) {
1221 if (bc1 == Solid::Voltage) {
1222 potential = volt[vol1[j]];
1223 } else if (bc1 == Solid::Charge) {
1224 chargeDensity = charge[vol1[j]];
1225 }
1226 } else {
1227 interfaceType = -1;
1228 }
1229 if (bc1 == Solid::Voltage && bc2 == Solid::Voltage) {
1230 const double v1 = volt[vol1[j]];
1231 const double v2 = volt[vol2[j]];
1232 if (fabs(v1 - v2) < 1.e-6 * (1. + fabs(v1) + fabs(v2))) {
1233 interfaceType = 0;
1234 }
1235 }
1236 } else if (bc1 == Solid::Dielectric || bc1 == Solid::DielectricCharge) {
1237 interfaceType = InterfaceType(bc2);
1238 // First volume is a dielectric.
1239 if (bc2 == Solid::Voltage) {
1240 potential = volt[vol2[j]];
1241 } else if (bc2 == Solid::Charge) {
1242 chargeDensity = charge[vol2[j]];
1243 } else if (bc2 == Solid::Dielectric ||
1244 bc2 == Solid::DielectricCharge) {
1245 const double eps1 = eps[vol1[j]];
1246 const double eps2 = eps[vol2[j]];
1247 if (fabs(eps1 - eps2) < 1.e-6 * (1. + fabs(eps1) + fabs(eps2))) {
1248 // Same epsilon. Skip.
1249 interfaceType = 0;
1250 } else {
1251 lambda = (eps1 - eps2) / (eps1 + eps2);
1252 }
1253 }
1254 }
1255 }
1256 if (m_debug) {
1257 if (interfaceType < 0) {
1258 std::cout << " Conflicting boundary conditions. Skip.\n";
1259 } else if (interfaceType < 1) {
1260 std::cout << " Trivial interface. Skip.\n";
1261 } else if (interfaceType > 5) {
1262 std::cout << " Interface type " << interfaceType
1263 << " is not implemented. Skip.\n";
1264 }
1265 }
1266 if (interfaceType < 0) {
1267 ++nConflicting;
1268 } else if (interfaceType < 1) {
1269 ++nTrivial;
1270 } else if (interfaceType > 5) {
1271 ++nNotImplemented;
1272 }
1273 if (interfaceType < 1 || interfaceType > 5) continue;
1274
1275 std::vector<Panel> panelsOut;
1276 // Reduce to rectangles and right-angle triangles.
1277 if (m_debug) std::cout << " Creating primitives.\n";
1278 MakePrimitives(newPanels[j], panelsOut);
1279 // Loop over the rectangles and triangles.
1280 for (auto& panel : panelsOut) {
1281 const auto& up = panel.xv;
1282 const auto& vp = panel.yv;
1283 const auto& wp = panel.zv;
1284 const unsigned int np = up.size();
1285 // Rotate.
1286 xp.assign(np, 0.);
1287 yp.assign(np, 0.);
1288 zp.assign(np, 0.);
1289 for (unsigned int k = 0; k < np; ++k) {
1290 xp[k] = rot[0][0] * up[k] + rot[1][0] * vp[k] + rot[2][0] * wp[k];
1291 yp[k] = rot[0][1] * up[k] + rot[1][1] * vp[k] + rot[2][1] * wp[k];
1292 zp[k] = rot[0][2] * up[k] + rot[1][2] * vp[k] + rot[2][2] * wp[k];
1293 }
1294 Primitive primitive;
1295 primitive.a = panel.a;
1296 primitive.b = panel.b;
1297 primitive.c = panel.c;
1298 primitive.xv = xp;
1299 primitive.yv = yp;
1300 primitive.zv = zp;
1301 primitive.v = potential;
1302 primitive.q = chargeDensity;
1303 primitive.lambda = lambda;
1304 primitive.interface = interfaceType;
1305 // Set the requested discretization level (target element size).
1306 primitive.elementSize = -1.;
1307 if (solids.find(vol1[j]) != solids.end()) {
1308 const auto solid = solids[vol1[j]];
1309 if (solid) {
1310 primitive.elementSize = solid->GetDiscretisationLevel(panel);
1311 }
1312 }
1313 primitive.vol1 = vol1[j];
1314 primitive.vol2 = vol2[j];
1315 m_primitives.push_back(std::move(primitive));
1316 }
1317 }
1318 }
1319
1320 // Add the wires.
1321 for (unsigned int i = 0; i < nSolids; ++i) {
1322 const auto solid = m_geometry->GetSolid(i);
1323 if (!solid) continue;
1324 if (!solid->IsWire()) continue;
1325 double x0 = 0., y0 = 0., z0 = 0.;
1326 solid->GetCentre(x0, y0, z0);
1327 double dx = 0., dy = 0., dz = 1.;
1328 solid->GetDirection(dx, dy, dz);
1329 const double dnorm = sqrt(dx * dx + dy * dy + dz * dz);
1330 if (dnorm < Small) {
1331 std::cerr << m_className << "::Initialise:\n"
1332 << " Wire has zero norm direction vector; skipped.\n";
1333 continue;
1334 }
1335 dx /= dnorm;
1336 dy /= dnorm;
1337 dz /= dnorm;
1338 const double h = solid->GetHalfLengthZ();
1339 Primitive primitive;
1340 primitive.a = solid->GetRadius();
1341 primitive.b = 0.;
1342 primitive.c = 0.;
1343 primitive.xv = {x0 - h * dx, x0 + h * dx};
1344 primitive.yv = {y0 - h * dy, y0 + h * dy};
1345 primitive.zv = {z0 - h * dz, z0 + h * dz};
1346 primitive.v = solid->GetBoundaryPotential();
1347 primitive.q = solid->GetBoundaryChargeDensity();
1348 primitive.lambda = 0.;
1349 primitive.interface = InterfaceType(solid->GetBoundaryConditionType());
1350 // Set the requested discretization level (target element size).
1351 Panel panel;
1352 primitive.elementSize = solid->GetDiscretisationLevel(panel);
1353 primitive.vol1 = solid->GetId();
1354 primitive.vol2 = -1;
1355 m_primitives.push_back(std::move(primitive));
1356 }
1357 // Print a warning if we have discarded some panels during the process.
1358 if (nTrivial > 0 || nConflicting > 0 || nNotImplemented > 0) {
1359 std::cerr << m_className << "::Initialise:\n";
1360 if (nConflicting > 0) {
1361 std::cerr << " Skipped " << nConflicting
1362 << " panels with conflicting boundary conditions.\n";
1363 }
1364 if (nNotImplemented > 0) {
1365 std::cerr << " Skipped " << nNotImplemented
1366 << " panels with not yet available boundary conditions.\n";
1367 }
1368 if (nTrivial > 0) {
1369 std::cerr << " Skipped " << nTrivial
1370 << " panels with trivial boundary conditions.\n";
1371 }
1372 }
1373 if (m_debug) {
1374 std::cout << m_className << "::Initialise:\n"
1375 << " Created " << m_primitives.size() << " primitives.\n";
1376 }
1377
1378 // Discretize the primitives.
1379 for (const auto& primitive : m_primitives) {
1380 const auto nVertices = primitive.xv.size();
1381 if (nVertices < 2 || nVertices > 4) continue;
1382 std::vector<Element> elements;
1383 // Get the target element size.
1384 double targetSize = primitive.elementSize;
1385 if (targetSize < MinDist) targetSize = m_targetElementSize;
1386 if (nVertices == 2) {
1387 DiscretizeWire(primitive, targetSize, elements);
1388 } else if (nVertices == 3) {
1389 DiscretizeTriangle(primitive, targetSize, elements);
1390 } else if (nVertices == 4) {
1391 DiscretizeRectangle(primitive, targetSize, elements);
1392 }
1393 for (auto& element : elements) {
1394 element.interface = primitive.interface;
1395 element.lambda = primitive.lambda;
1396 element.bc = primitive.v;
1397 element.assigned = primitive.q;
1398 element.solution = 0.;
1399 }
1400 m_elements.insert(m_elements.end(),
1401 std::make_move_iterator(elements.begin()),
1402 std::make_move_iterator(elements.end()));
1403 }
1404
1405 // Set the user options.
1406 // Number of threads and use of primary average properties
1407 neBEM::NbThreads = m_nThreads;
1408 neBEM::PrimAfter = m_primAfter;
1409 neBEM::WtFldPrimAfter = m_wtFldPrimAfter;
1410 neBEM::OptRmPrim = m_optRmPrim;
1411
1412 // Fast volume details (physical potential and field related)
1413 neBEM::OptFastVol = m_optFastVol;
1414 neBEM::OptCreateFastPF = m_optCreateFastPF;
1415 neBEM::OptReadFastPF = m_optReadFastPF;
1416 neBEM::VersionFV = m_versionFV;
1417 neBEM::NbBlocksFV = m_nbBlocksFV;
1418 // Weighting potential and field related Fast volume details
1419 // Since MAXWtFld is not a variable, we do not need neBEM::
1420 for (int id = 1; id < MAXWtFld; ++id) {
1421 neBEM::OptWtFldFastVol[id] = m_optWtFldFastVol[id];
1422 neBEM::OptCreateWtFldFastPF[id] = m_optCreateWtFldFastPF[id];
1423 neBEM::OptReadWtFldFastPF[id] = m_optReadWtFldFastPF[id];
1424 neBEM::VersionWtFldFV[id] = m_versionWtFldFV[id];
1425 neBEM::NbBlocksWtFldFV[id] = m_nbBlocksWtFldFV[id];
1426 }
1427 // Known charge options
1428 neBEM::OptKnCh = m_optKnownCharge;
1429 // Charging up options
1430 neBEM::OptChargingUp = m_optChargingUp;
1431
1432 // Initialize neBEM
1433 if (neBEM::neBEMInitialize() != 0) {
1434 std::cerr << m_className << "::Initialise:\n"
1435 << " Initialising neBEM failed.\n";
1436 return false;
1437 }
1438 gComponentNeBem3d = this;
1439
1440 // Discretization related
1441 neBEM::MinNbElementsOnLength = m_minNbElementsOnLength;
1442 neBEM::MaxNbElementsOnLength = m_maxNbElementsOnLength;
1443 neBEM::ElementLengthRqstd = m_targetElementSize * 0.01;
1444
1445 // New model / reuse existing model flag.
1446 neBEM::NewModel = m_newModel;
1447 neBEM::NewMesh = m_newMesh;
1448 neBEM::NewBC = m_newBC;
1449 neBEM::NewPP = m_newPP;
1450
1451 // Store and read options.
1452 neBEM::OptStoreInflMatrix = m_optStoreInflMatrix;
1453 neBEM::OptReadInflMatrix = m_optReadInflMatrix;
1454 neBEM::OptStoreInvMatrix = m_optStoreInvMatrix;
1455 neBEM::OptReadInvMatrix = m_optReadInvMatrix;
1456 neBEM::OptStorePrimitives = m_optStorePrimitives;
1457 neBEM::OptReadPrimitives = m_optReadPrimitives;
1458 neBEM::OptStoreElements = m_optStoreElements;
1459 neBEM::OptReadElements = m_optReadElements;
1460 neBEM::OptFormattedFile = m_optStoreFormatted;
1461 neBEM::OptUnformattedFile = m_optStoreUnformatted;
1462 neBEM::OptSystemChargeZero = m_optSystemChargeZero;
1463 neBEM::OptValidateSolution = m_optValidateSolution;
1464 neBEM::OptForceValidation = m_optForceValidation;
1465 neBEM::OptRepeatLHMatrix = m_optRepeatLHMatrix;
1466
1467 // Compute options
1468 neBEM::OptSystemChargeZero = m_optSystemChargeZero;
1469 neBEM::OptValidateSolution = m_optValidateSolution;
1470 neBEM::OptForceValidation = m_optForceValidation;
1471 neBEM::OptRepeatLHMatrix = m_optRepeatLHMatrix;
1472
1473 // Pass debug level.
1474 neBEM::DebugLevel = m_debug ? 101 : 0;
1475
1476 // Matrix inversion method (LU or SVD).
1477 if (m_inversion == Inversion::LU) {
1478 neBEM::OptInvMatProc = 0;
1479 } else {
1480 neBEM::OptInvMatProc = 1;
1481 }
1482 // Delete existing weighting fields (if any).
1483 if (neBEM::WtFieldChDen != NULL) {
1484 neBEM::neBEMDeleteAllWeightingFields();
1485 }
1486 // Transfer the geometry.
1487 if (neBEM::neBEMReadGeometry() != 0) {
1488 std::cerr << m_className << "::Initialise:\n"
1489 << " Transferring the geometry to neBEM failed.\n";
1490 return false;
1491 }
1492
1493 // Discretization.
1494 int** elementNbs = neBEM::imatrix(1, neBEM::NbPrimitives, 1, 2);
1495 for (int i = 1; i <= neBEM::NbPrimitives; ++i) {
1496 const int vol1 = neBEM::VolRef1[i];
1497 double size1 = -1.;
1498 if (solids.find(vol1) != solids.end()) {
1499 const auto solid = solids[vol1];
1500 if (solid) {
1501 Panel panel;
1502 panel.a = neBEM::XNorm[i];
1503 panel.b = neBEM::YNorm[i];
1504 panel.c = neBEM::ZNorm[i];
1505 const int nv = neBEM::NbVertices[i];
1506 panel.xv.resize(nv);
1507 panel.yv.resize(nv);
1508 panel.zv.resize(nv);
1509 for (int j = 0; j < nv; ++j) {
1510 panel.xv[j] = neBEM::XVertex[i][j];
1511 panel.yv[j] = neBEM::YVertex[i][j];
1512 panel.zv[j] = neBEM::ZVertex[i][j];
1513 }
1514 size1 = solid->GetDiscretisationLevel(panel);
1515 }
1516 }
1517 int vol2 = neBEM::VolRef2[i];
1518 double size2 = -1.;
1519 if (solids.find(vol2) != solids.end()) {
1520 const auto solid = solids[vol2];
1521 if (solid) {
1522 Panel panel;
1523 panel.a = -neBEM::XNorm[i];
1524 panel.b = -neBEM::YNorm[i];
1525 panel.c = -neBEM::ZNorm[i];
1526 const int nv = neBEM::NbVertices[i];
1527 panel.xv.resize(nv);
1528 panel.yv.resize(nv);
1529 panel.zv.resize(nv);
1530 for (int j = 0; j < nv; ++j) {
1531 panel.xv[j] = neBEM::XVertex[i][j];
1532 panel.yv[j] = neBEM::YVertex[i][j];
1533 panel.zv[j] = neBEM::ZVertex[i][j];
1534 }
1535 size2 = solid->GetDiscretisationLevel(panel);
1536 }
1537 }
1538 double size = m_targetElementSize;
1539 if (size1 > 0. && size2 > 0.) {
1540 size = std::min(size1, size2);
1541 } else if (size1 > 0.) {
1542 size = size1;
1543 } else if (size2 > 0.) {
1544 size = size2;
1545 }
1546 // Convert from cm to m.
1547 size *= 0.01;
1548
1549 // Work out the element dimensions.
1550 const double dx1 = neBEM::XVertex[i][0] - neBEM::XVertex[i][1];
1551 const double dy1 = neBEM::YVertex[i][0] - neBEM::YVertex[i][1];
1552 const double dz1 = neBEM::ZVertex[i][0] - neBEM::ZVertex[i][1];
1553 const double l1 = dx1 * dx1 + dy1 * dy1 + dz1 * dz1;
1554 int nb1 = (int)(sqrt(l1) / size);
1555 // Truncate to desired range.
1556 if (nb1 < neBEM::MinNbElementsOnLength) {
1557 nb1 = neBEM::MinNbElementsOnLength;
1558 } else if (nb1 > neBEM::MaxNbElementsOnLength) {
1559 nb1 = neBEM::MaxNbElementsOnLength;
1560 }
1561 int nb2 = 0;
1562 if (neBEM::NbVertices[i] > 2) {
1563 const double dx2 = neBEM::XVertex[i][2] - neBEM::XVertex[i][1];
1564 const double dy2 = neBEM::YVertex[i][2] - neBEM::YVertex[i][1];
1565 const double dz2 = neBEM::ZVertex[i][2] - neBEM::ZVertex[i][1];
1566 const double l2 = dx2 * dx2 + dy2 * dy2 + dz2 * dz2;
1567 nb2 = (int)(sqrt(l2) / size);
1568 if (nb2 < neBEM::MinNbElementsOnLength) {
1569 nb2 = neBEM::MinNbElementsOnLength;
1570 } else if (nb2 > neBEM::MaxNbElementsOnLength) {
1571 nb2 = neBEM::MaxNbElementsOnLength;
1572 }
1573 }
1574 elementNbs[i][1] = nb1;
1575 elementNbs[i][2] = nb2;
1576 }
1577
1578 if (neBEM::neBEMDiscretize(elementNbs) != 0) {
1579 std::cerr << m_className << "::Initialise: Discretization failed.\n";
1580 neBEM::free_imatrix(elementNbs, 1, neBEM::NbPrimitives, 1, 2);
1581 return false;
1582 }
1583 neBEM::free_imatrix(elementNbs, 1, neBEM::NbPrimitives, 1, 2);
1584 if (neBEM::neBEMBoundaryInitialConditions() != 0) {
1585 std::cerr << m_className << "::Initialise:\n"
1586 << " Setting the boundary and initial conditions failed.\n";
1587 return false;
1588 }
1589 if (neBEM::neBEMSolve() != 0) {
1590 std::cerr << m_className << "::Initialise: Solution failed.\n";
1591 return false;
1592 }
1593 // Now the weighting fields.
1594 std::set<std::string> labels;
1595 for (unsigned int i = 0; i < nSolids; ++i) {
1596 const auto solid = m_geometry->GetSolid(i);
1597 if (!solid) continue;
1598 const std::string label = solid->GetLabel();
1599 if (!label.empty()) labels.insert(label);
1600 }
1601 for (const auto& label : labels) {
1602 std::vector<int> primitives;
1603 for (unsigned int i = 0; i < nSolids; ++i) {
1604 const auto solid = m_geometry->GetSolid(i);
1605 if (!solid) continue;
1606 if (solid->GetLabel() != label) continue;
1607 const int id = solid->GetId();
1608 // Add the primitives associated to this solid to the list.
1609 for (int j = 1; j <= neBEM::NbPrimitives; ++j) {
1610 if (neBEM::VolRef1[j] == id || neBEM::VolRef2[j] == id) {
1611 primitives.push_back(j);
1612 }
1613 }
1614 }
1615 // Request the weighting field for this list of primitives.
1616 const int np = primitives.size();
1617 const int id = neBEM::neBEMPrepareWeightingField(np, primitives.data());
1618 if (id < 0) {
1619 std::cerr << m_className << "::Initialise:\n"
1620 << " Weighting field calculation for readout group \""
1621 << label << "\" failed.\n";
1622 continue;
1623 } else {
1624 std::cout << m_className << "::Initialise:\n"
1625 << " Prepared weighting field for readout group \"" << label
1626 << "\".\n";
1627 m_wfields[label] = id;
1628 }
1629 }
1630 // TODO! Not sure if we should call this here.
1631 // neBEM::neBEMEnd();
1632 m_ready = true;
1633 return true;
1634}
Garfield::Component::m_debug
bool m_debug
Switch on/off debugging messages.
Definition Component.hh:371
Garfield::gComponentNeBem3d
ComponentNeBem3d * gComponentNeBem3d
Definition ComponentNeBem3d.cc:380
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition DoubleAc.h:615
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition DoubleAc.cpp:314
MAXWtFld
#define MAXWtFld
Definition neBEM.h:30
InterfaceType
neBEMGLOBAL int * InterfaceType
Definition neBEM.h:75

Referenced by ElectricField().

◆ Reset()

void Garfield::ComponentNeBem3d::Reset ( )
overrideprotectedvirtual

Reset the component.

Implements Garfield::Component.

Definition at line 3456 of file ComponentNeBem3d.cc.

3456 {
3457 m_primitives.clear();
3458 m_elements.clear();
3459 m_ynplan.fill(false);
3460 m_coplan.fill(0.);
3461 m_vtplan.fill(0.);
3462 m_ready = false;
3463}

◆ SetChargingUpOptions()

void Garfield::ComponentNeBem3d::SetChargingUpOptions ( const unsigned int OptChargingUp)

Definition at line 822 of file ComponentNeBem3d.cc.

822 {
823 m_optChargingUp = OptChargingUp;
824}
OptChargingUp
neBEMGLOBAL int OptChargingUp
Definition neBEM.h:60

◆ SetComputeOptions()

void Garfield::ComponentNeBem3d::SetComputeOptions ( const unsigned int OptSystemChargeZero,
const unsigned int OptValidateSolution,
const unsigned int OptForceValidation,
const unsigned int OptRepeatLHMatrix )

Definition at line 760 of file ComponentNeBem3d.cc.

763 {
764 m_optSystemChargeZero = OptSystemChargeZero;
765 m_optValidateSolution = OptValidateSolution;
766 m_optForceValidation = OptForceVaildation;
767 m_optRepeatLHMatrix = OptRepeatLHMatrix;
768}
OptRepeatLHMatrix
neBEMGLOBAL int OptRepeatLHMatrix
Definition neBEM.h:58
OptSystemChargeZero
neBEMGLOBAL int OptSystemChargeZero
Definition neBEM.h:125
OptValidateSolution
neBEMGLOBAL int OptValidateSolution
Definition neBEM.h:44

◆ SetFastVolBlocks()

void Garfield::ComponentNeBem3d::SetFastVolBlocks ( const unsigned int NbBlocksFV)

Definition at line 785 of file ComponentNeBem3d.cc.

785 {
786 m_nbBlocksFV = NbBlocksFV;
787}
NbBlocksFV
neBEMGLOBAL int NbBlocksFV
Definition neBEM.h:423

◆ SetFastVolOptions()

void Garfield::ComponentNeBem3d::SetFastVolOptions ( const unsigned int OptFastVol,
const unsigned int OptCreateFastPF,
const unsigned int OptReadFastPF )

Definition at line 771 of file ComponentNeBem3d.cc.

773 {
774 m_optFastVol = OptFastVol;
775 m_optCreateFastPF = OptCreateFastPF;
776 m_optReadFastPF = OptReadFastPF;
777}
OptCreateFastPF
neBEMGLOBAL int OptCreateFastPF
Definition neBEM.h:420
OptFastVol
neBEMGLOBAL int OptFastVol
Definition neBEM.h:418
OptReadFastPF
neBEMGLOBAL int OptReadFastPF
Definition neBEM.h:421

◆ SetFastVolVersion()

void Garfield::ComponentNeBem3d::SetFastVolVersion ( const unsigned int VersionFV)

Definition at line 780 of file ComponentNeBem3d.cc.

780 {
781 m_versionFV = VersionFV;
782}
VersionFV
neBEMGLOBAL int VersionFV
Definition neBEM.h:422

◆ SetForceValidation()

void Garfield::ComponentNeBem3d::SetForceValidation ( const unsigned int OptForceValidation)

Definition at line 751 of file ComponentNeBem3d.cc.

752 {
753 m_optForceValidation = OptForceValidation;
754}
OptForceValidation
neBEMGLOBAL int OptForceValidation
Definition neBEM.h:45

◆ SetFormattedFile()

void Garfield::ComponentNeBem3d::SetFormattedFile ( const unsigned int OptFormattedFile)

Definition at line 706 of file ComponentNeBem3d.cc.

706 {
707 m_optStoreFormatted = OptFormattedFile;
708}
OptFormattedFile
neBEMGLOBAL int OptFormattedFile
Definition neBEM.h:56

◆ SetKnownChargeOptions()

void Garfield::ComponentNeBem3d::SetKnownChargeOptions ( const unsigned int OptKnownCharge)

Definition at line 816 of file ComponentNeBem3d.cc.

817 {
818 m_optKnownCharge = OptKnownCharge;
819}

◆ 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 639 of file ComponentNeBem3d.cc.

640 {
641 if (nmin == 0 || nmax == 0) {
642 std::cerr << m_className << "::SetMinMaxNumberOfElements:\n"
643 << " Values must be non-zero.\n";
644 return;
645 }
646 m_minNbElementsOnLength = std::min(nmin, nmax);
647 m_maxNbElementsOnLength = std::max(nmin, nmax);
648}

◆ SetMirrorPeriodicityX()

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

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

Definition at line 870 of file ComponentNeBem3d.cc.

870 {
871 if (s < Small) {
872 std::cerr << m_className << "::SetMirrorPeriodicityX:\n"
873 << " Periodic length must be greater than zero.\n";
874 return;
875 }
876 m_periodicLength[0] = s;
877 m_periodic[0] = false;
878 m_mirrorPeriodic[0] = true;
879 UpdatePeriodicity();
880}
Garfield::ComponentNeBem3d::UpdatePeriodicity
void UpdatePeriodicity() override
Verify periodicities.
Definition ComponentNeBem3d.cc:3465

◆ SetMirrorPeriodicityY()

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

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

Definition at line 882 of file ComponentNeBem3d.cc.

882 {
883 if (s < Small) {
884 std::cerr << m_className << "::SetMirrorPeriodicityY:\n"
885 << " Periodic length must be greater than zero.\n";
886 return;
887 }
888 m_periodicLength[1] = s;
889 m_periodic[1] = false;
890 m_mirrorPeriodic[1] = true;
891 UpdatePeriodicity();
892}

◆ SetMirrorPeriodicityZ()

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

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

Definition at line 894 of file ComponentNeBem3d.cc.

894 {
895 if (s < Small) {
896 std::cerr << m_className << "::SetMirrorPeriodicityZ:\n"
897 << " Periodic length must be greater than zero.\n";
898 return;
899 }
900 m_periodicLength[2] = s;
901 m_periodic[2] = false;
902 m_mirrorPeriodic[2] = true;
903 UpdatePeriodicity();
904}

◆ SetModelOptions()

void Garfield::ComponentNeBem3d::SetModelOptions ( const unsigned int NewModel,
const unsigned int NewMesh,
const unsigned int NewBC,
const unsigned int NewPP )

Definition at line 662 of file ComponentNeBem3d.cc.

665 {
666 m_newModel = NewModel;
667 m_newMesh = NewMesh;
668 m_newBC = NewBC;
669 m_newPP = NewPP;
670}
NewModel
neBEMGLOBAL int NewModel
Definition neBEM.h:180
NewMesh
neBEMGLOBAL int NewMesh
Definition neBEM.h:180
NewPP
neBEMGLOBAL int NewPP
Definition neBEM.h:180
NewBC
neBEMGLOBAL int NewBC
Definition neBEM.h:180

◆ SetNewBC()

void Garfield::ComponentNeBem3d::SetNewBC ( const unsigned int NewBC)

Definition at line 658 of file ComponentNeBem3d.cc.

658{ m_newBC = NewBC; }

◆ SetNewMesh()

void Garfield::ComponentNeBem3d::SetNewMesh ( const unsigned int NewMesh)

Definition at line 654 of file ComponentNeBem3d.cc.

654 {
655 m_newMesh = NewMesh;
656}

◆ SetNewModel()

void Garfield::ComponentNeBem3d::SetNewModel ( const unsigned int NewModel)

Definition at line 650 of file ComponentNeBem3d.cc.

650 {
651 m_newModel = NewModel;
652}

◆ SetNewPP()

void Garfield::ComponentNeBem3d::SetNewPP ( const unsigned int NewPP)

Definition at line 660 of file ComponentNeBem3d.cc.

660{ m_newPP = NewPP; }

◆ SetNumberOfThreads()

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

Set the number of threads to be used by neBEM.

Definition at line 183 of file ComponentNeBem3d.hh.

183{ m_nThreads = n > 0 ? n : 1; }

◆ SetOptRmPrim()

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

Set option related to removal of primitives.

Definition at line 196 of file ComponentNeBem3d.hh.

196{ m_optRmPrim = n; }

◆ 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 826 of file ComponentNeBem3d.cc.

828 {
829 m_nCopiesX = nx;
830 m_nCopiesY = ny;
831 m_nCopiesZ = nz;
832}

◆ SetPeriodicityX()

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

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

Definition at line 834 of file ComponentNeBem3d.cc.

834 {
835 if (s < Small) {
836 std::cerr << m_className << "::SetPeriodicityX:\n"
837 << " Periodic length must be greater than zero.\n";
838 return;
839 }
840 m_periodicLength[0] = s;
841 m_periodic[0] = true;
842 m_mirrorPeriodic[0] = false;
843 UpdatePeriodicity();
844}

◆ SetPeriodicityY()

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

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

Definition at line 846 of file ComponentNeBem3d.cc.

846 {
847 if (s < Small) {
848 std::cerr << m_className << "::SetPeriodicityY:\n"
849 << " Periodic length must be greater than zero.\n";
850 return;
851 }
852 m_periodicLength[1] = s;
853 m_periodic[1] = true;
854 m_mirrorPeriodic[1] = false;
855 UpdatePeriodicity();
856}

◆ SetPeriodicityZ()

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

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

Definition at line 858 of file ComponentNeBem3d.cc.

858 {
859 if (s < Small) {
860 std::cerr << m_className << "::SetPeriodicityZ:\n"
861 << " Periodic length must be greater than zero.\n";
862 return;
863 }
864 m_periodicLength[2] = s;
865 m_periodic[2] = true;
866 m_mirrorPeriodic[2] = false;
867 UpdatePeriodicity();
868}

◆ SetPrimAfter()

void Garfield::ComponentNeBem3d::SetPrimAfter ( const int n)
inline

Set the number of repetitions after which primitive properties are used for the physical field. A negative value (default) implies all the elements are always evaluated.

Definition at line 188 of file ComponentNeBem3d.hh.

188{ m_primAfter = n; }

◆ SetReadElements()

void Garfield::ComponentNeBem3d::SetReadElements ( const unsigned int OptReadElements)

Definition at line 702 of file ComponentNeBem3d.cc.

702 {
703 m_optReadElements = OptReadElements;
704}
OptReadElements
neBEMGLOBAL int OptReadElements
Definition neBEM.h:55

◆ SetReadInflMatrix()

void Garfield::ComponentNeBem3d::SetReadInflMatrix ( const unsigned int OptReadInflMatrix)

Definition at line 677 of file ComponentNeBem3d.cc.

677 {
678 m_optReadInflMatrix = OptReadInflMatrix;
679}
OptReadInflMatrix
neBEMGLOBAL int OptReadInflMatrix
Definition neBEM.h:48

◆ SetReadInvMatrix()

void Garfield::ComponentNeBem3d::SetReadInvMatrix ( const unsigned int OptReadInvMatrix)

Definition at line 685 of file ComponentNeBem3d.cc.

685 {
686 m_optReadInvMatrix = OptReadInvMatrix;
687}
OptReadInvMatrix
neBEMGLOBAL int OptReadInvMatrix
Definition neBEM.h:50

◆ SetReadPrimitives()

void Garfield::ComponentNeBem3d::SetReadPrimitives ( const unsigned int OptReadPrimitives)

Definition at line 694 of file ComponentNeBem3d.cc.

694 {
695 m_optReadPrimitives = OptReadPrimitives;
696}
OptReadPrimitives
neBEMGLOBAL int OptReadPrimitives
Definition neBEM.h:52

◆ SetRepeatLHMatrix()

void Garfield::ComponentNeBem3d::SetRepeatLHMatrix ( const unsigned int OptRepeatLHMatrix)

Definition at line 756 of file ComponentNeBem3d.cc.

756 {
757 m_optRepeatLHMatrix = OptRepeatLHMatrix;
758}

◆ SetReuseModel()

void Garfield::ComponentNeBem3d::SetReuseModel ( void )

Definition at line 734 of file ComponentNeBem3d.cc.

734 {
735 m_newModel = 0;
736 m_newMesh = 0;
737 m_newBC = 1;
738 m_optReadInvMatrix = 1;
739}

◆ SetStoreElements()

void Garfield::ComponentNeBem3d::SetStoreElements ( const unsigned int OptStoreElements)

Definition at line 698 of file ComponentNeBem3d.cc.

698 {
699 m_optStoreElements = OptStoreElements;
700}
OptStoreElements
neBEMGLOBAL int OptStoreElements
Definition neBEM.h:54

◆ SetStoreInflMatrix()

void Garfield::ComponentNeBem3d::SetStoreInflMatrix ( const unsigned int OptStoreInflMatrix)

Set storing options (OptStoreInflMatrix, OptStoreInvMatrix, OptStoreInvMatrix, OptStoreInvMatrix) OptStorePrimitives, OptStorePrimitives) OptStoreElements, OptStoreElements) OptFormattedFile, OptUnformattedFile)

Definition at line 672 of file ComponentNeBem3d.cc.

673 {
674 m_optStoreInflMatrix = OptStoreInflMatrix;
675}
OptStoreInflMatrix
neBEMGLOBAL int OptStoreInflMatrix
Definition neBEM.h:47

◆ SetStoreInvMatrix()

void Garfield::ComponentNeBem3d::SetStoreInvMatrix ( const unsigned int OptStoreInvMatrix)

Definition at line 681 of file ComponentNeBem3d.cc.

681 {
682 m_optStoreInvMatrix = OptStoreInvMatrix;
683}
OptStoreInvMatrix
neBEMGLOBAL int OptStoreInvMatrix
Definition neBEM.h:49

◆ SetStorePrimitives()

void Garfield::ComponentNeBem3d::SetStorePrimitives ( const unsigned int OptStorePrimitives)

Definition at line 689 of file ComponentNeBem3d.cc.

690 {
691 m_optStorePrimitives = OptStorePrimitives;
692}
OptStorePrimitives
neBEMGLOBAL int OptStorePrimitives
Definition neBEM.h:51

◆ SetStoreReadOptions()

void Garfield::ComponentNeBem3d::SetStoreReadOptions ( const unsigned int OptStoreInflMatrix,
const unsigned int OptReadInflMatrix,
const unsigned int OptStoreInvMatrix,
const unsigned int OptReadInvMatrix,
const unsigned int OptStorePrimitives,
const unsigned int OptReadPrimitives,
const unsigned int OptStoreElements,
const unsigned int OptReadElements,
const unsigned int OptFormattedFile,
const unsigned int OptUnformattedFile )

Definition at line 715 of file ComponentNeBem3d.cc.

721 {
722 m_optStoreInflMatrix = OptStoreInflMatrix;
723 m_optReadInflMatrix = OptReadInflMatrix;
724 m_optStoreInvMatrix = OptStoreInvMatrix;
725 m_optReadInvMatrix = OptReadInvMatrix;
726 m_optStorePrimitives = OptStorePrimitives;
727 m_optReadPrimitives = OptReadPrimitives;
728 m_optStoreElements = OptStoreElements;
729 m_optReadElements = OptReadElements;
730 m_optStoreFormatted = OptFormattedFile;
731 m_optStoreUnformatted = OptUnformattedFile;
732}
OptUnformattedFile
neBEMGLOBAL int OptUnformattedFile
Definition neBEM.h:57

◆ SetSystemChargeZero()

void Garfield::ComponentNeBem3d::SetSystemChargeZero ( const unsigned int OptSystemChargeZero)

Other functions to be, are void SetPlotOptions(OptGnuplot=0, OptGnuplotPrimitives=0, OptGnuplotElements=0, OptPrimitiveFiles=0, OptElementFiles=0)

Definition at line 741 of file ComponentNeBem3d.cc.

742 {
743 m_optSystemChargeZero = OptSystemChargeZero;
744}

◆ 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 630 of file ComponentNeBem3d.cc.

630 {
631 if (length < MinDist) {
632 std::cerr << m_className << "::SetTargetElementSize: Value must be > "
633 << MinDist << ".\n";
634 return;
635 }
636 m_targetElementSize = length;
637}

◆ SetUnformattedFile()

void Garfield::ComponentNeBem3d::SetUnformattedFile ( const unsigned int OptUnformattedFile)

Definition at line 710 of file ComponentNeBem3d.cc.

711 {
712 m_optStoreUnformatted = OptUnformattedFile;
713}

◆ SetValidateSolution()

void Garfield::ComponentNeBem3d::SetValidateSolution ( const unsigned int OptValidateSolution)

Definition at line 746 of file ComponentNeBem3d.cc.

747 {
748 m_optValidateSolution = OptValidateSolution;
749}

◆ SetWtFldFastVolBlocks()

void Garfield::ComponentNeBem3d::SetWtFldFastVolBlocks ( const unsigned int IdWtField,
const unsigned int NbBlocksWtFldFV )

Definition at line 809 of file ComponentNeBem3d.cc.

810 {
811 m_idWtField = IdWtField;
812 m_nbBlocksWtFldFV[IdWtField] = NbBlocksWtFldFV;
813}
NbBlocksWtFldFV
neBEMGLOBAL int NbBlocksWtFldFV[MAXWtFld]
Definition neBEM.h:497

◆ SetWtFldFastVolOptions()

void Garfield::ComponentNeBem3d::SetWtFldFastVolOptions ( const unsigned int IdWtField,
const unsigned int OptWtFldFastVol,
const unsigned int OptCreateWtFldFastPF,
const unsigned int OptReadWtFldFastPF )

Definition at line 791 of file ComponentNeBem3d.cc.

794 {
795 m_idWtField = IdWtField;
796 m_optWtFldFastVol[IdWtField] = OptWtFldFastVol;
797 m_optCreateWtFldFastPF[IdWtField] = OptCreateWtFldFastPF;
798 m_optReadWtFldFastPF[IdWtField] = OptReadWtFldFastPF;
799}
OptWtFldFastVol
neBEMGLOBAL int OptWtFldFastVol[MAXWtFld]
Definition neBEM.h:492
OptCreateWtFldFastPF
neBEMGLOBAL int OptCreateWtFldFastPF[MAXWtFld]
Definition neBEM.h:494
OptReadWtFldFastPF
neBEMGLOBAL int OptReadWtFldFastPF[MAXWtFld]
Definition neBEM.h:495

◆ SetWtFldFastVolVersion()

void Garfield::ComponentNeBem3d::SetWtFldFastVolVersion ( const unsigned int IdWtField,
const unsigned int VersionWtFldFV )

Definition at line 802 of file ComponentNeBem3d.cc.

803 {
804 m_idWtField = IdWtField;
805 m_versionWtFldFV[IdWtField] = VersionWtFldFV;
806}
VersionWtFldFV
neBEMGLOBAL int VersionWtFldFV[MAXWtFld]
Definition neBEM.h:496

◆ SetWtFldPrimAfter()

void Garfield::ComponentNeBem3d::SetWtFldPrimAfter ( const int n)
inline

Set the number of repetitions after which primitive properties are used for the weighting field. A negative value (default) implies all the elements are always evaluated.

Definition at line 193 of file ComponentNeBem3d.hh.

193{ m_wtFldPrimAfter = n; }

◆ UpdatePeriodicity()

void Garfield::ComponentNeBem3d::UpdatePeriodicity ( )
overrideprotectedvirtual

Verify periodicities.

Implements Garfield::Component.

Definition at line 3465 of file ComponentNeBem3d.cc.

3465 {
3466 for (unsigned int i = 0; i < 3; ++i) {
3467 // Cannot have regular and mirror periodicity at the same time.
3468 if (m_periodic[i] && m_mirrorPeriodic[i]) {
3469 std::cerr << m_className << "::UpdatePeriodicity:\n"
3470 << " Both simple and mirror periodicity requested. Reset.\n";
3471 m_periodic[i] = false;
3472 m_mirrorPeriodic[i] = false;
3473 continue;
3474 }
3475 if ((m_periodic[i] || m_mirrorPeriodic[i]) && m_periodicLength[i] < Small) {
3476 std::cerr << m_className << "::UpdatePeriodicity:\n"
3477 << " Periodic length is not set. Reset.\n";
3478 m_periodic[i] = m_mirrorPeriodic[i] = false;
3479 }
3480 }
3481
3482 if (m_axiallyPeriodic[0] || m_axiallyPeriodic[1] || m_axiallyPeriodic[2]) {
3483 std::cerr << m_className << "::UpdatePeriodicity:\n"
3484 << " Axial periodicity is not available.\n";
3485 m_axiallyPeriodic.fill(false);
3486 }
3487
3488 if (m_rotationSymmetric[0] || m_rotationSymmetric[1] ||
3489 m_rotationSymmetric[2]) {
3490 std::cerr << m_className << "::UpdatePeriodicity:\n"
3491 << " Rotation symmetry is not available.\n";
3492 m_rotationSymmetric.fill(false);
3493 }
3494}
Garfield::Component::m_rotationSymmetric
std::array< bool, 3 > m_rotationSymmetric
Rotation symmetry around x-axis, y-axis, z-axis.
Definition Component.hh:380
Garfield::Component::m_axiallyPeriodic
std::array< bool, 3 > m_axiallyPeriodic
Axial periodicity in x, y, z.
Definition Component.hh:378

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 147 of file ComponentNeBem3d.hh.

147{ m_inversion = Inversion::LU; }

◆ UseSVDInversion()

void Garfield::ComponentNeBem3d::UseSVDInversion ( )
inline

Invert the influence matrix using singular value decomposition.

Definition at line 149 of file ComponentNeBem3d.hh.

149{ 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 458 of file ComponentNeBem3d.cc.

460 {
461 wx = wy = wz = 0.;
462 if (m_wfields.count(label) == 0) return;
463 const int id = m_wfields[label];
464 neBEM::Point3D point;
465 point.X = 0.01 * x;
466 point.Y = 0.01 * y;
467 point.Z = 0.01 * z;
468 neBEM::Vector3D field;
469 if (neBEM::neBEMWeightingField(&point, &field, id) == DBL_MAX) {
470 std::cerr << m_className << "::WeightingField: Evaluation failed.\n";
471 return;
472 }
473 wx = 0.01 * field.X;
474 wy = 0.01 * field.Y;
475 wz = 0.01 * field.Z;
476}

◆ 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 478 of file ComponentNeBem3d.cc.

480 {
481 if (m_wfields.count(label) == 0) return 0.;
482 const int id = m_wfields[label];
483 neBEM::Point3D point;
484 point.X = 0.01 * x;
485 point.Y = 0.01 * y;
486 point.Z = 0.01 * z;
487 neBEM::Vector3D field;
488 const double v = neBEM::neBEMWeightingField(&point, &field, id);
489 return v == DBL_MAX ? 0. : v;
490}
The documentation for this class was generated from the following files: