Garfield++ v1r0
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)
 
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)
 
bool IsInBoundingBox (const double x, const double y, const double z)
 
void SetRange ()
 
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 ()
 
virtual ~ComponentFieldMap ()
 
virtual void SetRange ()
 
void PrintRange ()
 
virtual bool IsInBoundingBox (const double x, const double y, const double z)
 
virtual bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
 
bool GetVoltageRange (double &vmin, double &vmax)
 
void PrintMaterials ()
 
void DriftMedium (int imat)
 
void NotDriftMedium (int imat)
 
int GetNumberOfMaterials ()
 
double GetPermittivity (const int imat)
 
double GetConductivity (const int imat)
 
void SetMedium (const int imat, Medium *medium)
 
Medium * GetMedium (const unsigned int &i) const
 
Medium * GetMedium (const double &x, const double &y, const double &z)=0
 
int GetNumberOfMedia ()
 
int GetNumberOfElements () const
 
bool GetElement (const int i, double &vol, double &dmin, double &dmax)
 
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 ()
 
- Public Member Functions inherited from Garfield::ComponentBase
 ComponentBase ()
 
virtual ~ComponentBase ()
 
virtual void SetGeometry (GeometryBase *geo)
 
virtual void Clear ()
 
virtual Medium * GetMedium (const double &x, const double &y, const double &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
 
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)
 
virtual bool IsReady ()
 
virtual bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
 
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 (double x0, double y0, double z0, double &xw, double &yw, double &rw)
 
void EnablePeriodicityX ()
 
void DisablePeriodicityX ()
 
void EnablePeriodicityY ()
 
void DisablePeriodicityY ()
 
void EnablePeriodicityZ ()
 
void DisablePeriodicityZ ()
 
void EnableMirrorPeriodicityX ()
 
void DisableMirrorPeriodicityX ()
 
void EnableMirrorPeriodicityY ()
 
void DisableMirrorPeriodicityY ()
 
void EnableMirrorPeriodicityZ ()
 
void DisableMirrorPeriodicityZ ()
 
void EnableAxialPeriodicityX ()
 
void DisableAxialPeriodicityX ()
 
void EnableAxialPeriodicityY ()
 
void DisableAxialPeriodicityY ()
 
void EnableAxialPeriodicityZ ()
 
void DisableAxialPeriodicityZ ()
 
void EnableRotationSymmetryX ()
 
void DisableRotationSymmetryX ()
 
void EnableRotationSymmetryY ()
 
void DisableRotationSymmetryY ()
 
void EnableRotationSymmetryZ ()
 
void DisableRotationSymmetryZ ()
 
void EnableDebugging ()
 
void DisableDebugging ()
 

Protected Member Functions

void UpdatePeriodicity ()
 
double GetElementVolume (const int i)
 
void GetAspectRatio (const 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 ()
 
virtual void UpdatePeriodicity ()=0
 
void UpdatePeriodicity2d ()
 
void UpdatePeriodicityCommon ()
 
int Coordinates3 (double x, double y, double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det, int imap)
 
int Coordinates4 (double x, double y, double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det, int imap)
 
int Coordinates5 (double x, double y, double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det, int imap)
 
int Coordinates12 (double x, double y, double z, double &t1, double &t2, double &t3, double &t4, int imap)
 
int Coordinates13 (double x, double y, double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det, int imap)
 
int CoordinatesCube (double x, double y, double z, double &t1, double &t2, double &t3, TMatrixD *&jac, std::vector< TMatrixD * > &dN, int imap)
 
void Jacobian3 (int i, double u, double v, double w, double &det, double jac[4][4])
 
void Jacobian5 (int i, double u, double v, double &det, double jac[4][4])
 
void Jacobian13 (int i, double t, double u, double v, double w, double &det, double jac[4][4])
 
void JacobianCube (int i, double t1, double t2, double t3, TMatrixD *&jac, std::vector< TMatrixD * > &dN)
 
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)
 
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)
 
int FindElementCube (const double x, const double y, const double z, double &t1, double &t2, double &t3, TMatrixD *&jac, std::vector< TMatrixD * > &dN)
 
void MapCoordinates (double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
 
void UnmapFields (double &ex, double &ey, double &ez, double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation)
 
int ReadInteger (char *token, int def, bool &error)
 
double ReadDouble (char *token, double def, bool &error)
 
virtual double GetElementVolume (const int i)=0
 
virtual void GetAspectRatio (const int i, double &dmin, double &dmax)=0
 
void CalculateElementBoundingBoxes (void)
 
virtual void Reset ()=0
 
virtual void UpdatePeriodicity ()=0
 

Additional Inherited Members

- Protected Attributes inherited from Garfield::ComponentFieldMap
bool is3d
 
int nElements
 
std::vector< element > elements
 
int lastElement
 
bool cacheElemBoundingBoxes
 
int nNodes
 
std::vector< node > nodes
 
int 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 deleteBackground
 
bool checkMultipleElement
 
bool warning
 
- Protected Attributes inherited from Garfield::ComponentBase
std::string m_className
 
GeometryBase * theGeometry
 
bool ready
 
bool xPeriodic
 
bool yPeriodic
 
bool zPeriodic
 
bool xMirrorPeriodic
 
bool yMirrorPeriodic
 
bool zMirrorPeriodic
 
bool xAxiallyPeriodic
 
bool yAxiallyPeriodic
 
bool zAxiallyPeriodic
 
bool xRotationSymmetry
 
bool yRotationSymmetry
 
bool zRotationSymmetry
 
double bx0
 
double by0
 
double bz0
 
bool debug
 

Detailed Description

Definition at line 15 of file ComponentCST.hh.

Constructor & Destructor Documentation

◆ ComponentCST()

Garfield::ComponentCST::ComponentCST ( )

Definition at line 20 of file ComponentCST.cc.

20 : ComponentFieldMap() {
21
22 m_className = "ComponentCST";
23 ready = false;
24 // Default bounding box
25 zMinBoundingBox = -50.;
26 zMaxBoundingBox = 50.;
27 m_xlines.clear();
28 m_ylines.clear();
29 m_zlines.clear();
30 deleteBackground = false;
31 disableFieldComponent[0] = false;
32 disableFieldComponent[1] = false;
33 disableFieldComponent[2] = false;
34 doShaping = false;
35}

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

1092 {
1093 bool mirrored[3] = {false, false, false};
1094 double position_mapped[3] = {0., 0., 0.};
1095 double rcoordinate, rotation;
1096 return Coordinate2Index(x ,y, z, i, j, k, position_mapped, mirrored, rcoordinate, rotation);
1097}
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:1091

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

1111 {
1112 // Map the coordinates onto field map coordinates
1113 position_mapped[0] = xin;
1114 position_mapped[1] = yin;
1115 position_mapped[2] = zin;
1116 MapCoordinates(position_mapped[0], position_mapped[1], position_mapped[2],
1117 mirrored[0], mirrored[1], mirrored[2], rcoordinate, rotation);
1118
1119 std::vector<double>::iterator it_x, it_y, it_z;
1120 it_x = std::lower_bound(m_xlines.begin(),m_xlines.end(),position_mapped[0]);
1121 it_y = std::lower_bound(m_ylines.begin(),m_ylines.end(),position_mapped[1]);
1122 it_z = std::lower_bound(m_zlines.begin(),m_zlines.end(),position_mapped[2]);
1123 if(it_x == m_xlines.end() || it_y == m_ylines.end() || it_z == m_zlines.end() ||
1124 position_mapped[0] < m_xlines.at(0) || position_mapped[1] < m_ylines.at(0) || position_mapped[2] < m_zlines.at(0) ){
1125 if(debug){
1126 std::cerr << m_className << "::ElectricFieldBinary:" << std::endl;
1127 std::cerr << " Could not find the given coordinate!" << std::endl;
1128 std::cerr << " You ask for the following position: " << xin << ", " << yin << ", " << zin << std::endl;
1129 std::cerr << " The mapped position is: " << position_mapped[0] << ", " << position_mapped[1] << ", " << position_mapped[2] << std::endl;
1130 PrintRange();
1131 }
1132 return false;
1133 }
1134 /* Lower bound returns the next mesh line behind the position in question.
1135 * If the position in question is on a mesh line this mesh line is returned.
1136 * Since we are interested in the mesh line before the position in question we need to move the
1137 * iterator to the left except for the very first mesh line!
1138 */
1139 if(it_x == m_xlines.begin())
1140 i = 0;
1141 else
1142 i = std::distance(m_xlines.begin(),it_x-1);
1143 if(it_y == m_ylines.begin())
1144 j = 0;
1145 else
1146 j = std::distance(m_ylines.begin(),it_y-1);
1147 if(it_z == m_zlines.begin())
1148 k = 0;
1149 else
1150 k = std::distance(m_zlines.begin(),it_z-1);
1151 return true;
1152}
Garfield::ComponentFieldMap::MapCoordinates
void MapCoordinates(double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
Definition: ComponentFieldMap.cc:4036

◆ DisableShaping()

void Garfield::ComponentCST::DisableShaping ( )
inline

Definition at line 166 of file ComponentCST.hh.

166 {
167 doShaping = false;
168 }

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

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

◆ DisableYField()

void Garfield::ComponentCST::DisableYField ( )
inline

Definition at line 128 of file ComponentCST.hh.

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

◆ DisableZField()

void Garfield::ComponentCST::DisableZField ( )
inline

Definition at line 131 of file ComponentCST.hh.

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

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

891 {
892 ElectricFieldBinary(xin, yin, zin, ex, ey, ez, volt, m, status, true);
893}

◆ 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

Implements Garfield::ComponentFieldMap.

Definition at line 882 of file ComponentCST.cc.

884 {
885 double volt;
886 ElectricFieldBinary(xin, yin, zin, ex, ey, ez, volt, m, status);
887}

◆ 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 \f[E_1\] along the \f[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 163 of file ComponentCST.hh.

163 {
164 doShaping = true;
165 }

◆ GetAspectRatio()

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

Implements Garfield::ComponentFieldMap.

Definition at line 1160 of file ComponentCST.cc.

1160 {
1161
1162 if (element < 0 || element >= nElements) {
1163 dmin = dmax = 0.;
1164 return;
1165 }
1166 unsigned int i, j, k;
1167 Element2Index(element, i, j, k);
1168 std::vector<double> distances;
1169 distances.push_back(m_xlines.at(i+1) - m_xlines.at(i));
1170 distances.push_back(m_ylines.at(j+1) - m_ylines.at(j));
1171 distances.push_back(m_zlines.at(k+1) - m_zlines.at(k));
1172 std::sort(distances.begin(), distances.end());
1173 dmin = distances.at(0);
1174 dmax = distances.at(2);
1175}

◆ GetElementBoundaries()

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

Definition at line 1022 of file ComponentCST.cc.

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

◆ 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 int  i)
protectedvirtual

Implements Garfield::ComponentFieldMap.

Definition at line 1177 of file ComponentCST.cc.

1177 {
1178
1179 if (element < 0 || element >= nElements) return 0.;
1180 unsigned int i,j,k;
1181 Element2Index(element, i, j, k);
1182 const double volume =
1183 fabs((m_xlines.at(i+1) - m_xlines.at(i)) *
1184 (m_xlines.at(j+1) - m_ylines.at(j)) *
1185 (m_xlines.at(k+1) - m_zlines.at(k)));
1186 return volume;
1187}
fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:616

◆ GetMedium()

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

Implements Garfield::ComponentFieldMap.

Definition at line 1034 of file ComponentCST.cc.

1035 {
1036 unsigned int i, j, k;
1037 Coordinate2Index(xin,yin,zin,i,j,k);
1038 if(debug){
1039 std::cout << m_className << "::GetMedium:" << std::endl;
1040 std::cout << " Found position (" << xin << ", " << yin << ", " << zin << "): " << std:: endl;
1041 std::cout << " Indexes are: x: " << i << "/" << m_xlines.size()
1042 << "\t y: " << j << "/" << m_ylines.size()
1043 << "\t z: " << k << "/" << m_zlines.size() << std::endl;
1044 std::cout << " Element material index: " << Index2Element(i, j, k) << std::endl;
1045 std::cout << " Element index: " << (int)m_elementMaterial.at(Index2Element(i,j,k)) << std::endl;
1046 }
1047 return materials.at(m_elementMaterial.at(Index2Element(i,j,k))).medium;
1048}
Garfield::ComponentCST::Index2Element
int Index2Element(const unsigned int i, const unsigned int j, const unsigned int k)
Definition: ComponentCST.cc:1100
Garfield::ComponentFieldMap::materials
std::vector< material > materials
Definition: ComponentFieldMap.hh:124

◆ GetNumberOfMeshLines()

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

Definition at line 1016 of file ComponentCST.cc.

1016 {
1017 n_x = m_xlines.size();
1018 n_y = m_ylines.size();
1019 n_z = m_zlines.size();
1020}

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

1100 {
1101
1102 if (i>m_nx-2 || j>m_ny-2 || k>m_nz-2) {
1103 throw "FieldMap::ElementByIndex: Error. Element indexes out of bounds.";
1104 }
1105 return i+j*(m_nx-1)+k*(m_nx-1)*(m_ny-1);
1106}

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

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

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

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

◆ IsInBoundingBox()

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

Reimplemented from Garfield::ComponentFieldMap.

Definition at line 107 of file ComponentCST.hh.

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

◆ SetRange()

void Garfield::ComponentCST::SetRange ( )
virtual

Reimplemented from Garfield::ComponentFieldMap.

Definition at line 1050 of file ComponentCST.cc.

1050 {
1051 // Establish the ranges
1052 mapxmin = *m_xlines.begin();
1053 mapxmax = *(m_xlines.end()-1);
1054 mapymin = *m_ylines.begin();
1055 mapymax = *(m_ylines.end()-1);
1056 mapzmin = *m_zlines.begin();
1057 mapzmax = *(m_zlines.end()-1);
1058 mapvmin = *std::min_element(m_potential.begin(),m_potential.end());
1059 mapvmax = *std::max_element(m_potential.begin(),m_potential.end());
1060 // Set the periodicity length (maybe not needed).
1061 mapsx = fabs(mapxmax - mapxmin);
1062 mapsy = fabs(mapymax - mapymin);
1063 mapsz = fabs(mapzmax - mapzmin);
1064 // Set provisional cell dimensions.
1065 xMinBoundingBox = mapxmin;
1066 xMaxBoundingBox = mapxmax;
1067 yMinBoundingBox = mapymin;
1068 yMaxBoundingBox = mapymax;
1069 if (is3d) {
1070 zMinBoundingBox = mapzmin;
1071 zMaxBoundingBox = mapzmax;
1072 } else {
1073 mapzmin = zMinBoundingBox;
1074 mapzmax = zMaxBoundingBox;
1075 }
1076 hasBoundingBox = true;
1077
1078}

Referenced by Initialise(), and ShiftComponent().

◆ SetRangeZ()

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

Definition at line 1080 of file ComponentCST.cc.

1080 {
1081
1082 if (fabs(zmax - zmin) <= 0.) {
1083 std::cerr << m_className << "::SetRangeZ:" << std::endl;
1084 std::cerr << " Zero range is not permitted." << std::endl;
1085 return;
1086 }
1087 zMinBoundingBox = std::min(zmin, zmax);
1088 zMaxBoundingBox = std::max(zmin, zmax);
1089}

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

689 {
690 std::vector<float> potentials(nNodes);
691 if (!ready) {
692 std::cerr << m_className << "::SetWeightingField:" << std::endl;
693 std::cerr << " No valid field map is present." << std::endl;
694 std::cerr << " Weighting field cannot be added." << std::endl;
695 return false;
696 }
697
698 // Open the voltage list
699 std::ifstream fprnsol;
700 fprnsol.open(prnsol.c_str(), std::ios::in);
701 if (fprnsol.fail()) {
702 std::cerr << m_className << "::SetWeightingField:" << std::endl;
703 std::cerr << " Could not open potential file " << prnsol
704 << " for reading." << std::endl;
705 std::cerr << " The file perhaps does not exist." << std::endl;
706 return false;
707 }
708 // Check if a weighting field with the same label already exists.
709 std::map<std::string, std::vector<float> >::iterator it = m_weightingFields.find(label);
710 if (it != m_weightingFields.end()) {
711 std::cout << m_className << "::SetWeightingField:" << std::endl;
712 std::cout << " Replacing existing weighting field " << label << "."
713 << std::endl;
714 } else {
715 wfields.push_back(label);
716 wfieldsOk.push_back(false);
717 }
718
719 if(std::distance(m_weightingFields.begin(),it) != std::distance(wfields.begin(),find(wfields.begin(),wfields.end(),label))){
720 std::cerr << m_className << "::SetWeightingField:" << std::endl;
721 std::cerr << " Indexes of the weighting fields and the weighting field counter are not equal!" << std::endl;
722 return false;
723 }
724 unsigned int iField = std::distance(m_weightingFields.begin(),it);
725 int nread = 0;
726 bool ok = true;
727
728 if(isBinary) {
729 std::cout << m_className << "::SetWeightingField:" << std::endl;
730 std::cout << " Reading weighting field from binary file:" << prnsol.c_str() << std::endl;
731 FILE* f = fopen(prnsol.c_str(), "rb");
732 if (f == nullptr) {
733 std::cerr << m_className << "::Initilise:" << std::endl;
734 std::cerr << " Could not open file:" << prnsol.c_str() << std::endl;
735 return false;
736 }
737
738 struct stat fileStatus;
739 stat(prnsol.c_str(), &fileStatus);
740 int fileSize = fileStatus.st_size;
741
742 if (fileSize < 1000) {
743 fclose(f);
744 std::cerr << m_className << "::SetWeightingField:" << std::endl;
745 std::cerr << " Error. The file is extremely short and does not seem to contain a header or data." << std::endl;
746 ok = false;
747 }
748
749 char header[headerSize];
750 size_t result;
751 result = fread(header, sizeof(char), headerSize, f);
752 if (result != headerSize) {fputs ("Reading error while reading header.",stderr); exit (3);}
753
754 int nx = 0, ny = 0, nz = 0;
755 int m_x = 0, m_y = 0, m_z = 0;
756 int n_x = 0, n_y = 0, n_z = 0;
757 int e_s = 0, e_x = 0, e_y = 0, e_z = 0;
758 int e_m = 0;
759
760 int filled = 0;
761 filled = std::sscanf(header, (std::string("mesh_nx=%d mesh_ny=%d mesh_nz=%d\n")+
762 std::string("mesh_xlines=%d mesh_ylines=%d mesh_zlines=%d\n")+
763 std::string("nodes_scalar=%d nodes_vector_x=%d nodes_vector_y=%d nodes_vector_z=%d\n")+
764 std::string("elements_scalar=%d elements_vector_x=%d elements_vector_y=%d elements_vector_z=%d\n")+
765 std::string("elements_material=%d\n")+
766 std::string("n_materials=%d\n")).c_str(),
767 &nx, &ny, &nz,
768 &m_x, &m_y, &m_z,
769 &nread, &n_x, &n_y, &n_z,
770 &e_s, &e_x, &e_y, &e_z,
771 &e_m, &nMaterials);
772 if (filled != 16){
773 fclose(f);
774 std::cerr << m_className << "::SetWeightingField:" << std::endl;
775 std::cerr << " Error. File header of " << prnsol.c_str() << " is broken." << std::endl;
776 ok = false;
777 }
778 if (fileSize < 1000+(m_x+m_y+m_z)*8+(nread+n_x+n_y+n_z+e_s+e_x+e_y+e_z)*4+e_m*1+nMaterials*20) {
779 fclose(f);
780 ok = false;
781 }
782 if(debug){
783 std::cout << m_className << "::SetWeightingField:" << std::endl;
784 std::cout << " Information about the data stored in the given binary file:" << std::endl;
785 std::cout << " Mesh (nx): " << nx << "\t Mesh (ny): " << ny << "\t Mesh (nz): " << nz << std::endl;
786 std::cout << " Mesh (x_lines): " << m_x << "\t Mesh (y_lines): " << m_y << "\t Mesh (z_lines): " << m_z << std::endl;
787 std::cout << " Nodes (scalar): " << nread << "\t Nodes (x): " << n_x << "\t Nodes (y): " << n_y << "\t Nodes (z): " << n_z << std::endl;
788 std::cout << " Field (scalar): " << e_s << "\t Field (x): " << e_x << "\t Field (y): " << e_y << "\t Field (z): " << e_z << std::endl;
789 std::cout << " Elements: " << e_m << "\t Materials: " << nMaterials << std::endl;
790 }
791 // skip everything, but the potential
792 fseek(f, m_x*sizeof(double), SEEK_CUR);
793 fseek(f, m_y*sizeof(double), SEEK_CUR);
794 fseek(f, m_z*sizeof(double), SEEK_CUR);
795 result = fread(potentials.data(), sizeof(float), potentials.size(), f);
796 if (result != potentials.size()) {fputs ("Reading error while reading nodes.",stderr); exit (3);}
797 else if (result == 0) {fputs ("No wighting potentials are stored in the data file.",stderr); exit (3);}
798 fprnsol.close();
799 } else {
800 std::cout << m_className << "::SetWeightingField:" << std::endl;
801 std::cout << " Reading weighting field from text file:" << prnsol.c_str() << std::endl;
802 // Buffer for reading
803 const int size = 100;
804 char line[size];
805
806
807 // Read the voltage list
808 int il = 0;
809
810 bool readerror = false;
811 while (fprnsol.getline(line, size, '\n')) {
812 il++;
813 // Split the line in tokens
814 char* token = NULL;
815 token = strtok(line, " ");
816 // Skip blank lines and headers
817 if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
818 int(token[0]) == 10 || int(token[0]) == 13 ||
819 strcmp(token, "PRINT") == 0 || strcmp(token, "*****") == 0 ||
820 strcmp(token, "LOAD") == 0 || strcmp(token, "TIME=") == 0 ||
821 strcmp(token, "MAXIMUM") == 0 || strcmp(token, "VALUE") == 0 ||
822 strcmp(token, "NODE") == 0)
823 continue;
824 // Read the element
825 int inode = ReadInteger(token, -1, readerror);
826 token = strtok(NULL, " ");
827 double volt = ReadDouble(token, -1, readerror);
828 try{
829 potentials.at(inode-1) = volt;
830 nread++;
831 } catch (...){
832 std::cerr << m_className << "::SetWeightingField:" << std::endl;
833 std::cerr << " Node number " << inode << " out of range." << std::endl;
834 std::cerr << " on potential file " << prnsol << " (line " << il << ")." << std::endl;
835 std::cerr << " Size of the potential vector is: " << potentials.size() << std::endl;
836 ok = false;
837 }
838 }
839 // Close the file
840 fprnsol.close();
841 }
842 // Tell how many lines read
843 std::cout << m_className << "::SetWeightingField:" << std::endl;
844 std::cout << " Read " << nread << "/" << nNodes << " (expected) potentials from file " << prnsol << "."
845 << std::endl;
846 // Check number of nodes
847 if (nread != nNodes) {
848 std::cerr << m_className << "::SetWeightingField:" << std::endl;
849 std::cerr << " Number of nodes read (" << nread << ")"
850 << " on potential file (" << prnsol << ")" << std::endl;
851 std::cerr << " does not match the node list (" << nNodes << ")."
852 << std::endl;
853 ok = false;
854 }
855 if (!ok) {
856 std::cerr << m_className << "::SetWeightingField:" << std::endl;
857 std::cerr << " Field map could not be read "
858 << "and cannot be interpolated." << std::endl;
859 return false;
860 }
861
862 m_weightingFields[label] = potentials;
863
864 // Set the ready flag.
865 wfieldsOk[iField] = ok;
866 return true;
867}
Garfield::ComponentFieldMap::wfields
std::vector< std::string > wfields
Definition: ComponentFieldMap.hh:127

◆ ShiftComponent()

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

Definition at line 869 of file ComponentCST.cc.

869 {
870 std::transform(m_xlines.begin(), m_xlines.end(), m_xlines.begin(), std::bind1st(std::plus<double>(),xShift));
871 std::transform(m_ylines.begin(), m_ylines.end(), m_ylines.begin(), std::bind1st(std::plus<double>(),yShift));
872 std::transform(m_zlines.begin(), m_zlines.end(), m_zlines.begin(), std::bind1st(std::plus<double>(),zShift));
873 SetRange();
874 UpdatePeriodicity();
875
876 std::cout << m_className << "::ShiftComponent:" << std::endl;
877 std::cout << " Shifted component in x-direction: " << xShift
878 << "\t y-direction: " << yShift
879 << "\t z-direction: " << zShift << std::endl;
880}

◆ UpdatePeriodicity()

void Garfield::ComponentCST::UpdatePeriodicity ( )
protectedvirtual

Implements Garfield::ComponentFieldMap.

Definition at line 1154 of file ComponentCST.cc.

1154 {
1155
1156 UpdatePeriodicity2d();
1157 UpdatePeriodicityCommon();
1158}

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

Implements Garfield::ComponentFieldMap.

Definition at line 895 of file ComponentCST.cc.

897 {
898
899 // Initial values
900 wx = wy = wz = 0;
901
902 // Do not proceed if not properly initialised.
903 if (!ready) return;
904
905 // Look for the label.
906 std::map<std::string, std::vector<float> >::iterator it = m_weightingFields.find(label);
907 if(it == m_weightingFields.end()){
908 // Do not proceed if the requested weighting field does not exist.
909 std::cerr << "No weighting field named " << label.c_str() << " found!" << std::endl;
910 return;
911 }
912
913 // Check if the weighting field is properly initialised.
914 if (!wfieldsOk[std::distance(m_weightingFields.begin(),it)]) return;
915
916 // Copy the coordinates
917 double x = xin, y = yin, z = zin;
918
919 // Map the coordinates onto field map coordinates and get indexes
920 bool mirrored[3];
921 double rcoordinate, rotation;
922 unsigned int i,j,k;
923 double position_mapped[3] = {0., 0., 0.};
924 if(!Coordinate2Index(x, y, z, i, j, k, position_mapped, mirrored, rcoordinate, rotation))
925 return;
926
927 double rx = (position_mapped[0] - m_xlines.at(i))/(m_xlines.at(i+1) - m_xlines.at(i));
928 double ry = (position_mapped[1] - m_ylines.at(j))/(m_ylines.at(j+1) - m_ylines.at(j));
929 double rz = (position_mapped[2] - m_zlines.at(k))/(m_zlines.at(k+1) - m_zlines.at(k));
930
931 float fwx, fwy, fwz;
932 if(!disableFieldComponent[0])
933 fwx = GetFieldComponent(i, j, k, rx, ry, rz, 'x', &((*it).second));
934 if(!disableFieldComponent[1])
935 fwy = GetFieldComponent(i, j, k, rx, ry, rz, 'y', &((*it).second));
936 if(!disableFieldComponent[2])
937 fwz = GetFieldComponent(i, j, k, rx, ry, rz, 'z', &((*it).second));
938
939 if (m_elementMaterial.size()>0 && doShaping) {
940 ShapeField(fwx, fwy, fwz, rx, ry, rz, i, j, k, &((*it).second));
941 }
942 if(mirrored[0])
943 fwx *= -1.;
944 if(mirrored[1])
945 fwy *= -1.;
946 if(mirrored[2])
947 fwz *= -1.;
948 if (warning) {
949 std::cout << m_className << "::WeightingField:" << std::endl;
950 std::cout << " Warnings have been issued for this field map."
951 << std::endl;
952 }
953 if (materials.at(m_elementMaterial.at(Index2Element(i,j,k))).driftmedium) {
954 if (!disableFieldComponent[0]) wx = fwx;
955 if (!disableFieldComponent[1]) wy = fwy;
956 if (!disableFieldComponent[2]) wz = fwz;
957 }
958}

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

962 {
963
964 // Do not proceed if not properly initialised.
965 if (!ready) return 0.;
966
967
968 // Look for the label.
969 std::map<std::string, std::vector<float> >::iterator it = m_weightingFields.find(label);
970 if(it == m_weightingFields.end()){
971 // Do not proceed if the requested weighting field does not exist.
972 std::cerr << "No weighting field named " << label.c_str() << " found!" << std::endl;
973 return 0.;
974 }
975
976 // Check if the weighting field is properly initialised.
977 if (!wfieldsOk[std::distance(m_weightingFields.begin(),it)]) return 0.;
978
979 // Copy the coordinates
980 double x = xin, y = yin, z = zin;
981
982 // Map the coordinates onto field map coordinates
983 bool mirrored[3];
984 double rcoordinate, rotation;
985 unsigned int i,j,k;
986 double position_mapped[3] = {0., 0., 0.};
987 if(!Coordinate2Index(x, y, z, i, j, k, position_mapped, mirrored, rcoordinate, rotation))
988 return 0.;
989
990 double rx = (position_mapped[0] - m_xlines.at(i))/(m_xlines.at(i+1) - m_xlines.at(i));
991 double ry = (position_mapped[1] - m_ylines.at(j))/(m_ylines.at(j+1) - m_ylines.at(j));
992 double rz = (position_mapped[2] - m_zlines.at(k))/(m_zlines.at(k+1) - m_zlines.at(k));
993
994 double potential = GetPotential(i, j, k, rx, ry, rz, &((*it).second));
995
996 if (debug) {
997 std::cout << m_className << "::WeightingPotential:" << std::endl;
998 std::cout << " Global: (" << x << "," << y << "," << z << "),"
999 << std::endl;
1000 std::cout << " Local: (" << rx << "," << ry << "," << rz
1001 << ") in element with indexes: i=" << i << ", j=" << j << ", k=" << k << std::endl;
1002 std::cout << " Node xyzV:" << std::endl;
1003 std::cout << "Node 0 position: " << Index2Node(i+1, j , k ) << "\t potential: " << ((*it).second).at(Index2Node(i+1, j , k ))
1004 << "Node 1 position: " << Index2Node(i+1, j+1, k ) << "\t potential: " << ((*it).second).at(Index2Node(i+1, j+1, k ))
1005 << "Node 2 position: " << Index2Node(i , j+1, k ) << "\t potential: " << ((*it).second).at(Index2Node(i , j+1, k ))
1006 << "Node 3 position: " << Index2Node(i , j , k ) << "\t potential: " << ((*it).second).at(Index2Node(i , j , k ))
1007 << "Node 4 position: " << Index2Node(i+1, j , k+1) << "\t potential: " << ((*it).second).at(Index2Node(i+1, j , k+1))
1008 << "Node 5 position: " << Index2Node(i+1, j+1, k+1) << "\t potential: " << ((*it).second).at(Index2Node(i+1, j+1, k+1))
1009 << "Node 6 position: " << Index2Node(i , j+1, k+1) << "\t potential: " << ((*it).second).at(Index2Node(i , j+1, k+1))
1010 << "Node 7 position: " << Index2Node(i , j , k+1) << "\t potential: " << ((*it).second).at(Index2Node(i , j , k ))
1011 << std::endl;
1012 }
1013 return potential;
1014}
The documentation for this class was generated from the following files: