Garfield++ 3.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)
 
Medium * GetMedium (const double x, const double y, const double z) override
 Get the medium at a given location (x, y, z).
 
void GetNumberOfMeshLines (unsigned int &n_x, unsigned int &n_y, unsigned int &n_z)
 
void GetElementBoundaries (unsigned int element, double &xmin, double &xmax, double &ymin, double &ymax, double &zmin, double &zmax)
 
int GetElementMaterial (unsigned int element)
 
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 SetRange () override
 Calculate x, y, z, V and angular ranges.
 
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)
 
bool Coordinate2Index (const double x, const double y, const double z, unsigned int &i, unsigned int &j, unsigned int &k)
 
- Public Member Functions inherited from Garfield::ComponentFieldMap
 ComponentFieldMap ()
 Constructor.
 
virtual ~ComponentFieldMap ()
 Destructor.
 
virtual void SetRange ()
 Calculate x, y, z, V and angular ranges.
 
void PrintRange ()
 Show x, y, z, V and angular ranges.
 
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 GetVoltageRange (double &vmin, double &vmax) override
 Calculate the voltage range [V].
 
void PrintMaterials ()
 List all currently defined materials.
 
void DriftMedium (const unsigned int imat)
 Flag a field map material as a drift medium.
 
void NotDriftMedium (const unsigned int imat)
 Flag a field map materials as a non-drift medium.
 
unsigned int GetNumberOfMaterials () const
 Return the number of materials in the field map.
 
double GetPermittivity (const unsigned int imat) const
 Return the permittivity of a field map material.
 
double GetConductivity (const unsigned int imat) const
 Return the conductivity of a field map material.
 
void SetMedium (const unsigned int imat, Medium *medium)
 Associate a field map material with a Medium class.
 
Medium * GetMedium (const unsigned int i) const
 Return the Medium associated to a field map material.
 
unsigned int GetNumberOfMedia () const
 
int GetNumberOfElements () const
 Return the number of mesh elements.
 
bool GetElement (const unsigned int i, double &vol, double &dmin, double &dmax)
 Return the volume and aspect ratio of a mesh element.
 
void EnableCheckMapIndices ()
 
void DisableCheckMapIndices ()
 
void EnableDeleteBackgroundElements (const bool on=true)
 Option to eliminate mesh elements in conductors (default: on).
 
void EnableTetrahedralTreeForElementSearch (const bool on=true)
 
virtual Medium * GetMedium (const double x, const double y, const double z)
 Get the medium at a given location (x, y, z).
 
- Public Member Functions inherited from Garfield::ComponentBase
 ComponentBase ()
 Constructor.
 
virtual ~ComponentBase ()
 Destructor.
 
virtual void SetGeometry (GeometryBase *geo)
 Define the geometry.
 
virtual void Clear ()
 Reset.
 
virtual Medium * GetMedium (const double x, const double y, const double z)
 Get the medium at a given location (x, y, z).
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status)=0
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status)=0
 Calculate the drift field [V/cm] and potential [V] at (x, y, z).
 
virtual bool GetVoltageRange (double &vmin, double &vmax)=0
 Calculate the voltage range [V].
 
virtual void WeightingField (const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label)
 
virtual double WeightingPotential (const double x, const double y, const double z, const std::string &label)
 
virtual void DelayedWeightingField (const double x, const double y, const double z, const double t, double &wx, double &wy, double &wz, const std::string &label)
 
virtual 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.
 
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 IntegrateFlux (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)
 
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 DisablePeriodicityX ()
 
void EnablePeriodicityY (const bool on=true)
 Enable simple periodicity in the $y$ direction.
 
void DisablePeriodicityY ()
 
void EnablePeriodicityZ (const bool on=true)
 Enable simple periodicity in the $z$ direction.
 
void DisablePeriodicityZ ()
 
void EnableMirrorPeriodicityX (const bool on=true)
 Enable mirror periodicity in the $x$ direction.
 
void DisableMirrorPeriodicityX ()
 
void EnableMirrorPeriodicityY (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void DisableMirrorPeriodicityY ()
 
void EnableMirrorPeriodicityZ (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void DisableMirrorPeriodicityZ ()
 
void EnableAxialPeriodicityX (const bool on=true)
 Enable axial periodicity in the $x$ direction.
 
void DisableAxialPeriodicityX ()
 
void EnableAxialPeriodicityY (const bool on=true)
 Enable axial periodicity in the $y$ direction.
 
void DisableAxialPeriodicityY ()
 
void EnableAxialPeriodicityZ (const bool on=true)
 Enable axial periodicity in the $z$ direction.
 
void DisableAxialPeriodicityZ ()
 
void EnableRotationSymmetryX (const bool on=true)
 Enable rotation symmetry around the $x$ axis.
 
void DisableRotationSymmetryX ()
 
void EnableRotationSymmetryY (const bool on=true)
 Enable rotation symmetry around the $y$ axis.
 
void DisableRotationSymmetryY ()
 
void EnableRotationSymmetryZ (const bool on=true)
 Enable rotation symmetry around the $z$ axis.
 
void DisableRotationSymmetryZ ()
 
void EnableDebugging ()
 Switch on debugging messages.
 
void DisableDebugging ()
 Switch off debugging messages.
 
void ActivateTraps ()
 Request trapping to be taken care of by the component (for TCAD).
 
void DeactivateTraps ()
 
bool IsTrapActive ()
 
void ActivateVelocityMap ()
 Request velocity to be taken care of by the component (for TCAD).
 
void DectivateVelocityMap ()
 
bool IsVelocityActive ()
 
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 void ElectronVelocity (const double, const double, const double, double &vx, double &vy, double &vz, Medium *&, int &status)
 Get the electron drift velocity.
 
virtual void HoleVelocity (const double, const double, const double, double &vx, double &vy, double &vz, Medium *&, int &status)
 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.
 
double GetElementVolume (const unsigned int i) override
 
void GetAspectRatio (const unsigned int i, double &dmin, double &dmax) 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, double &rcoordinate, double &rotation)
 
- Protected Member Functions inherited from Garfield::ComponentFieldMap
void Reset () override
 Reset the component.
 
void UpdatePeriodicity2d ()
 
void UpdatePeriodicityCommon ()
 
int FindElement5 (const double x, const double y, const double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det)
 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)
 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)
 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, double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
 Move (ex, ey, ez) to global coordinates.
 
int ReadInteger (char *token, int def, bool &error)
 
double ReadDouble (char *token, double def, bool &error)
 
virtual double GetElementVolume (const unsigned int i)=0
 
virtual void GetAspectRatio (const unsigned int i, double &dmin, double &dmax)=0
 
void PrintWarning (const std::string &header)
 
void PrintNotReady (const std::string &header) 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 Element &element, const unsigned int n, const int iw=-1) const
 
virtual void Reset ()=0
 Reset the component.
 
virtual void UpdatePeriodicity ()=0
 Verify periodicities.
 

Additional Inherited Members

- Protected Attributes inherited from Garfield::ComponentFieldMap
bool m_is3d = true
 
int nElements = -1
 
std::vector< Element > elements
 
int nNodes = -1
 
std::vector< Node > nodes
 
unsigned int m_nMaterials = 0
 
std::vector< Material > materials
 
int nWeightingFields = 0
 
std::vector< std::string > wfields
 
std::vector< bool > wfieldsOk
 
bool hasBoundingBox = false
 
std::array< double, 3 > m_minBoundingBox
 
std::array< double, 3 > m_maxBoundingBox
 
std::array< double, 3 > m_mapmin
 
std::array< double, 3 > m_mapmax
 
std::array< double, 3 > m_mapamin
 
std::array< double, 3 > m_mapamax
 
std::array< double, 3 > m_mapna
 
std::array< double, 3 > m_cells
 
double m_mapvmin
 
double m_mapvmax
 
std::array< bool, 3 > m_setang
 
bool m_deleteBackground = true
 
bool m_warning = false
 
unsigned int m_nWarnings = 0
 
- Protected Attributes inherited from Garfield::ComponentBase
std::string m_className = "ComponentBase"
 Class name.
 
GeometryBase * m_geometry = nullptr
 Pointer to the geometry.
 
bool m_ready = false
 Ready for use?
 
bool m_activeTraps = false
 Does the component have traps?
 
bool m_hasVelocityMap = false
 Does the component have velocity maps?
 
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.
 
double m_bx0 = 0.
 
double m_by0 = 0.
 
double m_bz0 = 0.
 
bool m_debug = false
 Switch on/off debugging messages.
 

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 18 of file ComponentCST.cc.

18 : ComponentFieldMap() {
19 m_className = "ComponentCST";
20 // Default bounding box
21 m_minBoundingBox[2] = -50.;
22 m_maxBoundingBox[2] = 50.;
23 m_xlines.clear();
24 m_ylines.clear();
25 m_zlines.clear();
26 m_deleteBackground = false;
27 disableFieldComponent[0] = false;
28 disableFieldComponent[1] = false;
29 disableFieldComponent[2] = false;
30 doShaping = false;
31}
Garfield::ComponentBase::m_className
std::string m_className
Class name.
Definition: ComponentBase.hh:285
Garfield::ComponentFieldMap::m_maxBoundingBox
std::array< double, 3 > m_maxBoundingBox
Definition: ComponentFieldMap.hh:134
Garfield::ComponentFieldMap::m_minBoundingBox
std::array< double, 3 > m_minBoundingBox
Definition: ComponentFieldMap.hh:133

◆ ~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 
)

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 1210 of file ComponentCST.cc.

1212 {
1213 bool mirrored[3] = {false, false, false};
1214 double position_mapped[3] = {0., 0., 0.};
1215 double rcoordinate, rotation;
1216 return Coordinate2Index(x, y, z, i, j, k, position_mapped, mirrored,
1217 rcoordinate, rotation);
1218}
Garfield::ComponentCST::Coordinate2Index
bool Coordinate2Index(const double x, const double y, const double z, unsigned int &i, unsigned int &j, unsigned int &k)
Definition: ComponentCST.cc:1210

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,
double &  rcoordinate,
double &  rotation 
)
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.
rcoordinateInformation about rotation of the component. See ComponentFieldMap::MapCoordinates.
rotationInformation about rotation of the component. See ComponentFieldMap::MapCoordinates.

Definition at line 1228 of file ComponentCST.cc.

1232 {
1233 // Map the coordinates onto field map coordinates
1234 position_mapped[0] = xin;
1235 position_mapped[1] = yin;
1236 position_mapped[2] = zin;
1237 MapCoordinates(position_mapped[0], position_mapped[1], position_mapped[2],
1238 mirrored[0], mirrored[1], mirrored[2], rcoordinate, rotation);
1239
1240 auto it_x =
1241 std::lower_bound(m_xlines.begin(), m_xlines.end(), position_mapped[0]);
1242 auto it_y =
1243 std::lower_bound(m_ylines.begin(), m_ylines.end(), position_mapped[1]);
1244 auto it_z =
1245 std::lower_bound(m_zlines.begin(), m_zlines.end(), position_mapped[2]);
1246 if (it_x == m_xlines.end() || it_y == m_ylines.end() ||
1247 it_z == m_zlines.end() || position_mapped[0] < m_xlines.at(0) ||
1248 position_mapped[1] < m_ylines.at(0) ||
1249 position_mapped[2] < m_zlines.at(0)) {
1250 if (m_debug) {
1251 std::cerr << m_className << "::ElectricFieldBinary:" << std::endl;
1252 std::cerr << " Could not find the given coordinate!" << std::endl;
1253 std::cerr << " You ask for the following position: " << xin << ", "
1254 << yin << ", " << zin << std::endl;
1255 std::cerr << " The mapped position is: " << position_mapped[0] << ", "
1256 << position_mapped[1] << ", " << position_mapped[2]
1257 << std::endl;
1258 PrintRange();
1259 }
1260 return false;
1261 }
1262 /* Lower bound returns the next mesh line behind the position in question.
1263 * If the position in question is on a mesh line this mesh line is returned.
1264 * Since we are interested in the mesh line before the position in question we
1265 * need to move the
1266 * iterator to the left except for the very first mesh line!
1267 */
1268 if (it_x == m_xlines.begin())
1269 i = 0;
1270 else
1271 i = std::distance(m_xlines.begin(), it_x - 1);
1272 if (it_y == m_ylines.begin())
1273 j = 0;
1274 else
1275 j = std::distance(m_ylines.begin(), it_y - 1);
1276 if (it_z == m_zlines.begin())
1277 k = 0;
1278 else
1279 k = std::distance(m_zlines.begin(), it_z - 1);
1280 return true;
1281}
Garfield::ComponentBase::m_debug
bool m_debug
Switch on/off debugging messages.
Definition: ComponentBase.hh:310
Garfield::ComponentFieldMap::PrintRange
void PrintRange()
Show x, y, z, V and angular ranges.
Definition: ComponentFieldMap.cc:2031
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:2074

◆ DisableShaping()

void Garfield::ComponentCST::DisableShaping ( )
inline

Definition at line 177 of file ComponentCST.hh.

177{ 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 130 of file ComponentCST.hh.

130{ disableFieldComponent[0] = true; };

◆ DisableYField()

void Garfield::ComponentCST::DisableYField ( )
inline

Definition at line 131 of file ComponentCST.hh.

131{ disableFieldComponent[1] = true; };

◆ DisableZField()

void Garfield::ComponentCST::DisableZField ( )
inline

Definition at line 132 of file ComponentCST.hh.

132{ 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::ComponentBase.

Definition at line 993 of file ComponentCST.cc.

996 {
997 ElectricFieldBinary(xin, yin, zin, ex, ey, ez, volt, m, status, true);
998}

◆ 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::ComponentBase.

Definition at line 986 of file ComponentCST.cc.

988 {
989 double volt;
990 ElectricFieldBinary(xin, yin, zin, ex, ey, ez, volt, m, status);
991}

◆ 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 176 of file ComponentCST.hh.

176{ doShaping = true; }

◆ GetAspectRatio()

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

Implements Garfield::ComponentFieldMap.

Definition at line 1288 of file ComponentCST.cc.

1289 {
1290 if ((int)element >= nElements) {
1291 dmin = dmax = 0.;
1292 return;
1293 }
1294 unsigned int i, j, k;
1295 Element2Index(element, i, j, k);
1296 std::vector<double> distances;
1297 distances.push_back(m_xlines.at(i + 1) - m_xlines.at(i));
1298 distances.push_back(m_ylines.at(j + 1) - m_ylines.at(j));
1299 distances.push_back(m_zlines.at(k + 1) - m_zlines.at(k));
1300 std::sort(distances.begin(), distances.end());
1301 dmin = distances.at(0);
1302 dmax = distances.at(2);
1303}

◆ GetElementBoundaries()

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

Definition at line 1141 of file ComponentCST.cc.

1144 {
1145 unsigned int i, j, k;
1146 Element2Index(element, i, j, k);
1147 xmin = m_xlines.at(i);
1148 xmax = m_xlines.at(i + 1);
1149 ymin = m_ylines.at(j);
1150 ymax = m_ylines.at(j + 1);
1151 zmin = m_zlines.at(k);
1152 zmax = m_zlines.at(k + 1);
1153}

◆ GetElementMaterial()

int Garfield::ComponentCST::GetElementMaterial ( unsigned int  element)
inline

Definition at line 31 of file ComponentCST.hh.

31 {
32 return m_elementMaterial.at(element);
33 }

◆ GetElementVolume()

double Garfield::ComponentCST::GetElementVolume ( const unsigned int  i)
overrideprotectedvirtual

Implements Garfield::ComponentFieldMap.

Definition at line 1305 of file ComponentCST.cc.

1305 {
1306 if ((int)element >= nElements) return 0.;
1307 unsigned int i, j, k;
1308 Element2Index(element, i, j, k);
1309 const double volume = fabs((m_xlines.at(i + 1) - m_xlines.at(i)) *
1310 (m_xlines.at(j + 1) - m_ylines.at(j)) *
1311 (m_xlines.at(k + 1) - m_zlines.at(k)));
1312 return volume;
1313}
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615

◆ 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::ComponentBase.

Definition at line 1155 of file ComponentCST.cc.

1156 {
1157 unsigned int i, j, k;
1158 Coordinate2Index(xin, yin, zin, i, j, k);
1159 if (m_debug) {
1160 std::cout << m_className << "::GetMedium:" << std::endl;
1161 std::cout << " Found position (" << xin << ", " << yin << ", " << zin
1162 << "): " << std::endl;
1163 std::cout << " Indexes are: x: " << i << "/" << m_xlines.size()
1164 << "\t y: " << j << "/" << m_ylines.size() << "\t z: " << k << "/"
1165 << m_zlines.size() << std::endl;
1166 std::cout << " Element material index: " << Index2Element(i, j, k)
1167 << std::endl;
1168 std::cout << " Element index: "
1169 << (int)m_elementMaterial.at(Index2Element(i, j, k)) << std::endl;
1170 }
1171 return materials.at(m_elementMaterial.at(Index2Element(i, j, k))).medium;
1172}
Garfield::ComponentCST::Index2Element
int Index2Element(const unsigned int i, const unsigned int j, const unsigned int k)
Definition: ComponentCST.cc:1220
Garfield::ComponentFieldMap::materials
std::vector< Material > materials
Definition: ComponentFieldMap.hh:125

◆ GetNumberOfMeshLines()

void Garfield::ComponentCST::GetNumberOfMeshLines ( unsigned int &  n_x,
unsigned int &  n_y,
unsigned int &  n_z 
)

Definition at line 1134 of file ComponentCST.cc.

1135 {
1136 n_x = m_xlines.size();
1137 n_y = m_ylines.size();
1138 n_z = m_zlines.size();
1139}

◆ Index2Element()

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

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 1220 of file ComponentCST.cc.

1221 {
1222 if (i > m_nx - 2 || j > m_ny - 2 || k > m_nz - 2) {
1223 throw "FieldMap::ElementByIndex: Error. Element indexes out of bounds.";
1224 }
1225 return i + j * (m_nx - 1) + k * (m_nx - 1) * (m_ny - 1);
1226}

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 500 of file ComponentCST.cc.

500 {
501 m_ready = false;
502
503 // Keep track of the success
504 bool ok = true;
505 // Check the value of the unit
506 double funit;
507 if (strcmp(unit.c_str(), "mum") == 0 || strcmp(unit.c_str(), "micron") == 0 ||
508 strcmp(unit.c_str(), "micrometer") == 0) {
509 funit = 0.0001;
510 } else if (strcmp(unit.c_str(), "mm") == 0 ||
511 strcmp(unit.c_str(), "millimeter") == 0) {
512 funit = 0.1;
513 } else if (strcmp(unit.c_str(), "cm") == 0 ||
514 strcmp(unit.c_str(), "centimeter") == 0) {
515 funit = 1.0;
516 } else if (strcmp(unit.c_str(), "m") == 0 ||
517 strcmp(unit.c_str(), "meter") == 0) {
518 funit = 100.0;
519 } else {
520 std::cerr << m_className << "::Initialise:" << std::endl;
521 std::cerr << " Unknown length unit " << unit << "." << std::endl;
522 ok = false;
523 funit = 1.0;
524 }
525 if (m_debug) {
526 std::cout << m_className << "::Initialise:" << std::endl;
527 std::cout << " Unit scaling factor = " << funit << "." << std::endl;
528 }
529 FILE* f = fopen(dataFile.c_str(), "rb");
530 if (f == nullptr) {
531 std::cerr << m_className << "::Initilise:" << std::endl;
532 std::cerr << " Could not open file:" << dataFile.c_str() << std::endl;
533 return false;
534 }
535
536 struct stat fileStatus;
537 stat(dataFile.c_str(), &fileStatus);
538 int fileSize = fileStatus.st_size;
539
540 if (fileSize < 1000) {
541 fclose(f);
542 std::cerr << m_className << "::Initilise:" << std::endl;
543 std::cerr << " Error. The file is extremely short and does not seem to "
544 "contain a header or data."
545 << std::endl;
546 ok = false;
547 }
548
549 char header[headerSize];
550 size_t result;
551 result = fread(header, sizeof(char), headerSize, f);
552 if (result != headerSize) {
553 fputs("Reading error while reading header.", stderr);
554 exit(3);
555 }
556
557 int nx = 0, ny = 0, nz = 0;
558 int m_x = 0, m_y = 0, m_z = 0;
559 int n_s = 0, n_x = 0, n_y = 0, n_z = 0;
560 int e_s = 0, e_x = 0, e_y = 0, e_z = 0;
561 int e_m = 0;
562
563 int filled = 0;
564 filled = std::sscanf(
565 header,
566 (std::string("mesh_nx=%d mesh_ny=%d mesh_nz=%d\n") +
567 std::string("mesh_xlines=%d mesh_ylines=%d mesh_zlines=%d\n") +
568 std::string("nodes_scalar=%d nodes_vector_x=%d nodes_vector_y=%d "
569 "nodes_vector_z=%d\n") +
570 std::string("elements_scalar=%d elements_vector_x=%d "
571 "elements_vector_y=%d elements_vector_z=%d\n") +
572 std::string("elements_material=%d\n") + std::string("n_materials=%d\n"))
573 .c_str(),
574 &nx, &ny, &nz, &m_x, &m_y, &m_z, &n_s, &n_x, &n_y, &n_z, &e_s, &e_x, &e_y,
575 &e_z, &e_m, &m_nMaterials);
576 if (filled != 16) {
577 fclose(f);
578 std::cerr << m_className << "::Initilise:" << std::endl;
579 std::cerr << " Error. File header of " << dataFile.c_str()
580 << " is broken." << std::endl;
581 ok = false;
582 }
583 if (fileSize < 1000 + (m_x + m_y + m_z) * 8 +
584 (n_s + n_x + n_y + n_z + e_s + e_x + e_y + e_z) * 4 +
585 e_m * 1 + (int)m_nMaterials * 20) {
586 fclose(f);
587 ok = false;
588 }
589 if (m_debug) {
590 std::cout << m_className << "::Initialise:" << std::endl;
591 std::cout << " Information about the data stored in the given binary file:"
592 << std::endl;
593 std::cout << " Mesh (nx): " << nx << "\t Mesh (ny): " << ny
594 << "\t Mesh (nz): " << nz << std::endl;
595 std::cout << " Mesh (x_lines): " << m_x << "\t Mesh (y_lines): " << m_y
596 << "\t Mesh (z_lines): " << m_z << std::endl;
597 std::cout << " Nodes (scalar): " << n_s << "\t Nodes (x): " << n_x
598 << "\t Nodes (y): " << n_y << "\t Nodes (z): " << n_z
599 << std::endl;
600 std::cout << " Field (scalar): " << e_s << "\t Field (x): " << e_x
601 << "\t Field (y): " << e_y << "\t Field (z): " << e_z
602 << std::endl;
603 std::cout << " Elements: " << e_m << "\t Materials: " << m_nMaterials
604 << std::endl;
605 }
606 m_nx = m_x;
607 m_ny = m_y;
608 m_nz = m_z;
609 nNodes = m_nx * m_ny * m_nz;
610 nElements = (m_nx - 1) * (m_ny - 1) * (m_nz - 1);
611
612 m_xlines.resize(m_x);
613 m_ylines.resize(m_y);
614 m_zlines.resize(m_z);
615 m_potential.resize(n_s);
616 m_elementMaterial.resize(e_m);
617 // elements_scalar.resize(e_s);
618 materials.resize(m_nMaterials);
619 result = fread(m_xlines.data(), sizeof(double), m_xlines.size(), f);
620 if (result != m_xlines.size()) {
621 fputs("Reading error while reading xlines.", stderr);
622 exit(3);
623 } else if (result == 0) {
624 fputs("No xlines are stored in the data file.", stderr);
625 exit(3);
626 }
627 result = fread(m_ylines.data(), sizeof(double), m_ylines.size(), f);
628 if (result != m_ylines.size()) {
629 fputs("Reading error while reading ylines", stderr);
630 exit(3);
631 } else if (result == 0) {
632 fputs("No ylines are stored in the data file.", stderr);
633 exit(3);
634 }
635 result = fread(m_zlines.data(), sizeof(double), m_zlines.size(), f);
636 if (result != m_zlines.size()) {
637 fputs("Reading error while reasing zlines", stderr);
638 exit(3);
639 } else if (result == 0) {
640 fputs("No zlines are stored in the data file.", stderr);
641 exit(3);
642 }
643 result = fread(m_potential.data(), sizeof(float), m_potential.size(), f);
644 if (result != m_potential.size()) {
645 fputs("Reading error while reading nodes.", stderr);
646 exit(3);
647 } else if (result == 0) {
648 fputs("No potentials are stored in the data file.", stderr);
649 exit(3);
650 }
651 fseek(f, e_s * sizeof(float), SEEK_CUR);
652 // not needed in principle - thus it is ok if nothing is read
653 result = fread(m_elementMaterial.data(), sizeof(unsigned char),
654 m_elementMaterial.size(), f);
655 if (result != m_elementMaterial.size()) {
656 fputs("Reading error while reading element material", stderr);
657 exit(3);
658 }
659 std::stringstream st;
660 st << m_className << "::Initialise:" << std::endl;
661 /*
662 * The material vector is filled according to the material id!
663 * Thus material.at(0) is material with id 0.
664 */
665 for (unsigned int i = 0; i < materials.size(); i++) {
666 float id;
667 result = fread(&(id), sizeof(float), 1, f);
668 if (result != 1) {
669 fputs("Input error while reading material id.", stderr);
670 exit(3);
671 }
672 unsigned int description_size = 0;
673 result = fread(&(description_size), sizeof(int), 1, f);
674 if (result != 1) {
675 fputs("Input error while reading material description size.", stderr);
676 exit(3);
677 }
678 char* c = new char[description_size];
679 result = fread(c, sizeof(char), description_size, f);
680 if (result != description_size) {
681 fputs("Input error while reading material description.", stderr);
682 exit(3);
683 }
684 std::string name = c;
685 st << " Read material: " << name.c_str();
686 if (name.compare("gas") == 0) {
687 st << " (considered as drift medium)";
688 materials.at(id).driftmedium = true;
689 } else {
690 materials.at(id).driftmedium = false;
691 }
692 delete[] c;
693 float tmp_eps;
694 result = fread(&(tmp_eps), sizeof(float), 1, f);
695 materials.at(id).eps = tmp_eps;
696 if (result != 1) {
697 fputs("Reading error while reading eps.", stderr);
698 exit(3);
699 }
700 // float mue, rho;
701 // result = fread(&(mue), sizeof(float), 1, f);
702 // if (result != 1) {fputs ("Reading error while reading
703 //mue.",stderr);
704 // exit (3);}
705 // result = fread(&(rho), sizeof(float), 1, f);
706 // if (result != 1) {fputs ("Reading error while reading
707 //rho.",stderr);
708 // exit (3);}
709 st << "; eps is: " << materials.at(id).eps <<
710 // "\t mue is: " << mue <<
711 // "\t rho is: " << rho <<
712 "\t id is: " << id << std::endl;
713 // skip mue and rho
714 fseek(f, 2 * sizeof(float), SEEK_CUR);
715 // ToDo: Check if rho should be used to decide, which material is driftable
716 }
717 if (m_debug) {
718 std::cout << st.str();
719 for (auto it = materials.begin(), it_end = materials.end(); it != it_end;
720 it++) {
721 std::cout << "Material id: " << std::distance(materials.begin(), it)
722 << " \t driftable: " << (*it).driftmedium << std::endl;
723 }
724 }
725 // To be sure that they are sorted (should be already be the case)
726 std::sort(m_xlines.begin(), m_xlines.end());
727 std::sort(m_ylines.begin(), m_ylines.end());
728 std::sort(m_zlines.begin(), m_zlines.end());
729 if (funit != 1) {
730 std::transform(m_xlines.begin(), m_xlines.end(), m_xlines.begin(), [funit](double x) { return x * funit;});
731 std::transform(m_ylines.begin(), m_ylines.end(), m_ylines.begin(), [funit](double x) { return x * funit;});
732 std::transform(m_zlines.begin(), m_zlines.end(), m_zlines.begin(), [funit](double x) { return x * funit;});
733 }
734
735 std::cout << m_className << "::Initialise" << std::endl;
736 std::cout << " x range: " << *(m_xlines.begin()) << " - "
737 << *(m_xlines.end() - 1) << std::endl;
738 std::cout << " y range: " << *(m_ylines.begin()) << " - "
739 << *(m_ylines.end() - 1) << std::endl;
740 std::cout << " z range: " << *(m_zlines.begin()) << " - "
741 << *(m_zlines.end() - 1) << std::endl;
742 fclose(f);
743 // Set the ready flag
744 if (ok) {
745 m_ready = true;
746 } else {
747 std::cerr << m_className << "::Initialise:" << std::endl;
748 std::cerr << " Field map could not be read and cannot be interpolated."
749 << std::endl;
750 return false;
751 }
752
753 SetRange();
754 UpdatePeriodicity();
755 return true;
756}
Garfield::ComponentBase::m_ready
bool m_ready
Ready for use?
Definition: ComponentBase.hh:291
Garfield::ComponentCST::UpdatePeriodicity
void UpdatePeriodicity() override
Verify periodicities.
Definition: ComponentCST.cc:1283
Garfield::ComponentCST::SetRange
void SetRange() override
Calculate x, y, z, V and angular ranges.
Definition: ComponentCST.cc:1174

◆ 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 33 of file ComponentCST.cc.

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

◆ SetRange()

void Garfield::ComponentCST::SetRange ( )
overridevirtual

Calculate x, y, z, V and angular ranges.

Reimplemented from Garfield::ComponentFieldMap.

Definition at line 1174 of file ComponentCST.cc.

1174 {
1175 // Establish the ranges
1176 m_mapmin[0] = *m_xlines.begin();
1177 m_mapmax[0] = *(m_xlines.end() - 1);
1178 m_mapmin[1] = *m_ylines.begin();
1179 m_mapmax[1] = *(m_ylines.end() - 1);
1180 m_mapmin[2] = *m_zlines.begin();
1181 m_mapmax[2] = *(m_zlines.end() - 1);
1182 m_mapvmin = *std::min_element(m_potential.begin(), m_potential.end());
1183 m_mapvmax = *std::max_element(m_potential.begin(), m_potential.end());
1184
1185 // Set provisional cell dimensions.
1186 m_minBoundingBox[0] = m_mapmin[0];
1187 m_maxBoundingBox[0] = m_mapmax[0];
1188 m_minBoundingBox[1] = m_mapmin[1];
1189 m_maxBoundingBox[1] = m_mapmax[1];
1190 if (m_is3d) {
1191 m_minBoundingBox[2] = m_mapmin[2];
1192 m_maxBoundingBox[2] = m_mapmax[2];
1193 } else {
1194 m_mapmin[2] = m_minBoundingBox[2];
1195 m_mapmax[2] = m_maxBoundingBox[2];
1196 }
1197 hasBoundingBox = true;
1198}
Garfield::ComponentFieldMap::m_mapmin
std::array< double, 3 > m_mapmin
Definition: ComponentFieldMap.hh:137
Garfield::ComponentFieldMap::m_mapmax
std::array< double, 3 > m_mapmax
Definition: ComponentFieldMap.hh:138

Referenced by Initialise(), and ShiftComponent().

◆ SetRangeZ()

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

Definition at line 1200 of file ComponentCST.cc.

1200 {
1201 if (fabs(zmax - zmin) <= 0.) {
1202 std::cerr << m_className << "::SetRangeZ:" << std::endl;
1203 std::cerr << " Zero range is not permitted." << std::endl;
1204 return;
1205 }
1206 m_minBoundingBox[2] = std::min(zmin, zmax);
1207 m_maxBoundingBox[2] = std::max(zmin, zmax);
1208}

◆ 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 758 of file ComponentCST.cc.

759 {
760 std::vector<float> potentials(nNodes);
761 if (!m_ready) {
762 std::cerr << m_className << "::SetWeightingField:" << std::endl;
763 std::cerr << " No valid field map is present." << std::endl;
764 std::cerr << " Weighting field cannot be added." << std::endl;
765 return false;
766 }
767
768 // Open the voltage list
769 std::ifstream fprnsol;
770 fprnsol.open(prnsol.c_str(), std::ios::in);
771 if (fprnsol.fail()) {
772 std::cerr << m_className << "::SetWeightingField:" << std::endl;
773 std::cerr << " Could not open potential file " << prnsol
774 << " for reading." << std::endl;
775 std::cerr << " The file perhaps does not exist." << std::endl;
776 return false;
777 }
778 // Check if a weighting field with the same label already exists.
779 auto it = m_weightingFields.find(label);
780 if (it != m_weightingFields.end()) {
781 std::cout << m_className << "::SetWeightingField:" << std::endl;
782 std::cout << " Replacing existing weighting field " << label << "."
783 << std::endl;
784 } else {
785 wfields.push_back(label);
786 wfieldsOk.push_back(false);
787 }
788
789 if (std::distance(m_weightingFields.begin(), it) !=
790 std::distance(wfields.begin(),
791 find(wfields.begin(), wfields.end(), label))) {
792 std::cerr << m_className << "::SetWeightingField:" << std::endl;
793 std::cerr << " Indexes of the weighting fields and the weighting field "
794 "counter are not equal!"
795 << std::endl;
796 return false;
797 }
798 unsigned int iField = std::distance(m_weightingFields.begin(), it);
799 int nread = 0;
800 bool ok = true;
801
802 if (isBinary) {
803 std::cout << m_className << "::SetWeightingField:" << std::endl;
804 std::cout << " Reading weighting field from binary file:"
805 << prnsol.c_str() << std::endl;
806 FILE* f = fopen(prnsol.c_str(), "rb");
807 if (f == nullptr) {
808 std::cerr << m_className << "::Initilise:" << std::endl;
809 std::cerr << " Could not open file:" << prnsol.c_str() << std::endl;
810 return false;
811 }
812
813 struct stat fileStatus;
814 stat(prnsol.c_str(), &fileStatus);
815 int fileSize = fileStatus.st_size;
816
817 if (fileSize < 1000) {
818 fclose(f);
819 std::cerr << m_className << "::SetWeightingField:" << std::endl;
820 std::cerr << " Error. The file is extremely short and does not seem "
821 "to contain a header or data."
822 << std::endl;
823 ok = false;
824 }
825
826 char header[headerSize];
827 size_t result;
828 result = fread(header, sizeof(char), headerSize, f);
829 if (result != headerSize) {
830 fputs("Reading error while reading header.", stderr);
831 exit(3);
832 }
833
834 int nx = 0, ny = 0, nz = 0;
835 int m_x = 0, m_y = 0, m_z = 0;
836 int n_x = 0, n_y = 0, n_z = 0;
837 int e_s = 0, e_x = 0, e_y = 0, e_z = 0;
838 int e_m = 0;
839
840 int filled = 0;
841 filled = std::sscanf(
842 header,
843 (std::string("mesh_nx=%d mesh_ny=%d mesh_nz=%d\n") +
844 std::string("mesh_xlines=%d mesh_ylines=%d mesh_zlines=%d\n") +
845 std::string("nodes_scalar=%d nodes_vector_x=%d nodes_vector_y=%d "
846 "nodes_vector_z=%d\n") +
847 std::string("elements_scalar=%d elements_vector_x=%d "
848 "elements_vector_y=%d elements_vector_z=%d\n") +
849 std::string("elements_material=%d\n") +
850 std::string("n_materials=%d\n"))
851 .c_str(),
852 &nx, &ny, &nz, &m_x, &m_y, &m_z, &nread, &n_x, &n_y, &n_z, &e_s, &e_x,
853 &e_y, &e_z, &e_m, &m_nMaterials);
854 if (filled != 16) {
855 fclose(f);
856 std::cerr << m_className << "::SetWeightingField:" << std::endl;
857 std::cerr << " Error. File header of " << prnsol.c_str()
858 << " is broken." << std::endl;
859 ok = false;
860 }
861 if (fileSize < 1000 + (m_x + m_y + m_z) * 8 +
862 (nread + n_x + n_y + n_z + e_s + e_x + e_y + e_z) * 4 +
863 e_m * 1 + (int)m_nMaterials * 20) {
864 fclose(f);
865 ok = false;
866 }
867 if (m_debug) {
868 std::cout << m_className << "::SetWeightingField:" << std::endl;
869 std::cout
870 << " Information about the data stored in the given binary file:"
871 << std::endl;
872 std::cout << " Mesh (nx): " << nx << "\t Mesh (ny): " << ny
873 << "\t Mesh (nz): " << nz << std::endl;
874 std::cout << " Mesh (x_lines): " << m_x << "\t Mesh (y_lines): " << m_y
875 << "\t Mesh (z_lines): " << m_z << std::endl;
876 std::cout << " Nodes (scalar): " << nread << "\t Nodes (x): " << n_x
877 << "\t Nodes (y): " << n_y << "\t Nodes (z): " << n_z
878 << std::endl;
879 std::cout << " Field (scalar): " << e_s << "\t Field (x): " << e_x
880 << "\t Field (y): " << e_y << "\t Field (z): " << e_z
881 << std::endl;
882 std::cout << " Elements: " << e_m << "\t Materials: " << m_nMaterials
883 << std::endl;
884 }
885 // skip everything, but the potential
886 fseek(f, m_x * sizeof(double), SEEK_CUR);
887 fseek(f, m_y * sizeof(double), SEEK_CUR);
888 fseek(f, m_z * sizeof(double), SEEK_CUR);
889 result = fread(potentials.data(), sizeof(float), potentials.size(), f);
890 if (result != potentials.size()) {
891 fputs("Reading error while reading nodes.", stderr);
892 exit(3);
893 } else if (result == 0) {
894 fputs("No wighting potentials are stored in the data file.", stderr);
895 exit(3);
896 }
897 fprnsol.close();
898 } else {
899 std::cout << m_className << "::SetWeightingField:" << std::endl;
900 std::cout << " Reading weighting field from text file:" << prnsol.c_str()
901 << std::endl;
902 // Buffer for reading
903 const int size = 100;
904 char line[size];
905
906 // Read the voltage list
907 int il = 0;
908
909 bool readerror = false;
910 while (fprnsol.getline(line, size, '\n')) {
911 il++;
912 // Split the line in tokens
913 char* token = NULL;
914 token = strtok(line, " ");
915 // Skip blank lines and headers
916 if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
917 int(token[0]) == 10 || int(token[0]) == 13 ||
918 strcmp(token, "PRINT") == 0 || strcmp(token, "*****") == 0 ||
919 strcmp(token, "LOAD") == 0 || strcmp(token, "TIME=") == 0 ||
920 strcmp(token, "MAXIMUM") == 0 || strcmp(token, "VALUE") == 0 ||
921 strcmp(token, "NODE") == 0)
922 continue;
923 // Read the element
924 int inode = ReadInteger(token, -1, readerror);
925 token = strtok(NULL, " ");
926 double volt = ReadDouble(token, -1, readerror);
927 try {
928 potentials.at(inode - 1) = volt;
929 nread++;
930 } catch (...) {
931 std::cerr << m_className << "::SetWeightingField:" << std::endl;
932 std::cerr << " Node number " << inode << " out of range."
933 << std::endl;
934 std::cerr << " on potential file " << prnsol << " (line " << il
935 << ")." << std::endl;
936 std::cerr << " Size of the potential vector is: "
937 << potentials.size() << std::endl;
938 ok = false;
939 }
940 }
941 // Close the file
942 fprnsol.close();
943 }
944 // Tell how many lines read
945 std::cout << m_className << "::SetWeightingField:" << std::endl;
946 std::cout << " Read " << nread << "/" << nNodes
947 << " (expected) potentials from file " << prnsol << "."
948 << std::endl;
949 // Check number of nodes
950 if (nread != nNodes) {
951 std::cerr << m_className << "::SetWeightingField:" << std::endl;
952 std::cerr << " Number of nodes read (" << nread << ")"
953 << " on potential file (" << prnsol << ")" << std::endl;
954 std::cerr << " does not match the node list (" << nNodes << ")."
955 << std::endl;
956 ok = false;
957 }
958 if (!ok) {
959 std::cerr << m_className << "::SetWeightingField:" << std::endl;
960 std::cerr << " Field map could not be read "
961 << "and cannot be interpolated." << std::endl;
962 return false;
963 }
964
965 m_weightingFields[label] = potentials;
966
967 // Set the ready flag.
968 wfieldsOk[iField] = ok;
969 return true;
970}
Garfield::ComponentFieldMap::wfields
std::vector< std::string > wfields
Definition: ComponentFieldMap.hh:128

◆ ShiftComponent()

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

Definition at line 972 of file ComponentCST.cc.

973 {
974 std::transform(m_xlines.begin(), m_xlines.end(), m_xlines.begin(), [xShift](double x) { return x + xShift;});
975 std::transform(m_ylines.begin(), m_ylines.end(), m_ylines.begin(), [yShift](double x) { return x + yShift;});
976 std::transform(m_zlines.begin(), m_zlines.end(), m_zlines.begin(), [zShift](double x) { return x + zShift;});
977 SetRange();
978 UpdatePeriodicity();
979
980 std::cout << m_className << "::ShiftComponent:" << std::endl;
981 std::cout << " Shifted component in x-direction: " << xShift
982 << "\t y-direction: " << yShift << "\t z-direction: " << zShift
983 << std::endl;
984}

◆ UpdatePeriodicity()

void Garfield::ComponentCST::UpdatePeriodicity ( )
overrideprotectedvirtual

Verify periodicities.

Implements Garfield::ComponentBase.

Definition at line 1283 of file ComponentCST.cc.

1283 {
1284 UpdatePeriodicity2d();
1285 UpdatePeriodicityCommon();
1286}

Referenced by Initialise(), and ShiftComponent().

◆ 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::ComponentBase.

Definition at line 1000 of file ComponentCST.cc.

1002 {
1003 // Initial values
1004 wx = wy = wz = 0;
1005
1006 // Do not proceed if not properly initialised.
1007 if (!m_ready) return;
1008
1009 // Look for the label.
1010 auto it = m_weightingFields.find(label);
1011 if (it == m_weightingFields.end()) {
1012 // Do not proceed if the requested weighting field does not exist.
1013 std::cerr << "No weighting field named " << label.c_str() << " found!"
1014 << std::endl;
1015 return;
1016 }
1017
1018 // Check if the weighting field is properly initialised.
1019 if (!wfieldsOk[std::distance(m_weightingFields.begin(), it)]) return;
1020
1021 // Copy the coordinates
1022 double x = xin, y = yin, z = zin;
1023
1024 // Map the coordinates onto field map coordinates and get indexes
1025 bool mirrored[3];
1026 double rcoordinate, rotation;
1027 unsigned int i, j, k;
1028 double position_mapped[3] = {0., 0., 0.};
1029 if (!Coordinate2Index(x, y, z, i, j, k, position_mapped, mirrored,
1030 rcoordinate, rotation))
1031 return;
1032
1033 double rx = (position_mapped[0] - m_xlines.at(i)) /
1034 (m_xlines.at(i + 1) - m_xlines.at(i));
1035 double ry = (position_mapped[1] - m_ylines.at(j)) /
1036 (m_ylines.at(j + 1) - m_ylines.at(j));
1037 double rz = (position_mapped[2] - m_zlines.at(k)) /
1038 (m_zlines.at(k + 1) - m_zlines.at(k));
1039
1040 float fwx, fwy, fwz;
1041 if (!disableFieldComponent[0])
1042 fwx = GetFieldComponent(i, j, k, rx, ry, rz, 'x', &((*it).second));
1043 if (!disableFieldComponent[1])
1044 fwy = GetFieldComponent(i, j, k, rx, ry, rz, 'y', &((*it).second));
1045 if (!disableFieldComponent[2])
1046 fwz = GetFieldComponent(i, j, k, rx, ry, rz, 'z', &((*it).second));
1047
1048 if (m_elementMaterial.size() > 0 && doShaping) {
1049 ShapeField(fwx, fwy, fwz, rx, ry, rz, i, j, k, &((*it).second));
1050 }
1051 if (mirrored[0]) fwx *= -1.;
1052 if (mirrored[1]) fwy *= -1.;
1053 if (mirrored[2]) fwz *= -1.;
1054 if (m_warning) PrintWarning("WeightingField");
1055 if (materials.at(m_elementMaterial.at(Index2Element(i, j, k))).driftmedium) {
1056 if (!disableFieldComponent[0]) wx = fwx;
1057 if (!disableFieldComponent[1]) wy = fwy;
1058 if (!disableFieldComponent[2]) wz = fwz;
1059 }
1060}
Garfield::ComponentFieldMap::PrintWarning
void PrintWarning(const std::string &header)
Definition: ComponentFieldMap.hh:193

◆ 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::ComponentBase.

Definition at line 1062 of file ComponentCST.cc.

1064 {
1065 // Do not proceed if not properly initialised.
1066 if (!m_ready) return 0.;
1067
1068 // Look for the label.
1069 auto it = m_weightingFields.find(label);
1070 if (it == m_weightingFields.end()) {
1071 // Do not proceed if the requested weighting field does not exist.
1072 std::cerr << "No weighting field named " << label.c_str() << " found!"
1073 << std::endl;
1074 return 0.;
1075 }
1076
1077 // Check if the weighting field is properly initialised.
1078 if (!wfieldsOk[std::distance(m_weightingFields.begin(), it)]) return 0.;
1079
1080 // Copy the coordinates
1081 double x = xin, y = yin, z = zin;
1082
1083 // Map the coordinates onto field map coordinates
1084 bool mirrored[3];
1085 double rcoordinate, rotation;
1086 unsigned int i, j, k;
1087 double position_mapped[3] = {0., 0., 0.};
1088 if (!Coordinate2Index(x, y, z, i, j, k, position_mapped, mirrored,
1089 rcoordinate, rotation))
1090 return 0.;
1091
1092 double rx = (position_mapped[0] - m_xlines.at(i)) /
1093 (m_xlines.at(i + 1) - m_xlines.at(i));
1094 double ry = (position_mapped[1] - m_ylines.at(j)) /
1095 (m_ylines.at(j + 1) - m_ylines.at(j));
1096 double rz = (position_mapped[2] - m_zlines.at(k)) /
1097 (m_zlines.at(k + 1) - m_zlines.at(k));
1098
1099 double potential = GetPotential(i, j, k, rx, ry, rz, &((*it).second));
1100
1101 if (m_debug) {
1102 std::cout << m_className << "::WeightingPotential:" << std::endl;
1103 std::cout << " Global: (" << x << "," << y << "," << z << "),"
1104 << std::endl;
1105 std::cout << " Local: (" << rx << "," << ry << "," << rz
1106 << ") in element with indexes: i=" << i << ", j=" << j
1107 << ", k=" << k << std::endl;
1108 std::cout << " Node xyzV:" << std::endl;
1109 std::cout << "Node 0 position: " << Index2Node(i + 1, j, k)
1110 << "\t potential: " << ((*it).second).at(Index2Node(i + 1, j, k))
1111 << "Node 1 position: " << Index2Node(i + 1, j + 1, k)
1112 << "\t potential: "
1113 << ((*it).second).at(Index2Node(i + 1, j + 1, k))
1114 << "Node 2 position: " << Index2Node(i, j + 1, k)
1115 << "\t potential: " << ((*it).second).at(Index2Node(i, j + 1, k))
1116 << "Node 3 position: " << Index2Node(i, j, k)
1117 << "\t potential: " << ((*it).second).at(Index2Node(i, j, k))
1118 << "Node 4 position: " << Index2Node(i + 1, j, k + 1)
1119 << "\t potential: "
1120 << ((*it).second).at(Index2Node(i + 1, j, k + 1))
1121 << "Node 5 position: " << Index2Node(i + 1, j + 1, k + 1)
1122 << "\t potential: "
1123 << ((*it).second).at(Index2Node(i + 1, j + 1, k + 1))
1124 << "Node 6 position: " << Index2Node(i, j + 1, k + 1)
1125 << "\t potential: "
1126 << ((*it).second).at(Index2Node(i, j + 1, k + 1))
1127 << "Node 7 position: " << Index2Node(i, j, k + 1)
1128 << "\t potential: " << ((*it).second).at(Index2Node(i, j, k))
1129 << std::endl;
1130 }
1131 return potential;
1132}
The documentation for this class was generated from the following files: