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

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

#include <ComponentAnsys123.hh>

+ Inheritance diagram for Garfield::ComponentAnsys123:

Public Member Functions

 ComponentAnsys123 ()
 Constructor.
 
 ~ComponentAnsys123 ()
 Destructor.
 
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
 
Medium * GetMedium (const double x, const double y, const double z) override
 Get the medium at a given location (x, y, z).
 
bool Initialise (std::string elist="ELIST.lis", std::string nlist="NLIST.lis", std::string mplist="MPLIST.lis", std::string prnsol="PRNSOL.lis", std::string unit="cm")
 
bool SetWeightingField (std::string prnsol, std::string label)
 
- Public Member Functions inherited from Garfield::ComponentFieldMap
 ComponentFieldMap ()=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 three-dimensional ANSYS field maps.

Definition at line 10 of file ComponentAnsys123.hh.

Constructor & Destructor Documentation

◆ ComponentAnsys123()

Garfield::ComponentAnsys123::ComponentAnsys123 ( )

Constructor.

Definition at line 10 of file ComponentAnsys123.cc.

10: ComponentFieldMap("Ansys123") {}
Garfield::ComponentFieldMap::ComponentFieldMap
ComponentFieldMap()=delete
Default constructor.

◆ ~ComponentAnsys123()

Garfield::ComponentAnsys123::~ComponentAnsys123 ( )
inline

Destructor.

Definition at line 15 of file ComponentAnsys123.hh.

15{}

Member Function Documentation

◆ ElectricField() [1/2]

void Garfield::ComponentAnsys123::ElectricField ( const double  x,
const double  y,
const double  z,
double &  ex,
double &  ey,
double &  ez,
double &  v,
Medium *&  m,
int &  status 
)
overridevirtual

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

Implements Garfield::Component.

Definition at line 733 of file ComponentAnsys123.cc.

736 {
737 // Copy the coordinates
738 double x = xin, y = yin, z = zin;
739
740 // Map the coordinates onto field map coordinates
741 bool xmirr, ymirr, zmirr;
742 double rcoordinate, rotation;
743 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
744
745 // Initial values
746 ex = ey = ez = volt = 0.;
747 status = 0;
748 m = nullptr;
749
750 // Do not proceed if not properly initialised.
751 if (!m_ready) {
752 status = -10;
753 PrintNotReady("ElectricField");
754 return;
755 }
756
757 if (m_warning) PrintWarning("ElectricField");
758
759 // Find the element that contains this point
760 double t1, t2, t3, t4, jac[4][4], det;
761 const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
762 if (imap < 0) {
763 if (m_debug) {
764 std::cerr << m_className << "::ElectricField: Point (" << x << ", " << y
765 << ", " << z << ") not in the mesh.\n";
766 }
767 status = -6;
768 return;
769 }
770
771 const Element& element = m_elements[imap];
772 if (m_debug) {
773 PrintElement("ElectricField", x, y, z, t1, t2, t3, t4, element, 10);
774 }
775 const Node& n0 = m_nodes[element.emap[0]];
776 const Node& n1 = m_nodes[element.emap[1]];
777 const Node& n2 = m_nodes[element.emap[2]];
778 const Node& n3 = m_nodes[element.emap[3]];
779 const Node& n4 = m_nodes[element.emap[4]];
780 const Node& n5 = m_nodes[element.emap[5]];
781 const Node& n6 = m_nodes[element.emap[6]];
782 const Node& n7 = m_nodes[element.emap[7]];
783 const Node& n8 = m_nodes[element.emap[8]];
784 const Node& n9 = m_nodes[element.emap[9]];
785 // Shorthands.
786 const double fourt1 = 4 * t1;
787 const double fourt2 = 4 * t2;
788 const double fourt3 = 4 * t3;
789 const double fourt4 = 4 * t4;
790 const double invdet = 1. / det;
791 // Tetrahedral field
792 volt = n0.v * t1 * (2 * t1 - 1) + n1.v * t2 * (2 * t2 - 1) +
793 n2.v * t3 * (2 * t3 - 1) + n3.v * t4 * (2 * t4 - 1) +
794 n4.v * fourt1 * t2 + n5.v * fourt1 * t3 + n6.v * fourt1 * t4 +
795 n7.v * fourt2 * t3 + n8.v * fourt2 * t4 + n9.v * fourt3 * t4;
796 ex = -(n0.v * (fourt1 - 1) * jac[0][1] + n1.v * (fourt2 - 1) * jac[1][1] +
797 n2.v * (fourt3 - 1) * jac[2][1] + n3.v * (fourt4 - 1) * jac[3][1] +
798 n4.v * (fourt2 * jac[0][1] + fourt1 * jac[1][1]) +
799 n5.v * (fourt3 * jac[0][1] + fourt1 * jac[2][1]) +
800 n6.v * (fourt4 * jac[0][1] + fourt1 * jac[3][1]) +
801 n7.v * (fourt3 * jac[1][1] + fourt2 * jac[2][1]) +
802 n8.v * (fourt4 * jac[1][1] + fourt2 * jac[3][1]) +
803 n9.v * (fourt4 * jac[2][1] + fourt3 * jac[3][1])) *
804 invdet;
805
806 ey = -(n0.v * (fourt1 - 1) * jac[0][2] + n1.v * (fourt2 - 1) * jac[1][2] +
807 n2.v * (fourt3 - 1) * jac[2][2] + n3.v * (fourt4 - 1) * jac[3][2] +
808 n4.v * (fourt2 * jac[0][2] + fourt1 * jac[1][2]) +
809 n5.v * (fourt3 * jac[0][2] + fourt1 * jac[2][2]) +
810 n6.v * (fourt4 * jac[0][2] + fourt1 * jac[3][2]) +
811 n7.v * (fourt3 * jac[1][2] + fourt2 * jac[2][2]) +
812 n8.v * (fourt4 * jac[1][2] + fourt2 * jac[3][2]) +
813 n9.v * (fourt4 * jac[2][2] + fourt3 * jac[3][2])) *
814 invdet;
815
816 ez = -(n0.v * (fourt1 - 1) * jac[0][3] + n1.v * (fourt2 - 1) * jac[1][3] +
817 n2.v * (fourt3 - 1) * jac[2][3] + n3.v * (fourt4 - 1) * jac[3][3] +
818 n4.v * (fourt2 * jac[0][3] + fourt1 * jac[1][3]) +
819 n5.v * (fourt3 * jac[0][3] + fourt1 * jac[2][3]) +
820 n6.v * (fourt4 * jac[0][3] + fourt1 * jac[3][3]) +
821 n7.v * (fourt3 * jac[1][3] + fourt2 * jac[2][3]) +
822 n8.v * (fourt4 * jac[1][3] + fourt2 * jac[3][3]) +
823 n9.v * (fourt4 * jac[2][3] + fourt3 * jac[3][3])) *
824 invdet;
825
826 // Transform field to global coordinates
827 UnmapFields(ex, ey, ez, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
828
829 // Drift medium?
830 if (m_debug) {
831 std::cout << m_className << "::ElectricField:\n";
832 std::cout << " Material " << element.matmap << ", drift flag "
833 << m_materials[element.matmap].driftmedium << ".\n";
834 }
835 m = m_materials[element.matmap].medium;
836 status = -5;
837 if (m_materials[element.matmap].driftmedium) {
838 if (m && m->IsDriftable()) status = 0;
839 }
840}
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::FindElement13
int FindElement13(const double x, const double y, const double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det)
Find the element for a point in curved quadratic tetrahedra.
Definition: ComponentFieldMap.cc:337
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::ComponentAnsys123::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 726 of file ComponentAnsys123.cc.

728 {
729 double v = 0.;
730 ElectricField(x, y, z, ex, ey, ez, v, m, status);
731}
Garfield::ComponentAnsys123::ElectricField
void ElectricField(const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status) override
Definition: ComponentAnsys123.cc:726

Referenced by ElectricField().

◆ GetAspectRatio()

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

Implements Garfield::ComponentFieldMap.

Definition at line 1049 of file ComponentAnsys123.cc.

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

◆ GetElementVolume()

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

Implements Garfield::ComponentFieldMap.

Definition at line 1031 of file ComponentAnsys123.cc.

1031 {
1032 if (i >= m_elements.size()) return 0.;
1033 const Element& element = m_elements[i];
1034 const Node& n0 = m_nodes[element.emap[0]];
1035 const Node& n1 = m_nodes[element.emap[1]];
1036 const Node& n2 = m_nodes[element.emap[2]];
1037 const Node& n3 = m_nodes[element.emap[3]];
1038 const double vol =
1039 fabs((n3.x - n0.x) *
1040 ((n1.y - n0.y) * (n2.z - n0.z) - (n2.y - n0.y) * (n1.z - n0.z)) +
1041 (n3.y - n0.y) *
1042 ((n1.z - n0.z) * (n2.x - n0.x) - (n2.z - n0.z) * (n1.x - n0.x)) +
1043 (n3.z - n0.z) * ((n1.x - n0.x) * (n2.y - n0.y) -
1044 (n3.x - n0.x) * (n1.y - n0.y))) /
1045 6.;
1046 return vol;
1047}
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615

◆ GetMedium()

Medium * Garfield::ComponentAnsys123::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 984 of file ComponentAnsys123.cc.

985 {
986 // Copy the coordinates.
987 double x = xin, y = yin, z = zin;
988
989 // Map the coordinates onto field map coordinates.
990 bool xmirr, ymirr, zmirr;
991 double rcoordinate, rotation;
992 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
993
994 // Do not proceed if not properly initialised.
995 if (!m_ready) {
996 std::cerr << m_className << "::GetMedium:\n";
997 std::cerr << " Field map not available for interpolation.\n";
998 return nullptr;
999 }
1000 if (m_warning) PrintWarning("GetMedium");
1001
1002 // Find the element that contains this point.
1003 double t1, t2, t3, t4, jac[4][4], det;
1004 const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
1005 if (imap < 0) {
1006 if (m_debug) {
1007 std::cerr << m_className << "::GetMedium:\n";
1008 std::cerr << " Point (" << x << ", " << y << ", " << z
1009 << ") not in the mesh.\n";
1010 }
1011 return nullptr;
1012 }
1013 const Element& element = m_elements[imap];
1014 if (element.matmap >= m_materials.size()) {
1015 if (m_debug) {
1016 std::cerr << m_className << "::GetMedium:\n";
1017 std::cerr << " Point (" << x << ", " << y << ", " << z << ")"
1018 << " has out of range material number " << imap << ".\n";
1019 }
1020 return nullptr;
1021 }
1022
1023 if (m_debug) {
1024 PrintElement("GetMedium", x, y, z, t1, t2, t3, t4, element, 10);
1025 }
1026
1027 // Assign a medium.
1028 return m_materials[element.matmap].medium;
1029}

◆ Initialise()

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

Definition at line 12 of file ComponentAnsys123.cc.

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

◆ SetWeightingField()

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

Definition at line 606 of file ComponentAnsys123.cc.

607 {
608 if (!m_ready) {
609 PrintNotReady("SetWeightingField");
610 std::cerr << " Weighting field cannot be added.\n";
611 return false;
612 }
613
614 // Open the voltage list.
615 std::ifstream fprnsol;
616 fprnsol.open(prnsol.c_str(), std::ios::in);
617 if (fprnsol.fail()) {
618 PrintCouldNotOpen("SetWeightingField", prnsol);
619 return false;
620 }
621
622 // Check if a weighting field with the same label already exists.
623 const size_t iw = GetOrCreateWeightingFieldIndex(label);
624 if (iw + 1 != m_wfields.size()) {
625 std::cout << m_className << "::SetWeightingField:\n"
626 << " Replacing existing weighting field " << label << ".\n";
627 }
628 m_wfieldsOk[iw] = false;
629
630 // Buffer for reading
631 constexpr int size = 100;
632 char line[size];
633
634 bool ok = true;
635 // Read the voltage list.
636 int il = 0;
637 unsigned int nread = 0;
638 bool readerror = false;
639 while (fprnsol.getline(line, size, '\n')) {
640 il++;
641 // Skip page feed
642 if (strcmp(line, "1") == 0) {
643 fprnsol.getline(line, size, '\n');
644 il++;
645 fprnsol.getline(line, size, '\n');
646 il++;
647 fprnsol.getline(line, size, '\n');
648 il++;
649 fprnsol.getline(line, size, '\n');
650 il++;
651 fprnsol.getline(line, size, '\n');
652 il++;
653 continue;
654 }
655 // Skip page feed (Ansys > v15.x).
656 if (strstr(line, "***") != NULL) {
657 fprnsol.getline(line, size, '\n');
658 il++;
659 fprnsol.getline(line, size, '\n');
660 il++;
661 fprnsol.getline(line, size, '\n');
662 il++;
663 continue;
664 }
665 // Split the line in tokens.
666 char* token = NULL;
667 token = strtok(line, " ");
668 // Skip blank lines and headers.
669 if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
670 int(token[0]) == 10 || int(token[0]) == 13 ||
671 strcmp(token, "PRINT") == 0 || strcmp(token, "*****") == 0 ||
672 strcmp(token, "LOAD") == 0 || strcmp(token, "TIME=") == 0 ||
673 strcmp(token, "MAXIMUM") == 0 || strcmp(token, "VALUE") == 0 ||
674 strcmp(token, "NODE") == 0)
675 continue;
676 // Read the element.
677 int inode = ReadInteger(token, -1, readerror);
678 token = strtok(NULL, " ");
679 double volt = ReadDouble(token, -1, readerror);
680 // Check the syntax.
681 if (readerror) {
682 std::cerr << m_className << "::SetWeightingField:\n";
683 std::cerr << " Error reading file " << prnsol.c_str() << " (line "
684 << il << ").\n";
685 fprnsol.close();
686 return false;
687 }
688 // Check node number and store if OK.
689 if (inode < 1 || inode > (int)m_nodes.size()) {
690 std::cerr << m_className << "::SetWeightingField:\n";
691 std::cerr << " Node number " << inode << " out of range\n";
692 std::cerr << " on potential file " << prnsol.c_str() << " (line " << il
693 << ").\n";
694 ok = false;
695 } else {
696 m_nodes[inode - 1].w[iw] = volt;
697 nread++;
698 }
699 }
700 // Close the file.
701 fprnsol.close();
702
703 std::cout << m_className << "::SetWeightingField:\n"
704 << " Read " << nread << " potentials from file "
705 << prnsol << ".\n";
706 // Check the number of nodes.
707 if (nread != m_nodes.size()) {
708 std::cerr << m_className << "::SetWeightingField:\n"
709 << " Number of nodes read from potential file " << prnsol
710 << " (" << nread << ")\n does not match the node list ("
711 << m_nodes.size() << ").\n";
712 ok = false;
713 }
714
715 // Set the ready flag.
716 m_wfieldsOk[iw] = ok;
717 if (!ok) {
718 std::cerr << m_className << "::SetWeightingField:\n";
719 std::cerr << " Field map could not be read "
720 << "and cannot be interpolated.\n";
721 return false;
722 }
723 return true;
724}
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::ComponentAnsys123::UpdatePeriodicity ( )
inlineoverrideprotectedvirtual

Verify periodicities.

Implements Garfield::Component.

Definition at line 41 of file ComponentAnsys123.hh.

41{ UpdatePeriodicityCommon(); }

◆ WeightingField()

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

844 {
845 // Initial values
846 wx = wy = wz = 0;
847
848 // Do not proceed if not properly initialised.
849 if (!m_ready) return;
850
851 // Look for the label.
852 const size_t iw = GetWeightingFieldIndex(label);
853 // Do not proceed if the requested weighting field does not exist.
854 if (iw == m_wfields.size()) return;
855 // Check if the weighting field is properly initialised.
856 if (!m_wfieldsOk[iw]) return;
857
858 // Copy the coordinates.
859 double x = xin, y = yin, z = zin;
860
861 // Map the coordinates onto field map coordinates
862 bool xmirr, ymirr, zmirr;
863 double rcoordinate, rotation;
864 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
865
866 if (m_warning) PrintWarning("WeightingField");
867
868 // Find the element that contains this point.
869 double t1, t2, t3, t4, jac[4][4], det;
870 const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
871 // Check if the point is in the mesh.
872 if (imap < 0) return;
873
874 const Element& element = m_elements[imap];
875 if (m_debug) {
876 PrintElement("WeightingField", x, y, z, t1, t2, t3, t4, element, 10, iw);
877 }
878 const Node& n0 = m_nodes[element.emap[0]];
879 const Node& n1 = m_nodes[element.emap[1]];
880 const Node& n2 = m_nodes[element.emap[2]];
881 const Node& n3 = m_nodes[element.emap[3]];
882 const Node& n4 = m_nodes[element.emap[4]];
883 const Node& n5 = m_nodes[element.emap[5]];
884 const Node& n6 = m_nodes[element.emap[6]];
885 const Node& n7 = m_nodes[element.emap[7]];
886 const Node& n8 = m_nodes[element.emap[8]];
887 const Node& n9 = m_nodes[element.emap[9]];
888 // Tetrahedral field
889 const double invdet = 1. / det;
890 const double fourt1 = 4 * t1;
891 const double fourt2 = 4 * t2;
892 const double fourt3 = 4 * t3;
893 const double fourt4 = 4 * t4;
894 wx = -(n0.w[iw] * (fourt1 - 1) * jac[0][1] +
895 n1.w[iw] * (fourt2 - 1) * jac[1][1] +
896 n2.w[iw] * (fourt3 - 1) * jac[2][1] +
897 n3.w[iw] * (fourt4 - 1) * jac[3][1] +
898 n4.w[iw] * (fourt2 * jac[0][1] + fourt1 * jac[1][1]) +
899 n5.w[iw] * (fourt3 * jac[0][1] + fourt1 * jac[2][1]) +
900 n6.w[iw] * (fourt4 * jac[0][1] + fourt1 * jac[3][1]) +
901 n7.w[iw] * (fourt3 * jac[1][1] + fourt2 * jac[2][1]) +
902 n8.w[iw] * (fourt4 * jac[1][1] + fourt2 * jac[3][1]) +
903 n9.w[iw] * (fourt4 * jac[2][1] + fourt3 * jac[3][1])) *
904 invdet;
905
906 wy = -(n0.w[iw] * (fourt1 - 1) * jac[0][2] +
907 n1.w[iw] * (fourt2 - 1) * jac[1][2] +
908 n2.w[iw] * (fourt3 - 1) * jac[2][2] +
909 n3.w[iw] * (fourt4 - 1) * jac[3][2] +
910 n4.w[iw] * (fourt2 * jac[0][2] + fourt1 * jac[1][2]) +
911 n5.w[iw] * (fourt3 * jac[0][2] + fourt1 * jac[2][2]) +
912 n6.w[iw] * (fourt4 * jac[0][2] + fourt1 * jac[3][2]) +
913 n7.w[iw] * (fourt3 * jac[1][2] + fourt2 * jac[2][2]) +
914 n8.w[iw] * (fourt4 * jac[1][2] + fourt2 * jac[3][2]) +
915 n9.w[iw] * (fourt4 * jac[2][2] + fourt3 * jac[3][2])) *
916 invdet;
917
918 wz = -(n0.w[iw] * (fourt1 - 1) * jac[0][3] +
919 n1.w[iw] * (fourt2 - 1) * jac[1][3] +
920 n2.w[iw] * (fourt3 - 1) * jac[2][3] +
921 n3.w[iw] * (fourt4 - 1) * jac[3][3] +
922 n4.w[iw] * (fourt2 * jac[0][3] + fourt1 * jac[1][3]) +
923 n5.w[iw] * (fourt3 * jac[0][3] + fourt1 * jac[2][3]) +
924 n6.w[iw] * (fourt4 * jac[0][3] + fourt1 * jac[3][3]) +
925 n7.w[iw] * (fourt3 * jac[1][3] + fourt2 * jac[2][3]) +
926 n8.w[iw] * (fourt4 * jac[1][3] + fourt2 * jac[3][3]) +
927 n9.w[iw] * (fourt4 * jac[2][3] + fourt3 * jac[3][3])) *
928 invdet;
929
930 // Transform field to global coordinates
931 UnmapFields(wx, wy, wz, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
932}
Garfield::ComponentFieldMap::GetWeightingFieldIndex
size_t GetWeightingFieldIndex(const std::string &label) const
Definition: ComponentFieldMap.cc:2348

◆ WeightingPotential()

double Garfield::ComponentAnsys123::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 934 of file ComponentAnsys123.cc.

936 {
937 // Do not proceed if not properly initialised.
938 if (!m_ready) return 0.;
939
940 // Look for the label.
941 const size_t iw = GetWeightingFieldIndex(label);
942 // Do not proceed if the requested weighting field does not exist.
943 if (iw == m_wfields.size()) return 0.;
944 // Check if the weighting field is properly initialised.
945 if (!m_wfieldsOk[iw]) return 0.;
946
947 // Copy the coordinates.
948 double x = xin, y = yin, z = zin;
949
950 // Map the coordinates onto field map coordinates.
951 bool xmirr, ymirr, zmirr;
952 double rcoordinate, rotation;
953 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
954
955 if (m_warning) PrintWarning("WeightingPotential");
956
957 // Find the element that contains this point.
958 double t1, t2, t3, t4, jac[4][4], det;
959 const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
960 if (imap < 0) return 0.;
961
962 const Element& element = m_elements[imap];
963 if (m_debug) {
964 PrintElement("WeightingPotential", x, y, z, t1, t2, t3, t4, element, 10,
965 iw);
966 }
967 const Node& n0 = m_nodes[element.emap[0]];
968 const Node& n1 = m_nodes[element.emap[1]];
969 const Node& n2 = m_nodes[element.emap[2]];
970 const Node& n3 = m_nodes[element.emap[3]];
971 const Node& n4 = m_nodes[element.emap[4]];
972 const Node& n5 = m_nodes[element.emap[5]];
973 const Node& n6 = m_nodes[element.emap[6]];
974 const Node& n7 = m_nodes[element.emap[7]];
975 const Node& n8 = m_nodes[element.emap[8]];
976 const Node& n9 = m_nodes[element.emap[9]];
977 return n0.w[iw] * t1 * (2 * t1 - 1) + n1.w[iw] * t2 * (2 * t2 - 1) +
978 n2.w[iw] * t3 * (2 * t3 - 1) + n3.w[iw] * t4 * (2 * t4 - 1) +
979 4 * n4.w[iw] * t1 * t2 + 4 * n5.w[iw] * t1 * t3 +
980 4 * n6.w[iw] * t1 * t4 + 4 * n7.w[iw] * t2 * t3 +
981 4 * n8.w[iw] * t2 * t4 + 4 * n9.w[iw] * t3 * t4;
982}
The documentation for this class was generated from the following files: