Garfield++ v1r0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Garfield::ComponentAnsys121 Class Reference

#include <ComponentAnsys121.hh>

+ Inheritance diagram for Garfield::ComponentAnsys121:

Public Member Functions

 ComponentAnsys121 ()
 
 ~ComponentAnsys121 ()
 
MediumGetMedium (const double &x, const double &y, const double &z)
 
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="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)
 
bool IsInBoundingBox (const double x, const double y, const double z)
 
void SetRangeZ (const double zmin, const double zmax)
 
- 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)
 
MediumGetMedium (const unsigned int &i) const
 
MediumGetMedium (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 MediumGetMedium (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< elementelements
 
int lastElement
 
bool cacheElemBoundingBoxes
 
int nNodes
 
std::vector< nodenodes
 
int nMaterials
 
std::vector< materialmaterials
 
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
 
GeometryBasetheGeometry
 
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 ComponentAnsys121.hh.

Constructor & Destructor Documentation

◆ ComponentAnsys121()

Garfield::ComponentAnsys121::ComponentAnsys121 ( )

Definition at line 12 of file ComponentAnsys121.cc.

13
14 m_className = "ComponentAnsys121";
15 ready = false;
16 // Default bounding box
17 is3d = false;
18 zMinBoundingBox = -50;
19 zMaxBoundingBox = 50;
20}

◆ ~ComponentAnsys121()

Garfield::ComponentAnsys121::~ComponentAnsys121 ( )
inline

Definition at line 14 of file ComponentAnsys121.hh.

14{}

Member Function Documentation

◆ ElectricField() [1/2]

void Garfield::ComponentAnsys121::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 717 of file ComponentAnsys121.cc.

720 {
721
722 // Copy the coordinates
723 double x = xin, y = yin, z = 0.;
724
725 // Map the coordinates onto field map coordinates
726 bool xmirrored, ymirrored, zmirrored;
727 double rcoordinate, rotation;
728 MapCoordinates(x, y, z, xmirrored, ymirrored, zmirrored, rcoordinate,
729 rotation);
730
731 // Initial values
732 ex = ey = ez = volt = 0;
733 m = 0;
734 status = 0;
735
736 // Do not proceed if not properly initialised.
737 if (!ready) {
738 status = -10;
739 std::cerr << m_className << "::ElectricField:\n";
740 std::cerr << " Field map not available for interpolation.\n";
741 return;
742 }
743
744 if (warning) {
745 std::cerr << m_className << "::ElectricField:\n";
746 std::cerr << " Warnings have been issued for this field map.\n";
747 }
748
749 if (zin < zMinBoundingBox || zin > zMaxBoundingBox) {
750 status = -5;
751 return;
752 }
753
754 // Find the element that contains this point
755 double t1, t2, t3, t4, jac[4][4], det;
756 int imap = FindElement5(x, y, z, t1, t2, t3, t4, jac, det);
757 if (imap < 0) {
758 if (debug) {
759 std::cout << m_className << "::ElectricField:\n";
760 std::cout << " Point (" << x << ", " << y << ") not in the mesh.\n";
761 }
762 status = -6;
763 return;
764 }
765
766 if (debug) {
767 std::cout << m_className << "::ElectricField:\n";
768 std::cout << " Global: (" << x << ", " << y << "),\n";
769 std::cout << " Local: (" << t1 << ", " << t2 << ", " << t3 << ", " << t4
770 << ") in element " << imap
771 << " (degenerate: " << elements[imap].degenerate << ")\n";
772 std::cout
773 << " Node x y V\n";
774 for (int i = 0; i < 8; i++) {
775 printf(" %-5d %12g %12g %12g\n", elements[imap].emap[i],
776 nodes[elements[imap].emap[i]].x, nodes[elements[imap].emap[i]].y,
777 nodes[elements[imap].emap[i]].v);
778 }
779 }
780
781 // Calculate quadrilateral field, which can degenerate to a triangular field
782 if (elements[imap].degenerate) {
783 volt = nodes[elements[imap].emap[0]].v * t1 * (2 * t1 - 1) +
784 nodes[elements[imap].emap[1]].v * t2 * (2 * t2 - 1) +
785 nodes[elements[imap].emap[2]].v * t3 * (2 * t3 - 1) +
786 4 * nodes[elements[imap].emap[3]].v * t1 * t2 +
787 4 * nodes[elements[imap].emap[4]].v * t1 * t3 +
788 4 * nodes[elements[imap].emap[5]].v * t2 * t3;
789 ex = -(nodes[elements[imap].emap[0]].v * (4 * t1 - 1) * jac[0][1] +
790 nodes[elements[imap].emap[1]].v * (4 * t2 - 1) * jac[1][1] +
791 nodes[elements[imap].emap[2]].v * (4 * t3 - 1) * jac[2][1] +
792 nodes[elements[imap].emap[3]].v *
793 (4 * t2 * jac[0][1] + 4 * t1 * jac[1][1]) +
794 nodes[elements[imap].emap[4]].v *
795 (4 * t3 * jac[0][1] + 4 * t1 * jac[2][1]) +
796 nodes[elements[imap].emap[5]].v *
797 (4 * t3 * jac[1][1] + 4 * t2 * jac[2][1])) /
798 det;
799 ey = -(nodes[elements[imap].emap[0]].v * (4 * t1 - 1) * jac[0][2] +
800 nodes[elements[imap].emap[1]].v * (4 * t2 - 1) * jac[1][2] +
801 nodes[elements[imap].emap[2]].v * (4 * t3 - 1) * jac[2][2] +
802 nodes[elements[imap].emap[3]].v *
803 (4 * t2 * jac[0][2] + 4 * t1 * jac[1][2]) +
804 nodes[elements[imap].emap[4]].v *
805 (4 * t3 * jac[0][2] + 4 * t1 * jac[2][2]) +
806 nodes[elements[imap].emap[5]].v *
807 (4 * t3 * jac[1][2] + 4 * t2 * jac[2][2])) /
808 det;
809 } else {
810 volt =
811 -nodes[elements[imap].emap[0]].v * (1 - t1) * (1 - t2) * (1 + t1 + t2) /
812 4 -
813 nodes[elements[imap].emap[1]].v * (1 + t1) * (1 - t2) * (1 - t1 + t2) /
814 4 -
815 nodes[elements[imap].emap[2]].v * (1 + t1) * (1 + t2) * (1 - t1 - t2) /
816 4 -
817 nodes[elements[imap].emap[3]].v * (1 - t1) * (1 + t2) * (1 + t1 - t2) /
818 4 +
819 nodes[elements[imap].emap[4]].v * (1 - t1) * (1 + t1) * (1 - t2) / 2 +
820 nodes[elements[imap].emap[5]].v * (1 + t1) * (1 + t2) * (1 - t2) / 2 +
821 nodes[elements[imap].emap[6]].v * (1 - t1) * (1 + t1) * (1 + t2) / 2 +
822 nodes[elements[imap].emap[7]].v * (1 - t1) * (1 + t2) * (1 - t2) / 2;
823 ex = -(nodes[elements[imap].emap[0]].v *
824 ((1 - t2) * (2 * t1 + t2) * jac[0][0] +
825 (1 - t1) * (t1 + 2 * t2) * jac[1][0]) /
826 4 +
827 nodes[elements[imap].emap[1]].v *
828 ((1 - t2) * (2 * t1 - t2) * jac[0][0] -
829 (1 + t1) * (t1 - 2 * t2) * jac[1][0]) /
830 4 +
831 nodes[elements[imap].emap[2]].v *
832 ((1 + t2) * (2 * t1 + t2) * jac[0][0] +
833 (1 + t1) * (t1 + 2 * t2) * jac[1][0]) /
834 4 +
835 nodes[elements[imap].emap[3]].v *
836 ((1 + t2) * (2 * t1 - t2) * jac[0][0] -
837 (1 - t1) * (t1 - 2 * t2) * jac[1][0]) /
838 4 +
839 nodes[elements[imap].emap[4]].v *
840 (t1 * (t2 - 1) * jac[0][0] +
841 (t1 - 1) * (t1 + 1) * jac[1][0] / 2) +
842 nodes[elements[imap].emap[5]].v *
843 ((1 - t2) * (1 + t2) * jac[0][0] / 2 -
844 (1 + t1) * t2 * jac[1][0]) +
845 nodes[elements[imap].emap[6]].v *
846 (-t1 * (1 + t2) * jac[0][0] +
847 (1 - t1) * (1 + t1) * jac[1][0] / 2) +
848 nodes[elements[imap].emap[7]].v *
849 ((t2 - 1) * (t2 + 1) * jac[0][0] / 2 +
850 (t1 - 1) * t2 * jac[1][0])) /
851 det;
852 ey = -(nodes[elements[imap].emap[0]].v *
853 ((1 - t2) * (2 * t1 + t2) * jac[0][1] +
854 (1 - t1) * (t1 + 2 * t2) * jac[1][1]) /
855 4 +
856 nodes[elements[imap].emap[1]].v *
857 ((1 - t2) * (2 * t1 - t2) * jac[0][1] -
858 (1 + t1) * (t1 - 2 * t2) * jac[1][1]) /
859 4 +
860 nodes[elements[imap].emap[2]].v *
861 ((1 + t2) * (2 * t1 + t2) * jac[0][1] +
862 (1 + t1) * (t1 + 2 * t2) * jac[1][1]) /
863 4 +
864 nodes[elements[imap].emap[3]].v *
865 ((1 + t2) * (2 * t1 - t2) * jac[0][1] -
866 (1 - t1) * (t1 - 2 * t2) * jac[1][1]) /
867 4 +
868 nodes[elements[imap].emap[4]].v *
869 (t1 * (t2 - 1) * jac[0][1] +
870 (t1 - 1) * (t1 + 1) * jac[1][1] / 2) +
871 nodes[elements[imap].emap[5]].v *
872 ((1 - t2) * (1 + t2) * jac[0][1] / 2 -
873 (1 + t1) * t2 * jac[1][1]) +
874 nodes[elements[imap].emap[6]].v *
875 (-t1 * (1 + t2) * jac[0][1] +
876 (1 - t1) * (1 + t1) * jac[1][1] / 2) +
877 nodes[elements[imap].emap[7]].v *
878 ((t2 - 1) * (t2 + 1) * jac[0][1] / 2 +
879 (t1 - 1) * t2 * jac[1][1])) /
880 det;
881 }
882
883 // Transform field to global coordinates
884 UnmapFields(ex, ey, ez, x, y, z, xmirrored, ymirrored, zmirrored, rcoordinate,
885 rotation);
886
887 // Drift medium?
888 if (debug) {
889 std::cout << m_className << "::ElectricField:\n";
890 std::cout << " Material " << elements[imap].matmap << ", drift flag "
891 << materials[elements[imap].matmap].driftmedium << ".\n";
892 }
893 m = materials[elements[imap].matmap].medium;
894 status = -5;
895 if (materials[elements[imap].matmap].driftmedium) {
896 if (m != 0) {
897 if (m->IsDriftable()) status = 0;
898 }
899 }
900}
void MapCoordinates(double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
std::vector< material > materials
void UnmapFields(double &ex, double &ey, double &ez, double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation)
std::vector< element > elements
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)

◆ ElectricField() [2/2]

void Garfield::ComponentAnsys121::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 709 of file ComponentAnsys121.cc.

711 {
712
713 double v;
714 ElectricField(x, y, z, ex, ey, ez, v, m, status);
715}
void ElectricField(const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status)

Referenced by ElectricField().

◆ GetAspectRatio()

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

Implements Garfield::ComponentFieldMap.

Definition at line 1236 of file ComponentAnsys121.cc.

1237 {
1238
1239 if (i < 0 || i >= nElements) {
1240 dmin = dmax = 0.;
1241 return;
1242 }
1243
1244 const int np = 8;
1245 // Loop over all pairs of vertices.
1246 for (int j = 0; j < np - 1; ++j) {
1247 for (int k = j + 1; k < np; ++k) {
1248 // Compute distance.
1249 const double dist = sqrt(
1250 pow(nodes[elements[i].emap[j]].x - nodes[elements[i].emap[k]].x, 2) +
1251 pow(nodes[elements[i].emap[j]].y - nodes[elements[i].emap[k]].y, 2));
1252 if (k == 1) {
1253 dmin = dist;
1254 dmax = dist;
1255 } else {
1256 if (dist < dmin) dmin = dist;
1257 if (dist > dmax) dmax = dist;
1258 }
1259 }
1260 }
1261}
DoubleAc pow(const DoubleAc &f, double p)
Definition: DoubleAc.cpp:336
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:313

◆ GetElementVolume()

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

Implements Garfield::ComponentFieldMap.

Definition at line 1219 of file ComponentAnsys121.cc.

1219 {
1220
1221 if (i < 0 || i >= nElements) return 0.;
1222 const double surf =
1223 (fabs((nodes[elements[i].emap[1]].x - nodes[elements[i].emap[0]].x) *
1224 (nodes[elements[i].emap[2]].y - nodes[elements[i].emap[0]].y) -
1225 (nodes[elements[i].emap[2]].x - nodes[elements[i].emap[0]].x) *
1226 (nodes[elements[i].emap[1]].y - nodes[elements[i].emap[0]].y)) +
1227 fabs(
1228 (nodes[elements[i].emap[3]].x - nodes[elements[i].emap[0]].x) *
1229 (nodes[elements[i].emap[2]].y - nodes[elements[i].emap[0]].y) -
1230 (nodes[elements[i].emap[2]].x - nodes[elements[i].emap[0]].x) *
1231 (nodes[elements[i].emap[3]].y - nodes[elements[i].emap[0]].y))) /
1232 2.;
1233 return surf;
1234}
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:616

◆ GetMedium()

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

Implements Garfield::ComponentFieldMap.

Definition at line 1137 of file ComponentAnsys121.cc.

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

◆ Initialise()

bool Garfield::ComponentAnsys121::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 22 of file ComponentAnsys121.cc.

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

◆ IsInBoundingBox()

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

◆ SetRangeZ()

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

Definition at line 1202 of file ComponentAnsys121.cc.

1202 {
1203
1204 if (fabs(zmax - zmin) <= 0.) {
1205 std::cerr << m_className << "::SetRangeZ:\n";
1206 std::cerr << " Zero range is not permitted.\n";
1207 return;
1208 }
1209 zMinBoundingBox = mapzmin = std::min(zmin, zmax);
1210 zMaxBoundingBox = mapzmax = std::max(zmin, zmax);
1211}

◆ SetWeightingField()

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

Definition at line 594 of file ComponentAnsys121.cc.

595 {
596
597 if (!ready) {
598 std::cerr << m_className << "::SetWeightingField:\n";
599 std::cerr << " No valid field map is present.\n";
600 std::cerr << " Weighting field cannot be added.\n";
601 return false;
602 }
603
604 // Open the voltage list.
605 std::ifstream fprnsol;
606 fprnsol.open(prnsol.c_str(), std::ios::in);
607 if (fprnsol.fail()) {
608 std::cerr << m_className << "::SetWeightingField:\n";
609 std::cerr << " Could not open potential file " << prnsol
610 << " for reading.\n";
611 std::cerr << " The file perhaps does not exist.\n";
612 return false;
613 }
614
615 // Check if a weighting field with the same label already exists.
616 int iw = nWeightingFields;
617 for (int i = nWeightingFields; i--;) {
618 if (wfields[i] == label) {
619 iw = i;
620 break;
621 }
622 }
623 if (iw == nWeightingFields) {
627 for (int j = nNodes; j--;) {
628 nodes[j].w.resize(nWeightingFields);
629 }
630 } else {
631 std::cout << m_className << "::SetWeightingField:\n";
632 std::cout << " Replacing existing weighting field " << label << ".\n";
633 }
634 wfields[iw] = label;
635 wfieldsOk[iw] = false;
636
637 // Buffer for reading
638 const int size = 100;
639 char line[size];
640
641 bool ok = true;
642 // Read the voltage list.
643 int il = 0;
644 int nread = 0;
645 bool readerror = false;
646 while (fprnsol.getline(line, size, '\n')) {
647 il++;
648 // Split the line in tokens.
649 char* token = NULL;
650 token = strtok(line, " ");
651 // Skip blank lines and headers.
652 if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
653 int(token[0]) == 10 || int(token[0]) == 13 ||
654 strcmp(token, "PRINT") == 0 || strcmp(token, "*****") == 0 ||
655 strcmp(token, "LOAD") == 0 || strcmp(token, "TIME=") == 0 ||
656 strcmp(token, "MAXIMUM") == 0 || strcmp(token, "VALUE") == 0 ||
657 strcmp(token, "NODE") == 0)
658 continue;
659 // Read the element.
660 int inode = ReadInteger(token, -1, readerror);
661 token = strtok(NULL, " ");
662 double volt = ReadDouble(token, -1, readerror);
663 // Check the syntax.
664 if (readerror) {
665 std::cerr << m_className << "::SetWeightingField:\n";
666 std::cerr << " Error reading file " << prnsol << " (line " << il
667 << ").\n";
668 fprnsol.close();
669 return false;
670 }
671 // Check node number and store if OK.
672 if (inode < 1 || inode > nNodes) {
673 std::cerr << m_className << "::SetWeightingField:\n";
674 std::cerr << " Node number " << inode << " out of range\n";
675 std::cerr << " on potential file " << prnsol << " (line " << il
676 << ").\n";
677 ok = false;
678 } else {
679 nodes[inode - 1].w[iw] = volt;
680 nread++;
681 }
682 }
683 // Close the file.
684 fprnsol.close();
685 // Tell how many lines read.
686 std::cout << m_className << "::SetWeightingField:\n";
687 std::cout << " Read " << nread << " potentials from file " << prnsol
688 << ".\n";
689 // Check the number of nodes.
690 if (nread != nNodes) {
691 std::cerr << m_className << "::SetWeightingField:\n";
692 std::cerr << " Number of nodes read (" << nread << ") "
693 << " on potential file " << prnsol << "\n";
694 std::cerr << " does not match the node list (" << nNodes << ").\n";
695 ok = false;
696 }
697
698 // Set the ready flag.
699 wfieldsOk[iw] = ok;
700 if (!ok) {
701 std::cerr << m_className << "::SetWeightingField:\n";
702 std::cerr << " Field map could not be read "
703 << "and cannot be interpolated.\n";
704 return false;
705 }
706 return true;
707}

◆ UpdatePeriodicity()

void Garfield::ComponentAnsys121::UpdatePeriodicity ( )
protectedvirtual

Implements Garfield::ComponentFieldMap.

Definition at line 1213 of file ComponentAnsys121.cc.

Referenced by Initialise().

◆ WeightingField()

void Garfield::ComponentAnsys121::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 902 of file ComponentAnsys121.cc.

904 {
905
906 // Initial values
907 wx = wy = wz = 0;
908
909 // Do not proceed if not properly initialised.
910 if (!ready) return;
911
912 // Look for the label.
913 int iw = 0;
914 bool found = false;
915 for (int i = nWeightingFields; i--;) {
916 if (wfields[i] == label) {
917 iw = i;
918 found = true;
919 break;
920 }
921 }
922
923 // Do not proceed if the requested weighting field does not exist.
924 if (!found) return;
925 // Check if the weighting field is properly initialised.
926 if (!wfieldsOk[iw]) return;
927
928 // Copy the coordinates.
929 double x = xin, y = yin, z = zin;
930
931 // Map the coordinates onto field map coordinates
932 bool xmirrored, ymirrored, zmirrored;
933 double rcoordinate, rotation;
934 MapCoordinates(x, y, z, xmirrored, ymirrored, zmirrored, rcoordinate,
935 rotation);
936
937 if (warning) {
938 std::cerr << m_className << "::WeightingField:\n";
939 std::cerr << " Warnings have been issued for this field map.\n";
940 }
941
942 // Find the element that contains this point.
943 double t1, t2, t3, t4, jac[4][4], det;
944 int imap = FindElement5(x, y, z, t1, t2, t3, t4, jac, det);
945 // Check if the point is in the mesh.
946 if (imap < 0) return;
947
948 if (debug) {
949 std::cout << m_className << "::WeightingField:\n";
950 std::cout << " Global: (" << x << ", " << y << "),\n";
951 std::cout << " Local: (" << t1 << ", " << t2 << ", " << t3 << ", " << t4
952 << ") in element " << imap
953 << " (degenerate: " << elements[imap].degenerate << ")\n";
954 std::cout
955 << " Node x y V\n";
956 for (int i = 0; i < 8; i++) {
957 printf(" %-5d %12g %12g %12g\n", elements[imap].emap[i],
958 nodes[elements[imap].emap[i]].x, nodes[elements[imap].emap[i]].y,
959 nodes[elements[imap].emap[i]].w[iw]);
960 }
961 }
962
963 // Calculate quadrilateral field, which can degenerate to a triangular field
964 if (elements[imap].degenerate) {
965 wx = -(nodes[elements[imap].emap[0]].w[iw] * (4 * t1 - 1) * jac[0][1] +
966 nodes[elements[imap].emap[1]].w[iw] * (4 * t2 - 1) * jac[1][1] +
967 nodes[elements[imap].emap[2]].w[iw] * (4 * t3 - 1) * jac[2][1] +
968 nodes[elements[imap].emap[3]].w[iw] *
969 (4 * t2 * jac[0][1] + 4 * t1 * jac[1][1]) +
970 nodes[elements[imap].emap[4]].w[iw] *
971 (4 * t3 * jac[0][1] + 4 * t1 * jac[2][1]) +
972 nodes[elements[imap].emap[5]].w[iw] *
973 (4 * t3 * jac[1][1] + 4 * t2 * jac[2][1])) /
974 det;
975 wy = -(nodes[elements[imap].emap[0]].w[iw] * (4 * t1 - 1) * jac[0][2] +
976 nodes[elements[imap].emap[1]].w[iw] * (4 * t2 - 1) * jac[1][2] +
977 nodes[elements[imap].emap[2]].w[iw] * (4 * t3 - 1) * jac[2][2] +
978 nodes[elements[imap].emap[3]].w[iw] *
979 (4 * t2 * jac[0][2] + 4 * t1 * jac[1][2]) +
980 nodes[elements[imap].emap[4]].w[iw] *
981 (4 * t3 * jac[0][2] + 4 * t1 * jac[2][2]) +
982 nodes[elements[imap].emap[5]].w[iw] *
983 (4 * t3 * jac[1][2] + 4 * t2 * jac[2][2])) /
984 det;
985 } else {
986 wx = -(nodes[elements[imap].emap[0]].w[iw] *
987 ((1 - t2) * (2 * t1 + t2) * jac[0][0] +
988 (1 - t1) * (t1 + 2 * t2) * jac[1][0]) /
989 4 +
990 nodes[elements[imap].emap[1]].w[iw] *
991 ((1 - t2) * (2 * t1 - t2) * jac[0][0] -
992 (1 + t1) * (t1 - 2 * t2) * jac[1][0]) /
993 4 +
994 nodes[elements[imap].emap[2]].w[iw] *
995 ((1 + t2) * (2 * t1 + t2) * jac[0][0] +
996 (1 + t1) * (t1 + 2 * t2) * jac[1][0]) /
997 4 +
998 nodes[elements[imap].emap[3]].w[iw] *
999 ((1 + t2) * (2 * t1 - t2) * jac[0][0] -
1000 (1 - t1) * (t1 - 2 * t2) * jac[1][0]) /
1001 4 +
1002 nodes[elements[imap].emap[4]].w[iw] *
1003 (t1 * (t2 - 1) * jac[0][0] +
1004 (t1 - 1) * (t1 + 1) * jac[1][0] / 2) +
1005 nodes[elements[imap].emap[5]].w[iw] *
1006 ((1 - t2) * (1 + t2) * jac[0][0] / 2 -
1007 (1 + t1) * t2 * jac[1][0]) +
1008 nodes[elements[imap].emap[6]].w[iw] *
1009 (-t1 * (1 + t2) * jac[0][0] +
1010 (1 - t1) * (1 + t1) * jac[1][0] / 2) +
1011 nodes[elements[imap].emap[7]].w[iw] *
1012 ((t2 - 1) * (1 + t2) * jac[0][0] / 2 +
1013 (t1 - 1) * t2 * jac[1][0])) /
1014 det;
1015 wy = -(nodes[elements[imap].emap[0]].w[iw] *
1016 ((1 - t2) * (2 * t1 + t2) * jac[0][1] +
1017 (1 - t1) * (t1 + 2 * t2) * jac[1][1]) /
1018 4 +
1019 nodes[elements[imap].emap[1]].w[iw] *
1020 ((1 - t2) * (2 * t1 - t2) * jac[0][1] -
1021 (1 + t1) * (t1 - 2 * t2) * jac[1][1]) /
1022 4 +
1023 nodes[elements[imap].emap[2]].w[iw] *
1024 ((1 + t2) * (2 * t1 + t2) * jac[0][1] +
1025 (1 + t1) * (t1 + 2 * t2) * jac[1][1]) /
1026 4 +
1027 nodes[elements[imap].emap[3]].w[iw] *
1028 ((1 + t2) * (2 * t1 - t2) * jac[0][1] -
1029 (1 - t1) * (t1 - 2 * t2) * jac[1][1]) /
1030 4 +
1031 nodes[elements[imap].emap[4]].w[iw] *
1032 (t1 * (t2 - 1) * jac[0][1] +
1033 (t1 - 1) * (t1 + 1) * jac[1][1] / 2) +
1034 nodes[elements[imap].emap[5]].w[iw] *
1035 ((1 - t2) * (1 + t2) * jac[0][1] / 2 -
1036 (1 + t1) * t2 * jac[1][1]) +
1037 nodes[elements[imap].emap[6]].w[iw] *
1038 (-t1 * (1 + t2) * jac[0][1] +
1039 (1 - t1) * (1 + t1) * jac[1][1] / 2) +
1040 nodes[elements[imap].emap[7]].w[iw] *
1041 ((t2 - 1) * (t2 + 1) * jac[0][1] / 2 +
1042 (t1 - 1) * t2 * jac[1][1])) /
1043 det;
1044 }
1045
1046 // Transform field to global coordinates
1047 UnmapFields(wx, wy, wz, x, y, z, xmirrored, ymirrored, zmirrored, rcoordinate,
1048 rotation);
1049}

◆ WeightingPotential()

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

Implements Garfield::ComponentFieldMap.

Definition at line 1051 of file ComponentAnsys121.cc.

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

The documentation for this class was generated from the following files: