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

Component for importing and interpolating two-dimensional ANSYS field maps. More...

#include <ComponentAnsys121.hh>

+ Inheritance diagram for Garfield::ComponentAnsys121:

Public Member Functions

 ComponentAnsys121 ()
 Constructor.
 
 ~ComponentAnsys121 ()
 Destructor.
 
Medium * GetMedium (const double x, const double y, const double z)
 Get the medium at a given location (x, y, z).
 
void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status)
 
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)
 
virtual bool IsInBoundingBox (const double x, const double y, const double z) const
 
void SetRangeZ (const double zmin, const double zmax)
 
- Public Member Functions inherited from Garfield::ComponentFieldMap
 ComponentFieldMap ()
 Constructor.
 
virtual ~ComponentFieldMap ()
 Destructor.
 
virtual void SetRange ()
 Calculate x, y, z, V and angular ranges.
 
void PrintRange ()
 Show x, y, z, V and angular ranges.
 
virtual bool IsInBoundingBox (const double x, const double y, const double z) const
 
virtual bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
 Get the bounding box coordinates.
 
virtual bool GetVoltageRange (double &vmin, double &vmax)
 Calculate the voltage range [V].
 
void PrintMaterials ()
 List all currently defined materials.
 
void DriftMedium (const unsigned int imat)
 Flag a field map material as a drift medium.
 
void NotDriftMedium (const unsigned int imat)
 Flag a field map materials as a non-drift medium.
 
unsigned int GetNumberOfMaterials () const
 Return the number of materials in the field map.
 
double GetPermittivity (const unsigned int imat) const
 Return the permittivity of a field map material.
 
double GetConductivity (const unsigned int imat) const
 Return the conductivity of a field map material.
 
void SetMedium (const unsigned int imat, Medium *medium)
 Associate a field map material with a Medium class.
 
Medium * GetMedium (const unsigned int i) const
 Return the Medium associated to a field map material.
 
Medium * GetMedium (const double x, const double y, const double z)=0
 Get the medium at a given location (x, y, z).
 
unsigned int GetNumberOfMedia () const
 
int GetNumberOfElements () const
 Return the number of mesh elements.
 
bool GetElement (const unsigned int i, double &vol, double &dmin, double &dmax)
 Return the volume and aspect ratio of a mesh element.
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status)=0
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status)=0
 
virtual void WeightingField (const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label)=0
 
virtual double WeightingPotential (const double x, const double y, const double z, const std::string &label)=0
 
void EnableCheckMapIndices ()
 
void DisableCheckMapIndices ()
 
void EnableDeleteBackgroundElements ()
 
void DisableDeleteBackgroundElements ()
 
void EnableTetrahedralTreeForElementSearch (const bool on=true)
 
- Public Member Functions inherited from Garfield::ComponentBase
 ComponentBase ()
 Constructor.
 
virtual ~ComponentBase ()
 Destructor.
 
virtual void SetGeometry (GeometryBase *geo)
 Define the geometry.
 
virtual void Clear ()
 Reset.
 
virtual Medium * GetMedium (const double x, const double y, const double z)
 Get the medium at a given location (x, y, z).
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status)=0
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status)=0
 
virtual bool GetVoltageRange (double &vmin, double &vmax)=0
 Calculate the voltage range [V].
 
virtual void WeightingField (const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label)
 
virtual double WeightingPotential (const double x, const double y, const double z, const std::string &label)
 
virtual void MagneticField (const double x, const double y, const double z, double &bx, double &by, double &bz, int &status)
 
void SetMagneticField (const double bx, const double by, const double bz)
 Set a constant magnetic field.
 
virtual bool IsReady ()
 Ready for use?
 
virtual bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
 Get the bounding box coordinates.
 
virtual bool IsWireCrossed (const double x0, const double y0, const double z0, const double x1, const double y1, const double z1, double &xc, double &yc, double &zc)
 
virtual bool IsInTrapRadius (const double q0, const double x0, const double y0, const double z0, double &xw, double &yw, double &rw)
 
void EnablePeriodicityX (const bool on=true)
 Enable simple periodicity in the $x$ direction.
 
void DisablePeriodicityX ()
 
void EnablePeriodicityY (const bool on=true)
 Enable simple periodicity in the $y$ direction.
 
void DisablePeriodicityY ()
 
void EnablePeriodicityZ (const bool on=true)
 Enable simple periodicity in the $z$ direction.
 
void DisablePeriodicityZ ()
 
void EnableMirrorPeriodicityX (const bool on=true)
 Enable mirror periodicity in the $x$ direction.
 
void DisableMirrorPeriodicityX ()
 
void EnableMirrorPeriodicityY (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void DisableMirrorPeriodicityY ()
 
void EnableMirrorPeriodicityZ (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void DisableMirrorPeriodicityZ ()
 
void EnableAxialPeriodicityX (const bool on=true)
 Enable axial periodicity in the $x$ direction.
 
void DisableAxialPeriodicityX ()
 
void EnableAxialPeriodicityY (const bool on=true)
 Enable axial periodicity in the $y$ direction.
 
void DisableAxialPeriodicityY ()
 
void EnableAxialPeriodicityZ (const bool on=true)
 Enable axial periodicity in the $z$ direction.
 
void DisableAxialPeriodicityZ ()
 
void EnableRotationSymmetryX (const bool on=true)
 Enable rotation symmetry around the $x$ axis.
 
void DisableRotationSymmetryX ()
 
void EnableRotationSymmetryY (const bool on=true)
 Enable rotation symmetry around the $y$ axis.
 
void DisableRotationSymmetryY ()
 
void EnableRotationSymmetryZ (const bool on=true)
 Enable rotation symmetry around the $z$ axis.
 
void DisableRotationSymmetryZ ()
 
void EnableDebugging ()
 Switch on debugging messages.
 
void DisableDebugging ()
 Switch off debugging messages.
 
void ActivateTraps ()
 Request trapping to be taken care of by the component (for TCAD).
 
void DeactivateTraps ()
 
bool IsTrapActive ()
 
void ActivateVelocityMap ()
 Request velocity to be taken care of by the component (for TCAD).
 
void DectivateVelocityMap ()
 
bool IsVelocityActive ()
 
virtual bool ElectronAttachment (const double, const double, const double, double &eta)
 Get the electron attachment coefficient.
 
virtual bool HoleAttachment (const double, const double, const double, double &eta)
 Get the hole attachment coefficient.
 
virtual void ElectronVelocity (const double, const double, const double, double &vx, double &vy, double &vz, Medium *&, int &status)
 Get the electron drift velocity.
 
virtual void HoleVelocity (const double, const double, const double, double &vx, double &vy, double &vz, Medium *&, int &status)
 Get the hole drift velocity.
 
virtual bool GetElectronLifetime (const double, const double, const double, double &etau)
 
virtual bool GetHoleLifetime (const double, const double, const double, double &htau)
 

Protected Member Functions

void UpdatePeriodicity ()
 Verify periodicities.
 
double GetElementVolume (const unsigned int i)
 
void GetAspectRatio (const unsigned int i, double &dmin, double &dmax)
 
- Protected Member Functions inherited from Garfield::ComponentFieldMap
void Reset ()
 Geometry checks.
 
virtual void UpdatePeriodicity ()=0
 Verify periodicities.
 
void UpdatePeriodicity2d ()
 
void UpdatePeriodicityCommon ()
 
int FindElement5 (const double x, const double y, const double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det)
 Find the element for a point in curved quadratic quadrilaterals.
 
int FindElement13 (const double x, const double y, const double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det)
 Find the element for a point in curved quadratic tetrahedra.
 
int FindElementCube (const double x, const double y, const double z, double &t1, double &t2, double &t3, TMatrixD *&jac, std::vector< TMatrixD * > &dN)
 Find the element for a point in a cube.
 
void MapCoordinates (double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
 Move (xpos, ypos, zpos) to field map coordinates.
 
void UnmapFields (double &ex, double &ey, double &ez, double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
 Move (ex, ey, ez) to global coordinates.
 
int ReadInteger (char *token, int def, bool &error)
 
double ReadDouble (char *token, double def, bool &error)
 
virtual double GetElementVolume (const unsigned int i)=0
 
virtual void GetAspectRatio (const unsigned int i, double &dmin, double &dmax)=0
 
void PrintWarning (const std::string &header)
 
void PrintNotReady (const std::string &header) const
 
void PrintElement (const std::string &header, const double x, const double y, const double z, const double t1, const double t2, const double t3, const double t4, const unsigned int i, const unsigned int n, const int iw=-1) const
 
virtual void Reset ()=0
 Geometry checks.
 
virtual void UpdatePeriodicity ()=0
 Verify periodicities.
 

Additional Inherited Members

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

Detailed Description

Component for importing and interpolating two-dimensional ANSYS field maps.

Definition at line 10 of file ComponentAnsys121.hh.

Constructor & Destructor Documentation

◆ ComponentAnsys121()

Garfield::ComponentAnsys121::ComponentAnsys121 ( )

Constructor.

Definition at line 10 of file ComponentAnsys121.cc.

10 : ComponentFieldMap() {
11
12 m_className = "ComponentAnsys121";
13 // Default bounding box
14 m_is3d = false;
15 zMinBoundingBox = -50;
16 zMaxBoundingBox = 50;
17}
Garfield::ComponentBase::m_className
std::string m_className
Class name.
Definition: ComponentBase.hh:233

◆ ~ComponentAnsys121()

Garfield::ComponentAnsys121::~ComponentAnsys121 ( )
inline

Destructor.

Definition at line 16 of file ComponentAnsys121.hh.

16{}

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 xmirr, ymirr, zmirr;
727 double rcoordinate, rotation;
728 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
729
730 // Initial values
731 ex = ey = ez = volt = 0;
732 m = NULL;
733 status = 0;
734
735 // Do not proceed if not properly initialised.
736 if (!m_ready) {
737 status = -10;
738 PrintNotReady("ElectricField");
739 return;
740 }
741
742 if (m_warning) PrintWarning("ElectricField");
743
744 if (zin < zMinBoundingBox || zin > zMaxBoundingBox) {
745 status = -5;
746 return;
747 }
748
749 // Find the element that contains this point
750 double t1, t2, t3, t4, jac[4][4], det;
751 const int imap = FindElement5(x, y, z, t1, t2, t3, t4, jac, det);
752 if (imap < 0) {
753 if (m_debug) {
754 std::cout << m_className << "::ElectricField:\n";
755 std::cout << " Point (" << x << ", " << y << ") not in the mesh.\n";
756 }
757 status = -6;
758 return;
759 }
760
761 if (m_debug) {
762 PrintElement("ElectricField", x, y, z, t1, t2, t3, t4, imap, 8);
763 }
764
765 const Element& element = elements[imap];
766 const Node& n0 = nodes[element.emap[0]];
767 const Node& n1 = nodes[element.emap[1]];
768 const Node& n2 = nodes[element.emap[2]];
769 const Node& n3 = nodes[element.emap[3]];
770 const Node& n4 = nodes[element.emap[4]];
771 const Node& n5 = nodes[element.emap[5]];
772 // Calculate quadrilateral field, which can degenerate to a triangular field
773 const double invdet = 1. / det;
774 if (element.degenerate) {
775 volt = n0.v * t1 * (2 * t1 - 1) + n1.v * t2 * (2 * t2 - 1) +
776 n2.v * t3 * (2 * t3 - 1) + 4 * n3.v * t1 * t2 + 4 * n4.v * t1 * t3 +
777 4 * n5.v * t2 * t3;
778 ex = -(n0.v * (4 * t1 - 1) * jac[0][1] +
779 n1.v * (4 * t2 - 1) * jac[1][1] +
780 n2.v * (4 * t3 - 1) * jac[2][1] +
781 n3.v * (4 * t2 * jac[0][1] + 4 * t1 * jac[1][1]) +
782 n4.v * (4 * t3 * jac[0][1] + 4 * t1 * jac[2][1]) +
783 n5.v * (4 * t3 * jac[1][1] + 4 * t2 * jac[2][1])) * invdet;
784 ey = -(n0.v * (4 * t1 - 1) * jac[0][2] +
785 n1.v * (4 * t2 - 1) * jac[1][2] +
786 n2.v * (4 * t3 - 1) * jac[2][2] +
787 n3.v * (4 * t2 * jac[0][2] + 4 * t1 * jac[1][2]) +
788 n4.v * (4 * t3 * jac[0][2] + 4 * t1 * jac[2][2]) +
789 n5.v * (4 * t3 * jac[1][2] + 4 * t2 * jac[2][2])) * invdet;
790 } else {
791 const Node& n6 = nodes[element.emap[6]];
792 const Node& n7 = nodes[element.emap[7]];
793 volt = -n0.v * (1 - t1) * (1 - t2) * (1 + t1 + t2) * 0.25 -
794 n1.v * (1 + t1) * (1 - t2) * (1 - t1 + t2) * 0.25 -
795 n2.v * (1 + t1) * (1 + t2) * (1 - t1 - t2) * 0.25 -
796 n3.v * (1 - t1) * (1 + t2) * (1 + t1 - t2) * 0.25 +
797 n4.v * (1 - t1) * (1 + t1) * (1 - t2) * 0.5 +
798 n5.v * (1 + t1) * (1 + t2) * (1 - t2) * 0.5 +
799 n6.v * (1 - t1) * (1 + t1) * (1 + t2) * 0.5 +
800 n7.v * (1 - t1) * (1 + t2) * (1 - t2) * 0.5;
801 ex = -(n0.v * ((1 - t2) * (2 * t1 + t2) * jac[0][0] +
802 (1 - t1) * (t1 + 2 * t2) * jac[1][0]) * 0.25 +
803 n1.v * ((1 - t2) * (2 * t1 - t2) * jac[0][0] -
804 (1 + t1) * (t1 - 2 * t2) * jac[1][0]) * 0.25 +
805 n2.v * ((1 + t2) * (2 * t1 + t2) * jac[0][0] +
806 (1 + t1) * (t1 + 2 * t2) * jac[1][0]) * 0.25 +
807 n3.v * ((1 + t2) * (2 * t1 - t2) * jac[0][0] -
808 (1 - t1) * (t1 - 2 * t2) * jac[1][0]) * 0.25 +
809 n4.v * (t1 * (t2 - 1) * jac[0][0] +
810 (t1 - 1) * (t1 + 1) * jac[1][0] * 0.5) +
811 n5.v * ((1 - t2) * (1 + t2) * jac[0][0] * 0.5 -
812 (1 + t1) * t2 * jac[1][0]) +
813 n6.v * (-t1 * (1 + t2) * jac[0][0] +
814 (1 - t1) * (1 + t1) * jac[1][0] * 0.5) +
815 n7.v * ((t2 - 1) * (t2 + 1) * jac[0][0] * 0.5 +
816 (t1 - 1) * t2 * jac[1][0])) * invdet;
817 ey = -(n0.v * ((1 - t2) * (2 * t1 + t2) * jac[0][1] +
818 (1 - t1) * (t1 + 2 * t2) * jac[1][1]) * 0.25 +
819 n1.v * ((1 - t2) * (2 * t1 - t2) * jac[0][1] -
820 (1 + t1) * (t1 - 2 * t2) * jac[1][1]) * 0.25 +
821 n2.v * ((1 + t2) * (2 * t1 + t2) * jac[0][1] +
822 (1 + t1) * (t1 + 2 * t2) * jac[1][1]) * 0.25 +
823 n3.v * ((1 + t2) * (2 * t1 - t2) * jac[0][1] -
824 (1 - t1) * (t1 - 2 * t2) * jac[1][1]) * 0.25 +
825 n4.v * (t1 * (t2 - 1) * jac[0][1] +
826 (t1 - 1) * (t1 + 1) * jac[1][1] * 0.5) +
827 n5.v * ((1 - t2) * (1 + t2) * jac[0][1] * 0.5 -
828 (1 + t1) * t2 * jac[1][1]) +
829 n6.v * (-t1 * (1 + t2) * jac[0][1] +
830 (1 - t1) * (1 + t1) * jac[1][1] * 0.5) +
831 n7.v * ((t2 - 1) * (t2 + 1) * jac[0][1] * 0.5 +
832 (t1 - 1) * t2 * jac[1][1])) * invdet;
833 }
834
835 // Transform field to global coordinates
836 UnmapFields(ex, ey, ez, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
837
838 // Drift medium?
839 if (m_debug) {
840 std::cout << m_className << "::ElectricField:\n";
841 std::cout << " Material " << element.matmap << ", drift flag "
842 << materials[element.matmap].driftmedium << ".\n";
843 }
844 m = materials[element.matmap].medium;
845 status = -5;
846 if (materials[element.matmap].driftmedium) {
847 if (m && m->IsDriftable()) status = 0;
848 }
849}
Garfield::ComponentBase::m_ready
bool m_ready
Ready for use?
Definition: ComponentBase.hh:239
Garfield::ComponentBase::m_debug
bool m_debug
Switch on/off debugging messages.
Definition: ComponentBase.hh:274
Garfield::ComponentFieldMap::MapCoordinates
void MapCoordinates(double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
Move (xpos, ypos, zpos) to field map coordinates.
Definition: ComponentFieldMap.cc:2441
Garfield::ComponentFieldMap::PrintWarning
void PrintWarning(const std::string &header)
Definition: ComponentFieldMap.hh:201
Garfield::ComponentFieldMap::PrintElement
void PrintElement(const std::string &header, const double x, const double y, const double z, const double t1, const double t2, const double t3, const double t4, const unsigned int i, const unsigned int n, const int iw=-1) const
Definition: ComponentFieldMap.cc:2740
Garfield::ComponentFieldMap::FindElement5
int FindElement5(const double x, const double y, const double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det)
Find the element for a point in curved quadratic quadrilaterals.
Definition: ComponentFieldMap.cc:166
Garfield::ComponentFieldMap::materials
std::vector< Material > materials
Definition: ComponentFieldMap.hh:133
Garfield::ComponentFieldMap::elements
std::vector< Element > elements
Definition: ComponentFieldMap.hh:108
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) const
Move (ex, ey, ez) to global coordinates.
Definition: ComponentFieldMap.cc:2559
Garfield::ComponentFieldMap::PrintNotReady
void PrintNotReady(const std::string &header) const
Definition: ComponentFieldMap.hh:207

◆ 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

Calculate the drift field at given point.

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

Status flags: 

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

Implements Garfield::ComponentFieldMap.

Definition at line 709 of file ComponentAnsys121.cc.

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

Referenced by ElectricField().

◆ GetAspectRatio()

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

Implements Garfield::ComponentFieldMap.

Definition at line 1110 of file ComponentAnsys121.cc.

1111 {
1112
1113 if (i >= elements.size()) {
1114 dmin = dmax = 0.;
1115 return;
1116 }
1117
1118 const Element& element = elements[i];
1119 const int np = 8;
1120 // Loop over all pairs of vertices.
1121 for (int j = 0; j < np - 1; ++j) {
1122 const Node& nj = nodes[element.emap[j]];
1123 for (int k = j + 1; k < np; ++k) {
1124 const Node& nk = nodes[element.emap[k]];
1125 // Compute distance.
1126 const double dx = nj.x - nk.x;
1127 const double dy = nj.y - nk.y;
1128 const double dist = sqrt(dx * dx + dy * dy);
1129 if (k == 1) {
1130 dmin = dmax = dist;
1131 } else {
1132 if (dist < dmin) dmin = dist;
1133 if (dist > dmax) dmax = dist;
1134 }
1135 }
1136 }
1137}
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314

◆ GetElementVolume()

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

Implements Garfield::ComponentFieldMap.

Definition at line 1096 of file ComponentAnsys121.cc.

1096 {
1097
1098 if (i >= elements.size()) return 0.;
1099 const Element& element = elements[i];
1100 const Node& n0 = nodes[element.emap[0]];
1101 const Node& n1 = nodes[element.emap[1]];
1102 const Node& n2 = nodes[element.emap[2]];
1103 const Node& n3 = nodes[element.emap[3]];
1104 const double surf = 0.5 *
1105 (fabs((n1.x - n0.x) * (n2.y - n0.y) - (n2.x - n0.x) * (n1.y - n0.y)) +
1106 fabs((n3.x - n0.x) * (n2.y - n0.y) - (n2.x - n0.x) * (n3.y - n0.y)));
1107 return surf;
1108}
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615

◆ GetMedium()

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

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

Implements Garfield::ComponentFieldMap.

Definition at line 1029 of file ComponentAnsys121.cc.

1030 {
1031
1032 // Copy the coordinates.
1033 double x = xin, y = yin, z = 0.;
1034
1035 // Map the coordinates onto field map coordinates.
1036 bool xmirr, ymirr, zmirr;
1037 double rcoordinate, rotation;
1038 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
1039
1040 if (zin < zMinBoundingBox || z > zMaxBoundingBox) {
1041 return NULL;
1042 }
1043
1044 // Do not proceed if not properly initialised.
1045 if (!m_ready) {
1046 PrintNotReady("GetMedium");
1047 return NULL;
1048 }
1049 if (m_warning) PrintWarning("GetMedium");
1050
1051 // Find the element that contains this point.
1052 double t1, t2, t3, t4, jac[4][4], det;
1053 const int imap = FindElement5(x, y, z, t1, t2, t3, t4, jac, det);
1054 if (imap < 0) {
1055 if (m_debug) {
1056 std::cerr << m_className << "::GetMedium:\n";
1057 std::cerr << " Point (" << x << ", " << y << ") not in the mesh.\n";
1058 }
1059 return NULL;
1060 }
1061 const Element& element = elements[imap];
1062 if (element.matmap >= m_nMaterials) {
1063 if (m_debug) {
1064 std::cerr << m_className << "::GetMedium:\n";
1065 std::cerr << " Point (" << x << ", " << y << ")"
1066 << " has out of range material number " << imap << ".\n";
1067 }
1068 return NULL;
1069 }
1070
1071 if (m_debug) {
1072 PrintElement("GetMedium", x, y, z, t1, t2, t3, t4, imap, 8);
1073 }
1074
1075 // Assign a medium.
1076 return materials[element.matmap].medium;
1077}

◆ 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 19 of file ComponentAnsys121.cc.

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

◆ IsInBoundingBox()

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

Reimplemented from Garfield::ComponentFieldMap.

Definition at line 39 of file ComponentAnsys121.hh.

40 {
41 return x >= xMinBoundingBox && x <= xMaxBoundingBox &&
42 y >= yMinBoundingBox && y <= yMaxBoundingBox &&
43 z >= zMinBoundingBox && y <= zMaxBoundingBox;
44 }

◆ SetRangeZ()

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

Definition at line 1079 of file ComponentAnsys121.cc.

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

◆ SetWeightingField()

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

Definition at line 595 of file ComponentAnsys121.cc.

596 {
597
598 if (!m_ready) {
599 PrintNotReady("SetWeightingField");
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) {
624 ++nWeightingFields;
625 wfields.resize(nWeightingFields);
626 wfieldsOk.resize(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

Verify periodicities.

Implements Garfield::ComponentFieldMap.

Definition at line 1090 of file ComponentAnsys121.cc.

1090 {
1091
1092 UpdatePeriodicity2d();
1093 UpdatePeriodicityCommon();
1094}

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

Calculate the weighting field at a given point and for a given electrode.

Parameters
x,y,zcoordinates [cm].
wx,wy,wzcomponents of the weighting field [1/cm].
labelname of the electrode

Implements Garfield::ComponentFieldMap.

Definition at line 851 of file ComponentAnsys121.cc.

853 {
854
855 // Initial values
856 wx = wy = wz = 0;
857
858 // Do not proceed if not properly initialised.
859 if (!m_ready) return;
860
861 // Look for the label.
862 int iw = 0;
863 bool found = false;
864 for (int i = nWeightingFields; i--;) {
865 if (wfields[i] == label) {
866 iw = i;
867 found = true;
868 break;
869 }
870 }
871
872 // Do not proceed if the requested weighting field does not exist.
873 if (!found) return;
874 // Check if the weighting field is properly initialised.
875 if (!wfieldsOk[iw]) return;
876
877 // Copy the coordinates.
878 double x = xin, y = yin, z = zin;
879
880 // Map the coordinates onto field map coordinates
881 bool xmirr, ymirr, zmirr;
882 double rcoordinate, rotation;
883 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
884
885 if (m_warning) PrintWarning("WeightingField");
886
887 // Find the element that contains this point.
888 double t1, t2, t3, t4, jac[4][4], det;
889 const int imap = FindElement5(x, y, z, t1, t2, t3, t4, jac, det);
890 // Check if the point is in the mesh.
891 if (imap < 0) return;
892
893 if (m_debug) {
894 PrintElement("WeightingField", x, y, z, t1, t2, t3, t4, imap, 8, iw);
895 }
896
897 const Element& element = elements[imap];
898 const Node& n0 = nodes[element.emap[0]];
899 const Node& n1 = nodes[element.emap[1]];
900 const Node& n2 = nodes[element.emap[2]];
901 const Node& n3 = nodes[element.emap[3]];
902 const Node& n4 = nodes[element.emap[4]];
903 const Node& n5 = nodes[element.emap[5]];
904 // Calculate quadrilateral field, which can degenerate to a triangular field
905 const double invdet = 1. / det;
906 if (elements[imap].degenerate) {
907 wx = -(n0.w[iw] * (4 * t1 - 1) * jac[0][1] +
908 n1.w[iw] * (4 * t2 - 1) * jac[1][1] +
909 n2.w[iw] * (4 * t3 - 1) * jac[2][1] +
910 n3.w[iw] * (4 * t2 * jac[0][1] + 4 * t1 * jac[1][1]) +
911 n4.w[iw] * (4 * t3 * jac[0][1] + 4 * t1 * jac[2][1]) +
912 n5.w[iw] * (4 * t3 * jac[1][1] + 4 * t2 * jac[2][1])) * invdet;
913 wy = -(n0.w[iw] * (4 * t1 - 1) * jac[0][2] +
914 n1.w[iw] * (4 * t2 - 1) * jac[1][2] +
915 n2.w[iw] * (4 * t3 - 1) * jac[2][2] +
916 n3.w[iw] * (4 * t2 * jac[0][2] + 4 * t1 * jac[1][2]) +
917 n4.w[iw] * (4 * t3 * jac[0][2] + 4 * t1 * jac[2][2]) +
918 n5.w[iw] * (4 * t3 * jac[1][2] + 4 * t2 * jac[2][2])) * invdet;
919 } else {
920 const Node& n6 = nodes[element.emap[6]];
921 const Node& n7 = nodes[element.emap[7]];
922 wx = -(n0.w[iw] * ((1 - t2) * (2 * t1 + t2) * jac[0][0] +
923 (1 - t1) * (t1 + 2 * t2) * jac[1][0]) * 0.25 +
924 n1.w[iw] * ((1 - t2) * (2 * t1 - t2) * jac[0][0] -
925 (1 + t1) * (t1 - 2 * t2) * jac[1][0]) * 0.25 +
926 n2.w[iw] * ((1 + t2) * (2 * t1 + t2) * jac[0][0] +
927 (1 + t1) * (t1 + 2 * t2) * jac[1][0]) * 0.25 +
928 n3.w[iw] * ((1 + t2) * (2 * t1 - t2) * jac[0][0] -
929 (1 - t1) * (t1 - 2 * t2) * jac[1][0]) * 0.25 +
930 n4.w[iw] * (t1 * (t2 - 1) * jac[0][0] +
931 (t1 - 1) * (t1 + 1) * jac[1][0] * 0.5) +
932 n5.w[iw] * ((1 - t2) * (1 + t2) * jac[0][0] * 0.5 -
933 (1 + t1) * t2 * jac[1][0]) +
934 n6.w[iw] * (-t1 * (1 + t2) * jac[0][0] +
935 (1 - t1) * (1 + t1) * jac[1][0] * 0.5) +
936 n7.w[iw] * ((t2 - 1) * (1 + t2) * jac[0][0] * 0.5 +
937 (t1 - 1) * t2 * jac[1][0])) * invdet;
938 wy = -(n0.w[iw] * ((1 - t2) * (2 * t1 + t2) * jac[0][1] +
939 (1 - t1) * (t1 + 2 * t2) * jac[1][1]) * 0.25 +
940 n1.w[iw] * ((1 - t2) * (2 * t1 - t2) * jac[0][1] -
941 (1 + t1) * (t1 - 2 * t2) * jac[1][1]) * 0.25 +
942 n2.w[iw] * ((1 + t2) * (2 * t1 + t2) * jac[0][1] +
943 (1 + t1) * (t1 + 2 * t2) * jac[1][1]) * 0.25 +
944 n3.w[iw] * ((1 + t2) * (2 * t1 - t2) * jac[0][1] -
945 (1 - t1) * (t1 - 2 * t2) * jac[1][1]) * 0.25 +
946 n4.w[iw] * (t1 * (t2 - 1) * jac[0][1] +
947 (t1 - 1) * (t1 + 1) * jac[1][1] * 0.5) +
948 n5.w[iw] * ((1 - t2) * (1 + t2) * jac[0][1] * 0.5 -
949 (1 + t1) * t2 * jac[1][1]) +
950 n6.w[iw] * (-t1 * (1 + t2) * jac[0][1] +
951 (1 - t1) * (1 + t1) * jac[1][1] * 0.5) +
952 n7.w[iw] * ((t2 - 1) * (t2 + 1) * jac[0][1] * 0.5 +
953 (t1 - 1) * t2 * jac[1][1])) * invdet;
954 }
955
956 // Transform field to global coordinates
957 UnmapFields(wx, wy, wz, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
958}

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

962 {
963
964 // Do not proceed if not properly initialised.
965 if (!m_ready) return 0.;
966
967 // Look for the label.
968 int iw = 0;
969 bool found = false;
970 for (int i = nWeightingFields; i--;) {
971 if (wfields[i] == label) {
972 iw = i;
973 found = true;
974 break;
975 }
976 }
977
978 // Do not proceed if the requested weighting field does not exist.
979 if (!found) return 0.;
980 // Check if the weighting field is properly initialised.
981 if (!wfieldsOk[iw]) return 0.;
982
983 // Copy the coordinates.
984 double x = xin, y = yin, z = zin;
985
986 // Map the coordinates onto field map coordinates.
987 bool xmirr, ymirr, zmirr;
988 double rcoordinate, rotation;
989 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
990
991 if (m_warning) PrintWarning("WeightingPotential");
992
993 // Find the element that contains this point.
994 double t1, t2, t3, t4, jac[4][4], det;
995 const int imap = FindElement5(x, y, z, t1, t2, t3, t4, jac, det);
996 // Check if the point is in the mesh
997 if (imap < 0) return 0.;
998
999 if (m_debug) {
1000 PrintElement("WeightingPotential", x, y, z, t1, t2, t3, t4, imap, 8, iw);
1001 }
1002
1003 // Calculate quadrilateral field, which can degenerate to a triangular field
1004 const Element& element = elements[imap];
1005 const Node& n0 = nodes[element.emap[0]];
1006 const Node& n1 = nodes[element.emap[1]];
1007 const Node& n2 = nodes[element.emap[2]];
1008 const Node& n3 = nodes[element.emap[3]];
1009 const Node& n4 = nodes[element.emap[4]];
1010 const Node& n5 = nodes[element.emap[5]];
1011 if (element.degenerate) {
1012 return n0.w[iw] * t1 * (2 * t1 - 1) + n1.w[iw] * t2 * (2 * t2 - 1) +
1013 n2.w[iw] * t3 * (2 * t3 - 1) + 4 * n3.w[iw] * t1 * t2 +
1014 4 * n4.w[iw] * t1 * t3 + 4 * n5.w[iw] * t2 * t3;
1015 }
1016
1017 const Node& n6 = nodes[element.emap[6]];
1018 const Node& n7 = nodes[element.emap[7]];
1019 return -n0.w[iw] * (1 - t1) * (1 - t2) * (1 + t1 + t2) * 0.25 -
1020 n1.w[iw] * (1 + t1) * (1 - t2) * (1 - t1 + t2) * 0.25 -
1021 n2.w[iw] * (1 + t1) * (1 + t2) * (1 - t1 - t2) * 0.25 -
1022 n3.w[iw] * (1 - t1) * (1 + t2) * (1 + t1 - t2) * 0.25 +
1023 n4.w[iw] * (1 - t1) * (1 + t1) * (1 - t2) * 0.5 +
1024 n5.w[iw] * (1 + t1) * (1 + t2) * (1 - t2) * 0.5 +
1025 n6.w[iw] * (1 - t1) * (1 + t1) * (1 + t2) * 0.5 +
1026 n7.w[iw] * (1 - t1) * (1 + t2) * (1 - t2) * 0.5;
1027}
The documentation for this class was generated from the following files: