Garfield++ v2r0
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 ()
 
 ~ComponentCST ()
 
void ShiftComponent (const double xShift, const double yShift, const double zShift)
 
Medium * GetMedium (const double x, const double y, const double z)
 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)
 
void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status)
 
void WeightingField (const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label)
 
double WeightingPotential (const double x, const double y, const double z, const std::string &label)
 
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)
 
virtual bool IsInBoundingBox (const double x, const double y, const double z) const
 
void SetRange ()
 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.
 
virtual bool IsInBoundingBox (const double x, const double y, const double z) const
 
virtual bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
 Get the bounding box coordinates.
 
virtual bool GetVoltageRange (double &vmin, double &vmax)
 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.
 
Medium * GetMedium (const double x, const double y, const double z)=0
 Get the medium at a given location (x, y, z).
 
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.
 
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
 
virtual void WeightingField (const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label)=0
 
virtual double WeightingPotential (const double x, const double y, const double z, const std::string &label)=0
 
void EnableCheckMapIndices ()
 
void DisableCheckMapIndices ()
 
void EnableDeleteBackgroundElements ()
 
void DisableDeleteBackgroundElements ()
 
void EnableTetrahedralTreeForElementSearch (const bool on=true)
 
- 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
 
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 MagneticField (const double x, const double y, const double z, double &bx, double &by, double &bz, int &status)
 
void SetMagneticField (const double bx, const double by, const double bz)
 Set a constant magnetic field.
 
virtual bool IsReady ()
 Ready for use?
 
virtual bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
 Get the bounding box coordinates.
 
virtual bool 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)
 
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 ()
 Verify periodicities.
 
double GetElementVolume (const unsigned int i)
 
void GetAspectRatio (const unsigned int i, double &dmin, double &dmax)
 
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 ()
 Geometry checks.
 
virtual void UpdatePeriodicity ()=0
 Verify periodicities.
 
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 unsigned int i, const unsigned int n, const int iw=-1) const
 
virtual void Reset ()=0
 Geometry checks.
 
virtual void UpdatePeriodicity ()=0
 Verify periodicities.
 

Additional Inherited Members

- Protected Attributes inherited from Garfield::ComponentFieldMap
bool m_is3d
 
int nElements
 
std::vector< Element > elements
 
int nNodes
 
std::vector< Node > nodes
 
unsigned int m_nMaterials
 
std::vector< Material > materials
 
int nWeightingFields
 
std::vector< std::string > wfields
 
std::vector< bool > wfieldsOk
 
bool hasBoundingBox
 
double xMinBoundingBox
 
double yMinBoundingBox
 
double zMinBoundingBox
 
double xMaxBoundingBox
 
double yMaxBoundingBox
 
double zMaxBoundingBox
 
double mapxmin
 
double mapymin
 
double mapzmin
 
double mapxmax
 
double mapymax
 
double mapzmax
 
double mapxamin
 
double mapyamin
 
double mapzamin
 
double mapxamax
 
double mapyamax
 
double mapzamax
 
double mapvmin
 
double mapvmax
 
bool setangx
 
bool setangy
 
bool setangz
 
double mapsx
 
double mapsy
 
double mapsz
 
double cellsx
 
double cellsy
 
double cellsz
 
double mapnxa
 
double mapnya
 
double mapnza
 
bool m_deleteBackground
 
bool m_warning
 
unsigned int m_nWarnings
 
- Protected Attributes inherited from Garfield::ComponentBase
std::string m_className
 Class name.
 
GeometryBase * m_geometry
 Pointer to the geometry.
 
bool m_ready
 Ready for use?
 
bool m_activeTraps
 Does the component have traps?
 
bool m_hasVelocityMap
 Does the component have velocity maps?
 
bool m_xPeriodic
 Simple periodicity in x.
 
bool m_yPeriodic
 Simple periodicity in y.
 
bool m_zPeriodic
 Simple periodicity in z.
 
bool m_xMirrorPeriodic
 Mirror periodicity in x.
 
bool m_yMirrorPeriodic
 Mirror periodicity in y.
 
bool m_zMirrorPeriodic
 Mirror periodicity in z.
 
bool m_xAxiallyPeriodic
 Axial periodicity in x.
 
bool m_yAxiallyPeriodic
 Axial periodicity in y.
 
bool m_zAxiallyPeriodic
 Axial periodicity in z.
 
bool m_xRotationSymmetry
 Rotation symmetry around x-axis.
 
bool m_yRotationSymmetry
 Rotation symmetry around y-axis.
 
bool m_zRotationSymmetry
 Rotation symmetry around z-axis.
 
double m_bx0
 
double m_by0
 
double m_bz0
 
bool m_debug
 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 ( )

Definition at line 18 of file ComponentCST.cc.

18 : ComponentFieldMap() {
19
20 m_className = "ComponentCST";
21 // Default bounding box
22 zMinBoundingBox = -50.;
23 zMaxBoundingBox = 50.;
24 m_xlines.clear();
25 m_ylines.clear();
26 m_zlines.clear();
27 m_deleteBackground = false;
28 disableFieldComponent[0] = false;
29 disableFieldComponent[1] = false;
30 disableFieldComponent[2] = false;
31 doShaping = false;
32}
Garfield::ComponentBase::m_className
std::string m_className
Class name.
Definition: ComponentBase.hh:233

◆ ~ComponentCST()

Garfield::ComponentCST::~ComponentCST ( )
inline

Definition at line 21 of file ComponentCST.hh.

21{}

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

1228 {
1229 bool mirrored[3] = {false, false, false};
1230 double position_mapped[3] = {0., 0., 0.};
1231 double rcoordinate, rotation;
1232 return Coordinate2Index(x, y, z, i, j, k, position_mapped, mirrored,
1233 rcoordinate, rotation);
1234}
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:1226

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

1249 {
1250 // Map the coordinates onto field map coordinates
1251 position_mapped[0] = xin;
1252 position_mapped[1] = yin;
1253 position_mapped[2] = zin;
1254 MapCoordinates(position_mapped[0], position_mapped[1], position_mapped[2],
1255 mirrored[0], mirrored[1], mirrored[2], rcoordinate, rotation);
1256
1257 std::vector<double>::iterator it_x, it_y, it_z;
1258 it_x = std::lower_bound(m_xlines.begin(), m_xlines.end(), position_mapped[0]);
1259 it_y = std::lower_bound(m_ylines.begin(), m_ylines.end(), position_mapped[1]);
1260 it_z = std::lower_bound(m_zlines.begin(), m_zlines.end(), position_mapped[2]);
1261 if (it_x == m_xlines.end() || it_y == m_ylines.end() ||
1262 it_z == m_zlines.end() || position_mapped[0] < m_xlines.at(0) ||
1263 position_mapped[1] < m_ylines.at(0) ||
1264 position_mapped[2] < m_zlines.at(0)) {
1265 if (m_debug) {
1266 std::cerr << m_className << "::ElectricFieldBinary:" << std::endl;
1267 std::cerr << " Could not find the given coordinate!" << std::endl;
1268 std::cerr << " You ask for the following position: " << xin << ", "
1269 << yin << ", " << zin << std::endl;
1270 std::cerr << " The mapped position is: " << position_mapped[0] << ", "
1271 << position_mapped[1] << ", " << position_mapped[2]
1272 << std::endl;
1273 PrintRange();
1274 }
1275 return false;
1276 }
1277 /* Lower bound returns the next mesh line behind the position in question.
1278 * If the position in question is on a mesh line this mesh line is returned.
1279 * Since we are interested in the mesh line before the position in question we
1280 * need to move the
1281 * iterator to the left except for the very first mesh line!
1282 */
1283 if (it_x == m_xlines.begin())
1284 i = 0;
1285 else
1286 i = std::distance(m_xlines.begin(), it_x - 1);
1287 if (it_y == m_ylines.begin())
1288 j = 0;
1289 else
1290 j = std::distance(m_ylines.begin(), it_y - 1);
1291 if (it_z == m_zlines.begin())
1292 k = 0;
1293 else
1294 k = std::distance(m_zlines.begin(), it_z - 1);
1295 return true;
1296}
Garfield::ComponentBase::m_debug
bool m_debug
Switch on/off debugging messages.
Definition: ComponentBase.hh:274
Garfield::ComponentFieldMap::PrintRange
void PrintRange()
Show x, y, z, V and angular ranges.
Definition: ComponentFieldMap.cc:2354
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:2441

◆ DisableShaping()

void Garfield::ComponentCST::DisableShaping ( )
inline

Definition at line 167 of file ComponentCST.hh.

167 {
168 doShaping = false;
169 }

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

126 {
127 disableFieldComponent[0] = true;
128 };

◆ DisableYField()

void Garfield::ComponentCST::DisableYField ( )
inline

Definition at line 129 of file ComponentCST.hh.

129 {
130 disableFieldComponent[1] = true;
131 };

◆ DisableZField()

void Garfield::ComponentCST::DisableZField ( )
inline

Definition at line 132 of file ComponentCST.hh.

132 {
133 disableFieldComponent[2] = true;
134 };

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

Implements Garfield::ComponentFieldMap.

Definition at line 1001 of file ComponentCST.cc.

1004 {
1005 ElectricFieldBinary(xin, yin, zin, ex, ey, ez, volt, m, status, true);
1006}

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

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::ComponentFieldMap.

Definition at line 994 of file ComponentCST.cc.

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

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

164 {
165 doShaping = true;
166 }

◆ GetAspectRatio()

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

Implements Garfield::ComponentFieldMap.

Definition at line 1304 of file ComponentCST.cc.

1305 {
1306
1307 if ((int)element >= nElements) {
1308 dmin = dmax = 0.;
1309 return;
1310 }
1311 unsigned int i, j, k;
1312 Element2Index(element, i, j, k);
1313 std::vector<double> distances;
1314 distances.push_back(m_xlines.at(i + 1) - m_xlines.at(i));
1315 distances.push_back(m_ylines.at(j + 1) - m_ylines.at(j));
1316 distances.push_back(m_zlines.at(k + 1) - m_zlines.at(k));
1317 std::sort(distances.begin(), distances.end());
1318 dmin = distances.at(0);
1319 dmax = distances.at(2);
1320}

◆ GetElementBoundaries()

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

Definition at line 1153 of file ComponentCST.cc.

1156 {
1157 unsigned int i, j, k;
1158 Element2Index(element, i, j, k);
1159 xmin = m_xlines.at(i);
1160 xmax = m_xlines.at(i + 1);
1161 ymin = m_ylines.at(j);
1162 ymax = m_ylines.at(j + 1);
1163 zmin = m_zlines.at(k);
1164 zmax = m_zlines.at(k + 1);
1165}

◆ GetElementMaterial()

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

Definition at line 29 of file ComponentCST.hh.

29{return m_elementMaterial.at(element);}

◆ GetElementVolume()

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

Implements Garfield::ComponentFieldMap.

Definition at line 1322 of file ComponentCST.cc.

1322 {
1323
1324 if ((int)element >= nElements) return 0.;
1325 unsigned int i, j, k;
1326 Element2Index(element, i, j, k);
1327 const double volume = fabs((m_xlines.at(i + 1) - m_xlines.at(i)) *
1328 (m_xlines.at(j + 1) - m_ylines.at(j)) *
1329 (m_xlines.at(k + 1) - m_zlines.at(k)));
1330 return volume;
1331}
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 
)
virtual

Get the medium at a given location (x, y, z).

Implements Garfield::ComponentFieldMap.

Definition at line 1167 of file ComponentCST.cc.

1168 {
1169 unsigned int i, j, k;
1170 Coordinate2Index(xin, yin, zin, i, j, k);
1171 if (m_debug) {
1172 std::cout << m_className << "::GetMedium:" << std::endl;
1173 std::cout << " Found position (" << xin << ", " << yin << ", " << zin
1174 << "): " << std::endl;
1175 std::cout << " Indexes are: x: " << i << "/" << m_xlines.size()
1176 << "\t y: " << j << "/" << m_ylines.size() << "\t z: " << k << "/"
1177 << m_zlines.size() << std::endl;
1178 std::cout << " Element material index: " << Index2Element(i, j, k)
1179 << std::endl;
1180 std::cout << " Element index: "
1181 << (int)m_elementMaterial.at(Index2Element(i, j, k)) << std::endl;
1182 }
1183 return materials.at(m_elementMaterial.at(Index2Element(i, j, k))).medium;
1184}
Garfield::ComponentCST::Index2Element
int Index2Element(const unsigned int i, const unsigned int j, const unsigned int k)
Definition: ComponentCST.cc:1236
Garfield::ComponentFieldMap::materials
std::vector< Material > materials
Definition: ComponentFieldMap.hh:133

◆ GetNumberOfMeshLines()

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

Definition at line 1146 of file ComponentCST.cc.

1147 {
1148 n_x = m_xlines.size();
1149 n_y = m_ylines.size();
1150 n_z = m_zlines.size();
1151}

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

1237 {
1238
1239 if (i > m_nx - 2 || j > m_ny - 2 || k > m_nz - 2) {
1240 throw "FieldMap::ElementByIndex: Error. Element indexes out of bounds.";
1241 }
1242 return i + j * (m_nx - 1) + k * (m_nx - 1) * (m_ny - 1);
1243}

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

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

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

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

◆ IsInBoundingBox()

virtual bool Garfield::ComponentCST::IsInBoundingBox ( const double  x,
const double  y,
const double  z 
) const
inlinevirtual

Reimplemented from Garfield::ComponentFieldMap.

Definition at line 107 of file ComponentCST.hh.

108 {
109 return x >= xMinBoundingBox && x <= xMaxBoundingBox &&
110 y >= yMinBoundingBox && y <= yMaxBoundingBox &&
111 z >= zMinBoundingBox && z <= zMaxBoundingBox;
112 }

◆ SetRange()

void Garfield::ComponentCST::SetRange ( )
virtual

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

Reimplemented from Garfield::ComponentFieldMap.

Definition at line 1186 of file ComponentCST.cc.

1186 {
1187 // Establish the ranges
1188 mapxmin = *m_xlines.begin();
1189 mapxmax = *(m_xlines.end() - 1);
1190 mapymin = *m_ylines.begin();
1191 mapymax = *(m_ylines.end() - 1);
1192 mapzmin = *m_zlines.begin();
1193 mapzmax = *(m_zlines.end() - 1);
1194 mapvmin = *std::min_element(m_potential.begin(), m_potential.end());
1195 mapvmax = *std::max_element(m_potential.begin(), m_potential.end());
1196 // Set the periodicity length (maybe not needed).
1197 mapsx = fabs(mapxmax - mapxmin);
1198 mapsy = fabs(mapymax - mapymin);
1199 mapsz = fabs(mapzmax - mapzmin);
1200 // Set provisional cell dimensions.
1201 xMinBoundingBox = mapxmin;
1202 xMaxBoundingBox = mapxmax;
1203 yMinBoundingBox = mapymin;
1204 yMaxBoundingBox = mapymax;
1205 if (m_is3d) {
1206 zMinBoundingBox = mapzmin;
1207 zMaxBoundingBox = mapzmax;
1208 } else {
1209 mapzmin = zMinBoundingBox;
1210 mapzmax = zMaxBoundingBox;
1211 }
1212 hasBoundingBox = true;
1213}

Referenced by Initialise(), and ShiftComponent().

◆ SetRangeZ()

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

Definition at line 1215 of file ComponentCST.cc.

1215 {
1216
1217 if (fabs(zmax - zmin) <= 0.) {
1218 std::cerr << m_className << "::SetRangeZ:" << std::endl;
1219 std::cerr << " Zero range is not permitted." << std::endl;
1220 return;
1221 }
1222 zMinBoundingBox = std::min(zmin, zmax);
1223 zMaxBoundingBox = std::max(zmin, zmax);
1224}

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

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

◆ ShiftComponent()

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

Definition at line 977 of file ComponentCST.cc.

978 {
979 std::transform(m_xlines.begin(), m_xlines.end(), m_xlines.begin(),
980 std::bind1st(std::plus<double>(), xShift));
981 std::transform(m_ylines.begin(), m_ylines.end(), m_ylines.begin(),
982 std::bind1st(std::plus<double>(), yShift));
983 std::transform(m_zlines.begin(), m_zlines.end(), m_zlines.begin(),
984 std::bind1st(std::plus<double>(), zShift));
985 SetRange();
986 UpdatePeriodicity();
987
988 std::cout << m_className << "::ShiftComponent:" << std::endl;
989 std::cout << " Shifted component in x-direction: " << xShift
990 << "\t y-direction: " << yShift << "\t z-direction: " << zShift
991 << std::endl;
992}

◆ UpdatePeriodicity()

void Garfield::ComponentCST::UpdatePeriodicity ( )
protectedvirtual

Verify periodicities.

Implements Garfield::ComponentFieldMap.

Definition at line 1298 of file ComponentCST.cc.

1298 {
1299
1300 UpdatePeriodicity2d();
1301 UpdatePeriodicityCommon();
1302}

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 
)
virtual

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

Implements Garfield::ComponentFieldMap.

Definition at line 1008 of file ComponentCST.cc.

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

◆ WeightingPotential()

double Garfield::ComponentCST::WeightingPotential ( const double  x,
const double  y,
const double  z,
const std::string &  label 
)
virtual

Implements Garfield::ComponentFieldMap.

Definition at line 1072 of file ComponentCST.cc.

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