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

Component for importing field maps computed by Elmer. More...

#include <ComponentElmer.hh>

+ Inheritance diagram for Garfield::ComponentElmer:

Public Member Functions

 ComponentElmer ()
 Default constructor.
 
 ComponentElmer (const std::string &header, const std::string &elist, const std::string &nlist, const std::string &mplist, const std::string &volt, const std::string &unit)
 Constructor with a set of field map files, see Initialise().
 
 ~ComponentElmer ()
 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 (const std::string &header="mesh.header", const std::string &elist="mesh.elements", const std::string &nlist="mesh.nodes", const std::string &mplist="dielectrics.dat", const std::string &volt="out.result", const std::string &unit="cm")
 
bool SetWeightingField (std::string prnsol, std::string label)
 Import a list of voltages to be used as weighting field.
 
- 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 field maps computed by Elmer.

Definition at line 12 of file ComponentElmer.hh.

Constructor & Destructor Documentation

◆ ComponentElmer() [1/2]

Garfield::ComponentElmer::ComponentElmer ( )

Default constructor.

Definition at line 20 of file ComponentElmer.cc.

20: ComponentFieldMap("Elmer") {}
Garfield::ComponentFieldMap::ComponentFieldMap
ComponentFieldMap()=delete
Default constructor.

◆ ComponentElmer() [2/2]

Garfield::ComponentElmer::ComponentElmer ( const std::string &  header,
const std::string &  elist,
const std::string &  nlist,
const std::string &  mplist,
const std::string &  volt,
const std::string &  unit 
)

Constructor with a set of field map files, see Initialise().

Definition at line 22 of file ComponentElmer.cc.

27 : ComponentFieldMap("Elmer") {
28 Initialise(header, elist, nlist, mplist, volt, unit);
29}
Garfield::ComponentElmer::Initialise
bool Initialise(const std::string &header="mesh.header", const std::string &elist="mesh.elements", const std::string &nlist="mesh.nodes", const std::string &mplist="dielectrics.dat", const std::string &volt="out.result", const std::string &unit="cm")
Definition: ComponentElmer.cc:31

◆ ~ComponentElmer()

Garfield::ComponentElmer::~ComponentElmer ( )
inline

Destructor.

Definition at line 21 of file ComponentElmer.hh.

21{}

Member Function Documentation

◆ ElectricField() [1/2]

void Garfield::ComponentElmer::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 464 of file ComponentElmer.cc.

467 {
468 // Copy the coordinates
469 double x = xin, y = yin, z = zin;
470
471 // Map the coordinates onto field map coordinates
472 bool xmirr, ymirr, zmirr;
473 double rcoordinate, rotation;
474 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
475
476 // Initial values
477 ex = ey = ez = volt = 0.;
478 status = 0;
479 m = nullptr;
480
481 // Do not proceed if not properly initialised.
482 if (!m_ready) {
483 status = -10;
484 PrintNotReady("ElectricField");
485 return;
486 }
487
488 if (m_warning) PrintWarning("ElectricField");
489
490 // Find the element that contains this point
491 double t1, t2, t3, t4, jac[4][4], det;
492 const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
493 if (imap < 0) {
494 if (m_debug) {
495 std::cout << m_className << "::ElectricField:\n Point (" << x << ", "
496 << y << ", " << z << ") is not in the mesh.\n";
497 }
498 status = -6;
499 return;
500 }
501
502 const Element& element = m_elements[imap];
503 if (m_debug) {
504 PrintElement("ElectricField", x, y, z, t1, t2, t3, t4, element, 10);
505 }
506 const Node& n0 = m_nodes[element.emap[0]];
507 const Node& n1 = m_nodes[element.emap[1]];
508 const Node& n2 = m_nodes[element.emap[2]];
509 const Node& n3 = m_nodes[element.emap[3]];
510 const Node& n4 = m_nodes[element.emap[4]];
511 const Node& n5 = m_nodes[element.emap[5]];
512 const Node& n6 = m_nodes[element.emap[6]];
513 const Node& n7 = m_nodes[element.emap[7]];
514 const Node& n8 = m_nodes[element.emap[8]];
515 const Node& n9 = m_nodes[element.emap[9]];
516 // Shorthands.
517 const double fourt1 = 4 * t1;
518 const double fourt2 = 4 * t2;
519 const double fourt3 = 4 * t3;
520 const double fourt4 = 4 * t4;
521 const double invdet = 1. / det;
522 // Tetrahedral field
523 volt = n0.v * t1 * (2 * t1 - 1) + n1.v * t2 * (2 * t2 - 1) +
524 n2.v * t3 * (2 * t3 - 1) + n3.v * t4 * (2 * t4 - 1) +
525 4 * n4.v * t1 * t2 + 4 * n5.v * t1 * t3 + 4 * n6.v * t1 * t4 +
526 4 * n7.v * t2 * t3 + 4 * n8.v * t2 * t4 + 4 * n9.v * t3 * t4;
527 ex = -(n0.v * (fourt1 - 1) * jac[0][1] + n1.v * (fourt2 - 1) * jac[1][1] +
528 n2.v * (fourt3 - 1) * jac[2][1] + n3.v * (fourt4 - 1) * jac[3][1] +
529 n4.v * (fourt2 * jac[0][1] + fourt1 * jac[1][1]) +
530 n5.v * (fourt3 * jac[0][1] + fourt1 * jac[2][1]) +
531 n6.v * (fourt4 * jac[0][1] + fourt1 * jac[3][1]) +
532 n7.v * (fourt3 * jac[1][1] + fourt2 * jac[2][1]) +
533 n8.v * (fourt4 * jac[1][1] + fourt2 * jac[3][1]) +
534 n9.v * (fourt4 * jac[2][1] + fourt3 * jac[3][1])) *
535 invdet;
536 ey = -(n0.v * (fourt1 - 1) * jac[0][2] + n1.v * (fourt2 - 1) * jac[1][2] +
537 n2.v * (fourt3 - 1) * jac[2][2] + n3.v * (fourt4 - 1) * jac[3][2] +
538 n4.v * (fourt2 * jac[0][2] + fourt1 * jac[1][2]) +
539 n5.v * (fourt3 * jac[0][2] + fourt1 * jac[2][2]) +
540 n6.v * (fourt4 * jac[0][2] + fourt1 * jac[3][2]) +
541 n7.v * (fourt3 * jac[1][2] + fourt2 * jac[2][2]) +
542 n8.v * (fourt4 * jac[1][2] + fourt2 * jac[3][2]) +
543 n9.v * (fourt4 * jac[2][2] + fourt3 * jac[3][2])) *
544 invdet;
545 ez = -(n0.v * (fourt1 - 1) * jac[0][3] + n1.v * (fourt2 - 1) * jac[1][3] +
546 n2.v * (fourt3 - 1) * jac[2][3] + n3.v * (fourt4 - 1) * jac[3][3] +
547 n4.v * (fourt2 * jac[0][3] + fourt1 * jac[1][3]) +
548 n5.v * (fourt3 * jac[0][3] + fourt1 * jac[2][3]) +
549 n6.v * (fourt4 * jac[0][3] + fourt1 * jac[3][3]) +
550 n7.v * (fourt3 * jac[1][3] + fourt2 * jac[2][3]) +
551 n8.v * (fourt4 * jac[1][3] + fourt2 * jac[3][3]) +
552 n9.v * (fourt4 * jac[2][3] + fourt3 * jac[3][3])) *
553 invdet;
554
555 // Transform field to global coordinates
556 UnmapFields(ex, ey, ez, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
557
558 // Drift medium?
559 const Material& mat = m_materials[element.matmap];
560 if (m_debug) {
561 std::cout << m_className << "::ElectricField:\n Material "
562 << element.matmap << ", drift flag " << mat.driftmedium << ".\n";
563 }
564 m = mat.medium;
565 status = -5;
566 if (mat.driftmedium && m && m->IsDriftable()) status = 0;
567}
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::ComponentElmer::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 457 of file ComponentElmer.cc.

459 {
460 double v = 0.;
461 ElectricField(x, y, z, ex, ey, ez, v, m, status);
462}
Garfield::ComponentElmer::ElectricField
void ElectricField(const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status) override
Definition: ComponentElmer.cc:457

Referenced by ElectricField().

◆ GetAspectRatio()

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

Implements Garfield::ComponentFieldMap.

Definition at line 774 of file ComponentElmer.cc.

775 {
776 if (i >= m_elements.size()) {
777 dmin = dmax = 0.;
778 return;
779 }
780
781 const Element& element = m_elements[i];
782 const int np = 4;
783 // Loop over all pairs of vertices.
784 for (int j = 0; j < np - 1; ++j) {
785 const Node& nj = m_nodes[element.emap[j]];
786 for (int k = j + 1; k < np; ++k) {
787 const Node& nk = m_nodes[element.emap[k]];
788 // Compute distance.
789 const double dx = nj.x - nk.x;
790 const double dy = nj.y - nk.y;
791 const double dz = nj.z - nk.z;
792 const double dist = sqrt(dx * dx + dy * dy + dz * dz);
793 if (k == 1) {
794 dmin = dmax = dist;
795 } else {
796 if (dist < dmin) dmin = dist;
797 if (dist > dmax) dmax = dist;
798 }
799 }
800 }
801}
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314

◆ GetElementVolume()

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

Implements Garfield::ComponentFieldMap.

Definition at line 753 of file ComponentElmer.cc.

753 {
754 if (i >= m_elements.size()) return 0.;
755 const Element& element = m_elements[i];
756 const Node& n0 = m_nodes[element.emap[0]];
757 const Node& n1 = m_nodes[element.emap[1]];
758 const Node& n2 = m_nodes[element.emap[2]];
759 const Node& n3 = m_nodes[element.emap[3]];
760
761 // Uses formula V = |a (dot) b x c|/6
762 // with a => "3", b => "1", c => "2" and origin "0"
763 const double vol =
764 fabs((n3.x - n0.x) *
765 ((n1.y - n0.y) * (n2.z - n0.z) - (n2.y - n0.y) * (n1.z - n0.z)) +
766 (n3.y - n0.y) *
767 ((n1.z - n0.z) * (n2.x - n0.x) - (n2.z - n0.z) * (n1.x - n0.x)) +
768 (n3.z - n0.z) * ((n1.x - n0.x) * (n2.y - n0.y) -
769 (n3.x - n0.x) * (n1.y - n0.y))) /
770 6.;
771 return vol;
772}
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615

◆ GetMedium()

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

712 {
713 // Copy the coordinates
714 double x = xin, y = yin, z = zin;
715
716 // Map the coordinates onto field map coordinates
717 bool xmirr, ymirr, zmirr;
718 double rcoordinate, rotation;
719 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
720
721 // Do not proceed if not properly initialised.
722 if (!m_ready) {
723 PrintNotReady("GetMedium");
724 return nullptr;
725 }
726 if (m_warning) PrintWarning("GetMedium");
727
728 // Find the element that contains this point
729 double t1, t2, t3, t4, jac[4][4], det;
730 const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
731 if (imap < 0) {
732 if (m_debug) {
733 std::cout << m_className << "::GetMedium:\n Point (" << x << ", " << y
734 << ", " << z << ") is not in the mesh.\n";
735 }
736 return nullptr;
737 }
738 const Element& element = m_elements[imap];
739 if (element.matmap >= m_materials.size()) {
740 if (m_debug) {
741 std::cerr << m_className << "::GetMedium:\n Point (" << x << ", " << y
742 << ", " << z << ") has out of range material number " << imap
743 << ".\n";
744 }
745 return nullptr;
746 }
747
748 if (m_debug) PrintElement("GetMedium", x, y, z, t1, t2, t3, t4, element, 10);
749
750 return m_materials[element.matmap].medium;
751}

◆ Initialise()

bool Garfield::ComponentElmer::Initialise ( const std::string &  header = "mesh.header",
const std::string &  elist = "mesh.elements",
const std::string &  nlist = "mesh.nodes",
const std::string &  mplist = "dielectrics.dat",
const std::string &  volt = "out.result",
const std::string &  unit = "cm" 
)

Import a field map from a set of files.

Parameters
headername of the header file (contains the number of elements and nodes).
elistname of the file that contains the list of mesh elements
nlistname of the file that contains the list of mesh nodes
mplistname of the file that contains the material properties
voltoutput of the field solver (list of voltages)
unitlength unit to be used

Definition at line 31 of file ComponentElmer.cc.

36 {
37 const std::string hdr = m_className + "::Initialise:";
38 Reset();
39
40 // Keep track of the success.
41 bool ok = true;
42
43 // Buffer for reading
44 constexpr int size = 100;
45 char line[size];
46
47 // Open the header.
48 std::ifstream fheader;
49 fheader.open(header.c_str(), std::ios::in);
50 if (fheader.fail()) {
51 PrintCouldNotOpen("Initialise", header);
52 return false;
53 }
54
55 // Temporary variables for use in file reading
56 char* token = NULL;
57 bool readerror = false;
58 bool readstop = false;
59 int il = 0;
60
61 // Read the header to get the number of nodes and elements.
62 fheader.getline(line, size, '\n');
63 token = strtok(line, " ");
64 const int nNodes = ReadInteger(token, 0, readerror);
65 token = strtok(NULL, " ");
66 const int nElements = ReadInteger(token, 0, readerror);
67 std::cout << hdr << "\n Read " << nNodes << " nodes and " << nElements
68 << " elements from file " << header << ".\n";
69 if (readerror) {
70 PrintErrorReadingFile(hdr, header, il);
71 fheader.close();
72 return false;
73 }
74
75 // Close the header file.
76 fheader.close();
77
78 // Open the nodes list.
79 std::ifstream fnodes;
80 fnodes.open(nlist.c_str(), std::ios::in);
81 if (fnodes.fail()) {
82 PrintCouldNotOpen("Initialise", nlist);
83 return false;
84 }
85
86 // Check the value of the unit.
87 double funit = ScalingFactor(unit);
88 if (funit <= 0.) {
89 std::cerr << hdr << " Unknown length unit " << unit << ".\n";
90 ok = false;
91 funit = 1.0;
92 }
93 if (m_debug) std::cout << hdr << " Unit scaling factor = " << funit << ".\n";
94
95 // Read the nodes from the file.
96 for (il = 0; il < nNodes; il++) {
97 // Get a line from the nodes file.
98 fnodes.getline(line, size, '\n');
99
100 // Ignore the first two characters.
101 token = strtok(line, " ");
102 token = strtok(NULL, " ");
103
104 // Get the node coordinates.
105 token = strtok(NULL, " ");
106 double xnode = ReadDouble(token, -1, readerror);
107 token = strtok(NULL, " ");
108 double ynode = ReadDouble(token, -1, readerror);
109 token = strtok(NULL, " ");
110 double znode = ReadDouble(token, -1, readerror);
111 if (readerror) {
112 PrintErrorReadingFile(hdr, nlist, il);
113 fnodes.close();
114 return false;
115 }
116
117 // Set up and create a new node.
118 Node newNode;
119 newNode.w.clear();
120 newNode.x = xnode * funit;
121 newNode.y = ynode * funit;
122 newNode.z = znode * funit;
123 m_nodes.push_back(std::move(newNode));
124 }
125
126 // Close the nodes file.
127 fnodes.close();
128
129 // Open the potential file.
130 std::ifstream fvolt;
131 fvolt.open(volt.c_str(), std::ios::in);
132 if (fvolt.fail()) {
133 PrintCouldNotOpen("Initialise", volt);
134 return false;
135 }
136
137 // Reset the line counter.
138 il = 1;
139
140 // Read past the header.
141 while (!readstop && fvolt.getline(line, size, '\n')) {
142 token = strtok(line, " ");
143 if (strcmp(token, "Perm:") == 0) readstop = true;
144 il++;
145 }
146
147 // Should have stopped: if not, print error message.
148 if (!readstop) {
149 std::cerr << hdr << "\n Error reading past header of potentials file "
150 << volt << ".\n";
151 fvolt.close();
152 return false;
153 }
154
155 // Read past the permutation map (number of lines = nNodes).
156 for (int tl = 0; tl < nNodes; tl++) {
157 fvolt.getline(line, size, '\n');
158 il++;
159 }
160
161 // Read the potentials.
162 for (int tl = 0; tl < nNodes; tl++) {
163 fvolt.getline(line, size, '\n');
164 token = strtok(line, " ");
165 double v = ReadDouble(token, -1, readerror);
166 if (readerror) {
167 PrintErrorReadingFile(hdr, volt, il);
168 fvolt.close();
169 return false;
170 }
171 // Place the voltage in its appropriate node.
172 m_nodes[tl].v = v;
173 }
174
175 // Close the potentials file.
176 fvolt.close();
177
178 // Open the materials file.
179 std::ifstream fmplist;
180 fmplist.open(mplist.c_str(), std::ios::in);
181 if (fmplist.fail()) {
182 PrintCouldNotOpen("Initialise", mplist);
183 return false;
184 }
185
186 // Read the dielectric constants from the materials file.
187 fmplist.getline(line, size, '\n');
188 token = strtok(line, " ");
189 if (readerror) {
190 std::cerr << hdr << "\n Error reading number of materials from "
191 << mplist << ".\n";
192 fmplist.close();
193 return false;
194 }
195 const unsigned int nMaterials = ReadInteger(token, 0, readerror);
196 m_materials.resize(nMaterials);
197 for (auto& material : m_materials) {
198 material.ohm = -1;
199 material.eps = -1;
200 material.medium = nullptr;
201 }
202 for (il = 2; il < ((int)nMaterials + 2); il++) {
203 fmplist.getline(line, size, '\n');
204 token = strtok(line, " ");
205 ReadInteger(token, -1, readerror);
206 token = strtok(NULL, " ");
207 double dc = ReadDouble(token, -1.0, readerror);
208 if (readerror) {
209 PrintErrorReadingFile(hdr, mplist, il);
210 fmplist.close();
211 return false;
212 }
213 m_materials[il - 2].eps = dc;
214 std::cout << hdr << "\n Set material " << il - 2 << " of "
215 << nMaterials << " to eps " << dc << ".\n";
216 }
217
218 // Close the materials file.
219 fmplist.close();
220
221 // Find lowest epsilon, check for eps = 0, set default drift medium.
222 if (!SetDefaultDriftMedium()) ok = false;
223
224 // Open the elements file.
225 std::ifstream felems;
226 felems.open(elist.c_str(), std::ios::in);
227 if (felems.fail()) {
228 PrintCouldNotOpen("Initialise", elist);
229 return false;
230 }
231
232 // Read the elements and their material indices.
233 m_elements.clear();
234 int highestnode = 0;
235
236 for (il = 0; il < nElements; il++) {
237 // Get a line
238 felems.getline(line, size, '\n');
239
240 // Split into tokens.
241 token = strtok(line, " ");
242 // Read the 2nd-order element
243 // Note: Ordering of Elmer elements can be described in the
244 // ElmerSolver manual (appendix D. at the time of this comment)
245 // If the order read below is compared to the shape functions used
246 // eg. in ElectricField, the order is wrong, but note at the
247 // end of this function the order of elements 5,6,7 will change to
248 // 7,5,6 when actually recorded in newElement.emap to correct for this
249 token = strtok(NULL, " ");
250 int imat = ReadInteger(token, -1, readerror) - 1;
251 token = strtok(NULL, " ");
252 token = strtok(NULL, " ");
253 int in0 = ReadInteger(token, -1, readerror);
254 token = strtok(NULL, " ");
255 int in1 = ReadInteger(token, -1, readerror);
256 token = strtok(NULL, " ");
257 int in2 = ReadInteger(token, -1, readerror);
258 token = strtok(NULL, " ");
259 int in3 = ReadInteger(token, -1, readerror);
260 token = strtok(NULL, " ");
261 int in4 = ReadInteger(token, -1, readerror);
262 token = strtok(NULL, " ");
263 int in5 = ReadInteger(token, -1, readerror);
264 token = strtok(NULL, " ");
265 int in6 = ReadInteger(token, -1, readerror);
266 token = strtok(NULL, " ");
267 int in7 = ReadInteger(token, -1, readerror);
268 token = strtok(NULL, " ");
269 int in8 = ReadInteger(token, -1, readerror);
270 token = strtok(NULL, " ");
271 int in9 = ReadInteger(token, -1, readerror);
272
273 if (m_debug && il < 10) {
274 std::cout << " Read nodes " << in0 << ", " << in1 << ", " << in2
275 << ", " << in3 << ", ... from element " << il + 1 << " of "
276 << nElements << " with mat " << imat << ".\n";
277 }
278
279 // Check synchronisation.
280 if (readerror) {
281 PrintErrorReadingFile(hdr, elist, il);
282 felems.close();
283 return false;
284 }
285
286 // Check the material number and ensure that epsilon is non-negative.
287 if (imat < 0 || imat > (int)nMaterials) {
288 std::cerr << hdr << "\n Out-of-range material number on file " << elist
289 << " (line " << il << ").\n";
290 std::cerr << " Element: " << il << ", material: " << imat << "\n";
291 std::cerr << " nodes: (" << in0 << ", " << in1 << ", " << in2 << ", "
292 << in3 << ", " << in4 << ", " << in5 << ", " << in6 << ", "
293 << in7 << ", " << in8 << ", " << in9 << ")\n";
294 ok = false;
295 }
296 if (m_materials[imat].eps < 0) {
297 std::cerr << hdr << "\n Element " << il << " in element list " << elist
298 << "\n uses material " << imat
299 << " which has not been assigned a positive permittivity in "
300 << mplist << ".\n";
301 ok = false;
302 }
303
304 // Check the node numbers.
305 if (in0 < 1 || in1 < 1 || in2 < 1 || in3 < 1 || in4 < 1 || in5 < 1 ||
306 in6 < 1 || in7 < 1 || in8 < 1 || in9 < 1) {
307 std::cerr << hdr << "\n Found a node number < 1 on file " << elist
308 << " (line " << il << ").\n Element: " << il
309 << ", material: " << imat << "\n nodes: (" << in0 << ", "
310 << in1 << ", " << in2 << ", " << in3 << ", " << in4 << ", "
311 << in5 << ", " << in6 << ", " << in7 << ", " << in8 << ", "
312 << in9 << ")\n";
313 ok = false;
314 }
315 if (in0 > highestnode) highestnode = in0;
316 if (in1 > highestnode) highestnode = in1;
317 if (in2 > highestnode) highestnode = in2;
318 if (in3 > highestnode) highestnode = in3;
319 if (in4 > highestnode) highestnode = in4;
320 if (in5 > highestnode) highestnode = in5;
321 if (in6 > highestnode) highestnode = in6;
322 if (in7 > highestnode) highestnode = in7;
323 if (in8 > highestnode) highestnode = in8;
324 if (in9 > highestnode) highestnode = in9;
325
326 // These elements must not be degenerate.
327 if (in0 == in1 || in0 == in2 || in0 == in3 || in0 == in4 || in0 == in5 ||
328 in0 == in6 || in0 == in7 || in0 == in8 || in0 == in9 || in1 == in2 ||
329 in1 == in3 || in1 == in4 || in1 == in5 || in1 == in6 || in1 == in7 ||
330 in1 == in8 || in1 == in9 || in2 == in3 || in2 == in4 || in2 == in5 ||
331 in2 == in6 || in2 == in7 || in2 == in8 || in2 == in9 || in3 == in4 ||
332 in3 == in5 || in3 == in6 || in3 == in7 || in3 == in8 || in3 == in9 ||
333 in4 == in5 || in4 == in6 || in4 == in7 || in4 == in8 || in4 == in9 ||
334 in5 == in6 || in5 == in7 || in5 == in8 || in5 == in9 || in6 == in7 ||
335 in6 == in8 || in6 == in9 || in7 == in8 || in7 == in9 || in8 == in9) {
336 std::cerr << hdr << "\n Element " << il << " of file " << elist
337 << " is degenerate,\n"
338 << " no such elements are allowed in this type of map.\n";
339 ok = false;
340 }
341 Element newElement;
342 newElement.degenerate = false;
343
344 // Store the material reference.
345 newElement.matmap = imat;
346
347 // Node references
348 newElement.emap[0] = in0 - 1;
349 newElement.emap[1] = in1 - 1;
350 newElement.emap[2] = in2 - 1;
351 newElement.emap[3] = in3 - 1;
352 newElement.emap[4] = in4 - 1;
353 newElement.emap[7] = in5 - 1;
354 newElement.emap[5] = in6 - 1;
355 newElement.emap[6] = in7 - 1;
356 newElement.emap[8] = in8 - 1;
357 newElement.emap[9] = in9 - 1;
358 m_elements.push_back(std::move(newElement));
359 }
360
361 // Close the elements file.
362 felems.close();
363
364 // Set the ready flag.
365 if (ok) {
366 m_ready = true;
367 } else {
368 std::cerr << hdr << "\n Field map could not be "
369 << "read and cannot be interpolated.\n";
370 return false;
371 }
372
373 std::cout << hdr << " Finished.\n";
374
375 Prepare();
376 return true;
377}
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

Referenced by ComponentElmer().

◆ SetWeightingField()

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

Import a list of voltages to be used as weighting field.

Definition at line 379 of file ComponentElmer.cc.

379 {
380 const std::string hdr = m_className + "::SetWeightingField:";
381 if (!m_ready) {
382 PrintNotReady("SetWeightingField");
383 std::cerr << " Weighting field cannot be added.\n";
384 return false;
385 }
386
387 // Open the voltage list.
388 std::ifstream fwvolt;
389 fwvolt.open(wvolt.c_str(), std::ios::in);
390 if (fwvolt.fail()) {
391 PrintCouldNotOpen("SetWeightingField", wvolt);
392 return false;
393 }
394
395 // Check if a weighting field with the same label already exists.
396 const size_t iw = GetOrCreateWeightingFieldIndex(label);
397 if (iw + 1 != m_wfields.size()) {
398 std::cout << m_className << "::SetWeightingField:\n"
399 << " Replacing existing weighting field " << label << ".\n";
400 }
401 m_wfieldsOk[iw] = false;
402
403 // Temporary variables for use in file reading
404 constexpr int size = 100;
405 char line[size];
406 char* token = NULL;
407 bool readerror = false;
408 bool readstop = false;
409 int il = 1;
410
411 // Read past the header.
412 while (!readstop && fwvolt.getline(line, size, '\n')) {
413 token = strtok(line, " ");
414 if (strcmp(token, "Perm:") == 0) readstop = true;
415 il++;
416 }
417
418 // Should have stopped: if not, print error message.
419 if (!readstop) {
420 std::cerr << hdr << "\n Error reading past header of potentials file "
421 << wvolt << ".\n";
422 fwvolt.close();
423 return false;
424 }
425
426 // Read past the permutation map (number of lines = nNodes).
427 const int nNodes = m_nodes.size();
428 for (int tl = 0; tl < nNodes; tl++) {
429 fwvolt.getline(line, size, '\n');
430 il++;
431 }
432
433 // Read the potentials.
434 for (int tl = 0; tl < nNodes; tl++) {
435 double v;
436 fwvolt.getline(line, size, '\n');
437 token = strtok(line, " ");
438 v = ReadDouble(token, -1, readerror);
439 if (readerror) {
440 PrintErrorReadingFile(hdr, wvolt, il);
441 fwvolt.close();
442 return false;
443 }
444 // Place the weighting potential at its appropriate node and index.
445 m_nodes[tl].w[iw] = v;
446 }
447
448 // Close the potentials file.
449 fwvolt.close();
450 std::cout << hdr << "\n Read potentials from file " << wvolt << ".\n";
451
452 // Set the ready flag.
453 m_wfieldsOk[iw] = true;
454 return true;
455}
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::ComponentElmer::UpdatePeriodicity ( )
inlineoverrideprotectedvirtual

Verify periodicities.

Implements Garfield::Component.

Definition at line 57 of file ComponentElmer.hh.

57{ UpdatePeriodicityCommon(); }

◆ WeightingField()

void Garfield::ComponentElmer::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 569 of file ComponentElmer.cc.

571 {
572 // Initial values
573 wx = wy = wz = 0;
574
575 // Do not proceed if not properly initialised.
576 if (!m_ready) return;
577
578 // Look for the label.
579 const size_t iw = GetWeightingFieldIndex(label);
580 // Do not proceed if the requested weighting field does not exist.
581 if (iw == m_wfields.size()) return;
582 // Check if the weighting field is properly initialised.
583 if (!m_wfieldsOk[iw]) return;
584
585 // Copy the coordinates.
586 double x = xin, y = yin, z = zin;
587
588 // Map the coordinates onto field map coordinates
589 bool xmirr, ymirr, zmirr;
590 double rcoordinate, rotation;
591 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
592
593 if (m_warning) PrintWarning("WeightingField");
594
595 // Find the element that contains this point.
596 double t1, t2, t3, t4, jac[4][4], det;
597 const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
598 // Check if the point is in the mesh.
599 if (imap < 0) return;
600
601 const Element& element = m_elements[imap];
602 if (m_debug) {
603 PrintElement("WeightingField", x, y, z, t1, t2, t3, t4, element, 10, iw);
604 }
605 const Node& n0 = m_nodes[element.emap[0]];
606 const Node& n1 = m_nodes[element.emap[1]];
607 const Node& n2 = m_nodes[element.emap[2]];
608 const Node& n3 = m_nodes[element.emap[3]];
609 const Node& n4 = m_nodes[element.emap[4]];
610 const Node& n5 = m_nodes[element.emap[5]];
611 const Node& n6 = m_nodes[element.emap[6]];
612 const Node& n7 = m_nodes[element.emap[7]];
613 const Node& n8 = m_nodes[element.emap[8]];
614 const Node& n9 = m_nodes[element.emap[9]];
615 // Shorthands.
616 const double fourt1 = 4 * t1;
617 const double fourt2 = 4 * t2;
618 const double fourt3 = 4 * t3;
619 const double fourt4 = 4 * t4;
620 const double invdet = 1. / det;
621 // Tetrahedral field
622 wx = -(n0.w[iw] * (fourt1 - 1) * jac[0][1] +
623 n1.w[iw] * (fourt2 - 1) * jac[1][1] +
624 n2.w[iw] * (fourt3 - 1) * jac[2][1] +
625 n3.w[iw] * (fourt4 - 1) * jac[3][1] +
626 n4.w[iw] * (fourt2 * jac[0][1] + fourt1 * jac[1][1]) +
627 n5.w[iw] * (fourt3 * jac[0][1] + fourt1 * jac[2][1]) +
628 n6.w[iw] * (fourt4 * jac[0][1] + fourt1 * jac[3][1]) +
629 n7.w[iw] * (fourt3 * jac[1][1] + fourt2 * jac[2][1]) +
630 n8.w[iw] * (fourt4 * jac[1][1] + fourt2 * jac[3][1]) +
631 n9.w[iw] * (fourt4 * jac[2][1] + fourt3 * jac[3][1])) *
632 invdet;
633 wy = -(n0.w[iw] * (fourt1 - 1) * jac[0][2] +
634 n1.w[iw] * (fourt2 - 1) * jac[1][2] +
635 n2.w[iw] * (fourt3 - 1) * jac[2][2] +
636 n3.w[iw] * (fourt4 - 1) * jac[3][2] +
637 n4.w[iw] * (fourt2 * jac[0][2] + fourt1 * jac[1][2]) +
638 n5.w[iw] * (fourt3 * jac[0][2] + fourt1 * jac[2][2]) +
639 n6.w[iw] * (fourt4 * jac[0][2] + fourt1 * jac[3][2]) +
640 n7.w[iw] * (fourt3 * jac[1][2] + fourt2 * jac[2][2]) +
641 n8.w[iw] * (fourt4 * jac[1][2] + fourt2 * jac[3][2]) +
642 n9.w[iw] * (fourt4 * jac[2][2] + fourt3 * jac[3][2])) *
643 invdet;
644 wz = -(n0.w[iw] * (fourt1 - 1) * jac[0][3] +
645 n1.w[iw] * (fourt2 - 1) * jac[1][3] +
646 n2.w[iw] * (fourt3 - 1) * jac[2][3] +
647 n3.w[iw] * (fourt4 - 1) * jac[3][3] +
648 n4.w[iw] * (fourt2 * jac[0][3] + fourt1 * jac[1][3]) +
649 n5.w[iw] * (fourt3 * jac[0][3] + fourt1 * jac[2][3]) +
650 n6.w[iw] * (fourt4 * jac[0][3] + fourt1 * jac[3][3]) +
651 n7.w[iw] * (fourt3 * jac[1][3] + fourt2 * jac[2][3]) +
652 n8.w[iw] * (fourt4 * jac[1][3] + fourt2 * jac[3][3]) +
653 n9.w[iw] * (fourt4 * jac[2][3] + fourt3 * jac[3][3])) *
654 invdet;
655
656 // Transform field to global coordinates
657 UnmapFields(wx, wy, wz, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
658}
Garfield::ComponentFieldMap::GetWeightingFieldIndex
size_t GetWeightingFieldIndex(const std::string &label) const
Definition: ComponentFieldMap.cc:2348

◆ WeightingPotential()

double Garfield::ComponentElmer::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 660 of file ComponentElmer.cc.

662 {
663 // Do not proceed if not properly initialised.
664 if (!m_ready) return 0.;
665
666 // Look for the label.
667 const size_t iw = GetWeightingFieldIndex(label);
668 // Do not proceed if the requested weighting field does not exist.
669 if (iw == m_wfields.size()) return 0.;
670 // Check if the weighting field is properly initialised.
671 if (!m_wfieldsOk[iw]) return 0.;
672
673 // Copy the coordinates.
674 double x = xin, y = yin, z = zin;
675
676 // Map the coordinates onto field map coordinates.
677 bool xmirr, ymirr, zmirr;
678 double rcoordinate, rotation;
679 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
680
681 if (m_warning) PrintWarning("WeightingPotential");
682
683 // Find the element that contains this point.
684 double t1, t2, t3, t4, jac[4][4], det;
685 const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
686 if (imap < 0) return 0.;
687
688 const Element& element = m_elements[imap];
689 if (m_debug) {
690 PrintElement("WeightingPotential", x, y, z, t1, t2, t3, t4, element, 10,
691 iw);
692 }
693 const Node& n0 = m_nodes[element.emap[0]];
694 const Node& n1 = m_nodes[element.emap[1]];
695 const Node& n2 = m_nodes[element.emap[2]];
696 const Node& n3 = m_nodes[element.emap[3]];
697 const Node& n4 = m_nodes[element.emap[4]];
698 const Node& n5 = m_nodes[element.emap[5]];
699 const Node& n6 = m_nodes[element.emap[6]];
700 const Node& n7 = m_nodes[element.emap[7]];
701 const Node& n8 = m_nodes[element.emap[8]];
702 const Node& n9 = m_nodes[element.emap[9]];
703 // Tetrahedral field
704 return n0.w[iw] * t1 * (2 * t1 - 1) + n1.w[iw] * t2 * (2 * t2 - 1) +
705 n2.w[iw] * t3 * (2 * t3 - 1) + n3.w[iw] * t4 * (2 * t4 - 1) +
706 4 * n4.w[iw] * t1 * t2 + 4 * n5.w[iw] * t1 * t3 +
707 4 * n6.w[iw] * t1 * t4 + 4 * n7.w[iw] * t2 * t3 +
708 4 * n8.w[iw] * t2 * t4 + 4 * n9.w[iw] * t3 * t4;
709}
The documentation for this class was generated from the following files: