Garfield++ 4.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Garfield::ComponentTcadBase< N > Class Template Referenceabstract

Interpolation in a field map created by Sentaurus Device. More...

#include <ComponentTcadBase.hh>

+ Inheritance diagram for Garfield::ComponentTcadBase< N >:

Classes

struct  Defect
 
struct  Element
 
struct  Region
 

Public Member Functions

 ComponentTcadBase ()=delete
 Default constructor.
 
 ComponentTcadBase (const std::string &name)
 Constructor.
 
virtual ~ComponentTcadBase ()
 Destructor.
 
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
 
bool GetVoltageRange (double &vmin, double &vmax) override
 Calculate the voltage range [V].
 
bool Initialise (const std::string &gridfilename, const std::string &datafilename)
 
bool SetWeightingField (const std::string &datfile1, const std::string &datfile2, const double dv, const std::string &label)
 
bool SetWeightingFieldShift (const std::string &label, const double x, const double y, const double z)
 
void PrintRegions () const
 List all currently defined regions.
 
size_t GetNumberOfRegions () const
 Get the number of regions in the device.
 
void GetRegion (const size_t ireg, std::string &name, bool &active) const
 Get the name and "active volume" flag of a region.
 
void SetDriftRegion (const size_t ireg)
 Make a region active ("driftable").
 
void UnsetDriftRegion (const size_t ireg)
 Make a region inactive.
 
void SetMedium (const size_t ireg, Medium *m)
 Set the medium to be associated to a given region.
 
size_t GetNumberOfElements () const
 Get the number of elements in the mesh.
 
size_t GetNumberOfNodes () const
 Get the number of vertices in the mesh.
 
void EnableVelocityMap (const bool on)
 Switch use of the imported velocity map on/off.
 
bool HasVelocityMap () const override
 Does the component have velocity maps?
 
bool ElectronVelocity (const double x, const double y, const double z, double &vx, double &vy, double &vz) override
 Get the electron drift velocity.
 
bool HoleVelocity (const double x, const double y, const double z, double &vx, double &vy, double &vz) override
 Get the hole drift velocity.
 
size_t GetNumberOfDonors ()
 Get the number of donor states found in the map.
 
size_t GetNumberOfAcceptors ()
 Get the number of acceptor states found in the map.
 
bool SetDonor (const size_t donorNumber, const double exsec, const double hxsec, const double concentration)
 
bool SetAcceptor (const size_t acceptorNumber, const double exsec, const double hxsec, const double concentration)
 Set the properties of an acceptor-type defect state.
 
void EnableAttachmentMap (const bool on)
 Switch use of the imported trapping map on/off.
 
bool HasAttachmentMap () const override
 Does the component have attachment maps?
 
bool ElectronAttachment (const double x, const double y, const double z, double &eta) override
 Get the electron attachment coefficient.
 
bool HoleAttachment (const double x, const double y, const double z, double &eta) override
 Get the hole attachment coefficient.
 
bool GetElectronLifetime (const double x, const double y, const double z, double &etau) override
 
bool GetHoleLifetime (const double x, const double y, const double z, double &htau) override
 
bool GetElectronMobility (const double x, const double y, const double z, double &mob)
 
bool GetHoleMobility (const double x, const double y, const double z, double &mob)
 
- 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 MediumGetMedium (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 UpdatePeriodicity () override
 Verify periodicities.
 
void Cleanup ()
 
virtual bool Interpolate (const double x, const double y, const double z, const std::vector< double > &field, double &f)=0
 
virtual bool Interpolate (const double x, const double y, const double z, const std::vector< std::array< double, N > > &field, double &fx, double &fy, double &fz)=0
 
virtual void FillTree ()=0
 
size_t FindRegion (const std::string &name) const
 
void MapCoordinates (std::array< double, N > &x, std::array< bool, N > &mirr) const
 
bool InBoundingBox (const std::array< double, N > &x) const
 
void UpdateAttachment ()
 
bool LoadGrid (const std::string &gridfilename)
 
bool LoadData (const std::string &datafilename)
 
bool ReadDataset (std::ifstream &datafile, const std::string &dataset)
 
bool LoadWeightingField (const std::string &datafilename, std::vector< std::array< double, N > > &wf, std::vector< double > &wp)
 
- Protected Member Functions inherited from Garfield::Component
virtual void Reset ()=0
 Reset the component.
 
virtual void UpdatePeriodicity ()=0
 Verify periodicities.
 

Static Protected Member Functions

static unsigned int ElementVertices (const Element &element)
 

Protected Attributes

std::vector< Regionm_regions
 
std::vector< std::array< double, N > > m_vertices
 
std::vector< Elementm_elements
 
std::vector< double > m_epot
 
std::vector< std::array< double, N > > m_efield
 
std::vector< std::array< double, N > > m_wfield
 
std::vector< double > m_wpot
 
std::vector< std::string > m_wlabel
 
std::vector< std::array< double, 3 > > m_wshift
 
std::vector< std::array< double, N > > m_eVelocity
 
std::vector< std::array< double, N > > m_hVelocity
 
std::vector< double > m_eMobility
 
std::vector< double > m_hMobility
 
std::vector< double > m_eLifetime
 
std::vector< double > m_hLifetime
 
std::vector< std::vector< float > > m_donorOcc
 
std::vector< std::vector< float > > m_acceptorOcc
 
std::vector< double > m_eAttachment
 
std::vector< double > m_hAttachment
 
std::vector< Defectm_donors
 
std::vector< Defectm_acceptors
 
bool m_useVelocityMap = false
 
bool m_useAttachmentMap = false
 
std::array< double, 3 > m_bbMin = {{0., 0., 0.}}
 
std::array< double, 3 > m_bbMax = {{0., 0., 0.}}
 
double m_pMin = 0.
 
double m_pMax = 0.
 
- Protected Attributes inherited from Garfield::Component
std::string m_className = "Component"
 Class name.
 
Geometrym_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.
 

Static Protected Attributes

static constexpr size_t nMaxVertices = 4
 

Detailed Description

template<size_t N>
class Garfield::ComponentTcadBase< N >

Interpolation in a field map created by Sentaurus Device.

Definition at line 14 of file ComponentTcadBase.hh.

Constructor & Destructor Documentation

◆ ComponentTcadBase() [1/2]

template<size_t N>
Garfield::ComponentTcadBase< N >::ComponentTcadBase ( )
delete

Default constructor.

◆ ComponentTcadBase() [2/2]

template<size_t N>
Garfield::ComponentTcadBase< N >::ComponentTcadBase ( const std::string &  name)
inline

Constructor.

Definition at line 19 of file ComponentTcadBase.hh.

19 : Component(name) {
20 m_regions.reserve(10);
21 m_vertices.reserve(10000);
22 m_elements.reserve(10000);
23 }
std::vector< std::array< double, N > > m_vertices
std::vector< Element > m_elements
std::vector< Region > m_regions
Component()=delete
Default constructor.

◆ ~ComponentTcadBase()

template<size_t N>
virtual Garfield::ComponentTcadBase< N >::~ComponentTcadBase ( )
inlinevirtual

Destructor.

Definition at line 25 of file ComponentTcadBase.hh.

25{}

Member Function Documentation

◆ Cleanup()

template<size_t N>
void Garfield::ComponentTcadBase< N >::Cleanup
protected

Definition at line 1518 of file ComponentTcadBase.cc.

1518 {
1519 // Vertices
1520 m_vertices.clear();
1521 // Elements
1522 m_elements.clear();
1523 // Regions
1524 m_regions.clear();
1525 // Potential and electric field.
1526 m_epot.clear();
1527 m_efield.clear();
1528 // Weighting potential and field.
1529 m_wpot.clear();
1530 m_wfield.clear();
1531 m_wlabel.clear();
1532 m_wshift.clear();
1533
1534 // Other data.
1535 m_eVelocity.clear();
1536 m_hVelocity.clear();
1537 m_eMobility.clear();
1538 m_hMobility.clear();
1539 m_eLifetime.clear();
1540 m_hLifetime.clear();
1541 m_donors.clear();
1542 m_acceptors.clear();
1543 m_donorOcc.clear();
1544 m_acceptorOcc.clear();
1545 m_eAttachment.clear();
1546 m_hAttachment.clear();
1547}
std::vector< double > m_hMobility
std::vector< Defect > m_acceptors
std::vector< std::vector< float > > m_acceptorOcc
std::vector< std::string > m_wlabel
std::vector< std::array< double, 3 > > m_wshift
std::vector< std::vector< float > > m_donorOcc
std::vector< double > m_eAttachment
std::vector< std::array< double, N > > m_hVelocity
std::vector< double > m_epot
std::vector< std::array< double, N > > m_efield
std::vector< std::array< double, N > > m_eVelocity
std::vector< double > m_hLifetime
std::vector< double > m_wpot
std::vector< Defect > m_donors
std::vector< double > m_eMobility
std::vector< std::array< double, N > > m_wfield
std::vector< double > m_hAttachment
std::vector< double > m_eLifetime

◆ ElectronAttachment()

template<size_t N>
bool Garfield::ComponentTcadBase< N >::ElectronAttachment ( const double  ,
const double  ,
const double  ,
double &  eta 
)
overridevirtual

Get the electron attachment coefficient.

Reimplemented from Garfield::Component.

Definition at line 1438 of file ComponentTcadBase.cc.

1439 {
1440 Interpolate(x, y, z, m_eAttachment, eta);
1441 return true;
1442}
virtual bool Interpolate(const double x, const double y, const double z, const std::vector< double > &field, double &f)=0

◆ ElectronVelocity()

template<size_t N>
bool Garfield::ComponentTcadBase< N >::ElectronVelocity ( const double  ,
const double  ,
const double  ,
double &  vx,
double &  vy,
double &  vz 
)
overridevirtual

Get the electron drift velocity.

Reimplemented from Garfield::Component.

Definition at line 1452 of file ComponentTcadBase.cc.

1454 {
1455 return Interpolate(x, y, z, m_eVelocity, vx, vy, vz);
1456}

◆ ElementVertices()

template<size_t N>
static unsigned int Garfield::ComponentTcadBase< N >::ElementVertices ( const Element element)
inlinestaticprotected

Definition at line 226 of file ComponentTcadBase.hh.

226 {
227 return std::min(element.type + 1, 4U);
228 }

◆ EnableAttachmentMap()

template<size_t N>
void Garfield::ComponentTcadBase< N >::EnableAttachmentMap ( const bool  on)
inline

Switch use of the imported trapping map on/off.

Definition at line 107 of file ComponentTcadBase.hh.

◆ EnableVelocityMap()

template<size_t N>
void Garfield::ComponentTcadBase< N >::EnableVelocityMap ( const bool  on)

Switch use of the imported velocity map on/off.

Definition at line 376 of file ComponentTcadBase.cc.

376 {
377 m_useVelocityMap = on;
378 if (m_ready && (m_eVelocity.empty() && m_hVelocity.empty())) {
379 std::cout << m_className << "::EnableVelocityMap:\n"
380 << " Warning: current map does not include velocity data.\n";
381 }
382}
std::string m_className
Class name.
Definition: Component.hh:329
bool m_ready
Ready for use?
Definition: Component.hh:338

◆ FillTree()

template<size_t N>
virtual void Garfield::ComponentTcadBase< N >::FillTree ( )
protectedpure virtual

◆ FindRegion()

template<size_t N>
size_t Garfield::ComponentTcadBase< N >::FindRegion ( const std::string &  name) const
protected

Definition at line 1573 of file ComponentTcadBase.cc.

1573 {
1574 const auto nRegions = m_regions.size();
1575 for (size_t j = 0; j < nRegions; ++j) {
1576 if (name == m_regions[j].name) return j;
1577 }
1578 return m_regions.size();
1579}

◆ GetElectronLifetime()

template<size_t N>
bool Garfield::ComponentTcadBase< N >::GetElectronLifetime ( const double  x,
const double  y,
const double  z,
double &  etau 
)
overridevirtual

Reimplemented from Garfield::Component.

Definition at line 1466 of file ComponentTcadBase.cc.

1467 {
1468 return Interpolate(x, y, z, m_eLifetime, tau);
1469}

◆ GetElectronMobility()

template<size_t N>
bool Garfield::ComponentTcadBase< N >::GetElectronMobility ( const double  x,
const double  y,
const double  z,
double &  mob 
)

Definition at line 1478 of file ComponentTcadBase.cc.

1479 {
1480 return Interpolate(x, y, z, m_eMobility, mob);
1481}

◆ GetHoleLifetime()

template<size_t N>
bool Garfield::ComponentTcadBase< N >::GetHoleLifetime ( const double  x,
const double  y,
const double  z,
double &  htau 
)
overridevirtual

Reimplemented from Garfield::Component.

Definition at line 1472 of file ComponentTcadBase.cc.

1473 {
1474 return Interpolate(x, y, z, m_hLifetime, tau);
1475}

◆ GetHoleMobility()

template<size_t N>
bool Garfield::ComponentTcadBase< N >::GetHoleMobility ( const double  x,
const double  y,
const double  z,
double &  mob 
)

Definition at line 1484 of file ComponentTcadBase.cc.

1485 {
1486 return Interpolate(x, y, z, m_hMobility, mob);
1487}

◆ GetNumberOfAcceptors()

template<size_t N>
size_t Garfield::ComponentTcadBase< N >::GetNumberOfAcceptors ( )
inline

Get the number of acceptor states found in the map.

Definition at line 92 of file ComponentTcadBase.hh.

92{ return m_acceptors.size(); }

◆ GetNumberOfDonors()

template<size_t N>
size_t Garfield::ComponentTcadBase< N >::GetNumberOfDonors ( )
inline

Get the number of donor states found in the map.

Definition at line 90 of file ComponentTcadBase.hh.

90{ return m_donors.size(); }

◆ GetNumberOfElements()

template<size_t N>
size_t Garfield::ComponentTcadBase< N >::GetNumberOfElements ( ) const
inline

Get the number of elements in the mesh.

Definition at line 75 of file ComponentTcadBase.hh.

75{ return m_elements.size(); }

◆ GetNumberOfNodes()

template<size_t N>
size_t Garfield::ComponentTcadBase< N >::GetNumberOfNodes ( ) const
inline

Get the number of vertices in the mesh.

Definition at line 77 of file ComponentTcadBase.hh.

77{ return m_vertices.size(); }

◆ GetNumberOfRegions()

template<size_t N>
size_t Garfield::ComponentTcadBase< N >::GetNumberOfRegions ( ) const
inline

Get the number of regions in the device.

Definition at line 64 of file ComponentTcadBase.hh.

64{ return m_regions.size(); }

Referenced by main().

◆ GetRegion()

template<size_t N>
void Garfield::ComponentTcadBase< N >::GetRegion ( const size_t  ireg,
std::string &  name,
bool &  active 
) const

Get the name and "active volume" flag of a region.

Definition at line 1363 of file ComponentTcadBase.cc.

1364 {
1365 if (i >= m_regions.size()) {
1366 std::cerr << m_className << "::GetRegion: Index out of range.\n";
1367 return;
1368 }
1369 name = m_regions[i].name;
1370 active = m_regions[i].drift;
1371}

◆ GetVoltageRange()

template<size_t N>
bool Garfield::ComponentTcadBase< N >::GetVoltageRange ( double &  vmin,
double &  vmax 
)
overridevirtual

Calculate the voltage range [V].

Implements Garfield::Component.

Definition at line 269 of file ComponentTcadBase.cc.

269 {
270 if (!m_ready) return false;
271 vmin = m_pMin;
272 vmax = m_pMax;
273 return true;
274}

◆ HasAttachmentMap()

template<size_t N>
bool Garfield::ComponentTcadBase< N >::HasAttachmentMap ( ) const
inlineoverridevirtual

Does the component have attachment maps?

Reimplemented from Garfield::Component.

Definition at line 108 of file ComponentTcadBase.hh.

108 {
109 return (m_useAttachmentMap && !(m_acceptors.empty() && m_donors.empty()));
110 }

◆ HasVelocityMap()

template<size_t N>
bool Garfield::ComponentTcadBase< N >::HasVelocityMap ( ) const
inlineoverridevirtual

Does the component have velocity maps?

Reimplemented from Garfield::Component.

Definition at line 81 of file ComponentTcadBase.hh.

81 {
82 return m_useVelocityMap && !(m_eVelocity.empty() && m_hVelocity.empty());
83 }

◆ HoleAttachment()

template<size_t N>
bool Garfield::ComponentTcadBase< N >::HoleAttachment ( const double  ,
const double  ,
const double  ,
double &  eta 
)
overridevirtual

Get the hole attachment coefficient.

Reimplemented from Garfield::Component.

Definition at line 1445 of file ComponentTcadBase.cc.

1446 {
1447 Interpolate(x, y, z, m_hAttachment, eta);
1448 return true;
1449}

◆ HoleVelocity()

template<size_t N>
bool Garfield::ComponentTcadBase< N >::HoleVelocity ( const double  ,
const double  ,
const double  ,
double &  vx,
double &  vy,
double &  vz 
)
overridevirtual

Get the hole drift velocity.

Reimplemented from Garfield::Component.

Definition at line 1459 of file ComponentTcadBase.cc.

1461 {
1462 return Interpolate(x, y, z, m_hVelocity, vx, vy, vz);
1463}

◆ InBoundingBox()

template<size_t N>
bool Garfield::ComponentTcadBase< N >::InBoundingBox ( const std::array< double, N > &  x) const
inlineprotected

Definition at line 239 of file ComponentTcadBase.hh.

239 {
240 for (size_t i = 0; i < N; ++i) {
241 if (x[i] < m_bbMin[i] || x[i] > m_bbMax[i]) return false;
242 }
243 return true;
244 }
std::array< double, 3 > m_bbMax
std::array< double, 3 > m_bbMin

◆ Initialise()

template<size_t N>
bool Garfield::ComponentTcadBase< N >::Initialise ( const std::string &  gridfilename,
const std::string &  datafilename 
)

Import mesh and field map from files.

Parameters
gridfilenamename of the .grd file containing the mesh
datafilenamename of the .dat file containing the nodal solution

Definition at line 92 of file ComponentTcadBase.cc.

93 {
94
95 m_ready = false;
96 Cleanup();
97 // Import mesh data from .grd file.
98 if (!LoadGrid(gridfilename)) {
99 std::cerr << m_className << "::Initialise:\n"
100 << " Importing mesh data failed.\n";
101 Cleanup();
102 return false;
103 }
104 // Import electric field, potential and other data from .dat file.
105 if (!LoadData(datafilename)) {
106 std::cerr << m_className << "::Initialise:\n"
107 << " Importing electric field and potential failed.\n";
108 Cleanup();
109 return false;
110 }
111
112 // Find min./max. coordinates and potentials.
113 for (size_t i = 0; i < N; ++i) {
114 m_bbMax[i] = m_vertices[m_elements[0].vertex[0]][i];
115 m_bbMin[i] = m_bbMax[i];
116 }
117 const size_t nElements = m_elements.size();
118 for (size_t i = 0; i < nElements; ++i) {
119 Element& element = m_elements[i];
120 std::array<double, N> xmin = m_vertices[element.vertex[0]];
121 std::array<double, N> xmax = m_vertices[element.vertex[0]];
122 const auto nV = ElementVertices(m_elements[i]);
123 for (unsigned int j = 0; j < nV; ++j) {
124 const auto& v = m_vertices[m_elements[i].vertex[j]];
125 for (size_t k = 0; k < N; ++k) {
126 xmin[k] = std::min(xmin[k], v[k]);
127 xmax[k] = std::max(xmax[k], v[k]);
128 }
129 }
130 constexpr double tol = 1.e-6;
131 for (size_t k = 0; k < N; ++k) {
132 m_elements[i].bbMin[k] = xmin[k] - tol;
133 m_elements[i].bbMax[k] = xmax[k] + tol;
134 m_bbMin[k] = std::min(m_bbMin[k], xmin[k]);
135 m_bbMax[k] = std::max(m_bbMax[k], xmax[k]);
136 }
137 }
138 m_pMin = *std::min_element(m_epot.begin(), m_epot.end());
139 m_pMax = *std::max_element(m_epot.begin(), m_epot.end());
140
141 std::cout << m_className << "::Initialise:\n"
142 << " Available data:\n";
143 if (!m_epot.empty()) std::cout << " Electrostatic potential\n";
144 if (!m_efield.empty()) std::cout << " Electric field\n";
145 if (!m_eMobility.empty()) std::cout << " Electron mobility\n";
146 if (!m_hMobility.empty()) std::cout << " Hole mobility\n";
147 if (!m_eVelocity.empty()) std::cout << " Electron velocity\n";
148 if (!m_hVelocity.empty()) std::cout << " Hole velocity\n";
149 if (!m_eLifetime.empty()) std::cout << " Electron lifetime\n";
150 if (!m_hLifetime.empty()) std::cout << " Hole lifetime\n";
151 if (!m_donors.empty()) {
152 std::cout << " " << m_donors.size() << " donor-type traps\n";
153 }
154 if (!m_acceptors.empty()) {
155 std::cout << " " << m_acceptors.size() << " acceptor-type traps\n";
156 }
157 const std::array<std::string, 3> axes = {"x", "y", "z"};
158 std::cout << " Bounding box:\n";
159 for (size_t i = 0; i < N; ++i) {
160 std::cout << " " << m_bbMin[i] << " < " << axes[i] << " [cm] < "
161 << m_bbMax[i] << "\n";
162 }
163 std::cout << " Voltage range:\n"
164 << " " << m_pMin << " < V < " << m_pMax << "\n";
165
166 bool ok = true;
167
168 // Count the number of elements belonging to a region.
169 const auto nRegions = m_regions.size();
170 std::vector<size_t> nElementsByRegion(nRegions, 0);
171 // Keep track of elements that are not part of any region.
172 std::vector<size_t> looseElements;
173
174 // Count the different element shapes.
175 std::map<int, unsigned int> nElementsByShape;
176 if (N == 2) {
177 nElementsByShape = {{0, 0}, {1, 0}, {2, 0}, {3, 0}};
178 } else {
179 nElementsByShape = {{2, 0}, {5, 0}};
180 }
181 unsigned int nElementsOther = 0;
182
183 // Keep track of degenerate elements.
184 std::vector<size_t> degenerateElements;
185
186 for (size_t i = 0; i < nElements; ++i) {
187 const Element& element = m_elements[i];
188 if (element.region < nRegions) {
189 ++nElementsByRegion[element.region];
190 } else {
191 looseElements.push_back(i);
192 }
193 if (nElementsByShape.count(element.type) == 0) {
194 ++nElementsOther;
195 continue;
196 }
197 nElementsByShape[element.type] += 1;
198 bool degenerate = false;
199 const auto nV = ElementVertices(m_elements[i]);
200 for (unsigned int j = 0; j < nV; ++j) {
201 for (unsigned int k = j + 1; k < nV; ++k) {
202 if (element.vertex[j] == element.vertex[k]) {
203 degenerate = true;
204 break;
205 }
206 }
207 if (degenerate) break;
208 }
209 if (degenerate) {
210 degenerateElements.push_back(i);
211 }
212 }
213
214 if (!degenerateElements.empty()) {
215 std::cerr << m_className << "::Initialise:\n"
216 << " The following elements are degenerate:\n";
217 for (size_t i : degenerateElements) std::cerr << " " << i << "\n";
218 ok = false;
219 }
220
221 if (!looseElements.empty()) {
222 std::cerr << m_className << "::Initialise:\n"
223 << " The following elements are not part of any region:\n";
224 for (size_t i : looseElements) std::cerr << " " << i << "\n";
225 ok = false;
226 }
227
228 std::cout << m_className << "::Initialise:\n"
229 << " Number of regions: " << nRegions << "\n";
230 for (size_t i = 0; i < nRegions; ++i) {
231 std::cout << " " << i << ": " << m_regions[i].name << ", "
232 << nElementsByRegion[i] << " elements\n";
233 }
234
235 std::map<int, std::string> shapes = {
236 {0, "points"}, {1, "lines"}, {2, "triangles"}, {3, "rectangles"},
237 {5, "tetrahedra"}};
238
239 std::cout << " Number of elements: " << nElements << "\n";
240 for (const auto& n : nElementsByShape) {
241 if (n.second > 0) {
242 std::cout << " " << n.second << " " << shapes[n.first] << "\n";
243 }
244 }
245 if (nElementsOther > 0) {
246 std::cerr << " " << nElementsOther << " elements of unknown type.\n"
247 << " Program bug!\n";
248 m_ready = false;
249 Cleanup();
250 return false;
251 }
252
253 std::cout << " Number of vertices: " << m_vertices.size() << "\n";
254 if (!ok) {
255 m_ready = false;
256 Cleanup();
257 return false;
258 }
259
260 FillTree();
261
262 m_ready = true;
264 std::cout << m_className << "::Initialise: Initialisation finished.\n";
265 return true;
266}
bool LoadGrid(const std::string &gridfilename)
virtual void FillTree()=0
static unsigned int ElementVertices(const Element &element)
bool LoadData(const std::string &datafilename)
void UpdatePeriodicity() override
Verify periodicities.
Definition: neBEM.h:155

Referenced by main().

◆ Interpolate() [1/2]

template<size_t N>
virtual bool Garfield::ComponentTcadBase< N >::Interpolate ( const double  x,
const double  y,
const double  z,
const std::vector< double > &  field,
double &  f 
)
protectedpure virtual

◆ Interpolate() [2/2]

template<size_t N>
virtual bool Garfield::ComponentTcadBase< N >::Interpolate ( const double  x,
const double  y,
const double  z,
const std::vector< std::array< double, N > > &  field,
double &  fx,
double &  fy,
double &  fz 
)
protectedpure virtual

◆ LoadData()

template<size_t N>
bool Garfield::ComponentTcadBase< N >::LoadData ( const std::string &  datafilename)
protected

Definition at line 917 of file ComponentTcadBase.cc.

917 {
918
919 std::ifstream datafile;
920 datafile.open(filename, std::ios::in);
921 if (!datafile) {
922 std::cerr << m_className << "::LoadData:\n"
923 << " Could not open file " << filename << ".\n";
924 return false;
925 }
926 const size_t nVertices = m_vertices.size();
927 std::vector<unsigned int> fillCount(nVertices, 0);
928
929 std::array<double, N> zeroVector{};
930 // Read the file line by line.
931 std::string line;
932 while (std::getline(datafile, line)) {
933 // Strip white space from the beginning of the line.
934 ltrim(line);
935 // Find data section.
936 if (line.substr(0, 8) != "function") continue;
937 // Read type of data set.
938 const auto pEq = line.find('=');
939 if (pEq == std::string::npos) {
940 // No "=" found.
941 std::cerr << m_className << "::LoadData:\n"
942 << " Error reading file " << filename << ".\n"
943 << " Line:\n " << line << "\n";
944 return false;
945 }
946 line = line.substr(pEq + 1);
947 std::string dataset;
948 std::istringstream data;
949 data.str(line);
950 data >> dataset;
951 data.clear();
952 if (dataset == "ElectrostaticPotential") {
953 m_epot.assign(nVertices, 0.);
954 if (!ReadDataset(datafile, dataset)) {
955 m_epot.clear();
956 return false;
957 }
958 } else if (dataset == "ElectricField") {
959 m_efield.assign(nVertices, zeroVector);
960 if (!ReadDataset(datafile, dataset)) {
961 m_efield.clear();
962 return false;
963 }
964 } else if (dataset == "eDriftVelocity") {
965 m_eVelocity.assign(nVertices, zeroVector);
966 if (!ReadDataset(datafile, dataset)) {
967 m_eVelocity.clear();
968 return false;
969 }
970 } else if (dataset == "hDriftVelocity") {
971 m_hVelocity.assign(nVertices, zeroVector);
972 if (!ReadDataset(datafile, dataset)) {
973 m_hVelocity.clear();
974 return false;
975 }
976 } else if (dataset == "eMobility") {
977 m_eMobility.assign(nVertices, 0.);
978 if (!ReadDataset(datafile, dataset)) {
979 m_eMobility.clear();
980 return false;
981 }
982 } else if (dataset == "hMobility") {
983 m_hMobility.assign(nVertices, 0.);
984 if (!ReadDataset(datafile, dataset)) {
985 m_hMobility.clear();
986 return false;
987 }
988 } else if (dataset == "eLifetime") {
989 m_eLifetime.assign(nVertices, 0.);
990 if (!ReadDataset(datafile, dataset)) {
991 m_eLifetime.clear();
992 return false;
993 }
994 } else if (dataset == "hLifetime") {
995 m_hLifetime.assign(nVertices, 0.);
996 if (!ReadDataset(datafile, dataset)) {
997 m_hLifetime.clear();
998 return false;
999 }
1000 } else if (dataset.substr(0, 14) == "TrapOccupation" &&
1001 dataset.substr(17, 2) == "Do") {
1002 if (!ReadDataset(datafile, dataset)) return false;
1003 Defect donor;
1004 donor.xsece = -1.;
1005 donor.xsech = -1.;
1006 donor.conc = -1.;
1007 m_donors.push_back(donor);
1008 } else if (dataset.substr(0, 14) == "TrapOccupation" &&
1009 dataset.substr(17, 2) == "Ac") {
1010 if (!ReadDataset(datafile, dataset)) return false;
1011 Defect acceptor;
1012 acceptor.xsece = -1.;
1013 acceptor.xsech = -1.;
1014 acceptor.conc = -1.;
1015 m_acceptors.push_back(acceptor);
1016 }
1017 }
1018 if (datafile.fail() && !datafile.eof()) {
1019 std::cerr << m_className << "::LoadData:\n"
1020 << " Error reading file " << filename << "\n";
1021 return false;
1022 }
1023 return true;
1024}
bool ReadDataset(std::ifstream &datafile, const std::string &dataset)
void ltrim(std::string &line)
Definition: Utilities.hh:10

◆ LoadGrid()

template<size_t N>
bool Garfield::ComponentTcadBase< N >::LoadGrid ( const std::string &  gridfilename)
protected

Definition at line 385 of file ComponentTcadBase.cc.

385 {
386 // Open the file containing the mesh description.
387 std::ifstream gridfile;
388 gridfile.open(filename, std::ios::in);
389 if (!gridfile) {
390 std::cerr << m_className << "::LoadGrid:\n"
391 << " Could not open file " << filename << ".\n";
392 return false;
393 }
394 // Delete existing mesh information.
395 Cleanup();
396 // Count line numbers.
397 unsigned int iLine = 0;
398 // Get the number of regions.
399 size_t nRegions = 0;
400 // Read the file line by line.
401 std::string line;
402 while (std::getline(gridfile, line)) {
403 ++iLine;
404 // Strip white space from the beginning of the line.
405 ltrim(line);
406 // Find entry 'nb_regions'.
407 if (line.substr(0, 10) != "nb_regions") continue;
408 const auto pEq = line.find('=');
409 if (pEq == std::string::npos) {
410 // No "=" sign found.
411 std::cerr << m_className << "::LoadGrid:\n"
412 << " Could not read number of regions.\n";
413 return false;
414 }
415 line = line.substr(pEq + 1);
416 std::istringstream data;
417 data.str(line);
418 data >> nRegions;
419 break;
420 }
421 if (gridfile.eof()) {
422 // Reached end of file.
423 std::cerr << m_className << "::LoadGrid:\n"
424 << " Could not find entry 'nb_regions' in file\n"
425 << " " << filename << ".\n";
426 return false;
427 } else if (gridfile.fail()) {
428 // Error reading from the file.
429 PrintError(m_className + "::LoadGrid", filename, iLine);
430 return false;
431 }
432 m_regions.resize(nRegions);
433 for (size_t j = 0; j < nRegions; ++j) {
434 m_regions[j].name = "";
435 m_regions[j].drift = false;
436 m_regions[j].medium = nullptr;
437 }
438 if (m_debug) {
439 std::cout << m_className << "::LoadGrid:\n"
440 << " Found " << nRegions << " regions.\n";
441 }
442 // Get the region names.
443 while (std::getline(gridfile, line)) {
444 ++iLine;
445 ltrim(line);
446 // Find entry 'regions'.
447 if (line.substr(0, 7) != "regions") continue;
448 // Get region names (given in brackets).
449 if (!ExtractFromSquareBrackets(line)) {
450 std::cerr << m_className << "::LoadGrid:\n"
451 << " Could not read region names.\n";
452 return false;
453 }
454 std::istringstream data;
455 data.str(line);
456 for (size_t j = 0; j < nRegions; ++j) {
457 data >> m_regions[j].name;
458 data.clear();
459 // Assume by default that all regions are active.
460 m_regions[j].drift = true;
461 m_regions[j].medium = nullptr;
462 }
463 break;
464 }
465 if (gridfile.eof()) {
466 // Reached end of file.
467 std::cerr << m_className << "::LoadGrid:\n"
468 << " Could not find entry 'regions' in file\n"
469 << " " << filename << ".\n";
470 return false;
471 } else if (gridfile.fail()) {
472 // Error reading from the file.
473 PrintError(m_className + "::LoadGrid", filename, iLine);
474 return false;
475 }
476
477 // Get the vertices.
478 size_t nVertices = 0;
479 while (std::getline(gridfile, line)) {
480 ++iLine;
481 ltrim(line);
482 // Find section 'Vertices'.
483 if (line.substr(0, 8) != "Vertices") continue;
484 // Get number of vertices (given in brackets).
485 if (!ExtractFromBrackets(line)) {
486 std::cerr << m_className << "::LoadGrid:\n"
487 << " Could not read number of vertices.\n";
488 return false;
489 }
490 std::istringstream data;
491 data.str(line);
492 data >> nVertices;
493 m_vertices.resize(nVertices);
494 // Get the coordinates of every vertex.
495 for (size_t j = 0; j < nVertices; ++j) {
496 for (size_t k = 0; k < N; ++k) {
497 gridfile >> m_vertices[j][k];
498 // Change units from micron to cm.
499 m_vertices[j][k] *= 1.e-4;
500 }
501 ++iLine;
502 }
503 break;
504 }
505 if (gridfile.eof()) {
506 std::cerr << m_className << "::LoadGrid:\n"
507 << " Could not find section 'Vertices' in file\n"
508 << " " << filename << ".\n";
509 return false;
510 } else if (gridfile.fail()) {
511 PrintError(m_className + "::LoadGrid", filename, iLine);
512 return false;
513 }
514
515 // Get the "edges" (lines connecting two vertices).
516 size_t nEdges = 0;
517 // Temporary arrays for storing edge points.
518 std::vector<unsigned > edgeP1;
519 std::vector<unsigned > edgeP2;
520 while (std::getline(gridfile, line)) {
521 ++iLine;
522 ltrim(line);
523 // Find section 'Edges'.
524 if (line.substr(0, 5) != "Edges") continue;
525 // Get the number of edges (given in brackets).
526 if (!ExtractFromBrackets(line)) {
527 std::cerr << m_className << "::LoadGrid:\n"
528 << " Could not read number of edges.\n";
529 return false;
530 }
531 std::istringstream data;
532 data.str(line);
533 data >> nEdges;
534 edgeP1.resize(nEdges);
535 edgeP2.resize(nEdges);
536 // Get the indices of the two endpoints.
537 for (size_t j = 0; j < nEdges; ++j) {
538 gridfile >> edgeP1[j] >> edgeP2[j];
539 ++iLine;
540 }
541 break;
542 }
543 if (gridfile.eof()) {
544 std::cerr << m_className << "::LoadGrid:\n"
545 << " Could not find section 'Edges' in file\n"
546 << " " << filename << ".\n";
547 return false;
548 } else if (gridfile.fail()) {
549 PrintError(m_className + "::LoadGrid", filename, iLine);
550 return false;
551 }
552
553 for (size_t i = 0; i < nEdges; ++i) {
554 // Make sure the indices of the edge endpoints are not out of range.
555 if (edgeP1[i] >= nVertices || edgeP2[i] >= nVertices) {
556 std::cerr << m_className << "::LoadGrid:\n"
557 << " Vertex index of edge " << i << " out of range.\n";
558 return false;
559 }
560 // Make sure the edge is non-degenerate.
561 if (edgeP1[i] == edgeP2[i]) {
562 std::cerr << m_className << "::LoadGrid:\n"
563 << " Edge " << i << " is degenerate.\n";
564 return false;
565 }
566 }
567
568 // Get the "faces" (only for 3D maps).
569 struct Face {
570 // Indices of edges
571 int edge[4];
572 int type;
573 };
574 size_t nFaces = 0;
575 std::vector<Face> faces;
576 if (N == 3) {
577 while (std::getline(gridfile, line)) {
578 ++iLine;
579 ltrim(line);
580 // Find section 'Faces'.
581 if (line.substr(0, 5) != "Faces") continue;
582 // Get the number of faces (given in brackets).
583 if (!ExtractFromBrackets(line)) {
584 std::cerr << m_className << "::LoadGrid:\n"
585 << " Could not read number of faces.\n";
586 return false;
587 }
588 std::istringstream data;
589 data.str(line);
590 data >> nFaces;
591 faces.resize(nFaces);
592 // Get the indices of the edges constituting this face.
593 for (size_t j = 0; j < nFaces; ++j) {
594 gridfile >> faces[j].type;
595 if (faces[j].type != 3 && faces[j].type != 4) {
596 std::cerr << m_className << "::LoadGrid:\n"
597 << " Face with index " << j
598 << " has invalid number of edges, " << faces[j].type << ".\n";
599 return false;
600 }
601 for (int k = 0; k < faces[j].type; ++k) {
602 gridfile >> faces[j].edge[k];
603 }
604 }
605 iLine += nFaces - 1;
606 break;
607 }
608 if (gridfile.eof()) {
609 std::cerr << m_className << "::LoadGrid:\n"
610 << " Could not find section 'Faces' in file\n"
611 << " " << filename << ".\n";
612 return false;
613 } else if (gridfile.fail()) {
614 PrintError(m_className + "::LoadGrid", filename, iLine);
615 return false;
616 }
617 }
618
619 // Get the elements.
620 size_t nElements = 0;
621 while (std::getline(gridfile, line)) {
622 ++iLine;
623 ltrim(line);
624 // Find section 'Elements'.
625 if (line.substr(0, 8) != "Elements") continue;
626 // Get number of elements (given in brackets).
627 if (!ExtractFromBrackets(line)) {
628 std::cerr << m_className << "::LoadGrid:\n"
629 << " Could not read number of elements.\n";
630 return false;
631 }
632 std::istringstream data;
633 data.str(line);
634 data >> nElements;
635 data.clear();
636 // Resize the list of elements.
637 m_elements.resize(nElements);
638 // Get type and constituting edges of each element.
639 for (size_t j = 0; j < nElements; ++j) {
640 ++iLine;
641 unsigned int type = 0;
642 gridfile >> type;
643 if (N == 2) {
644 if (type == 0) {
645 // Point
646 unsigned int p = 0;
647 gridfile >> p;
648 // Make sure the index is not out of range.
649 if (p >= nVertices) {
650 PrintError(m_className + "::LoadGrid", filename, iLine);
651 std::cerr << " Vertex index out of range.\n";
652 return false;
653 }
654 m_elements[j].vertex[0] = p;
655 } else if (type == 1) {
656 // Line
657 for (size_t k = 0; k < 2; ++k) {
658 int p = 0;
659 gridfile >> p;
660 if (p < 0) p = -p - 1;
661 // Make sure the index is not out of range.
662 if (p >= (int)nVertices) {
663 PrintError(m_className + "::LoadGrid", filename, iLine);
664 std::cerr << " Vertex index out of range.\n";
665 return false;
666 }
667 m_elements[j].vertex[k] = p;
668 }
669 } else if (type == 2) {
670 // Triangle
671 int p0 = 0, p1 = 0, p2 = 0;
672 gridfile >> p0 >> p1 >> p2;
673 // Negative edge index means that the sequence of the two points
674 // is supposed to be inverted.
675 // The actual index is then given by "-index - 1".
676 if (p0 < 0) p0 = -p0 - 1;
677 if (p1 < 0) p1 = -p1 - 1;
678 if (p2 < 0) p2 = -p2 - 1;
679 // Make sure the indices are not out of range.
680 if (p0 >= (int)nEdges || p1 >= (int)nEdges || p2 >= (int)nEdges) {
681 PrintError(m_className + "::LoadGrid", filename, iLine);
682 std::cerr << " Edge index out of range.\n";
683 return false;
684 }
685 m_elements[j].vertex[0] = edgeP1[p0];
686 m_elements[j].vertex[1] = edgeP2[p0];
687 if (edgeP1[p1] != m_elements[j].vertex[0] &&
688 edgeP1[p1] != m_elements[j].vertex[1]) {
689 m_elements[j].vertex[2] = edgeP1[p1];
690 } else {
691 m_elements[j].vertex[2] = edgeP2[p1];
692 }
693 // Rearrange vertices such that point 0 is on the left.
694 while (m_vertices[m_elements[j].vertex[0]][0] >
695 m_vertices[m_elements[j].vertex[1]][0] ||
696 m_vertices[m_elements[j].vertex[0]][0] >
697 m_vertices[m_elements[j].vertex[2]][0]) {
698 const int tmp = m_elements[j].vertex[0];
699 m_elements[j].vertex[0] = m_elements[j].vertex[1];
700 m_elements[j].vertex[1] = m_elements[j].vertex[2];
701 m_elements[j].vertex[2] = tmp;
702 }
703 } else if (type == 3) {
704 // Rectangle
705 for (size_t k = 0; k < 4; ++k) {
706 int p = 0;
707 gridfile >> p;
708 // Make sure the index is not out of range.
709 if (p >= (int)nEdges || -p - 1 >= (int)nEdges) {
710 PrintError(m_className + "::LoadGrid", filename, iLine);
711 std::cerr << " Edge index out of range.\n";
712 return false;
713 }
714 if (p >= 0) {
715 m_elements[j].vertex[k] = edgeP1[p];
716 } else {
717 m_elements[j].vertex[k] = edgeP2[-p - 1];
718 }
719 }
720 // Rearrange vertices such that point 0 is on the left.
721 while (m_vertices[m_elements[j].vertex[0]][0] >
722 m_vertices[m_elements[j].vertex[1]][0] ||
723 m_vertices[m_elements[j].vertex[0]][0] >
724 m_vertices[m_elements[j].vertex[2]][0] ||
725 m_vertices[m_elements[j].vertex[0]][0] >
726 m_vertices[m_elements[j].vertex[3]][0]) {
727 const int tmp = m_elements[j].vertex[0];
728 m_elements[j].vertex[0] = m_elements[j].vertex[1];
729 m_elements[j].vertex[1] = m_elements[j].vertex[2];
730 m_elements[j].vertex[2] = m_elements[j].vertex[3];
731 m_elements[j].vertex[3] = tmp;
732 }
733 } else {
734 // Other element types are not permitted for 2d grids.
735 PrintError(m_className + "::LoadGrid", filename, iLine);
736 std::cerr << " Invalid element type (" << type
737 << ") for 2d mesh.\n";
738 return false;
739 }
740 } else if (N == 3) {
741 if (type == 2) {
742 // Triangle
743 int edge0, edge1, edge2;
744 gridfile >> edge0 >> edge1 >> edge2;
745 // Get the vertices.
746 // Negative edge index means that the sequence of the two points
747 // is supposed to be inverted.
748 // The actual index is then given by "-index - 1".
749 // For our purposes, the orientation does not matter.
750 if (edge0 < 0) edge0 = -edge0 - 1;
751 if (edge1 < 0) edge1 = -edge1 - 1;
752 if (edge2 < 0) edge2 = -edge2 - 1;
753 // Make sure the indices are not out of range.
754 if (edge0 >= (int)nEdges || edge1 >= (int)nEdges ||
755 edge2 >= (int)nEdges) {
756 PrintError(m_className + "::LoadGrid", filename, iLine);
757 std::cerr << " Edge index out of range.\n";
758 return false;
759 }
760 m_elements[j].vertex[0] = edgeP1[edge0];
761 m_elements[j].vertex[1] = edgeP2[edge0];
762 if (edgeP1[edge1] != m_elements[j].vertex[0] &&
763 edgeP1[edge1] != m_elements[j].vertex[1]) {
764 m_elements[j].vertex[2] = edgeP1[edge1];
765 } else {
766 m_elements[j].vertex[2] = edgeP2[edge1];
767 }
768 } else if (type == 5) {
769 // Tetrahedron
770 // Get the faces.
771 // Negative face index means that the sequence of the edges
772 // is supposed to be inverted.
773 // For our purposes, the orientation does not matter.
774 int face0, face1, face2, face3;
775 gridfile >> face0 >> face1 >> face2 >> face3;
776 if (face0 < 0) face0 = -face0 - 1;
777 if (face1 < 0) face1 = -face1 - 1;
778 if (face2 < 0) face2 = -face2 - 1;
779 if (face3 < 0) face3 = -face3 - 1;
780 // Make sure the face indices are not out of range.
781 if (face0 >= (int)nFaces || face1 >= (int)nFaces ||
782 face2 >= (int)nFaces || face3 >= (int)nFaces) {
783 PrintError(m_className + "::LoadGrid", filename, iLine);
784 std::cerr << " Face index out of range.\n";
785 return false;
786 }
787 // Get the edges of the first face.
788 int edge0 = faces[face0].edge[0];
789 int edge1 = faces[face0].edge[1];
790 int edge2 = faces[face0].edge[2];
791 if (edge0 < 0) edge0 = -edge0 - 1;
792 if (edge1 < 0) edge1 = -edge1 - 1;
793 if (edge2 < 0) edge2 = -edge2 - 1;
794 // Make sure the edge indices are not out of range.
795 if (edge0 >= (int)nEdges || edge1 >= (int)nEdges ||
796 edge2 >= (int)nEdges) {
797 PrintError(m_className + "::LoadGrid", filename, iLine);
798 std::cerr << " Edge index out of range.\n";
799 return false;
800 }
801 // Get the first three vertices.
802 m_elements[j].vertex[0] = edgeP1[edge0];
803 m_elements[j].vertex[1] = edgeP2[edge0];
804 if (edgeP1[edge1] != m_elements[j].vertex[0] &&
805 edgeP1[edge1] != m_elements[j].vertex[1]) {
806 m_elements[j].vertex[2] = edgeP1[edge1];
807 } else {
808 m_elements[j].vertex[2] = edgeP2[edge1];
809 }
810 // Get the fourth vertex from face 1.
811 edge0 = faces[face1].edge[0];
812 edge1 = faces[face1].edge[1];
813 edge2 = faces[face1].edge[2];
814 if (edge0 < 0) edge0 = -edge0 - 1;
815 if (edge1 < 0) edge1 = -edge1 - 1;
816 if (edge2 < 0) edge2 = -edge2 - 1;
817 const auto v0 = m_elements[j].vertex[0];
818 const auto v1 = m_elements[j].vertex[1];
819 const auto v2 = m_elements[j].vertex[2];
820 if (edgeP1[edge0] != v0 && edgeP1[edge0] != v1 && edgeP1[edge0] != v2) {
821 m_elements[j].vertex[3] = edgeP1[edge0];
822 } else if (edgeP2[edge0] != v0 && edgeP2[edge0] != v1 &&
823 edgeP2[edge0] != v2) {
824 m_elements[j].vertex[3] = edgeP2[edge0];
825 } else if (edgeP1[edge1] != v0 &&
826 edgeP1[edge1] != v1 &&
827 edgeP1[edge1] != v2) {
828 m_elements[j].vertex[3] = edgeP1[edge1];
829 } else if (edgeP2[edge1] != v0 &&
830 edgeP2[edge1] != v1 &&
831 edgeP2[edge1] != v2) {
832 m_elements[j].vertex[3] = edgeP2[edge1];
833 } else {
834 PrintError(m_className + "::LoadGrid", filename, iLine);
835 std::cerr << " Face 1 of element " << j << " is degenerate.\n";
836 return false;
837 }
838 } else {
839 // Other element types are not allowed.
840 PrintError(m_className + "::LoadGrid", filename, iLine);
841 if (type == 0 || type == 1) {
842 std::cerr << " Invalid element type (" << type
843 << ") for 3d mesh.\n";
844 } else {
845 std::cerr << " Element type " << type << " is not supported.\n"
846 << " Remesh with option -t to create only"
847 << " triangles and tetrahedra.\n";
848 }
849 return false;
850 }
851 }
852 m_elements[j].type = type;
853 m_elements[j].region = m_regions.size();
854 }
855 break;
856 }
857 if (gridfile.eof()) {
858 std::cerr << m_className << "::LoadGrid:\n"
859 << " Could not find section 'Elements' in file\n"
860 << " " << filename << ".\n";
861 return false;
862 } else if (gridfile.fail()) {
863 PrintError(m_className + "::LoadGrid", filename, iLine);
864 return false;
865 }
866
867 // Assign regions to elements.
868 while (std::getline(gridfile, line)) {
869 ltrim(line);
870 // Find section 'Region'.
871 if (line.substr(0, 6) != "Region") continue;
872 // Get region name (given in brackets).
873 if (!ExtractFromBrackets(line)) {
874 std::cerr << m_className << "::LoadGrid:\n"
875 << " Could not read region name.\n";
876 return false;
877 }
878 std::istringstream data;
879 data.str(line);
880 std::string name;
881 data >> name;
882 data.clear();
883 const size_t index = FindRegion(name);
884 if (index >= m_regions.size()) {
885 // Specified region name is not in the list.
886 std::cerr << m_className << "::LoadGrid:\n"
887 << " Unknown region " << name << ".\n";
888 continue;
889 }
890 std::getline(gridfile, line);
891 std::getline(gridfile, line);
892 if (!ExtractFromBrackets(line)) {
893 std::cerr << m_className << "::LoadGrid:\n"
894 << " Could not read number of elements in region "
895 << name << ".\n";
896 return false;
897 }
898 int nElementsRegion;
899 int iElement;
900 data.str(line);
901 data >> nElementsRegion;
902 data.clear();
903 for (int j = 0; j < nElementsRegion; ++j) {
904 gridfile >> iElement;
905 m_elements[iElement].region = index;
906 }
907 }
908 if (gridfile.fail() && !gridfile.eof()) {
909 std::cerr << m_className << "::LoadGrid:\n"
910 << " Error reading file " << filename << ".\n";
911 return false;
912 }
913 return true;
914}
size_t FindRegion(const std::string &name) const
bool m_debug
Switch on/off debugging messages.
Definition: Component.hh:341

◆ LoadWeightingField()

template<size_t N>
bool Garfield::ComponentTcadBase< N >::LoadWeightingField ( const std::string &  datafilename,
std::vector< std::array< double, N > > &  wf,
std::vector< double > &  wp 
)
protected

Definition at line 1200 of file ComponentTcadBase.cc.

1202 {
1203
1204 std::ifstream datafile;
1205 datafile.open(filename.c_str(), std::ios::in);
1206 if (!datafile) {
1207 std::cerr << m_className << "::LoadWeightingField:\n"
1208 << " Could not open file " << filename << ".\n";
1209 return false;
1210 }
1211 const size_t nVertices = m_vertices.size();
1212 const std::array<double, N> zeroVector{};
1213 bool ok = true;
1214 // Read the file line by line.
1215 std::string line;
1216 while (std::getline(datafile, line)) {
1217 // Strip white space from the beginning of the line.
1218 ltrim(line);
1219 // Find data section.
1220 if (line.substr(0, 8) != "function") continue;
1221 // Read type of data set.
1222 const auto pEq = line.find('=');
1223 if (pEq == std::string::npos) {
1224 // No "=" found.
1225 std::cerr << m_className << "::LoadWeightingField:\n"
1226 << " Error reading file " << filename << ".\n"
1227 << " Line:\n " << line << "\n";
1228 return false;
1229 }
1230 line = line.substr(pEq + 1);
1231 std::string dataset;
1232 std::istringstream data;
1233 data.str(line);
1234 data >> dataset;
1235 data.clear();
1236 if (dataset != "ElectrostaticPotential" && dataset != "ElectricField") {
1237 continue;
1238 }
1239 bool field = false;
1240 if (dataset == "ElectricField") {
1241 wf.assign(nVertices, zeroVector);
1242 field = true;
1243 } else {
1244 wp.assign(nVertices, 0.);
1245 }
1246 std::getline(datafile, line);
1247 std::getline(datafile, line);
1248 std::getline(datafile, line);
1249 std::getline(datafile, line);
1250 // Get the region name (given in brackets).
1251 if (!ExtractFromSquareBrackets(line)) {
1252 std::cerr << m_className << "::LoadWeightingField:\n"
1253 << " Cannot extract region name.\n"
1254 << " Line:\n " << line << "\n";
1255 ok = false;
1256 break;
1257 }
1258 std::string name;
1259 data.str(line);
1260 data >> name;
1261 data.clear();
1262 // Check if the region name matches one from the mesh file.
1263 const auto index = FindRegion(name);
1264 if (index >= m_regions.size()) {
1265 std::cerr << m_className << "::LoadWeightingField:\n"
1266 << " Unknown region " << name << ".\n";
1267 ok = false;
1268 break;
1269 }
1270 // Get the number of values.
1271 std::getline(datafile, line);
1272 if (!ExtractFromBrackets(line)) {
1273 std::cerr << m_className << "::LoadWeightingField:\n"
1274 << " Cannot extract number of values to be read.\n"
1275 << " Line:\n " << line << "\n";
1276 ok = false;
1277 break;
1278 }
1279 int nValues;
1280 data.str(line);
1281 data >> nValues;
1282 data.clear();
1283 if (field) nValues /= N;
1284 // Mark the vertices belonging to this region.
1285 std::vector<bool> isInRegion(nVertices, false);
1286 const size_t nElements = m_elements.size();
1287 for (size_t j = 0; j < nElements; ++j) {
1288 if (m_elements[j].region != index) continue;
1289 const unsigned int nV = ElementVertices(m_elements[j]);
1290 for (unsigned int k = 0; k < nV; ++k) {
1291 isInRegion[m_elements[j].vertex[k]] = true;
1292 }
1293 }
1294 unsigned int ivertex = 0;
1295 for (int j = 0; j < nValues; ++j) {
1296 // Read the next value.
1297 std::array<double, N> val;
1298 if (field) {
1299 for (size_t k = 0; k < N; ++k) datafile >> val[k];
1300 } else {
1301 datafile >> val[0];
1302 }
1303 // Find the next vertex belonging to the region.
1304 while (ivertex < nVertices) {
1305 if (isInRegion[ivertex]) break;
1306 ++ivertex;
1307 }
1308 // Check if there is a mismatch between the number of vertices
1309 // and the number of values.
1310 if (ivertex >= nVertices) {
1311 std::cerr << m_className << "::LoadWeightingField:\n"
1312 << " Dataset " << dataset
1313 << " has more values than vertices in region " << name << "\n";
1314 ok = false;
1315 break;
1316 }
1317 if (field) {
1318 wf[ivertex] = val;
1319 } else {
1320 wp[ivertex] = val[0];
1321 }
1322 ++ivertex;
1323 }
1324 }
1325
1326 if (!ok || (datafile.fail() && !datafile.eof())) {
1327 std::cerr << m_className << "::LoadWeightingField:\n"
1328 << " Error reading file " << filename << "\n";
1329 return false;
1330 }
1331 return true;
1332}

◆ MapCoordinates()

template<size_t N>
void Garfield::ComponentTcadBase< N >::MapCoordinates ( std::array< double, N > &  x,
std::array< bool, N > &  mirr 
) const
protected

Definition at line 1550 of file ComponentTcadBase.cc.

1551 {
1552 mirr.fill(false);
1553 for (size_t i = 0; i < N; ++i) {
1554 // In case of periodicity, reduce to the cell volume.
1555 const double cellsx = m_bbMax[i] - m_bbMin[i];
1556 if (m_periodic[i]) {
1557 x[i] = m_bbMin[i] + fmod(x[i] - m_bbMin[i], cellsx);
1558 if (x[i] < m_bbMin[i]) x[i] += cellsx;
1559 } else if (m_mirrorPeriodic[i]) {
1560 double xNew = m_bbMin[i] + fmod(x[i] - m_bbMin[i], cellsx);
1561 if (xNew < m_bbMin[i]) xNew += cellsx;
1562 const int nx = int(floor(0.5 + (xNew - x[i]) / cellsx));
1563 if (nx != 2 * (nx / 2)) {
1564 xNew = m_bbMin[i] + m_bbMax[i] - xNew;
1565 mirr[i] = true;
1566 }
1567 x[i] = xNew;
1568 }
1569 }
1570}
std::array< bool, 3 > m_mirrorPeriodic
Mirror periodicity in x, y, z.
Definition: Component.hh:346
std::array< bool, 3 > m_periodic
Simple periodicity in x, y, z.
Definition: Component.hh:344

◆ PrintRegions()

template<size_t N>
void Garfield::ComponentTcadBase< N >::PrintRegions

List all currently defined regions.

Definition at line 1335 of file ComponentTcadBase.cc.

1335 {
1336
1337 if (m_regions.empty()) {
1338 std::cerr << m_className << "::PrintRegions:\n"
1339 << " No regions are currently defined.\n";
1340 return;
1341 }
1342
1343 const size_t nRegions = m_regions.size();
1344 std::cout << m_className << "::PrintRegions:\n"
1345 << " Currently " << nRegions << " regions are defined.\n"
1346 << " Index Name Medium\n";
1347 for (size_t i = 0; i < nRegions; ++i) {
1348 std::cout << " " << i << " " << m_regions[i].name;
1349 if (!m_regions[i].medium) {
1350 std::cout << " none ";
1351 } else {
1352 std::cout << " " << m_regions[i].medium->GetName();
1353 }
1354 if (m_regions[i].drift) {
1355 std::cout << " (active region)\n";
1356 } else {
1357 std::cout << "\n";
1358 }
1359 }
1360}

◆ ReadDataset()

template<size_t N>
bool Garfield::ComponentTcadBase< N >::ReadDataset ( std::ifstream &  datafile,
const std::string &  dataset 
)
protected

Definition at line 1027 of file ComponentTcadBase.cc.

1028 {
1029
1030 if (!datafile.is_open()) return false;
1031 enum DataSet {
1032 ElectrostaticPotential,
1033 EField,
1034 eDriftVelocity,
1035 hDriftVelocity,
1036 eMobility,
1037 hMobility,
1038 eLifetime,
1039 hLifetime,
1040 DonorTrapOccupation,
1041 AcceptorTrapOccupation,
1042 Unknown
1043 };
1044 DataSet ds = Unknown;
1045 if (dataset == "ElectrostaticPotential") {
1046 ds = ElectrostaticPotential;
1047 } else if (dataset == "ElectricField") {
1048 ds = EField;
1049 } else if (dataset == "eDriftVelocity") {
1050 ds = eDriftVelocity;
1051 } else if (dataset == "hDriftVelocity") {
1052 ds = hDriftVelocity;
1053 } else if (dataset == "eMobility") {
1054 ds = eMobility;
1055 } else if (dataset == "hMobility") {
1056 ds = hMobility;
1057 } else if (dataset == "eLifetime") {
1058 ds = eLifetime;
1059 } else if (dataset == "hLifetime") {
1060 ds = hLifetime;
1061 } else if (dataset.substr(0, 14) == "TrapOccupation") {
1062 if (dataset.substr(17, 2) == "Do") {
1063 ds = DonorTrapOccupation;
1064 } else if (dataset.substr(17, 2) == "Ac") {
1065 ds = AcceptorTrapOccupation;
1066 }
1067 } else {
1068 std::cerr << m_className << "::ReadDataset:\n"
1069 << " Unexpected dataset " << dataset << ".\n";
1070 return false;
1071 }
1072 bool isVector = false;
1073 if (ds == EField || ds == eDriftVelocity || ds == hDriftVelocity) {
1074 isVector = true;
1075 }
1076
1077 std::string line;
1078 std::getline(datafile, line);
1079 std::getline(datafile, line);
1080 std::getline(datafile, line);
1081 std::getline(datafile, line);
1082 // Get the region name (given in brackets).
1083 if (!ExtractFromSquareBrackets(line)) {
1084 std::cerr << m_className << "::ReadDataset:\n"
1085 << " Cannot extract region name.\n"
1086 << " Line:\n " << line << "\n";
1087 return false;
1088 }
1089 std::string name;
1090 std::istringstream data;
1091 data.str(line);
1092 data >> name;
1093 data.clear();
1094 // Check if the region name matches one from the mesh file.
1095 const size_t index = FindRegion(name);
1096 if (index >= m_regions.size()) {
1097 std::cerr << m_className << "::ReadDataset:\n"
1098 << " Unknown region " << name << ".\n";
1099 return false;
1100 }
1101 // Get the number of values.
1102 std::getline(datafile, line);
1103 if (!ExtractFromBrackets(line)) {
1104 std::cerr << m_className << "::ReadDataset:\n"
1105 << " Cannot extract number of values to be read.\n"
1106 << " Line:\n " << line << "\n";
1107 return false;
1108 }
1109 int nValues;
1110 data.str(line);
1111 data >> nValues;
1112 if (isVector) nValues /= N;
1113 // Mark the vertices belonging to this region.
1114 const size_t nVertices = m_vertices.size();
1115 std::vector<bool> isInRegion(nVertices, false);
1116 const size_t nElements = m_elements.size();
1117 for (size_t j = 0; j < nElements; ++j) {
1118 if (m_elements[j].region != index) continue;
1119 const unsigned int nV = ElementVertices(m_elements[j]);
1120 for (unsigned int k = 0; k < nV; ++k) {
1121 isInRegion[m_elements[j].vertex[k]] = true;
1122 }
1123 }
1124
1125 unsigned int ivertex = 0;
1126 for (int j = 0; j < nValues; ++j) {
1127 // Read the next value.
1128 std::array<double, N> val;
1129 if (isVector) {
1130 for (size_t k = 0; k < N; ++k) datafile >> val[k];
1131 } else {
1132 datafile >> val[0];
1133 }
1134 // Find the next vertex belonging to the region.
1135 while (ivertex < nVertices) {
1136 if (isInRegion[ivertex]) break;
1137 ++ivertex;
1138 }
1139 // Check if there is a mismatch between the number of m_vertices
1140 // and the number of potential values.
1141 if (ivertex >= nVertices) {
1142 std::cerr << m_className << "::ReadDataset:\n"
1143 << " Dataset " << dataset << " has more values than "
1144 << "there are vertices in region " << name << "\n";
1145 return false;
1146 }
1147
1148 switch (ds) {
1149 case ElectrostaticPotential:
1150 m_epot[ivertex] = val[0];
1151 break;
1152 case EField:
1153 for (size_t k = 0; k < N; ++k) m_efield[ivertex][k] = val[k];
1154 break;
1155 case eDriftVelocity:
1156 // Scale from cm/s to cm/ns.
1157 for (size_t k = 0; k < N; ++k) {
1158 m_eVelocity[ivertex][k] = val[k] * 1.e-9;
1159 }
1160 break;
1161 case hDriftVelocity:
1162 // Scale from cm/s to cm/ns.
1163 for (size_t k = 0; k < N; ++k) {
1164 m_hVelocity[ivertex][k] = val[k] * 1.e-9;
1165 }
1166 break;
1167 case eMobility:
1168 // Convert from cm2 / (V s) to cm2 / (V ns).
1169 m_eMobility[ivertex] = val[0] * 1.e-9;
1170 break;
1171 case hMobility:
1172 // Convert from cm2 / (V s) to cm2 / (V ns).
1173 m_hMobility[ivertex] = val[0] * 1.e-9;
1174 break;
1175 case eLifetime:
1176 // Convert from s to ns.
1177 m_eLifetime[ivertex] = val[0] * 1.e9;
1178 break;
1179 case hLifetime:
1180 // Convert from s to ns.
1181 m_hLifetime[ivertex] = val[0] * 1.e9;
1182 break;
1183 case DonorTrapOccupation:
1184 m_donorOcc[ivertex].push_back(val[0]);
1185 break;
1186 case AcceptorTrapOccupation:
1187 m_acceptorOcc[ivertex].push_back(val[0]);
1188 break;
1189 default:
1190 std::cerr << m_className << "::ReadDataset:\n"
1191 << " Unexpected dataset (" << ds << "). Program bug!\n";
1192 return false;
1193 }
1194 ++ivertex;
1195 }
1196 return true;
1197}

◆ SetAcceptor()

template<size_t N>
bool Garfield::ComponentTcadBase< N >::SetAcceptor ( const size_t  acceptorNumber,
const double  exsec,
const double  hxsec,
const double  concentration 
)

Set the properties of an acceptor-type defect state.

Definition at line 1422 of file ComponentTcadBase.cc.

1424 {
1425 if (acceptorNumber >= m_acceptors.size()) {
1426 std::cerr << m_className << "::SetAcceptor: Index out of range.\n";
1427 return false;
1428 }
1429 m_acceptors[acceptorNumber].xsece = eXsec;
1430 m_acceptors[acceptorNumber].xsech = hXsec;
1431 m_acceptors[acceptorNumber].conc = conc;
1432
1434 return true;
1435}

◆ SetDonor()

template<size_t N>
bool Garfield::ComponentTcadBase< N >::SetDonor ( const size_t  donorNumber,
const double  exsec,
const double  hxsec,
const double  concentration 
)

Set the properties of a donor-type defect state.

Parameters
donorNumberindex of the donor
exseccross-section [cm2] for electrons
hxseccross-section [cm2] for holes
concentrationdefect density [cm-3]

Definition at line 1406 of file ComponentTcadBase.cc.

1408 {
1409 if (donorNumber >= m_donors.size()) {
1410 std::cerr << m_className << "::SetDonor: Index out of range.\n";
1411 return false;
1412 }
1413 m_donors[donorNumber].xsece = eXsec;
1414 m_donors[donorNumber].xsech = hXsec;
1415 m_donors[donorNumber].conc = conc;
1416
1418 return true;
1419}

◆ SetDriftRegion()

template<size_t N>
void Garfield::ComponentTcadBase< N >::SetDriftRegion ( const size_t  ireg)

Make a region active ("driftable").

Definition at line 1374 of file ComponentTcadBase.cc.

1374 {
1375 if (i >= m_regions.size()) {
1376 std::cerr << m_className << "::SetDriftRegion: Index out of range.\n";
1377 return;
1378 }
1379 m_regions[i].drift = true;
1380}

◆ SetMedium()

template<size_t N>
void Garfield::ComponentTcadBase< N >::SetMedium ( const size_t  ireg,
Medium m 
)

Set the medium to be associated to a given region.

Definition at line 1392 of file ComponentTcadBase.cc.

1392 {
1393 if (i >= m_regions.size()) {
1394 std::cerr << m_className << "::SetMedium: Index out of range.\n";
1395 return;
1396 }
1397
1398 if (!medium) {
1399 std::cerr << m_className << "::SetMedium: Null pointer.\n";
1400 return;
1401 }
1402 m_regions[i].medium = medium;
1403}

Referenced by main().

◆ SetWeightingField()

template<size_t N>
bool Garfield::ComponentTcadBase< N >::SetWeightingField ( const std::string &  datfile1,
const std::string &  datfile2,
const double  dv,
const std::string &  label 
)

Import field maps defining the weighting field and potential.

Parameters
datfile1.dat file containing the field map at nominal bias.
datfile2.dat file containing the field map for a configuration with the potential at the electrode to be read out increased by a small voltage dv.
dvincrease in electrode potential between the two field maps.
labelname of the electrode

The field maps must use the same mesh as the drift field.

Definition at line 277 of file ComponentTcadBase.cc.

280 {
281
282 if (!m_ready) {
283 std::cerr << m_className << "::SetWeightingField:\n"
284 << " Mesh is not available. Call Initialise first.\n";
285 return false;
286 }
287 if (dv < Small) {
288 std::cerr << m_className << "::SetWeightingField:\n"
289 << " Voltage difference must be > 0.\n";
290 return false;
291 }
292 const double s = 1. / dv;
293 m_wfield.clear();
294 m_wpot.clear();
295 m_wlabel.clear();
296 m_wshift.clear();
297
298 // Load first the field/potential at nominal bias.
299 std::vector<std::array<double, N> > wf1;
300 std::vector<double> wp1;
301 if (!LoadWeightingField(datfile1, wf1, wp1)) {
302 std::cerr << m_className << "::SetWeightingField:\n"
303 << " Could not import data from " << datfile1 << ".\n";
304 return false;
305 }
306
307 // Then load the field/potential for the configuration with the potential
308 // at the electrode to be read out increased by small voltage dv.
309 std::vector<std::array<double, N> > wf2;
310 std::vector<double> wp2;
311 if (!LoadWeightingField(datfile2, wf2, wp2)) {
312 std::cerr << m_className << "::SetWeightingField:\n"
313 << " Could not import data from " << datfile2 << ".\n";
314 return false;
315 }
316 const size_t nVertices = m_vertices.size();
317 bool foundField = true;
318 if (wf1.size() != nVertices || wf2.size() != nVertices) {
319 foundField = false;
320 std::cerr << m_className << "::SetWeightingField:\n"
321 << " Could not load electric field values.\n";
322 }
323 bool foundPotential = true;
324 if (wp1.size() != nVertices || wp2.size() != nVertices) {
325 foundPotential = false;
326 std::cerr << m_className << "::SetWeightingField:\n"
327 << " Could not load electrostatic potentials.\n";
328 }
329 if (!foundField && !foundPotential) return false;
330 if (foundField) {
331 m_wfield.resize(nVertices);
332 for (size_t i = 0; i < nVertices; ++i) {
333 for (size_t j = 0; j < N; ++j) {
334 m_wfield[i][j] = (wf2[i][j] - wf1[i][j]) * s;
335 }
336 }
337 }
338 if (foundPotential) {
339 m_wpot.assign(nVertices, 0.);
340 for (size_t i = 0; i < nVertices; ++i) {
341 m_wpot[i] = (wp2[i] - wp1[i]) * s;
342 }
343 }
344 m_wlabel.push_back(label);
345 m_wshift.push_back({0., 0., 0.});
346 return true;
347}
bool LoadWeightingField(const std::string &datafilename, std::vector< std::array< double, N > > &wf, std::vector< double > &wp)

◆ SetWeightingFieldShift()

template<size_t N>
bool Garfield::ComponentTcadBase< N >::SetWeightingFieldShift ( const std::string &  label,
const double  x,
const double  y,
const double  z 
)

Shift the maps of weighting field/potential for a given electrode with respect to the original mesh. If the electrode does not exist yet, a new one will be added to the list.

Definition at line 350 of file ComponentTcadBase.cc.

351 {
352 if (m_wlabel.empty()) {
353 std::cerr << m_className << "::SetWeightingFieldShift:\n"
354 << " No map of weighting potentials/fields loaded.\n";
355 return false;
356 }
357 const size_t n = m_wlabel.size();
358 for (size_t i = 0; i < n; ++i) {
359 if (m_wlabel[i] == label) {
360 m_wshift[i] = {x, y, z};
361 std::cout << m_className << "::SetWeightingFieldShift:\n"
362 << " Changing offset of electrode \'" << label
363 << "\' to (" << x << ", " << y << ", " << z << ") cm.\n";
364 return true;
365 }
366 }
367 m_wlabel.push_back(label);
368 m_wshift.push_back({x, y, z});
369 std::cout << m_className << "::SetWeightingFieldShift:\n"
370 << " Adding electrode \'" << label << "\' with offset ("
371 << x << ", " << y << ", " << z << ") cm.\n";
372 return true;
373}

◆ UnsetDriftRegion()

template<size_t N>
void Garfield::ComponentTcadBase< N >::UnsetDriftRegion ( const size_t  ireg)

Make a region inactive.

Definition at line 1383 of file ComponentTcadBase.cc.

1383 {
1384 if (i >= m_regions.size()) {
1385 std::cerr << m_className << "::UnsetDriftRegion: Index out of range.\n";
1386 return;
1387 }
1388 m_regions[i].drift = false;
1389}

◆ UpdateAttachment()

template<size_t N>
void Garfield::ComponentTcadBase< N >::UpdateAttachment
protected

Definition at line 1582 of file ComponentTcadBase.cc.

1582 {
1583
1584 if (m_vertices.empty()) return;
1585 const size_t nVertices = m_vertices.size();
1586 m_eAttachment.assign(nVertices, 0.);
1587 m_hAttachment.assign(nVertices, 0.);
1588
1589 const size_t nAcceptors = m_acceptors.size();
1590 for (size_t i = 0; i < nAcceptors; ++i) {
1591 const auto& defect = m_acceptors[i];
1592 if (defect.conc < 0.) continue;
1593 for (size_t j = 0; j < nVertices; ++j) {
1594 // Get the occupation probability.
1595 const double f = m_acceptorOcc[j][i];
1596 if (defect.xsece > 0.) {
1597 m_eAttachment[j] += defect.conc * defect.xsece * (1. - f);
1598 }
1599 if (defect.xsech > 0.) {
1600 m_hAttachment[j] += defect.conc * defect.xsech * f;
1601 }
1602 }
1603 }
1604 const size_t nDonors = m_donors.size();
1605 for (size_t i = 0; i < nDonors; ++i) {
1606 const auto& defect = m_donors[i];
1607 if (defect.conc < 0.) continue;
1608 for (size_t j = 0; j < nVertices; ++j) {
1609 const double f = m_donorOcc[j][i];
1610 if (defect.xsece > 0.) {
1611 m_eAttachment[j] += defect.conc * defect.xsece * f;
1612 }
1613 if (defect.xsech > 0.) {
1614 m_hAttachment[j] += defect.conc * defect.xsech * (1. - f);
1615 }
1616 }
1617 }
1618}

◆ UpdatePeriodicity()

template<size_t N>
void Garfield::ComponentTcadBase< N >::UpdatePeriodicity ( )
overrideprotectedvirtual

Verify periodicities.

Implements Garfield::Component.

Definition at line 1490 of file ComponentTcadBase.cc.

1490 {
1491 if (!m_ready) {
1492 std::cerr << m_className << "::UpdatePeriodicity:\n"
1493 << " Field map not available.\n";
1494 return;
1495 }
1496
1497 // Check for conflicts.
1498 for (size_t i = 0; i < 3; ++i) {
1499 if (m_periodic[i] && m_mirrorPeriodic[i]) {
1500 std::cerr << m_className << "::UpdatePeriodicity:\n"
1501 << " Both simple and mirror periodicity requested. Reset.\n";
1502 m_periodic[i] = m_mirrorPeriodic[i] = false;
1503 }
1504 if (m_axiallyPeriodic[i]) {
1505 std::cerr << m_className << "::UpdatePeriodicity:\n"
1506 << " Axial symmetry is not supported. Reset.\n";
1507 m_axiallyPeriodic.fill(false);
1508 }
1509 if (m_rotationSymmetric[i]) {
1510 std::cerr << m_className << "::UpdatePeriodicity:\n"
1511 << " Rotation symmetry is not supported. Reset.\n";
1512 m_rotationSymmetric.fill(false);
1513 }
1514 }
1515}
std::array< bool, 3 > m_rotationSymmetric
Rotation symmetry around x-axis, y-axis, z-axis.
Definition: Component.hh:350
std::array< bool, 3 > m_axiallyPeriodic
Axial periodicity in x, y, z.
Definition: Component.hh:348

◆ WeightingField()

template<size_t N>
void Garfield::ComponentTcadBase< N >::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 48 of file ComponentTcadBase.cc.

50 {
51 wx = wy = wz = 0.;
52 if (m_wfield.empty()) {
53 std::cerr << m_className << "::WeightingField: Not available.\n";
54 return;
55 }
56 const size_t n = m_wlabel.size();
57 for (size_t i = 0; i < n; ++i) {
58 if (label == m_wlabel[i]) {
59 const double x = xin - m_wshift[i][0];
60 const double y = yin - m_wshift[i][1];
61 const double z = zin - m_wshift[i][2];
62 Interpolate(x, y, z, m_wfield, wx, wy, wz);
63 return;
64 }
65 }
66}

◆ WeightingPotential()

template<size_t N>
double Garfield::ComponentTcadBase< N >::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 69 of file ComponentTcadBase.cc.

71 {
72
73 if (m_wpot.empty()) {
74 std::cerr << m_className << "::WeightingPotential: Not available.\n";
75 return 0.;
76 }
77 const size_t n = m_wlabel.size();
78 for (size_t i = 0; i < n; ++i) {
79 if (label == m_wlabel[i]) {
80 const double x = xin - m_wshift[i][0];
81 const double y = yin - m_wshift[i][1];
82 const double z = zin - m_wshift[i][2];
83 double v = 0.;
84 Interpolate(x, y, z, m_wpot, v);
85 return v;
86 }
87 }
88 return 0.;
89}

Member Data Documentation

◆ m_acceptorOcc

template<size_t N>
std::vector<std::vector<float> > Garfield::ComponentTcadBase< N >::m_acceptorOcc
protected

Definition at line 193 of file ComponentTcadBase.hh.

◆ m_acceptors

template<size_t N>
std::vector<Defect> Garfield::ComponentTcadBase< N >::m_acceptors
protected

◆ m_bbMax

template<size_t N>
std::array<double, 3> Garfield::ComponentTcadBase< N >::m_bbMax = {{0., 0., 0.}}
protected

◆ m_bbMin

template<size_t N>
std::array<double, 3> Garfield::ComponentTcadBase< N >::m_bbMin = {{0., 0., 0.}}
protected

◆ m_donorOcc

template<size_t N>
std::vector<std::vector<float> > Garfield::ComponentTcadBase< N >::m_donorOcc
protected

Definition at line 192 of file ComponentTcadBase.hh.

◆ m_donors

template<size_t N>
std::vector<Defect> Garfield::ComponentTcadBase< N >::m_donors
protected

◆ m_eAttachment

template<size_t N>
std::vector<double> Garfield::ComponentTcadBase< N >::m_eAttachment
protected

Definition at line 195 of file ComponentTcadBase.hh.

◆ m_efield

template<size_t N>
std::vector<std::array<double, N> > Garfield::ComponentTcadBase< N >::m_efield
protected

Definition at line 173 of file ComponentTcadBase.hh.

◆ m_elements

template<size_t N>
std::vector<Element> Garfield::ComponentTcadBase< N >::m_elements
protected

◆ m_eLifetime

template<size_t N>
std::vector<double> Garfield::ComponentTcadBase< N >::m_eLifetime
protected

Definition at line 189 of file ComponentTcadBase.hh.

◆ m_eMobility

template<size_t N>
std::vector<double> Garfield::ComponentTcadBase< N >::m_eMobility
protected

Definition at line 186 of file ComponentTcadBase.hh.

◆ m_epot

template<size_t N>
std::vector<double> Garfield::ComponentTcadBase< N >::m_epot
protected

Definition at line 171 of file ComponentTcadBase.hh.

◆ m_eVelocity

template<size_t N>
std::vector<std::array<double, N> > Garfield::ComponentTcadBase< N >::m_eVelocity
protected

◆ m_hAttachment

template<size_t N>
std::vector<double> Garfield::ComponentTcadBase< N >::m_hAttachment
protected

Definition at line 196 of file ComponentTcadBase.hh.

◆ m_hLifetime

template<size_t N>
std::vector<double> Garfield::ComponentTcadBase< N >::m_hLifetime
protected

Definition at line 190 of file ComponentTcadBase.hh.

◆ m_hMobility

template<size_t N>
std::vector<double> Garfield::ComponentTcadBase< N >::m_hMobility
protected

Definition at line 187 of file ComponentTcadBase.hh.

◆ m_hVelocity

template<size_t N>
std::vector<std::array<double, N> > Garfield::ComponentTcadBase< N >::m_hVelocity
protected

◆ m_pMax

template<size_t N>
double Garfield::ComponentTcadBase< N >::m_pMax = 0.
protected

Definition at line 220 of file ComponentTcadBase.hh.

◆ m_pMin

template<size_t N>
double Garfield::ComponentTcadBase< N >::m_pMin = 0.
protected

Definition at line 219 of file ComponentTcadBase.hh.

◆ m_regions

template<size_t N>
std::vector<Region> Garfield::ComponentTcadBase< N >::m_regions
protected

◆ m_useAttachmentMap

template<size_t N>
bool Garfield::ComponentTcadBase< N >::m_useAttachmentMap = false
protected

◆ m_useVelocityMap

template<size_t N>
bool Garfield::ComponentTcadBase< N >::m_useVelocityMap = false
protected

◆ m_vertices

template<size_t N>
std::vector<std::array<double, N> > Garfield::ComponentTcadBase< N >::m_vertices
protected

◆ m_wfield

template<size_t N>
std::vector<std::array<double, N> > Garfield::ComponentTcadBase< N >::m_wfield
protected

Definition at line 176 of file ComponentTcadBase.hh.

◆ m_wlabel

template<size_t N>
std::vector<std::string> Garfield::ComponentTcadBase< N >::m_wlabel
protected

Definition at line 179 of file ComponentTcadBase.hh.

◆ m_wpot

template<size_t N>
std::vector<double> Garfield::ComponentTcadBase< N >::m_wpot
protected

Definition at line 177 of file ComponentTcadBase.hh.

◆ m_wshift

template<size_t N>
std::vector<std::array<double, 3> > Garfield::ComponentTcadBase< N >::m_wshift
protected

Definition at line 180 of file ComponentTcadBase.hh.

◆ nMaxVertices

template<size_t N>
constexpr size_t Garfield::ComponentTcadBase< N >::nMaxVertices = 4
staticconstexprprotected

Definition at line 128 of file ComponentTcadBase.hh.


The documentation for this class was generated from the following files: