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::ComponentCST Class Reference

#include <ComponentCST.hh>

+ Inheritance diagram for Garfield::ComponentCST:

Public Member Functions

 ComponentCST ()
 Constructor.
 
 ~ComponentCST ()
 Destructor.
 
void ShiftComponent (const double xShift, const double yShift, const double zShift)
 
void GetNumberOfMeshLines (unsigned int &nx, unsigned int &ny, unsigned int &nz) const
 
size_t GetNumberOfElements () const override
 Return the number of mesh elements.
 
bool GetElement (const size_t i, size_t &mat, bool &drift, std::vector< size_t > &nodes) const override
 Return the material and node indices of a mesh element.
 
size_t GetNumberOfNodes () const override
 
bool GetNode (const size_t i, double &x, double &y, double &z) const override
 
void GetElementBoundaries (unsigned int element, double &xmin, double &xmax, double &ymin, double &ymax, double &zmin, double &zmax) const
 
Medium * GetMedium (const double x, const double y, const double z) override
 Get the medium at a given location (x, y, z).
 
void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status) override
 
void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status) override
 Calculate the drift field [V/cm] and potential [V] at (x, y, z).
 
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 Initialise (std::string elist, std::string nlist, std::string mplist, std::string prnsol, std::string unit="cm")
 
bool Initialise (std::string dataFile, std::string unit="cm")
 
bool SetWeightingField (std::string prnsol, std::string label, bool isBinary=true)
 
void SetRangeZ (const double zmin, const double zmax)
 
void DisableXField ()
 
void DisableYField ()
 
void DisableZField ()
 
void EnableShaping ()
 
void DisableShaping ()
 
int Index2Element (const unsigned int i, const unsigned int j, const unsigned int k) const
 
bool Coordinate2Index (const double x, const double y, const double z, unsigned int &i, unsigned int &j, unsigned int &k) const
 
- Public Member Functions inherited from Garfield::ComponentFieldMap
 ComponentFieldMap ()=delete
 Default constructor.
 
 ComponentFieldMap (const std::string &name)
 Constructor.
 
virtual ~ComponentFieldMap ()
 Destructor.
 
bool Check ()
 Check element aspect ratio.
 
void PrintRange ()
 Show x, y, z, V and angular ranges.
 
void PrintMaterials ()
 List all currently defined materials.
 
void DriftMedium (const size_t imat)
 Flag a field map material as a drift medium.
 
void NotDriftMedium (const size_t imat)
 Flag a field map materials as a non-drift medium.
 
size_t GetNumberOfMaterials () const
 Return the number of materials in the field map.
 
double GetPermittivity (const size_t imat) const
 Return the relative permittivity of a field map material.
 
double GetConductivity (const size_t imat) const
 Return the conductivity of a field map material.
 
void SetMedium (const size_t imat, Medium *medium)
 Associate a field map material with a Medium object.
 
Medium * GetMedium (const size_t imat) const
 Return the Medium associated to a field map material.
 
void SetGas (Medium *medium)
 
bool GetElement (const size_t i, double &vol, double &dmin, double &dmax) const
 Return the volume and aspect ratio of a mesh element.
 
double GetPotential (const size_t i) const
 
void EnableCheckMapIndices (const bool on=true)
 
void EnableDeleteBackgroundElements (const bool on=true)
 Option to eliminate mesh elements in conductors (default: on).
 
void EnableTetrahedralTreeForElementSearch (const bool on=true)
 
void EnableConvergenceWarnings (const bool on=true)
 
Medium * GetMedium (const double x, const double y, const double z) override
 Get the medium at a given location (x, y, z).
 
void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status) override
 
void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status) override
 Calculate the drift field [V/cm] and potential [V] at (x, y, z).
 
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
 
double DelayedWeightingPotential (double x, double y, double z, const double t, const std::string &label) override
 
bool IsInBoundingBox (const double x, const double y, const double z) const
 
bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax) override
 Get the bounding box coordinates.
 
bool GetElementaryCell (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax) override
 Get the coordinates of the elementary cell.
 
bool GetVoltageRange (double &vmin, double &vmax) override
 Calculate the voltage range [V].
 
void CopyWeightingPotential (const std::string &label, const std::string &labelSource, const double x, const double y, const double z, const double alpha, const double beta, const double gamma)
 
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 void MagneticField (const double x, const double y, const double z, double &bx, double &by, double &bz, int &status)
 
void SetMagneticField (const double bx, const double by, const double bz)
 Set a constant magnetic field.
 
virtual bool IsReady ()
 Ready for use?
 
double IntegrateFluxCircle (const double xc, const double yc, const double r, const unsigned int nI=50)
 
double IntegrateFluxSphere (const double xc, const double yc, const double zc, const double r, const unsigned int nI=20)
 
double IntegrateFluxParallelogram (const double x0, const double y0, const double z0, const double dx1, const double dy1, const double dz1, const double dx2, const double dy2, const double dz2, const unsigned int nU=20, const unsigned int nV=20)
 
double IntegrateWeightingFluxParallelogram (const std::string &label, const double x0, const double y0, const double z0, const double dx1, const double dy1, const double dz1, const double dx2, const double dy2, const double dz2, const unsigned int nU=20, const unsigned int nV=20)
 
double IntegrateFluxLine (const double x0, const double y0, const double z0, const double x1, const double y1, const double z1, const double xp, const double yp, const double zp, const unsigned int nI, const int isign=0)
 
virtual bool 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 SetRange () override
 
double GetElementVolume (const size_t i) const override
 
void GetAspectRatio (const size_t i, double &dmin, double &dmax) const override
 
bool Coordinate2Index (const double x, const double y, const double z, unsigned int &i, unsigned int &j, unsigned int &k, double *position_mapped, bool *mirrored) const
 
- Protected Member Functions inherited from Garfield::ComponentFieldMap
void Reset () override
 Reset the component.
 
void Prepare ()
 
void UpdatePeriodicity () override
 Verify periodicities.
 
void UpdatePeriodicity2d ()
 
void UpdatePeriodicityCommon ()
 
bool SetDefaultDriftMedium ()
 Find lowest epsilon, check for eps = 0, set default drift media flags.
 
int Field (const double x, const double y, const double z, double &fx, double &fy, double &fz, int &iel, const std::vector< double > &potentials) const
 Compute the electric/weighting field.
 
double Potential (const double x, const double y, const double z, const std::vector< double > &potentials) const
 Compute the electrostatic/weighting potential.
 
int FindElement5 (const double x, const double y, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det) const
 Find the element for a point in curved quadratic quadrilaterals.
 
int FindElement13 (const double x, const double y, const double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det) const
 Find the element for a point in curved quadratic tetrahedra.
 
int FindElementCube (const double x, const double y, const double z, double &t1, double &t2, double &t3, TMatrixD *&jac, std::vector< TMatrixD * > &dN) const
 Find the element for a point in a cube.
 
void MapCoordinates (double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
 Move (xpos, ypos, zpos) to field map coordinates.
 
void UnmapFields (double &ex, double &ey, double &ez, const double xpos, const double ypos, const double zpos, const bool xmirrored, const bool ymirrored, const bool zmirrored, const double rcoordinate, const double rotation) const
 Move (ex, ey, ez) to global coordinates.
 
void PrintWarning (const std::string &header)
 
void PrintNotReady (const std::string &header) const
 
void PrintCouldNotOpen (const std::string &header, const std::string &filename) const
 
void PrintElement (const std::string &header, const double x, const double y, const double z, const double t1, const double t2, const double t3, const double t4, const size_t i, const std::vector< double > &potential) const
 
void TimeInterpolation (const double t, double &f0, double &f1, int &i0, int &i1)
 Interpolation of potential between two time slices.
 

Additional Inherited Members

- Protected Types inherited from Garfield::ComponentFieldMap
enum class  ElementType { Unknown = 0 , Serendipity = 5 , CurvedTetrahedron = 13 }
 
- Static Protected Member Functions inherited from Garfield::ComponentFieldMap
static double ScalingFactor (std::string unit)
 
static double Potential3 (const std::array< double, 6 > &v, const std::array< double, 3 > &t)
 Interpolate the potential in a triangle.
 
static void Field3 (const std::array< double, 6 > &v, const std::array< double, 3 > &t, double jac[4][4], const double det, double &ex, double &ey)
 Interpolate the field in a triangle.
 
static double Potential5 (const std::array< double, 8 > &v, const std::array< double, 2 > &t)
 Interpolate the potential in a curved quadrilateral.
 
static void Field5 (const std::array< double, 8 > &v, const std::array< double, 2 > &t, double jac[4][4], const double det, double &ex, double &ey)
 Interpolate the field in a curved quadrilateral.
 
static double Potential13 (const std::array< double, 10 > &v, const std::array< double, 4 > &t)
 Interpolate the potential in a curved quadratic tetrahedron.
 
static void Field13 (const std::array< double, 10 > &v, const std::array< double, 4 > &t, double jac[4][4], const double det, double &ex, double &ey, double &ez)
 Interpolate the field in a curved quadratic tetrahedron.
 
static int ReadInteger (char *token, int def, bool &error)
 
static double ReadDouble (char *token, double def, bool &error)
 
- Protected Attributes inherited from Garfield::ComponentFieldMap
bool m_is3d = true
 
ElementType m_elementType = ElementType::CurvedTetrahedron
 
std::vector< Element > m_elements
 
std::vector< int > m_elementIndices
 
std::vector< bool > m_degenerate
 
std::vector< std::array< double, 3 > > m_bbMin
 
std::vector< std::array< double, 3 > > m_bbMax
 
std::vector< std::array< std::array< double, 3 >, 4 > > m_w12
 
std::vector< Node > m_nodes
 
std::vector< double > m_pot
 
std::map< std::string, std::vector< double > > m_wpot
 
std::map< std::string, std::vector< std::vector< double > > > m_dwpot
 
std::vector< Material > m_materials
 
std::map< std::string, WeightingFieldCopy > m_wfieldCopies
 
std::vector< double > m_wdtimes
 
bool m_hasBoundingBox = false
 
std::array< double, 3 > m_minBoundingBox = {{0., 0., 0.}}
 
std::array< double, 3 > m_maxBoundingBox = {{0., 0., 0.}}
 
std::array< double, 3 > m_mapmin = {{0., 0., 0.}}
 
std::array< double, 3 > m_mapmax = {{0., 0., 0.}}
 
std::array< double, 3 > m_mapamin = {{0., 0., 0.}}
 
std::array< double, 3 > m_mapamax = {{0., 0., 0.}}
 
std::array< double, 3 > m_mapna = {{0., 0., 0.}}
 
std::array< double, 3 > m_cells = {{0., 0., 0.}}
 
double m_mapvmin = 0.
 
double m_mapvmax = 0.
 
std::array< bool, 3 > m_setang
 
bool m_deleteBackground = true
 
bool m_warning = false
 
unsigned int m_nWarnings = 0
 
bool m_printConvergenceWarnings = true
 
- 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

Component for importing and interpolating field maps from CST. This interface assumes a certain format of the ascii files Please find the tools to extract the field data from CST in the correct way here: http://www.desy.de/~zenker/garfieldpp.html

Definition at line 15 of file ComponentCST.hh.

Constructor & Destructor Documentation

◆ ComponentCST()

Garfield::ComponentCST::ComponentCST ( )

Constructor.

Definition at line 80 of file ComponentCST.cc.

80 : ComponentFieldMap("CST") {
81 // Default bounding box
82 m_minBoundingBox[2] = -50.;
83 m_maxBoundingBox[2] = 50.;
84 m_deleteBackground = false;
85}
Garfield::ComponentFieldMap::m_maxBoundingBox
std::array< double, 3 > m_maxBoundingBox
Definition ComponentFieldMap.hh:206
Garfield::ComponentFieldMap::m_minBoundingBox
std::array< double, 3 > m_minBoundingBox
Definition ComponentFieldMap.hh:205
Garfield::ComponentFieldMap::ComponentFieldMap
ComponentFieldMap()=delete
Default constructor.

◆ ~ComponentCST()

Garfield::ComponentCST::~ComponentCST ( )
inline

Destructor.

Definition at line 20 of file ComponentCST.hh.

20{}

Member Function Documentation

◆ Coordinate2Index() [1/2]

bool Garfield::ComponentCST::Coordinate2Index ( const double x,
const double y,
const double z,
unsigned int & i,
unsigned int & j,
unsigned int & k ) const

Find the positions in the x/y/z position vectors (m_xlines, m_ylines, m_zlines) for a given point. The operator used for the comparison is <=. Therefore, the last entry in the vector will never be returned for a point inside the mesh.

Definition at line 1055 of file ComponentCST.cc.

1057 {
1058 bool mirrored[3] = {false, false, false};
1059 double pos[3] = {0., 0., 0.};
1060 return Coordinate2Index(x, y, z, i, j, k, pos, mirrored);
1061}
Garfield::ComponentCST::Coordinate2Index
bool Coordinate2Index(const double x, const double y, const double z, unsigned int &i, unsigned int &j, unsigned int &k) const
Definition ComponentCST.cc:1055

Referenced by Coordinate2Index(), GetMedium(), WeightingField(), and WeightingPotential().

◆ Coordinate2Index() [2/2]

bool Garfield::ComponentCST::Coordinate2Index ( const double x,
const double y,
const double z,
unsigned int & i,
unsigned int & j,
unsigned int & k,
double * position_mapped,
bool * mirrored ) const
protected

Calculate the index in the vectors m_xlines, m_ylines, m_zlines, which is before the given coordinate.

Remarks
x, y, z need to be mapped before using this function!
Parameters
xThe x coordinate mapped to the basic cell.
yThe y coordinate mapped to the basic cell.
zThe z coordinate mapped to the basic cell.
iThe index of the m_xlines vector, where m_xlines.at(i) < x < m_xlines.at(i+1).
jThe index of the m_ylines vector, where m_ylines.at(j) < y < m_ylines.at(j+1).
kThe index of the m_zlines vector, where m_zlines.at(k) < z < m_zlines.at(k+1).
position_mappedThe calculated mapped position (x,y,z) -> (x_mapped, y_mapped, z_mapped)
mirroredInformation if x, y, or z direction is mirrored.

Definition at line 1064 of file ComponentCST.cc.

1067 {
1068 // Map the coordinates onto field map coordinates
1069 pos[0] = xin;
1070 pos[1] = yin;
1071 pos[2] = zin;
1072 double rcoordinate = 0.;
1073 double rotation = 0.;
1074 MapCoordinates(pos[0], pos[1], pos[2],
1075 mirrored[0], mirrored[1], mirrored[2], rcoordinate, rotation);
1076
1077 auto it_x =
1078 std::lower_bound(m_xlines.begin(), m_xlines.end(), pos[0]);
1079 auto it_y =
1080 std::lower_bound(m_ylines.begin(), m_ylines.end(), pos[1]);
1081 auto it_z =
1082 std::lower_bound(m_zlines.begin(), m_zlines.end(), pos[2]);
1083 if (it_x == m_xlines.end() || it_y == m_ylines.end() ||
1084 it_z == m_zlines.end() || pos[0] < m_xlines.at(0) ||
1085 pos[1] < m_ylines.at(0) ||
1086 pos[2] < m_zlines.at(0)) {
1087 if (m_debug) {
1088 std::cerr << m_className << "::ElectricFieldBinary:" << std::endl;
1089 std::cerr << " Could not find the given coordinate!" << std::endl;
1090 std::cerr << " You ask for the following position: " << xin << ", "
1091 << yin << ", " << zin << std::endl;
1092 std::cerr << " The mapped position is: " << pos[0] << ", "
1093 << pos[1] << ", " << pos[2]
1094 << std::endl;
1095 }
1096 return false;
1097 }
1098 /* Lower bound returns the next mesh line behind the position in question.
1099 * If the position in question is on a mesh line this mesh line is returned.
1100 * Since we are interested in the mesh line before the position in question we
1101 * need to move the
1102 * iterator to the left except for the very first mesh line!
1103 */
1104 if (it_x == m_xlines.begin())
1105 i = 0;
1106 else
1107 i = std::distance(m_xlines.begin(), it_x - 1);
1108 if (it_y == m_ylines.begin())
1109 j = 0;
1110 else
1111 j = std::distance(m_ylines.begin(), it_y - 1);
1112 if (it_z == m_zlines.begin())
1113 k = 0;
1114 else
1115 k = std::distance(m_zlines.begin(), it_z - 1);
1116 return true;
1117}
Garfield::ComponentFieldMap::MapCoordinates
void MapCoordinates(double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
Move (xpos, ypos, zpos) to field map coordinates.
Definition ComponentFieldMap.cc:2446
Garfield::Component::m_debug
bool m_debug
Switch on/off debugging messages.
Definition Component.hh:371
Garfield::Component::m_className
std::string m_className
Class name.
Definition Component.hh:359

◆ DisableShaping()

void Garfield::ComponentCST::DisableShaping ( )
inline

Definition at line 181 of file ComponentCST.hh.

181{ doShaping = false; }

◆ DisableXField()

void Garfield::ComponentCST::DisableXField ( )
inline

Use these functions to disable a certain field component. Is a field component is disabled ElectricField and WeightingField will return 0 for this component. This is useful if you want to have calculated global field distortions and you want to add the field of a GEM. If you would simply add both components the field component in drift direction would be added twice!

Definition at line 134 of file ComponentCST.hh.

134{ disableFieldComponent[0] = true; }

◆ DisableYField()

void Garfield::ComponentCST::DisableYField ( )
inline

Definition at line 135 of file ComponentCST.hh.

135{ disableFieldComponent[1] = true; }

◆ DisableZField()

void Garfield::ComponentCST::DisableZField ( )
inline

Definition at line 136 of file ComponentCST.hh.

136{ disableFieldComponent[2] = true; }

◆ ElectricField() [1/2]

void Garfield::ComponentCST::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 808 of file ComponentCST.cc.

811 {
812 ElectricFieldBinary(xin, yin, zin, ex, ey, ez, volt, m, status, true);
813}

◆ ElectricField() [2/2]

void Garfield::ComponentCST::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 801 of file ComponentCST.cc.

803 {
804 double volt;
805 ElectricFieldBinary(xin, yin, zin, ex, ey, ez, volt, m, status);
806}

◆ EnableShaping()

void Garfield::ComponentCST::EnableShaping ( )
inline

If you calculate the electric field component in $x$ direction along a line in x direction this field component will be constant inside mesh elements (by construction). This can be observed by plotting $E_x$ in $x$ direction. If you plot $E_x$ in y direction the field will be smooth (also by construction). Yuri Piadyk proposed also to shape the electric field. This is done as follows. The field component calculated as described above is assumed to appear in the center of the mesh element.

| M1 P | M2 | x direction |-—x—x-|--—x---—|-->

__________ ____________

element 1 element 2

Lets consider only the $x$ direction and we want to calculate $E_x(P)$. The field in the center of the element containing $P$ is $E_x(M_1) = E_1$. Without shaping it is $E_1$ along the $x$ direction in everywhere in element 1. The idea of the shaping is to do a linear interpolation of the $E_x$ between the field $E_1$ and $E_x(M_2)=E_2$. This results in a smooth electric field $E_x$ also in $x$ direction. If P would be left from $M_1$ the field in the left neighboring element would be considered. In addition it is also checked if the material in both elements used for the interpolation is the same. Else no interpolation is done.

Remarks
This shaping gives you a nice and smooth field, but you introduce additional information. This information is not coming from the CST simulation, but from the assumption that the field between elements changes in a linear way, which might be wrong! So you might consider to increase the number of mesh cells used in the simulation rather than using this smoothing.

Definition at line 180 of file ComponentCST.hh.

180{ doShaping = true; }

◆ GetAspectRatio()

void Garfield::ComponentCST::GetAspectRatio ( const size_t i,
double & dmin,
double & dmax ) const
overrideprotectedvirtual

Reimplemented from Garfield::ComponentFieldMap.

Definition at line 1127 of file ComponentCST.cc.

1128 {
1129 if (element >= m_nElements) {
1130 dmin = dmax = 0.;
1131 return;
1132 }
1133 unsigned int i, j, k;
1134 Element2Index(element, i, j, k);
1135 const double dx = fabs(m_xlines.at(i + 1) - m_xlines.at(i));
1136 const double dy = fabs(m_ylines.at(j + 1) - m_ylines.at(j));
1137 const double dz = fabs(m_zlines.at(k + 1) - m_zlines.at(k));
1138 dmin = std::min({dx, dy, dz});
1139 dmax = std::max({dx, dy, dz});
1140}
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition DoubleAc.h:615

◆ GetElement()

bool Garfield::ComponentCST::GetElement ( const size_t i,
size_t & mat,
bool & drift,
std::vector< size_t > & nodes ) const
overridevirtual

Return the material and node indices of a mesh element.

Reimplemented from Garfield::ComponentFieldMap.

Definition at line 948 of file ComponentCST.cc.

949 {
950 if (element >= m_nElements || element >= m_elementMaterial.size()) {
951 std::cerr << m_className << "::GetElement: Index out of range.\n";
952 return false;
953 }
954 mat = m_elementMaterial[element];
955 drift = m_materials[mat].driftmedium;
956 nodes.clear();
957 unsigned int i0 = 0, j0 = 0, k0 = 0;
958 Element2Index(element, i0, j0, k0);
959 const auto i1 = i0 + 1;
960 const auto j1 = j0 + 1;
961 const auto k1 = k0 + 1;
962 nodes.push_back(Index2Node(i0, j0, k0));
963 nodes.push_back(Index2Node(i1, j0, k0));
964 nodes.push_back(Index2Node(i0, j1, k0));
965 nodes.push_back(Index2Node(i1, j1, k0));
966 nodes.push_back(Index2Node(i0, j0, k1));
967 nodes.push_back(Index2Node(i1, j0, k1));
968 nodes.push_back(Index2Node(i0, j1, k1));
969 nodes.push_back(Index2Node(i1, j1, k1));
970 return true;
971}
Garfield::ComponentFieldMap::m_materials
std::vector< Material > m_materials
Definition ComponentFieldMap.hh:188

◆ GetElementBoundaries()

void Garfield::ComponentCST::GetElementBoundaries ( unsigned int element,
double & xmin,
double & xmax,
double & ymin,
double & ymax,
double & zmin,
double & zmax ) const

Definition at line 987 of file ComponentCST.cc.

990 {
991 unsigned int i, j, k;
992 Element2Index(element, i, j, k);
993 xmin = m_xlines.at(i);
994 xmax = m_xlines.at(i + 1);
995 ymin = m_ylines.at(j);
996 ymax = m_ylines.at(j + 1);
997 zmin = m_zlines.at(k);
998 zmax = m_zlines.at(k + 1);
999}

◆ GetElementVolume()

double Garfield::ComponentCST::GetElementVolume ( const size_t i) const
overrideprotectedvirtual

Reimplemented from Garfield::ComponentFieldMap.

Definition at line 1142 of file ComponentCST.cc.

1142 {
1143 if (element >= m_nElements) return 0.;
1144 unsigned int i, j, k;
1145 Element2Index(element, i, j, k);
1146 const double dx = fabs(m_xlines.at(i + 1) - m_xlines.at(i));
1147 const double dy = fabs(m_ylines.at(j + 1) - m_ylines.at(j));
1148 const double dz = fabs(m_zlines.at(k + 1) - m_zlines.at(k));
1149 return dx * dy * dz;
1150}

◆ GetMedium()

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

1002 {
1003 unsigned int i, j, k;
1004 Coordinate2Index(x, y, z, i, j, k);
1005 if (m_debug) {
1006 std::cout << m_className << "::GetMedium:\n"
1007 << " Position (" << x << ", " << y << ", " << z << "):\n"
1008 << " Indices are: x: " << i << "/" << m_xlines.size()
1009 << "\t y: " << j << "/" << m_ylines.size()
1010 << "\t z: " << k << "/" << m_zlines.size() << std::endl;
1011 const auto element = Index2Element(i, j, k);
1012 std::cout << " Element index: " << element << std::endl
1013 << " Material index: "
1014 << (int)m_elementMaterial.at(element) << std::endl;
1015 }
1016 return m_materials.at(m_elementMaterial.at(Index2Element(i, j, k))).medium;
1017}
Garfield::ComponentCST::Index2Element
int Index2Element(const unsigned int i, const unsigned int j, const unsigned int k) const
Definition ComponentCST.cc:1119

◆ GetNode()

bool Garfield::ComponentCST::GetNode ( const size_t i,
double & x,
double & y,
double & z ) const
overridevirtual

Reimplemented from Garfield::ComponentFieldMap.

Definition at line 973 of file ComponentCST.cc.

974 {
975 if (node >= m_nNodes) {
976 std::cerr << m_className << "::GetNode: Index out of range.\n";
977 return false;
978 }
979 unsigned int i = 0, j = 0, k = 0;
980 Node2Index(node, i, j, k);
981 x = m_xlines[i];
982 y = m_ylines[j];
983 z = m_zlines[k];
984 return true;
985}

◆ GetNumberOfElements()

size_t Garfield::ComponentCST::GetNumberOfElements ( ) const
inlineoverridevirtual

Return the number of mesh elements.

Reimplemented from Garfield::ComponentFieldMap.

Definition at line 27 of file ComponentCST.hh.

27{ return m_nElements; }

◆ GetNumberOfMeshLines()

void Garfield::ComponentCST::GetNumberOfMeshLines ( unsigned int & nx,
unsigned int & ny,
unsigned int & nz ) const

Definition at line 941 of file ComponentCST.cc.

942 {
943 n_x = m_xlines.size();
944 n_y = m_ylines.size();
945 n_z = m_zlines.size();
946}

◆ GetNumberOfNodes()

size_t Garfield::ComponentCST::GetNumberOfNodes ( ) const
inlineoverridevirtual

Reimplemented from Garfield::ComponentFieldMap.

Definition at line 30 of file ComponentCST.hh.

30{ return m_nNodes; }

◆ Index2Element()

int Garfield::ComponentCST::Index2Element ( const unsigned int i,
const unsigned int j,
const unsigned int k ) const

Calculate the element index from the position in the x/y/z position vectors (m_xlines, m_ylines, m_zlines). This is public since it is used in ViewFEMesh::DrawCST.

Definition at line 1119 of file ComponentCST.cc.

1120 {
1121 if (i > m_nx - 2 || j > m_ny - 2 || k > m_nz - 2) {
1122 throw "ComponentCST::Index2Element: Error. Element indices out of bounds.";
1123 }
1124 return i + j * (m_nx - 1) + k * (m_nx - 1) * (m_ny - 1);
1125}

Referenced by GetMedium(), and WeightingField().

◆ Initialise() [1/2]

bool Garfield::ComponentCST::Initialise ( std::string dataFile,
std::string unit = "cm" )

Import of field data based on binary files. See http://www.desy.de/~zenker/garfieldpp.html to get information about the binary files export from CST.

Parameters
dataFileThe binary file containing the field data exported from CST.
unitThe units used in the binary file. They are not necessarily equal to CST units.

Definition at line 472 of file ComponentCST.cc.

472 {
473 m_ready = false;
474
475 // Keep track of the success
476 bool ok = true;
477 // Check the value of the unit
478 double funit = ScalingFactor(unit);
479 if (funit <= 0.) {
480 std::cerr << m_className << "::Initialise:\n"
481 << " Unknown length unit " << unit << ".\n";
482 ok = false;
483 funit = 1.0;
484 }
485 if (m_debug) {
486 std::cout << m_className << "::Initialise: Unit scaling factor = "
487 << funit << ".\n";
488 }
489 FILE* f = fopen(dataFile.c_str(), "rb");
490 if (f == nullptr) {
491 PrintCouldNotOpen("Initialise", dataFile);
492 return false;
493 }
494
495 struct stat fileStatus;
496 stat(dataFile.c_str(), &fileStatus);
497 int fileSize = fileStatus.st_size;
498
499 int nLinesX = 0, nLinesY = 0, nLinesZ = 0;
500 int nNS = 0, nES = 0, nEM = 0;
501 int nMaterials = 0;
502 if (!ReadHeader(f, fileSize, m_debug, nLinesX, nLinesY, nLinesZ,
503 nNS, nES, nEM, nMaterials)) {
504 fclose(f);
505 return false;
506 }
507 m_nx = nLinesX;
508 m_ny = nLinesY;
509 m_nz = nLinesZ;
510 m_nNodes = m_nx * m_ny * m_nz;
511 m_nElements = (m_nx - 1) * (m_ny - 1) * (m_nz - 1);
512
513 m_xlines.resize(nLinesX);
514 m_ylines.resize(nLinesY);
515 m_zlines.resize(nLinesZ);
516 m_potential.resize(nNS);
517 m_elementMaterial.resize(nEM);
518 m_materials.resize(nMaterials);
519 auto result = fread(m_xlines.data(), sizeof(double), m_xlines.size(), f);
520 if (result != m_xlines.size()) {
521 fputs("Reading error while reading xlines.", stderr);
522 exit(3);
523 } else if (result == 0) {
524 fputs("No xlines are stored in the data file.", stderr);
525 exit(3);
526 }
527 result = fread(m_ylines.data(), sizeof(double), m_ylines.size(), f);
528 if (result != m_ylines.size()) {
529 fputs("Reading error while reading ylines", stderr);
530 exit(3);
531 } else if (result == 0) {
532 fputs("No ylines are stored in the data file.", stderr);
533 exit(3);
534 }
535 result = fread(m_zlines.data(), sizeof(double), m_zlines.size(), f);
536 if (result != m_zlines.size()) {
537 fputs("Reading error while reading zlines", stderr);
538 exit(3);
539 } else if (result == 0) {
540 fputs("No zlines are stored in the data file.", stderr);
541 exit(3);
542 }
543 result = fread(m_potential.data(), sizeof(float), m_potential.size(), f);
544 if (result != m_potential.size()) {
545 fputs("Reading error while reading nodes.", stderr);
546 exit(3);
547 } else if (result == 0) {
548 fputs("No potentials are stored in the data file.", stderr);
549 exit(3);
550 }
551 fseek(f, nES * sizeof(float), SEEK_CUR);
552 // not needed in principle - thus it is ok if nothing is read
553 result = fread(m_elementMaterial.data(), sizeof(unsigned char),
554 m_elementMaterial.size(), f);
555 if (result != m_elementMaterial.size()) {
556 fputs("Reading error while reading element material", stderr);
557 exit(3);
558 }
559 std::stringstream st;
560 st << m_className << "::Initialise:" << std::endl;
561 /*
562 * The material vector is filled according to the material id!
563 * Thus material.at(0) is material with id 0.
564 */
565 for (unsigned int i = 0; i < m_materials.size(); i++) {
566 float id;
567 result = fread(&(id), sizeof(float), 1, f);
568 if (result != 1) {
569 fputs("Input error while reading material id.", stderr);
570 exit(3);
571 }
572 // const unsigned int index = id;
573 const unsigned int index = i;
574 unsigned int description_size = 0;
575 result = fread(&(description_size), sizeof(int), 1, f);
576 if (result != 1) {
577 fputs("Input error while reading material description size.", stderr);
578 exit(3);
579 }
580 char* c = new char[description_size];
581 result = fread(c, sizeof(char), description_size, f);
582 if (result != description_size) {
583 fputs("Input error while reading material description.", stderr);
584 exit(3);
585 }
586 std::string name = c;
587 st << " Read material: " << name;
588 if (name.compare("gas") == 0) {
589 st << " (considered as drift medium)";
590 m_materials.at(index).driftmedium = true;
591 } else {
592 m_materials.at(index).driftmedium = false;
593 }
594 delete[] c;
595 float eps;
596 result = fread(&(eps), sizeof(float), 1, f);
597 m_materials.at(index).eps = eps;
598 if (result != 1) {
599 fputs("Reading error while reading eps.", stderr);
600 exit(3);
601 }
602 st << "; eps is: " << m_materials.at(index).eps;
603 // float mue;
604 // result = fread(&(mue), sizeof(float), 1, f);
605 // if (result != 1) {
606 // fputs ("Reading error while reading mue.", stderr);
607 // exit (3);
608 // }
609 // st << "\t mue is: " << mue;
610 // float rho;
611 // result = fread(&(rho), sizeof(float), 1, f);
612 // if (result != 1) {
613 // fputs ("Reading error while reading rho.", stderr);
614 // exit (3);
615 // }
616 // st << "\t rho is: " << rho;
617 st << "\t id is: " << id << std::endl;
618 // Skip mue and rho
619 fseek(f, 2 * sizeof(float), SEEK_CUR);
620 // ToDo: Check if rho should be used to decide, which material is driftable
621 }
622 // To be sure that they are sorted (should be already be the case)
623 std::sort(m_xlines.begin(), m_xlines.end());
624 std::sort(m_ylines.begin(), m_ylines.end());
625 std::sort(m_zlines.begin(), m_zlines.end());
626 if (funit != 1) {
627 std::transform(m_xlines.begin(), m_xlines.end(), m_xlines.begin(), [funit](double x) { return x * funit;});
628 std::transform(m_ylines.begin(), m_ylines.end(), m_ylines.begin(), [funit](double x) { return x * funit;});
629 std::transform(m_zlines.begin(), m_zlines.end(), m_zlines.begin(), [funit](double x) { return x * funit;});
630 }
631
632 std::cout << m_className << "::Initialise" << std::endl;
633 std::cout << " x range: " << *(m_xlines.begin()) << " - "
634 << *(m_xlines.end() - 1) << std::endl;
635 std::cout << " y range: " << *(m_ylines.begin()) << " - "
636 << *(m_ylines.end() - 1) << std::endl;
637 std::cout << " z range: " << *(m_zlines.begin()) << " - "
638 << *(m_zlines.end() - 1) << std::endl;
639 fclose(f);
640 // Set the ready flag
641 if (ok) {
642 m_ready = true;
643 } else {
644 std::cerr << m_className << "::Initialise:" << std::endl;
645 std::cerr << " Field map could not be read and cannot be interpolated."
646 << std::endl;
647 return false;
648 }
649
650 SetRange();
651 UpdatePeriodicity();
652 return true;
653}
Garfield::ComponentCST::SetRange
void SetRange() override
Definition ComponentCST.cc:1019
Garfield::ComponentFieldMap::ScalingFactor
static double ScalingFactor(std::string unit)
Definition ComponentFieldMap.cc:2637
Garfield::ComponentFieldMap::UpdatePeriodicity
void UpdatePeriodicity() override
Verify periodicities.
Definition ComponentFieldMap.hh:244
Garfield::ComponentFieldMap::PrintCouldNotOpen
void PrintCouldNotOpen(const std::string &header, const std::string &filename) const
Definition ComponentFieldMap.cc:2785
Garfield::Component::m_ready
bool m_ready
Ready for use?
Definition Component.hh:368

◆ Initialise() [2/2]

bool Garfield::ComponentCST::Initialise ( std::string elist,
std::string nlist,
std::string mplist,
std::string prnsol,
std::string unit = "cm" )

Deprecated version of the interface based on text file import of field data.

Parameters
elistInformation about the element material of mesh cells. Each line contains the element number and the material index:
0 3
...
nlistInformation about the mesh like this:
xmax 136 ymax 79 zmax 425
x−l i n e s
0
8 . 9 2 8 5 7 e −07
1 . 7 8 5 7 1 e −06
...
y−l i n e s
0
8 . 9 2 8 5 7 e −07
1 . 7 8 5 7 1 e −06
...
z−l i n e s
0.0027
0.00270674
...
mplistInformation about material properties used in the simulation:
Materials 4
Material 1 PERX 1 . 0 0 0 0 0 0
Material 2 RSVX 0 . 0 0 0 0 0 0 PERX 0 . 1 0 0 0 0 0 0E+11
Material 3 PERX 3 . 5 0 0 0 0 0
Material 4 PERX 4 . 8 0 0 0 0 0
prnsolInformation about the node potentials. Each line contains the node number and the potential:
0 1000.00
...
unitThe units used in the nlist input file

Definition at line 87 of file ComponentCST.cc.

89 {
90 Reset();
91 // Keep track of the success
92 bool ok = true;
93
94 // Buffer for reading
95 const int size = 200;
96 char line[size];
97 // Open the material list
98 std::ifstream fmplist(mplist);
99 if (!fmplist) {
100 PrintCouldNotOpen("Initialise", mplist);
101 return false;
102 }
103
104 // Read the material list
105 int il = 0;
106 bool readerror = false;
107 while (fmplist.getline(line, size, '\n')) {
108 il++;
109 // Split the line in tokens
110 char* token = strtok(line, " ");
111 // Skip blank lines and headers
112 if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
113 int(token[0]) == 10 || int(token[0]) == 13)
114 continue;
115 // Read number of materials,
116 // ensure it does not exceed the maximum and initialize the list
117 if (strcmp(token, "Materials") == 0) {
118 token = strtok(nullptr, " ");
119 const int nMaterials = ReadInteger(token, -1, readerror);
120 if (readerror) {
121 std::cerr << m_className << "::Initialise:\n"
122 << " Error reading file " << mplist << " (line " << il
123 << ")." << std::endl;
124 fmplist.close();
125 return false;
126 }
127 m_materials.resize(nMaterials);
128 for (auto& material : m_materials) {
129 material.ohm = -1;
130 material.eps = -1;
131 material.medium = nullptr;
132 }
133 if (m_debug) {
134 std::cout << m_className << "::Initialise:" << std::endl;
135 std::cout << " Number of materials: " << nMaterials << "\n";
136 }
137 } else if (strcmp(token, "Material") == 0) {
138 token = strtok(nullptr, " ");
139 int imat = ReadInteger(token, -1, readerror);
140 if (readerror) {
141 std::cerr << m_className << "::Initialise:" << std::endl;
142 std::cerr << " Error reading file " << mplist << " (line " << il
143 << ").\n";
144 fmplist.close();
145 return false;
146 } else if (imat < 1 || imat > (int)m_materials.size()) {
147 std::cerr << m_className << "::Initialise:\n"
148 << " Found out-of-range material index " << imat << "in\n"
149 << " material properties file " << mplist << ".\n";
150 ok = false;
151 } else {
152 token = strtok(nullptr, " ");
153 int itype = 0;
154 if (strcmp(token, "PERX") == 0) {
155 itype = 1;
156 } else if (strcmp(token, "RSVX") == 0) {
157 itype = 2;
158 } else {
159 std::cerr << m_className << "::Initialise:\n"
160 << " Unknown material property flag " << token << "\n"
161 << " in material properties file " << mplist
162 << " (line " << il << ").\n";
163 ok = false;
164 }
165 token = strtok(nullptr, " ");
166 if (itype == 1) {
167 m_materials[imat - 1].eps = ReadDouble(token, -1, readerror);
168 } else if (itype == 2) {
169 m_materials[imat - 1].ohm = ReadDouble(token, -1, readerror);
170 token = strtok(nullptr, " ");
171 if (strcmp(token, "PERX") != 0) {
172 std::cerr << m_className << "::Initialise:\n"
173 << " Unknown material property flag " << token << "\n"
174 << " in material file " << mplist << " (material "
175 << imat << ").\n";
176 ok = false;
177 } else {
178 token = strtok(nullptr, " ");
179 m_materials[imat - 1].eps = ReadDouble(token, -1, readerror);
180 }
181 }
182 if (readerror) {
183 std::cerr << m_className << "::Initialise:\n"
184 << " Error reading file " << mplist
185 << " (line " << il << ")." << std::endl;
186 fmplist.close();
187 return false;
188 }
189 if (m_debug) {
190 std::cout << m_className << "::Initialise:" << std::endl;
191 std::cout << " Read material properties for material "
192 << (imat - 1) << "" << std::endl;
193 if (itype == 2) {
194 std::cout << " eps = " << m_materials[imat - 1].eps
195 << " ohm = " << m_materials[imat - 1].ohm << ""
196 << std::endl;
197 } else {
198 std::cout << " eps = " << m_materials[imat - 1].eps << ""
199 << std::endl;
200 }
201 }
202 }
203 }
204 }
205 // Close the file
206 fmplist.close();
207
208 // Find lowest epsilon, check for eps = 0, set default drift medium.
209 if (!SetDefaultDriftMedium()) ok = false;
210
211 // Tell how many lines read
212 std::cout << m_className << "::Initialise:\n"
213 << " Read properties of " << m_materials.size() << " materials\n"
214 << " from file " << mplist << "." << std::endl;
215 if (m_debug) PrintMaterials();
216
217 // Check the value of the unit
218 double funit = ScalingFactor(unit);
219 if (funit <= 0.) {
220 std::cerr << m_className << "::Initialise:\n"
221 << " Unknown length unit " << unit << ".\n";
222 ok = false;
223 funit = 1.0;
224 }
225 if (m_debug) {
226 std::cout << m_className << "::Initialise: Unit scaling factor = "
227 << funit << ".\n";
228 }
229
230 // Open the node list
231 std::ifstream fnlist(nlist);
232 if (!fnlist) {
233 PrintCouldNotOpen("Initialise", nlist);
234 return false;
235 }
236 // Read the node list
237 m_nNodes = 0;
238 il = 0;
239 int xlines = 0, ylines = 0, zlines = 0;
240 int lines_type = -1;
241 double line_tmp;
242 while (fnlist.getline(line, size, '\n')) {
243 il++;
244 // Split the line in tokens
245 char* token = strtok(line, " ");
246 // Skip blank lines and headers
247 if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
248 int(token[0]) == 10 || int(token[0]) == 13)
249 continue;
250 // Read max sizes
251 if (strcmp(token, "xmax") == 0) {
252 token = strtok(nullptr, " ");
253 xlines = ReadInteger(token, -1, readerror);
254 token = strtok(nullptr, " ");
255 token = strtok(nullptr, " ");
256 ylines = ReadInteger(token, -1, readerror);
257 token = strtok(nullptr, " ");
258 token = strtok(nullptr, " ");
259 zlines = ReadInteger(token, -1, readerror);
260 if (readerror) break;
261 continue;
262 }
263 if (strcmp(token, "x-lines\n") == 0 || strcmp(token, "x-lines") == 0) {
264 lines_type = 1;
265 if (m_debug) {
266 std::cout << m_className << "::Initialise:\n"
267 << " Reading x-lines from file " << nlist << ".\n";
268 }
269 continue;
270 }
271 if (strcmp(token, "y-lines\n") == 0 || strcmp(token, "y-lines") == 0) {
272 lines_type = 2;
273 if (m_debug) {
274 std::cout << m_className << "::Initialise:\n"
275 << " Reading y-lines from file " << nlist << ".\n";
276 }
277 continue;
278 }
279 if (strcmp(token, "z-lines\n") == 0 || strcmp(token, "z-lines") == 0) {
280 lines_type = 3;
281 if (m_debug) {
282 std::cout << m_className << "::Initialise:\n"
283 << " Reading z-lines from file " << nlist << ".\n";
284 }
285 continue;
286 }
287 line_tmp = ReadDouble(token, -1, readerror);
288 if (lines_type == 1)
289 m_xlines.push_back(line_tmp * funit);
290 else if (lines_type == 2)
291 m_ylines.push_back(line_tmp * funit);
292 else if (lines_type == 3)
293 m_zlines.push_back(line_tmp * funit);
294 else {
295 std::cerr << m_className << "::Initialise:" << std::endl;
296 std::cerr << " Line type was not set in " << nlist << " (line " << il
297 << ", token = " << token << ")." << std::endl;
298 std::cerr << " Maybe it is in the wrong format" << std::endl;
299 std::cerr << " e.g. missing tailing space after x-lines." << std::endl;
300 ok = false;
301 break;
302 }
303 if (readerror) break;
304 }
305 // Check syntax
306 if (readerror) {
307 std::cerr << m_className << "::Initialise:\n"
308 << " Error reading file " << nlist
309 << " (line " << il << ").\n";
310 fnlist.close();
311 return false;
312 }
313 // Close the file
314 fnlist.close();
315
316 if ((unsigned)xlines == m_xlines.size() &&
317 (unsigned)ylines == m_ylines.size() &&
318 (unsigned)zlines == m_zlines.size()) {
319 std::cout << m_className << "::Initialise:" << std::endl;
320 std::cout << " Found in file " << nlist << "\n " << xlines
321 << " x-lines\n " << ylines << " y-lines\n " << zlines
322 << " z-lines" << std::endl;
323 } else {
324 std::cerr << m_className << "::Initialise:" << std::endl;
325 std::cerr << " There should be " << xlines << " x-lines, " << ylines
326 << " y-lines and " << zlines << " z-lines in file " << nlist
327 << " but I found :\n " << m_xlines.size() << " x-lines, "
328 << m_ylines.size() << " x-lines, " << m_zlines.size()
329 << " z-lines." << std::endl;
330 }
331 m_nx = m_xlines.size();
332 m_ny = m_ylines.size();
333 m_nz = m_zlines.size();
334 m_nNodes = m_nx * m_ny * m_nz;
335 m_nElements = (m_nx - 1) * (m_ny - 1) * (m_nz - 1);
336
337 // Tell how many lines read
338 std::cout << m_className << "::Initialise:" << std::endl;
339 std::cout << " Read " << m_nNodes << " nodes from file " << nlist << "."
340 << std::endl;
341
342 // Open the element list
343 std::ifstream felist(elist);
344 if (!felist) {
345 PrintCouldNotOpen("Initialise", elist);
346 return false;
347 }
348 // Read the element list
349 m_elementMaterial.resize(m_nElements);
350 il = 0;
351 while (felist.getline(line, size, '\n')) {
352 il++;
353 // Split the line in tokens
354 char* token = strtok(line, " ");
355 // Skip blank lines and headers
356 if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
357 int(token[0]) == 10 || int(token[0]) == 13 ||
358 strcmp(token, "LIST") == 0 || strcmp(token, "ELEM") == 0)
359 continue;
360 // Read the element
361 int ielem = ReadInteger(token, -1, readerror);
362 token = strtok(nullptr, " ");
363 if (!token) continue;
364 unsigned char imat = atoi(token);
365 // construct node numbers
366 std::vector<int> node_nb;
367 try {
368 // Read element material - the number of the material is stored (1, 2,
369 // ...) but we need the index (0, 1, ...)
370 m_elementMaterial.at(ielem) = (imat - 1);
371 } catch (...) {
372 std::cerr << m_className << "::Initialise:" << std::endl;
373 std::cerr << " Error reading file " << elist << " (line " << il << ")."
374 << std::endl;
375 std::cerr << " The element index (" << ielem
376 << ") is not in the expected range: 0 - " << m_nElements
377 << std::endl;
378 ok = false;
379 }
380 // Check the material number and ensure that epsilon is non-negative
381 // int check_mat = imat;
382 if (imat < 1 || imat > m_materials.size()) {
383 std::cerr << m_className << "::Initialise:" << std::endl;
384 std::cerr << " Out-of-range material number on file " << elist
385 << " (line " << il << ")." << std::endl;
386 std::cerr << " Element: " << ielem << ", material: " << imat
387 << std::endl;
388 ok = false;
389 }
390 if (m_materials[imat - 1].eps < 0) {
391 std::cerr << m_className << "::Initialise:" << std::endl;
392 std::cerr << " Element " << ielem << " in element list " << elist
393 << " uses material " << imat << " which" << std::endl;
394 std::cerr << " has not been assigned a positive permittivity"
395 << std::endl;
396 std::cerr << " in material list " << mplist << "." << std::endl;
397 ok = false;
398 }
399 }
400 // Close the file
401 felist.close();
402 // Tell how many lines read
403 std::cout << m_className << "::Initialise:" << std::endl;
404 std::cout << " Read " << m_nElements << " elements.\n";
405
406 // Open the voltage list
407 m_potential.resize(m_nNodes);
408 std::ifstream fprnsol(prnsol);
409 if (!fprnsol) {
410 PrintCouldNotOpen("Initialise", prnsol);
411 return false;
412 }
413 // Read the voltage list
414 il = 0;
415 int nread = 0;
416 readerror = false;
417 while (fprnsol.getline(line, size, '\n')) {
418 il++;
419 // Split the line in tokens
420 char* token = strtok(line, " ");
421 // Skip blank lines and headers
422 if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
423 int(token[0]) == 10 || int(token[0]) == 13 || strcmp(token, "Max") == 0)
424 continue;
425 // Read node potential (in prnsol node id starts with 1 and here we will use
426 // 0 as first node...)
427 int inode = ReadInteger(token, -1, readerror);
428 token = strtok(nullptr, " ");
429 double volt = ReadDouble(token, -1, readerror);
430
431 try {
432 m_potential.at(inode - 1) = volt;
433 nread++;
434 } catch (...) {
435 std::cerr << m_className << "::Initialise:" << std::endl;
436 std::cerr << " Error reading file " << prnsol << " (line " << il
437 << ")." << std::endl;
438 std::cerr << " The node index (" << inode - 1
439 << ") is not in the expected range: 0 - " << m_nNodes
440 << std::endl;
441 ok = false;
442 }
443 }
444 // Close the file
445 fprnsol.close();
446 // Tell how many lines read
447 std::cout << m_className << "::Initialise:" << std::endl;
448 std::cout << " Read " << nread << "/" << m_nNodes
449 << " (expected) potentials from file " << prnsol << "."
450 << std::endl;
451 // Check number of nodes
452 if (nread != (int)m_nNodes) {
453 std::cerr << m_className << "::Initialise:" << std::endl;
454 std::cerr << " Number of nodes read (" << nread << ") on potential file "
455 << prnsol << " does not" << std::endl;
456 std::cerr << " match the node list (" << m_nNodes << ").\n";
457 ok = false;
458 }
459 // Set the ready flag
460 if (ok) {
461 m_ready = true;
462 } else {
463 std::cerr << m_className << "::Initialise:" << std::endl;
464 std::cerr << " Field map could not be read and cannot be interpolated."
465 << std::endl;
466 return false;
467 }
468 Prepare();
469 return true;
470}
Garfield::ComponentFieldMap::PrintMaterials
void PrintMaterials()
List all currently defined materials.
Definition ComponentFieldMap.cc:511
Garfield::ComponentFieldMap::SetDefaultDriftMedium
bool SetDefaultDriftMedium()
Find lowest epsilon, check for eps = 0, set default drift media flags.
Definition ComponentFieldMap.cc:755
Garfield::ComponentFieldMap::ReadDouble
static double ReadDouble(char *token, double def, bool &error)
Definition ComponentFieldMap.cc:2660
Garfield::ComponentFieldMap::ReadInteger
static int ReadInteger(char *token, int def, bool &error)
Definition ComponentFieldMap.cc:2651
Garfield::ComponentFieldMap::Reset
void Reset() override
Reset the component.
Definition ComponentFieldMap.cc:2099

◆ SetRange()

void Garfield::ComponentCST::SetRange ( )
overrideprotectedvirtual

Reimplemented from Garfield::ComponentFieldMap.

Definition at line 1019 of file ComponentCST.cc.

1019 {
1020 // Establish the ranges
1021 m_mapmin[0] = m_xlines.front();
1022 m_mapmax[0] = m_xlines.back();
1023 m_mapmin[1] = m_ylines.front();
1024 m_mapmax[1] = m_ylines.back();
1025 m_mapmin[2] = m_zlines.front();
1026 m_mapmax[2] = m_zlines.back();
1027 m_mapvmin = *std::min_element(m_potential.begin(), m_potential.end());
1028 m_mapvmax = *std::max_element(m_potential.begin(), m_potential.end());
1029
1030 // Set provisional cell dimensions.
1031 m_minBoundingBox[0] = m_mapmin[0];
1032 m_maxBoundingBox[0] = m_mapmax[0];
1033 m_minBoundingBox[1] = m_mapmin[1];
1034 m_maxBoundingBox[1] = m_mapmax[1];
1035 if (m_is3d) {
1036 m_minBoundingBox[2] = m_mapmin[2];
1037 m_maxBoundingBox[2] = m_mapmax[2];
1038 } else {
1039 m_mapmin[2] = m_minBoundingBox[2];
1040 m_mapmax[2] = m_maxBoundingBox[2];
1041 }
1042 m_hasBoundingBox = true;
1043}
Garfield::ComponentFieldMap::m_mapmin
std::array< double, 3 > m_mapmin
Definition ComponentFieldMap.hh:209
Garfield::ComponentFieldMap::m_mapmax
std::array< double, 3 > m_mapmax
Definition ComponentFieldMap.hh:210

Referenced by Initialise(), and ShiftComponent().

◆ SetRangeZ()

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

Definition at line 1045 of file ComponentCST.cc.

1045 {
1046 if (fabs(zmax - zmin) <= 0.) {
1047 std::cerr << m_className << "::SetRangeZ:" << std::endl;
1048 std::cerr << " Zero range is not permitted." << std::endl;
1049 return;
1050 }
1051 m_minBoundingBox[2] = std::min(zmin, zmax);
1052 m_maxBoundingBox[2] = std::max(zmin, zmax);
1053}

◆ SetWeightingField()

bool Garfield::ComponentCST::SetWeightingField ( std::string prnsol,
std::string label,
bool isBinary = true )

Initialise a weighting field. This function can handle the deprecated text based file format (see Initialise( std::string elist...) for the expected file format, which is similar to prnsol. It also also handles binary files including the weighting field.

Parameters
prnsolThe input file (binary/text file)
labelThe name of the weighting field to be added. If a weighting field with same name already exist it is replaced by the new one.
isBinaryDepending on the file type you use, adapt this switch.

Definition at line 655 of file ComponentCST.cc.

656 {
657 std::vector<float> potentials(m_nNodes);
658 if (!m_ready) {
659 std::cerr << m_className << "::SetWeightingField:" << std::endl;
660 std::cerr << " No valid field map is present." << std::endl;
661 std::cerr << " Weighting field cannot be added." << std::endl;
662 return false;
663 }
664
665 // Open the voltage list
666 std::ifstream fprnsol(prnsol);
667 if (!fprnsol) {
668 PrintCouldNotOpen("SetWeightingField", prnsol);
669 return false;
670 }
671 // Check if a weighting field with the same label already exists.
672 auto it = m_weightingFields.find(label);
673 if (it != m_weightingFields.end()) {
674 std::cout << m_className << "::SetWeightingField:" << std::endl;
675 std::cout << " Replacing existing weighting field " << label << "."
676 << std::endl;
677 }
678
679 int nread = 0;
680 bool ok = true;
681
682 if (isBinary) {
683 std::cout << m_className << "::SetWeightingField:" << std::endl;
684 std::cout << " Reading weighting field from binary file:"
685 << prnsol << std::endl;
686 FILE* f = fopen(prnsol.c_str(), "rb");
687 if (f == nullptr) {
688 PrintCouldNotOpen("SetWeightingField", prnsol);
689 return false;
690 }
691
692 struct stat fileStatus;
693 stat(prnsol.c_str(), &fileStatus);
694 int fileSize = fileStatus.st_size;
695
696 int nLinesX = 0, nLinesY = 0, nLinesZ = 0;
697 int nES = 0, nEM = 0;
698 int nMaterials = 0;
699 if (!ReadHeader(f, fileSize, m_debug, nLinesX, nLinesY, nLinesZ,
700 nread, nES, nEM, nMaterials)) {
701 fclose(f);
702 return false;
703 }
704 // Skip everything, but the potential
705 fseek(f, nLinesX * sizeof(double), SEEK_CUR);
706 fseek(f, nLinesY * sizeof(double), SEEK_CUR);
707 fseek(f, nLinesZ * sizeof(double), SEEK_CUR);
708 auto result = fread(potentials.data(), sizeof(float), potentials.size(), f);
709 if (result != potentials.size()) {
710 fputs("Reading error while reading nodes.", stderr);
711 exit(3);
712 } else if (result == 0) {
713 fputs("No weighting potentials are stored in the data file.", stderr);
714 exit(3);
715 }
716 fprnsol.close();
717 } else {
718 std::cout << m_className << "::SetWeightingField:" << std::endl;
719 std::cout << " Reading weighting field from text file:" << prnsol
720 << std::endl;
721 // Buffer for reading
722 const int size = 100;
723 char line[size];
724
725 // Read the voltage list
726 int il = 0;
727
728 bool readerror = false;
729 while (fprnsol.getline(line, size, '\n')) {
730 il++;
731 // Split the line in tokens
732 char* token = strtok(line, " ");
733 // Skip blank lines and headers
734 if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
735 int(token[0]) == 10 || int(token[0]) == 13 ||
736 strcmp(token, "PRINT") == 0 || strcmp(token, "*****") == 0 ||
737 strcmp(token, "LOAD") == 0 || strcmp(token, "TIME=") == 0 ||
738 strcmp(token, "MAXIMUM") == 0 || strcmp(token, "VALUE") == 0 ||
739 strcmp(token, "NODE") == 0)
740 continue;
741 // Read the element
742 int inode = ReadInteger(token, -1, readerror);
743 token = strtok(nullptr, " ");
744 double volt = ReadDouble(token, -1, readerror);
745 try {
746 potentials.at(inode - 1) = volt;
747 nread++;
748 } catch (...) {
749 std::cerr << m_className << "::SetWeightingField:" << std::endl;
750 std::cerr << " Node number " << inode << " out of range."
751 << std::endl;
752 std::cerr << " on potential file " << prnsol << " (line " << il
753 << ")." << std::endl;
754 std::cerr << " Size of the potential vector is: "
755 << potentials.size() << std::endl;
756 ok = false;
757 }
758 }
759 // Close the file
760 fprnsol.close();
761 }
762 // Tell how many lines read
763 std::cout << m_className << "::SetWeightingField:" << std::endl;
764 std::cout << " Read " << nread << "/" << m_nNodes
765 << " (expected) potentials from file " << prnsol << "."
766 << std::endl;
767 // Check number of nodes
768 if (nread != (int)m_nNodes) {
769 std::cerr << m_className << "::SetWeightingField:" << std::endl;
770 std::cerr << " Number of nodes read (" << nread << ")"
771 << " on potential file (" << prnsol << ")" << std::endl;
772 std::cerr << " does not match the node list (" << m_nNodes << ")."
773 << std::endl;
774 ok = false;
775 }
776 if (!ok) {
777 std::cerr << m_className << "::SetWeightingField:" << std::endl;
778 std::cerr << " Field map could not be read "
779 << "and cannot be interpolated." << std::endl;
780 return false;
781 }
782
783 m_weightingFields[label] = potentials;
784 return true;
785}

◆ ShiftComponent()

void Garfield::ComponentCST::ShiftComponent ( const double xShift,
const double yShift,
const double zShift )

Definition at line 787 of file ComponentCST.cc.

788 {
789 std::transform(m_xlines.begin(), m_xlines.end(), m_xlines.begin(), [xShift](double x) { return x + xShift;});
790 std::transform(m_ylines.begin(), m_ylines.end(), m_ylines.begin(), [yShift](double x) { return x + yShift;});
791 std::transform(m_zlines.begin(), m_zlines.end(), m_zlines.begin(), [zShift](double x) { return x + zShift;});
792 SetRange();
793 UpdatePeriodicity();
794
795 std::cout << m_className << "::ShiftComponent:" << std::endl;
796 std::cout << " Shifted component in x-direction: " << xShift
797 << "\t y-direction: " << yShift << "\t z-direction: " << zShift
798 << std::endl;
799}

◆ WeightingField()

void Garfield::ComponentCST::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 815 of file ComponentCST.cc.

817 {
818 // Initial values
819 wx = wy = wz = 0;
820
821 // Do not proceed if not properly initialised.
822 if (!m_ready) return;
823
824 // Look for the label.
825 auto it = m_weightingFields.find(label);
826 if (it == m_weightingFields.end()) {
827 // Do not proceed if the requested weighting field does not exist.
828 std::cerr << "No weighting field named " << label << " found!\n";
829 return;
830 }
831
832 // Check if the weighting field is properly initialised.
833 if (m_weightingFields[label].empty()) return;
834
835 // Copy the coordinates
836 double x = xin, y = yin, z = zin;
837
838 // Map the coordinates onto field map coordinates and get indexes
839 bool mirrored[3];
840 unsigned int i, j, k;
841 double pos[3] = {0., 0., 0.};
842 if (!Coordinate2Index(x, y, z, i, j, k, pos, mirrored)) {
843 return;
844 }
845
846 double rx = (pos[0] - m_xlines.at(i)) / (m_xlines.at(i + 1) - m_xlines.at(i));
847 double ry = (pos[1] - m_ylines.at(j)) / (m_ylines.at(j + 1) - m_ylines.at(j));
848 double rz = (pos[2] - m_zlines.at(k)) / (m_zlines.at(k + 1) - m_zlines.at(k));
849
850 float fwx = 0., fwy = 0., fwz = 0.;
851 if (!disableFieldComponent[0])
852 fwx = GetFieldComponent(i, j, k, rx, ry, rz, 'x', (*it).second);
853 if (!disableFieldComponent[1])
854 fwy = GetFieldComponent(i, j, k, rx, ry, rz, 'y', (*it).second);
855 if (!disableFieldComponent[2])
856 fwz = GetFieldComponent(i, j, k, rx, ry, rz, 'z', (*it).second);
857
858 if (m_elementMaterial.size() > 0 && doShaping) {
859 ShapeField(fwx, fwy, fwz, rx, ry, rz, i, j, k, (*it).second);
860 }
861 if (mirrored[0]) fwx *= -1.f;
862 if (mirrored[1]) fwy *= -1.f;
863 if (mirrored[2]) fwz *= -1.f;
864 if (m_warning) PrintWarning("WeightingField");
865 if (m_materials.at(m_elementMaterial.at(Index2Element(i, j, k))).driftmedium) {
866 if (!disableFieldComponent[0]) wx = fwx;
867 if (!disableFieldComponent[1]) wy = fwy;
868 if (!disableFieldComponent[2]) wz = fwz;
869 }
870}
Garfield::ComponentFieldMap::PrintWarning
void PrintWarning(const std::string &header)
Definition ComponentFieldMap.cc:2773

◆ WeightingPotential()

double Garfield::ComponentCST::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 872 of file ComponentCST.cc.

874 {
875 // Do not proceed if not properly initialised.
876 if (!m_ready) return 0.;
877
878 // Look for the label.
879 auto it = m_weightingFields.find(label);
880 if (it == m_weightingFields.end()) {
881 // Do not proceed if the requested weighting field does not exist.
882 std::cerr << "No weighting field named " << label << " found!\n";
883 return 0.;
884 }
885
886 // Check if the weighting field is properly initialised.
887 if (m_weightingFields[label].empty()) return 0.;
888
889 // Copy the coordinates
890 double x = xin, y = yin, z = zin;
891
892 // Map the coordinates onto field map coordinates
893 bool mirrored[3];
894 unsigned int i, j, k;
895 double pos[3] = {0., 0., 0.};
896 if (!Coordinate2Index(x, y, z, i, j, k, pos, mirrored)) {
897 return 0.;
898 }
899 double rx = (pos[0] - m_xlines.at(i)) /
900 (m_xlines.at(i + 1) - m_xlines.at(i));
901 double ry = (pos[1] - m_ylines.at(j)) /
902 (m_ylines.at(j + 1) - m_ylines.at(j));
903 double rz = (pos[2] - m_zlines.at(k)) /
904 (m_zlines.at(k + 1) - m_zlines.at(k));
905
906 double potential = GetPotential(i, j, k, rx, ry, rz, (*it).second);
907
908 if (m_debug) {
909 std::cout << m_className << "::WeightingPotential:" << std::endl;
910 std::cout << " Global: (" << x << "," << y << "," << z << "),"
911 << std::endl;
912 std::cout << " Local: (" << rx << "," << ry << "," << rz
913 << ") in element with indexes: i=" << i << ", j=" << j
914 << ", k=" << k << std::endl;
915 std::cout << " Node xyzV:" << std::endl;
916 std::cout << "Node 0 position: " << Index2Node(i + 1, j, k)
917 << "\t potential: " << ((*it).second).at(Index2Node(i + 1, j, k))
918 << "Node 1 position: " << Index2Node(i + 1, j + 1, k)
919 << "\t potential: "
920 << ((*it).second).at(Index2Node(i + 1, j + 1, k))
921 << "Node 2 position: " << Index2Node(i, j + 1, k)
922 << "\t potential: " << ((*it).second).at(Index2Node(i, j + 1, k))
923 << "Node 3 position: " << Index2Node(i, j, k)
924 << "\t potential: " << ((*it).second).at(Index2Node(i, j, k))
925 << "Node 4 position: " << Index2Node(i + 1, j, k + 1)
926 << "\t potential: "
927 << ((*it).second).at(Index2Node(i + 1, j, k + 1))
928 << "Node 5 position: " << Index2Node(i + 1, j + 1, k + 1)
929 << "\t potential: "
930 << ((*it).second).at(Index2Node(i + 1, j + 1, k + 1))
931 << "Node 6 position: " << Index2Node(i, j + 1, k + 1)
932 << "\t potential: "
933 << ((*it).second).at(Index2Node(i, j + 1, k + 1))
934 << "Node 7 position: " << Index2Node(i, j, k + 1)
935 << "\t potential: " << ((*it).second).at(Index2Node(i, j, k))
936 << std::endl;
937 }
938 return potential;
939}
The documentation for this class was generated from the following files: