Garfield++ 4.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Public Member Functions | Protected Member Functions | List of all members
Garfield::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) override
 Get the medium at a given location (x, y, z).
 
void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status) override
 
void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status) override
 Calculate the drift field [V/cm] and potential [V] at (x, y, z).
 
void WeightingField (const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label) override
 
double WeightingPotential (const double x, const double y, const double z, const std::string &label) override
 
bool Initialise (std::string elist="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)
 Import a weighting field map.
 
void SetRangeZ (const double zmin, const double zmax)
 Set the limits of the active region along z.
 
- Public Member Functions inherited from Garfield::ComponentFieldMap
 ComponentFieldMap ()=delete
 Default constructor.
 
 ComponentFieldMap (const std::string &name)
 Constructor.
 
virtual ~ComponentFieldMap ()
 Destructor.
 
virtual void SetRange ()
 Calculate x, y, z, V and angular ranges.
 
void PrintRange ()
 Show x, y, z, V and angular ranges.
 
bool IsInBoundingBox (const double x, const double y, const double z) const
 
bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax) override
 Get the bounding box coordinates.
 
bool GetElementaryCell (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax) override
 Get the coordinates of the elementary cell.
 
bool GetVoltageRange (double &vmin, double &vmax) override
 Calculate the voltage range [V].
 
void PrintMaterials ()
 List all currently defined materials.
 
void DriftMedium (const unsigned int imat)
 Flag a field map material as a drift medium.
 
void NotDriftMedium (const unsigned int imat)
 Flag a field map materials as a non-drift medium.
 
size_t GetNumberOfMaterials () const
 Return the number of materials in the field map.
 
double GetPermittivity (const unsigned int imat) const
 Return the permittivity of a field map material.
 
double GetConductivity (const unsigned int imat) const
 Return the conductivity of a field map material.
 
void SetMedium (const unsigned int imat, Medium *medium)
 Associate a field map material with a Medium class.
 
Medium * GetMedium (const unsigned int i) const
 Return the Medium associated to a field map material.
 
virtual size_t GetNumberOfElements () const
 Return the number of mesh elements.
 
bool GetElement (const size_t i, double &vol, double &dmin, double &dmax)
 Return the volume and aspect ratio of a mesh element.
 
virtual bool GetElement (const size_t i, size_t &mat, bool &drift, std::vector< size_t > &nodes) const
 Return the material and node indices of a mesh element.
 
virtual size_t GetNumberOfNodes () const
 
virtual bool GetNode (const size_t i, double &x, double &y, double &z) const
 
void EnableCheckMapIndices (const bool on=true)
 
void EnableDeleteBackgroundElements (const bool on=true)
 Option to eliminate mesh elements in conductors (default: on).
 
void EnableTetrahedralTreeForElementSearch (const bool on=true)
 
void EnableConvergenceWarnings (const bool on=true)
 
virtual Medium * GetMedium (const double x, const double y, const double z)
 Get the medium at a given location (x, y, z).
 
- Public Member Functions inherited from Garfield::Component
 Component ()=delete
 Default constructor.
 
 Component (const std::string &name)
 Constructor.
 
virtual ~Component ()
 Destructor.
 
virtual void SetGeometry (Geometry *geo)
 Define the geometry.
 
virtual void Clear ()
 Reset.
 
virtual Medium * GetMedium (const double x, const double y, const double z)
 Get the medium at a given location (x, y, z).
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status)=0
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status)=0
 Calculate the drift field [V/cm] and potential [V] at (x, y, z).
 
virtual bool GetVoltageRange (double &vmin, double &vmax)=0
 Calculate the voltage range [V].
 
virtual void WeightingField (const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label)
 
virtual double WeightingPotential (const double x, const double y, const double z, const std::string &label)
 
virtual void DelayedWeightingField (const double x, const double y, const double z, const double t, double &wx, double &wy, double &wz, const std::string &label)
 
virtual double DelayedWeightingPotential (const double x, const double y, const double z, const double t, const std::string &label)
 
virtual void MagneticField (const double x, const double y, const double z, double &bx, double &by, double &bz, int &status)
 
void SetMagneticField (const double bx, const double by, const double bz)
 Set a constant magnetic field.
 
virtual bool IsReady ()
 Ready for use?
 
virtual bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
 Get the bounding box coordinates.
 
virtual bool GetElementaryCell (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
 Get the coordinates of the elementary cell.
 
double IntegrateFluxCircle (const double xc, const double yc, const double r, const unsigned int nI=50)
 
double IntegrateFluxSphere (const double xc, const double yc, const double zc, const double r, const unsigned int nI=20)
 
double IntegrateFluxParallelogram (const double x0, const double y0, const double z0, const double dx1, const double dy1, const double dz1, const double dx2, const double dy2, const double dz2, const unsigned int nU=20, const unsigned int nV=20)
 
double IntegrateWeightingFluxParallelogram (const std::string &label, const double x0, const double y0, const double z0, const double dx1, const double dy1, const double dz1, const double dx2, const double dy2, const double dz2, const unsigned int nU=20, const unsigned int nV=20)
 
double IntegrateFluxLine (const double x0, const double y0, const double z0, const double x1, const double y1, const double z1, const double xp, const double yp, const double zp, const unsigned int nI, const int isign=0)
 
virtual bool IsWireCrossed (const double x0, const double y0, const double z0, const double x1, const double y1, const double z1, double &xc, double &yc, double &zc, const bool centre, double &rc)
 
virtual bool IsInTrapRadius (const double q0, const double x0, const double y0, const double z0, double &xw, double &yw, double &rw)
 
void EnablePeriodicityX (const bool on=true)
 Enable simple periodicity in the $x$ direction.
 
void EnablePeriodicityY (const bool on=true)
 Enable simple periodicity in the $y$ direction.
 
void EnablePeriodicityZ (const bool on=true)
 Enable simple periodicity in the $z$ direction.
 
void IsPeriodic (bool &perx, bool &pery, bool &perz)
 Return periodicity flags.
 
void EnableMirrorPeriodicityX (const bool on=true)
 Enable mirror periodicity in the $x$ direction.
 
void EnableMirrorPeriodicityY (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void EnableMirrorPeriodicityZ (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void IsMirrorPeriodic (bool &perx, bool &pery, bool &perz)
 Return mirror periodicity flags.
 
void EnableAxialPeriodicityX (const bool on=true)
 Enable axial periodicity in the $x$ direction.
 
void EnableAxialPeriodicityY (const bool on=true)
 Enable axial periodicity in the $y$ direction.
 
void EnableAxialPeriodicityZ (const bool on=true)
 Enable axial periodicity in the $z$ direction.
 
void IsAxiallyPeriodic (bool &perx, bool &pery, bool &perz)
 Return axial periodicity flags.
 
void EnableRotationSymmetryX (const bool on=true)
 Enable rotation symmetry around the $x$ axis.
 
void EnableRotationSymmetryY (const bool on=true)
 Enable rotation symmetry around the $y$ axis.
 
void EnableRotationSymmetryZ (const bool on=true)
 Enable rotation symmetry around the $z$ axis.
 
void IsRotationSymmetric (bool &rotx, bool &roty, bool &rotz)
 Return rotation symmetry flags.
 
void EnableDebugging ()
 Switch on debugging messages.
 
void DisableDebugging ()
 Switch off debugging messages.
 
virtual bool HasAttachmentMap () const
 Does the component have attachment maps?
 
virtual bool HasVelocityMap () const
 Does the component have velocity maps?
 
virtual bool ElectronAttachment (const double, const double, const double, double &eta)
 Get the electron attachment coefficient.
 
virtual bool HoleAttachment (const double, const double, const double, double &eta)
 Get the hole attachment coefficient.
 
virtual bool ElectronVelocity (const double, const double, const double, double &vx, double &vy, double &vz)
 Get the electron drift velocity.
 
virtual bool HoleVelocity (const double, const double, const double, double &vx, double &vy, double &vz)
 Get the hole drift velocity.
 
virtual bool GetElectronLifetime (const double, const double, const double, double &etau)
 
virtual bool GetHoleLifetime (const double, const double, const double, double &htau)
 

Protected Member Functions

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

Additional Inherited Members

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

Detailed Description

Component for importing and interpolating 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("Ansys121") {
11 m_is3d = false;
12 // Default bounding box
13 m_minBoundingBox[2] = -50;
14 m_maxBoundingBox[2] = 50;
15}
Garfield::ComponentFieldMap::m_maxBoundingBox
std::array< double, 3 > m_maxBoundingBox
Definition: ComponentFieldMap.hh:149
Garfield::ComponentFieldMap::m_minBoundingBox
std::array< double, 3 > m_minBoundingBox
Definition: ComponentFieldMap.hh:148
Garfield::ComponentFieldMap::ComponentFieldMap
ComponentFieldMap()=delete
Default constructor.

◆ ~ComponentAnsys121()

Garfield::ComponentAnsys121::~ComponentAnsys121 ( )
inline

Destructor.

Definition at line 15 of file ComponentAnsys121.hh.

15{}

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 
)
overridevirtual

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

Implements Garfield::Component.

Definition at line 659 of file ComponentAnsys121.cc.

662 {
663 // Copy the coordinates
664 double x = xin, y = yin, z = 0.;
665
666 // Map the coordinates onto field map coordinates
667 bool xmirr, ymirr, zmirr;
668 double rcoordinate, rotation;
669 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
670
671 // Initial values
672 ex = ey = ez = volt = 0;
673 m = nullptr;
674 status = 0;
675
676 // Do not proceed if not properly initialised.
677 if (!m_ready) {
678 status = -10;
679 PrintNotReady("ElectricField");
680 return;
681 }
682
683 if (m_warning) PrintWarning("ElectricField");
684
685 if (zin < m_minBoundingBox[2] || zin > m_maxBoundingBox[2]) {
686 status = -5;
687 return;
688 }
689
690 // Find the element that contains this point
691 double t1, t2, t3, t4, jac[4][4], det;
692 const int imap = FindElement5(x, y, z, t1, t2, t3, t4, jac, det);
693 if (imap < 0) {
694 if (m_debug) {
695 std::cout << m_className << "::ElectricField:\n";
696 std::cout << " Point (" << x << ", " << y << ") not in the mesh.\n";
697 }
698 status = -6;
699 return;
700 }
701
702 const Element& element = m_elements[imap];
703 if (m_debug) {
704 PrintElement("ElectricField", x, y, z, t1, t2, t3, t4, element, 8);
705 }
706
707 const Node& n0 = m_nodes[element.emap[0]];
708 const Node& n1 = m_nodes[element.emap[1]];
709 const Node& n2 = m_nodes[element.emap[2]];
710 const Node& n3 = m_nodes[element.emap[3]];
711 const Node& n4 = m_nodes[element.emap[4]];
712 const Node& n5 = m_nodes[element.emap[5]];
713 // Calculate quadrilateral field, which can degenerate to a triangular field
714 const double invdet = 1. / det;
715 if (element.degenerate) {
716 volt = n0.v * t1 * (2 * t1 - 1) + n1.v * t2 * (2 * t2 - 1) +
717 n2.v * t3 * (2 * t3 - 1) + 4 * n3.v * t1 * t2 + 4 * n4.v * t1 * t3 +
718 4 * n5.v * t2 * t3;
719 ex = -(n0.v * (4 * t1 - 1) * jac[0][1] + n1.v * (4 * t2 - 1) * jac[1][1] +
720 n2.v * (4 * t3 - 1) * jac[2][1] +
721 n3.v * (4 * t2 * jac[0][1] + 4 * t1 * jac[1][1]) +
722 n4.v * (4 * t3 * jac[0][1] + 4 * t1 * jac[2][1]) +
723 n5.v * (4 * t3 * jac[1][1] + 4 * t2 * jac[2][1])) *
724 invdet;
725 ey = -(n0.v * (4 * t1 - 1) * jac[0][2] + n1.v * (4 * t2 - 1) * jac[1][2] +
726 n2.v * (4 * t3 - 1) * jac[2][2] +
727 n3.v * (4 * t2 * jac[0][2] + 4 * t1 * jac[1][2]) +
728 n4.v * (4 * t3 * jac[0][2] + 4 * t1 * jac[2][2]) +
729 n5.v * (4 * t3 * jac[1][2] + 4 * t2 * jac[2][2])) *
730 invdet;
731 } else {
732 const Node& n6 = m_nodes[element.emap[6]];
733 const Node& n7 = m_nodes[element.emap[7]];
734 volt = -n0.v * (1 - t1) * (1 - t2) * (1 + t1 + t2) * 0.25 -
735 n1.v * (1 + t1) * (1 - t2) * (1 - t1 + t2) * 0.25 -
736 n2.v * (1 + t1) * (1 + t2) * (1 - t1 - t2) * 0.25 -
737 n3.v * (1 - t1) * (1 + t2) * (1 + t1 - t2) * 0.25 +
738 n4.v * (1 - t1) * (1 + t1) * (1 - t2) * 0.5 +
739 n5.v * (1 + t1) * (1 + t2) * (1 - t2) * 0.5 +
740 n6.v * (1 - t1) * (1 + t1) * (1 + t2) * 0.5 +
741 n7.v * (1 - t1) * (1 + t2) * (1 - t2) * 0.5;
742 ex = -(n0.v * ((1 - t2) * (2 * t1 + t2) * jac[0][0] +
743 (1 - t1) * (t1 + 2 * t2) * jac[1][0]) *
744 0.25 +
745 n1.v * ((1 - t2) * (2 * t1 - t2) * jac[0][0] -
746 (1 + t1) * (t1 - 2 * t2) * jac[1][0]) *
747 0.25 +
748 n2.v * ((1 + t2) * (2 * t1 + t2) * jac[0][0] +
749 (1 + t1) * (t1 + 2 * t2) * jac[1][0]) *
750 0.25 +
751 n3.v * ((1 + t2) * (2 * t1 - t2) * jac[0][0] -
752 (1 - t1) * (t1 - 2 * t2) * jac[1][0]) *
753 0.25 +
754 n4.v * (t1 * (t2 - 1) * jac[0][0] +
755 (t1 - 1) * (t1 + 1) * jac[1][0] * 0.5) +
756 n5.v * ((1 - t2) * (1 + t2) * jac[0][0] * 0.5 -
757 (1 + t1) * t2 * jac[1][0]) +
758 n6.v * (-t1 * (1 + t2) * jac[0][0] +
759 (1 - t1) * (1 + t1) * jac[1][0] * 0.5) +
760 n7.v * ((t2 - 1) * (t2 + 1) * jac[0][0] * 0.5 +
761 (t1 - 1) * t2 * jac[1][0])) *
762 invdet;
763 ey = -(n0.v * ((1 - t2) * (2 * t1 + t2) * jac[0][1] +
764 (1 - t1) * (t1 + 2 * t2) * jac[1][1]) *
765 0.25 +
766 n1.v * ((1 - t2) * (2 * t1 - t2) * jac[0][1] -
767 (1 + t1) * (t1 - 2 * t2) * jac[1][1]) *
768 0.25 +
769 n2.v * ((1 + t2) * (2 * t1 + t2) * jac[0][1] +
770 (1 + t1) * (t1 + 2 * t2) * jac[1][1]) *
771 0.25 +
772 n3.v * ((1 + t2) * (2 * t1 - t2) * jac[0][1] -
773 (1 - t1) * (t1 - 2 * t2) * jac[1][1]) *
774 0.25 +
775 n4.v * (t1 * (t2 - 1) * jac[0][1] +
776 (t1 - 1) * (t1 + 1) * jac[1][1] * 0.5) +
777 n5.v * ((1 - t2) * (1 + t2) * jac[0][1] * 0.5 -
778 (1 + t1) * t2 * jac[1][1]) +
779 n6.v * (-t1 * (1 + t2) * jac[0][1] +
780 (1 - t1) * (1 + t1) * jac[1][1] * 0.5) +
781 n7.v * ((t2 - 1) * (t2 + 1) * jac[0][1] * 0.5 +
782 (t1 - 1) * t2 * jac[1][1])) *
783 invdet;
784 }
785
786 // Transform field to global coordinates
787 UnmapFields(ex, ey, ez, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
788
789 // Drift medium?
790 if (m_debug) {
791 std::cout << m_className << "::ElectricField:\n";
792 std::cout << " Material " << element.matmap << ", drift flag "
793 << m_materials[element.matmap].driftmedium << ".\n";
794 }
795 m = m_materials[element.matmap].medium;
796 status = -5;
797 if (m_materials[element.matmap].driftmedium) {
798 if (m && m->IsDriftable()) status = 0;
799 }
800}
Garfield::ComponentFieldMap::MapCoordinates
void MapCoordinates(double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
Move (xpos, ypos, zpos) to field map coordinates.
Definition: ComponentFieldMap.cc:2021
Garfield::ComponentFieldMap::PrintWarning
void PrintWarning(const std::string &header)
Definition: ComponentFieldMap.cc:2375
Garfield::ComponentFieldMap::m_materials
std::vector< Material > m_materials
Definition: ComponentFieldMap.hh:138
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:190
Garfield::ComponentFieldMap::m_elements
std::vector< Element > m_elements
Definition: ComponentFieldMap.hh:113
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 Element &element, const unsigned int n, const int iw=-1) const
Definition: ComponentFieldMap.cc:2395
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:2138
Garfield::ComponentFieldMap::PrintNotReady
void PrintNotReady(const std::string &header) const
Definition: ComponentFieldMap.cc:2383
Garfield::Component::m_debug
bool m_debug
Switch on/off debugging messages.
Definition: Component.hh:341
Garfield::Component::m_className
std::string m_className
Class name.
Definition: Component.hh:329
Garfield::Component::m_ready
bool m_ready
Ready for use?
Definition: Component.hh:338
Element
Definition: neBEM.h:155

◆ 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 
)
overridevirtual

Calculate the drift field at given point.

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

Status flags: 

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

Implements Garfield::Component.

Definition at line 652 of file ComponentAnsys121.cc.

654 {
655 double v;
656 ElectricField(x, y, z, ex, ey, ez, v, m, status);
657}
Garfield::ComponentAnsys121::ElectricField
void ElectricField(const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status) override
Definition: ComponentAnsys121.cc:652

Referenced by ElectricField().

◆ GetAspectRatio()

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

Implements Garfield::ComponentFieldMap.

Definition at line 1047 of file ComponentAnsys121.cc.

1048 {
1049 if (i >= m_elements.size()) {
1050 dmin = dmax = 0.;
1051 return;
1052 }
1053
1054 const Element& element = m_elements[i];
1055 const int np = 8;
1056 // Loop over all pairs of vertices.
1057 for (int j = 0; j < np - 1; ++j) {
1058 const Node& nj = m_nodes[element.emap[j]];
1059 for (int k = j + 1; k < np; ++k) {
1060 const Node& nk = m_nodes[element.emap[k]];
1061 // Compute distance.
1062 const double dx = nj.x - nk.x;
1063 const double dy = nj.y - nk.y;
1064 const double dist = sqrt(dx * dx + dy * dy);
1065 if (k == 1) {
1066 dmin = dmax = dist;
1067 } else {
1068 if (dist < dmin) dmin = dist;
1069 if (dist > dmax) dmax = dist;
1070 }
1071 }
1072 }
1073}
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314

◆ GetElementVolume()

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

Implements Garfield::ComponentFieldMap.

Definition at line 1033 of file ComponentAnsys121.cc.

1033 {
1034 if (i >= m_elements.size()) return 0.;
1035 const Element& element = m_elements[i];
1036 const Node& n0 = m_nodes[element.emap[0]];
1037 const Node& n1 = m_nodes[element.emap[1]];
1038 const Node& n2 = m_nodes[element.emap[2]];
1039 const Node& n3 = m_nodes[element.emap[3]];
1040 const double surf =
1041 0.5 *
1042 (fabs((n1.x - n0.x) * (n2.y - n0.y) - (n2.x - n0.x) * (n1.y - n0.y)) +
1043 fabs((n3.x - n0.x) * (n2.y - n0.y) - (n2.x - n0.x) * (n3.y - n0.y)));
1044 return surf;
1045}
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 
)
overridevirtual

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

Reimplemented from Garfield::Component.

Definition at line 970 of file ComponentAnsys121.cc.

971 {
972 // Copy the coordinates.
973 double x = xin, y = yin, z = 0.;
974
975 // Map the coordinates onto field map coordinates.
976 bool xmirr, ymirr, zmirr;
977 double rcoordinate, rotation;
978 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
979
980 if (zin < m_minBoundingBox[2] || z > m_maxBoundingBox[2]) {
981 return nullptr;
982 }
983
984 // Do not proceed if not properly initialised.
985 if (!m_ready) {
986 PrintNotReady("GetMedium");
987 return nullptr;
988 }
989 if (m_warning) PrintWarning("GetMedium");
990
991 // Find the element that contains this point.
992 double t1, t2, t3, t4, jac[4][4], det;
993 const int imap = FindElement5(x, y, z, t1, t2, t3, t4, jac, det);
994 if (imap < 0) {
995 if (m_debug) {
996 std::cerr << m_className << "::GetMedium:\n";
997 std::cerr << " Point (" << x << ", " << y << ") not in the mesh.\n";
998 }
999 return nullptr;
1000 }
1001 const Element& element = m_elements[imap];
1002 if (element.matmap >= m_materials.size()) {
1003 if (m_debug) {
1004 std::cerr << m_className << "::GetMedium:\n";
1005 std::cerr << " Point (" << x << ", " << y << ")"
1006 << " has out of range material number " << imap << ".\n";
1007 }
1008 return nullptr;
1009 }
1010
1011 if (m_debug) {
1012 PrintElement("GetMedium", x, y, z, t1, t2, t3, t4, element, 8);
1013 }
1014
1015 // Assign a medium.
1016 return m_materials[element.matmap].medium;
1017}

◆ 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" 
)

Import a field map.

Parameters
elistname of the file containing the list of elements
nlistname of the file containing the list of nodes
mplistname of the file containing the list of materials
prnsolname of the file containing the nodal solutions
unitlength unit

Definition at line 17 of file ComponentAnsys121.cc.

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

◆ SetRangeZ()

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

Set the limits of the active region along z.

Definition at line 1019 of file ComponentAnsys121.cc.

1019 {
1020 if (fabs(zmax - zmin) <= 0.) {
1021 std::cerr << m_className << "::SetRangeZ: Zero range is not permitted.\n";
1022 return;
1023 }
1024 m_minBoundingBox[2] = m_mapmin[2] = std::min(zmin, zmax);
1025 m_maxBoundingBox[2] = m_mapmax[2] = std::max(zmin, zmax);
1026}
Garfield::ComponentFieldMap::m_mapmin
std::array< double, 3 > m_mapmin
Definition: ComponentFieldMap.hh:152
Garfield::ComponentFieldMap::m_mapmax
std::array< double, 3 > m_mapmax
Definition: ComponentFieldMap.hh:153

◆ SetWeightingField()

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

Import a weighting field map.

Definition at line 549 of file ComponentAnsys121.cc.

550 {
551 if (!m_ready) {
552 PrintNotReady("SetWeightingField");
553 std::cerr << " Weighting field cannot be added.\n";
554 return false;
555 }
556
557 // Open the voltage list.
558 std::ifstream fprnsol;
559 fprnsol.open(prnsol.c_str(), std::ios::in);
560 if (fprnsol.fail()) {
561 PrintCouldNotOpen("SetWeightingField", prnsol);
562 return false;
563 }
564
565 // Check if a weighting field with the same label already exists.
566 const size_t iw = GetOrCreateWeightingFieldIndex(label);
567 if (iw + 1 != m_wfields.size()) {
568 std::cout << m_className << "::SetWeightingField:\n"
569 << " Replacing existing weighting field " << label << ".\n";
570 }
571
572 // Buffer for reading
573 constexpr int size = 100;
574 char line[size];
575
576 bool ok = true;
577 // Read the voltage list.
578 int il = 0;
579 unsigned int nread = 0;
580 bool readerror = false;
581 while (fprnsol.getline(line, size, '\n')) {
582 il++;
583 // Skip page feed in batch page title
584 if (strstr(line,"VERSION") != NULL) {
585 fprnsol.getline(line, size, '\n');
586 il++;
587 fprnsol.getline(line, size, '\n');
588 il++;
589 continue;
590 }
591 // Split the line in tokens.
592 char* token = NULL;
593 token = strtok(line, " ");
594 // Skip blank lines and headers.
595 if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
596 int(token[0]) == 10 || int(token[0]) == 13 ||
597 strcmp(token, "PRINT") == 0 || strcmp(token, "*****") == 0 ||
598 strcmp(token, "LOAD") == 0 || strcmp(token, "TIME=") == 0 ||
599 strcmp(token, "MAXIMUM") == 0 || strcmp(token, "VALUE") == 0 ||
600 strcmp(token, "NODE") == 0)
601 continue;
602 // Read the element.
603 int inode = ReadInteger(token, -1, readerror);
604 token = strtok(NULL, " ");
605 double volt = ReadDouble(token, -1, readerror);
606 // Check the syntax.
607 if (readerror) {
608 std::cerr << m_className << "::SetWeightingField:\n";
609 std::cerr << " Error reading file " << prnsol << " (line " << il
610 << ").\n";
611 fprnsol.close();
612 return false;
613 }
614 // Check node number and store if OK.
615 if (inode < 1 || inode > (int)m_nodes.size()) {
616 std::cerr << m_className << "::SetWeightingField:\n";
617 std::cerr << " Node number " << inode << " out of range\n";
618 std::cerr << " on potential file " << prnsol << " (line " << il
619 << ").\n";
620 ok = false;
621 } else {
622 m_nodes[inode - 1].w[iw] = volt;
623 nread++;
624 }
625 }
626 // Close the file.
627 fprnsol.close();
628 // Tell how many lines read.
629 std::cout << m_className << "::SetWeightingField:\n";
630 std::cout << " Read " << nread << " potentials from file " << prnsol
631 << ".\n";
632 // Check the number of nodes.
633 if (nread != m_nodes.size()) {
634 std::cerr << m_className << "::SetWeightingField:\n"
635 << " Number of nodes read from potential file " << prnsol
636 << " (" << nread << ")\n does not match the node list ("
637 << m_nodes.size() << ").\n";
638 ok = false;
639 }
640
641 // Set the ready flag.
642 m_wfieldsOk[iw] = ok;
643 if (!ok) {
644 std::cerr << m_className << "::SetWeightingField:\n";
645 std::cerr << " Field map could not be read "
646 << "and cannot be interpolated.\n";
647 return false;
648 }
649 return true;
650}
Garfield::ComponentFieldMap::m_wfieldsOk
std::vector< bool > m_wfieldsOk
Definition: ComponentFieldMap.hh:141
Garfield::ComponentFieldMap::m_wfields
std::vector< std::string > m_wfields
Definition: ComponentFieldMap.hh:140
Garfield::ComponentFieldMap::GetOrCreateWeightingFieldIndex
size_t GetOrCreateWeightingFieldIndex(const std::string &label)
Definition: ComponentFieldMap.cc:2357

◆ UpdatePeriodicity()

void Garfield::ComponentAnsys121::UpdatePeriodicity ( )
overrideprotectedvirtual

Verify periodicities.

Implements Garfield::Component.

Definition at line 1028 of file ComponentAnsys121.cc.

1028 {
1029 UpdatePeriodicity2d();
1030 UpdatePeriodicityCommon();
1031}

◆ WeightingField()

void Garfield::ComponentAnsys121::WeightingField ( const double  x,
const double  y,
const double  z,
double &  wx,
double &  wy,
double &  wz,
const std::string &  label 
)
overridevirtual

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

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

Reimplemented from Garfield::Component.

Definition at line 802 of file ComponentAnsys121.cc.

804 {
805 // Initial values
806 wx = wy = wz = 0;
807
808 // Do not proceed if not properly initialised.
809 if (!m_ready) return;
810
811 // Look for the label.
812 const size_t iw = GetWeightingFieldIndex(label);
813 // Do not proceed if the requested weighting field does not exist.
814 if (iw == m_wfields.size()) return;
815 // Check if the weighting field is properly initialised.
816 if (!m_wfieldsOk[iw]) return;
817
818 // Copy the coordinates.
819 double x = xin, y = yin, z = zin;
820
821 // Map the coordinates onto field map coordinates
822 bool xmirr, ymirr, zmirr;
823 double rcoordinate, rotation;
824 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
825
826 if (m_warning) PrintWarning("WeightingField");
827
828 // Find the element that contains this point.
829 double t1, t2, t3, t4, jac[4][4], det;
830 const int imap = FindElement5(x, y, z, t1, t2, t3, t4, jac, det);
831 // Check if the point is in the mesh.
832 if (imap < 0) return;
833
834 const Element& element = m_elements[imap];
835 if (m_debug) {
836 PrintElement("WeightingField", x, y, z, t1, t2, t3, t4, element, 8, iw);
837 }
838 const Node& n0 = m_nodes[element.emap[0]];
839 const Node& n1 = m_nodes[element.emap[1]];
840 const Node& n2 = m_nodes[element.emap[2]];
841 const Node& n3 = m_nodes[element.emap[3]];
842 const Node& n4 = m_nodes[element.emap[4]];
843 const Node& n5 = m_nodes[element.emap[5]];
844 // Calculate quadrilateral field, which can degenerate to a triangular field
845 const double invdet = 1. / det;
846 if (m_elements[imap].degenerate) {
847 wx = -(n0.w[iw] * (4 * t1 - 1) * jac[0][1] +
848 n1.w[iw] * (4 * t2 - 1) * jac[1][1] +
849 n2.w[iw] * (4 * t3 - 1) * jac[2][1] +
850 n3.w[iw] * (4 * t2 * jac[0][1] + 4 * t1 * jac[1][1]) +
851 n4.w[iw] * (4 * t3 * jac[0][1] + 4 * t1 * jac[2][1]) +
852 n5.w[iw] * (4 * t3 * jac[1][1] + 4 * t2 * jac[2][1])) *
853 invdet;
854 wy = -(n0.w[iw] * (4 * t1 - 1) * jac[0][2] +
855 n1.w[iw] * (4 * t2 - 1) * jac[1][2] +
856 n2.w[iw] * (4 * t3 - 1) * jac[2][2] +
857 n3.w[iw] * (4 * t2 * jac[0][2] + 4 * t1 * jac[1][2]) +
858 n4.w[iw] * (4 * t3 * jac[0][2] + 4 * t1 * jac[2][2]) +
859 n5.w[iw] * (4 * t3 * jac[1][2] + 4 * t2 * jac[2][2])) *
860 invdet;
861 } else {
862 const Node& n6 = m_nodes[element.emap[6]];
863 const Node& n7 = m_nodes[element.emap[7]];
864 wx = -(n0.w[iw] * ((1 - t2) * (2 * t1 + t2) * jac[0][0] +
865 (1 - t1) * (t1 + 2 * t2) * jac[1][0]) *
866 0.25 +
867 n1.w[iw] * ((1 - t2) * (2 * t1 - t2) * jac[0][0] -
868 (1 + t1) * (t1 - 2 * t2) * jac[1][0]) *
869 0.25 +
870 n2.w[iw] * ((1 + t2) * (2 * t1 + t2) * jac[0][0] +
871 (1 + t1) * (t1 + 2 * t2) * jac[1][0]) *
872 0.25 +
873 n3.w[iw] * ((1 + t2) * (2 * t1 - t2) * jac[0][0] -
874 (1 - t1) * (t1 - 2 * t2) * jac[1][0]) *
875 0.25 +
876 n4.w[iw] * (t1 * (t2 - 1) * jac[0][0] +
877 (t1 - 1) * (t1 + 1) * jac[1][0] * 0.5) +
878 n5.w[iw] * ((1 - t2) * (1 + t2) * jac[0][0] * 0.5 -
879 (1 + t1) * t2 * jac[1][0]) +
880 n6.w[iw] * (-t1 * (1 + t2) * jac[0][0] +
881 (1 - t1) * (1 + t1) * jac[1][0] * 0.5) +
882 n7.w[iw] * ((t2 - 1) * (1 + t2) * jac[0][0] * 0.5 +
883 (t1 - 1) * t2 * jac[1][0])) *
884 invdet;
885 wy = -(n0.w[iw] * ((1 - t2) * (2 * t1 + t2) * jac[0][1] +
886 (1 - t1) * (t1 + 2 * t2) * jac[1][1]) *
887 0.25 +
888 n1.w[iw] * ((1 - t2) * (2 * t1 - t2) * jac[0][1] -
889 (1 + t1) * (t1 - 2 * t2) * jac[1][1]) *
890 0.25 +
891 n2.w[iw] * ((1 + t2) * (2 * t1 + t2) * jac[0][1] +
892 (1 + t1) * (t1 + 2 * t2) * jac[1][1]) *
893 0.25 +
894 n3.w[iw] * ((1 + t2) * (2 * t1 - t2) * jac[0][1] -
895 (1 - t1) * (t1 - 2 * t2) * jac[1][1]) *
896 0.25 +
897 n4.w[iw] * (t1 * (t2 - 1) * jac[0][1] +
898 (t1 - 1) * (t1 + 1) * jac[1][1] * 0.5) +
899 n5.w[iw] * ((1 - t2) * (1 + t2) * jac[0][1] * 0.5 -
900 (1 + t1) * t2 * jac[1][1]) +
901 n6.w[iw] * (-t1 * (1 + t2) * jac[0][1] +
902 (1 - t1) * (1 + t1) * jac[1][1] * 0.5) +
903 n7.w[iw] * ((t2 - 1) * (t2 + 1) * jac[0][1] * 0.5 +
904 (t1 - 1) * t2 * jac[1][1])) *
905 invdet;
906 }
907
908 // Transform field to global coordinates
909 UnmapFields(wx, wy, wz, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
910}
Garfield::ComponentFieldMap::GetWeightingFieldIndex
size_t GetWeightingFieldIndex(const std::string &label) const
Definition: ComponentFieldMap.cc:2348

◆ WeightingPotential()

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

Calculate the weighting potential at a given point.

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

Reimplemented from Garfield::Component.

Definition at line 912 of file ComponentAnsys121.cc.

914 {
915 // Do not proceed if not properly initialised.
916 if (!m_ready) return 0.;
917
918 // Look for the label.
919 const size_t iw = GetWeightingFieldIndex(label);
920 // Do not proceed if the requested weighting field does not exist.
921 if (iw == m_wfields.size()) return 0.;
922 // Check if the weighting field is properly initialised.
923 if (!m_wfieldsOk[iw]) return 0.;
924
925 // Copy the coordinates.
926 double x = xin, y = yin, z = zin;
927
928 // Map the coordinates onto field map coordinates.
929 bool xmirr, ymirr, zmirr;
930 double rcoordinate, rotation;
931 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
932
933 if (m_warning) PrintWarning("WeightingPotential");
934
935 // Find the element that contains this point.
936 double t1, t2, t3, t4, jac[4][4], det;
937 const int imap = FindElement5(x, y, z, t1, t2, t3, t4, jac, det);
938 // Check if the point is in the mesh
939 if (imap < 0) return 0.;
940
941 const Element& element = m_elements[imap];
942 if (m_debug) {
943 PrintElement("WeightingPotential", x, y, z, t1, t2, t3, t4, element, 8, iw);
944 }
945 // Calculate quadrilateral field, which can degenerate to a triangular field
946 const Node& n0 = m_nodes[element.emap[0]];
947 const Node& n1 = m_nodes[element.emap[1]];
948 const Node& n2 = m_nodes[element.emap[2]];
949 const Node& n3 = m_nodes[element.emap[3]];
950 const Node& n4 = m_nodes[element.emap[4]];
951 const Node& n5 = m_nodes[element.emap[5]];
952 if (element.degenerate) {
953 return n0.w[iw] * t1 * (2 * t1 - 1) + n1.w[iw] * t2 * (2 * t2 - 1) +
954 n2.w[iw] * t3 * (2 * t3 - 1) + 4 * n3.w[iw] * t1 * t2 +
955 4 * n4.w[iw] * t1 * t3 + 4 * n5.w[iw] * t2 * t3;
956 }
957
958 const Node& n6 = m_nodes[element.emap[6]];
959 const Node& n7 = m_nodes[element.emap[7]];
960 return -n0.w[iw] * (1 - t1) * (1 - t2) * (1 + t1 + t2) * 0.25 -
961 n1.w[iw] * (1 + t1) * (1 - t2) * (1 - t1 + t2) * 0.25 -
962 n2.w[iw] * (1 + t1) * (1 + t2) * (1 - t1 - t2) * 0.25 -
963 n3.w[iw] * (1 - t1) * (1 + t2) * (1 + t1 - t2) * 0.25 +
964 n4.w[iw] * (1 - t1) * (1 + t1) * (1 - t2) * 0.5 +
965 n5.w[iw] * (1 + t1) * (1 + t2) * (1 - t2) * 0.5 +
966 n6.w[iw] * (1 - t1) * (1 + t1) * (1 + t2) * 0.5 +
967 n7.w[iw] * (1 - t1) * (1 + t2) * (1 - t2) * 0.5;
968}
The documentation for this class was generated from the following files: