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

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

#include <ComponentElmer2D.hh>

+ Inheritance diagram for Garfield::ComponentElmer2D:

Public Member Functions

 ComponentElmer2D ()
 Default constructor.
 
 ComponentElmer2D (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().
 
 ~ComponentElmer2D ()
 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.
 
void SetRangeZ (const double zmin, const double zmax)
 
- 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 ComponentElmer2D.hh.

Constructor & Destructor Documentation

◆ ComponentElmer2D() [1/2]

Garfield::ComponentElmer2D::ComponentElmer2D ( )

Default constructor.

Definition at line 20 of file ComponentElmer2D.cc.

20 : ComponentFieldMap("Elmer2D") {
21 m_is3d = false;
22
23 // Default bounding box
24 m_minBoundingBox[2] = -50;
25 m_maxBoundingBox[2] = 50;
26}
Garfield::ComponentFieldMap::m_maxBoundingBox
std::array< double, 3 > m_maxBoundingBox
Definition: ComponentFieldMap.hh:149
Garfield::ComponentFieldMap::m_minBoundingBox
std::array< double, 3 > m_minBoundingBox
Definition: ComponentFieldMap.hh:148
Garfield::ComponentFieldMap::ComponentFieldMap
ComponentFieldMap()=delete
Default constructor.

◆ ComponentElmer2D() [2/2]

Garfield::ComponentElmer2D::ComponentElmer2D ( 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 28 of file ComponentElmer2D.cc.

34 : ComponentElmer2D() {
35
36 Initialise(header, elist, nlist, mplist, volt, unit);
37}
Garfield::ComponentElmer2D::ComponentElmer2D
ComponentElmer2D()
Default constructor.
Definition: ComponentElmer2D.cc:20
Garfield::ComponentElmer2D::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: ComponentElmer2D.cc:39

◆ ~ComponentElmer2D()

Garfield::ComponentElmer2D::~ComponentElmer2D ( )
inline

Destructor.

Definition at line 21 of file ComponentElmer2D.hh.

21{}

Member Function Documentation

◆ ElectricField() [1/2]

void Garfield::ComponentElmer2D::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 482 of file ComponentElmer2D.cc.

485 {
486 // Copy the coordinates
487 double x = xin, y = yin, z = 0.;
488
489 // Map the coordinates onto field map coordinates
490 bool xmirr, ymirr, zmirr;
491 double rcoordinate, rotation;
492 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
493
494 // Initial values
495 ex = ey = ez = volt = 0.;
496 status = 0;
497 m = nullptr;
498
499 // Do not proceed if not properly initialised.
500 if (!m_ready) {
501 status = -10;
502 PrintNotReady("ElectricField");
503 return;
504 }
505
506 if (m_warning) PrintWarning("ElectricField");
507
508 // Check that the input z-coordinate is within the established bounding box.
509 if (zin < m_minBoundingBox[2] || zin > m_maxBoundingBox[2]) {
510 status = -5;
511 return;
512 }
513
514 // Find the element that contains this point
515 double t1, t2, t3, t4, jac[4][4], det;
516 const int imap = FindElement5(x, y, z, t1, t2, t3, t4, jac, det);
517 if (imap < 0) {
518 if (m_debug) {
519 std::cout << m_className << "::ElectricField:\n Point (" << x << ", "
520 << y << ", " << z << ") is not in the mesh.\n";
521 }
522 status = -6;
523 return;
524 }
525
526 const Element& element = m_elements[imap];
527 if (m_debug) {
528 PrintElement("ElectricField", x, y, z, t1, t2, t3, t4, element, 10);
529 }
530 const Node& n0 = m_nodes[element.emap[0]];
531 const Node& n1 = m_nodes[element.emap[1]];
532 const Node& n2 = m_nodes[element.emap[2]];
533 const Node& n3 = m_nodes[element.emap[3]];
534 const Node& n4 = m_nodes[element.emap[4]];
535 const Node& n5 = m_nodes[element.emap[5]];
536 const Node& n6 = m_nodes[element.emap[6]];
537 const Node& n7 = m_nodes[element.emap[7]];
538
539 // Shorthands.
540 const double invdet = 1. / det;
541 volt = -n0.v * (1 - t1) * (1 - t2) * (1 + t1 + t2) * 0.25 -
542 n1.v * (1 + t1) * (1 - t2) * (1 - t1 + t2) * 0.25 -
543 n2.v * (1 + t1) * (1 + t2) * (1 - t1 - t2) * 0.25 -
544 n3.v * (1 - t1) * (1 + t2) * (1 + t1 - t2) * 0.25 +
545 n4.v * (1 - t1) * (1 + t1) * (1 - t2) * 0.5 +
546 n5.v * (1 + t1) * (1 + t2) * (1 - t2) * 0.5 +
547 n6.v * (1 - t1) * (1 + t1) * (1 + t2) * 0.5 +
548 n7.v * (1 - t1) * (1 + t2) * (1 - t2) * 0.5;
549 ex = -(n0.v * ((1 - t2) * (2 * t1 + t2) * jac[0][0] +
550 (1 - t1) * (t1 + 2 * t2) * jac[1][0]) *
551 0.25 +
552 n1.v * ((1 - t2) * (2 * t1 - t2) * jac[0][0] -
553 (1 + t1) * (t1 - 2 * t2) * jac[1][0]) *
554 0.25 +
555 n2.v * ((1 + t2) * (2 * t1 + t2) * jac[0][0] +
556 (1 + t1) * (t1 + 2 * t2) * jac[1][0]) *
557 0.25 +
558 n3.v * ((1 + t2) * (2 * t1 - t2) * jac[0][0] -
559 (1 - t1) * (t1 - 2 * t2) * jac[1][0]) *
560 0.25 +
561 n4.v * (t1 * (t2 - 1) * jac[0][0] +
562 (t1 - 1) * (t1 + 1) * jac[1][0] * 0.5) +
563 n5.v * ((1 - t2) * (1 + t2) * jac[0][0] * 0.5 -
564 (1 + t1) * t2 * jac[1][0]) +
565 n6.v * (-t1 * (1 + t2) * jac[0][0] +
566 (1 - t1) * (1 + t1) * jac[1][0] * 0.5) +
567 n7.v * ((t2 - 1) * (t2 + 1) * jac[0][0] * 0.5 +
568 (t1 - 1) * t2 * jac[1][0])) *
569 invdet;
570 ey = -(n0.v * ((1 - t2) * (2 * t1 + t2) * jac[0][1] +
571 (1 - t1) * (t1 + 2 * t2) * jac[1][1]) *
572 0.25 +
573 n1.v * ((1 - t2) * (2 * t1 - t2) * jac[0][1] -
574 (1 + t1) * (t1 - 2 * t2) * jac[1][1]) *
575 0.25 +
576 n2.v * ((1 + t2) * (2 * t1 + t2) * jac[0][1] +
577 (1 + t1) * (t1 + 2 * t2) * jac[1][1]) *
578 0.25 +
579 n3.v * ((1 + t2) * (2 * t1 - t2) * jac[0][1] -
580 (1 - t1) * (t1 - 2 * t2) * jac[1][1]) *
581 0.25 +
582 n4.v * (t1 * (t2 - 1) * jac[0][1] +
583 (t1 - 1) * (t1 + 1) * jac[1][1] * 0.5) +
584 n5.v * ((1 - t2) * (1 + t2) * jac[0][1] * 0.5 -
585 (1 + t1) * t2 * jac[1][1]) +
586 n6.v * (-t1 * (1 + t2) * jac[0][1] +
587 (1 - t1) * (1 + t1) * jac[1][1] * 0.5) +
588 n7.v * ((t2 - 1) * (t2 + 1) * jac[0][1] * 0.5 +
589 (t1 - 1) * t2 * jac[1][1])) *
590 invdet;
591
592 // Transform field to global coordinates
593 UnmapFields(ex, ey, ez, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
594
595 // Drift medium?
596 const Material& mat = m_materials[element.matmap];
597 if (m_debug) {
598 std::cout << m_className << "::ElectricField:\n Material "
599 << element.matmap << ", drift flag " << mat.driftmedium << ".\n";
600 }
601 m = mat.medium;
602 status = -5;
603 if (mat.driftmedium && m && m->IsDriftable()) status = 0;
604}
Garfield::ComponentFieldMap::MapCoordinates
void MapCoordinates(double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
Move (xpos, ypos, zpos) to field map coordinates.
Definition: ComponentFieldMap.cc:2021
Garfield::ComponentFieldMap::PrintWarning
void PrintWarning(const std::string &header)
Definition: ComponentFieldMap.cc:2375
Garfield::ComponentFieldMap::m_materials
std::vector< Material > m_materials
Definition: ComponentFieldMap.hh:138
Garfield::ComponentFieldMap::FindElement5
int FindElement5(const double x, const double y, const double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det)
Find the element for a point in curved quadratic quadrilaterals.
Definition: ComponentFieldMap.cc:190
Garfield::ComponentFieldMap::m_elements
std::vector< Element > m_elements
Definition: ComponentFieldMap.hh:113
Garfield::ComponentFieldMap::PrintElement
void PrintElement(const std::string &header, const double x, const double y, const double z, const double t1, const double t2, const double t3, const double t4, const Element &element, const unsigned int n, const int iw=-1) const
Definition: ComponentFieldMap.cc:2395
Garfield::ComponentFieldMap::UnmapFields
void UnmapFields(double &ex, double &ey, double &ez, double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
Move (ex, ey, ez) to global coordinates.
Definition: ComponentFieldMap.cc:2138
Garfield::ComponentFieldMap::PrintNotReady
void PrintNotReady(const std::string &header) const
Definition: ComponentFieldMap.cc:2383
Garfield::Component::m_debug
bool m_debug
Switch on/off debugging messages.
Definition: Component.hh:341
Garfield::Component::m_className
std::string m_className
Class name.
Definition: Component.hh:329
Garfield::Component::m_ready
bool m_ready
Ready for use?
Definition: Component.hh:338
Element
Definition: neBEM.h:155

◆ ElectricField() [2/2]

void Garfield::ComponentElmer2D::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 475 of file ComponentElmer2D.cc.

477 {
478 double v = 0.;
479 ElectricField(x, y, z, ex, ey, ez, v, m, status);
480}
Garfield::ComponentElmer2D::ElectricField
void ElectricField(const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status) override
Definition: ComponentElmer2D.cc:475

Referenced by ElectricField().

◆ GetAspectRatio()

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

Implements Garfield::ComponentFieldMap.

Definition at line 836 of file ComponentElmer2D.cc.

837 {
838 if (i >= m_elements.size()) {
839 dmin = dmax = 0.;
840 return;
841 }
842
843 const Element& element = m_elements[i];
844 const int np = 8;
845 // Loop over all pairs of vertices.
846 for (int j = 0; j < np - 1; ++j) {
847 const Node& nj = m_nodes[element.emap[j]];
848 for (int k = j + 1; k < np; ++k) {
849 const Node& nk = m_nodes[element.emap[k]];
850 // Compute distance.
851 const double dx = nj.x - nk.x;
852 const double dy = nj.y - nk.y;
853 const double dist = sqrt(dx * dx + dy * dy);
854 if (k == 1) {
855 dmin = dmax = dist;
856 } else {
857 if (dist < dmin) dmin = dist;
858 if (dist > dmax) dmax = dist;
859 }
860 }
861 }
862}
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314

◆ GetElementVolume()

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

Implements Garfield::ComponentFieldMap.

Definition at line 822 of file ComponentElmer2D.cc.

822 {
823 if (i >= m_elements.size()) return 0.;
824 const Element& element = m_elements[i];
825 const Node& n0 = m_nodes[element.emap[0]];
826 const Node& n1 = m_nodes[element.emap[1]];
827 const Node& n2 = m_nodes[element.emap[2]];
828 const Node& n3 = m_nodes[element.emap[3]];
829 const double surf =
830 0.5 *
831 (fabs((n1.x - n0.x) * (n2.y - n0.y) - (n2.x - n0.x) * (n1.y - n0.y)) +
832 fabs((n3.x - n0.x) * (n2.y - n0.y) - (n2.x - n0.x) * (n3.y - n0.y)));
833 return surf;
834}
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615

◆ GetMedium()

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

761 {
762 // Copy the coordinates
763 double x = xin, y = yin, z = 0.;
764
765 // Map the coordinates onto field map coordinates
766 bool xmirr, ymirr, zmirr;
767 double rcoordinate, rotation;
768 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
769
770 if (zin < m_minBoundingBox[2] || z > m_maxBoundingBox[2]) {
771 return nullptr;
772 }
773
774 // Do not proceed if not properly initialised.
775 if (!m_ready) {
776 PrintNotReady("GetMedium");
777 return nullptr;
778 }
779 if (m_warning) PrintWarning("GetMedium");
780
781 // Find the element that contains this point.
782 double t1, t2, t3, t4, jac[4][4], det;
783 const int imap = FindElement5(x, y, z, t1, t2, t3, t4, jac, det);
784 if (imap < 0) {
785 if (m_debug) {
786 std::cout << m_className << "::GetMedium:\n Point (" << x << ", " << y
787 << ", " << z << ") is not in the mesh.\n";
788 }
789 return nullptr;
790 }
791 const Element& element = m_elements[imap];
792 if (element.matmap >= m_materials.size()) {
793 if (m_debug) {
794 std::cerr << m_className << "::GetMedium:\n Point (" << x << ", " << y
795 << ", " << z << ") has out of range material number " << imap
796 << ".\n";
797 }
798 return nullptr;
799 }
800
801 if (m_debug) PrintElement("GetMedium", x, y, z, t1, t2, t3, t4, element, 10);
802
803 return m_materials[element.matmap].medium;
804}

◆ Initialise()

bool Garfield::ComponentElmer2D::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 39 of file ComponentElmer2D.cc.

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

Referenced by ComponentElmer2D().

◆ SetRangeZ()

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

Definition at line 807 of file ComponentElmer2D.cc.

807 {
808 if (fabs(zmax - zmin) <= 0.) {
809 std::cerr << m_className << "::SetRangeZ: Zero range is not permitted.\n";
810 return;
811 }
812 m_minBoundingBox[2] = m_mapmin[2] = std::min(zmin, zmax);
813 m_maxBoundingBox[2] = m_mapmax[2] = std::max(zmin, zmax);
814}
Garfield::ComponentFieldMap::m_mapmin
std::array< double, 3 > m_mapmin
Definition: ComponentFieldMap.hh:152
Garfield::ComponentFieldMap::m_mapmax
std::array< double, 3 > m_mapmax
Definition: ComponentFieldMap.hh:153

Referenced by main().

◆ SetWeightingField()

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

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

Definition at line 397 of file ComponentElmer2D.cc.

397 {
398 const std::string hdr = m_className + "::SetWeightingField:";
399 if (!m_ready) {
400 PrintNotReady("SetWeightingField");
401 std::cerr << " Weighting field cannot be added.\n";
402 return false;
403 }
404
405 // Open the voltage list.
406 std::ifstream fwvolt;
407 fwvolt.open(wvolt.c_str(), std::ios::in);
408 if (fwvolt.fail()) {
409 PrintCouldNotOpen("SetWeightingField", wvolt);
410 return false;
411 }
412
413 // Check if a weighting field with the same label already exists.
414 const size_t iw = GetOrCreateWeightingFieldIndex(label);
415 if (iw + 1 != m_wfields.size()) {
416 std::cout << m_className << "::SetWeightingField:\n"
417 << " Replacing existing weighting field " << label << ".\n";
418 }
419 m_wfieldsOk[iw] = false;
420
421 // Temporary variables for use in file reading
422 constexpr int size = 100;
423 char line[size];
424 char* token = NULL;
425 bool readerror = false;
426 bool readstop = false;
427 int il = 1;
428
429 // Read past the header.
430 while (!readstop && fwvolt.getline(line, size, '\n')) {
431 token = strtok(line, " ");
432 if (strcmp(token, "Perm:") == 0) readstop = true;
433 il++;
434 }
435
436 // Should have stopped: if not, print error message.
437 if (!readstop) {
438 std::cerr << hdr << "\n Error reading past header of potentials file "
439 << wvolt << ".\n";
440 fwvolt.close();
441 return false;
442 }
443
444 // Read past the permutation map (number of lines = nNodes).
445 const int nNodes = m_nodes.size();
446 for (int tl = 0; tl < nNodes; tl++) {
447 fwvolt.getline(line, size, '\n');
448 il++;
449 }
450
451 // Read the potentials.
452 for (int tl = 0; tl < nNodes; tl++) {
453 double v;
454 fwvolt.getline(line, size, '\n');
455 token = strtok(line, " ");
456 v = ReadDouble(token, -1, readerror);
457 if (readerror) {
458 PrintErrorReadingFile(hdr, wvolt, il);
459 fwvolt.close();
460 return false;
461 }
462 // Place the weighting potential at its appropriate node and index.
463 m_nodes[tl].w[iw] = v;
464 }
465
466 // Close the potentials file.
467 fwvolt.close();
468 std::cout << hdr << "\n Read potentials from file " << wvolt << ".\n";
469
470 // Set the ready flag.
471 m_wfieldsOk[iw] = true;
472 return true;
473}
Garfield::ComponentFieldMap::GetOrCreateWeightingFieldIndex
size_t GetOrCreateWeightingFieldIndex(const std::string &label)
Definition: ComponentFieldMap.cc:2357

◆ UpdatePeriodicity()

void Garfield::ComponentElmer2D::UpdatePeriodicity ( )
overrideprotectedvirtual

Verify periodicities.

Implements Garfield::Component.

Definition at line 816 of file ComponentElmer2D.cc.

816 {
817 UpdatePeriodicity2d();
818 UpdatePeriodicityCommon();
819}

◆ WeightingField()

void Garfield::ComponentElmer2D::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 606 of file ComponentElmer2D.cc.

608 {
609 // Initial values
610 wx = wy = wz = 0;
611
612 // Do not proceed if not properly initialised.
613 if (!m_ready) return;
614
615 // Look for the label.
616 const size_t iw = GetWeightingFieldIndex(label);
617 // Do not proceed if the requested weighting field does not exist.
618 if (iw == m_wfields.size()) return;
619 // Check if the weighting field is properly initialised.
620 if (!m_wfieldsOk[iw]) return;
621
622 // Copy the coordinates.
623 double x = xin, y = yin, z = 0.;
624
625 // Map the coordinates onto field map coordinates
626 bool xmirr, ymirr, zmirr;
627 double rcoordinate, rotation;
628 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
629
630 if (m_warning) PrintWarning("WeightingField");
631
632 // Check that the input z-coordinate is within the established bounding box.
633 if (zin < m_minBoundingBox[2] || zin > m_maxBoundingBox[2]) {
634 return;
635 }
636
637 // Find the element that contains this point.
638 double t1, t2, t3, t4, jac[4][4], det;
639 const int imap = FindElement5(x, y, z, t1, t2, t3, t4, jac, det);
640 // Check if the point is in the mesh.
641 if (imap < 0) return;
642
643 const Element& element = m_elements[imap];
644 if (m_debug) {
645 PrintElement("WeightingField", x, y, z, t1, t2, t3, t4, element, 10, iw);
646 }
647 const Node& n0 = m_nodes[element.emap[0]];
648 const Node& n1 = m_nodes[element.emap[1]];
649 const Node& n2 = m_nodes[element.emap[2]];
650 const Node& n3 = m_nodes[element.emap[3]];
651 const Node& n4 = m_nodes[element.emap[4]];
652 const Node& n5 = m_nodes[element.emap[5]];
653 const Node& n6 = m_nodes[element.emap[6]];
654 const Node& n7 = m_nodes[element.emap[7]];
655 // Shorthands.
656 const double invdet = 1. / det;
657 wx = -(n0.w[iw] * ((1 - t2) * (2 * t1 + t2) * jac[0][0] +
658 (1 - t1) * (t1 + 2 * t2) * jac[1][0]) *
659 0.25 +
660 n1.w[iw] * ((1 - t2) * (2 * t1 - t2) * jac[0][0] -
661 (1 + t1) * (t1 - 2 * t2) * jac[1][0]) *
662 0.25 +
663 n2.w[iw] * ((1 + t2) * (2 * t1 + t2) * jac[0][0] +
664 (1 + t1) * (t1 + 2 * t2) * jac[1][0]) *
665 0.25 +
666 n3.w[iw] * ((1 + t2) * (2 * t1 - t2) * jac[0][0] -
667 (1 - t1) * (t1 - 2 * t2) * jac[1][0]) *
668 0.25 +
669 n4.w[iw] * (t1 * (t2 - 1) * jac[0][0] +
670 (t1 - 1) * (t1 + 1) * jac[1][0] * 0.5) +
671 n5.w[iw] * ((1 - t2) * (1 + t2) * jac[0][0] * 0.5 -
672 (1 + t1) * t2 * jac[1][0]) +
673 n6.w[iw] * (-t1 * (1 + t2) * jac[0][0] +
674 (1 - t1) * (1 + t1) * jac[1][0] * 0.5) +
675 n7.w[iw] * ((t2 - 1) * (1 + t2) * jac[0][0] * 0.5 +
676 (t1 - 1) * t2 * jac[1][0])) *
677 invdet;
678 wy = -(n0.w[iw] * ((1 - t2) * (2 * t1 + t2) * jac[0][1] +
679 (1 - t1) * (t1 + 2 * t2) * jac[1][1]) *
680 0.25 +
681 n1.w[iw] * ((1 - t2) * (2 * t1 - t2) * jac[0][1] -
682 (1 + t1) * (t1 - 2 * t2) * jac[1][1]) *
683 0.25 +
684 n2.w[iw] * ((1 + t2) * (2 * t1 + t2) * jac[0][1] +
685 (1 + t1) * (t1 + 2 * t2) * jac[1][1]) *
686 0.25 +
687 n3.w[iw] * ((1 + t2) * (2 * t1 - t2) * jac[0][1] -
688 (1 - t1) * (t1 - 2 * t2) * jac[1][1]) *
689 0.25 +
690 n4.w[iw] * (t1 * (t2 - 1) * jac[0][1] +
691 (t1 - 1) * (t1 + 1) * jac[1][1] * 0.5) +
692 n5.w[iw] * ((1 - t2) * (1 + t2) * jac[0][1] * 0.5 -
693 (1 + t1) * t2 * jac[1][1]) +
694 n6.w[iw] * (-t1 * (1 + t2) * jac[0][1] +
695 (1 - t1) * (1 + t1) * jac[1][1] * 0.5) +
696 n7.w[iw] * ((t2 - 1) * (t2 + 1) * jac[0][1] * 0.5 +
697 (t1 - 1) * t2 * jac[1][1])) *
698 invdet;
699
700 // Transform field to global coordinates
701 UnmapFields(wx, wy, wz, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
702}
Garfield::ComponentFieldMap::GetWeightingFieldIndex
size_t GetWeightingFieldIndex(const std::string &label) const
Definition: ComponentFieldMap.cc:2348

◆ WeightingPotential()

double Garfield::ComponentElmer2D::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 704 of file ComponentElmer2D.cc.

706 {
707 // Do not proceed if not properly initialised.
708 if (!m_ready) return 0.;
709
710 // Look for the label.
711 const size_t iw = GetWeightingFieldIndex(label);
712 // Do not proceed if the requested weighting field does not exist.
713 if (iw == m_wfields.size()) return 0.;
714 // Check if the weighting field is properly initialised.
715 if (!m_wfieldsOk[iw]) return 0.;
716
717 // Copy the coordinates.
718 double x = xin, y = yin, z = 0.;
719
720 // Map the coordinates onto field map coordinates.
721 bool xmirr, ymirr, zmirr;
722 double rcoordinate, rotation;
723 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
724
725 if (m_warning) PrintWarning("WeightingPotential");
726
727 // Check that the input z-coordinate is within the established bounding box.
728 if (zin < m_minBoundingBox[2] || zin > m_maxBoundingBox[2]) {
729 return 0.;
730 }
731
732 // Find the element that contains this point.
733 double t1, t2, t3, t4, jac[4][4], det;
734 const int imap = FindElement5(x, y, z, t1, t2, t3, t4, jac, det);
735 if (imap < 0) return 0.;
736
737 const Element& element = m_elements[imap];
738 if (m_debug) {
739 PrintElement("WeightingPotential", x, y, z, t1, t2, t3, t4, element, 10,
740 iw);
741 }
742 const Node& n0 = m_nodes[element.emap[0]];
743 const Node& n1 = m_nodes[element.emap[1]];
744 const Node& n2 = m_nodes[element.emap[2]];
745 const Node& n3 = m_nodes[element.emap[3]];
746 const Node& n4 = m_nodes[element.emap[4]];
747 const Node& n5 = m_nodes[element.emap[5]];
748 const Node& n6 = m_nodes[element.emap[6]];
749 const Node& n7 = m_nodes[element.emap[7]];
750 return -n0.w[iw] * (1 - t1) * (1 - t2) * (1 + t1 + t2) * 0.25 -
751 n1.w[iw] * (1 + t1) * (1 - t2) * (1 - t1 + t2) * 0.25 -
752 n2.w[iw] * (1 + t1) * (1 + t2) * (1 - t1 - t2) * 0.25 -
753 n3.w[iw] * (1 - t1) * (1 + t2) * (1 + t1 - t2) * 0.25 +
754 n4.w[iw] * (1 - t1) * (1 + t1) * (1 - t2) * 0.5 +
755 n5.w[iw] * (1 + t1) * (1 + t2) * (1 - t2) * 0.5 +
756 n6.w[iw] * (1 - t1) * (1 + t1) * (1 + t2) * 0.5 +
757 n7.w[iw] * (1 - t1) * (1 + t2) * (1 - t2) * 0.5;
758}
The documentation for this class was generated from the following files: