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)
void	SetRangeZ (const double zmin, const double zmax)
Public Member Functions inherited from Garfield::ComponentFieldMap
	ComponentFieldMap ()
	Constructor.
virtual	~ComponentFieldMap ()
	Destructor.
virtual void	SetRange ()
	Calculate x, y, z, V and angular ranges.
void	PrintRange ()
	Show x, y, z, V and angular ranges.
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	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.
unsigned int	GetNumberOfMaterials () const
	Return the number of materials in the field map.
double	GetPermittivity (const unsigned int imat) const
	Return the permittivity of a field map material.
double	GetConductivity (const unsigned int imat) const
	Return the conductivity of a field map material.
void	SetMedium (const unsigned int imat, Medium *medium)
	Associate a field map material with a Medium class.
Medium *	GetMedium (const unsigned int i) const
	Return the Medium associated to a field map material.
unsigned int	GetNumberOfMedia () const
int	GetNumberOfElements () const
	Return the number of mesh elements.
bool	GetElement (const unsigned int i, double &vol, double &dmin, double &dmax)
	Return the volume and aspect ratio of a mesh element.
void	EnableCheckMapIndices ()
void	DisableCheckMapIndices ()
void	EnableDeleteBackgroundElements (const bool on=true)
	Option to eliminate mesh elements in conductors (default: on).
void	EnableTetrahedralTreeForElementSearch (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::ComponentBase
	ComponentBase ()
	Constructor.
virtual	~ComponentBase ()
	Destructor.
virtual void	SetGeometry (GeometryBase *geo)
	Define the geometry.
virtual void	Clear ()
	Reset.
virtual Medium *	GetMedium (const double x, const double y, const double z)
	Get the medium at a given location (x, y, z).
virtual void	ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status)=0
virtual void	ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status)=0
	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 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.
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	IntegrateFlux (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)
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 direction.
void	DisablePeriodicityX ()
void	EnablePeriodicityY (const bool on=true)
	Enable simple periodicity in the direction.
void	DisablePeriodicityY ()
void	EnablePeriodicityZ (const bool on=true)
	Enable simple periodicity in the direction.
void	DisablePeriodicityZ ()
void	EnableMirrorPeriodicityX (const bool on=true)
	Enable mirror periodicity in the direction.
void	DisableMirrorPeriodicityX ()
void	EnableMirrorPeriodicityY (const bool on=true)
	Enable mirror periodicity in the direction.
void	DisableMirrorPeriodicityY ()
void	EnableMirrorPeriodicityZ (const bool on=true)
	Enable mirror periodicity in the direction.
void	DisableMirrorPeriodicityZ ()
void	EnableAxialPeriodicityX (const bool on=true)
	Enable axial periodicity in the direction.
void	DisableAxialPeriodicityX ()
void	EnableAxialPeriodicityY (const bool on=true)
	Enable axial periodicity in the direction.
void	DisableAxialPeriodicityY ()
void	EnableAxialPeriodicityZ (const bool on=true)
	Enable axial periodicity in the direction.
void	DisableAxialPeriodicityZ ()
void	EnableRotationSymmetryX (const bool on=true)
	Enable rotation symmetry around the axis.
void	DisableRotationSymmetryX ()
void	EnableRotationSymmetryY (const bool on=true)
	Enable rotation symmetry around the axis.
void	DisableRotationSymmetryY ()
void	EnableRotationSymmetryZ (const bool on=true)
	Enable rotation symmetry around the axis.
void	DisableRotationSymmetryZ ()
void	EnableDebugging ()
	Switch on debugging messages.
void	DisableDebugging ()
	Switch off debugging messages.
void	ActivateTraps ()
	Request trapping to be taken care of by the component (for TCAD).
void	DeactivateTraps ()
bool	IsTrapActive ()
void	ActivateVelocityMap ()
	Request velocity to be taken care of by the component (for TCAD).
void	DectivateVelocityMap ()
bool	IsVelocityActive ()
virtual bool	ElectronAttachment (const double, const double, const double, double &eta)
	Get the electron attachment coefficient.
virtual bool	HoleAttachment (const double, const double, const double, double &eta)
	Get the hole attachment coefficient.
virtual void	ElectronVelocity (const double, const double, const double, double &vx, double &vy, double &vz, Medium *&, int &status)
	Get the electron drift velocity.
virtual void	HoleVelocity (const double, const double, const double, double &vx, double &vy, double &vz, Medium *&, int &status)
	Get the hole drift velocity.
virtual bool	GetElectronLifetime (const double, const double, const double, double &etau)
virtual bool	GetHoleLifetime (const double, const double, const double, double &htau)

Protected Member Functions
void	UpdatePeriodicity () 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	UpdatePeriodicity2d ()
void	UpdatePeriodicityCommon ()
int	FindElement5 (const double x, const double y, const double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det)
	Find the element for a point in curved quadratic quadrilaterals.
int	FindElement13 (const double x, const double y, const double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det)
	Find the element for a point in curved quadratic tetrahedra.
int	FindElementCube (const double x, const double y, const double z, double &t1, double &t2, double &t3, TMatrixD *&jac, std::vector< TMatrixD * > &dN)
	Find the element for a point in a cube.
void	MapCoordinates (double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
	Move (xpos, ypos, zpos) to field map coordinates.
void	UnmapFields (double &ex, double &ey, double &ez, double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
	Move (ex, ey, ez) to global coordinates.
int	ReadInteger (char *token, int def, bool &error)
double	ReadDouble (char *token, double def, bool &error)
virtual double	GetElementVolume (const unsigned int i)=0
virtual void	GetAspectRatio (const unsigned int i, double &dmin, double &dmax)=0
void	PrintWarning (const std::string &header)
void	PrintNotReady (const std::string &header) const
void	PrintElement (const std::string &header, const double x, const double y, const double z, const double t1, const double t2, const double t3, const double t4, const 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
Protected Attributes inherited from Garfield::ComponentFieldMap
bool	m_is3d = true
int	nElements = -1
std::vector< Element >	elements
int	nNodes = -1
std::vector< Node >	nodes
unsigned int	m_nMaterials = 0
std::vector< Material >	materials
int	nWeightingFields = 0
std::vector< std::string >	wfields
std::vector< bool >	wfieldsOk
bool	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
double	m_mapvmax
std::array< bool, 3 >	m_setang
bool	m_deleteBackground = true
bool	m_warning = false
unsigned int	m_nWarnings = 0
Protected Attributes inherited from Garfield::ComponentBase
std::string	m_className = "ComponentBase"
	Class name.
GeometryBase *	m_geometry = nullptr
	Pointer to the geometry.
bool	m_ready = false
	Ready for use?
bool	m_activeTraps = false
	Does the component have traps?
bool	m_hasVelocityMap = false
	Does the component have velocity maps?
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.
double	m_bx0 = 0.
double	m_by0 = 0.
double	m_bz0 = 0.
bool	m_debug = false
	Switch on/off debugging messages.

Detailed Description

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

Definition at line 10 of file ComponentAnsys121.hh.

Constructor & Destructor Documentation

ComponentAnsys121()

Garfield::ComponentAnsys121::ComponentAnsys121

(

)

Constructor.

Definition at line 10 of file ComponentAnsys121.cc.

                                     : ComponentFieldMap() {
  m_className = "ComponentAnsys121";
  m_is3d = false;
  // Default bounding box
  m_minBoundingBox[2] = -50;
  m_maxBoundingBox[2] = 50;
}

~ComponentAnsys121()

Garfield::ComponentAnsys121::~ComponentAnsys121 ( )

inline

Destructor.

Definition at line 15 of file ComponentAnsys121.hh.

{}

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::ComponentBase.

Definition at line 712 of file ComponentAnsys121.cc.

                                                   {
  // Copy the coordinates
  double x = xin, y = yin, z = 0.;
 
  // Map the coordinates onto field map coordinates
  bool xmirr, ymirr, zmirr;
  double rcoordinate, rotation;
  MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
 
  // Initial values
  ex = ey = ez = volt = 0;
  m = nullptr;
  status = 0;
 
  // Do not proceed if not properly initialised.
  if (!m_ready) {
    status = -10;
    PrintNotReady("ElectricField");
    return;
  }
 
  if (m_warning) PrintWarning("ElectricField");
 
  if (zin < m_minBoundingBox[2] || zin > m_maxBoundingBox[2]) {
    status = -5;
    return;
  }
 
  // Find the element that contains this point
  double t1, t2, t3, t4, jac[4][4], det;
  const int imap = FindElement5(x, y, z, t1, t2, t3, t4, jac, det);
  if (imap < 0) {
    if (m_debug) {
      std::cout << m_className << "::ElectricField:\n";
      std::cout << "    Point (" << x << ", " << y << ") not in the mesh.\n";
    }
    status = -6;
    return;
  }
 
  const Element& element = elements[imap];
  if (m_debug) {
    PrintElement("ElectricField", x, y, z, t1, t2, t3, t4, element, 8);
  }
 
  const Node& n0 = nodes[element.emap[0]];
  const Node& n1 = nodes[element.emap[1]];
  const Node& n2 = nodes[element.emap[2]];
  const Node& n3 = nodes[element.emap[3]];
  const Node& n4 = nodes[element.emap[4]];
  const Node& n5 = nodes[element.emap[5]];
  // Calculate quadrilateral field, which can degenerate to a triangular field
  const double invdet = 1. / det;
  if (element.degenerate) {
    volt = n0.v * t1 * (2 * t1 - 1) + n1.v * t2 * (2 * t2 - 1) +
           n2.v * t3 * (2 * t3 - 1) + 4 * n3.v * t1 * t2 + 4 * n4.v * t1 * t3 +
           4 * n5.v * t2 * t3;
    ex = -(n0.v * (4 * t1 - 1) * jac[0][1] + n1.v * (4 * t2 - 1) * jac[1][1] +
           n2.v * (4 * t3 - 1) * jac[2][1] +
           n3.v * (4 * t2 * jac[0][1] + 4 * t1 * jac[1][1]) +
           n4.v * (4 * t3 * jac[0][1] + 4 * t1 * jac[2][1]) +
           n5.v * (4 * t3 * jac[1][1] + 4 * t2 * jac[2][1])) *
         invdet;
    ey = -(n0.v * (4 * t1 - 1) * jac[0][2] + n1.v * (4 * t2 - 1) * jac[1][2] +
           n2.v * (4 * t3 - 1) * jac[2][2] +
           n3.v * (4 * t2 * jac[0][2] + 4 * t1 * jac[1][2]) +
           n4.v * (4 * t3 * jac[0][2] + 4 * t1 * jac[2][2]) +
           n5.v * (4 * t3 * jac[1][2] + 4 * t2 * jac[2][2])) *
         invdet;
  } else {
    const Node& n6 = nodes[element.emap[6]];
    const Node& n7 = nodes[element.emap[7]];
    volt = -n0.v * (1 - t1) * (1 - t2) * (1 + t1 + t2) * 0.25 -
           n1.v * (1 + t1) * (1 - t2) * (1 - t1 + t2) * 0.25 -
           n2.v * (1 + t1) * (1 + t2) * (1 - t1 - t2) * 0.25 -
           n3.v * (1 - t1) * (1 + t2) * (1 + t1 - t2) * 0.25 +
           n4.v * (1 - t1) * (1 + t1) * (1 - t2) * 0.5 +
           n5.v * (1 + t1) * (1 + t2) * (1 - t2) * 0.5 +
           n6.v * (1 - t1) * (1 + t1) * (1 + t2) * 0.5 +
           n7.v * (1 - t1) * (1 + t2) * (1 - t2) * 0.5;
    ex = -(n0.v * ((1 - t2) * (2 * t1 + t2) * jac[0][0] +
                   (1 - t1) * (t1 + 2 * t2) * jac[1][0]) *
               0.25 +
           n1.v * ((1 - t2) * (2 * t1 - t2) * jac[0][0] -
                   (1 + t1) * (t1 - 2 * t2) * jac[1][0]) *
               0.25 +
           n2.v * ((1 + t2) * (2 * t1 + t2) * jac[0][0] +
                   (1 + t1) * (t1 + 2 * t2) * jac[1][0]) *
               0.25 +
           n3.v * ((1 + t2) * (2 * t1 - t2) * jac[0][0] -
                   (1 - t1) * (t1 - 2 * t2) * jac[1][0]) *
               0.25 +
           n4.v * (t1 * (t2 - 1) * jac[0][0] +
                   (t1 - 1) * (t1 + 1) * jac[1][0] * 0.5) +
           n5.v * ((1 - t2) * (1 + t2) * jac[0][0] * 0.5 -
                   (1 + t1) * t2 * jac[1][0]) +
           n6.v * (-t1 * (1 + t2) * jac[0][0] +
                   (1 - t1) * (1 + t1) * jac[1][0] * 0.5) +
           n7.v * ((t2 - 1) * (t2 + 1) * jac[0][0] * 0.5 +
                   (t1 - 1) * t2 * jac[1][0])) *
         invdet;
    ey = -(n0.v * ((1 - t2) * (2 * t1 + t2) * jac[0][1] +
                   (1 - t1) * (t1 + 2 * t2) * jac[1][1]) *
               0.25 +
           n1.v * ((1 - t2) * (2 * t1 - t2) * jac[0][1] -
                   (1 + t1) * (t1 - 2 * t2) * jac[1][1]) *
               0.25 +
           n2.v * ((1 + t2) * (2 * t1 + t2) * jac[0][1] +
                   (1 + t1) * (t1 + 2 * t2) * jac[1][1]) *
               0.25 +
           n3.v * ((1 + t2) * (2 * t1 - t2) * jac[0][1] -
                   (1 - t1) * (t1 - 2 * t2) * jac[1][1]) *
               0.25 +
           n4.v * (t1 * (t2 - 1) * jac[0][1] +
                   (t1 - 1) * (t1 + 1) * jac[1][1] * 0.5) +
           n5.v * ((1 - t2) * (1 + t2) * jac[0][1] * 0.5 -
                   (1 + t1) * t2 * jac[1][1]) +
           n6.v * (-t1 * (1 + t2) * jac[0][1] +
                   (1 - t1) * (1 + t1) * jac[1][1] * 0.5) +
           n7.v * ((t2 - 1) * (t2 + 1) * jac[0][1] * 0.5 +
                   (t1 - 1) * t2 * jac[1][1])) *
         invdet;
  }
 
  // Transform field to global coordinates
  UnmapFields(ex, ey, ez, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
 
  // Drift medium?
  if (m_debug) {
    std::cout << m_className << "::ElectricField:\n";
    std::cout << "    Material " << element.matmap << ", drift flag "
              << materials[element.matmap].driftmedium << ".\n";
  }
  m = materials[element.matmap].medium;
  status = -5;
  if (materials[element.matmap].driftmedium) {
    if (m && m->IsDriftable()) status = 0;
  }
}

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,z	coordinates [cm].
ex,ey,ez	components of the electric field [V/cm].
m	pointer to the medium at this location.
status	status 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::ComponentBase.

Definition at line 705 of file ComponentAnsys121.cc.

                                                                           {
  double v;
  ElectricField(x, y, z, ex, ey, ez, v, m, status);
}

Referenced by ElectricField().

GetAspectRatio()

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

overrideprotectedvirtual

Implements Garfield::ComponentFieldMap.

Definition at line 1118 of file ComponentAnsys121.cc.

                                                     {
  if (i >= elements.size()) {
    dmin = dmax = 0.;
    return;
  }
 
  const Element& element = elements[i];
  const int np = 8;
  // Loop over all pairs of vertices.
  for (int j = 0; j < np - 1; ++j) {
    const Node& nj = nodes[element.emap[j]];
    for (int k = j + 1; k < np; ++k) {
      const Node& nk = nodes[element.emap[k]];
      // Compute distance.
      const double dx = nj.x - nk.x;
      const double dy = nj.y - nk.y;
      const double dist = sqrt(dx * dx + dy * dy);
      if (k == 1) {
        dmin = dmax = dist;
      } else {
        if (dist < dmin) dmin = dist;
        if (dist > dmax) dmax = dist;
      }
    }
  }
}

GetElementVolume()

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

overrideprotectedvirtual

Implements Garfield::ComponentFieldMap.

Definition at line 1104 of file ComponentAnsys121.cc.

                                                               {
  if (i >= elements.size()) return 0.;
  const Element& element = elements[i];
  const Node& n0 = nodes[element.emap[0]];
  const Node& n1 = nodes[element.emap[1]];
  const Node& n2 = nodes[element.emap[2]];
  const Node& n3 = nodes[element.emap[3]];
  const double surf =
      0.5 *
      (fabs((n1.x - n0.x) * (n2.y - n0.y) - (n2.x - n0.x) * (n1.y - n0.y)) +
       fabs((n3.x - n0.x) * (n2.y - n0.y) - (n2.x - n0.x) * (n3.y - n0.y)));
  return surf;
}

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::ComponentBase.

Definition at line 1041 of file ComponentAnsys121.cc.

                                                       {
  // Copy the coordinates.
  double x = xin, y = yin, z = 0.;
 
  // Map the coordinates onto field map coordinates.
  bool xmirr, ymirr, zmirr;
  double rcoordinate, rotation;
  MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
 
  if (zin < m_minBoundingBox[2] || z > m_maxBoundingBox[2]) {
    return nullptr;
  }
 
  // Do not proceed if not properly initialised.
  if (!m_ready) {
    PrintNotReady("GetMedium");
    return nullptr;
  }
  if (m_warning) PrintWarning("GetMedium");
 
  // Find the element that contains this point.
  double t1, t2, t3, t4, jac[4][4], det;
  const int imap = FindElement5(x, y, z, t1, t2, t3, t4, jac, det);
  if (imap < 0) {
    if (m_debug) {
      std::cerr << m_className << "::GetMedium:\n";
      std::cerr << "    Point (" << x << ", " << y << ") not in the mesh.\n";
    }
    return nullptr;
  }
  const Element& element = elements[imap];
  if (element.matmap >= m_nMaterials) {
    if (m_debug) {
      std::cerr << m_className << "::GetMedium:\n";
      std::cerr << "    Point (" << x << ", " << y << ")"
                << " has out of range material number " << imap << ".\n";
    }
    return nullptr;
  }
 
  if (m_debug) {
    PrintElement("GetMedium", x, y, z, t1, t2, t3, t4, element, 8);
  }
 
  // Assign a medium.
  return materials[element.matmap].medium;
}

Initialise()

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

Definition at line 18 of file ComponentAnsys121.cc.

                                                   {
  m_ready = false;
  m_warning = false;
  m_nWarnings = 0;
  // Keep track of the success.
  bool ok = true;
 
  // Buffer for reading
  const int size = 100;
  char line[size];
 
  // Open the material list.
  std::ifstream fmplist;
  fmplist.open(mplist.c_str(), std::ios::in);
  if (fmplist.fail()) {
    std::cerr << m_className << "::Initialise:\n"
              << "    Could not open material file " << mplist
              << " for reading. The file perhaps does not exist.\n";
    return false;
  }
 
  // Read the material list.
  m_nMaterials = 0;
  int il = 0;
  unsigned int icurrmat = 0;
  bool readerror = false;
  while (fmplist.getline(line, size, '\n')) {
    il++;
    // Skip page feed
    if (strcmp(line, "1") == 0) {
      fmplist.getline(line, size, '\n');
      il++;
      fmplist.getline(line, size, '\n');
      il++;
      fmplist.getline(line, size, '\n');
      il++;
      fmplist.getline(line, size, '\n');
      il++;
      fmplist.getline(line, size, '\n');
      il++;
      continue;
    }
    // Split the line in tokens
    char* token = NULL;
    token = strtok(line, " ");
    // Skip blank lines and headers
    if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
        strcmp(token, "TEMPERATURE") == 0 || strcmp(token, "PROPERTY=") == 0 ||
        int(token[0]) == 10 || int(token[0]) == 13)
      continue;
    // Read number of materials,
    // ensure it does not exceed the maximum and initialise the list
    if (strcmp(token, "LIST") == 0) {
      token = strtok(NULL, " ");
      token = strtok(NULL, " ");
      token = strtok(NULL, " ");
      token = strtok(NULL, " ");
      m_nMaterials = ReadInteger(token, -1, readerror);
      if (readerror) {
        std::cerr << m_className << "::Initialise:\n";
        std::cerr << "    Error reading file " << mplist << " (line " << il
                  << ").\n";
        fmplist.close();
        ok = false;
        return false;
      }
      materials.resize(m_nMaterials);
      for (unsigned int i = 0; i < m_nMaterials; ++i) {
        materials[i].ohm = -1;
        materials[i].eps = -1;
        materials[i].medium = nullptr;
      }
      if (m_debug) {
        std::cout << m_className << "::Initialise:\n";
        std::cout << "    Number of materials: " << m_nMaterials << "\n";
      }
    } else if (strcmp(token, "MATERIAL") == 0) {
      // Version 12 format: read material number
      token = strtok(NULL, " ");
      token = strtok(NULL, " ");
      const int imat = ReadInteger(token, -1, readerror);
      if (readerror || imat < 0) {
        std::cerr << m_className << "::Initialise:\n";
        std::cerr << "    Error reading file " << mplist << " (line " << il
                  << ").\n";
        fmplist.close();
        ok = false;
        return false;
      }
      icurrmat = imat;
    } else if (strcmp(token, "TEMP") == 0) {
      // Version 12 format: read property tag and value
      token = strtok(NULL, " ");
      int itype = 0;
      if (strncmp(token, "PERX", 4) == 0) {
        itype = 1;
      } else if (strncmp(token, "RSVX", 4) == 0) {
        itype = 2;
      } else {
        std::cerr << m_className << "::Initialise:\n";
        std::cerr << "    Found unknown material property flag " << token
                  << "\n";
        std::cerr << "    on material properties file " << mplist << " (line "
                  << il << ").\n";
        ok = false;
      }
      fmplist.getline(line, size, '\n');
      il++;
      token = NULL;
      token = strtok(line, " ");
      if (icurrmat < 1 || icurrmat > m_nMaterials) {
        std::cerr << m_className << "::Initialise:\n";
        std::cerr << "    Found out-of-range current material index "
                  << icurrmat << "\n";
        std::cerr << "    in material properties file " << mplist << ".\n";
        ok = false;
        readerror = false;
      } else if (itype == 1) {
        materials[icurrmat - 1].eps = ReadDouble(token, -1, readerror);
      } else if (itype == 2) {
        materials[icurrmat - 1].ohm = ReadDouble(token, -1, readerror);
      }
      if (readerror) {
        std::cerr << m_className << "::Initialise:\n";
        std::cerr << "    Error reading file " << mplist << " (line " << il
                  << ").\n";
        fmplist.close();
        ok = false;
        return false;
      }
    } else if (strcmp(token, "PROPERTY") == 0) {
      // Version 11 format
      token = strtok(NULL, " ");
      token = strtok(NULL, " ");
      int itype = 0;
      if (strcmp(token, "PERX") == 0) {
        itype = 1;
      } else if (strcmp(token, "RSVX") == 0) {
        itype = 2;
      } else {
        std::cerr << m_className << "::Initialise:\n";
        std::cerr << m_className << "::Initialise:\n";
        std::cerr << "    Found unknown material property flag " << token
                  << "\n";
        std::cerr << "    on material properties file " << mplist << " (line "
                  << il << ").\n";
        ok = false;
      }
      token = strtok(NULL, " ");
      token = strtok(NULL, " ");
      int imat = ReadInteger(token, -1, readerror);
      if (readerror) {
        std::cerr << m_className << "::Initialise:\n";
        std::cerr << "    Error reading file " << mplist << " (line " << il
                  << ").\n";
        fmplist.close();
        ok = false;
        return false;
      } else if (imat < 1 || imat > (int)m_nMaterials) {
        std::cerr << m_className << "::Initialise:\n";
        std::cerr << "    Found out-of-range current material index " << imat
                  << "\n";
        std::cerr << "    in material properties file " << mplist << ".\n";
        ok = false;
      } else {
        fmplist.getline(line, size, '\n');
        il++;
        fmplist.getline(line, size, '\n');
        il++;
        token = NULL;
        token = strtok(line, " ");
        token = strtok(NULL, " ");
        if (itype == 1) {
          materials[imat - 1].eps = ReadDouble(token, -1, readerror);
        } else if (itype == 2) {
          materials[imat - 1].ohm = ReadDouble(token, -1, readerror);
        }
        if (readerror) {
          std::cerr << m_className << "::Initialise:\n";
          std::cerr << "    Error reading file " << mplist << " (line " << il
                    << ").\n";
          fmplist.close();
          ok = false;
          return false;
        }
      }
    }
  }
 
  // Close the file
  fmplist.close();
 
  // Find the lowest epsilon, check for eps = 0, set default drift media
  double epsmin = -1;
  unsigned int iepsmin = 0;
  for (unsigned int imat = 0; imat < m_nMaterials; ++imat) {
    if (materials[imat].eps < 0) continue;
    if (materials[imat].eps == 0) {
      std::cerr << m_className << "::Initialise:\n";
      std::cerr << "    Material " << imat
                << " has been assigned a permittivity\n";
      std::cerr << "    equal to zero in " << mplist << ".\n";
      ok = false;
    } else if (epsmin < 0. || epsmin > materials[imat].eps) {
      epsmin = materials[imat].eps;
      iepsmin = imat;
    }
  }
 
  if (epsmin < 0.) {
    std::cerr << m_className << "::Initialise:\n";
    std::cerr << "    No material with positive permittivity found \n";
    std::cerr << "    in material list " << mplist << ".\n";
    ok = false;
  } else {
    for (unsigned int imat = 0; imat < m_nMaterials; ++imat) {
      if (imat == iepsmin) {
        materials[imat].driftmedium = true;
      } else {
        materials[imat].driftmedium = false;
      }
    }
  }
 
  // Tell how many lines read
  std::cout << m_className << "::Initialise:\n";
  std::cout << "    Read properties of " << m_nMaterials
            << " materials from file " << mplist << ".\n";
  if (m_debug) PrintMaterials();
 
  // Open the element list
  std::ifstream felist;
  felist.open(elist.c_str(), std::ios::in);
  if (felist.fail()) {
    std::cerr << m_className << "::Initialise:\n";
    std::cerr << "    Could not open element file " << elist
              << " for reading.\n";
    std::cerr << "    The file perhaps does not exist.\n";
    return false;
  }
  // Read the element list
  elements.clear();
  nElements = 0;
  Element newElement;
  int ndegenerate = 0;
  int nbackground = 0;
  il = 0;
  int highestnode = 0;
  while (felist.getline(line, size, '\n')) {
    il++;
    // Split the line in tokens
    char* token = NULL;
    // Split into tokens
    token = strtok(line, " ");
    // Skip blank lines and headers
    if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
        int(token[0]) == 10 || int(token[0]) == 13 ||
        strcmp(token, "LIST") == 0 || strcmp(token, "ELEM") == 0)
      continue;
    // Read the element
    int ielem = ReadInteger(token, -1, readerror);
    token = strtok(NULL, " ");
    int imat = ReadInteger(token, -1, readerror);
    token = strtok(NULL, " ");
    token = strtok(NULL, " ");
    token = strtok(NULL, " ");
    token = strtok(NULL, " ");
    token = strtok(NULL, " ");
    if (!token) std::cerr << "No token available\n";
    int in0 = ReadInteger(token, -1, readerror);
    token = strtok(NULL, " ");
    int in1 = ReadInteger(token, -1, readerror);
    token = strtok(NULL, " ");
    int in2 = ReadInteger(token, -1, readerror);
    token = strtok(NULL, " ");
    int in3 = ReadInteger(token, -1, readerror);
    token = strtok(NULL, " ");
    int in4 = ReadInteger(token, -1, readerror);
    token = strtok(NULL, " ");
    int in5 = ReadInteger(token, -1, readerror);
    token = strtok(NULL, " ");
    int in6 = ReadInteger(token, -1, readerror);
    token = strtok(NULL, " ");
    int in7 = ReadInteger(token, -1, readerror);
 
    // Check synchronisation
    if (readerror) {
      std::cerr << m_className << "::Initialise:\n";
      std::cerr << "    Error reading file " << elist << " (line " << il
                << ").\n";
      felist.close();
      ok = false;
      return false;
    } else if (ielem - 1 != nElements + nbackground) {
      std::cerr << m_className << "::Initialise:\n";
      std::cerr << "    Synchronisation lost on file " << elist << " (line "
                << il << ").\n";
      std::cerr << "    Element: " << ielem << " (expected " << nElements
                << "), material: " << imat << ", nodes: (" << in0 << ", " << in1
                << ", " << in2 << ", " << in3 << ", " << in4 << ", " << in5
                << ", " << in6 << ", " << in7 << ")\n";
      ok = false;
    }
    // Check the material number and ensure that epsilon is non-negative
    if (imat < 1 || imat > (int)m_nMaterials) {
      std::cerr << m_className << "::Initialise:\n";
      std::cerr << "   Out-of-range material number on file " << elist
                << " (line " << il << ").\n";
      std::cerr << "    Element: " << ielem << ", material: " << imat
                << ", nodes: (" << in0 << ", " << in1 << ", " << in2 << ", "
                << in3 << ", " << in4 << ", " << in5 << ", " << in6 << ", "
                << in7 << ")\n";
      ok = false;
    }
    if (materials[imat - 1].eps < 0) {
      std::cerr << m_className << "::Initialise:\n";
      std::cerr << "    Element " << ielem << " in element list " << elist
                << " uses material " << imat << " which\n",
          std::cerr << "    has not been assigned a positive permittivity\n";
      std::cerr << "    in material list " << mplist << ".\n";
      ok = false;
    }
    // Check the node numbers
    if (in0 < 1 || in1 < 1 || in2 < 1 || in3 < 1 || in4 < 1 || in5 < 1 ||
        in6 < 1 || in7 < 1) {
      std::cerr << m_className << "::Initialise:\n";
      std::cerr << "    Found a node number < 1 on file " << elist << " (line "
                << il << ").\n";
      std::cerr << "    Element: " << ielem << ", material: " << imat << "\n";
      ok = false;
    }
    if (in0 > highestnode) highestnode = in0;
    if (in1 > highestnode) highestnode = in1;
    if (in2 > highestnode) highestnode = in2;
    if (in3 > highestnode) highestnode = in3;
    if (in4 > highestnode) highestnode = in4;
    if (in5 > highestnode) highestnode = in5;
    if (in6 > highestnode) highestnode = in6;
    if (in7 > highestnode) highestnode = in7;
    // Skip quadrilaterals which are background.
    if (m_deleteBackground && materials[imat - 1].ohm == 0) {
      nbackground++;
      continue;
    }
    // Store the element, degeneracy
    if (in2 == in3 && in3 == in6) {
      ndegenerate++;
      newElement.degenerate = true;
    } else {
      newElement.degenerate = false;
    }
    // Store the material reference
    newElement.matmap = imat - 1;
    // Node references
    if (newElement.degenerate) {
      newElement.emap[0] = in0 - 1;
      newElement.emap[1] = in1 - 1;
      newElement.emap[2] = in2 - 1;
      newElement.emap[3] = in4 - 1;
      newElement.emap[4] = in7 - 1;
      newElement.emap[5] = in5 - 1;
      newElement.emap[6] = in3 - 1;
      newElement.emap[7] = in6 - 1;
    } else {
      newElement.emap[0] = in0 - 1;
      newElement.emap[1] = in1 - 1;
      newElement.emap[2] = in2 - 1;
      newElement.emap[3] = in3 - 1;
      newElement.emap[4] = in4 - 1;
      newElement.emap[5] = in5 - 1;
      newElement.emap[6] = in6 - 1;
      newElement.emap[7] = in7 - 1;
    }
    elements.push_back(newElement);
    ++nElements;
  }
  // Close the file
  felist.close();
  // Tell how many lines read
  std::cout << m_className << "::Initialise:\n";
  std::cout << "    Read " << nElements << " elements from file " << elist
            << ",\n";
  std::cout << "    highest node number: " << highestnode << ",\n";
  std::cout << "    degenerate elements: " << ndegenerate << ",\n";
  std::cout << "    background elements skipped: " << nbackground << ".\n";
  // Check the value of the unit
  double funit;
  if (strcmp(unit.c_str(), "mum") == 0 || strcmp(unit.c_str(), "micron") == 0 ||
      strcmp(unit.c_str(), "micrometer") == 0) {
    funit = 0.0001;
  } else if (strcmp(unit.c_str(), "mm") == 0 ||
             strcmp(unit.c_str(), "millimeter") == 0) {
    funit = 0.1;
  } else if (strcmp(unit.c_str(), "cm") == 0 ||
             strcmp(unit.c_str(), "centimeter") == 0) {
    funit = 1.0;
  } else if (strcmp(unit.c_str(), "m") == 0 ||
             strcmp(unit.c_str(), "meter") == 0) {
    funit = 100.0;
  } else {
    std::cerr << m_className << "::Initialise:\n";
    std::cerr << "    Unknown length unit " << unit << ".\n";
    ok = false;
    funit = 1.0;
  }
  if (m_debug) {
    std::cout << m_className << "::Initialise:\n";
    std::cout << "    Unit scaling factor = " << funit << ".\n";
  }
 
  // Open the node list
  std::ifstream fnlist;
  fnlist.open(nlist.c_str(), std::ios::in);
  if (fnlist.fail()) {
    std::cerr << m_className << "::Initialise:\n";
    std::cerr << "    Could not open nodes file " << nlist << " for reading.\n";
    std::cerr << "    The file perhaps does not exist.\n";
    return false;
  }
  // Read the node list
  nodes.clear();
  nNodes = 0;
  Node newNode;
  newNode.w.clear();
  il = 0;
  while (fnlist.getline(line, size, '\n')) {
    il++;
    // Split the line in tokens
    char* token = NULL;
    // Split into tokens
    token = strtok(line, " ");
    // Skip blank lines and headers
    if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
        int(token[0]) == 10 || int(token[0]) == 13 ||
        strcmp(token, "LIST") == 0 || strcmp(token, "NODE") == 0)
      continue;
    // Read the element
    int inode = ReadInteger(token, -1, readerror);
    token = strtok(NULL, " ");
    double xnode = ReadDouble(token, -1, readerror);
    token = strtok(NULL, " ");
    double ynode = ReadDouble(token, -1, readerror);
    token = strtok(NULL, " ");
    double znode = ReadDouble(token, -1, readerror);
    // Check syntax
    if (readerror) {
      std::cerr << m_className << "::Initialise:\n";
      std::cerr << "    Error reading file " << nlist << " (line " << il
                << ").\n";
      fnlist.close();
      ok = false;
      return false;
    }
    // Check synchronisation
    if (inode - 1 != nNodes) {
      std::cerr << m_className << "::Initialise:\n";
      std::cerr << "    Synchronisation lost on file " << nlist << " (line "
                << il << ").\n";
      std::cerr << "    Node: " << inode << " (expected " << nNodes
                << "), (x,y,z) = (" << xnode << ", " << ynode << ", " << znode
                << ")\n";
      ok = false;
    }
 
    // Store the point coordinates
    newNode.x = xnode * funit;
    newNode.y = ynode * funit;
    newNode.z = znode * funit;
    nodes.push_back(newNode);
    ++nNodes;
  }
  // Close the file
  fnlist.close();
  // Tell how many lines read
  std::cout << m_className << "::Initialise:\n";
  std::cout << "    Read " << nNodes << " nodes from file " << nlist << ".\n";
  // Check number of nodes
  if (nNodes != highestnode) {
    std::cerr << m_className << "::Initialise:\n";
    std::cerr << "    Number of nodes read (" << nNodes << ") on " << nlist
              << "\n";
    std::cerr << "    does not match element list (" << highestnode << ").\n";
    ok = false;
  }
 
  // Open the voltage list
  std::ifstream fprnsol;
  fprnsol.open(prnsol.c_str(), std::ios::in);
  if (fprnsol.fail()) {
    std::cerr << m_className << "::Initialise:\n";
    std::cerr << "    Could not open potential file " << prnsol
              << " for reading.\n";
    std::cerr << "    The file perhaps does not exist.\n";
    return false;
  }
  // Read the voltage list
  il = 0;
  int nread = 0;
  readerror = false;
  while (fprnsol.getline(line, size, '\n')) {
    il++;
    // Split the line in tokens
    char* token = NULL;
    token = strtok(line, " ");
    // Skip blank lines and headers
    if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
        int(token[0]) == 10 || int(token[0]) == 13 ||
        strcmp(token, "PRINT") == 0 || strcmp(token, "*****") == 0 ||
        strcmp(token, "LOAD") == 0 || strcmp(token, "TIME=") == 0 ||
        strcmp(token, "MAXIMUM") == 0 || strcmp(token, "VALUE") == 0 ||
        strcmp(token, "NODE") == 0)
      continue;
    // Read the element
    int inode = ReadInteger(token, -1, readerror);
    token = strtok(NULL, " ");
    double volt = ReadDouble(token, -1, readerror);
    // Check syntax
    if (readerror) {
      std::cerr << m_className << "::Initialise:\n";
      std::cerr << "    Error reading file " << prnsol << " (line " << il
                << ").\n";
      fprnsol.close();
      ok = false;
      return false;
    }
    // Check node number and store if OK
    if (inode < 1 || inode > nNodes) {
      std::cerr << m_className << "::Initialise:\n";
      std::cerr << "    Node number " << inode << " out of range\n";
      std::cerr << "    on potential file " << prnsol << " (line " << il
                << ").\n";
      ok = false;
    } else {
      nodes[inode - 1].v = volt;
      nread++;
    }
  }
  // Close the file
  fprnsol.close();
  // Tell how many lines read
  std::cout << m_className << "::Initialise:\n";
  std::cout << "    Read " << nread << " potentials from file " << prnsol
            << ".\n";
  // Check number of nodes
  if (nread != nNodes) {
    std::cerr << m_className << "::Initialise:\n";
    std::cerr << "    Number of nodes read (" << nread << ") on potential file "
              << prnsol << " does not\n",
        std::cerr << "    match the node list (" << nNodes << ").\n";
    ok = false;
  }
  // Set the ready flag
  if (ok) {
    m_ready = true;
  } else {
    std::cerr << m_className << "::Initialise:\n";
    std::cerr
        << "    Field map could not be read and cannot be interpolated.\n";
    return false;
  }
 
  // Remove weighting fields (if any).
  wfields.clear();
  wfieldsOk.clear();
  nWeightingFields = 0;
 
  // Establish the ranges
  SetRange();
  UpdatePeriodicity();
  return true;
}

SetRangeZ()

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

Definition at line 1090 of file ComponentAnsys121.cc.

                                                                      {
  if (fabs(zmax - zmin) <= 0.) {
    std::cerr << m_className << "::SetRangeZ: Zero range is not permitted.\n";
    return;
  }
  m_minBoundingBox[2] = m_mapmin[2] = std::min(zmin, zmax);
  m_maxBoundingBox[2] = m_mapmax[2] = std::max(zmin, zmax);
}

SetWeightingField()

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

Definition at line 592 of file ComponentAnsys121.cc.

                                                           {
  if (!m_ready) {
    PrintNotReady("SetWeightingField");
    std::cerr << "    Weighting field cannot be added.\n";
    return false;
  }
 
  // Open the voltage list.
  std::ifstream fprnsol;
  fprnsol.open(prnsol.c_str(), std::ios::in);
  if (fprnsol.fail()) {
    std::cerr << m_className << "::SetWeightingField:\n";
    std::cerr << "    Could not open potential file " << prnsol
              << " for reading.\n";
    std::cerr << "    The file perhaps does not exist.\n";
    return false;
  }
 
  // Check if a weighting field with the same label already exists.
  int iw = nWeightingFields;
  for (int i = nWeightingFields; i--;) {
    if (wfields[i] == label) {
      iw = i;
      break;
    }
  }
  if (iw == nWeightingFields) {
    ++nWeightingFields;
    wfields.resize(nWeightingFields);
    wfieldsOk.resize(nWeightingFields);
    for (int j = nNodes; j--;) {
      nodes[j].w.resize(nWeightingFields);
    }
  } else {
    std::cout << m_className << "::SetWeightingField:\n";
    std::cout << "    Replacing existing weighting field " << label << ".\n";
  }
  wfields[iw] = label;
  wfieldsOk[iw] = false;
 
  // Buffer for reading
  const int size = 100;
  char line[size];
 
  bool ok = true;
  // Read the voltage list.
  int il = 0;
  int nread = 0;
  bool readerror = false;
  while (fprnsol.getline(line, size, '\n')) {
    il++;
    // Split the line in tokens.
    char* token = NULL;
    token = strtok(line, " ");
    // Skip blank lines and headers.
    if (!token || strcmp(token, " ") == 0 || strcmp(token, "\n") == 0 ||
        int(token[0]) == 10 || int(token[0]) == 13 ||
        strcmp(token, "PRINT") == 0 || strcmp(token, "*****") == 0 ||
        strcmp(token, "LOAD") == 0 || strcmp(token, "TIME=") == 0 ||
        strcmp(token, "MAXIMUM") == 0 || strcmp(token, "VALUE") == 0 ||
        strcmp(token, "NODE") == 0)
      continue;
    // Read the element.
    int inode = ReadInteger(token, -1, readerror);
    token = strtok(NULL, " ");
    double volt = ReadDouble(token, -1, readerror);
    // Check the syntax.
    if (readerror) {
      std::cerr << m_className << "::SetWeightingField:\n";
      std::cerr << "    Error reading file " << prnsol << " (line " << il
                << ").\n";
      fprnsol.close();
      return false;
    }
    // Check node number and store if OK.
    if (inode < 1 || inode > nNodes) {
      std::cerr << m_className << "::SetWeightingField:\n";
      std::cerr << "    Node number " << inode << " out of range\n";
      std::cerr << "    on potential file " << prnsol << " (line " << il
                << ").\n";
      ok = false;
    } else {
      nodes[inode - 1].w[iw] = volt;
      nread++;
    }
  }
  // Close the file.
  fprnsol.close();
  // Tell how many lines read.
  std::cout << m_className << "::SetWeightingField:\n";
  std::cout << "    Read " << nread << " potentials from file " << prnsol
            << ".\n";
  // Check the number of nodes.
  if (nread != nNodes) {
    std::cerr << m_className << "::SetWeightingField:\n";
    std::cerr << "    Number of nodes read (" << nread << ") "
              << "    on potential file " << prnsol << "\n";
    std::cerr << "    does not match the node list (" << nNodes << ").\n";
    ok = false;
  }
 
  // Set the ready flag.
  wfieldsOk[iw] = ok;
  if (!ok) {
    std::cerr << m_className << "::SetWeightingField:\n";
    std::cerr << "    Field map could not be read "
              << "and cannot be interpolated.\n";
    return false;
  }
  return true;
}

UpdatePeriodicity()

void Garfield::ComponentAnsys121::UpdatePeriodicity ( )

overrideprotectedvirtual

Verify periodicities.

Implements Garfield::ComponentBase.

Definition at line 1099 of file ComponentAnsys121.cc.

                                          {
  UpdatePeriodicity2d();
  UpdatePeriodicityCommon();
}

Referenced by Initialise().

WeightingField()

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

overridevirtual

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

Parameters

x,y,z	coordinates [cm].
wx,wy,wz	components of the weighting field [1/cm].
label	name of the electrode

Reimplemented from Garfield::ComponentBase.

Definition at line 855 of file ComponentAnsys121.cc.

                                                                           {
  // Initial values
  wx = wy = wz = 0;
 
  // Do not proceed if not properly initialised.
  if (!m_ready) return;
 
  // Look for the label.
  int iw = 0;
  bool found = false;
  for (int i = nWeightingFields; i--;) {
    if (wfields[i] == label) {
      iw = i;
      found = true;
      break;
    }
  }
 
  // Do not proceed if the requested weighting field does not exist.
  if (!found) return;
  // Check if the weighting field is properly initialised.
  if (!wfieldsOk[iw]) return;
 
  // Copy the coordinates.
  double x = xin, y = yin, z = zin;
 
  // Map the coordinates onto field map coordinates
  bool xmirr, ymirr, zmirr;
  double rcoordinate, rotation;
  MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
 
  if (m_warning) PrintWarning("WeightingField");
 
  // Find the element that contains this point.
  double t1, t2, t3, t4, jac[4][4], det;
  const int imap = FindElement5(x, y, z, t1, t2, t3, t4, jac, det);
  // Check if the point is in the mesh.
  if (imap < 0) return;
 
  const Element& element = elements[imap];
  if (m_debug) {
    PrintElement("WeightingField", x, y, z, t1, t2, t3, t4, element, 8, iw);
  }
  const Node& n0 = nodes[element.emap[0]];
  const Node& n1 = nodes[element.emap[1]];
  const Node& n2 = nodes[element.emap[2]];
  const Node& n3 = nodes[element.emap[3]];
  const Node& n4 = nodes[element.emap[4]];
  const Node& n5 = nodes[element.emap[5]];
  // Calculate quadrilateral field, which can degenerate to a triangular field
  const double invdet = 1. / det;
  if (elements[imap].degenerate) {
    wx = -(n0.w[iw] * (4 * t1 - 1) * jac[0][1] +
           n1.w[iw] * (4 * t2 - 1) * jac[1][1] +
           n2.w[iw] * (4 * t3 - 1) * jac[2][1] +
           n3.w[iw] * (4 * t2 * jac[0][1] + 4 * t1 * jac[1][1]) +
           n4.w[iw] * (4 * t3 * jac[0][1] + 4 * t1 * jac[2][1]) +
           n5.w[iw] * (4 * t3 * jac[1][1] + 4 * t2 * jac[2][1])) *
         invdet;
    wy = -(n0.w[iw] * (4 * t1 - 1) * jac[0][2] +
           n1.w[iw] * (4 * t2 - 1) * jac[1][2] +
           n2.w[iw] * (4 * t3 - 1) * jac[2][2] +
           n3.w[iw] * (4 * t2 * jac[0][2] + 4 * t1 * jac[1][2]) +
           n4.w[iw] * (4 * t3 * jac[0][2] + 4 * t1 * jac[2][2]) +
           n5.w[iw] * (4 * t3 * jac[1][2] + 4 * t2 * jac[2][2])) *
         invdet;
  } else {
    const Node& n6 = nodes[element.emap[6]];
    const Node& n7 = nodes[element.emap[7]];
    wx = -(n0.w[iw] * ((1 - t2) * (2 * t1 + t2) * jac[0][0] +
                       (1 - t1) * (t1 + 2 * t2) * jac[1][0]) *
               0.25 +
           n1.w[iw] * ((1 - t2) * (2 * t1 - t2) * jac[0][0] -
                       (1 + t1) * (t1 - 2 * t2) * jac[1][0]) *
               0.25 +
           n2.w[iw] * ((1 + t2) * (2 * t1 + t2) * jac[0][0] +
                       (1 + t1) * (t1 + 2 * t2) * jac[1][0]) *
               0.25 +
           n3.w[iw] * ((1 + t2) * (2 * t1 - t2) * jac[0][0] -
                       (1 - t1) * (t1 - 2 * t2) * jac[1][0]) *
               0.25 +
           n4.w[iw] * (t1 * (t2 - 1) * jac[0][0] +
                       (t1 - 1) * (t1 + 1) * jac[1][0] * 0.5) +
           n5.w[iw] * ((1 - t2) * (1 + t2) * jac[0][0] * 0.5 -
                       (1 + t1) * t2 * jac[1][0]) +
           n6.w[iw] * (-t1 * (1 + t2) * jac[0][0] +
                       (1 - t1) * (1 + t1) * jac[1][0] * 0.5) +
           n7.w[iw] * ((t2 - 1) * (1 + t2) * jac[0][0] * 0.5 +
                       (t1 - 1) * t2 * jac[1][0])) *
         invdet;
    wy = -(n0.w[iw] * ((1 - t2) * (2 * t1 + t2) * jac[0][1] +
                       (1 - t1) * (t1 + 2 * t2) * jac[1][1]) *
               0.25 +
           n1.w[iw] * ((1 - t2) * (2 * t1 - t2) * jac[0][1] -
                       (1 + t1) * (t1 - 2 * t2) * jac[1][1]) *
               0.25 +
           n2.w[iw] * ((1 + t2) * (2 * t1 + t2) * jac[0][1] +
                       (1 + t1) * (t1 + 2 * t2) * jac[1][1]) *
               0.25 +
           n3.w[iw] * ((1 + t2) * (2 * t1 - t2) * jac[0][1] -
                       (1 - t1) * (t1 - 2 * t2) * jac[1][1]) *
               0.25 +
           n4.w[iw] * (t1 * (t2 - 1) * jac[0][1] +
                       (t1 - 1) * (t1 + 1) * jac[1][1] * 0.5) +
           n5.w[iw] * ((1 - t2) * (1 + t2) * jac[0][1] * 0.5 -
                       (1 + t1) * t2 * jac[1][1]) +
           n6.w[iw] * (-t1 * (1 + t2) * jac[0][1] +
                       (1 - t1) * (1 + t1) * jac[1][1] * 0.5) +
           n7.w[iw] * ((t2 - 1) * (t2 + 1) * jac[0][1] * 0.5 +
                       (t1 - 1) * t2 * jac[1][1])) *
         invdet;
  }
 
  // Transform field to global coordinates
  UnmapFields(wx, wy, wz, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
}

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,z	coordinates [cm].
label	name of the electrode.

Returns

weighting potential [dimensionless].

Reimplemented from Garfield::ComponentBase.

Definition at line 974 of file ComponentAnsys121.cc.

                                                                     {
  // Do not proceed if not properly initialised.
  if (!m_ready) return 0.;
 
  // Look for the label.
  int iw = 0;
  bool found = false;
  for (int i = nWeightingFields; i--;) {
    if (wfields[i] == label) {
      iw = i;
      found = true;
      break;
    }
  }
 
  // Do not proceed if the requested weighting field does not exist.
  if (!found) return 0.;
  // Check if the weighting field is properly initialised.
  if (!wfieldsOk[iw]) return 0.;
 
  // Copy the coordinates.
  double x = xin, y = yin, z = zin;
 
  // Map the coordinates onto field map coordinates.
  bool xmirr, ymirr, zmirr;
  double rcoordinate, rotation;
  MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
 
  if (m_warning) PrintWarning("WeightingPotential");
 
  // Find the element that contains this point.
  double t1, t2, t3, t4, jac[4][4], det;
  const int imap = FindElement5(x, y, z, t1, t2, t3, t4, jac, det);
  // Check if the point is in the mesh
  if (imap < 0) return 0.;
 
  const Element& element = elements[imap];
  if (m_debug) {
    PrintElement("WeightingPotential", x, y, z, t1, t2, t3, t4, element, 8, iw);
  }
  // Calculate quadrilateral field, which can degenerate to a triangular field
  const Node& n0 = nodes[element.emap[0]];
  const Node& n1 = nodes[element.emap[1]];
  const Node& n2 = nodes[element.emap[2]];
  const Node& n3 = nodes[element.emap[3]];
  const Node& n4 = nodes[element.emap[4]];
  const Node& n5 = nodes[element.emap[5]];
  if (element.degenerate) {
    return n0.w[iw] * t1 * (2 * t1 - 1) + n1.w[iw] * t2 * (2 * t2 - 1) +
           n2.w[iw] * t3 * (2 * t3 - 1) + 4 * n3.w[iw] * t1 * t2 +
           4 * n4.w[iw] * t1 * t3 + 4 * n5.w[iw] * t2 * t3;
  }
 
  const Node& n6 = nodes[element.emap[6]];
  const Node& n7 = nodes[element.emap[7]];
  return -n0.w[iw] * (1 - t1) * (1 - t2) * (1 + t1 + t2) * 0.25 -
         n1.w[iw] * (1 + t1) * (1 - t2) * (1 - t1 + t2) * 0.25 -
         n2.w[iw] * (1 + t1) * (1 + t2) * (1 - t1 - t2) * 0.25 -
         n3.w[iw] * (1 - t1) * (1 + t2) * (1 + t1 - t2) * 0.25 +
         n4.w[iw] * (1 - t1) * (1 + t1) * (1 - t2) * 0.5 +
         n5.w[iw] * (1 + t1) * (1 + t2) * (1 - t2) * 0.5 +
         n6.w[iw] * (1 - t1) * (1 + t1) * (1 + t2) * 0.5 +
         n7.w[iw] * (1 - t1) * (1 + t2) * (1 - t2) * 0.5;
}

---

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

• ComponentAnsys121.hh
• ComponentAnsys121.cc