Garfield++ 4.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Public Member Functions | Protected Member Functions | List of all members
Garfield::ComponentCST Class Reference

#include <ComponentCST.hh>

+ Inheritance diagram for Garfield::ComponentCST:

Public Member Functions

 ComponentCST ()
 Constructor.
 
 ~ComponentCST ()
 Destructor.
 
void ShiftComponent (const double xShift, const double yShift, const double zShift)
 
void GetNumberOfMeshLines (unsigned int &nx, unsigned int &ny, unsigned int &nz) const
 
size_t GetNumberOfElements () const override
 Return the number of mesh elements.
 
bool GetElement (const size_t i, size_t &mat, bool &drift, std::vector< size_t > &nodes) const override
 Return the material and node indices of a mesh element.
 
size_t GetNumberOfNodes () const override
 
bool GetNode (const size_t i, double &x, double &y, double &z) const override
 
void GetElementBoundaries (unsigned int element, double &xmin, double &xmax, double &ymin, double &ymax, double &zmin, double &zmax) const
 
Medium * GetMedium (const double x, const double y, const double z) override
 Get the medium at a given location (x, y, z).
 
void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status) override
 
void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status) override
 Calculate the drift field [V/cm] and potential [V] at (x, y, z).
 
void WeightingField (const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label) override
 
double WeightingPotential (const double x, const double y, const double z, const std::string &label) override
 
bool Initialise (std::string elist, std::string nlist, std::string mplist, std::string prnsol, std::string unit="cm")
 
bool Initialise (std::string dataFile, std::string unit="cm")
 
bool SetWeightingField (std::string prnsol, std::string label, bool isBinary=true)
 
void 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) const
 
bool Coordinate2Index (const double x, const double y, const double z, unsigned int &i, unsigned int &j, unsigned int &k) const
 
- Public Member Functions inherited from Garfield::ComponentFieldMap
 ComponentFieldMap ()=delete
 Default constructor.
 
 ComponentFieldMap (const std::string &name)
 Constructor.
 
virtual ~ComponentFieldMap ()
 Destructor.
 
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 GetElementaryCell (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax) override
 Get the coordinates of the elementary cell.
 
bool GetVoltageRange (double &vmin, double &vmax) override
 Calculate the voltage range [V].
 
void 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.
 
size_t 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.
 
virtual size_t GetNumberOfElements () const
 Return the number of mesh elements.
 
bool GetElement (const size_t i, double &vol, double &dmin, double &dmax)
 Return the volume and aspect ratio of a mesh element.
 
virtual bool GetElement (const size_t i, size_t &mat, bool &drift, std::vector< size_t > &nodes) const
 Return the material and node indices of a mesh element.
 
virtual size_t GetNumberOfNodes () const
 
virtual bool GetNode (const size_t i, double &x, double &y, double &z) const
 
void EnableCheckMapIndices (const bool on=true)
 
void EnableDeleteBackgroundElements (const bool on=true)
 Option to eliminate mesh elements in conductors (default: on).
 
void EnableTetrahedralTreeForElementSearch (const bool on=true)
 
void EnableConvergenceWarnings (const bool on=true)
 
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::Component
 Component ()=delete
 Default constructor.
 
 Component (const std::string &name)
 Constructor.
 
virtual ~Component ()
 Destructor.
 
virtual void SetGeometry (Geometry *geo)
 Define the geometry.
 
virtual void Clear ()
 Reset.
 
virtual 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 double DelayedWeightingPotential (const double x, const double y, const double z, const double t, const std::string &label)
 
virtual void MagneticField (const double x, const double y, const double z, double &bx, double &by, double &bz, int &status)
 
void SetMagneticField (const double bx, const double by, const double bz)
 Set a constant magnetic field.
 
virtual bool IsReady ()
 Ready for use?
 
virtual bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
 Get the bounding box coordinates.
 
virtual bool GetElementaryCell (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
 Get the coordinates of the elementary cell.
 
double IntegrateFluxCircle (const double xc, const double yc, const double r, const unsigned int nI=50)
 
double IntegrateFluxSphere (const double xc, const double yc, const double zc, const double r, const unsigned int nI=20)
 
double IntegrateFluxParallelogram (const double x0, const double y0, const double z0, const double dx1, const double dy1, const double dz1, const double dx2, const double dy2, const double dz2, const unsigned int nU=20, const unsigned int nV=20)
 
double IntegrateWeightingFluxParallelogram (const std::string &label, const double x0, const double y0, const double z0, const double dx1, const double dy1, const double dz1, const double dx2, const double dy2, const double dz2, const unsigned int nU=20, const unsigned int nV=20)
 
double IntegrateFluxLine (const double x0, const double y0, const double z0, const double x1, const double y1, const double z1, const double xp, const double yp, const double zp, const unsigned int nI, const int isign=0)
 
virtual bool IsWireCrossed (const double x0, const double y0, const double z0, const double x1, const double y1, const double z1, double &xc, double &yc, double &zc, const bool centre, double &rc)
 
virtual bool IsInTrapRadius (const double q0, const double x0, const double y0, const double z0, double &xw, double &yw, double &rw)
 
void EnablePeriodicityX (const bool on=true)
 Enable simple periodicity in the $x$ direction.
 
void EnablePeriodicityY (const bool on=true)
 Enable simple periodicity in the $y$ direction.
 
void EnablePeriodicityZ (const bool on=true)
 Enable simple periodicity in the $z$ direction.
 
void IsPeriodic (bool &perx, bool &pery, bool &perz)
 Return periodicity flags.
 
void EnableMirrorPeriodicityX (const bool on=true)
 Enable mirror periodicity in the $x$ direction.
 
void EnableMirrorPeriodicityY (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void EnableMirrorPeriodicityZ (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void IsMirrorPeriodic (bool &perx, bool &pery, bool &perz)
 Return mirror periodicity flags.
 
void EnableAxialPeriodicityX (const bool on=true)
 Enable axial periodicity in the $x$ direction.
 
void EnableAxialPeriodicityY (const bool on=true)
 Enable axial periodicity in the $y$ direction.
 
void EnableAxialPeriodicityZ (const bool on=true)
 Enable axial periodicity in the $z$ direction.
 
void IsAxiallyPeriodic (bool &perx, bool &pery, bool &perz)
 Return axial periodicity flags.
 
void EnableRotationSymmetryX (const bool on=true)
 Enable rotation symmetry around the $x$ axis.
 
void EnableRotationSymmetryY (const bool on=true)
 Enable rotation symmetry around the $y$ axis.
 
void EnableRotationSymmetryZ (const bool on=true)
 Enable rotation symmetry around the $z$ axis.
 
void IsRotationSymmetric (bool &rotx, bool &roty, bool &rotz)
 Return rotation symmetry flags.
 
void EnableDebugging ()
 Switch on debugging messages.
 
void DisableDebugging ()
 Switch off debugging messages.
 
virtual bool HasAttachmentMap () const
 Does the component have attachment maps?
 
virtual bool HasVelocityMap () const
 Does the component have velocity maps?
 
virtual bool ElectronAttachment (const double, const double, const double, double &eta)
 Get the electron attachment coefficient.
 
virtual bool HoleAttachment (const double, const double, const double, double &eta)
 Get the hole attachment coefficient.
 
virtual bool ElectronVelocity (const double, const double, const double, double &vx, double &vy, double &vz)
 Get the electron drift velocity.
 
virtual bool HoleVelocity (const double, const double, const double, double &vx, double &vy, double &vz)
 Get the hole drift velocity.
 
virtual bool GetElectronLifetime (const double, const double, const double, double &etau)
 
virtual bool GetHoleLifetime (const double, const double, const double, double &htau)
 

Protected Member Functions

void UpdatePeriodicity () override
 Verify periodicities.
 
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) const
 
- Protected Member Functions inherited from Garfield::ComponentFieldMap
void Reset () override
 Reset the component.
 
void Prepare ()
 
void UpdatePeriodicity2d ()
 
void UpdatePeriodicityCommon ()
 
bool SetDefaultDriftMedium ()
 Find lowest epsilon, check for eps = 0, set default drift media flags.
 
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
 
size_t GetWeightingFieldIndex (const std::string &label) const
 
size_t GetOrCreateWeightingFieldIndex (const std::string &label)
 
void PrintWarning (const std::string &header)
 
void PrintNotReady (const std::string &header) const
 
void PrintCouldNotOpen (const std::string &header, const std::string &filename) const
 
void PrintElement (const std::string &header, const double x, const double y, const double z, const double t1, const double t2, const double t3, const double t4, const 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

- Static Protected Member Functions inherited from Garfield::ComponentFieldMap
static double ScalingFactor (std::string unit)
 
- Protected Attributes inherited from Garfield::ComponentFieldMap
bool m_is3d = true
 
std::vector< Element > m_elements
 
std::vector< Node > m_nodes
 
std::vector< Material > m_materials
 
std::vector< std::string > m_wfields
 
std::vector< bool > m_wfieldsOk
 
std::vector< bool > m_dwfieldsOk
 
std::vector< double > m_wdtimes
 
bool m_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 = 0.
 
double m_mapvmax = 0.
 
std::array< bool, 3 > m_setang
 
bool m_deleteBackground = true
 
bool m_warning = false
 
unsigned int m_nWarnings = 0
 
bool m_printConvergenceWarnings = true
 
- Protected Attributes inherited from Garfield::Component
std::string m_className = "Component"
 Class name.
 
Geometry * m_geometry = nullptr
 Pointer to the geometry.
 
std::array< double, 3 > m_b0 = {{0., 0., 0.}}
 Constant magnetic field.
 
bool m_ready = false
 Ready for use?
 
bool m_debug = false
 Switch on/off debugging messages.
 
std::array< bool, 3 > m_periodic = {{false, false, false}}
 Simple periodicity in x, y, z.
 
std::array< bool, 3 > m_mirrorPeriodic = {{false, false, false}}
 Mirror periodicity in x, y, z.
 
std::array< bool, 3 > m_axiallyPeriodic = {{false, false, false}}
 Axial periodicity in x, y, z.
 
std::array< bool, 3 > m_rotationSymmetric = {{false, false, false}}
 Rotation symmetry around x-axis, y-axis, z-axis.
 

Detailed Description

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

Definition at line 15 of file ComponentCST.hh.

Constructor & Destructor Documentation

◆ ComponentCST()

Garfield::ComponentCST::ComponentCST ( )

Constructor.

Definition at line 80 of file ComponentCST.cc.

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

◆ ~ComponentCST()

Garfield::ComponentCST::~ComponentCST ( )
inline

Destructor.

Definition at line 20 of file ComponentCST.hh.

20{}

Member Function Documentation

◆ Coordinate2Index() [1/2]

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

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

Definition at line 1084 of file ComponentCST.cc.

1086 {
1087 bool mirrored[3] = {false, false, false};
1088 double pos[3] = {0., 0., 0.};
1089 return Coordinate2Index(x, y, z, i, j, k, pos, mirrored);
1090}
Garfield::ComponentCST::Coordinate2Index
bool Coordinate2Index(const double x, const double y, const double z, unsigned int &i, unsigned int &j, unsigned int &k) const
Definition: ComponentCST.cc:1084

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

◆ Coordinate2Index() [2/2]

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

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

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

Definition at line 1093 of file ComponentCST.cc.

1096 {
1097 // Map the coordinates onto field map coordinates
1098 pos[0] = xin;
1099 pos[1] = yin;
1100 pos[2] = zin;
1101 double rcoordinate = 0.;
1102 double rotation = 0.;
1103 MapCoordinates(pos[0], pos[1], pos[2],
1104 mirrored[0], mirrored[1], mirrored[2], rcoordinate, rotation);
1105
1106 auto it_x =
1107 std::lower_bound(m_xlines.begin(), m_xlines.end(), pos[0]);
1108 auto it_y =
1109 std::lower_bound(m_ylines.begin(), m_ylines.end(), pos[1]);
1110 auto it_z =
1111 std::lower_bound(m_zlines.begin(), m_zlines.end(), pos[2]);
1112 if (it_x == m_xlines.end() || it_y == m_ylines.end() ||
1113 it_z == m_zlines.end() || pos[0] < m_xlines.at(0) ||
1114 pos[1] < m_ylines.at(0) ||
1115 pos[2] < m_zlines.at(0)) {
1116 if (m_debug) {
1117 std::cerr << m_className << "::ElectricFieldBinary:" << std::endl;
1118 std::cerr << " Could not find the given coordinate!" << std::endl;
1119 std::cerr << " You ask for the following position: " << xin << ", "
1120 << yin << ", " << zin << std::endl;
1121 std::cerr << " The mapped position is: " << pos[0] << ", "
1122 << pos[1] << ", " << pos[2]
1123 << std::endl;
1124 }
1125 return false;
1126 }
1127 /* Lower bound returns the next mesh line behind the position in question.
1128 * If the position in question is on a mesh line this mesh line is returned.
1129 * Since we are interested in the mesh line before the position in question we
1130 * need to move the
1131 * iterator to the left except for the very first mesh line!
1132 */
1133 if (it_x == m_xlines.begin())
1134 i = 0;
1135 else
1136 i = std::distance(m_xlines.begin(), it_x - 1);
1137 if (it_y == m_ylines.begin())
1138 j = 0;
1139 else
1140 j = std::distance(m_ylines.begin(), it_y - 1);
1141 if (it_z == m_zlines.begin())
1142 k = 0;
1143 else
1144 k = std::distance(m_zlines.begin(), it_z - 1);
1145 return true;
1146}
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:2021
Garfield::Component::m_debug
bool m_debug
Switch on/off debugging messages.
Definition: Component.hh:341
Garfield::Component::m_className
std::string m_className
Class name.
Definition: Component.hh:329

◆ DisableShaping()

void Garfield::ComponentCST::DisableShaping ( )
inline

Definition at line 181 of file ComponentCST.hh.

181{ doShaping = false; }

◆ DisableXField()

void Garfield::ComponentCST::DisableXField ( )
inline

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

Definition at line 134 of file ComponentCST.hh.

134{ disableFieldComponent[0] = true; }

◆ DisableYField()

void Garfield::ComponentCST::DisableYField ( )
inline

Definition at line 135 of file ComponentCST.hh.

135{ disableFieldComponent[1] = true; }

◆ DisableZField()

void Garfield::ComponentCST::DisableZField ( )
inline

Definition at line 136 of file ComponentCST.hh.

136{ disableFieldComponent[2] = true; }

◆ ElectricField() [1/2]

void Garfield::ComponentCST::ElectricField ( const double  x,
const double  y,
const double  z,
double &  ex,
double &  ey,
double &  ez,
double &  v,
Medium *&  m,
int &  status 
)
overridevirtual

Calculate the drift field [V/cm] and potential [V] at (x, y, z).

Implements Garfield::Component.

Definition at line 835 of file ComponentCST.cc.

838 {
839 ElectricFieldBinary(xin, yin, zin, ex, ey, ez, volt, m, status, true);
840}

◆ ElectricField() [2/2]

void Garfield::ComponentCST::ElectricField ( const double  x,
const double  y,
const double  z,
double &  ex,
double &  ey,
double &  ez,
Medium *&  m,
int &  status 
)
overridevirtual

Calculate the drift field at given point.

Parameters
x,y,zcoordinates [cm].
ex,ey,ezcomponents of the electric field [V/cm].
mpointer to the medium at this location.
statusstatus flag

Status flags: 

        0: Inside an active medium
      > 0: Inside a wire of type X
-4 ... -1: On the side of a plane where no wires are
       -5: Inside the mesh but not in an active medium
       -6: Outside the mesh
      -10: Unknown potential type (should not occur)
    other: Other cases (should not occur)

Implements Garfield::Component.

Definition at line 828 of file ComponentCST.cc.

830 {
831 double volt;
832 ElectricFieldBinary(xin, yin, zin, ex, ey, ez, volt, m, status);
833}

◆ EnableShaping()

void Garfield::ComponentCST::EnableShaping ( )
inline

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

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

__________ ____________

element 1 element 2

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

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

Definition at line 180 of file ComponentCST.hh.

180{ doShaping = true; }

◆ GetAspectRatio()

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

Implements Garfield::ComponentFieldMap.

Definition at line 1161 of file ComponentCST.cc.

1162 {
1163 if (element >= m_nElements) {
1164 dmin = dmax = 0.;
1165 return;
1166 }
1167 unsigned int i, j, k;
1168 Element2Index(element, i, j, k);
1169 const double dx = fabs(m_xlines.at(i + 1) - m_xlines.at(i));
1170 const double dy = fabs(m_ylines.at(j + 1) - m_ylines.at(j));
1171 const double dz = fabs(m_zlines.at(k + 1) - m_zlines.at(k));
1172 dmin = std::min({dx, dy, dz});
1173 dmax = std::max({dx, dy, dz});
1174}
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615

◆ GetElement()

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

Return the material and node indices of a mesh element.

Reimplemented from Garfield::ComponentFieldMap.

Definition at line 977 of file ComponentCST.cc.

978 {
979 if (element >= m_nElements || element >= m_elementMaterial.size()) {
980 std::cerr << m_className << "::GetElement: Index out of range.\n";
981 return false;
982 }
983 mat = m_elementMaterial[element];
984 drift = m_materials[mat].driftmedium;
985 nodes.clear();
986 unsigned int i0 = 0, j0 = 0, k0 = 0;
987 Element2Index(element, i0, j0, k0);
988 const auto i1 = i0 + 1;
989 const auto j1 = j0 + 1;
990 const auto k1 = k0 + 1;
991 nodes.push_back(Index2Node(i0, j0, k0));
992 nodes.push_back(Index2Node(i1, j0, k0));
993 nodes.push_back(Index2Node(i0, j1, k0));
994 nodes.push_back(Index2Node(i1, j1, k0));
995 nodes.push_back(Index2Node(i0, j0, k1));
996 nodes.push_back(Index2Node(i1, j0, k1));
997 nodes.push_back(Index2Node(i0, j1, k1));
998 nodes.push_back(Index2Node(i1, j1, k1));
999 return true;
1000}
Garfield::ComponentFieldMap::m_materials
std::vector< Material > m_materials
Definition: ComponentFieldMap.hh:138

◆ GetElementBoundaries()

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

Definition at line 1016 of file ComponentCST.cc.

1019 {
1020 unsigned int i, j, k;
1021 Element2Index(element, i, j, k);
1022 xmin = m_xlines.at(i);
1023 xmax = m_xlines.at(i + 1);
1024 ymin = m_ylines.at(j);
1025 ymax = m_ylines.at(j + 1);
1026 zmin = m_zlines.at(k);
1027 zmax = m_zlines.at(k + 1);
1028}

◆ GetElementVolume()

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

Implements Garfield::ComponentFieldMap.

Definition at line 1176 of file ComponentCST.cc.

1176 {
1177 if (element >= m_nElements) return 0.;
1178 unsigned int i, j, k;
1179 Element2Index(element, i, j, k);
1180 const double dx = fabs(m_xlines.at(i + 1) - m_xlines.at(i));
1181 const double dy = fabs(m_ylines.at(j + 1) - m_ylines.at(j));
1182 const double dz = fabs(m_zlines.at(k + 1) - m_zlines.at(k));
1183 return dx * dy * dz;
1184}

◆ GetMedium()

Medium * Garfield::ComponentCST::GetMedium ( const double  x,
const double  y,
const double  z 
)
overridevirtual

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

Reimplemented from Garfield::Component.

Definition at line 1030 of file ComponentCST.cc.

1031 {
1032 unsigned int i, j, k;
1033 Coordinate2Index(x, y, z, i, j, k);
1034 if (m_debug) {
1035 std::cout << m_className << "::GetMedium:\n"
1036 << " Position (" << x << ", " << y << ", " << z << "):\n"
1037 << " Indices are: x: " << i << "/" << m_xlines.size()
1038 << "\t y: " << j << "/" << m_ylines.size()
1039 << "\t z: " << k << "/" << m_zlines.size() << std::endl;
1040 const auto element = Index2Element(i, j, k);
1041 std::cout << " Element index: " << element << std::endl
1042 << " Material index: "
1043 << (int)m_elementMaterial.at(element) << std::endl;
1044 }
1045 return m_materials.at(m_elementMaterial.at(Index2Element(i, j, k))).medium;
1046}
Garfield::ComponentCST::Index2Element
int Index2Element(const unsigned int i, const unsigned int j, const unsigned int k) const
Definition: ComponentCST.cc:1148

◆ GetNode()

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

Reimplemented from Garfield::ComponentFieldMap.

Definition at line 1002 of file ComponentCST.cc.

1003 {
1004 if (node >= m_nNodes) {
1005 std::cerr << m_className << "::GetNode: Index out of range.\n";
1006 return false;
1007 }
1008 unsigned int i = 0, j = 0, k = 0;
1009 Node2Index(node, i, j, k);
1010 x = m_xlines[i];
1011 y = m_ylines[j];
1012 z = m_zlines[k];
1013 return true;
1014}

◆ GetNumberOfElements()

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

Return the number of mesh elements.

Reimplemented from Garfield::ComponentFieldMap.

Definition at line 27 of file ComponentCST.hh.

27{ return m_nElements; }

◆ GetNumberOfMeshLines()

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

Definition at line 970 of file ComponentCST.cc.

971 {
972 n_x = m_xlines.size();
973 n_y = m_ylines.size();
974 n_z = m_zlines.size();
975}

◆ GetNumberOfNodes()

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

Reimplemented from Garfield::ComponentFieldMap.

Definition at line 30 of file ComponentCST.hh.

30{ return m_nNodes; }

◆ Index2Element()

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

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

Definition at line 1148 of file ComponentCST.cc.

1149 {
1150 if (i > m_nx - 2 || j > m_ny - 2 || k > m_nz - 2) {
1151 throw "ComponentCST::Index2Element: Error. Element indices out of bounds.";
1152 }
1153 return i + j * (m_nx - 1) + k * (m_nx - 1) * (m_ny - 1);
1154}

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

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

◆ Initialise() [2/2]

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

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

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

Definition at line 87 of file ComponentCST.cc.

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

◆ SetRange()

void Garfield::ComponentCST::SetRange ( )
overridevirtual

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

Reimplemented from Garfield::ComponentFieldMap.

Definition at line 1048 of file ComponentCST.cc.

1048 {
1049 // Establish the ranges
1050 m_mapmin[0] = m_xlines.front();
1051 m_mapmax[0] = m_xlines.back();
1052 m_mapmin[1] = m_ylines.front();
1053 m_mapmax[1] = m_ylines.back();
1054 m_mapmin[2] = m_zlines.front();
1055 m_mapmax[2] = m_zlines.back();
1056 m_mapvmin = *std::min_element(m_potential.begin(), m_potential.end());
1057 m_mapvmax = *std::max_element(m_potential.begin(), m_potential.end());
1058
1059 // Set provisional cell dimensions.
1060 m_minBoundingBox[0] = m_mapmin[0];
1061 m_maxBoundingBox[0] = m_mapmax[0];
1062 m_minBoundingBox[1] = m_mapmin[1];
1063 m_maxBoundingBox[1] = m_mapmax[1];
1064 if (m_is3d) {
1065 m_minBoundingBox[2] = m_mapmin[2];
1066 m_maxBoundingBox[2] = m_mapmax[2];
1067 } else {
1068 m_mapmin[2] = m_minBoundingBox[2];
1069 m_mapmax[2] = m_maxBoundingBox[2];
1070 }
1071 m_hasBoundingBox = true;
1072}
Garfield::ComponentFieldMap::m_mapmin
std::array< double, 3 > m_mapmin
Definition: ComponentFieldMap.hh:152
Garfield::ComponentFieldMap::m_mapmax
std::array< double, 3 > m_mapmax
Definition: ComponentFieldMap.hh:153

Referenced by Initialise(), and ShiftComponent().

◆ SetRangeZ()

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

Definition at line 1074 of file ComponentCST.cc.

1074 {
1075 if (fabs(zmax - zmin) <= 0.) {
1076 std::cerr << m_className << "::SetRangeZ:" << std::endl;
1077 std::cerr << " Zero range is not permitted." << std::endl;
1078 return;
1079 }
1080 m_minBoundingBox[2] = std::min(zmin, zmax);
1081 m_maxBoundingBox[2] = std::max(zmin, zmax);
1082}

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

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

◆ ShiftComponent()

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

Definition at line 814 of file ComponentCST.cc.

815 {
816 std::transform(m_xlines.begin(), m_xlines.end(), m_xlines.begin(), [xShift](double x) { return x + xShift;});
817 std::transform(m_ylines.begin(), m_ylines.end(), m_ylines.begin(), [yShift](double x) { return x + yShift;});
818 std::transform(m_zlines.begin(), m_zlines.end(), m_zlines.begin(), [zShift](double x) { return x + zShift;});
819 SetRange();
820 UpdatePeriodicity();
821
822 std::cout << m_className << "::ShiftComponent:" << std::endl;
823 std::cout << " Shifted component in x-direction: " << xShift
824 << "\t y-direction: " << yShift << "\t z-direction: " << zShift
825 << std::endl;
826}

◆ UpdatePeriodicity()

void Garfield::ComponentCST::UpdatePeriodicity ( )
overrideprotectedvirtual

Verify periodicities.

Implements Garfield::Component.

Definition at line 1156 of file ComponentCST.cc.

1156 {
1157 UpdatePeriodicity2d();
1158 UpdatePeriodicityCommon();
1159}

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

Definition at line 842 of file ComponentCST.cc.

844 {
845 // Initial values
846 wx = wy = wz = 0;
847
848 // Do not proceed if not properly initialised.
849 if (!m_ready) return;
850
851 // Look for the label.
852 auto it = m_weightingFields.find(label);
853 if (it == m_weightingFields.end()) {
854 // Do not proceed if the requested weighting field does not exist.
855 std::cerr << "No weighting field named " << label.c_str() << " found!"
856 << std::endl;
857 return;
858 }
859
860 // Check if the weighting field is properly initialised.
861 if (!m_wfieldsOk[std::distance(m_weightingFields.begin(), it)]) return;
862
863 // Copy the coordinates
864 double x = xin, y = yin, z = zin;
865
866 // Map the coordinates onto field map coordinates and get indexes
867 bool mirrored[3];
868 unsigned int i, j, k;
869 double pos[3] = {0., 0., 0.};
870 if (!Coordinate2Index(x, y, z, i, j, k, pos, mirrored)) {
871 return;
872 }
873
874 double rx = (pos[0] - m_xlines.at(i)) / (m_xlines.at(i + 1) - m_xlines.at(i));
875 double ry = (pos[1] - m_ylines.at(j)) / (m_ylines.at(j + 1) - m_ylines.at(j));
876 double rz = (pos[2] - m_zlines.at(k)) / (m_zlines.at(k + 1) - m_zlines.at(k));
877
878 float fwx = 0., fwy = 0., fwz = 0.;
879 if (!disableFieldComponent[0])
880 fwx = GetFieldComponent(i, j, k, rx, ry, rz, 'x', (*it).second);
881 if (!disableFieldComponent[1])
882 fwy = GetFieldComponent(i, j, k, rx, ry, rz, 'y', (*it).second);
883 if (!disableFieldComponent[2])
884 fwz = GetFieldComponent(i, j, k, rx, ry, rz, 'z', (*it).second);
885
886 if (m_elementMaterial.size() > 0 && doShaping) {
887 ShapeField(fwx, fwy, fwz, rx, ry, rz, i, j, k, (*it).second);
888 }
889 if (mirrored[0]) fwx *= -1.f;
890 if (mirrored[1]) fwy *= -1.f;
891 if (mirrored[2]) fwz *= -1.f;
892 if (m_warning) PrintWarning("WeightingField");
893 if (m_materials.at(m_elementMaterial.at(Index2Element(i, j, k))).driftmedium) {
894 if (!disableFieldComponent[0]) wx = fwx;
895 if (!disableFieldComponent[1]) wy = fwy;
896 if (!disableFieldComponent[2]) wz = fwz;
897 }
898}
Garfield::ComponentFieldMap::PrintWarning
void PrintWarning(const std::string &header)
Definition: ComponentFieldMap.cc:2375

◆ WeightingPotential()

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

Calculate the weighting potential at a given point.

Parameters
x,y,zcoordinates [cm].
labelname of the electrode.
Returns
weighting potential [dimensionless].

Reimplemented from Garfield::Component.

Definition at line 900 of file ComponentCST.cc.

902 {
903 // Do not proceed if not properly initialised.
904 if (!m_ready) return 0.;
905
906 // Look for the label.
907 auto it = m_weightingFields.find(label);
908 if (it == m_weightingFields.end()) {
909 // Do not proceed if the requested weighting field does not exist.
910 std::cerr << "No weighting field named " << label.c_str() << " found!"
911 << std::endl;
912 return 0.;
913 }
914
915 // Check if the weighting field is properly initialised.
916 if (!m_wfieldsOk[std::distance(m_weightingFields.begin(), it)]) return 0.;
917
918 // Copy the coordinates
919 double x = xin, y = yin, z = zin;
920
921 // Map the coordinates onto field map coordinates
922 bool mirrored[3];
923 unsigned int i, j, k;
924 double pos[3] = {0., 0., 0.};
925 if (!Coordinate2Index(x, y, z, i, j, k, pos, mirrored)) {
926 return 0.;
927 }
928 double rx = (pos[0] - m_xlines.at(i)) /
929 (m_xlines.at(i + 1) - m_xlines.at(i));
930 double ry = (pos[1] - m_ylines.at(j)) /
931 (m_ylines.at(j + 1) - m_ylines.at(j));
932 double rz = (pos[2] - m_zlines.at(k)) /
933 (m_zlines.at(k + 1) - m_zlines.at(k));
934
935 double potential = GetPotential(i, j, k, rx, ry, rz, (*it).second);
936
937 if (m_debug) {
938 std::cout << m_className << "::WeightingPotential:" << std::endl;
939 std::cout << " Global: (" << x << "," << y << "," << z << "),"
940 << std::endl;
941 std::cout << " Local: (" << rx << "," << ry << "," << rz
942 << ") in element with indexes: i=" << i << ", j=" << j
943 << ", k=" << k << std::endl;
944 std::cout << " Node xyzV:" << std::endl;
945 std::cout << "Node 0 position: " << Index2Node(i + 1, j, k)
946 << "\t potential: " << ((*it).second).at(Index2Node(i + 1, j, k))
947 << "Node 1 position: " << Index2Node(i + 1, j + 1, k)
948 << "\t potential: "
949 << ((*it).second).at(Index2Node(i + 1, j + 1, k))
950 << "Node 2 position: " << Index2Node(i, j + 1, k)
951 << "\t potential: " << ((*it).second).at(Index2Node(i, j + 1, k))
952 << "Node 3 position: " << Index2Node(i, j, k)
953 << "\t potential: " << ((*it).second).at(Index2Node(i, j, k))
954 << "Node 4 position: " << Index2Node(i + 1, j, k + 1)
955 << "\t potential: "
956 << ((*it).second).at(Index2Node(i + 1, j, k + 1))
957 << "Node 5 position: " << Index2Node(i + 1, j + 1, k + 1)
958 << "\t potential: "
959 << ((*it).second).at(Index2Node(i + 1, j + 1, k + 1))
960 << "Node 6 position: " << Index2Node(i, j + 1, k + 1)
961 << "\t potential: "
962 << ((*it).second).at(Index2Node(i, j + 1, k + 1))
963 << "Node 7 position: " << Index2Node(i, j, k + 1)
964 << "\t potential: " << ((*it).second).at(Index2Node(i, j, k))
965 << std::endl;
966 }
967 return potential;
968}
The documentation for this class was generated from the following files: