Garfield++ v2r0
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)
 
void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status)
 
void WeightingField (const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label)
 
double WeightingPotential (const double x, const double y, const double z, const std::string &label)
 
Medium * GetMedium (const double x, const double y, const double z)
 Get the medium at a given location (x, y, z).
 
virtual bool IsInBoundingBox (const double x, const double y, const double z) const
 
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 ()
 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.
 
virtual bool IsInBoundingBox (const double x, const double y, const double z) const
 
virtual bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
 Get the bounding box coordinates.
 
virtual bool GetVoltageRange (double &vmin, double &vmax)
 Calculate the voltage range [V].
 
void PrintMaterials ()
 List all currently defined materials.
 
void DriftMedium (const unsigned int imat)
 Flag a field map material as a drift medium.
 
void NotDriftMedium (const unsigned int imat)
 Flag a field map materials as a non-drift medium.
 
unsigned int GetNumberOfMaterials () const
 Return the number of materials in the field map.
 
double GetPermittivity (const unsigned int imat) const
 Return the permittivity of a field map material.
 
double GetConductivity (const unsigned int imat) const
 Return the conductivity of a field map material.
 
void SetMedium (const unsigned int imat, Medium *medium)
 Associate a field map material with a Medium class.
 
Medium * GetMedium (const unsigned int i) const
 Return the Medium associated to a field map material.
 
Medium * GetMedium (const double x, const double y, const double z)=0
 Get the medium at a given location (x, y, z).
 
unsigned int GetNumberOfMedia () const
 
int GetNumberOfElements () const
 Return the number of mesh elements.
 
bool GetElement (const unsigned int i, double &vol, double &dmin, double &dmax)
 Return the volume and aspect ratio of a mesh element.
 
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
 
virtual void WeightingField (const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label)=0
 
virtual double WeightingPotential (const double x, const double y, const double z, const std::string &label)=0
 
void EnableCheckMapIndices ()
 
void DisableCheckMapIndices ()
 
void EnableDeleteBackgroundElements ()
 
void DisableDeleteBackgroundElements ()
 
void EnableTetrahedralTreeForElementSearch (const bool on=true)
 
- Public Member Functions inherited from Garfield::ComponentBase
 ComponentBase ()
 Constructor.
 
virtual ~ComponentBase ()
 Destructor.
 
virtual void SetGeometry (GeometryBase *geo)
 Define the geometry.
 
virtual void Clear ()
 Reset.
 
virtual Medium * GetMedium (const double x, const double y, const double z)
 Get the medium at a given location (x, y, z).
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status)=0
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status)=0
 
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 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 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)
 
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 DisablePeriodicityX ()
 
void EnablePeriodicityY (const bool on=true)
 Enable simple periodicity in the $y$ direction.
 
void DisablePeriodicityY ()
 
void EnablePeriodicityZ (const bool on=true)
 Enable simple periodicity in the $z$ direction.
 
void DisablePeriodicityZ ()
 
void EnableMirrorPeriodicityX (const bool on=true)
 Enable mirror periodicity in the $x$ direction.
 
void DisableMirrorPeriodicityX ()
 
void EnableMirrorPeriodicityY (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void DisableMirrorPeriodicityY ()
 
void EnableMirrorPeriodicityZ (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void DisableMirrorPeriodicityZ ()
 
void EnableAxialPeriodicityX (const bool on=true)
 Enable axial periodicity in the $x$ direction.
 
void DisableAxialPeriodicityX ()
 
void EnableAxialPeriodicityY (const bool on=true)
 Enable axial periodicity in the $y$ direction.
 
void DisableAxialPeriodicityY ()
 
void EnableAxialPeriodicityZ (const bool on=true)
 Enable axial periodicity in the $z$ direction.
 
void DisableAxialPeriodicityZ ()
 
void EnableRotationSymmetryX (const bool on=true)
 Enable rotation symmetry around the $x$ axis.
 
void DisableRotationSymmetryX ()
 
void EnableRotationSymmetryY (const bool on=true)
 Enable rotation symmetry around the $y$ axis.
 
void DisableRotationSymmetryY ()
 
void EnableRotationSymmetryZ (const bool on=true)
 Enable rotation symmetry around the $z$ axis.
 
void DisableRotationSymmetryZ ()
 
void EnableDebugging ()
 Switch on debugging messages.
 
void DisableDebugging ()
 Switch off debugging messages.
 
void ActivateTraps ()
 Request trapping to be taken care of by the component (for TCAD).
 
void DeactivateTraps ()
 
bool IsTrapActive ()
 
void ActivateVelocityMap ()
 Request velocity to be taken care of by the component (for TCAD).
 
void DectivateVelocityMap ()
 
bool IsVelocityActive ()
 
virtual bool ElectronAttachment (const double, const double, const double, double &eta)
 Get the electron attachment coefficient.
 
virtual bool HoleAttachment (const double, const double, const double, double &eta)
 Get the hole attachment coefficient.
 
virtual void ElectronVelocity (const double, const double, const double, double &vx, double &vy, double &vz, Medium *&, int &status)
 Get the electron drift velocity.
 
virtual void HoleVelocity (const double, const double, const double, double &vx, double &vy, double &vz, Medium *&, int &status)
 Get the hole drift velocity.
 
virtual bool GetElectronLifetime (const double, const double, const double, double &etau)
 
virtual bool GetHoleLifetime (const double, const double, const double, double &htau)
 

Protected Member Functions

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

Additional Inherited Members

- Protected Attributes inherited from Garfield::ComponentFieldMap
bool m_is3d
 
int nElements
 
std::vector< Element > elements
 
int nNodes
 
std::vector< Node > nodes
 
unsigned int m_nMaterials
 
std::vector< Material > materials
 
int nWeightingFields
 
std::vector< std::string > wfields
 
std::vector< bool > wfieldsOk
 
bool hasBoundingBox
 
double xMinBoundingBox
 
double yMinBoundingBox
 
double zMinBoundingBox
 
double xMaxBoundingBox
 
double yMaxBoundingBox
 
double zMaxBoundingBox
 
double mapxmin
 
double mapymin
 
double mapzmin
 
double mapxmax
 
double mapymax
 
double mapzmax
 
double mapxamin
 
double mapyamin
 
double mapzamin
 
double mapxamax
 
double mapyamax
 
double mapzamax
 
double mapvmin
 
double mapvmax
 
bool setangx
 
bool setangy
 
bool setangz
 
double mapsx
 
double mapsy
 
double mapsz
 
double cellsx
 
double cellsy
 
double cellsz
 
double mapnxa
 
double mapnya
 
double mapnza
 
bool m_deleteBackground
 
bool m_warning
 
unsigned int m_nWarnings
 
- Protected Attributes inherited from Garfield::ComponentBase
std::string m_className
 Class name.
 
GeometryBase * m_geometry
 Pointer to the geometry.
 
bool m_ready
 Ready for use?
 
bool m_activeTraps
 Does the component have traps?
 
bool m_hasVelocityMap
 Does the component have velocity maps?
 
bool m_xPeriodic
 Simple periodicity in x.
 
bool m_yPeriodic
 Simple periodicity in y.
 
bool m_zPeriodic
 Simple periodicity in z.
 
bool m_xMirrorPeriodic
 Mirror periodicity in x.
 
bool m_yMirrorPeriodic
 Mirror periodicity in y.
 
bool m_zMirrorPeriodic
 Mirror periodicity in z.
 
bool m_xAxiallyPeriodic
 Axial periodicity in x.
 
bool m_yAxiallyPeriodic
 Axial periodicity in y.
 
bool m_zAxiallyPeriodic
 Axial periodicity in z.
 
bool m_xRotationSymmetry
 Rotation symmetry around x-axis.
 
bool m_yRotationSymmetry
 Rotation symmetry around y-axis.
 
bool m_zRotationSymmetry
 Rotation symmetry around z-axis.
 
double m_bx0
 
double m_by0
 
double m_bz0
 
bool m_debug
 Switch on/off debugging messages.
 

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 30 of file ComponentElmer.cc.

30 : ComponentFieldMap() {
31
32 m_className = "ComponentElmer";
33}
Garfield::ComponentBase::m_className
std::string m_className
Class name.
Definition: ComponentBase.hh:233

◆ 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 35 of file ComponentElmer.cc.

41 : ComponentFieldMap() {
42
43 m_className = "ComponentElmer";
44 Initialise(header, elist, nlist, mplist, volt, unit);
45}
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:47

◆ ~ComponentElmer()

Garfield::ComponentElmer::~ComponentElmer ( )
inline

Destructor.

Definition at line 22 of file ComponentElmer.hh.

22{}

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

Implements Garfield::ComponentFieldMap.

Definition at line 567 of file ComponentElmer.cc.

570 {
571
572 // Copy the coordinates
573 double x = xin, y = yin, z = zin;
574
575 // Map the coordinates onto field map coordinates
576 bool xmirr, ymirr, zmirr;
577 double rcoordinate, rotation;
578 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
579
580 // Initial values
581 ex = ey = ez = volt = 0.;
582 status = 0;
583 m = NULL;
584
585 // Do not proceed if not properly initialised.
586 if (!m_ready) {
587 status = -10;
588 PrintNotReady("ElectricField");
589 return;
590 }
591
592 if (m_warning) PrintWarning("ElectricField");
593
594 // Find the element that contains this point
595 double t1, t2, t3, t4, jac[4][4], det;
596 const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
597 if (imap < 0) {
598 if (m_debug) {
599 std::cout << m_className << "::ElectricField:\n Point ("
600 << x << ", " << y << ", " << z << ") is not in the mesh.\n";
601 }
602 status = -6;
603 return;
604 }
605
606 if (m_debug) {
607 PrintElement("ElectricField", x, y, z, t1, t2, t3, t4, imap, 10);
608 }
609
610 const Element& element = elements[imap];
611 const Node& n0 = nodes[element.emap[0]];
612 const Node& n1 = nodes[element.emap[1]];
613 const Node& n2 = nodes[element.emap[2]];
614 const Node& n3 = nodes[element.emap[3]];
615 const Node& n4 = nodes[element.emap[4]];
616 const Node& n5 = nodes[element.emap[5]];
617 const Node& n6 = nodes[element.emap[6]];
618 const Node& n7 = nodes[element.emap[7]];
619 const Node& n8 = nodes[element.emap[8]];
620 const Node& n9 = nodes[element.emap[9]];
621 // Tetrahedral field
622 volt = n0.v * t1 * (2 * t1 - 1) + n1.v * t2 * (2 * t2 - 1) +
623 n2.v * t3 * (2 * t3 - 1) + n3.v * t4 * (2 * t4 - 1) +
624 4 * n4.v * t1 * t2 + 4 * n5.v * t1 * t3 + 4 * n6.v * t1 * t4 +
625 4 * n7.v * t2 * t3 + 4 * n8.v * t2 * t4 + 4 * n9.v * t3 * t4;
626 const double invdet = 1. / det;
627 ex = -(n0.v * (4 * t1 - 1) * jac[0][1] + n1.v * (4 * t2 - 1) * jac[1][1] +
628 n2.v * (4 * t3 - 1) * jac[2][1] + n3.v * (4 * t4 - 1) * jac[3][1] +
629 n4.v * (4 * t2 * jac[0][1] + 4 * t1 * jac[1][1]) +
630 n5.v * (4 * t3 * jac[0][1] + 4 * t1 * jac[2][1]) +
631 n6.v * (4 * t4 * jac[0][1] + 4 * t1 * jac[3][1]) +
632 n7.v * (4 * t3 * jac[1][1] + 4 * t2 * jac[2][1]) +
633 n8.v * (4 * t4 * jac[1][1] + 4 * t2 * jac[3][1]) +
634 n9.v * (4 * t4 * jac[2][1] + 4 * t3 * jac[3][1])) * invdet;
635 ey = -(n0.v * (4 * t1 - 1) * jac[0][2] + n1.v * (4 * t2 - 1) * jac[1][2] +
636 n2.v * (4 * t3 - 1) * jac[2][2] + n3.v * (4 * t4 - 1) * jac[3][2] +
637 n4.v * (4 * t2 * jac[0][2] + 4 * t1 * jac[1][2]) +
638 n5.v * (4 * t3 * jac[0][2] + 4 * t1 * jac[2][2]) +
639 n6.v * (4 * t4 * jac[0][2] + 4 * t1 * jac[3][2]) +
640 n7.v * (4 * t3 * jac[1][2] + 4 * t2 * jac[2][2]) +
641 n8.v * (4 * t4 * jac[1][2] + 4 * t2 * jac[3][2]) +
642 n9.v * (4 * t4 * jac[2][2] + 4 * t3 * jac[3][2])) * invdet;
643 ez = -(n0.v * (4 * t1 - 1) * jac[0][3] + n1.v * (4 * t2 - 1) * jac[1][3] +
644 n2.v * (4 * t3 - 1) * jac[2][3] + n3.v * (4 * t4 - 1) * jac[3][3] +
645 n4.v * (4 * t2 * jac[0][3] + 4 * t1 * jac[1][3]) +
646 n5.v * (4 * t3 * jac[0][3] + 4 * t1 * jac[2][3]) +
647 n6.v * (4 * t4 * jac[0][3] + 4 * t1 * jac[3][3]) +
648 n7.v * (4 * t3 * jac[1][3] + 4 * t2 * jac[2][3]) +
649 n8.v * (4 * t4 * jac[1][3] + 4 * t2 * jac[3][3]) +
650 n9.v * (4 * t4 * jac[2][3] + 4 * t3 * jac[3][3])) * invdet;
651
652 // Transform field to global coordinates
653 UnmapFields(ex, ey, ez, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
654
655 // Drift medium?
656 const Material& mat = materials[element.matmap];
657 if (m_debug) {
658 std::cout << m_className << "::ElectricField:\n Material "
659 << element.matmap << ", drift flag " << mat.driftmedium << ".\n";
660 }
661 m = mat.medium;
662 status = -5;
663 if (mat.driftmedium && m && m->IsDriftable()) status = 0;
664}
Garfield::ComponentBase::m_ready
bool m_ready
Ready for use?
Definition: ComponentBase.hh:239
Garfield::ComponentBase::m_debug
bool m_debug
Switch on/off debugging messages.
Definition: ComponentBase.hh:274
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:2441
Garfield::ComponentFieldMap::PrintWarning
void PrintWarning(const std::string &header)
Definition: ComponentFieldMap.hh:201
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 unsigned int i, const unsigned int n, const int iw=-1) const
Definition: ComponentFieldMap.cc:2740
Garfield::ComponentFieldMap::materials
std::vector< Material > materials
Definition: ComponentFieldMap.hh:133
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:331
Garfield::ComponentFieldMap::elements
std::vector< Element > elements
Definition: ComponentFieldMap.hh:108
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:2559
Garfield::ComponentFieldMap::PrintNotReady
void PrintNotReady(const std::string &header) const
Definition: ComponentFieldMap.hh:207

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

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

Definition at line 559 of file ComponentElmer.cc.

561 {
562
563 double v = 0.;
564 ElectricField(x, y, z, ex, ey, ez, v, m, status);
565}
Garfield::ComponentElmer::ElectricField
void ElectricField(const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status)
Definition: ComponentElmer.cc:559

Referenced by ElectricField().

◆ GetAspectRatio()

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

Implements Garfield::ComponentFieldMap.

Definition at line 886 of file ComponentElmer.cc.

887 {
888
889 if (i >= elements.size()) {
890 dmin = dmax = 0.;
891 return;
892 }
893
894 const Element& element = elements[i];
895 const int np = 4;
896 // Loop over all pairs of vertices.
897 for (int j = 0; j < np - 1; ++j) {
898 const Node& nj = nodes[element.emap[j]];
899 for (int k = j + 1; k < np; ++k) {
900 const Node& nk = nodes[element.emap[k]];
901 // Compute distance.
902 const double dx = nj.x - nk.x;
903 const double dy = nj.y - nk.y;
904 const double dz = nj.z - nk.z;
905 const double dist = sqrt(dx * dx + dy * dy + dz * dz);
906 if (k == 1) {
907 dmin = dmax = dist;
908 } else {
909 if (dist < dmin) dmin = dist;
910 if (dist > dmax) dmax = dist;
911 }
912 }
913 }
914}
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314

◆ GetElementVolume()

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

Implements Garfield::ComponentFieldMap.

Definition at line 864 of file ComponentElmer.cc.

864 {
865
866 if (i >= elements.size()) return 0.;
867 const Element& element = elements[i];
868 const Node& n0 = nodes[element.emap[0]];
869 const Node& n1 = nodes[element.emap[1]];
870 const Node& n2 = nodes[element.emap[2]];
871 const Node& n3 = nodes[element.emap[3]];
872
873 // Uses formula V = |a (dot) b x c|/6
874 // with a => "3", b => "1", c => "2" and origin "0"
875 const double vol =
876 fabs((n3.x - n0.x) *
877 ((n1.y - n0.y) * (n2.z - n0.z) - (n2.y - n0.y) * (n1.z - n0.z)) +
878 (n3.y - n0.y) *
879 ((n1.z - n0.z) * (n2.x - n0.x) - (n2.z - n0.z) * (n1.x - n0.x)) +
880 (n3.z - n0.z) * ((n1.x - n0.x) * (n2.y - n0.y) -
881 (n3.x - n0.x) * (n1.y - n0.y))) /
882 6.;
883 return vol;
884}
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 
)
virtual

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

Implements Garfield::ComponentFieldMap.

Definition at line 821 of file ComponentElmer.cc.

822 {
823
824 // Copy the coordinates
825 double x = xin, y = yin, z = zin;
826
827 // Map the coordinates onto field map coordinates
828 bool xmirr, ymirr, zmirr;
829 double rcoordinate, rotation;
830 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
831
832 // Do not proceed if not properly initialised.
833 if (!m_ready) {
834 PrintNotReady("GetMedium");
835 return NULL;
836 }
837 if (m_warning) PrintWarning("GetMedium");
838
839 // Find the element that contains this point
840 double t1, t2, t3, t4, jac[4][4], det;
841 const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
842 if (imap < 0) {
843 if (m_debug) {
844 std::cout << m_className << "::GetMedium:\n Point ("
845 << x << ", " << y << ", " << z << ") is not in the mesh.\n";
846 }
847 return NULL;
848 }
849 const Element& element = elements[imap];
850 if (element.matmap >= m_nMaterials) {
851 if (m_debug) {
852 std::cerr << m_className << "::GetMedium:\n Point ("
853 << x << ", " << y << ", " << z
854 << ") has out of range material number " << imap << ".\n";
855 }
856 return NULL;
857 }
858
859 if (m_debug) PrintElement("GetMedium", x, y, z, t1, t2, t3, t4, imap, 10);
860
861 return materials[element.matmap].medium;
862}

◆ 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 47 of file ComponentElmer.cc.

52 {
53
54 m_debug = false;
55 m_ready = false;
56 m_warning = false;
57 m_nWarnings = 0;
58
59 // Keep track of the success.
60 bool ok = true;
61
62 // Buffer for reading
63 const int size = 100;
64 char line[size];
65
66 // Open the header.
67 std::ifstream fheader;
68 fheader.open(header.c_str(), std::ios::in);
69 if (fheader.fail()) {
70 printErrorOpeningFile(m_className + "::Initialise", "header", header);
71 }
72
73 // Temporary variables for use in file reading
74 char* token = NULL;
75 bool readerror = false;
76 bool readstop = false;
77 int il = 0;
78
79 // Read the header to get the number of nodes and elements.
80 fheader.getline(line, size, '\n');
81 token = strtok(line, " ");
82 nNodes = ReadInteger(token, 0, readerror);
83 token = strtok(NULL, " ");
84 nElements = ReadInteger(token, 0, readerror);
85 std::cout << m_className << "::Initialise:\n Read " << nNodes
86 << " nodes and " << nElements
87 << " elements from file " << header << ".\n";
88 if (readerror) {
89 printErrorReadingFile(m_className + "::Initialise", header, il);
90 fheader.close();
91 return false;
92 }
93
94 // Close the header file.
95 fheader.close();
96
97 // Open the nodes list.
98 std::ifstream fnodes;
99 fnodes.open(nlist.c_str(), std::ios::in);
100 if (fnodes.fail()) {
101 printErrorOpeningFile(m_className + "::Initialise", "nodes", nlist);
102 }
103
104 // Check the value of the unit.
105 double funit;
106 if (strcmp(unit.c_str(), "mum") == 0 || strcmp(unit.c_str(), "micron") == 0 ||
107 strcmp(unit.c_str(), "micrometer") == 0) {
108 funit = 0.0001;
109 } else if (strcmp(unit.c_str(), "mm") == 0 ||
110 strcmp(unit.c_str(), "millimeter") == 0) {
111 funit = 0.1;
112 } else if (strcmp(unit.c_str(), "cm") == 0 ||
113 strcmp(unit.c_str(), "centimeter") == 0) {
114 funit = 1.0;
115 } else if (strcmp(unit.c_str(), "m") == 0 ||
116 strcmp(unit.c_str(), "meter") == 0) {
117 funit = 100.0;
118 } else {
119 std::cerr << m_className << "::Initialise:\n"
120 << " Unknown length unit " << unit << ".\n";
121 ok = false;
122 funit = 1.0;
123 }
124 if (m_debug) {
125 std::cout << m_className << "::Initialise:\n"
126 << " Unit scaling factor = " << funit << ".\n";
127 }
128
129 // Read the nodes from the file.
130 Node newNode;
131 newNode.w.clear();
132 for (il = 0; il < nNodes; il++) {
133
134 // Get a line from the nodes file.
135 fnodes.getline(line, size, '\n');
136
137 // Ignore the first two characters.
138 token = strtok(line, " ");
139 token = strtok(NULL, " ");
140
141 // Get the node coordinates.
142 token = strtok(NULL, " ");
143 double xnode = ReadDouble(token, -1, readerror);
144 token = strtok(NULL, " ");
145 double ynode = ReadDouble(token, -1, readerror);
146 token = strtok(NULL, " ");
147 double znode = ReadDouble(token, -1, readerror);
148 if (readerror) {
149 printErrorReadingFile(m_className + "::Initialise", nlist, il);
150 fnodes.close();
151 return false;
152 }
153
154 // Set up and create a new node.
155 newNode.x = xnode * funit;
156 newNode.y = ynode * funit;
157 newNode.z = znode * funit;
158 nodes.push_back(newNode);
159 }
160
161 // Close the nodes file.
162 fnodes.close();
163
164 // Open the potential file.
165 std::ifstream fvolt;
166 fvolt.open(volt.c_str(), std::ios::in);
167 if (fvolt.fail()) {
168 printErrorOpeningFile(m_className + "::Initialise", "potentials", volt);
169 }
170
171 // Reset the line counter.
172 il = 1;
173
174 // Read past the header.
175 while (!readstop && fvolt.getline(line, size, '\n')) {
176 token = strtok(line, " ");
177 if (strcmp(token, "Perm:") == 0) readstop = true;
178 il++;
179 }
180
181 // Should have stopped: if not, print error message.
182 if (!readstop) {
183 std::cerr << m_className << "::Initialise:\n Error reading past header "
184 << "of potentials file " << volt << ".\n";
185 fvolt.close();
186 return false;
187 }
188
189 // Read past the permutation map (number of lines = nNodes).
190 for (int tl = 0; tl < nNodes; tl++) {
191 fvolt.getline(line, size, '\n');
192 il++;
193 }
194
195 // Read the potentials.
196 for (int tl = 0; tl < nNodes; tl++) {
197 double v;
198 fvolt.getline(line, size, '\n');
199 token = strtok(line, " ");
200 v = ReadDouble(token, -1, readerror);
201 if (readerror) {
202 printErrorReadingFile(m_className + "::Initialise", volt, il);
203 fvolt.close();
204 return false;
205 }
206 // Place the voltage in its appropriate node.
207 nodes[tl].v = v;
208 }
209
210 // Close the potentials file.
211 fvolt.close();
212
213 // Open the materials file.
214 std::ifstream fmplist;
215 fmplist.open(mplist.c_str(), std::ios::in);
216 if (fmplist.fail()) {
217 printErrorOpeningFile(m_className + "::Initialise", "materials", mplist);
218 }
219
220 // Read the dielectric constants from the materials file.
221 fmplist.getline(line, size, '\n');
222 token = strtok(line, " ");
223 if (readerror) {
224 std::cerr << m_className << "::Initialise:\n "
225 << "Error reading number of materials from " << mplist << ".\n";
226 fmplist.close();
227 ok = false;
228 return false;
229 }
230 m_nMaterials = ReadInteger(token, 0, readerror);
231 materials.resize(m_nMaterials);
232 for (unsigned int i = 0; i < m_nMaterials; ++i) {
233 materials[i].ohm = -1;
234 materials[i].eps = -1;
235 materials[i].medium = NULL;
236 }
237 for (il = 2; il < ((int)m_nMaterials + 2); il++) {
238 fmplist.getline(line, size, '\n');
239 token = strtok(line, " ");
240 ReadInteger(token, -1, readerror);
241 token = strtok(NULL, " ");
242 double dc = ReadDouble(token, -1.0, readerror);
243 if (readerror) {
244 printErrorReadingFile(m_className + "::Initialise", mplist, il);
245 fmplist.close();
246 ok = false;
247 return false;
248 }
249 materials[il - 2].eps = dc;
250 std::cout << m_className << "::Initialise:\n Set material " << il - 2
251 << " of " << m_nMaterials << " to eps " << dc << ".\n";
252 }
253
254 // Close the materials file.
255 fmplist.close();
256
257 // Find the lowest epsilon, check for eps = 0, set default drift media.
258 double epsmin = -1.;
259 unsigned int iepsmin = 0;
260 for (unsigned int imat = 0; imat < m_nMaterials; ++imat) {
261 if (materials[imat].eps < 0) continue;
262 if (materials[imat].eps == 0) {
263 std::cerr << m_className << "::Initialise:\n Material " << imat
264 << " has been assigned a permittivity equal to zero in\n "
265 << mplist << ".\n";
266 ok = false;
267 } else if (epsmin < 0. || epsmin > materials[imat].eps) {
268 epsmin = materials[imat].eps;
269 iepsmin = imat;
270 }
271 }
272
273 if (epsmin < 0.) {
274 std::cerr << m_className << "::Initialise:\n"
275 << " No material with positive permittivity found \n"
276 << " in material list " << mplist << ".\n";
277 ok = false;
278 } else {
279 for (unsigned int imat = 0; imat < m_nMaterials; ++imat) {
280 if (imat == iepsmin) {
281 materials[imat].driftmedium = true;
282 } else {
283 materials[imat].driftmedium = false;
284 }
285 }
286 }
287
288 // Open the elements file.
289 std::ifstream felems;
290 felems.open(elist.c_str(), std::ios::in);
291 if (felems.fail()) {
292 printErrorOpeningFile(m_className + "::Initialise", "elements", elist);
293 }
294
295 // Read the elements and their material indices.
296 elements.clear();
297 int highestnode = 0;
298 Element newElement;
299 for (il = 0; il < nElements; il++) {
300
301 // Get a line
302 felems.getline(line, size, '\n');
303
304 // Split into tokens.
305 token = strtok(line, " ");
306 // Read the 2nd-order element
307 // Note: Ordering of Elmer elements can be described in the
308 // ElmerSolver manual (appendix D. at the time of this comment)
309 // If the order read below is compared to the shape functions used
310 // eg. in ElectricField, the order is wrong, but note at the
311 // end of this function the order of elements 5,6,7 will change to
312 // 7,5,6 when actually recorded in newElement.emap to correct for this
313 token = strtok(NULL, " ");
314 int imat = ReadInteger(token, -1, readerror) - 1;
315 token = strtok(NULL, " ");
316 token = strtok(NULL, " ");
317 int in0 = ReadInteger(token, -1, readerror);
318 token = strtok(NULL, " ");
319 int in1 = ReadInteger(token, -1, readerror);
320 token = strtok(NULL, " ");
321 int in2 = ReadInteger(token, -1, readerror);
322 token = strtok(NULL, " ");
323 int in3 = ReadInteger(token, -1, readerror);
324 token = strtok(NULL, " ");
325 int in4 = ReadInteger(token, -1, readerror);
326 token = strtok(NULL, " ");
327 int in5 = ReadInteger(token, -1, readerror);
328 token = strtok(NULL, " ");
329 int in6 = ReadInteger(token, -1, readerror);
330 token = strtok(NULL, " ");
331 int in7 = ReadInteger(token, -1, readerror);
332 token = strtok(NULL, " ");
333 int in8 = ReadInteger(token, -1, readerror);
334 token = strtok(NULL, " ");
335 int in9 = ReadInteger(token, -1, readerror);
336
337 if (m_debug && il < 10) {
338 std::cout << " Read nodes " << in0 << ", " << in1 << ", " << in2
339 << ", " << in3 << ", ... from element " << il + 1 << " of "
340 << nElements << " with mat " << imat << ".\n";
341 }
342
343 // Check synchronisation.
344 if (readerror) {
345 printErrorReadingFile(m_className + "::Initialise", elist, il);
346 felems.close();
347 ok = false;
348 return false;
349 }
350
351 // Check the material number and ensure that epsilon is non-negative.
352 if (imat < 0 || imat > (int)m_nMaterials) {
353 std::cerr << m_className << "::Initialise:\n";
354 std::cerr << " Out-of-range material number on file " << elist
355 << " (line " << il << ").\n";
356 std::cerr << " Element: " << il << ", material: " << imat << "\n";
357 std::cerr << " nodes: (" << in0 << ", " << in1 << ", " << in2 << ", "
358 << in3 << ", " << in4 << ", " << in5 << ", " << in6 << ", "
359 << in7 << ", " << in8 << ", " << in9 << ")\n";
360 ok = false;
361 }
362 if (materials[imat].eps < 0) {
363 std::cerr << m_className << "::Initialise:\n";
364 std::cerr << " Element " << il << " in element list " << elist << "\n";
365 std::cerr << " uses material " << imat
366 << " which has not been assigned\n";
367 std::cerr << " a positive permittivity in material list " << mplist
368 << ".\n";
369 ok = false;
370 }
371
372 // Check the node numbers.
373 if (in0 < 1 || in1 < 1 || in2 < 1 || in3 < 1 || in4 < 1 || in5 < 1 ||
374 in6 < 1 || in7 < 1 || in8 < 1 || in9 < 1) {
375 std::cerr << m_className << "::Initialise:\n";
376 std::cerr << " Found a node number < 1 on file " << elist << " (line "
377 << il << ").\n";
378 std::cerr << " Element: " << il << ", material: " << imat << "\n";
379 std::cerr << " nodes: (" << in0 << ", " << in1 << ", " << in2 << ", "
380 << in3 << ", " << in4 << ", " << in5 << ", " << in6 << ", "
381 << in7 << ", " << in8 << ", " << in9 << ")\n";
382 ok = false;
383 }
384 if (in0 > highestnode) highestnode = in0;
385 if (in1 > highestnode) highestnode = in1;
386 if (in2 > highestnode) highestnode = in2;
387 if (in3 > highestnode) highestnode = in3;
388 if (in4 > highestnode) highestnode = in4;
389 if (in5 > highestnode) highestnode = in5;
390 if (in6 > highestnode) highestnode = in6;
391 if (in7 > highestnode) highestnode = in7;
392 if (in8 > highestnode) highestnode = in8;
393 if (in9 > highestnode) highestnode = in9;
394
395 // These elements must not be degenerate.
396 if (in0 == in1 || in0 == in2 || in0 == in3 || in0 == in4 || in0 == in5 ||
397 in0 == in6 || in0 == in7 || in0 == in8 || in0 == in9 || in1 == in2 ||
398 in1 == in3 || in1 == in4 || in1 == in5 || in1 == in6 || in1 == in7 ||
399 in1 == in8 || in1 == in9 || in2 == in3 || in2 == in4 || in2 == in5 ||
400 in2 == in6 || in2 == in7 || in2 == in8 || in2 == in9 || in3 == in4 ||
401 in3 == in5 || in3 == in6 || in3 == in7 || in3 == in8 || in3 == in9 ||
402 in4 == in5 || in4 == in6 || in4 == in7 || in4 == in8 || in4 == in9 ||
403 in5 == in6 || in5 == in7 || in5 == in8 || in5 == in9 || in6 == in7 ||
404 in6 == in8 || in6 == in9 || in7 == in8 || in7 == in9 || in8 == in9) {
405 std::cerr << m_className << "::Initialise:\n Element " << il
406 << " of file " << elist << " is degenerate,\n"
407 << " no such elements are allowed in this type of map.\n";
408 ok = false;
409 }
410
411 newElement.degenerate = false;
412
413 // Store the material reference.
414 newElement.matmap = imat;
415
416 // Node references
417 newElement.emap[0] = in0 - 1;
418 newElement.emap[1] = in1 - 1;
419 newElement.emap[2] = in2 - 1;
420 newElement.emap[3] = in3 - 1;
421 newElement.emap[4] = in4 - 1;
422 newElement.emap[7] = in5 - 1;
423 newElement.emap[5] = in6 - 1;
424 newElement.emap[6] = in7 - 1;
425 newElement.emap[8] = in8 - 1;
426 newElement.emap[9] = in9 - 1;
427 elements.push_back(newElement);
428 }
429
430 // Close the elements file.
431 felems.close();
432
433 // Set the ready flag.
434 if (ok) {
435 m_ready = true;
436 } else {
437 std::cerr << m_className << "::Initialise:\n Field map could not be "
438 << "read and cannot be interpolated.\n";
439 return false;
440 }
441
442 std::cout << m_className << "::Initialise:\n Finished.\n";
443
444 // Remove weighting fields (if any).
445 wfields.clear();
446 wfieldsOk.clear();
447 nWeightingFields = 0;
448
449 // Establish the ranges.
450 SetRange();
451 UpdatePeriodicity();
452 return true;
453}
Garfield::ComponentElmer::UpdatePeriodicity
void UpdatePeriodicity()
Verify periodicities.
Definition: ComponentElmer.hh:65
Garfield::ComponentFieldMap::ReadDouble
double ReadDouble(char *token, double def, bool &error)
Definition: ComponentFieldMap.cc:2645
Garfield::ComponentFieldMap::ReadInteger
int ReadInteger(char *token, int def, bool &error)
Definition: ComponentFieldMap.cc:2635
Garfield::ComponentFieldMap::SetRange
virtual void SetRange()
Calculate x, y, z, V and angular ranges.
Definition: ComponentFieldMap.cc:2239
Garfield::ComponentFieldMap::wfields
std::vector< std::string > wfields
Definition: ComponentFieldMap.hh:136

Referenced by ComponentElmer().

◆ IsInBoundingBox()

virtual bool Garfield::ComponentElmer::IsInBoundingBox ( const double  x,
const double  y,
const double  z 
) const
inlinevirtual

Reimplemented from Garfield::ComponentFieldMap.

Definition at line 38 of file ComponentElmer.hh.

39 {
40 return x >= xMinBoundingBox && x <= xMaxBoundingBox &&
41 y >= yMinBoundingBox && y <= yMaxBoundingBox &&
42 z >= zMinBoundingBox && y <= zMaxBoundingBox;
43 }

◆ 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 455 of file ComponentElmer.cc.

455 {
456
457 if (!m_ready) {
458 PrintNotReady("SetWeightingField");
459 std::cerr << " Weighting field cannot be added.\n";
460 return false;
461 }
462
463 // Keep track of the success.
464 bool ok = true;
465
466 // Open the voltage list.
467 std::ifstream fwvolt;
468 fwvolt.open(wvolt.c_str(), std::ios::in);
469 if (fwvolt.fail()) {
470 printErrorOpeningFile(m_className + "::SetWeightingField",
471 "potential", wvolt);
472
473 return false;
474 }
475
476 // Check if a weighting field with the same label already exists.
477 int iw = nWeightingFields;
478 for (int i = nWeightingFields; i--;) {
479 if (wfields[i] == label) {
480 iw = i;
481 break;
482 }
483 }
484 if (iw == nWeightingFields) {
485 ++nWeightingFields;
486 wfields.resize(nWeightingFields);
487 wfieldsOk.resize(nWeightingFields);
488 for (int j = nNodes; j--;) {
489 nodes[j].w.resize(nWeightingFields);
490 }
491 } else {
492 std::cout << m_className << "::SetWeightingField:\n"
493 << " Replacing existing weighting field " << label << ".\n";
494 }
495 wfields[iw] = label;
496 wfieldsOk[iw] = false;
497
498 // Temporary variables for use in file reading
499 const int size = 100;
500 char line[size];
501 char* token = NULL;
502 bool readerror = false;
503 bool readstop = false;
504 int il = 1;
505
506 // Read past the header.
507 while (!readstop && fwvolt.getline(line, size, '\n')) {
508 token = strtok(line, " ");
509 if (strcmp(token, "Perm:") == 0) readstop = true;
510 il++;
511 }
512
513 // Should have stopped: if not, print error message.
514 if (!readstop) {
515 std::cerr << m_className << "::SetWeightingField:\n Error reading past "
516 << "header of potentials file " << wvolt << ".\n";
517 fwvolt.close();
518 ok = false;
519 return false;
520 }
521
522 // Read past the permutation map (number of lines = nNodes).
523 for (int tl = 0; tl < nNodes; tl++) {
524 fwvolt.getline(line, size, '\n');
525 il++;
526 }
527
528 // Read the potentials.
529 for (int tl = 0; tl < nNodes; tl++) {
530 double v;
531 fwvolt.getline(line, size, '\n');
532 token = strtok(line, " ");
533 v = ReadDouble(token, -1, readerror);
534 if (readerror) {
535 printErrorReadingFile(m_className + "::SetWeightingField", wvolt, il);
536 fwvolt.close();
537 ok = false;
538 return false;
539 }
540 // Place the weighting potential at its appropriate node and index.
541 nodes[tl].w[iw] = v;
542 }
543
544 // Close the potentials file.
545 fwvolt.close();
546 std::cout << m_className << "::SetWeightingField:\n"
547 << " Read potentials from file " << wvolt.c_str() << ".\n";
548
549 // Set the ready flag.
550 wfieldsOk[iw] = ok;
551 if (!ok) {
552 std::cerr << m_className << "::SetWeightingField:\n Field map could not "
553 << "be read and cannot be interpolated.\n";
554 return false;
555 }
556 return true;
557}

◆ UpdatePeriodicity()

void Garfield::ComponentElmer::UpdatePeriodicity ( )
inlineprotectedvirtual

Verify periodicities.

Implements Garfield::ComponentFieldMap.

Definition at line 65 of file ComponentElmer.hh.

65{ UpdatePeriodicityCommon(); }

Referenced by Initialise().

◆ WeightingField()

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

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

Implements Garfield::ComponentFieldMap.

Definition at line 666 of file ComponentElmer.cc.

668 {
669
670 // Initial values
671 wx = wy = wz = 0;
672
673 // Do not proceed if not properly initialised.
674 if (!m_ready) return;
675
676 // Look for the label.
677 int iw = 0;
678 bool found = false;
679 for (int i = nWeightingFields; i--;) {
680 if (wfields[i] == label) {
681 iw = i;
682 found = true;
683 break;
684 }
685 }
686
687 // Do not proceed if the requested weighting field does not exist.
688 if (!found) return;
689 // Check if the weighting field is properly initialised.
690 if (!wfieldsOk[iw]) return;
691
692 // Copy the coordinates.
693 double x = xin, y = yin, z = zin;
694
695 // Map the coordinates onto field map coordinates
696 bool xmirr, ymirr, zmirr;
697 double rcoordinate, rotation;
698 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
699
700 if (m_warning) PrintWarning("WeightingField");
701
702 // Find the element that contains this point.
703 double t1, t2, t3, t4, jac[4][4], det;
704 const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
705 // Check if the point is in the mesh.
706 if (imap < 0) return;
707
708 if (m_debug) {
709 PrintElement("WeightingField", x, y, z, t1, t2, t3, t4, imap, 10, iw);
710 }
711
712 const Element& element = elements[imap];
713 const Node& n0 = nodes[element.emap[0]];
714 const Node& n1 = nodes[element.emap[1]];
715 const Node& n2 = nodes[element.emap[2]];
716 const Node& n3 = nodes[element.emap[3]];
717 const Node& n4 = nodes[element.emap[4]];
718 const Node& n5 = nodes[element.emap[5]];
719 const Node& n6 = nodes[element.emap[6]];
720 const Node& n7 = nodes[element.emap[7]];
721 const Node& n8 = nodes[element.emap[8]];
722 const Node& n9 = nodes[element.emap[9]];
723 // Tetrahedral field
724 const double invdet = 1. / det;
725 wx = -(n0.w[iw] * (4 * t1 - 1) * jac[0][1] +
726 n1.w[iw] * (4 * t2 - 1) * jac[1][1] +
727 n2.w[iw] * (4 * t3 - 1) * jac[2][1] +
728 n3.w[iw] * (4 * t4 - 1) * jac[3][1] +
729 n4.w[iw] * (4 * t2 * jac[0][1] + 4 * t1 * jac[1][1]) +
730 n5.w[iw] * (4 * t3 * jac[0][1] + 4 * t1 * jac[2][1]) +
731 n6.w[iw] * (4 * t4 * jac[0][1] + 4 * t1 * jac[3][1]) +
732 n7.w[iw] * (4 * t3 * jac[1][1] + 4 * t2 * jac[2][1]) +
733 n8.w[iw] * (4 * t4 * jac[1][1] + 4 * t2 * jac[3][1]) +
734 n9.w[iw] * (4 * t4 * jac[2][1] + 4 * t3 * jac[3][1])) * invdet;
735 wy = -(n0.w[iw] * (4 * t1 - 1) * jac[0][2] +
736 n1.w[iw] * (4 * t2 - 1) * jac[1][2] +
737 n2.w[iw] * (4 * t3 - 1) * jac[2][2] +
738 n3.w[iw] * (4 * t4 - 1) * jac[3][2] +
739 n4.w[iw] * (4 * t2 * jac[0][2] + 4 * t1 * jac[1][2]) +
740 n5.w[iw] * (4 * t3 * jac[0][2] + 4 * t1 * jac[2][2]) +
741 n6.w[iw] * (4 * t4 * jac[0][2] + 4 * t1 * jac[3][2]) +
742 n7.w[iw] * (4 * t3 * jac[1][2] + 4 * t2 * jac[2][2]) +
743 n8.w[iw] * (4 * t4 * jac[1][2] + 4 * t2 * jac[3][2]) +
744 n9.w[iw] * (4 * t4 * jac[2][2] + 4 * t3 * jac[3][2])) * invdet;
745 wz = -(n0.w[iw] * (4 * t1 - 1) * jac[0][3] +
746 n1.w[iw] * (4 * t2 - 1) * jac[1][3] +
747 n2.w[iw] * (4 * t3 - 1) * jac[2][3] +
748 n3.w[iw] * (4 * t4 - 1) * jac[3][3] +
749 n4.w[iw] * (4 * t2 * jac[0][3] + 4 * t1 * jac[1][3]) +
750 n5.w[iw] * (4 * t3 * jac[0][3] + 4 * t1 * jac[2][3]) +
751 n6.w[iw] * (4 * t4 * jac[0][3] + 4 * t1 * jac[3][3]) +
752 n7.w[iw] * (4 * t3 * jac[1][3] + 4 * t2 * jac[2][3]) +
753 n8.w[iw] * (4 * t4 * jac[1][3] + 4 * t2 * jac[3][3]) +
754 n9.w[iw] * (4 * t4 * jac[2][3] + 4 * t3 * jac[3][3])) * invdet;
755
756 // Transform field to global coordinates
757 UnmapFields(wx, wy, wz, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
758}

◆ WeightingPotential()

double Garfield::ComponentElmer::WeightingPotential ( const double  x,
const double  y,
const double  z,
const std::string &  label 
)
virtual

Implements Garfield::ComponentFieldMap.

Definition at line 760 of file ComponentElmer.cc.

762 {
763
764 // Do not proceed if not properly initialised.
765 if (!m_ready) return 0.;
766
767 // Look for the label.
768 int iw = 0;
769 bool found = false;
770 for (int i = nWeightingFields; i--;) {
771 if (wfields[i] == label) {
772 iw = i;
773 found = true;
774 break;
775 }
776 }
777
778 // Do not proceed if the requested weighting field does not exist.
779 if (!found) return 0.;
780 // Check if the weighting field is properly initialised.
781 if (!wfieldsOk[iw]) return 0.;
782
783 // Copy the coordinates.
784 double x = xin, y = yin, z = zin;
785
786 // Map the coordinates onto field map coordinates.
787 bool xmirr, ymirr, zmirr;
788 double rcoordinate, rotation;
789 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
790
791 if (m_warning) PrintWarning("WeightingPotential");
792
793 // Find the element that contains this point.
794 double t1, t2, t3, t4, jac[4][4], det;
795 const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
796 if (imap < 0) return 0.;
797
798 if (m_debug) {
799 PrintElement("WeightingPotential", x, y, z, t1, t2, t3, t4, imap, 10, iw);
800 }
801
802 const Element& element = elements[imap];
803 const Node& n0 = nodes[element.emap[0]];
804 const Node& n1 = nodes[element.emap[1]];
805 const Node& n2 = nodes[element.emap[2]];
806 const Node& n3 = nodes[element.emap[3]];
807 const Node& n4 = nodes[element.emap[4]];
808 const Node& n5 = nodes[element.emap[5]];
809 const Node& n6 = nodes[element.emap[6]];
810 const Node& n7 = nodes[element.emap[7]];
811 const Node& n8 = nodes[element.emap[8]];
812 const Node& n9 = nodes[element.emap[9]];
813 // Tetrahedral field
814 return n0.w[iw] * t1 * (2 * t1 - 1) + n1.w[iw] * t2 * (2 * t2 - 1) +
815 n2.w[iw] * t3 * (2 * t3 - 1) + n3.w[iw] * t4 * (2 * t4 - 1) +
816 4 * n4.w[iw] * t1 * t2 + 4 * n5.w[iw] * t1 * t3 +
817 4 * n6.w[iw] * t1 * t4 + 4 * n7.w[iw] * t2 * t3 +
818 4 * n8.w[iw] * t2 * t4 + 4 * n9.w[iw] * t3 * t4;
819}
The documentation for this class was generated from the following files: