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

#include <ComponentAnsys123.hh>

+ Inheritance diagram for Garfield::ComponentAnsys123:

Public Member Functions

 ComponentAnsys123 ()
 
 ~ComponentAnsys123 ()
 
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)
 
Medium * GetMedium (const double &x, const double &y, const double &z)
 
bool IsInBoundingBox (const double x, const double y, const double z)
 
bool Initialise (std::string elist="ELIST.lis", std::string nlist="NLIST.lis", std::string mplist="MPLIST.lis", std::string prnsol="PRNSOL.lis", std::string unit="cm")
 
bool SetWeightingField (std::string prnsol, std::string label)
 
- 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)
 
- 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 8 of file ComponentAnsys123.hh.

Constructor & Destructor Documentation

◆ ComponentAnsys123()

Garfield::ComponentAnsys123::ComponentAnsys123 ( )

Definition at line 12 of file ComponentAnsys123.cc.

12 : ComponentFieldMap() {
13
14 m_className = "ComponentAnsys123";
15 ready = false;
16}

◆ ~ComponentAnsys123()

Garfield::ComponentAnsys123::~ComponentAnsys123 ( )
inline

Definition at line 14 of file ComponentAnsys123.hh.

14{}

Member Function Documentation

◆ ElectricField() [1/2]

void Garfield::ComponentAnsys123::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 795 of file ComponentAnsys123.cc.

798 {
799
800 // Copy the coordinates
801 double x = xin, y = yin, z = zin;
802
803 // Map the coordinates onto field map coordinates
804 bool xmirrored, ymirrored, zmirrored;
805 double rcoordinate, rotation;
806 MapCoordinates(x, y, z, xmirrored, ymirrored, zmirrored, rcoordinate,
807 rotation);
808
809 // Initial values
810 ex = ey = ez = volt = 0.;
811 status = 0;
812 m = 0;
813
814 // Do not proceed if not properly initialised.
815 if (!ready) {
816 status = -10;
817 std::cerr << m_className << "::ElectricField:\n";
818 std::cerr << " Field map not available for interpolation.\n";
819 return;
820 }
821
822 if (warning) {
823 std::cerr << m_className << "::ElectricField:\n";
824 std::cerr << " Warnings have been issued for this field map.\n";
825 }
826
827 // Find the element that contains this point
828 double t1, t2, t3, t4, jac[4][4], det;
829 int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
830 if (imap < 0) {
831 if (debug) {
832 std::cerr << m_className << "::ElectricField:\n";
833 std::cerr << " Point (" << x << ", " << y << ", " << z
834 << ") not in the mesh.\n";
835 }
836 status = -6;
837 return;
838 }
839
840 if (debug) {
841 std::cout << m_className << "::ElectricField:\n";
842 std::cout << " Global: (" << x << ", " << y << ", " << z << "),\n";
843 std::cout << " Local: (" << t1 << ", " << t2 << ", " << t3 << ", " << t4
844 << ") in element " << imap << ".\n";
845 std::cout
846 << " Node x y z V\n";
847 for (int i = 0; i < 10; i++) {
848 printf(" %-5d %12g %12g %12g %12g\n", elements[imap].emap[i],
849 nodes[elements[imap].emap[i]].x, nodes[elements[imap].emap[i]].y,
850 nodes[elements[imap].emap[i]].z, nodes[elements[imap].emap[i]].v);
851 }
852 }
853
854 // Tetrahedral field
855 volt = nodes[elements[imap].emap[0]].v * t1 * (2 * t1 - 1) +
856 nodes[elements[imap].emap[1]].v * t2 * (2 * t2 - 1) +
857 nodes[elements[imap].emap[2]].v * t3 * (2 * t3 - 1) +
858 nodes[elements[imap].emap[3]].v * t4 * (2 * t4 - 1) +
859 4 * nodes[elements[imap].emap[4]].v * t1 * t2 +
860 4 * nodes[elements[imap].emap[5]].v * t1 * t3 +
861 4 * nodes[elements[imap].emap[6]].v * t1 * t4 +
862 4 * nodes[elements[imap].emap[7]].v * t2 * t3 +
863 4 * nodes[elements[imap].emap[8]].v * t2 * t4 +
864 4 * nodes[elements[imap].emap[9]].v * t3 * t4;
865
866 ex = -(nodes[elements[imap].emap[0]].v * (4 * t1 - 1) * jac[0][1] +
867 nodes[elements[imap].emap[1]].v * (4 * t2 - 1) * jac[1][1] +
868 nodes[elements[imap].emap[2]].v * (4 * t3 - 1) * jac[2][1] +
869 nodes[elements[imap].emap[3]].v * (4 * t4 - 1) * jac[3][1] +
870 nodes[elements[imap].emap[4]].v *
871 (4 * t2 * jac[0][1] + 4 * t1 * jac[1][1]) +
872 nodes[elements[imap].emap[5]].v *
873 (4 * t3 * jac[0][1] + 4 * t1 * jac[2][1]) +
874 nodes[elements[imap].emap[6]].v *
875 (4 * t4 * jac[0][1] + 4 * t1 * jac[3][1]) +
876 nodes[elements[imap].emap[7]].v *
877 (4 * t3 * jac[1][1] + 4 * t2 * jac[2][1]) +
878 nodes[elements[imap].emap[8]].v *
879 (4 * t4 * jac[1][1] + 4 * t2 * jac[3][1]) +
880 nodes[elements[imap].emap[9]].v *
881 (4 * t4 * jac[2][1] + 4 * t3 * jac[3][1])) /
882 det;
883
884 ey = -(nodes[elements[imap].emap[0]].v * (4 * t1 - 1) * jac[0][2] +
885 nodes[elements[imap].emap[1]].v * (4 * t2 - 1) * jac[1][2] +
886 nodes[elements[imap].emap[2]].v * (4 * t3 - 1) * jac[2][2] +
887 nodes[elements[imap].emap[3]].v * (4 * t4 - 1) * jac[3][2] +
888 nodes[elements[imap].emap[4]].v *
889 (4 * t2 * jac[0][2] + 4 * t1 * jac[1][2]) +
890 nodes[elements[imap].emap[5]].v *
891 (4 * t3 * jac[0][2] + 4 * t1 * jac[2][2]) +
892 nodes[elements[imap].emap[6]].v *
893 (4 * t4 * jac[0][2] + 4 * t1 * jac[3][2]) +
894 nodes[elements[imap].emap[7]].v *
895 (4 * t3 * jac[1][2] + 4 * t2 * jac[2][2]) +
896 nodes[elements[imap].emap[8]].v *
897 (4 * t4 * jac[1][2] + 4 * t2 * jac[3][2]) +
898 nodes[elements[imap].emap[9]].v *
899 (4 * t4 * jac[2][2] + 4 * t3 * jac[3][2])) /
900 det;
901
902 ez = -(nodes[elements[imap].emap[0]].v * (4 * t1 - 1) * jac[0][3] +
903 nodes[elements[imap].emap[1]].v * (4 * t2 - 1) * jac[1][3] +
904 nodes[elements[imap].emap[2]].v * (4 * t3 - 1) * jac[2][3] +
905 nodes[elements[imap].emap[3]].v * (4 * t4 - 1) * jac[3][3] +
906 nodes[elements[imap].emap[4]].v *
907 (4 * t2 * jac[0][3] + 4 * t1 * jac[1][3]) +
908 nodes[elements[imap].emap[5]].v *
909 (4 * t3 * jac[0][3] + 4 * t1 * jac[2][3]) +
910 nodes[elements[imap].emap[6]].v *
911 (4 * t4 * jac[0][3] + 4 * t1 * jac[3][3]) +
912 nodes[elements[imap].emap[7]].v *
913 (4 * t3 * jac[1][3] + 4 * t2 * jac[2][3]) +
914 nodes[elements[imap].emap[8]].v *
915 (4 * t4 * jac[1][3] + 4 * t2 * jac[3][3]) +
916 nodes[elements[imap].emap[9]].v *
917 (4 * t4 * jac[2][3] + 4 * t3 * jac[3][3])) /
918 det;
919
920 // Transform field to global coordinates
921 UnmapFields(ex, ey, ez, x, y, z, xmirrored, ymirrored, zmirrored, rcoordinate,
922 rotation);
923
924 // Drift medium?
925 if (debug) {
926 std::cout << m_className << "::ElectricField:\n";
927 std::cout << " Material " << elements[imap].matmap << ", drift flag "
928 << materials[elements[imap].matmap].driftmedium << ".\n";
929 }
930 m = materials[elements[imap].matmap].medium;
931 status = -5;
932 if (materials[elements[imap].matmap].driftmedium) {
933 if (m != 0) {
934 if (m->IsDriftable()) status = 0;
935 }
936 }
937}
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
Garfield::ComponentFieldMap::materials
std::vector< material > materials
Definition: ComponentFieldMap.hh:124
Garfield::ComponentFieldMap::UnmapFields
void UnmapFields(double &ex, double &ey, double &ez, double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation)
Definition: ComponentFieldMap.cc:4154
Garfield::ComponentFieldMap::elements
std::vector< element > elements
Definition: ComponentFieldMap.hh:96
Garfield::ComponentFieldMap::FindElement13
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)
Definition: ComponentFieldMap.cc:356

◆ ElectricField() [2/2]

void Garfield::ComponentAnsys123::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 787 of file ComponentAnsys123.cc.

789 {
790
791 double v = 0.;
792 ElectricField(x, y, z, ex, ey, ez, v, m, status);
793}
Garfield::ComponentAnsys123::ElectricField
void ElectricField(const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status)
Definition: ComponentAnsys123.cc:787

Referenced by ElectricField().

◆ GetAspectRatio()

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

Implements Garfield::ComponentFieldMap.

Definition at line 1222 of file ComponentAnsys123.cc.

1223 {
1224
1225 if (i < 0 || i >= nElements) {
1226 dmin = dmax = 0.;
1227 return;
1228 }
1229
1230 const int np = 4;
1231 // Loop over all pairs of vertices.
1232 for (int j = 0; j < np - 1; ++j) {
1233 for (int k = j + 1; k < np; ++k) {
1234 // Compute distance.
1235 const double dist = sqrt(
1236 pow(nodes[elements[i].emap[j]].x - nodes[elements[i].emap[k]].x, 2) +
1237 pow(nodes[elements[i].emap[j]].y - nodes[elements[i].emap[k]].y, 2) +
1238 pow(nodes[elements[i].emap[j]].z - nodes[elements[i].emap[k]].z, 2));
1239 if (k == 1) {
1240 dmin = dist;
1241 dmax = dist;
1242 } else {
1243 if (dist < dmin) dmin = dist;
1244 if (dist > dmax) dmax = dist;
1245 }
1246 }
1247 }
1248}
pow
DoubleAc pow(const DoubleAc &f, double p)
Definition: DoubleAc.cpp:336
sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:313

◆ GetElementVolume()

double Garfield::ComponentAnsys123::GetElementVolume ( const int  i)
protectedvirtual

Implements Garfield::ComponentFieldMap.

Definition at line 1192 of file ComponentAnsys123.cc.

1192 {
1193
1194 if (i < 0 || i >= nElements) return 0.;
1195
1196 const double vol =
1197 fabs((nodes[elements[i].emap[3]].x - nodes[elements[i].emap[0]].x) *
1198 ((nodes[elements[i].emap[1]].y - nodes[elements[i].emap[0]].y) *
1199 (nodes[elements[i].emap[2]].z -
1200 nodes[elements[i].emap[0]].z) -
1201 (nodes[elements[i].emap[2]].y - nodes[elements[i].emap[0]].y) *
1202 (nodes[elements[i].emap[1]].z -
1203 nodes[elements[i].emap[0]].z)) +
1204 (nodes[elements[i].emap[3]].y - nodes[elements[i].emap[0]].y) *
1205 ((nodes[elements[i].emap[1]].z - nodes[elements[i].emap[0]].z) *
1206 (nodes[elements[i].emap[2]].x -
1207 nodes[elements[i].emap[0]].x) -
1208 (nodes[elements[i].emap[2]].z - nodes[elements[i].emap[0]].z) *
1209 (nodes[elements[i].emap[1]].x -
1210 nodes[elements[i].emap[0]].x)) +
1211 (nodes[elements[i].emap[3]].z - nodes[elements[i].emap[0]].z) *
1212 ((nodes[elements[i].emap[1]].x - nodes[elements[i].emap[0]].x) *
1213 (nodes[elements[i].emap[2]].y -
1214 nodes[elements[i].emap[0]].y) -
1215 (nodes[elements[i].emap[3]].x - nodes[elements[i].emap[0]].x) *
1216 (nodes[elements[i].emap[1]].y -
1217 nodes[elements[i].emap[0]].y))) /
1218 6.;
1219 return vol;
1220}
fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:616

◆ GetMedium()

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

Implements Garfield::ComponentFieldMap.

Definition at line 1130 of file ComponentAnsys123.cc.

1131 {
1132
1133 // Copy the coordinates.
1134 double x = xin, y = yin, z = zin;
1135
1136 // Map the coordinates onto field map coordinates.
1137 bool xmirrored, ymirrored, zmirrored;
1138 double rcoordinate, rotation;
1139 MapCoordinates(x, y, z, xmirrored, ymirrored, zmirrored, rcoordinate,
1140 rotation);
1141
1142 // Do not proceed if not properly initialised.
1143 if (!ready) {
1144 std::cerr << m_className << "::GetMedium:\n";
1145 std::cerr << " Field map not available for interpolation.\n";
1146 return NULL;
1147 }
1148 if (warning) {
1149 std::cerr << m_className << "::GetMedium:\n";
1150 std::cerr << " Warnings have been issued for this field map.\n";
1151 }
1152
1153 // Find the element that contains this point.
1154 double t1, t2, t3, t4, jac[4][4], det;
1155 int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
1156 if (imap < 0) {
1157 if (debug) {
1158 std::cerr << m_className << "::GetMedium:\n";
1159 std::cerr << " Point (" << x << ", " << y << ", " << z
1160 << ") not in the mesh.\n";
1161 }
1162 return NULL;
1163 }
1164 if (elements[imap].matmap < 0 || elements[imap].matmap >= nMaterials) {
1165 if (debug) {
1166 std::cerr << m_className << "::GetMedium:\n";
1167 std::cerr << " Point (" << x << ", " << y << ", " << z << ")"
1168 << " has out of range material number " << imap << ".\n";
1169 }
1170 return NULL;
1171 }
1172
1173 if (debug) {
1174 std::cout << m_className << "::GetMedium:\n";
1175 std::cout << " Global: (" << x << ", " << y << ", " << z << ")\n";
1176 std::cout << " Local: (" << t1 << ", " << t2 << ", " << t3 << ", " << t4
1177 << ") in element " << imap << "\n";
1178 std::cout
1179 << " Node x y z V\n";
1180 for (int ii = 0; ii < 10; ++ii) {
1181 printf(" %-5d %12g %12g %12g %12g\n", elements[imap].emap[ii],
1182 nodes[elements[imap].emap[ii]].x, nodes[elements[imap].emap[ii]].y,
1183 nodes[elements[imap].emap[ii]].z,
1184 nodes[elements[imap].emap[ii]].v);
1185 }
1186 }
1187
1188 // Assign a medium.
1189 return materials[elements[imap].matmap].medium;
1190}

◆ Initialise()

bool Garfield::ComponentAnsys123::Initialise ( std::string  elist = "ELIST.lis",
std::string  nlist = "NLIST.lis",
std::string  mplist = "MPLIST.lis",
std::string  prnsol = "PRNSOL.lis",
std::string  unit = "cm" 
)

Definition at line 18 of file ComponentAnsys123.cc.

20 {
21
22 ready = false;
23 // Keep track of the success.
24 bool ok = true;
25
26 // Buffer for reading
27 const int size = 100;
28 char line[size];
29
30 // Open the material list.
31 std::ifstream fmplist;
32 fmplist.open(mplist.c_str(), std::ios::in);
33 if (fmplist.fail()) {
34 std::cerr << m_className << "::Initialise:\n";
35 std::cerr << " Could not open material file " << mplist
36 << " for reading.\n",
37 std::cerr << " The file perhaps does not exist.\n";
38 return false;
39 }
40
41 // Read the material list.
42 nMaterials = 0;
43 int il = 0, icurrmat = -1;
44 bool readerror = false;
45 while (fmplist.getline(line, size, '\n')) {
46 il++;
47 // Skip page feed
48 if (strcmp(line, "1") == 0) {
49 fmplist.getline(line, size, '\n');
50 il++;
51 fmplist.getline(line, size, '\n');
52 il++;
53 fmplist.getline(line, size, '\n');
54 il++;
55 fmplist.getline(line, size, '\n');
56 il++;
57 fmplist.getline(line, size, '\n');
58 il++;
59 continue;
60 }
61 // Split the line in tokens
62 char* token = NULL;
63 token = strtok(line, " ");
64 // Skip blank lines and headers
65 if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
66 strcmp(token, "TEMPERATURE") == 0 || strcmp(token, "PROPERTY=") == 0 ||
67 int(token[0]) == 10 || int(token[0]) == 13)
68 continue;
69 // Read number of materials,
70 // ensure it does not exceed the maximum and initialise the list
71 if (strcmp(token, "LIST") == 0) {
72 token = strtok(NULL, " ");
73 token = strtok(NULL, " ");
74 token = strtok(NULL, " ");
75 token = strtok(NULL, " ");
76 nMaterials = ReadInteger(token, -1, readerror);
77 if (readerror) {
78 std::cerr << m_className << "::Initialise:\n";
79 std::cerr << " Error reading file " << mplist << " (line " << il
80 << ").\n";
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:\n";
93 std::cout << " Number of materials: " << nMaterials << "\n";
94 }
95 } else if (strcmp(token, "MATERIAL") == 0) {
96 // Version 12 format: read material number
97 token = strtok(NULL, " ");
98 token = strtok(NULL, " ");
99 icurrmat = ReadInteger(token, -1, readerror);
100 if (readerror) {
101 std::cerr << m_className << "::Initialise:\n";
102 std::cerr << " Error reading file " << mplist << " (line " << il
103 << ").\n";
104 fmplist.close();
105 ok = false;
106 return false;
107 }
108 } else if (strcmp(token, "TEMP") == 0) {
109 // Version 12 format: read property tag and value
110 token = strtok(NULL, " ");
111 int itype = 0;
112 if (strncmp(token, "PERX", 4) == 0) {
113 itype = 1;
114 } else if (strncmp(token, "RSVX", 4) == 0) {
115 itype = 2;
116 } else {
117 std::cerr << m_className << "::Initialise:\n";
118 std::cerr << " Found unknown material property flag " << token
119 << "\n";
120 std::cerr << " on material properties file " << mplist << " (line "
121 << il << ").\n";
122 ok = false;
123 }
124 fmplist.getline(line, size, '\n');
125 il++;
126 token = NULL;
127 token = strtok(line, " ");
128 if (icurrmat < 1 || icurrmat > nMaterials) {
129 std::cerr << m_className << "::Initialise:\n";
130 std::cerr << " Found out-of-range current material index "
131 << icurrmat << "\n";
132 std::cerr << " in material properties file " << mplist << ".\n";
133 ok = false;
134 readerror = false;
135 } else if (itype == 1) {
136 materials[icurrmat - 1].eps = ReadDouble(token, -1, readerror);
137 } else if (itype == 2) {
138 materials[icurrmat - 1].ohm = ReadDouble(token, -1, readerror);
139 }
140 if (readerror) {
141 std::cerr << m_className << "::Initialise:\n";
142 std::cerr << " Error reading file " << mplist << " line " << il
143 << ").\n";
144 fmplist.close();
145 ok = false;
146 return false;
147 }
148 } else if (strcmp(token, "PROPERTY") == 0) {
149 // Version 11 format
150 token = strtok(NULL, " ");
151 token = strtok(NULL, " ");
152 int itype = 0;
153 if (strcmp(token, "PERX") == 0) {
154 itype = 1;
155 } else if (strcmp(token, "RSVX") == 0) {
156 itype = 2;
157 } else {
158 std::cerr << m_className << "::Initialise:\n";
159 std::cerr << " Found unknown material property flag " << token
160 << "\n";
161 std::cerr << " on material properties file " << mplist << " (line "
162 << il << ").\n";
163 ok = false;
164 }
165 token = strtok(NULL, " ");
166 token = strtok(NULL, " ");
167 int imat = ReadInteger(token, -1, readerror);
168 if (readerror) {
169 std::cerr << m_className << "::Initialise:\n";
170 std::cerr << " Error reading file " << mplist << " (line " << il
171 << ").\n";
172 fmplist.close();
173 ok = false;
174 return false;
175 } else if (imat < 1 || imat > nMaterials) {
176 std::cerr << m_className << "::Initialise:\n";
177 std::cerr << " Found out-of-range material index " << imat << "\n";
178 std::cerr << " in material properties file " << mplist << ".\n";
179 ok = false;
180 } else {
181 fmplist.getline(line, size, '\n');
182 il++;
183 fmplist.getline(line, size, '\n');
184 il++;
185 token = NULL;
186 token = strtok(line, " ");
187 token = strtok(NULL, " ");
188 if (itype == 1) {
189 materials[imat - 1].eps = ReadDouble(token, -1, readerror);
190 } else if (itype == 2) {
191 materials[imat - 1].ohm = ReadDouble(token, -1, readerror);
192 }
193 if (readerror) {
194 std::cerr << m_className << "::Initialise:\n";
195 std::cerr << " Error reading file " << mplist << " (line " << il
196 << ").\n";
197 fmplist.close();
198 ok = false;
199 return false;
200 }
201 }
202 }
203 }
204
205 // Close the file
206 fmplist.close();
207
208 // Find the lowest epsilon, check for eps = 0, set default drift media
209 double epsmin = -1;
210 int iepsmin = -1;
211 for (int imat = 0; imat < nMaterials; ++imat) {
212 if (materials[imat].eps < 0) continue;
213 if (materials[imat].eps == 0) {
214 std::cerr << m_className << "::Initialise:\n";
215 std::cerr << " Material " << imat
216 << " has been assigned a permittivity\n";
217 std::cerr << " equal to zero in " << mplist << ".\n";
218 ok = false;
219 } else if (iepsmin < 0 || epsmin > materials[imat].eps) {
220 epsmin = materials[imat].eps;
221 iepsmin = imat;
222 }
223 }
224
225 if (iepsmin < 0) {
226 std::cerr << m_className << "::Initialise:\n";
227 std::cerr << " No material with positive permittivity found \n";
228 std::cerr << " in material list " << mplist << ".\n";
229 ok = false;
230 } else {
231 for (int imat = 0; imat < nMaterials; ++imat) {
232 if (imat == iepsmin) {
233 materials[imat].driftmedium = true;
234 } else {
235 materials[imat].driftmedium = false;
236 }
237 }
238 }
239
240 // Tell how many lines read
241 std::cout << m_className << "::Initialise:\n";
242 std::cout << " Read properties of " << nMaterials
243 << " materials from file " << mplist << ".\n";
244 if (debug) PrintMaterials();
245
246 // Open the element list
247 std::ifstream felist;
248 felist.open(elist.c_str(), std::ios::in);
249 if (felist.fail()) {
250 std::cerr << m_className << "::Initialise:\n";
251 std::cerr << " Could not open element file " << elist
252 << " for reading.\n";
253 std::cerr << " The file perhaps does not exist.\n";
254 return false;
255 }
256
257 // Read the element list
258 elements.clear();
259 nElements = 0;
260 element newElement;
261 int nbackground = 0;
262 il = 0;
263 int highestnode = 0;
264 while (felist.getline(line, size, '\n')) {
265 il++;
266 // Skip page feed
267 if (strcmp(line, "1") == 0) {
268 felist.getline(line, size, '\n');
269 il++;
270 felist.getline(line, size, '\n');
271 il++;
272 felist.getline(line, size, '\n');
273 il++;
274 felist.getline(line, size, '\n');
275 il++;
276 felist.getline(line, size, '\n');
277 il++;
278 continue;
279 }
280 // Split the line in tokens
281 char* token = NULL;
282 // Split into tokens
283 token = strtok(line, " ");
284 // Skip blank lines and headers
285 if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
286 int(token[0]) == 10 || int(token[0]) == 13 ||
287 strcmp(token, "LIST") == 0 || strcmp(token, "ELEM") == 0) {
288 continue;
289 }
290 // Read the element
291 int ielem = ReadInteger(token, -1, readerror);
292 token = strtok(NULL, " ");
293 int imat = ReadInteger(token, -1, readerror);
294 token = strtok(NULL, " ");
295 token = strtok(NULL, " ");
296 token = strtok(NULL, " ");
297 token = strtok(NULL, " ");
298 token = strtok(NULL, " ");
299 int in0 = ReadInteger(token, -1, readerror);
300 token = strtok(NULL, " ");
301 int in1 = ReadInteger(token, -1, readerror);
302 token = strtok(NULL, " ");
303 int in2 = ReadInteger(token, -1, readerror);
304 token = strtok(NULL, " ");
305 int in3 = ReadInteger(token, -1, readerror);
306 token = strtok(NULL, " ");
307 int in4 = ReadInteger(token, -1, readerror);
308 token = strtok(NULL, " ");
309 int in5 = ReadInteger(token, -1, readerror);
310 token = strtok(NULL, " ");
311 int in6 = ReadInteger(token, -1, readerror);
312 token = strtok(NULL, " ");
313 int in7 = ReadInteger(token, -1, readerror);
314 if (!felist.getline(line, size, '\n')) {
315 std::cerr << m_className << "::Initialise:\n";
316 std::cerr << " Error reading element " << ielem << ".\n";
317 ok = false;
318 break;
319 }
320 token = NULL;
321 token = strtok(line, " ");
322 int in8 = ReadInteger(token, -1, readerror);
323 token = strtok(NULL, " ");
324 int in9 = ReadInteger(token, -1, readerror);
325
326 // Check synchronisation
327 if (readerror) {
328 std::cerr << m_className << "::Initialise:\n";
329 std::cerr << " Error reading file " << elist << " (line " << il
330 << ").\n";
331 felist.close();
332 ok = false;
333 return false;
334 } else if (ielem - 1 != nElements + nbackground) {
335 std::cerr << m_className << "::Initialise:\n";
336 std::cerr << " Synchronisation lost on file " << elist << " (line "
337 << il << ").\n";
338 std::cerr << " Element: " << ielem << " (expected " << nElements
339 << "), material: " << imat << ",\n";
340 std::cerr << " nodes: (" << in0 << ", " << in1 << ", " << in2 << ", "
341 << in3 << ", " << in4 << ", " << in5 << ", " << in6 << ", "
342 << in7 << ", " << in8 << ", " << in9 << ")\n";
343 ok = false;
344 }
345
346 // Check the material number and ensure that epsilon is non-negative
347 if (imat < 1 || imat > nMaterials) {
348 std::cerr << m_className << "::Initialise:\n";
349 std::cerr << " Out-of-range material number on file " << elist
350 << " (line " << il << ").\n";
351 std::cerr << " Element: " << ielem << ", material: " << imat << ",\n";
352 std::cerr << " nodes: (" << in0 << ", " << in1 << ", " << in2 << ", "
353 << in3 << ", " << in4 << ", " << in5 << ", " << in6 << ", "
354 << in7 << ", " << in8 << ", " << in9 << ")\n";
355 ok = false;
356 }
357 if (materials[imat - 1].eps < 0) {
358 std::cerr << m_className << "::Initialise:\n";
359 std::cerr << " Element " << ielem << " in element list " << elist
360 << "\n";
361 std::cerr << " uses material " << imat
362 << " which has not been assigned\n";
363 std::cerr << " a positive permittivity in material list " << mplist
364 << ".\n";
365 ok = false;
366 }
367
368 // Check the node numbers
369 if (in0 < 1 || in1 < 1 || in2 < 1 || in3 < 1 || in4 < 1 || in5 < 1 ||
370 in6 < 1 || in7 < 1 || in8 < 1 || in9 < 1) {
371 std::cerr << m_className << "::Initialise:\n";
372 std::cerr << " Found a node number < 1 on file " << elist << " (line "
373 << il << ").\n";
374 std::cerr << " Element: " << ielem << ", material: " << imat << ",\n";
375 std::cerr << " nodes: (" << in0 << ", " << in1 << ", " << in2 << ", "
376 << in3 << ", " << in4 << ", " << in5 << ", " << in6 << ", "
377 << in7 << ", " << in8 << ", " << in9 << ")\n";
378 ok = false;
379 }
380 if (in0 > highestnode) highestnode = in0;
381 if (in1 > highestnode) highestnode = in1;
382 if (in2 > highestnode) highestnode = in2;
383 if (in3 > highestnode) highestnode = in3;
384 if (in4 > highestnode) highestnode = in4;
385 if (in5 > highestnode) highestnode = in5;
386 if (in6 > highestnode) highestnode = in6;
387 if (in7 > highestnode) highestnode = in7;
388 if (in8 > highestnode) highestnode = in8;
389 if (in9 > highestnode) highestnode = in9;
390
391 // Skip quadrilaterals which are background.
392 if (deleteBackground && materials[imat - 1].ohm == 0) {
393 nbackground++;
394 continue;
395 }
396
397 // These elements must not be degenerate.
398 if (in0 == in1 || in0 == in2 || in0 == in3 || in0 == in4 || in0 == in5 ||
399 in0 == in6 || in0 == in7 || in0 == in8 || in0 == in9 || in1 == in2 ||
400 in1 == in3 || in1 == in4 || in1 == in5 || in1 == in6 || in1 == in7 ||
401 in1 == in8 || in1 == in9 || in2 == in3 || in2 == in4 || in2 == in5 ||
402 in2 == in6 || in2 == in7 || in2 == in8 || in2 == in9 || in3 == in4 ||
403 in3 == in5 || in3 == in6 || in3 == in7 || in3 == in8 || in3 == in9 ||
404 in4 == in5 || in4 == in6 || in4 == in7 || in4 == in8 || in4 == in9 ||
405 in5 == in6 || in5 == in7 || in5 == in8 || in5 == in9 || in6 == in7 ||
406 in6 == in8 || in6 == in9 || in7 == in8 || in7 == in9 || in8 == in9) {
407 std::cerr << m_className << "::Initialise:\n";
408 std::cerr << " Element " << ielem << " of file " << elist
409 << " is degenerate,\n";
410 std::cerr << " no such elements allowed in this type of map.\n";
411 ok = false;
412 }
413
414 newElement.degenerate = false;
415
416 // Store the material reference
417 newElement.matmap = imat - 1;
418
419 // Node references
420 newElement.emap[0] = in0 - 1;
421 newElement.emap[1] = in1 - 1;
422 newElement.emap[2] = in2 - 1;
423 newElement.emap[3] = in3 - 1;
424 newElement.emap[4] = in4 - 1;
425 newElement.emap[7] = in5 - 1;
426 newElement.emap[5] = in6 - 1;
427 newElement.emap[6] = in7 - 1;
428 newElement.emap[8] = in8 - 1;
429 newElement.emap[9] = in9 - 1;
430 elements.push_back(newElement);
431 nElements++;
432 }
433 // Close the file
434 felist.close();
435
436 // Tell how many lines read.
437 std::cout << m_className << "::Initialise:\n";
438 std::cout << " Read " << nElements << " elements from file " << elist
439 << ",\n";
440 std::cout << " highest node number: " << highestnode << ",\n";
441 std::cout << " background elements skipped: " << nbackground << "\n";
442 // Check the value of the unit
443 double funit;
444 if (strcmp(unit.c_str(), "mum") == 0 || strcmp(unit.c_str(), "micron") == 0 ||
445 strcmp(unit.c_str(), "micrometer") == 0) {
446 funit = 0.0001;
447 } else if (strcmp(unit.c_str(), "mm") == 0 ||
448 strcmp(unit.c_str(), "millimeter") == 0) {
449 funit = 0.1;
450 } else if (strcmp(unit.c_str(), "cm") == 0 ||
451 strcmp(unit.c_str(), "centimeter") == 0) {
452 funit = 1.0;
453 } else if (strcmp(unit.c_str(), "m") == 0 ||
454 strcmp(unit.c_str(), "meter") == 0) {
455 funit = 100.0;
456 } else {
457 std::cerr << m_className << "::Initialise:\n";
458 std::cerr << " Unknown length unit " << unit << ".\n";
459 ok = false;
460 funit = 1.0;
461 }
462 if (debug) {
463 std::cout << m_className << ":Initialise:\n";
464 std::cout << " Unit scaling factor = " << funit << ".\n";
465 }
466
467 // Open the node list
468 std::ifstream fnlist;
469 fnlist.open(nlist.c_str(), std::ios::in);
470 if (fnlist.fail()) {
471 std::cerr << m_className << "::Initialise:\n";
472 std::cerr << " Could not open nodes file " << nlist << " for reading.\n";
473 std::cerr << " The file perhaps does not exist.\n";
474 return false;
475 }
476
477 // Read the node list
478 nodes.clear();
479 nNodes = 0;
480 node newNode;
481 newNode.w.clear();
482 il = 0;
483 while (fnlist.getline(line, size, '\n')) {
484 il++;
485 // Skip page feed
486 if (strcmp(line, "1") == 0) {
487 fnlist.getline(line, size, '\n');
488 il++;
489 fnlist.getline(line, size, '\n');
490 il++;
491 fnlist.getline(line, size, '\n');
492 il++;
493 fnlist.getline(line, size, '\n');
494 il++;
495 fnlist.getline(line, size, '\n');
496 il++;
497 continue;
498 }
499 // Split the line in tokens
500 char* token = NULL;
501 token = strtok(line, " ");
502 // Skip blank lines and headers
503 if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
504 int(token[0]) == 10 || int(token[0]) == 13 ||
505 strcmp(token, "LIST") == 0 || strcmp(token, "NODE") == 0)
506 continue;
507 // Read the element
508 int inode = ReadInteger(token, -1, readerror);
509 token = strtok(NULL, " ");
510 double xnode = ReadDouble(token, -1, readerror);
511 token = strtok(NULL, " ");
512 double ynode = ReadDouble(token, -1, readerror);
513 token = strtok(NULL, " ");
514 double znode = ReadDouble(token, -1, readerror);
515 // Check syntax
516 if (readerror) {
517 std::cerr << m_className << "::Initialise:\n";
518 std::cerr << " Error reading file " << nlist << " (line " << il
519 << ").\n";
520 fnlist.close();
521 ok = false;
522 return false;
523 }
524 // Check synchronisation
525 if (inode - 1 != nNodes) {
526 std::cerr << m_className << "::Initialise:\n";
527 std::cerr << " Synchronisation lost on file " << nlist << " (line "
528 << il << ").\n";
529 std::cerr << " Node: " << inode << " (expected " << nNodes
530 << "), (x,y,z) = (" << xnode << ", " << ynode << ", " << znode
531 << ")\n";
532 ok = false;
533 }
534 // Store the point coordinates
535 newNode.x = xnode * funit;
536 newNode.y = ynode * funit;
537 newNode.z = znode * funit;
538 nodes.push_back(newNode);
539 nNodes++;
540 }
541 // Close the file
542 fnlist.close();
543 // Tell how many lines read
544 std::cout << m_className << "::Initialise:\n";
545 std::cout << " Read " << nNodes << " nodes from file " << nlist << ".\n";
546 // Check number of nodes
547 if (nNodes != highestnode) {
548 std::cerr << m_className << "::Initialise:\n";
549 std::cerr << " Number of nodes read (" << nNodes << ") on " << nlist
550 << "\n";
551 std::cerr << " does not match element list (" << highestnode << ").\n";
552 ok = false;
553 }
554
555 // Open the voltage list
556 std::ifstream fprnsol;
557 fprnsol.open(prnsol.c_str(), std::ios::in);
558 if (fprnsol.fail()) {
559 std::cerr << m_className << "::Initialise:\n";
560 std::cerr << " Could not open potential file " << prnsol
561 << " for reading.\n";
562 std::cerr << " The file perhaps does not exist.\n";
563 return false;
564 }
565
566 // Read the voltage list
567 il = 0;
568 int nread = 0;
569 while (fprnsol.getline(line, size, '\n')) {
570 il++;
571 // Skip page feed
572 if (strcmp(line, "1") == 0) {
573 fprnsol.getline(line, size, '\n');
574 il++;
575 fprnsol.getline(line, size, '\n');
576 il++;
577 fprnsol.getline(line, size, '\n');
578 il++;
579 fprnsol.getline(line, size, '\n');
580 il++;
581 fprnsol.getline(line, size, '\n');
582 il++;
583 continue;
584 }
585 // Split the line in tokens
586 char* token = NULL;
587 token = strtok(line, " ");
588 // Skip blank lines and headers
589 if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
590 int(token[0]) == 10 || int(token[0]) == 13 ||
591 strcmp(token, "PRINT") == 0 || strcmp(token, "*****") == 0 ||
592 strcmp(token, "LOAD") == 0 || strcmp(token, "TIME=") == 0 ||
593 strcmp(token, "MAXIMUM") == 0 || strcmp(token, "VALUE") == 0 ||
594 strcmp(token, "NODE") == 0)
595 continue;
596 // Read the element
597 int inode = ReadInteger(token, -1, readerror);
598 token = strtok(NULL, " ");
599 double volt = ReadDouble(token, -1, readerror);
600 // Check syntax
601 if (readerror) {
602 std::cerr << m_className << "::Initialise:\n";
603 std::cerr << " Error reading file " << prnsol << " (line << " << il
604 << ").\n";
605 fprnsol.close();
606 ok = false;
607 return false;
608 }
609 // Check node number and store if OK
610 if (inode < 1 || inode > highestnode) {
611 std::cerr << m_className << "::Initialise:\n";
612 std::cerr << " Node number " << inode << " out of range\n";
613 std::cerr << " on potential file " << prnsol << " (line " << il
614 << ").\n";
615 ok = false;
616 } else {
617 nodes[inode - 1].v = volt;
618 nread++;
619 }
620 }
621 // Close the file
622 fprnsol.close();
623 // Tell how many lines read
624 std::cout << m_className << "::Initialise:\n";
625 std::cout << " Read " << nread << " potentials from file " << prnsol
626 << ".\n";
627 // Check number of nodes
628 if (nread != nNodes) {
629 std::cerr << m_className << "::Initialise:\n";
630 std::cerr << " Number of nodes read (" << nread << ") on potential file "
631 << prnsol << "\n";
632 std::cerr << " does not match the node list (" << nNodes << ").\n";
633 ok = false;
634 }
635
636 // Set the ready flag
637 if (ok) {
638 ready = true;
639 } else {
640 std::cerr << m_className << "::Initialise:\n";
641 std::cerr
642 << " Field map could not be read and can not be interpolated.\n";
643 return false;
644 }
645
646 // Remove weighting fields (if any).
647 wfields.clear();
648 wfieldsOk.clear();
649 nWeightingFields = 0;
650
651 // Establish the ranges
652 SetRange();
653 UpdatePeriodicity();
654 return true;
655}
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
Garfield::ComponentFieldMap::wfields
std::vector< std::string > wfields
Definition: ComponentFieldMap.hh:127

◆ IsInBoundingBox()

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

Reimplemented from Garfield::ComponentFieldMap.

Definition at line 31 of file ComponentAnsys123.hh.

31 {
32 return x >= xMinBoundingBox && x <= xMaxBoundingBox &&
33 y >= yMinBoundingBox && y <= yMaxBoundingBox &&
34 z >= zMinBoundingBox && y <= zMaxBoundingBox;
35 }

◆ SetWeightingField()

bool Garfield::ComponentAnsys123::SetWeightingField ( std::string  prnsol,
std::string  label 
)

Definition at line 657 of file ComponentAnsys123.cc.

658 {
659
660 if (!ready) {
661 std::cerr << m_className << "::SetWeightingField:\n";
662 std::cerr << " No valid field map is present.\n";
663 std::cerr << " Weighting field cannot be added.\n";
664 return false;
665 }
666
667 // Open the voltage list.
668 std::ifstream fprnsol;
669 fprnsol.open(prnsol.c_str(), std::ios::in);
670 if (fprnsol.fail()) {
671 std::cerr << m_className << "::SetWeightingField:\n";
672 std::cerr << " Could not open potential file " << prnsol
673 << " for reading.\n";
674 std::cerr << " The file perhaps does not exist.\n";
675 return false;
676 }
677
678 // Check if a weighting field with the same label already exists.
679 int iw = nWeightingFields;
680 for (int i = nWeightingFields; i--;) {
681 if (wfields[i] == label) {
682 iw = i;
683 break;
684 }
685 }
686 if (iw == nWeightingFields) {
687 ++nWeightingFields;
688 wfields.resize(nWeightingFields);
689 wfieldsOk.resize(nWeightingFields);
690 for (int j = nNodes; j--;) {
691 nodes[j].w.resize(nWeightingFields);
692 }
693 } else {
694 std::cout << m_className << "::SetWeightingField:\n";
695 std::cout << " Replacing existing weighting field " << label << ".\n";
696 }
697 wfields[iw] = label;
698 wfieldsOk[iw] = false;
699
700 // Buffer for reading
701 const int size = 100;
702 char line[size];
703
704 bool ok = true;
705 // Read the voltage list.
706 int il = 0;
707 int nread = 0;
708 bool readerror = false;
709
710 while (fprnsol.getline(line, size, '\n')) {
711 il++;
712 // Skip page feed
713 if (strcmp(line, "1") == 0) {
714 fprnsol.getline(line, size, '\n');
715 il++;
716 fprnsol.getline(line, size, '\n');
717 il++;
718 fprnsol.getline(line, size, '\n');
719 il++;
720 fprnsol.getline(line, size, '\n');
721 il++;
722 fprnsol.getline(line, size, '\n');
723 il++;
724 continue;
725 }
726 // Split the line in tokens.
727 char* token = NULL;
728 token = strtok(line, " ");
729 // Skip blank lines and headers.
730 if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
731 int(token[0]) == 10 || int(token[0]) == 13 ||
732 strcmp(token, "PRINT") == 0 || strcmp(token, "*****") == 0 ||
733 strcmp(token, "LOAD") == 0 || strcmp(token, "TIME=") == 0 ||
734 strcmp(token, "MAXIMUM") == 0 || strcmp(token, "VALUE") == 0 ||
735 strcmp(token, "NODE") == 0)
736 continue;
737 // Read the element.
738 int inode = ReadInteger(token, -1, readerror);
739 token = strtok(NULL, " ");
740 double volt = ReadDouble(token, -1, readerror);
741 // Check the syntax.
742 if (readerror) {
743 std::cerr << m_className << "::SetWeightingField:\n";
744 std::cerr << " Error reading file " << prnsol.c_str() << " (line "
745 << il << ").\n";
746 fprnsol.close();
747 return false;
748 }
749 // Check node number and store if OK.
750 if (inode < 1 || inode > nNodes) {
751 std::cerr << m_className << "::SetWeightingField:\n";
752 std::cerr << " Node number " << inode << " out of range\n";
753 std::cerr << " on potential file " << prnsol.c_str() << " (line " << il
754 << ").\n";
755 ok = false;
756 } else {
757 nodes[inode - 1].w[iw] = volt;
758 nread++;
759 }
760 }
761 // Close the file.
762 fprnsol.close();
763
764 std::cout << m_className << "::SetWeightingField:\n";
765 std::cout << " Read " << nread << " potentials from file "
766 << prnsol.c_str() << ".\n";
767 // Check the number of nodes.
768 if (nread != nNodes) {
769 std::cerr << m_className << "::SetWeightingField:\n";
770 std::cerr << " Number of nodes read (" << nread << ") "
771 << " on potential file " << prnsol.c_str() << "\n";
772 std::cerr << " does not match the node list (" << nNodes << ").\n";
773 ok = false;
774 }
775
776 // Set the ready flag.
777 wfieldsOk[iw] = ok;
778 if (!ok) {
779 std::cerr << m_className << "::SetWeightingField:\n";
780 std::cerr << " Field map could not be read "
781 << "and cannot be interpolated.\n";
782 return false;
783 }
784 return true;
785}

◆ UpdatePeriodicity()

void Garfield::ComponentAnsys123::UpdatePeriodicity ( )
inlineprotectedvirtual

Implements Garfield::ComponentFieldMap.

Definition at line 46 of file ComponentAnsys123.hh.

46{ UpdatePeriodicityCommon(); }

Referenced by Initialise().

◆ WeightingField()

void Garfield::ComponentAnsys123::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 939 of file ComponentAnsys123.cc.

941 {
942
943 // Initial values
944 wx = wy = wz = 0;
945
946 // Do not proceed if not properly initialised.
947 if (!ready) return;
948
949 // Look for the label.
950 int iw = 0;
951 bool found = false;
952 for (int i = nWeightingFields; i--;) {
953 if (wfields[i] == label) {
954 iw = i;
955 found = true;
956 break;
957 }
958 }
959
960 // Do not proceed if the requested weighting field does not exist.
961 if (!found) return;
962 // Check if the weighting field is properly initialised.
963 if (!wfieldsOk[iw]) return;
964
965 // Copy the coordinates.
966 double x = xin, y = yin, z = zin;
967
968 // Map the coordinates onto field map coordinates
969 bool xmirrored, ymirrored, zmirrored;
970 double rcoordinate, rotation;
971 MapCoordinates(x, y, z, xmirrored, ymirrored, zmirrored, rcoordinate,
972 rotation);
973
974 if (warning) {
975 std::cerr << m_className << "::WeightingField:\n";
976 std::cerr << " Warnings have been issued for this field map.\n";
977 }
978
979 // Find the element that contains this point.
980 double t1, t2, t3, t4, jac[4][4], det;
981 int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
982 // Check if the point is in the mesh.
983 if (imap < 0) return;
984
985 if (debug) {
986 std::cout << m_className << "::WeightingField:\n";
987 std::cout << " Global: (" << x << ", " << y << ", " << z << ")\n";
988 std::cout << " Local: (" << t1 << ", " << t2 << ", " << t3 << ", " << t4
989 << ") in element " << imap << "\n",
990 std::cout << " Node x y z "
991 " V\n";
992 for (int i = 0; i < 10; i++) {
993 printf(" %-5d %12g %12g %12g %12g\n", elements[imap].emap[i],
994 nodes[elements[imap].emap[i]].x, nodes[elements[imap].emap[i]].y,
995 nodes[elements[imap].emap[i]].z,
996 nodes[elements[imap].emap[i]].w[iw]);
997 }
998 }
999
1000 // Tetrahedral field
1001 wx = -(nodes[elements[imap].emap[0]].w[iw] * (4 * t1 - 1) * jac[0][1] +
1002 nodes[elements[imap].emap[1]].w[iw] * (4 * t2 - 1) * jac[1][1] +
1003 nodes[elements[imap].emap[2]].w[iw] * (4 * t3 - 1) * jac[2][1] +
1004 nodes[elements[imap].emap[3]].w[iw] * (4 * t4 - 1) * jac[3][1] +
1005 nodes[elements[imap].emap[4]].w[iw] *
1006 (4 * t2 * jac[0][1] + 4 * t1 * jac[1][1]) +
1007 nodes[elements[imap].emap[5]].w[iw] *
1008 (4 * t3 * jac[0][1] + 4 * t1 * jac[2][1]) +
1009 nodes[elements[imap].emap[6]].w[iw] *
1010 (4 * t4 * jac[0][1] + 4 * t1 * jac[3][1]) +
1011 nodes[elements[imap].emap[7]].w[iw] *
1012 (4 * t3 * jac[1][1] + 4 * t2 * jac[2][1]) +
1013 nodes[elements[imap].emap[8]].w[iw] *
1014 (4 * t4 * jac[1][1] + 4 * t2 * jac[3][1]) +
1015 nodes[elements[imap].emap[9]].w[iw] *
1016 (4 * t4 * jac[2][1] + 4 * t3 * jac[3][1])) /
1017 det;
1018
1019 wy = -(nodes[elements[imap].emap[0]].w[iw] * (4 * t1 - 1) * jac[0][2] +
1020 nodes[elements[imap].emap[1]].w[iw] * (4 * t2 - 1) * jac[1][2] +
1021 nodes[elements[imap].emap[2]].w[iw] * (4 * t3 - 1) * jac[2][2] +
1022 nodes[elements[imap].emap[3]].w[iw] * (4 * t4 - 1) * jac[3][2] +
1023 nodes[elements[imap].emap[4]].w[iw] *
1024 (4 * t2 * jac[0][2] + 4 * t1 * jac[1][2]) +
1025 nodes[elements[imap].emap[5]].w[iw] *
1026 (4 * t3 * jac[0][2] + 4 * t1 * jac[2][2]) +
1027 nodes[elements[imap].emap[6]].w[iw] *
1028 (4 * t4 * jac[0][2] + 4 * t1 * jac[3][2]) +
1029 nodes[elements[imap].emap[7]].w[iw] *
1030 (4 * t3 * jac[1][2] + 4 * t2 * jac[2][2]) +
1031 nodes[elements[imap].emap[8]].w[iw] *
1032 (4 * t4 * jac[1][2] + 4 * t2 * jac[3][2]) +
1033 nodes[elements[imap].emap[9]].w[iw] *
1034 (4 * t4 * jac[2][2] + 4 * t3 * jac[3][2])) /
1035 det;
1036
1037 wz = -(nodes[elements[imap].emap[0]].w[iw] * (4 * t1 - 1) * jac[0][3] +
1038 nodes[elements[imap].emap[1]].w[iw] * (4 * t2 - 1) * jac[1][3] +
1039 nodes[elements[imap].emap[2]].w[iw] * (4 * t3 - 1) * jac[2][3] +
1040 nodes[elements[imap].emap[3]].w[iw] * (4 * t4 - 1) * jac[3][3] +
1041 nodes[elements[imap].emap[4]].w[iw] *
1042 (4 * t2 * jac[0][3] + 4 * t1 * jac[1][3]) +
1043 nodes[elements[imap].emap[5]].w[iw] *
1044 (4 * t3 * jac[0][3] + 4 * t1 * jac[2][3]) +
1045 nodes[elements[imap].emap[6]].w[iw] *
1046 (4 * t4 * jac[0][3] + 4 * t1 * jac[3][3]) +
1047 nodes[elements[imap].emap[7]].w[iw] *
1048 (4 * t3 * jac[1][3] + 4 * t2 * jac[2][3]) +
1049 nodes[elements[imap].emap[8]].w[iw] *
1050 (4 * t4 * jac[1][3] + 4 * t2 * jac[3][3]) +
1051 nodes[elements[imap].emap[9]].w[iw] *
1052 (4 * t4 * jac[2][3] + 4 * t3 * jac[3][3])) /
1053 det;
1054
1055 // Transform field to global coordinates
1056 UnmapFields(wx, wy, wz, x, y, z, xmirrored, ymirrored, zmirrored, rcoordinate,
1057 rotation);
1058}

◆ WeightingPotential()

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

Implements Garfield::ComponentFieldMap.

Definition at line 1060 of file ComponentAnsys123.cc.

1062 {
1063
1064 // Do not proceed if not properly initialised.
1065 if (!ready) return 0.;
1066
1067 // Look for the label.
1068 int iw = 0;
1069 bool found = false;
1070 for (int i = nWeightingFields; i--;) {
1071 if (wfields[i] == label) {
1072 iw = i;
1073 found = true;
1074 break;
1075 }
1076 }
1077
1078 // Do not proceed if the requested weighting field does not exist.
1079 if (!found) return 0.;
1080 // Check if the weighting field is properly initialised.
1081 if (!wfieldsOk[iw]) return 0.;
1082
1083 // Copy the coordinates.
1084 double x = xin, y = yin, z = zin;
1085
1086 // Map the coordinates onto field map coordinates.
1087 bool xmirrored, ymirrored, zmirrored;
1088 double rcoordinate, rotation;
1089 MapCoordinates(x, y, z, xmirrored, ymirrored, zmirrored, rcoordinate,
1090 rotation);
1091
1092 if (warning) {
1093 std::cerr << m_className << "::WeightingPotential:\n";
1094 std::cerr << " Warnings have been issued for this field map.\n";
1095 }
1096
1097 // Find the element that contains this point.
1098 double t1, t2, t3, t4, jac[4][4], det;
1099 int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
1100 if (imap < 0) return 0.;
1101
1102 if (debug) {
1103 std::cerr << m_className << "::WeightingPotential:\n";
1104 std::cout << " Global: (" << x << ", " << y << ", "
1105 << ", " << z << ")\n";
1106 std::cout << " Local: (" << t1 << ", " << t2 << ", " << t3 << ", " << t4
1107 << ") in element " << imap << "\n";
1108 std::cout
1109 << " Node x y z V\n";
1110 for (int i = 0; i < 10; i++) {
1111 printf(" %-5d %12g %12g %12g %12g\n", elements[imap].emap[i],
1112 nodes[elements[imap].emap[i]].x, nodes[elements[imap].emap[i]].y,
1113 nodes[elements[imap].emap[i]].z, nodes[elements[imap].emap[i]].v);
1114 }
1115 }
1116
1117 // Tetrahedral field
1118 return nodes[elements[imap].emap[0]].w[iw] * t1 * (2 * t1 - 1) +
1119 nodes[elements[imap].emap[1]].w[iw] * t2 * (2 * t2 - 1) +
1120 nodes[elements[imap].emap[2]].w[iw] * t3 * (2 * t3 - 1) +
1121 nodes[elements[imap].emap[3]].w[iw] * t4 * (2 * t4 - 1) +
1122 4 * nodes[elements[imap].emap[4]].w[iw] * t1 * t2 +
1123 4 * nodes[elements[imap].emap[5]].w[iw] * t1 * t3 +
1124 4 * nodes[elements[imap].emap[6]].w[iw] * t1 * t4 +
1125 4 * nodes[elements[imap].emap[7]].w[iw] * t2 * t3 +
1126 4 * nodes[elements[imap].emap[8]].w[iw] * t2 * t4 +
1127 4 * nodes[elements[imap].emap[9]].w[iw] * t3 * t4;
1128}
The documentation for this class was generated from the following files: