Garfield++ v1r0
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

#include <ComponentElmer.hh>

+ Inheritance diagram for Garfield::ComponentElmer:

Public Member Functions

 ComponentElmer ()
 
 ComponentElmer (std::string header, std::string elist, std::string nlist, std::string mplist, std::string volt, std::string unit)
 
 ~ComponentElmer ()
 
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)
 
bool IsInBoundingBox (const double x, const double y, const double z)
 
bool Initialise (std::string header="mesh.header", std::string elist="mesh.elements", std::string nlist="mesh.nodes", std::string mplist="dielectrics.dat", std::string volt="out.result", std::string unit="cm")
 
bool SetWeightingField (std::string prnsol, std::string label)
 
- Public Member Functions inherited from Garfield::ComponentFieldMap
 ComponentFieldMap ()
 
virtual ~ComponentFieldMap ()
 
virtual void SetRange ()
 
void PrintRange ()
 
virtual bool IsInBoundingBox (const double x, const double y, const double z)
 
virtual bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
 
bool GetVoltageRange (double &vmin, double &vmax)
 
void PrintMaterials ()
 
void DriftMedium (int imat)
 
void NotDriftMedium (int imat)
 
int GetNumberOfMaterials ()
 
double GetPermittivity (const int imat)
 
double GetConductivity (const int imat)
 
void SetMedium (const int imat, Medium *medium)
 
Medium * GetMedium (const unsigned int &i) const
 
Medium * GetMedium (const double &x, const double &y, const double &z)=0
 
int GetNumberOfMedia ()
 
int GetNumberOfElements () const
 
bool GetElement (const int i, double &vol, double &dmin, double &dmax)
 
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 ()
 
- Public Member Functions inherited from Garfield::ComponentBase
 ComponentBase ()
 
virtual ~ComponentBase ()
 
virtual void SetGeometry (GeometryBase *geo)
 
virtual void Clear ()
 
virtual Medium * GetMedium (const double &x, const double &y, const double &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
 
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)
 
virtual bool IsReady ()
 
virtual bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
 
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 (double x0, double y0, double z0, double &xw, double &yw, double &rw)
 
void EnablePeriodicityX ()
 
void DisablePeriodicityX ()
 
void EnablePeriodicityY ()
 
void DisablePeriodicityY ()
 
void EnablePeriodicityZ ()
 
void DisablePeriodicityZ ()
 
void EnableMirrorPeriodicityX ()
 
void DisableMirrorPeriodicityX ()
 
void EnableMirrorPeriodicityY ()
 
void DisableMirrorPeriodicityY ()
 
void EnableMirrorPeriodicityZ ()
 
void DisableMirrorPeriodicityZ ()
 
void EnableAxialPeriodicityX ()
 
void DisableAxialPeriodicityX ()
 
void EnableAxialPeriodicityY ()
 
void DisableAxialPeriodicityY ()
 
void EnableAxialPeriodicityZ ()
 
void DisableAxialPeriodicityZ ()
 
void EnableRotationSymmetryX ()
 
void DisableRotationSymmetryX ()
 
void EnableRotationSymmetryY ()
 
void DisableRotationSymmetryY ()
 
void EnableRotationSymmetryZ ()
 
void DisableRotationSymmetryZ ()
 
void EnableDebugging ()
 
void DisableDebugging ()
 

Protected Member Functions

void UpdatePeriodicity ()
 
double GetElementVolume (const int i)
 
void GetAspectRatio (const int i, double &dmin, double &dmax)
 
- Protected Member Functions inherited from Garfield::ComponentFieldMap
void Reset ()
 
virtual void UpdatePeriodicity ()=0
 
void UpdatePeriodicity2d ()
 
void UpdatePeriodicityCommon ()
 
int Coordinates3 (double x, double y, double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det, int imap)
 
int Coordinates4 (double x, double y, double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det, int imap)
 
int Coordinates5 (double x, double y, double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det, int imap)
 
int Coordinates12 (double x, double y, double z, double &t1, double &t2, double &t3, double &t4, int imap)
 
int Coordinates13 (double x, double y, double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det, int imap)
 
int CoordinatesCube (double x, double y, double z, double &t1, double &t2, double &t3, TMatrixD *&jac, std::vector< TMatrixD * > &dN, int imap)
 
void Jacobian3 (int i, double u, double v, double w, double &det, double jac[4][4])
 
void Jacobian5 (int i, double u, double v, double &det, double jac[4][4])
 
void Jacobian13 (int i, double t, double u, double v, double w, double &det, double jac[4][4])
 
void JacobianCube (int i, double t1, double t2, double t3, TMatrixD *&jac, std::vector< TMatrixD * > &dN)
 
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)
 
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)
 
int FindElementCube (const double x, const double y, const double z, double &t1, double &t2, double &t3, TMatrixD *&jac, std::vector< TMatrixD * > &dN)
 
void MapCoordinates (double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
 
void UnmapFields (double &ex, double &ey, double &ez, double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation)
 
int ReadInteger (char *token, int def, bool &error)
 
double ReadDouble (char *token, double def, bool &error)
 
virtual double GetElementVolume (const int i)=0
 
virtual void GetAspectRatio (const int i, double &dmin, double &dmax)=0
 
void CalculateElementBoundingBoxes (void)
 
virtual void Reset ()=0
 
virtual void UpdatePeriodicity ()=0
 

Additional Inherited Members

- Protected Attributes inherited from Garfield::ComponentFieldMap
bool is3d
 
int nElements
 
std::vector< element > elements
 
int lastElement
 
bool cacheElemBoundingBoxes
 
int nNodes
 
std::vector< node > nodes
 
int 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 deleteBackground
 
bool checkMultipleElement
 
bool warning
 
- Protected Attributes inherited from Garfield::ComponentBase
std::string m_className
 
GeometryBase * theGeometry
 
bool ready
 
bool xPeriodic
 
bool yPeriodic
 
bool zPeriodic
 
bool xMirrorPeriodic
 
bool yMirrorPeriodic
 
bool zMirrorPeriodic
 
bool xAxiallyPeriodic
 
bool yAxiallyPeriodic
 
bool zAxiallyPeriodic
 
bool xRotationSymmetry
 
bool yRotationSymmetry
 
bool zRotationSymmetry
 
double bx0
 
double by0
 
double bz0
 
bool debug
 

Detailed Description

Definition at line 10 of file ComponentElmer.hh.

Constructor & Destructor Documentation

◆ ComponentElmer() [1/2]

Garfield::ComponentElmer::ComponentElmer ( )

Definition at line 14 of file ComponentElmer.cc.

14 : ComponentFieldMap() {
15
16 m_className = "ComponentElmer";
17 ready = false;
18}

◆ ComponentElmer() [2/2]

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

Definition at line 20 of file ComponentElmer.cc.

23 : ComponentFieldMap() {
24
25 m_className = "ComponentElmer";
26 Initialise(header, elist, nlist, mplist, volt, unit);
27}
Garfield::ComponentElmer::Initialise
bool Initialise(std::string header="mesh.header", std::string elist="mesh.elements", std::string nlist="mesh.nodes", std::string mplist="dielectrics.dat", std::string volt="out.result", std::string unit="cm")
Definition: ComponentElmer.cc:29

◆ ~ComponentElmer()

Garfield::ComponentElmer::~ComponentElmer ( )
inline

Definition at line 18 of file ComponentElmer.hh.

18{}

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

582 {
583
584 // Copy the coordinates
585 double x = xin, y = yin, z = zin;
586
587 // Map the coordinates onto field map coordinates
588 bool xmirrored, ymirrored, zmirrored;
589 double rcoordinate, rotation;
590 MapCoordinates(x, y, z, xmirrored, ymirrored, zmirrored, rcoordinate,
591 rotation);
592
593 // Initial values
594 ex = ey = ez = volt = 0.;
595 status = 0;
596 m = 0;
597
598 // Do not proceed if not properly initialised.
599 if (!ready) {
600 status = -10;
601 std::cerr << m_className << "::ElectricField:\n";
602 std::cerr << " Field map not available for interpolation.\n";
603 return;
604 }
605
606 if (warning) {
607 std::cerr << m_className << "::ElectricField:\n";
608 std::cerr << " Warnings have been issued for this field map.\n";
609 }
610
611 // Find the element that contains this point
612 double t1, t2, t3, t4, jac[4][4], det;
613 int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
614 if (imap < 0) {
615 if (debug) {
616 std::cout << m_className << "::ElectricField:\n";
617 std::cout << " Point (" << x << ", " << y << ", " << z
618 << " not in the mesh.\n";
619 }
620 status = -6;
621 return;
622 }
623
624 if (debug) {
625 std::cout << m_className << "::ElectricField:\n";
626 std::cout << " Global: (" << x << ", " << y << ", " << z << "),\n";
627 std::cout << " Local: (" << t1 << ", " << t2 << ", " << t3 << ", " << t4
628 << " in element " << imap << "\n";
629 std::cout
630 << " Node x y z V\n";
631 for (int i = 0; i < 10; i++) {
632 printf(" %-5d %12g %12g %12g %12g\n", elements[imap].emap[i],
633 nodes[elements[imap].emap[i]].x, nodes[elements[imap].emap[i]].y,
634 nodes[elements[imap].emap[i]].z, nodes[elements[imap].emap[i]].v);
635 }
636 }
637
638 // Tetrahedral field
639 volt = nodes[elements[imap].emap[0]].v * t1 * (2 * t1 - 1) +
640 nodes[elements[imap].emap[1]].v * t2 * (2 * t2 - 1) +
641 nodes[elements[imap].emap[2]].v * t3 * (2 * t3 - 1) +
642 nodes[elements[imap].emap[3]].v * t4 * (2 * t4 - 1) +
643 4 * nodes[elements[imap].emap[4]].v * t1 * t2 +
644 4 * nodes[elements[imap].emap[5]].v * t1 * t3 +
645 4 * nodes[elements[imap].emap[6]].v * t1 * t4 +
646 4 * nodes[elements[imap].emap[7]].v * t2 * t3 +
647 4 * nodes[elements[imap].emap[8]].v * t2 * t4 +
648 4 * nodes[elements[imap].emap[9]].v * t3 * t4;
649 ex = -(nodes[elements[imap].emap[0]].v * (4 * t1 - 1) * jac[0][1] +
650 nodes[elements[imap].emap[1]].v * (4 * t2 - 1) * jac[1][1] +
651 nodes[elements[imap].emap[2]].v * (4 * t3 - 1) * jac[2][1] +
652 nodes[elements[imap].emap[3]].v * (4 * t4 - 1) * jac[3][1] +
653 nodes[elements[imap].emap[4]].v *
654 (4 * t2 * jac[0][1] + 4 * t1 * jac[1][1]) +
655 nodes[elements[imap].emap[5]].v *
656 (4 * t3 * jac[0][1] + 4 * t1 * jac[2][1]) +
657 nodes[elements[imap].emap[6]].v *
658 (4 * t4 * jac[0][1] + 4 * t1 * jac[3][1]) +
659 nodes[elements[imap].emap[7]].v *
660 (4 * t3 * jac[1][1] + 4 * t2 * jac[2][1]) +
661 nodes[elements[imap].emap[8]].v *
662 (4 * t4 * jac[1][1] + 4 * t2 * jac[3][1]) +
663 nodes[elements[imap].emap[9]].v *
664 (4 * t4 * jac[2][1] + 4 * t3 * jac[3][1])) /
665 det;
666 ey = -(nodes[elements[imap].emap[0]].v * (4 * t1 - 1) * jac[0][2] +
667 nodes[elements[imap].emap[1]].v * (4 * t2 - 1) * jac[1][2] +
668 nodes[elements[imap].emap[2]].v * (4 * t3 - 1) * jac[2][2] +
669 nodes[elements[imap].emap[3]].v * (4 * t4 - 1) * jac[3][2] +
670 nodes[elements[imap].emap[4]].v *
671 (4 * t2 * jac[0][2] + 4 * t1 * jac[1][2]) +
672 nodes[elements[imap].emap[5]].v *
673 (4 * t3 * jac[0][2] + 4 * t1 * jac[2][2]) +
674 nodes[elements[imap].emap[6]].v *
675 (4 * t4 * jac[0][2] + 4 * t1 * jac[3][2]) +
676 nodes[elements[imap].emap[7]].v *
677 (4 * t3 * jac[1][2] + 4 * t2 * jac[2][2]) +
678 nodes[elements[imap].emap[8]].v *
679 (4 * t4 * jac[1][2] + 4 * t2 * jac[3][2]) +
680 nodes[elements[imap].emap[9]].v *
681 (4 * t4 * jac[2][2] + 4 * t3 * jac[3][2])) /
682 det;
683 ez = -(nodes[elements[imap].emap[0]].v * (4 * t1 - 1) * jac[0][3] +
684 nodes[elements[imap].emap[1]].v * (4 * t2 - 1) * jac[1][3] +
685 nodes[elements[imap].emap[2]].v * (4 * t3 - 1) * jac[2][3] +
686 nodes[elements[imap].emap[3]].v * (4 * t4 - 1) * jac[3][3] +
687 nodes[elements[imap].emap[4]].v *
688 (4 * t2 * jac[0][3] + 4 * t1 * jac[1][3]) +
689 nodes[elements[imap].emap[5]].v *
690 (4 * t3 * jac[0][3] + 4 * t1 * jac[2][3]) +
691 nodes[elements[imap].emap[6]].v *
692 (4 * t4 * jac[0][3] + 4 * t1 * jac[3][3]) +
693 nodes[elements[imap].emap[7]].v *
694 (4 * t3 * jac[1][3] + 4 * t2 * jac[2][3]) +
695 nodes[elements[imap].emap[8]].v *
696 (4 * t4 * jac[1][3] + 4 * t2 * jac[3][3]) +
697 nodes[elements[imap].emap[9]].v *
698 (4 * t4 * jac[2][3] + 4 * t3 * jac[3][3])) /
699 det;
700
701 // Transform field to global coordinates
702 UnmapFields(ex, ey, ez, x, y, z, xmirrored, ymirrored, zmirrored, rcoordinate,
703 rotation);
704
705 // Drift medium?
706 if (debug) {
707 std::cout << m_className << "::ElectricField:\n";
708 std::cout << " Material " << elements[imap].matmap << ", drift flag "
709 << materials[elements[imap].matmap].driftmedium << "\n";
710 }
711 m = materials[elements[imap].matmap].medium;
712 status = -5;
713 if (materials[elements[imap].matmap].driftmedium) {
714 if (m != 0) {
715 if (m->IsDriftable()) status = 0;
716 }
717 }
718}
Garfield::ComponentFieldMap::MapCoordinates
void MapCoordinates(double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
Definition: ComponentFieldMap.cc:4036
Garfield::ComponentFieldMap::materials
std::vector< material > materials
Definition: ComponentFieldMap.hh:124
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)
Definition: ComponentFieldMap.cc:4154
Garfield::ComponentFieldMap::elements
std::vector< element > elements
Definition: ComponentFieldMap.hh:96
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)
Definition: ComponentFieldMap.cc:356

◆ 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

Implements Garfield::ComponentFieldMap.

Definition at line 571 of file ComponentElmer.cc.

573 {
574
575 double v = 0.;
576 ElectricField(x, y, z, ex, ey, ez, v, m, status);
577}
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:571

Referenced by ElectricField().

◆ GetAspectRatio()

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

Implements Garfield::ComponentFieldMap.

Definition at line 1003 of file ComponentElmer.cc.

1003 {
1004
1005 if (i < 0 || i >= nElements) {
1006 dmin = dmax = 0.;
1007 return;
1008 }
1009
1010 const int np = 4;
1011 // Loop over all pairs of vertices.
1012 for (int j = 0; j < np - 1; ++j) {
1013 for (int k = j + 1; k < np; ++k) {
1014 // Compute distance.
1015 const double dist = sqrt(
1016 pow(nodes[elements[i].emap[j]].x - nodes[elements[i].emap[k]].x, 2) +
1017 pow(nodes[elements[i].emap[j]].y - nodes[elements[i].emap[k]].y, 2) +
1018 pow(nodes[elements[i].emap[j]].z - nodes[elements[i].emap[k]].z, 2));
1019 if (k == 1) {
1020 dmin = dist;
1021 dmax = dist;
1022 } else {
1023 if (dist < dmin) dmin = dist;
1024 if (dist > dmax) dmax = dist;
1025 }
1026 }
1027 }
1028}
pow
DoubleAc pow(const DoubleAc &f, double p)
Definition: DoubleAc.cpp:336
sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:313

◆ GetElementVolume()

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

Implements Garfield::ComponentFieldMap.

Definition at line 971 of file ComponentElmer.cc.

971 {
972
973 if (i < 0 || i >= nElements) return 0.;
974
975 // Uses formula V = |a (dot) b x c|/6
976 // with a => "3", b => "1", c => "2" and origin "0"
977 const double vol =
978 fabs((nodes[elements[i].emap[3]].x - nodes[elements[i].emap[0]].x) *
979 ((nodes[elements[i].emap[1]].y - nodes[elements[i].emap[0]].y) *
980 (nodes[elements[i].emap[2]].z -
981 nodes[elements[i].emap[0]].z) -
982 (nodes[elements[i].emap[2]].y - nodes[elements[i].emap[0]].y) *
983 (nodes[elements[i].emap[1]].z -
984 nodes[elements[i].emap[0]].z)) +
985 (nodes[elements[i].emap[3]].y - nodes[elements[i].emap[0]].y) *
986 ((nodes[elements[i].emap[1]].z - nodes[elements[i].emap[0]].z) *
987 (nodes[elements[i].emap[2]].x -
988 nodes[elements[i].emap[0]].x) -
989 (nodes[elements[i].emap[2]].z - nodes[elements[i].emap[0]].z) *
990 (nodes[elements[i].emap[1]].x -
991 nodes[elements[i].emap[0]].x)) +
992 (nodes[elements[i].emap[3]].z - nodes[elements[i].emap[0]].z) *
993 ((nodes[elements[i].emap[1]].x - nodes[elements[i].emap[0]].x) *
994 (nodes[elements[i].emap[2]].y -
995 nodes[elements[i].emap[0]].y) -
996 (nodes[elements[i].emap[3]].x - nodes[elements[i].emap[0]].x) *
997 (nodes[elements[i].emap[1]].y -
998 nodes[elements[i].emap[0]].y))) /
999 6.;
1000 return vol;
1001}
fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:616

◆ GetMedium()

Medium * Garfield::ComponentElmer::GetMedium ( const double &  x,
const double &  y,
const double &  z 
)
virtual

Implements Garfield::ComponentFieldMap.

Definition at line 910 of file ComponentElmer.cc.

911 {
912
913 // Copy the coordinates
914 double x = xin, y = yin, z = zin;
915
916 // Map the coordinates onto field map coordinates
917 bool xmirrored, ymirrored, zmirrored;
918 double rcoordinate, rotation;
919 MapCoordinates(x, y, z, xmirrored, ymirrored, zmirrored, rcoordinate,
920 rotation);
921
922 // Do not proceed if not properly initialised.
923 if (!ready) {
924 std::cerr << m_className << "::GetMedium:\n";
925 std::cerr << " Field map not available for interpolation.\n";
926 return NULL;
927 }
928 if (warning) {
929 std::cerr << m_className << "::GetMedium:\n";
930 std::cerr << " Warnings have been issued for this field map.\n";
931 }
932
933 // Find the element that contains this point
934 double t1, t2, t3, t4, jac[4][4], det;
935 int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
936 if (imap < 0) {
937 if (debug) {
938 std::cout << m_className << "::GetMedium:\n";
939 std::cout << " Point (" << x << ", " << y << ", " << z
940 << ") not in the mesh.\n";
941 }
942 return NULL;
943 }
944 if (elements[imap].matmap < 0 || elements[imap].matmap >= nMaterials) {
945 if (debug) {
946 std::cerr << m_className << "::GetMedium:\n";
947 std::cerr << " Point (" << x << ", " << y
948 << ") has out of range material number " << imap << ".\n";
949 }
950 return NULL;
951 }
952
953 if (debug) {
954 std::cout << m_className << "::GetMedium:\n";
955 std::cout << " Global: (" << x << ", " << y << ", " << z << "),\n";
956 std::cout << " Local: (" << t1 << ", " << t2 << ", " << t3 << ", " << t4
957 << " in element " << imap << "\n";
958 std::cout
959 << " Node x y z V\n";
960 for (int ii = 0; ii < 10; ++ii) {
961 printf(" %-5d %12g %12g %12g %12g\n", elements[imap].emap[ii],
962 nodes[elements[imap].emap[ii]].x, nodes[elements[imap].emap[ii]].y,
963 nodes[elements[imap].emap[ii]].z,
964 nodes[elements[imap].emap[ii]].v);
965 }
966 }
967
968 return materials[elements[imap].matmap].medium;
969}

◆ Initialise()

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

Definition at line 29 of file ComponentElmer.cc.

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

Referenced by ComponentElmer().

◆ IsInBoundingBox()

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

Reimplemented from Garfield::ComponentFieldMap.

Definition at line 35 of file ComponentElmer.hh.

35 {
36 return x >= xMinBoundingBox && x <= xMaxBoundingBox &&
37 y >= yMinBoundingBox && y <= yMaxBoundingBox &&
38 z >= zMinBoundingBox && y <= zMaxBoundingBox;
39 }

◆ SetWeightingField()

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

Definition at line 461 of file ComponentElmer.cc.

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

◆ UpdatePeriodicity()

void Garfield::ComponentElmer::UpdatePeriodicity ( )
inlineprotectedvirtual

Implements Garfield::ComponentFieldMap.

Definition at line 51 of file ComponentElmer.hh.

51{ 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

Implements Garfield::ComponentFieldMap.

Definition at line 720 of file ComponentElmer.cc.

722 {
723
724 // Initial values
725 wx = wy = wz = 0;
726
727 // Do not proceed if not properly initialised.
728 if (!ready) return;
729
730 // Look for the label.
731 int iw = 0;
732 bool found = false;
733 for (int i = nWeightingFields; i--;) {
734 if (wfields[i] == label) {
735 iw = i;
736 found = true;
737 break;
738 }
739 }
740
741 // Do not proceed if the requested weighting field does not exist.
742 if (!found) return;
743 // Check if the weighting field is properly initialised.
744 if (!wfieldsOk[iw]) return;
745
746 // Copy the coordinates.
747 double x = xin, y = yin, z = zin;
748
749 // Map the coordinates onto field map coordinates
750 bool xmirrored, ymirrored, zmirrored;
751 double rcoordinate, rotation;
752 MapCoordinates(x, y, z, xmirrored, ymirrored, zmirrored, rcoordinate,
753 rotation);
754
755 if (warning) {
756 std::cerr << m_className << "::WeightingField:\n";
757 std::cerr << " Warnings have been issued for this field map.\n";
758 }
759
760 // Find the element that contains this point.
761 double t1, t2, t3, t4, jac[4][4], det;
762 int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
763 // Check if the point is in the mesh.
764 if (imap < 0) return;
765
766 if (debug) {
767 std::cout << m_className << "::WeightingField:\n";
768 std::cout << " Global: (" << x << ", " << y << ", " << z << "),\n";
769 std::cout << " Local: (" << t1 << ", " << t2 << ", " << t3 << ", " << t4
770 << " in element " << imap << "\n";
771 std::cout
772 << " Node x y z V\n";
773 for (int i = 0; i < 10; i++) {
774 printf(" %-5d %12g %12g %12g %12g\n", elements[imap].emap[i],
775 nodes[elements[imap].emap[i]].x, nodes[elements[imap].emap[i]].y,
776 nodes[elements[imap].emap[i]].z,
777 nodes[elements[imap].emap[i]].w[iw]);
778 }
779 }
780
781 // Tetrahedral field
782 wx = -(nodes[elements[imap].emap[0]].w[iw] * (4 * t1 - 1) * jac[0][1] +
783 nodes[elements[imap].emap[1]].w[iw] * (4 * t2 - 1) * jac[1][1] +
784 nodes[elements[imap].emap[2]].w[iw] * (4 * t3 - 1) * jac[2][1] +
785 nodes[elements[imap].emap[3]].w[iw] * (4 * t4 - 1) * jac[3][1] +
786 nodes[elements[imap].emap[4]].w[iw] *
787 (4 * t2 * jac[0][1] + 4 * t1 * jac[1][1]) +
788 nodes[elements[imap].emap[5]].w[iw] *
789 (4 * t3 * jac[0][1] + 4 * t1 * jac[2][1]) +
790 nodes[elements[imap].emap[6]].w[iw] *
791 (4 * t4 * jac[0][1] + 4 * t1 * jac[3][1]) +
792 nodes[elements[imap].emap[7]].w[iw] *
793 (4 * t3 * jac[1][1] + 4 * t2 * jac[2][1]) +
794 nodes[elements[imap].emap[8]].w[iw] *
795 (4 * t4 * jac[1][1] + 4 * t2 * jac[3][1]) +
796 nodes[elements[imap].emap[9]].w[iw] *
797 (4 * t4 * jac[2][1] + 4 * t3 * jac[3][1])) /
798 det;
799
800 wy = -(nodes[elements[imap].emap[0]].w[iw] * (4 * t1 - 1) * jac[0][2] +
801 nodes[elements[imap].emap[1]].w[iw] * (4 * t2 - 1) * jac[1][2] +
802 nodes[elements[imap].emap[2]].w[iw] * (4 * t3 - 1) * jac[2][2] +
803 nodes[elements[imap].emap[3]].w[iw] * (4 * t4 - 1) * jac[3][2] +
804 nodes[elements[imap].emap[4]].w[iw] *
805 (4 * t2 * jac[0][2] + 4 * t1 * jac[1][2]) +
806 nodes[elements[imap].emap[5]].w[iw] *
807 (4 * t3 * jac[0][2] + 4 * t1 * jac[2][2]) +
808 nodes[elements[imap].emap[6]].w[iw] *
809 (4 * t4 * jac[0][2] + 4 * t1 * jac[3][2]) +
810 nodes[elements[imap].emap[7]].w[iw] *
811 (4 * t3 * jac[1][2] + 4 * t2 * jac[2][2]) +
812 nodes[elements[imap].emap[8]].w[iw] *
813 (4 * t4 * jac[1][2] + 4 * t2 * jac[3][2]) +
814 nodes[elements[imap].emap[9]].w[iw] *
815 (4 * t4 * jac[2][2] + 4 * t3 * jac[3][2])) /
816 det;
817
818 wz = -(nodes[elements[imap].emap[0]].w[iw] * (4 * t1 - 1) * jac[0][3] +
819 nodes[elements[imap].emap[1]].w[iw] * (4 * t2 - 1) * jac[1][3] +
820 nodes[elements[imap].emap[2]].w[iw] * (4 * t3 - 1) * jac[2][3] +
821 nodes[elements[imap].emap[3]].w[iw] * (4 * t4 - 1) * jac[3][3] +
822 nodes[elements[imap].emap[4]].w[iw] *
823 (4 * t2 * jac[0][3] + 4 * t1 * jac[1][3]) +
824 nodes[elements[imap].emap[5]].w[iw] *
825 (4 * t3 * jac[0][3] + 4 * t1 * jac[2][3]) +
826 nodes[elements[imap].emap[6]].w[iw] *
827 (4 * t4 * jac[0][3] + 4 * t1 * jac[3][3]) +
828 nodes[elements[imap].emap[7]].w[iw] *
829 (4 * t3 * jac[1][3] + 4 * t2 * jac[2][3]) +
830 nodes[elements[imap].emap[8]].w[iw] *
831 (4 * t4 * jac[1][3] + 4 * t2 * jac[3][3]) +
832 nodes[elements[imap].emap[9]].w[iw] *
833 (4 * t4 * jac[2][3] + 4 * t3 * jac[3][3])) /
834 det;
835
836 // Transform field to global coordinates
837 UnmapFields(wx, wy, wz, x, y, z, xmirrored, ymirrored, zmirrored, rcoordinate,
838 rotation);
839}

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

843 {
844
845 // Do not proceed if not properly initialised.
846 if (!ready) return 0.;
847
848 // Look for the label.
849 int iw = 0;
850 bool found = false;
851 for (int i = nWeightingFields; i--;) {
852 if (wfields[i] == label) {
853 iw = i;
854 found = true;
855 break;
856 }
857 }
858
859 // Do not proceed if the requested weighting field does not exist.
860 if (!found) return 0.;
861 // Check if the weighting field is properly initialised.
862 if (!wfieldsOk[iw]) return 0.;
863
864 // Copy the coordinates.
865 double x = xin, y = yin, z = zin;
866
867 // Map the coordinates onto field map coordinates.
868 bool xmirrored, ymirrored, zmirrored;
869 double rcoordinate, rotation;
870 MapCoordinates(x, y, z, xmirrored, ymirrored, zmirrored, rcoordinate,
871 rotation);
872
873 if (warning) {
874 std::cerr << m_className << "::WeightingPotential:\n";
875 std::cerr << " Warnings have been issued for this field map.\n";
876 }
877
878 // Find the element that contains this point.
879 double t1, t2, t3, t4, jac[4][4], det;
880 int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
881 if (imap < 0) return 0.;
882
883 if (debug) {
884 std::cout << m_className << "::WeightingPotential:\n";
885 std::cout << " Global: (" << x << ", " << y << ", " << z << "),\n";
886 std::cout << " Local: (" << t1 << ", " << t2 << ", " << t3 << ", " << t4
887 << " in element " << imap << "\n";
888 std::cout
889 << " Node x y z V\n";
890 for (int i = 0; i < 10; i++) {
891 printf(" %-5d %12g %12g %12g %12g\n", elements[imap].emap[i],
892 nodes[elements[imap].emap[i]].x, nodes[elements[imap].emap[i]].y,
893 nodes[elements[imap].emap[i]].z, nodes[elements[imap].emap[i]].v);
894 }
895 }
896
897 // Tetrahedral field
898 return nodes[elements[imap].emap[0]].w[iw] * t1 * (2 * t1 - 1) +
899 nodes[elements[imap].emap[1]].w[iw] * t2 * (2 * t2 - 1) +
900 nodes[elements[imap].emap[2]].w[iw] * t3 * (2 * t3 - 1) +
901 nodes[elements[imap].emap[3]].w[iw] * t4 * (2 * t4 - 1) +
902 4 * nodes[elements[imap].emap[4]].w[iw] * t1 * t2 +
903 4 * nodes[elements[imap].emap[5]].w[iw] * t1 * t3 +
904 4 * nodes[elements[imap].emap[6]].w[iw] * t1 * t4 +
905 4 * nodes[elements[imap].emap[7]].w[iw] * t2 * t3 +
906 4 * nodes[elements[imap].emap[8]].w[iw] * t2 * t4 +
907 4 * nodes[elements[imap].emap[9]].w[iw] * t3 * t4;
908}
The documentation for this class was generated from the following files: