Garfield++ 5.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Classes | Public Member Functions | Protected Types | Protected Member Functions | Static Protected Member Functions | Protected Attributes | Friends | List of all members
Garfield::ComponentFieldMap Class Reference

Base class for components based on finite-element field maps. More...

#include <ComponentFieldMap.hh>

+ Inheritance diagram for Garfield::ComponentFieldMap:

Classes

struct  Element
 
struct  Material
 
struct  Node
 
struct  WeightingFieldCopy
 

Public Member Functions

 ComponentFieldMap ()=delete
 Default constructor.
 
 ComponentFieldMap (const std::string &name)
 Constructor.
 
virtual ~ComponentFieldMap ()
 Destructor.
 
bool Check ()
 Check element aspect ratio.
 
void PrintRange ()
 Show x, y, z, V and angular ranges.
 
void PrintMaterials ()
 List all currently defined materials.
 
void DriftMedium (const size_t imat)
 Flag a field map material as a drift medium.
 
void NotDriftMedium (const size_t imat)
 Flag a field map materials as a non-drift medium.
 
size_t GetNumberOfMaterials () const
 Return the number of materials in the field map.
 
double GetPermittivity (const size_t imat) const
 Return the relative permittivity of a field map material.
 
double GetConductivity (const size_t imat) const
 Return the conductivity of a field map material.
 
void SetMedium (const size_t imat, Medium *medium)
 Associate a field map material with a Medium object.
 
Medium * GetMedium (const size_t imat) const
 Return the Medium associated to a field map material.
 
void SetGas (Medium *medium)
 
virtual size_t GetNumberOfElements () const
 Return the number of mesh elements.
 
bool GetElement (const size_t i, double &vol, double &dmin, double &dmax) const
 Return the volume and aspect ratio of a mesh element.
 
virtual bool GetElement (const size_t i, size_t &mat, bool &drift, std::vector< size_t > &nodes) const
 Return the material and node indices of a mesh element.
 
virtual size_t GetNumberOfNodes () const
 
virtual bool GetNode (const size_t i, double &x, double &y, double &z) const
 
double GetPotential (const size_t i) const
 
void EnableCheckMapIndices (const bool on=true)
 
void EnableDeleteBackgroundElements (const bool on=true)
 Option to eliminate mesh elements in conductors (default: on).
 
void EnableTetrahedralTreeForElementSearch (const bool on=true)
 
void EnableConvergenceWarnings (const bool on=true)
 
Medium * GetMedium (const double x, const double y, const double z) override
 Get the medium at a given location (x, y, z).
 
void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status) override
 
void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status) override
 Calculate the drift field [V/cm] and potential [V] at (x, y, z).
 
void WeightingField (const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label) override
 
double WeightingPotential (const double x, const double y, const double z, const std::string &label) override
 
double DelayedWeightingPotential (double x, double y, double z, const double t, const std::string &label) override
 
bool IsInBoundingBox (const double x, const double y, const double z) const
 
bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax) override
 Get the bounding box coordinates.
 
bool GetElementaryCell (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax) override
 Get the coordinates of the elementary cell.
 
bool GetVoltageRange (double &vmin, double &vmax) override
 Calculate the voltage range [V].
 
void CopyWeightingPotential (const std::string &label, const std::string &labelSource, const double x, const double y, const double z, const double alpha, const double beta, const double gamma)
 
std::array< double, 3 > ElectricField (const double x, const double y, const double z)
 Calculate the drift field [V/cm] at (x, y, z).
 
- Public Member Functions inherited from Garfield::Component
 Component ()=delete
 Default constructor.
 
 Component (const std::string &name)
 Constructor.
 
virtual ~Component ()
 Destructor.
 
virtual void SetGeometry (Geometry *geo)
 Define the geometry.
 
virtual void Clear ()
 Reset.
 
std::array< double, 3 > ElectricField (const double x, const double y, const double z)
 Calculate the drift field [V/cm] at (x, y, z).
 
virtual double ElectricPotential (const double x, const double y, const double z)
 Calculate the (drift) electrostatic potential [V] at (x, y, z).
 
virtual void DelayedWeightingField (const double x, const double y, const double z, const double t, double &wx, double &wy, double &wz, const std::string &label)
 
virtual void MagneticField (const double x, const double y, const double z, double &bx, double &by, double &bz, int &status)
 
void SetMagneticField (const double bx, const double by, const double bz)
 Set a constant magnetic field.
 
virtual bool IsReady ()
 Ready for use?
 
double IntegrateFluxCircle (const double xc, const double yc, const double r, const unsigned int nI=50)
 
double IntegrateFluxSphere (const double xc, const double yc, const double zc, const double r, const unsigned int nI=20)
 
double IntegrateFluxParallelogram (const double x0, const double y0, const double z0, const double dx1, const double dy1, const double dz1, const double dx2, const double dy2, const double dz2, const unsigned int nU=20, const unsigned int nV=20)
 
double IntegrateWeightingFluxParallelogram (const std::string &label, const double x0, const double y0, const double z0, const double dx1, const double dy1, const double dz1, const double dx2, const double dy2, const double dz2, const unsigned int nU=20, const unsigned int nV=20)
 
double IntegrateFluxLine (const double x0, const double y0, const double z0, const double x1, const double y1, const double z1, const double xp, const double yp, const double zp, const unsigned int nI, const int isign=0)
 
virtual bool CrossedWire (const double x0, const double y0, const double z0, const double x1, const double y1, const double z1, double &xc, double &yc, double &zc, const bool centre, double &rc)
 
virtual bool InTrapRadius (const double q0, const double x0, const double y0, const double z0, double &xw, double &yw, double &rw)
 
virtual bool CrossedPlane (const double x0, const double y0, const double z0, const double x1, const double y1, const double z1, double &xc, double &yc, double &zc)
 
void EnablePeriodicityX (const bool on=true)
 Enable simple periodicity in the $x$ direction.
 
void EnablePeriodicityY (const bool on=true)
 Enable simple periodicity in the $y$ direction.
 
void EnablePeriodicityZ (const bool on=true)
 Enable simple periodicity in the $z$ direction.
 
void IsPeriodic (bool &perx, bool &pery, bool &perz)
 Return periodicity flags.
 
void EnableMirrorPeriodicityX (const bool on=true)
 Enable mirror periodicity in the $x$ direction.
 
void EnableMirrorPeriodicityY (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void EnableMirrorPeriodicityZ (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void IsMirrorPeriodic (bool &perx, bool &pery, bool &perz)
 Return mirror periodicity flags.
 
void EnableAxialPeriodicityX (const bool on=true)
 Enable axial periodicity in the $x$ direction.
 
void EnableAxialPeriodicityY (const bool on=true)
 Enable axial periodicity in the $y$ direction.
 
void EnableAxialPeriodicityZ (const bool on=true)
 Enable axial periodicity in the $z$ direction.
 
void IsAxiallyPeriodic (bool &perx, bool &pery, bool &perz)
 Return axial periodicity flags.
 
void EnableRotationSymmetryX (const bool on=true)
 Enable rotation symmetry around the $x$ axis.
 
void EnableRotationSymmetryY (const bool on=true)
 Enable rotation symmetry around the $y$ axis.
 
void EnableRotationSymmetryZ (const bool on=true)
 Enable rotation symmetry around the $z$ axis.
 
void IsRotationSymmetric (bool &rotx, bool &roty, bool &rotz)
 Return rotation symmetry flags.
 
void EnableDebugging ()
 Switch on debugging messages.
 
void DisableDebugging ()
 Switch off debugging messages.
 
virtual bool HasMagneticField () const
 Does the component have a non-zero magnetic field?
 
virtual bool HasTownsendMap () const
 Does the component have maps of the Townsend coefficient?
 
virtual bool HasAttachmentMap () const
 Does the component have attachment maps?
 
virtual bool HasVelocityMap () const
 Does the component have velocity maps?
 
virtual bool ElectronAttachment (const double, const double, const double, double &eta)
 Get the electron attachment coefficient.
 
virtual bool HoleAttachment (const double, const double, const double, double &eta)
 Get the hole attachment coefficient.
 
virtual bool ElectronTownsend (const double, const double, const double, double &alpha)
 Get the electron Townsend coefficient.
 
virtual bool HoleTownsend (const double, const double, const double, double &alpha)
 Get the hole Townsend coefficient.
 
virtual bool ElectronVelocity (const double, const double, const double, double &vx, double &vy, double &vz)
 Get the electron drift velocity.
 
virtual bool HoleVelocity (const double, const double, const double, double &vx, double &vy, double &vz)
 Get the hole drift velocity.
 
virtual bool GetElectronLifetime (const double, const double, const double, double &etau)
 
virtual bool GetHoleLifetime (const double, const double, const double, double &htau)
 
virtual double StepSizeHint ()
 

Protected Types

enum class  ElementType { Unknown = 0 , Serendipity = 5 , CurvedTetrahedron = 13 }
 

Protected Member Functions

void Reset () override
 Reset the component.
 
void Prepare ()
 
virtual void SetRange ()
 
void UpdatePeriodicity () override
 Verify periodicities.
 
void UpdatePeriodicity2d ()
 
void UpdatePeriodicityCommon ()
 
bool SetDefaultDriftMedium ()
 Find lowest epsilon, check for eps = 0, set default drift media flags.
 
int Field (const double x, const double y, const double z, double &fx, double &fy, double &fz, int &iel, const std::vector< double > &potentials) const
 Compute the electric/weighting field.
 
double Potential (const double x, const double y, const double z, const std::vector< double > &potentials) const
 Compute the electrostatic/weighting potential.
 
int FindElement5 (const double x, const double y, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det) const
 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) const
 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) const
 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, const double xpos, const double ypos, const double zpos, const bool xmirrored, const bool ymirrored, const bool zmirrored, const double rcoordinate, const double rotation) const
 Move (ex, ey, ez) to global coordinates.
 
virtual double GetElementVolume (const size_t i) const
 
virtual void GetAspectRatio (const size_t i, double &dmin, double &dmax) const
 
void PrintWarning (const std::string &header)
 
void PrintNotReady (const std::string &header) const
 
void PrintCouldNotOpen (const std::string &header, const std::string &filename) const
 
void PrintElement (const std::string &header, const double x, const double y, const double z, const double t1, const double t2, const double t3, const double t4, const size_t i, const std::vector< double > &potential) const
 
void TimeInterpolation (const double t, double &f0, double &f1, int &i0, int &i1)
 Interpolation of potential between two time slices.
 

Static Protected Member Functions

static double ScalingFactor (std::string unit)
 
static double Potential3 (const std::array< double, 6 > &v, const std::array< double, 3 > &t)
 Interpolate the potential in a triangle.
 
static void Field3 (const std::array< double, 6 > &v, const std::array< double, 3 > &t, double jac[4][4], const double det, double &ex, double &ey)
 Interpolate the field in a triangle.
 
static double Potential5 (const std::array< double, 8 > &v, const std::array< double, 2 > &t)
 Interpolate the potential in a curved quadrilateral.
 
static void Field5 (const std::array< double, 8 > &v, const std::array< double, 2 > &t, double jac[4][4], const double det, double &ex, double &ey)
 Interpolate the field in a curved quadrilateral.
 
static double Potential13 (const std::array< double, 10 > &v, const std::array< double, 4 > &t)
 Interpolate the potential in a curved quadratic tetrahedron.
 
static void Field13 (const std::array< double, 10 > &v, const std::array< double, 4 > &t, double jac[4][4], const double det, double &ex, double &ey, double &ez)
 Interpolate the field in a curved quadratic tetrahedron.
 
static int ReadInteger (char *token, int def, bool &error)
 
static double ReadDouble (char *token, double def, bool &error)
 

Protected Attributes

bool m_is3d = true
 
ElementType m_elementType = ElementType::CurvedTetrahedron
 
std::vector< Element > m_elements
 
std::vector< int > m_elementIndices
 
std::vector< bool > m_degenerate
 
std::vector< std::array< double, 3 > > m_bbMin
 
std::vector< std::array< double, 3 > > m_bbMax
 
std::vector< std::array< std::array< double, 3 >, 4 > > m_w12
 
std::vector< Node > m_nodes
 
std::vector< double > m_pot
 
std::map< std::string, std::vector< double > > m_wpot
 
std::map< std::string, std::vector< std::vector< double > > > m_dwpot
 
std::vector< Material > m_materials
 
std::map< std::string, WeightingFieldCopy > m_wfieldCopies
 
std::vector< double > m_wdtimes
 
bool m_hasBoundingBox = false
 
std::array< double, 3 > m_minBoundingBox = {{0., 0., 0.}}
 
std::array< double, 3 > m_maxBoundingBox = {{0., 0., 0.}}
 
std::array< double, 3 > m_mapmin = {{0., 0., 0.}}
 
std::array< double, 3 > m_mapmax = {{0., 0., 0.}}
 
std::array< double, 3 > m_mapamin = {{0., 0., 0.}}
 
std::array< double, 3 > m_mapamax = {{0., 0., 0.}}
 
std::array< double, 3 > m_mapna = {{0., 0., 0.}}
 
std::array< double, 3 > m_cells = {{0., 0., 0.}}
 
double m_mapvmin = 0.
 
double m_mapvmax = 0.
 
std::array< bool, 3 > m_setang
 
bool m_deleteBackground = true
 
bool m_warning = false
 
unsigned int m_nWarnings = 0
 
bool m_printConvergenceWarnings = true
 
- Protected Attributes inherited from Garfield::Component
std::string m_className = "Component"
 Class name.
 
Geometry * m_geometry = nullptr
 Pointer to the geometry.
 
std::array< double, 3 > m_b0 = {{0., 0., 0.}}
 Constant magnetic field.
 
bool m_ready = false
 Ready for use?
 
bool m_debug = false
 Switch on/off debugging messages.
 
std::array< bool, 3 > m_periodic = {{false, false, false}}
 Simple periodicity in x, y, z.
 
std::array< bool, 3 > m_mirrorPeriodic = {{false, false, false}}
 Mirror periodicity in x, y, z.
 
std::array< bool, 3 > m_axiallyPeriodic = {{false, false, false}}
 Axial periodicity in x, y, z.
 
std::array< bool, 3 > m_rotationSymmetric = {{false, false, false}}
 Rotation symmetry around x-axis, y-axis, z-axis.
 

Friends

class ViewFEMesh
 

Detailed Description

Base class for components based on finite-element field maps.

Definition at line 19 of file ComponentFieldMap.hh.

Member Enumeration Documentation

◆ ElementType

enum class Garfield::ComponentFieldMap::ElementType
strongprotected
Enumerator
Unknown 
Serendipity 
CurvedTetrahedron 

Definition at line 140 of file ComponentFieldMap.hh.

140 {
141 Unknown = 0,
142 Serendipity = 5,
143 CurvedTetrahedron = 13
144 };

Constructor & Destructor Documentation

◆ ComponentFieldMap() [1/2]

Garfield::ComponentFieldMap::ComponentFieldMap ( )
delete

Default constructor.

Referenced by Garfield::ComponentAnsys121::ComponentAnsys121(), Garfield::ComponentAnsys123::ComponentAnsys123(), Garfield::ComponentComsol::ComponentComsol(), Garfield::ComponentCST::ComponentCST(), Garfield::ComponentElmer::ComponentElmer(), Garfield::ComponentElmer::ComponentElmer(), and Garfield::ComponentElmer2d::ComponentElmer2d().

◆ ComponentFieldMap() [2/2]

Garfield::ComponentFieldMap::ComponentFieldMap ( const std::string & name)

Constructor.

Definition at line 19 of file ComponentFieldMap.cc.

20 : Component(name) {}
Garfield::Component::Component
Component()=delete
Default constructor.

◆ ~ComponentFieldMap()

Garfield::ComponentFieldMap::~ComponentFieldMap ( )
virtual

Destructor.

Definition at line 22 of file ComponentFieldMap.cc.

22{}

Member Function Documentation

◆ Check()

bool Garfield::ComponentFieldMap::Check ( )

Check element aspect ratio.

Definition at line 386 of file ComponentFieldMap.cc.

386 {
387 // MAPCHK
388 // Ensure there are some mesh elements.
389 if (!m_ready) {
390 PrintNotReady("Check");
391 return false;
392 }
393 // Compute the range of volumes.
394 const size_t nElements = m_elements.size();
395 double vmin = 0., vmax = 0.;
396 for (size_t i = 0; i < nElements; ++i) {
397 const double v = GetElementVolume(i);
398 if (i == 0) {
399 vmin = vmax = v;
400 } else {
401 vmin = std::min(vmin, v);
402 vmax = std::max(vmax, v);
403 }
404 }
405 // Number of bins.
406 constexpr int nBins = 100;
407 double scale = 1.;
408 std::string unit = "cm";
409 if (m_is3d) {
410 if (vmax < 1.e-9) {
411 unit = "um";
412 scale = 1.e12;
413 } else if (vmax < 1.e-3) {
414 unit = "mm";
415 scale = 1.e3;
416 }
417 } else {
418 if (vmax < 1.e-6) {
419 unit = "um";
420 scale = 1.e8;
421 } else if (vmax < 1.e-2) {
422 unit = "mm";
423 scale = 1.e2;
424 }
425 }
426 vmin *= scale;
427 vmax *= scale;
428 // Check we do have a range and round it.
429 vmin = std::max(0., vmin - 0.1 * (vmax - vmin));
430 vmax = vmax + 0.1 * (vmax - vmin);
431 if (vmin == vmax) {
432 vmin -= 1. + std::abs(vmin);
433 vmax += 1. + std::abs(vmax);
434 }
435 // CALL ROUND(SMIN,SMAX,NCHA,'LARGER,COARSER',STEP)
436 std::string title = m_is3d ? ";volume [" : ";surface [";
437 if (unit == "um") {
438 title += "#mum";
439 } else {
440 title += unit;
441 }
442 if (m_is3d) {
443 title += "^{3}];";
444 } else {
445 title += "^{2}];";
446 }
447 TH1F hElementVolume("hElementVolume", title.c_str(), nBins, vmin, vmax);
448
449 TH1F hAspectRatio("hAspectRatio", ";largest / smallest vertex distance;",
450 nBins, 0., 100.);
451
452 // Loop over all mesh elements.
453 size_t nZero = 0;
454 double rmin = 0., rmax = 0.;
455 for (size_t i = 0; i < nElements; ++i) {
456 double v = 0., dmin = 0., dmax = 0.;
457 if (!GetElement(i, v, dmin, dmax)) return false;
458 // Check for null-sizes.
459 if (dmin <= 0. && !m_degenerate[i]) {
460 std::cerr << m_className << "::Check:\n"
461 << " Found element with zero-length vertex separation.\n";
462 return false;
463 }
464 const double r = dmax / dmin;
465 hAspectRatio.Fill(r);
466 if (v <= 0.) ++nZero;
467 v *= scale;
468 hElementVolume.Fill(v);
469 // Update maxima and minima.
470 if (i == 0) {
471 vmin = vmax = v;
472 rmin = rmax = r;
473 } else {
474 vmin = std::min(vmin, v);
475 vmax = std::max(vmax, v);
476 rmin = std::min(rmin, r);
477 rmax = std::max(rmax, r);
478 }
479 }
480 if (nZero > 0) {
481 std::cerr << m_className << "::Check:\n";
482 if (m_is3d) {
483 std::cerr << " Found " << nZero << " element(s) with zero volume.\n";
484 } else {
485 std::cerr << " Found " << nZero << " element(s) with zero surface.\n";
486 }
487 }
488 TCanvas* c1 = new TCanvas("cAspectRatio", "Aspect ratio", 600, 600);
489 c1->cd();
490 hAspectRatio.DrawCopy();
491 c1->Update();
492 TCanvas* c2 = new TCanvas("cElementVolume", "Element measure", 600, 600);
493 c2->cd();
494 hElementVolume.DrawCopy();
495 c2->Update();
496
497 // Printout.
498 std::cout << m_className << "::Check:\n"
499 << " Smallest Largest\n";
500 std::printf(" Aspect ratios: %15.8f %15.8f\n", rmin, rmax);
501 if (m_is3d) {
502 std::printf(" Volumes [%s3]: %15.8f %15.8f\n", unit.c_str(), vmin,
503 vmax);
504 } else {
505 std::printf(" Surfaces [%s2]: %15.8f %15.8f\n", unit.c_str(), vmin,
506 vmax);
507 }
508 return true;
509}
Garfield::ComponentFieldMap::GetElementVolume
virtual double GetElementVolume(const size_t i) const
Definition ComponentFieldMap.cc:640
Garfield::ComponentFieldMap::GetElement
bool GetElement(const size_t i, double &vol, double &dmin, double &dmax) const
Return the volume and aspect ratio of a mesh element.
Definition ComponentFieldMap.cc:628
Garfield::ComponentFieldMap::m_elements
std::vector< Element > m_elements
Definition ComponentFieldMap.hh:154
Garfield::ComponentFieldMap::PrintNotReady
void PrintNotReady(const std::string &header) const
Definition ComponentFieldMap.cc:2780
Garfield::Component::m_className
std::string m_className
Class name.
Definition Component.hh:359
Garfield::Component::m_ready
bool m_ready
Ready for use?
Definition Component.hh:368

◆ CopyWeightingPotential()

void Garfield::ComponentFieldMap::CopyWeightingPotential ( const std::string & label,
const std::string & labelSource,
const double x,
const double y,
const double z,
const double alpha,
const double beta,
const double gamma )

Makes a weighting potential copy of a imported map which can be translated and rotated.

Parameters
labelname of new electrode
labelSourcename of the source electrode that will be copied
xtranslation in the x-direction.
ytranslation in the y-direction.
ztranslation in the z-direction.
alpharotation around the x-axis.
betarotation around the y-axis.
gammarotation around the z-axis.

Definition at line 2823 of file ComponentFieldMap.cc.

2826 {
2827 // Check if a weighting field with the same label already exists.
2828 if (m_wpot.count(label) > 0) {
2829 std::cout << m_className << "::CopyWeightingPotential:\n"
2830 << " Electrode " << label << " exists already.\n";
2831 return;
2832 }
2833 if (m_wfieldCopies.count(label) > 0) {
2834 std::cout << m_className << "::CopyWeightingPotential:\n"
2835 << " A copy named " << label << " exists already.\n";
2836 return;
2837 }
2838
2839 if (m_wpot.count(labelSource) == 0) {
2840 std::cout << m_className << "::CopyWeightingPotential:\n"
2841 << " Source electrode " << labelSource
2842 << " does not exist.\n";
2843 return;
2844 }
2845
2846 WeightingFieldCopy wfieldCopy;
2847 wfieldCopy.source = labelSource;
2848
2849 TMatrixD Rx(3, 3); // Rotation around the y-axis.
2850 Rx(0, 0) = 1;
2851 Rx(1, 1) = TMath::Cos(-alpha);
2852 Rx(1, 2) = -TMath::Sin(-alpha);
2853 Rx(2, 1) = TMath::Sin(-alpha);
2854 Rx(2, 2) = TMath::Cos(-alpha);
2855
2856 TMatrixD Ry(3, 3); // Rotation around the y-axis.
2857 Ry(1, 1) = 1;
2858 Ry(0, 0) = TMath::Cos(-beta);
2859 Ry(2, 0) = -TMath::Sin(-beta);
2860 Ry(0, 2) = TMath::Sin(-beta);
2861 Ry(2, 2) = TMath::Cos(-beta);
2862
2863 TMatrixD Rz(3, 3); // Rotation around the z-axis.
2864 Rz(2, 2) = 1;
2865 Rz(0, 0) = TMath::Cos(-gamma);
2866 Rz(0, 1) = -TMath::Sin(-gamma);
2867 Rz(1, 0) = TMath::Sin(-gamma);
2868 Rz(1, 1) = TMath::Cos(-gamma);
2869
2870 TVectorD trans(3);
2871 trans(0) = -x;
2872 trans(1) = -y;
2873 trans(2) = -z;
2874 wfieldCopy.rot = Rx * Ry * Rz;
2875 wfieldCopy.trans = trans;
2876 m_wfieldCopies[label] = wfieldCopy;
2877
2878 std::cout << m_className << "::CopyWeightingPotential:\n"
2879 << " Copy named " << label << " of weighting potential "
2880 << labelSource << " made.\n";
2881}
Garfield::ComponentFieldMap::m_wfieldCopies
std::map< std::string, WeightingFieldCopy > m_wfieldCopies
Definition ComponentFieldMap.hh:199
Garfield::ComponentFieldMap::m_wpot
std::map< std::string, std::vector< double > > m_wpot
Definition ComponentFieldMap.hh:174

◆ DelayedWeightingPotential()

double Garfield::ComponentFieldMap::DelayedWeightingPotential ( double x,
double y,
double z,
const double t,
const std::string & label )
overridevirtual

Calculate the delayed weighting potential at a given point and time and for a given electrode.

Parameters
x,y,zcoordinates [cm].
ttime [ns].
labelname of the electrode

Reimplemented from Garfield::Component.

Definition at line 234 of file ComponentFieldMap.cc.

237 {
238 if (m_wdtimes.empty()) return 0.;
239 // Assume no weighting field for times outside the range of available maps.
240 if (tin < m_wdtimes.front()) return 0.;
241 double t = tin;
242 if (tin > m_wdtimes.back()) t = m_wdtimes.back();
243
244 // Do not proceed if not properly initialised.
245 if (!m_ready) return 0.;
246
247 std::string label = label0;
248 if (m_wfieldCopies.count(label0) > 0) {
249 label = m_wfieldCopies[label0].source;
250 TVectorD pos(3);
251 pos(0) = xin;
252 pos(1) = yin;
253 pos(2) = zin;
254 pos = m_wfieldCopies[label0].rot * pos + m_wfieldCopies[label0].trans;
255 xin = pos(0);
256 yin = pos(1);
257 zin = pos(2);
258 }
259
260 // Do not proceed if the requested weighting field does not exist.
261 if (m_dwpot.count(label) == 0) return 0.;
262 if (m_dwpot[label].empty()) return 0.;
263
264 // Copy the coordinates.
265 double x = xin, y = yin, z = zin;
266
267 // Map the coordinates onto field map coordinates.
268 bool xmirr, ymirr, zmirr;
269 double rcoordinate, rotation;
270 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
271
272 if (m_warning) PrintWarning("DelayedWeightingPotential");
273
274 // Find the element that contains this point.
275 double t1, t2, t3, t4, jac[4][4], det;
276
277 int imap = -1;
278 if (m_elementType == ElementType::Serendipity) {
279 imap = FindElement5(x, y, t1, t2, t3, t4, jac, det);
280 } else if (m_elementType == ElementType::CurvedTetrahedron) {
281 imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
282 }
283 if (imap < 0) return 0.;
284
285 // Linear interpolation between time slices
286 int i0;
287 int i1;
288 double f0;
289 double f1;
290
291 TimeInterpolation(t, f0, f1, i0, i1);
292
293 // Get potential value.
294 double dp0 = 0;
295 double dp1 = 0;
296 const Element& element = m_elements[imap];
297 if (m_elementType == ElementType::Serendipity) {
298 if (m_degenerate[imap]) {
299 std::array<double, 6> v0, v1;
300 for (size_t i = 0; i < 6; ++i) {
301 v0[i] = m_dwpot[label][element.emap[i]][i0];
302 v1[i] = m_dwpot[label][element.emap[i]][i1];
303 }
304 dp0 = Potential3(v0, {t1, t2, t3});
305 dp1 = Potential3(v1, {t1, t2, t3});
306 } else {
307 std::array<double, 8> v0, v1;
308 for (size_t i = 0; i < 8; ++i) {
309 v0[i] = m_dwpot[label][element.emap[i]][i0];
310 v1[i] = m_dwpot[label][element.emap[i]][i1];
311 }
312 dp0 = Potential5(v0, {t1, t2});
313 dp1 = Potential5(v1, {t1, t2});
314 }
315 } else if (m_elementType == ElementType::CurvedTetrahedron) {
316 std::array<double, 10> v0, v1;
317 for (size_t i = 0; i < 10; ++i) {
318 v0[i] = m_dwpot[label][element.emap[i]][i0];
319 v1[i] = m_dwpot[label][element.emap[i]][i1];
320 }
321 dp0 = Potential13(v0, {t1, t2, t3, t4});
322 dp1 = Potential13(v1, {t1, t2, t3, t4});
323 }
324
325 return f0 * dp0 + f1 * dp1;
326}
Garfield::ComponentFieldMap::Potential3
static double Potential3(const std::array< double, 6 > &v, const std::array< double, 3 > &t)
Interpolate the potential in a triangle.
Definition ComponentFieldMap.cc:782
Garfield::ComponentFieldMap::FindElement5
int FindElement5(const double x, const double y, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det) const
Find the element for a point in curved quadratic quadrilaterals.
Definition ComponentFieldMap.cc:871
Garfield::ComponentFieldMap::TimeInterpolation
void TimeInterpolation(const double t, double &f0, double &f1, int &i0, int &i1)
Interpolation of potential between two time slices.
Definition ComponentFieldMap.cc:2883
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) const
Find the element for a point in curved quadratic tetrahedra.
Definition ComponentFieldMap.cc:976
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:2446
Garfield::ComponentFieldMap::PrintWarning
void PrintWarning(const std::string &header)
Definition ComponentFieldMap.cc:2773
Garfield::ComponentFieldMap::Potential5
static double Potential5(const std::array< double, 8 > &v, const std::array< double, 2 > &t)
Interpolate the potential in a curved quadrilateral.
Definition ComponentFieldMap.cc:804
Garfield::ComponentFieldMap::Potential13
static double Potential13(const std::array< double, 10 > &v, const std::array< double, 4 > &t)
Interpolate the potential in a curved quadratic tetrahedron.
Definition ComponentFieldMap.cc:839
Garfield::ComponentFieldMap::m_wdtimes
std::vector< double > m_wdtimes
Definition ComponentFieldMap.hh:201
Garfield::ComponentFieldMap::m_dwpot
std::map< std::string, std::vector< std::vector< double > > > m_dwpot
Definition ComponentFieldMap.hh:176

◆ DriftMedium()

void Garfield::ComponentFieldMap::DriftMedium ( const size_t imat)

Flag a field map material as a drift medium.

Definition at line 541 of file ComponentFieldMap.cc.

541 {
542 // Do not proceed if not properly initialised.
543 if (!m_ready) PrintNotReady("DriftMedium");
544
545 // Check value
546 if (imat >= m_materials.size()) {
547 std::cerr << m_className << "::DriftMedium: Index out of range.\n";
548 return;
549 }
550
551 // Make drift medium
552 m_materials[imat].driftmedium = true;
553}
Garfield::ComponentFieldMap::m_materials
std::vector< Material > m_materials
Definition ComponentFieldMap.hh:188

◆ ElectricField() [1/3]

std::array< double, 3 > Garfield::Component::ElectricField ( const double x,
const double y,
const double z )

Calculate the drift field [V/cm] at (x, y, z).

Definition at line 55 of file Component.cc.

43 {
44 double ex = 0., ey = 0., ez = 0.;
45 Medium* medium = nullptr;
46 int status = 0;
47 ElectricField(x, y, z, ex, ey, ez, medium, status);
48 std::array<double, 3> efield = {ex, ey, ez};
49 return efield;
50}
Garfield::ComponentFieldMap::ElectricField
void ElectricField(const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status) override
Definition ComponentFieldMap.cc:31

◆ ElectricField() [2/3]

void Garfield::ComponentFieldMap::ElectricField ( const double x,
const double y,
const double z,
double & ex,
double & ey,
double & ez,
double & v,
Medium *& m,
int & status )
overridevirtual

Calculate the drift field [V/cm] and potential [V] at (x, y, z).

Implements Garfield::Component.

Definition at line 24 of file ComponentFieldMap.cc.

26 {
27 ElectricField(x, y, z, ex, ey, ez, m, status);
28 volt = Potential(x, y, z, m_pot);
29}
Garfield::ComponentFieldMap::Potential
double Potential(const double x, const double y, const double z, const std::vector< double > &potentials) const
Compute the electrostatic/weighting potential.
Definition ComponentFieldMap.cc:133

◆ ElectricField() [3/3]

void Garfield::ComponentFieldMap::ElectricField ( const double x,
const double y,
const double z,
double & ex,
double & ey,
double & ez,
Medium *& m,
int & status )
overridevirtual

Calculate the drift field at given point.

Parameters
x,y,zcoordinates [cm].
ex,ey,ezcomponents of the electric field [V/cm].
mpointer to the medium at this location.
statusstatus flag

Status flags: 

        0: Inside an active medium
      > 0: Inside a wire of type X
-4 ... -1: On the side of a plane where no wires are
       -5: Inside the mesh but not in an active medium
       -6: Outside the mesh
      -10: Unknown potential type (should not occur)
    other: Other cases (should not occur)

Implements Garfield::Component.

Definition at line 31 of file ComponentFieldMap.cc.

33 {
34 // Initial values
35 ex = ey = ez = 0.;
36 m = nullptr;
37 int iel = -1;
38 status = Field(xin, yin, zin, ex, ey, ez, iel, m_pot);
39 if (status < 0 || iel < 0) {
40 if (status == -10) PrintNotReady("ElectricField");
41 return;
42 }
43
44 const auto& element = m_elements[iel];
45 // Drift medium?
46 if (element.matmap >= m_materials.size()) {
47 if (m_debug) {
48 std::cout << m_className << "::ElectricField: "
49 << "Out-of-range material number.\n";
50 }
51 status = -5;
52 return;
53 }
54
55 const auto& mat = m_materials[element.matmap];
56 if (m_debug) {
57 std::cout << " Material " << element.matmap << ", drift flag "
58 << mat.driftmedium << ".\n";
59 }
60 m = mat.medium;
61 status = -5;
62 if (mat.driftmedium) {
63 if (m && m->IsDriftable()) status = 0;
64 }
65}
Garfield::ComponentFieldMap::Field
int Field(const double x, const double y, const double z, double &fx, double &fy, double &fz, int &iel, const std::vector< double > &potentials) const
Compute the electric/weighting field.
Definition ComponentFieldMap.cc:67
Garfield::Component::m_debug
bool m_debug
Switch on/off debugging messages.
Definition Component.hh:371

Referenced by ElectricField().

◆ EnableCheckMapIndices()

void Garfield::ComponentFieldMap::EnableCheckMapIndices ( const bool on = true)
inline

Definition at line 68 of file ComponentFieldMap.hh.

68 {
69 m_checkMultipleElement = on;
70 }

◆ EnableConvergenceWarnings()

void Garfield::ComponentFieldMap::EnableConvergenceWarnings ( const bool on = true)
inline

Enable or disable warnings that the calculation of the local coordinates did not achieve the requested precision.

Definition at line 84 of file ComponentFieldMap.hh.

84 {
85 m_printConvergenceWarnings = on;
86 }

◆ EnableDeleteBackgroundElements()

void Garfield::ComponentFieldMap::EnableDeleteBackgroundElements ( const bool on = true)
inline

Option to eliminate mesh elements in conductors (default: on).

Definition at line 72 of file ComponentFieldMap.hh.

72 {
73 m_deleteBackground = on;
74 }

◆ EnableTetrahedralTreeForElementSearch()

void Garfield::ComponentFieldMap::EnableTetrahedralTreeForElementSearch ( const bool on = true)
inline

Enable or disable the usage of the tetrahedral tree for searching the element in the mesh.

Definition at line 78 of file ComponentFieldMap.hh.

78 {
79 m_useTetrahedralTree = on;
80 }

◆ Field()

int Garfield::ComponentFieldMap::Field ( const double x,
const double y,
const double z,
double & fx,
double & fy,
double & fz,
int & iel,
const std::vector< double > & potentials ) const
protected

Compute the electric/weighting field.

Definition at line 67 of file ComponentFieldMap.cc.

70 {
71
72 // Do not proceed if not properly initialised.
73 if (!m_ready) return -10;
74
75 // Copy the coordinates.
76 double x = xin, y = yin;
77 double z = m_is3d ? zin : 0.;
78
79 // Map the coordinates onto field map coordinates
80 bool xmirr, ymirr, zmirr;
81 double rcoordinate, rotation;
82 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
83
84 if (!m_is3d) {
85 if (zin < m_minBoundingBox[2] || zin > m_maxBoundingBox[2]) {
86 return -5;
87 }
88 }
89
90 // Find the element that contains this point.
91 double t1 = 0., t2 = 0., t3 = 0., t4 = 0.;
92 double jac[4][4];
93 double det = 0.;
94 imap = -1;
95 if (m_elementType == ElementType::Serendipity) {
96 imap = FindElement5(x, y, t1, t2, t3, t4, jac, det);
97 } else if (m_elementType == ElementType::CurvedTetrahedron) {
98 imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
99 }
100 // Stop if the point is not in the mesh.
101 if (imap < 0) {
102 if (m_debug) {
103 std::cerr << m_className << "::Field: ("
104 << x << ", " << y << ", " << z << ") is not in the mesh.\n";
105 }
106 return -6;
107 }
108
109 const Element& element = m_elements[imap];
110 if (m_elementType == ElementType::Serendipity) {
111 if (m_degenerate[imap]) {
112 std::array<double, 6> v;
113 for (size_t i = 0; i < 6; ++i) v[i] = pot[element.emap[i]];
114 Field3(v, {t1, t2, t3}, jac, det, fx, fy);
115 } else {
116 std::array<double, 8> v;
117 for (size_t i = 0; i < 8; ++i) v[i] = pot[element.emap[i]];
118 Field5(v, {t1, t2}, jac, det, fx, fy);
119 }
120 } else if (m_elementType == ElementType::CurvedTetrahedron) {
121 std::array<double, 10> v;
122 for (size_t i = 0; i < 10; ++i) v[i] = pot[element.emap[i]];
123 Field13(v, {t1, t2, t3, t4}, jac, 4 * det, fx, fy, fz);
124 }
125 if (m_debug) {
126 PrintElement("Field", x, y, z, t1, t2, t3, t4, imap, pot);
127 }
128 // Transform field to global coordinates.
129 UnmapFields(fx, fy, fz, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
130 return 0;
131}
Garfield::ComponentFieldMap::Field5
static void Field5(const std::array< double, 8 > &v, const std::array< double, 2 > &t, double jac[4][4], const double det, double &ex, double &ey)
Interpolate the field in a curved quadrilateral.
Definition ComponentFieldMap.cc:816
Garfield::ComponentFieldMap::m_maxBoundingBox
std::array< double, 3 > m_maxBoundingBox
Definition ComponentFieldMap.hh:206
Garfield::ComponentFieldMap::m_minBoundingBox
std::array< double, 3 > m_minBoundingBox
Definition ComponentFieldMap.hh:205
Garfield::ComponentFieldMap::Field3
static void Field3(const std::array< double, 6 > &v, const std::array< double, 3 > &t, double jac[4][4], const double det, double &ex, double &ey)
Interpolate the field in a triangle.
Definition ComponentFieldMap.cc:792
Garfield::ComponentFieldMap::UnmapFields
void UnmapFields(double &ex, double &ey, double &ez, const double xpos, const double ypos, const double zpos, const bool xmirrored, const bool ymirrored, const bool zmirrored, const double rcoordinate, const double rotation) const
Move (ex, ey, ez) to global coordinates.
Definition ComponentFieldMap.cc:2563
Garfield::ComponentFieldMap::Field13
static void Field13(const std::array< double, 10 > &v, const std::array< double, 4 > &t, double jac[4][4], const double det, double &ex, double &ey, double &ez)
Interpolate the field in a curved quadratic tetrahedron.
Definition ComponentFieldMap.cc:851
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 size_t i, const std::vector< double > &potential) const
Definition ComponentFieldMap.cc:2792

Referenced by ElectricField(), and WeightingField().

◆ Field13()

void Garfield::ComponentFieldMap::Field13 ( const std::array< double, 10 > & v,
const std::array< double, 4 > & t,
double jac[4][4],
const double det,
double & ex,
double & ey,
double & ez )
staticprotected

Interpolate the field in a curved quadratic tetrahedron.

Definition at line 851 of file ComponentFieldMap.cc.

854 {
855 std::array<double, 4> g;
856 g[0] = v[0] * (t[0] - 0.25) + v[4] * t[1] + v[5] * t[2] + v[6] * t[3];
857 g[1] = v[1] * (t[1] - 0.25) + v[4] * t[0] + v[7] * t[2] + v[8] * t[3];
858 g[2] = v[2] * (t[2] - 0.25) + v[5] * t[0] + v[7] * t[1] + v[9] * t[3];
859 g[3] = v[3] * (t[3] - 0.25) + v[6] * t[0] + v[8] * t[1] + v[9] * t[2];
860 std::array<double, 3> f = {0., 0., 0.};
861 for (size_t j = 0; j < 4; ++j) {
862 for (size_t i = 0; i < 3; ++i) {
863 f[i] += g[j] * jac[j][i + 1];
864 }
865 }
866 ex = -f[0] * det;
867 ey = -f[1] * det;
868 ez = -f[2] * det;
869}

Referenced by Field().

◆ Field3()

void Garfield::ComponentFieldMap::Field3 ( const std::array< double, 6 > & v,
const std::array< double, 3 > & t,
double jac[4][4],
const double det,
double & ex,
double & ey )
staticprotected

Interpolate the field in a triangle.

Definition at line 792 of file ComponentFieldMap.cc.

794 {
795 std::array<double, 3> g;
796 g[0] = v[0] * (4 * t[0] - 1) + v[3] * 4 * t[1] + v[4] * 4 * t[2];
797 g[1] = v[1] * (4 * t[1] - 1) + v[3] * 4 * t[0] + v[5] * 4 * t[2];
798 g[2] = v[2] * (4 * t[2] - 1) + v[4] * 4 * t[0] + v[5] * 4 * t[1];
799 const double invdet = 1. / det;
800 ex = -(jac[0][1] * g[0] + jac[1][1] * g[1] + jac[2][1] * g[2]) * invdet;
801 ey = -(jac[0][2] * g[0] + jac[1][2] * g[1] + jac[2][2] * g[2]) * invdet;
802}

Referenced by Field().

◆ Field5()

void Garfield::ComponentFieldMap::Field5 ( const std::array< double, 8 > & v,
const std::array< double, 2 > & t,
double jac[4][4],
const double det,
double & ex,
double & ey )
staticprotected

Interpolate the field in a curved quadrilateral.

Definition at line 816 of file ComponentFieldMap.cc.

818 {
819 std::array<double, 2> g;
820 g[0] = (v[0] * (1 - t[1]) * (2 * t[0] + t[1]) +
821 v[1] * (1 - t[1]) * (2 * t[0] - t[1]) +
822 v[2] * (1 + t[1]) * (2 * t[0] + t[1]) +
823 v[3] * (1 + t[1]) * (2 * t[0] - t[1])) *
824 0.25 +
825 v[4] * t[0] * (t[1] - 1) + v[5] * (1 - t[1]) * (1 + t[1]) * 0.5 -
826 v[6] * t[0] * (1 + t[1]) + v[7] * (t[1] - 1) * (t[1] + 1) * 0.5;
827 g[1] = (v[0] * (1 - t[0]) * (t[0] + 2 * t[1]) -
828 v[1] * (1 + t[0]) * (t[0] - 2 * t[1]) +
829 v[2] * (1 + t[0]) * (t[0] + 2 * t[1]) -
830 v[3] * (1 - t[0]) * (t[0] - 2 * t[1])) *
831 0.25 +
832 v[4] * (t[0] - 1) * (t[0] + 1) * 0.5 - v[5] * (1 + t[0]) * t[1] +
833 v[6] * (1 - t[0]) * (1 + t[0]) * 0.5 + v[7] * (t[0] - 1) * t[1];
834 const double invdet = 1. / det;
835 ex = -(g[0] * jac[0][0] + g[1] * jac[1][0]) * invdet;
836 ey = -(g[0] * jac[0][1] + g[1] * jac[1][1]) * invdet;
837}

Referenced by Field().

◆ FindElement13()

int Garfield::ComponentFieldMap::FindElement13 ( const double x,
const double y,
const double z,
double & t1,
double & t2,
double & t3,
double & t4,
double jac[4][4],
double & det ) const
protected

Find the element for a point in curved quadratic tetrahedra.

Definition at line 976 of file ComponentFieldMap.cc.

979 {
980
981 // Backup
982 double jacbak[4][4];
983 double detbak = 1.;
984 double t1bak = 0., t2bak = 0., t3bak = 0., t4bak = 0.;
985 int imapbak = -1;
986
987 // Initial values.
988 t1 = t2 = t3 = t4 = 0.;
989
990 // Verify the count of volumes that contain the point.
991 int nfound = 0;
992 int imap = -1;
993 std::array<double, 10> xn;
994 std::array<double, 10> yn;
995 std::array<double, 10> zn;
996 const auto& elements = (m_useTetrahedralTree && m_octree) ?
997 m_octree->GetElementsInBlock(Vec3(x, y, z)) : m_elementIndices;
998 for (const auto i : elements) {
999 if (x < m_bbMin[i][0] || y < m_bbMin[i][1] || z < m_bbMin[i][2] ||
1000 x > m_bbMax[i][0] || y > m_bbMax[i][1] || z > m_bbMax[i][2]) {
1001 continue;
1002 }
1003 for (size_t j = 0; j < 10; ++j) {
1004 const auto& node = m_nodes[m_elements[i].emap[j]];
1005 xn[j] = node.x;
1006 yn[j] = node.y;
1007 zn[j] = node.z;
1008 }
1009 if (Coordinates13(x, y, z, t1, t2, t3, t4, jac, det, xn, yn, zn, m_w12[i]) != 0) {
1010 continue;
1011 }
1012 if (t1 < 0 || t1 > 1 || t2 < 0 || t2 > 1 || t3 < 0 || t3 > 1 || t4 < 0 ||
1013 t4 > 1) {
1014 continue;
1015 }
1016 if (m_debug) {
1017 std::cout << m_className << "::FindElement13:\n"
1018 << " Found matching element " << i << ".\n";
1019 }
1020 if (!m_checkMultipleElement) return i;
1021 ++nfound;
1022 imap = i;
1023 for (int j = 0; j < 4; ++j) {
1024 for (int k = 0; k < 4; ++k) jacbak[j][k] = jac[j][k];
1025 }
1026 detbak = det;
1027 t1bak = t1;
1028 t2bak = t2;
1029 t3bak = t3;
1030 t4bak = t4;
1031 imapbak = imap;
1032 }
1033
1034 // In checking mode, verify the tetrahedron/triangle count.
1035 if (m_checkMultipleElement) {
1036 if (nfound < 1) {
1037 if (m_debug) {
1038 std::cout << m_className << "::FindElement13:\n"
1039 << " No element matching point ("
1040 << x << ", " << y << ", " << z << ") found.\n";
1041 }
1042 return -1;
1043 }
1044 if (nfound > 1) {
1045 std::cerr << m_className << "::FindElement13:\n"
1046 << " Found << " << nfound << " elements matching point ("
1047 << x << ", " << y << ", " << z << ").\n";
1048 }
1049 for (int j = 0; j < 4; ++j) {
1050 for (int k = 0; k < 4; ++k) jac[j][k] = jacbak[j][k];
1051 }
1052 det = detbak;
1053 t1 = t1bak;
1054 t2 = t2bak;
1055 t3 = t3bak;
1056 t4 = t4bak;
1057 imap = imapbak;
1058 return imap;
1059 }
1060 if (m_debug) {
1061 std::cout << m_className << "::FindElement13:\n"
1062 << " No element matching point (" << x << ", " << y << ", "
1063 << z << ") found.\n";
1064 }
1065 return -1;
1066}
Garfield::ComponentFieldMap::m_bbMax
std::vector< std::array< double, 3 > > m_bbMax
Definition ComponentFieldMap.hh:160
Garfield::ComponentFieldMap::m_w12
std::vector< std::array< std::array< double, 3 >, 4 > > m_w12
Definition ComponentFieldMap.hh:162
Garfield::ComponentFieldMap::m_elementIndices
std::vector< int > m_elementIndices
Definition ComponentFieldMap.hh:155
Garfield::ComponentFieldMap::m_bbMin
std::vector< std::array< double, 3 > > m_bbMin
Definition ComponentFieldMap.hh:159

Referenced by DelayedWeightingPotential(), Field(), GetMedium(), and Potential().

◆ FindElement5()

int Garfield::ComponentFieldMap::FindElement5 ( const double x,
const double y,
double & t1,
double & t2,
double & t3,
double & t4,
double jac[4][4],
double & det ) const
protected

Find the element for a point in curved quadratic quadrilaterals.

Definition at line 871 of file ComponentFieldMap.cc.

874 {
875 // Backup
876 double jacbak[4][4], detbak = 1.;
877 double t1bak = 0., t2bak = 0., t3bak = 0., t4bak = 0.;
878 int imapbak = -1;
879
880 // Initial values.
881 t1 = t2 = t3 = t4 = 0;
882
883 // Verify the count of volumes that contain the point.
884 int nfound = 0;
885 int imap = -1;
886 std::array<double, 8> xn;
887 std::array<double, 8> yn;
888 const auto& elements = (m_useTetrahedralTree && m_octree) ?
889 m_octree->GetElementsInBlock(Vec3(x, y, 0.)) :
890 m_elementIndices;
891 for (const auto i : elements) {
892 if (x < m_bbMin[i][0] || y < m_bbMin[i][1] ||
893 x > m_bbMax[i][0] || y > m_bbMax[i][1]) continue;
894 const Element& element = m_elements[i];
895 if (m_degenerate[i]) {
896 // Degenerate element
897 for (size_t j = 0; j < 6; ++j) {
898 const auto& node = m_nodes[element.emap[j]];
899 xn[j] = node.x;
900 yn[j] = node.y;
901 }
902 if (Coordinates3(x, y, t1, t2, t3, t4, jac, det, xn, yn) != 0) {
903 continue;
904 }
905 if (t1 < 0 || t1 > 1 || t2 < 0 || t2 > 1 || t3 < 0 || t3 > 1) continue;
906 } else {
907 // Non-degenerate element
908 for (size_t j = 0; j < 8; ++j) {
909 const auto& node = m_nodes[element.emap[j]];
910 xn[j] = node.x;
911 yn[j] = node.y;
912 }
913 if (Coordinates5(x, y, t1, t2, t3, t4, jac, det, xn, yn) != 0) {
914 continue;
915 }
916 if (t1 < -1 || t1 > 1 || t2 < -1 || t2 > 1) continue;
917 }
918 ++nfound;
919 imap = i;
920 if (m_debug) {
921 std::cout << m_className << "::FindElement5:\n";
922 if (m_degenerate[i]) {
923 std::cout << " Found matching degenerate element ";
924 } else {
925 std::cout << " Found matching non-degenerate element ";
926 }
927 std::cout << i << ".\n";
928 }
929 if (!m_checkMultipleElement) return i;
930 for (int j = 0; j < 4; ++j) {
931 for (int k = 0; k < 4; ++k) jacbak[j][k] = jac[j][k];
932 }
933 detbak = det;
934 t1bak = t1;
935 t2bak = t2;
936 t3bak = t3;
937 t4bak = t4;
938 imapbak = imap;
939 }
940
941 // In checking mode, verify the element count.
942 if (m_checkMultipleElement) {
943 if (nfound < 1) {
944 if (m_debug) {
945 std::cout << m_className << "::FindElement5:\n"
946 << " No element matching point (" << x << ", " << y
947 << ") found.\n";
948 }
949 return -1;
950 }
951 if (nfound > 1) {
952 std::cout << m_className << "::FindElement5:\n"
953 << " Found " << nfound << " elements matching point ("
954 << x << ", " << y << ").\n";
955 }
956 for (int j = 0; j < 4; ++j) {
957 for (int k = 0; k < 4; ++k) jac[j][k] = jacbak[j][k];
958 }
959 det = detbak;
960 t1 = t1bak;
961 t2 = t2bak;
962 t3 = t3bak;
963 t4 = t4bak;
964 imap = imapbak;
965 return imap;
966 }
967
968 if (m_debug) {
969 std::cout << m_className << "::FindElement5:\n"
970 << " No element matching point (" << x << ", " << y
971 << ") found.\n";
972 }
973 return -1;
974}

Referenced by DelayedWeightingPotential(), Field(), GetMedium(), and Potential().

◆ FindElementCube()

int Garfield::ComponentFieldMap::FindElementCube ( const double x,
const double y,
const double z,
double & t1,
double & t2,
double & t3,
TMatrixD *& jac,
std::vector< TMatrixD * > & dN ) const
protected

Find the element for a point in a cube.

Definition at line 1068 of file ComponentFieldMap.cc.

1071 {
1072 int imap = -1;
1073 const size_t nElements = m_elements.size();
1074 for (size_t i = 0; i < nElements; ++i) {
1075 const Element& element = m_elements[i];
1076 const Node& n3 = m_nodes[element.emap[3]];
1077 if (x < n3.x || y < n3.y || z < n3.z) continue;
1078 const Node& n0 = m_nodes[element.emap[0]];
1079 const Node& n2 = m_nodes[element.emap[2]];
1080 const Node& n7 = m_nodes[element.emap[7]];
1081 if (x < n0.x && y < n2.y && z < n7.z) {
1082 imap = i;
1083 break;
1084 }
1085 }
1086
1087 if (imap < 0) {
1088 if (m_debug) {
1089 std::cout << m_className << "::FindElementCube:\n"
1090 << " Point (" << x << "," << y << "," << z
1091 << ") not in the mesh, it is background or PEC.\n";
1092 const Node& first0 = m_nodes[m_elements.front().emap[0]];
1093 const Node& first2 = m_nodes[m_elements.front().emap[2]];
1094 const Node& first3 = m_nodes[m_elements.front().emap[3]];
1095 const Node& first7 = m_nodes[m_elements.front().emap[7]];
1096 std::cout << " First node (" << first3.x << "," << first3.y << ","
1097 << first3.z << ") in the mesh.\n";
1098 std::cout << " dx= " << (first0.x - first3.x)
1099 << ", dy= " << (first2.y - first3.y)
1100 << ", dz= " << (first7.z - first3.z) << "\n";
1101 const Node& last0 = m_nodes[m_elements.back().emap[0]];
1102 const Node& last2 = m_nodes[m_elements.back().emap[2]];
1103 const Node& last3 = m_nodes[m_elements.back().emap[3]];
1104 const Node& last5 = m_nodes[m_elements.back().emap[5]];
1105 const Node& last7 = m_nodes[m_elements.back().emap[7]];
1106 std::cout << " Last node (" << last5.x << "," << last5.y << ","
1107 << last5.z << ") in the mesh.\n";
1108 std::cout << " dx= " << (last0.x - last3.x)
1109 << ", dy= " << (last2.y - last3.y)
1110 << ", dz= " << (last7.z - last3.z) << "\n";
1111 }
1112 return -1;
1113 }
1114 CoordinatesCube(x, y, z, t1, t2, t3, jac, dN, m_elements[imap]);
1115 return imap;
1116}

◆ GetAspectRatio()

void Garfield::ComponentFieldMap::GetAspectRatio ( const size_t i,
double & dmin,
double & dmax ) const
protectedvirtual

Reimplemented in Garfield::ComponentCST.

Definition at line 671 of file ComponentFieldMap.cc.

672 {
673 if (i >= m_elements.size()) {
674 dmin = dmax = 0.;
675 return;
676 }
677
678 const Element& element = m_elements[i];
679 if (m_elementType == ElementType::CurvedTetrahedron) {
680 const int np = 4;
681 // Loop over all pairs of vertices.
682 for (int j = 0; j < np - 1; ++j) {
683 const Node& nj = m_nodes[element.emap[j]];
684 for (int k = j + 1; k < np; ++k) {
685 const Node& nk = m_nodes[element.emap[k]];
686 // Compute distance.
687 const double dx = nj.x - nk.x;
688 const double dy = nj.y - nk.y;
689 const double dz = nj.z - nk.z;
690 const double dist = sqrt(dx * dx + dy * dy + dz * dz);
691 if (k == 1) {
692 dmin = dmax = dist;
693 } else {
694 if (dist < dmin) dmin = dist;
695 if (dist > dmax) dmax = dist;
696 }
697 }
698 }
699 } else if (m_elementType == ElementType::Serendipity) {
700 const int np = 8;
701 // Loop over all pairs of vertices.
702 for (int j = 0; j < np - 1; ++j) {
703 const Node& nj = m_nodes[element.emap[j]];
704 for (int k = j + 1; k < np; ++k) {
705 const Node& nk = m_nodes[element.emap[k]];
706 // Compute distance.
707 const double dx = nj.x - nk.x;
708 const double dy = nj.y - nk.y;
709 const double dist = sqrt(dx * dx + dy * dy);
710 if (k == 1) {
711 dmin = dmax = dist;
712 } else {
713 if (dist < dmin) dmin = dist;
714 if (dist > dmax) dmax = dist;
715 }
716 }
717 }
718 }
719}
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition DoubleAc.cpp:314

Referenced by GetElement().

◆ GetBoundingBox()

bool Garfield::ComponentFieldMap::GetBoundingBox ( double & xmin,
double & ymin,
double & zmin,
double & xmax,
double & ymax,
double & zmax )
overridevirtual

Get the bounding box coordinates.

Reimplemented from Garfield::Component.

Definition at line 2419 of file ComponentFieldMap.cc.

2421 {
2422 if (!m_ready) return false;
2423
2424 xmin = m_minBoundingBox[0];
2425 xmax = m_maxBoundingBox[0];
2426 ymin = m_minBoundingBox[1];
2427 ymax = m_maxBoundingBox[1];
2428 zmin = m_minBoundingBox[2];
2429 zmax = m_maxBoundingBox[2];
2430 return true;
2431}

◆ GetConductivity()

double Garfield::ComponentFieldMap::GetConductivity ( const size_t imat) const

Return the conductivity of a field map material.

Definition at line 577 of file ComponentFieldMap.cc.

577 {
578 if (imat >= m_materials.size()) {
579 std::cerr << m_className << "::GetConductivity: Index out of range.\n";
580 return -1.;
581 }
582 return m_materials[imat].ohm;
583}

◆ GetElement() [1/2]

bool Garfield::ComponentFieldMap::GetElement ( const size_t i,
double & vol,
double & dmin,
double & dmax ) const

Return the volume and aspect ratio of a mesh element.

Definition at line 628 of file ComponentFieldMap.cc.

629 {
630 if (i >= m_elements.size()) {
631 std::cerr << m_className << "::GetElement: Index out of range.\n";
632 return false;
633 }
634
635 vol = GetElementVolume(i);
636 GetAspectRatio(i, dmin, dmax);
637 return true;
638}
Garfield::ComponentFieldMap::GetAspectRatio
virtual void GetAspectRatio(const size_t i, double &dmin, double &dmax) const
Definition ComponentFieldMap.cc:671

Referenced by Check().

◆ GetElement() [2/2]

bool Garfield::ComponentFieldMap::GetElement ( const size_t i,
size_t & mat,
bool & drift,
std::vector< size_t > & nodes ) const
virtual

Return the material and node indices of a mesh element.

Reimplemented in Garfield::ComponentCST.

Definition at line 721 of file ComponentFieldMap.cc.

722 {
723 if (i >= m_elements.size()) {
724 std::cerr << m_className << "::GetElement: Index out of range.\n";
725 return false;
726 }
727 const auto& element = m_elements[i];
728 mat = element.matmap;
729 drift = m_materials[mat].driftmedium;
730 size_t nNodes = 4;
731 if (m_elementType == ElementType::Serendipity && m_degenerate[i]) {
732 nNodes = 3;
733 }
734 nodes.resize(nNodes);
735 for (size_t j = 0; j < nNodes; ++j) nodes[j] = element.emap[j];
736 return true;
737}

◆ GetElementaryCell()

bool Garfield::ComponentFieldMap::GetElementaryCell ( double & xmin,
double & ymin,
double & zmin,
double & xmax,
double & ymax,
double & zmax )
overridevirtual

Get the coordinates of the elementary cell.

Reimplemented from Garfield::Component.

Definition at line 2433 of file ComponentFieldMap.cc.

2435 {
2436 if (!m_ready) return false;
2437 xmin = m_mapmin[0];
2438 xmax = m_mapmax[0];
2439 ymin = m_mapmin[1];
2440 ymax = m_mapmax[1];
2441 zmin = m_mapmin[2];
2442 zmax = m_mapmax[2];
2443 return true;
2444}
Garfield::ComponentFieldMap::m_mapmin
std::array< double, 3 > m_mapmin
Definition ComponentFieldMap.hh:209
Garfield::ComponentFieldMap::m_mapmax
std::array< double, 3 > m_mapmax
Definition ComponentFieldMap.hh:210

◆ GetElementVolume()

double Garfield::ComponentFieldMap::GetElementVolume ( const size_t i) const
protectedvirtual

Reimplemented in Garfield::ComponentCST.

Definition at line 640 of file ComponentFieldMap.cc.

640 {
641 if (i >= m_elements.size()) return 0.;
642
643 const Element& element = m_elements[i];
644 if (m_elementType == ElementType::CurvedTetrahedron) {
645 const Node& n0 = m_nodes[element.emap[0]];
646 const Node& n1 = m_nodes[element.emap[1]];
647 const Node& n2 = m_nodes[element.emap[2]];
648 const Node& n3 = m_nodes[element.emap[3]];
649 // Uses formula V = |a (dot) b x c|/6
650 // with a => "3", b => "1", c => "2" and origin "0"
651 const double vol = (n3.x - n0.x) * ((n1.y - n0.y) * (n2.z - n0.z) -
652 (n2.y - n0.y) * (n1.z - n0.z)) +
653 (n3.y - n0.y) * ((n1.z - n0.z) * (n2.x - n0.x) -
654 (n2.z - n0.z) * (n1.x - n0.x)) +
655 (n3.z - n0.z) * ((n1.x - n0.x) * (n2.y - n0.y) -
656 (n3.x - n0.x) * (n1.y - n0.y));
657 return fabs(vol) / 6.;
658 } else if (m_elementType == ElementType::Serendipity) {
659 const Node& n0 = m_nodes[element.emap[0]];
660 const Node& n1 = m_nodes[element.emap[1]];
661 const Node& n2 = m_nodes[element.emap[2]];
662 const Node& n3 = m_nodes[element.emap[3]];
663 const double surf = 0.5 *
664 (fabs((n1.x - n0.x) * (n2.y - n0.y) - (n2.x - n0.x) * (n1.y - n0.y)) +
665 fabs((n3.x - n0.x) * (n2.y - n0.y) - (n2.x - n0.x) * (n3.y - n0.y)));
666 return surf;
667 }
668 return 0.;
669}
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition DoubleAc.h:615

Referenced by Check(), and GetElement().

◆ GetMedium() [1/2]

Medium * Garfield::ComponentFieldMap::GetMedium ( const double x,
const double y,
const double z )
overridevirtual

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

Reimplemented from Garfield::Component.

Definition at line 328 of file ComponentFieldMap.cc.

329 {
330 // Copy the coordinates.
331 double x = xin, y = yin;
332 double z = m_is3d ? zin : 0.;
333
334 // Map the coordinates onto field map coordinates.
335 bool xmirr, ymirr, zmirr;
336 double rcoordinate, rotation;
337 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
338
339 if (!m_is3d) {
340 if (zin < m_minBoundingBox[2] || z > m_maxBoundingBox[2]) {
341 return nullptr;
342 }
343 }
344
345 // Do not proceed if not properly initialised.
346 if (!m_ready) {
347 PrintNotReady("GetMedium");
348 return nullptr;
349 }
350 if (m_warning) PrintWarning("GetMedium");
351
352 // Find the element that contains this point.
353 double t1 = 0., t2 = 0., t3 = 0., t4 = 0.;
354 double jac[4][4];
355 double det = 0.;
356 int imap = -1;
357 if (m_elementType == ElementType::Serendipity) {
358 imap = FindElement5(x, y, t1, t2, t3, t4, jac, det);
359 } else if (m_elementType == ElementType::CurvedTetrahedron) {
360 imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
361 }
362 if (imap < 0) {
363 if (m_debug) {
364 std::cerr << m_className << "::GetMedium: (" << x << ", " << y << ", "
365 << z << ") is not in the mesh.\n";
366 }
367 return nullptr;
368 }
369 const Element& element = m_elements[imap];
370 if (element.matmap >= m_materials.size()) {
371 if (m_debug) {
372 std::cerr << m_className << "::GetMedium: Element " << imap
373 << " has out-of-range material number " << element.matmap
374 << ".\n";
375 }
376 return nullptr;
377 }
378 if (m_debug) {
379 PrintElement("GetMedium", x, y, z, t1, t2, t3, t4, imap, m_pot);
380 }
381
382 // Assign a medium.
383 return m_materials[element.matmap].medium;
384}

◆ GetMedium() [2/2]

Medium * Garfield::ComponentFieldMap::GetMedium ( const size_t imat) const

Return the Medium associated to a field map material.

Definition at line 601 of file ComponentFieldMap.cc.

601 {
602 if (imat >= m_materials.size()) {
603 std::cerr << m_className << "::GetMedium: Index out of range.\n";
604 return nullptr;
605 }
606 return m_materials[imat].medium;
607}

◆ GetNode()

bool Garfield::ComponentFieldMap::GetNode ( const size_t i,
double & x,
double & y,
double & z ) const
virtual

Reimplemented in Garfield::ComponentCST.

Definition at line 739 of file ComponentFieldMap.cc.

740 {
741 if (i >= m_nodes.size()) {
742 std::cerr << m_className << "::GetNode: Index out of range.\n";
743 return false;
744 }
745 x = m_nodes[i].x;
746 y = m_nodes[i].y;
747 z = m_nodes[i].z;
748 return true;
749}

◆ GetNumberOfElements()

virtual size_t Garfield::ComponentFieldMap::GetNumberOfElements ( ) const
inlinevirtual

Return the number of mesh elements.

Reimplemented in Garfield::ComponentCST.

Definition at line 56 of file ComponentFieldMap.hh.

56{ return m_elements.size(); }

◆ GetNumberOfMaterials()

size_t Garfield::ComponentFieldMap::GetNumberOfMaterials ( ) const
inline

Return the number of materials in the field map.

Definition at line 41 of file ComponentFieldMap.hh.

41{ return m_materials.size(); }

◆ GetNumberOfNodes()

virtual size_t Garfield::ComponentFieldMap::GetNumberOfNodes ( ) const
inlinevirtual

Reimplemented in Garfield::ComponentCST.

Definition at line 63 of file ComponentFieldMap.hh.

63{ return m_nodes.size(); }

◆ GetPermittivity()

double Garfield::ComponentFieldMap::GetPermittivity ( const size_t imat) const

Return the relative permittivity of a field map material.

Definition at line 569 of file ComponentFieldMap.cc.

569 {
570 if (imat >= m_materials.size()) {
571 std::cerr << m_className << "::GetPermittivity: Index out of range.\n";
572 return -1.;
573 }
574 return m_materials[imat].eps;
575}

◆ GetPotential()

double Garfield::ComponentFieldMap::GetPotential ( const size_t i) const

Definition at line 751 of file ComponentFieldMap.cc.

751 {
752 return i >= m_pot.size() ? 0. : m_pot[i];
753}

◆ GetVoltageRange()

bool Garfield::ComponentFieldMap::GetVoltageRange ( double & vmin,
double & vmax )
inlineoverridevirtual

Calculate the voltage range [V].

Implements Garfield::Component.

Definition at line 116 of file ComponentFieldMap.hh.

116 {
117 vmin = m_mapvmin;
118 vmax = m_mapvmax;
119 return true;
120 }

◆ IsInBoundingBox()

bool Garfield::ComponentFieldMap::IsInBoundingBox ( const double x,
const double y,
const double z ) const
inline

Definition at line 105 of file ComponentFieldMap.hh.

105 {
106 return x >= m_minBoundingBox[0] && x <= m_maxBoundingBox[0] &&
107 y >= m_minBoundingBox[1] && y <= m_maxBoundingBox[1] &&
108 z >= m_minBoundingBox[2] && y <= m_maxBoundingBox[2];
109 }

◆ MapCoordinates()

void Garfield::ComponentFieldMap::MapCoordinates ( double & xpos,
double & ypos,
double & zpos,
bool & xmirrored,
bool & ymirrored,
bool & zmirrored,
double & rcoordinate,
double & rotation ) const
protected

Move (xpos, ypos, zpos) to field map coordinates.

Definition at line 2446 of file ComponentFieldMap.cc.

2449 {
2450 // Initial values
2451 rotation = 0;
2452
2453 // If chamber is periodic, reduce to the cell volume.
2454 xmirrored = false;
2455 if (m_periodic[0]) {
2456 const double xrange = m_mapmax[0] - m_mapmin[0];
2457 xpos = m_mapmin[0] + fmod(xpos - m_mapmin[0], xrange);
2458 if (xpos < m_mapmin[0]) xpos += xrange;
2459 } else if (m_mirrorPeriodic[0]) {
2460 const double xrange = m_mapmax[0] - m_mapmin[0];
2461 double xnew = m_mapmin[0] + fmod(xpos - m_mapmin[0], xrange);
2462 if (xnew < m_mapmin[0]) xnew += xrange;
2463 int nx = int(floor(0.5 + (xnew - xpos) / xrange));
2464 if (nx != 2 * (nx / 2)) {
2465 xnew = m_mapmin[0] + m_mapmax[0] - xnew;
2466 xmirrored = true;
2467 }
2468 xpos = xnew;
2469 }
2470 if (m_axiallyPeriodic[0] && (zpos != 0 || ypos != 0)) {
2471 const double auxr = sqrt(zpos * zpos + ypos * ypos);
2472 double auxphi = atan2(zpos, ypos);
2473 const double phirange = m_mapamax[0] - m_mapamin[0];
2474 const double phim = 0.5 * (m_mapamin[0] + m_mapamax[0]);
2475 rotation = phirange * floor(0.5 + (auxphi - phim) / phirange);
2476 if (auxphi - rotation < m_mapamin[0]) rotation -= phirange;
2477 if (auxphi - rotation > m_mapamax[0]) rotation += phirange;
2478 auxphi = auxphi - rotation;
2479 ypos = auxr * cos(auxphi);
2480 zpos = auxr * sin(auxphi);
2481 }
2482
2483 ymirrored = false;
2484 if (m_periodic[1]) {
2485 const double yrange = m_mapmax[1] - m_mapmin[1];
2486 ypos = m_mapmin[1] + fmod(ypos - m_mapmin[1], yrange);
2487 if (ypos < m_mapmin[1]) ypos += yrange;
2488 } else if (m_mirrorPeriodic[1]) {
2489 const double yrange = m_mapmax[1] - m_mapmin[1];
2490 double ynew = m_mapmin[1] + fmod(ypos - m_mapmin[1], yrange);
2491 if (ynew < m_mapmin[1]) ynew += yrange;
2492 int ny = int(floor(0.5 + (ynew - ypos) / yrange));
2493 if (ny != 2 * (ny / 2)) {
2494 ynew = m_mapmin[1] + m_mapmax[1] - ynew;
2495 ymirrored = true;
2496 }
2497 ypos = ynew;
2498 }
2499 if (m_axiallyPeriodic[1] && (xpos != 0 || zpos != 0)) {
2500 const double auxr = sqrt(xpos * xpos + zpos * zpos);
2501 double auxphi = atan2(xpos, zpos);
2502 const double phirange = (m_mapamax[1] - m_mapamin[1]);
2503 const double phim = 0.5 * (m_mapamin[1] + m_mapamax[1]);
2504 rotation = phirange * floor(0.5 + (auxphi - phim) / phirange);
2505 if (auxphi - rotation < m_mapamin[1]) rotation -= phirange;
2506 if (auxphi - rotation > m_mapamax[1]) rotation += phirange;
2507 auxphi = auxphi - rotation;
2508 zpos = auxr * cos(auxphi);
2509 xpos = auxr * sin(auxphi);
2510 }
2511
2512 zmirrored = false;
2513 if (m_periodic[2]) {
2514 const double zrange = m_mapmax[2] - m_mapmin[2];
2515 zpos = m_mapmin[2] + fmod(zpos - m_mapmin[2], zrange);
2516 if (zpos < m_mapmin[2]) zpos += zrange;
2517 } else if (m_mirrorPeriodic[2]) {
2518 const double zrange = m_mapmax[2] - m_mapmin[2];
2519 double znew = m_mapmin[2] + fmod(zpos - m_mapmin[2], zrange);
2520 if (znew < m_mapmin[2]) znew += zrange;
2521 int nz = int(floor(0.5 + (znew - zpos) / zrange));
2522 if (nz != 2 * (nz / 2)) {
2523 znew = m_mapmin[2] + m_mapmax[2] - znew;
2524 zmirrored = true;
2525 }
2526 zpos = znew;
2527 }
2528 if (m_axiallyPeriodic[2] && (ypos != 0 || xpos != 0)) {
2529 const double auxr = sqrt(ypos * ypos + xpos * xpos);
2530 double auxphi = atan2(ypos, xpos);
2531 const double phirange = m_mapamax[2] - m_mapamin[2];
2532 const double phim = 0.5 * (m_mapamin[2] + m_mapamax[2]);
2533 rotation = phirange * floor(0.5 + (auxphi - phim) / phirange);
2534 if (auxphi - rotation < m_mapamin[2]) rotation -= phirange;
2535 if (auxphi - rotation > m_mapamax[2]) rotation += phirange;
2536 auxphi = auxphi - rotation;
2537 xpos = auxr * cos(auxphi);
2538 ypos = auxr * sin(auxphi);
2539 }
2540
2541 // If we have a rotationally symmetric field map, store coordinates.
2542 rcoordinate = 0;
2543 double zcoordinate = 0;
2544 if (m_rotationSymmetric[0]) {
2545 rcoordinate = sqrt(ypos * ypos + zpos * zpos);
2546 zcoordinate = xpos;
2547 } else if (m_rotationSymmetric[1]) {
2548 rcoordinate = sqrt(xpos * xpos + zpos * zpos);
2549 zcoordinate = ypos;
2550 } else if (m_rotationSymmetric[2]) {
2551 rcoordinate = sqrt(xpos * xpos + ypos * ypos);
2552 zcoordinate = zpos;
2553 }
2554
2555 if (m_rotationSymmetric[0] || m_rotationSymmetric[1] ||
2556 m_rotationSymmetric[2]) {
2557 xpos = rcoordinate;
2558 ypos = zcoordinate;
2559 zpos = 0;
2560 }
2561}
Garfield::ComponentFieldMap::m_mapamin
std::array< double, 3 > m_mapamin
Definition ComponentFieldMap.hh:211
Garfield::ComponentFieldMap::m_mapamax
std::array< double, 3 > m_mapamax
Definition ComponentFieldMap.hh:212
Garfield::Component::m_rotationSymmetric
std::array< bool, 3 > m_rotationSymmetric
Rotation symmetry around x-axis, y-axis, z-axis.
Definition Component.hh:380
Garfield::Component::m_mirrorPeriodic
std::array< bool, 3 > m_mirrorPeriodic
Mirror periodicity in x, y, z.
Definition Component.hh:376
Garfield::Component::m_periodic
std::array< bool, 3 > m_periodic
Simple periodicity in x, y, z.
Definition Component.hh:374
Garfield::Component::m_axiallyPeriodic
std::array< bool, 3 > m_axiallyPeriodic
Axial periodicity in x, y, z.
Definition Component.hh:378
Heed::cos
DoubleAc cos(const DoubleAc &f)
Definition DoubleAc.cpp:432
Heed::sin
DoubleAc sin(const DoubleAc &f)
Definition DoubleAc.cpp:384

Referenced by Garfield::ComponentCST::Coordinate2Index(), DelayedWeightingPotential(), Field(), GetMedium(), and Potential().

◆ NotDriftMedium()

void Garfield::ComponentFieldMap::NotDriftMedium ( const size_t imat)

Flag a field map materials as a non-drift medium.

Definition at line 555 of file ComponentFieldMap.cc.

555 {
556 // Do not proceed if not properly initialised.
557 if (!m_ready) PrintNotReady("NotDriftMedium");
558
559 // Check value
560 if (imat >= m_materials.size()) {
561 std::cerr << m_className << "::NotDriftMedium: Index out of range.\n";
562 return;
563 }
564
565 // Make drift medium
566 m_materials[imat].driftmedium = false;
567}

◆ Potential()

double Garfield::ComponentFieldMap::Potential ( const double x,
const double y,
const double z,
const std::vector< double > & potentials ) const
protected

Compute the electrostatic/weighting potential.

Definition at line 133 of file ComponentFieldMap.cc.

135 {
136
137 // Do not proceed if not properly initialised.
138 if (!m_ready) return 0.;
139
140 // Copy the coordinates.
141 double x = xin, y = yin;
142 double z = m_is3d ? zin : 0.;
143
144 // Map the coordinates onto field map coordinates.
145 bool xmirr, ymirr, zmirr;
146 double rcoordinate, rotation;
147 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
148
149 if (!m_is3d) {
150 if (zin < m_minBoundingBox[2] || zin > m_maxBoundingBox[2]) {
151 return 0.;
152 }
153 }
154
155 // Find the element that contains this point.
156 double t1 = 0., t2 = 0., t3 = 0., t4 = 0.;
157 double jac[4][4];
158 double det = 0.;
159 int imap = -1;
160 if (m_elementType == ElementType::Serendipity) {
161 imap = FindElement5(x, y, t1, t2, t3, t4, jac, det);
162 } else if (m_elementType == ElementType::CurvedTetrahedron) {
163 imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
164 }
165 if (imap < 0) return 0.;
166
167 double volt = 0.;
168 const Element& element = m_elements[imap];
169 if (m_elementType == ElementType::Serendipity) {
170 if (m_degenerate[imap]) {
171 std::array<double, 6> v;
172 for (size_t i = 0; i < 6; ++i) v[i] = pot[element.emap[i]];
173 volt = Potential3(v, {t1, t2, t3});
174 } else {
175 std::array<double, 8> v;
176 for (size_t i = 0; i < 8; ++i) v[i] = pot[element.emap[i]];
177 volt = Potential5(v, {t1, t2});
178 }
179 } else if (m_elementType == ElementType::CurvedTetrahedron) {
180 std::array<double, 10> v;
181 for (size_t i = 0; i < 10; ++i) v[i] = pot[element.emap[i]];
182 volt = Potential13(v, {t1, t2, t3, t4});
183 }
184 if (m_debug) {
185 PrintElement("Potential", x, y, z, t1, t2, t3, t4, imap, pot);
186 }
187 return volt;
188}

Referenced by ElectricField(), and WeightingPotential().

◆ Potential13()

double Garfield::ComponentFieldMap::Potential13 ( const std::array< double, 10 > & v,
const std::array< double, 4 > & t )
staticprotected

Interpolate the potential in a curved quadratic tetrahedron.

Definition at line 839 of file ComponentFieldMap.cc.

840 {
841 double sum = 0.;
842 for (size_t i = 0; i < 4; ++i) {
843 sum += v[i] * t[i] * (t[i] - 0.5);
844 }
845 sum *= 2;
846 sum += 4 * (v[4] * t[0] * t[1] + v[5] * t[0] * t[2] + v[6] * t[0] * t[3] +
847 v[7] * t[1] * t[2] + v[8] * t[1] * t[3] + v[9] * t[2] * t[3]);
848 return sum;
849}

Referenced by DelayedWeightingPotential(), and Potential().

◆ Potential3()

double Garfield::ComponentFieldMap::Potential3 ( const std::array< double, 6 > & v,
const std::array< double, 3 > & t )
staticprotected

Interpolate the potential in a triangle.

Definition at line 782 of file ComponentFieldMap.cc.

783 {
784 double sum = 0.;
785 for (size_t i = 0; i < 3; ++i) {
786 sum += v[i] * t[i] * (2 * t[i] - 1);
787 }
788 sum += 4 * (v[3] * t[0] * t[1] + v[4] * t[0] * t[2] + v[5] * t[1] * t[2]);
789 return sum;
790}

Referenced by DelayedWeightingPotential(), and Potential().

◆ Potential5()

double Garfield::ComponentFieldMap::Potential5 ( const std::array< double, 8 > & v,
const std::array< double, 2 > & t )
staticprotected

Interpolate the potential in a curved quadrilateral.

Definition at line 804 of file ComponentFieldMap.cc.

805 {
806 return -v[0] * (1 - t[0]) * (1 - t[1]) * (1 + t[0] + t[1]) * 0.25 -
807 v[1] * (1 + t[0]) * (1 - t[1]) * (1 - t[0] + t[1]) * 0.25 -
808 v[2] * (1 + t[0]) * (1 + t[1]) * (1 - t[0] - t[1]) * 0.25 -
809 v[3] * (1 - t[0]) * (1 + t[1]) * (1 + t[0] - t[1]) * 0.25 +
810 v[4] * (1 - t[0]) * (1 + t[0]) * (1 - t[1]) * 0.5 +
811 v[5] * (1 + t[0]) * (1 + t[1]) * (1 - t[1]) * 0.5 +
812 v[6] * (1 - t[0]) * (1 + t[0]) * (1 + t[1]) * 0.5 +
813 v[7] * (1 - t[0]) * (1 + t[1]) * (1 - t[1]) * 0.5;
814}

Referenced by DelayedWeightingPotential(), and Potential().

◆ Prepare()

void Garfield::ComponentFieldMap::Prepare ( )
protected

Definition at line 2121 of file ComponentFieldMap.cc.

2121 {
2122 // Establish the ranges.
2123 SetRange();
2124 UpdatePeriodicity();
2125 std::cout << m_className << "::Prepare:\n"
2126 << " Caching the bounding boxes of all elements...";
2127 CalculateElementBoundingBoxes();
2128 std::cout << " done.\n";
2129 // Initialize the tetrahedral tree.
2130 if (InitializeTetrahedralTree()) {
2131 std::cout << " Initialized tetrahedral tree.\n";
2132 }
2133 m_elementIndices.resize(m_elements.size());
2134 std::iota(m_elementIndices.begin(), m_elementIndices.end(), 0);
2135 // Precompute terms for interpolation in linear tetrahedra.
2136 if (m_elementType == ElementType::CurvedTetrahedron) {
2137 std::array<double, 10> xn;
2138 std::array<double, 10> yn;
2139 std::array<double, 10> zn;
2140 for (const auto& element : m_elements) {
2141 for (size_t j = 0; j < 10; ++j) {
2142 const auto& node = m_nodes[element.emap[j]];
2143 xn[j] = node.x;
2144 yn[j] = node.y;
2145 zn[j] = node.z;
2146 }
2147 m_w12.emplace_back(Weights12(xn, yn, zn));
2148 }
2149 }
2150}
Garfield::ComponentFieldMap::UpdatePeriodicity
void UpdatePeriodicity() override
Verify periodicities.
Definition ComponentFieldMap.hh:244

Referenced by Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentComsol::Initialise(), Garfield::ComponentCST::Initialise(), Garfield::ComponentElmer2d::Initialise(), and Garfield::ComponentElmer::Initialise().

◆ PrintCouldNotOpen()

void Garfield::ComponentFieldMap::PrintCouldNotOpen ( const std::string & header,
const std::string & filename ) const
protected

Definition at line 2785 of file ComponentFieldMap.cc.

2786 {
2787 std::cerr << m_className << "::" << header << ":\n"
2788 << " Could not open file " << filename << " for reading.\n"
2789 << " The file perhaps does not exist.\n";
2790}

Referenced by Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentComsol::Initialise(), Garfield::ComponentCST::Initialise(), Garfield::ComponentCST::Initialise(), Garfield::ComponentElmer2d::Initialise(), Garfield::ComponentElmer::Initialise(), Garfield::ComponentComsol::SetDynamicWeightingPotential(), and Garfield::ComponentCST::SetWeightingField().

◆ PrintElement()

void Garfield::ComponentFieldMap::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 size_t i,
const std::vector< double > & potential ) const
protected

Definition at line 2792 of file ComponentFieldMap.cc.

2797 {
2798 std::cout << m_className << "::" << header << ":\n"
2799 << " Global = (" << x << ", " << y << ", " << z << ")\n"
2800 << " Local = (" << t1 << ", " << t2 << ", " << t3 << ", " << t4
2801 << ")\n";
2802 if (m_degenerate[i]) std::cout << " Element is degenerate.\n";
2803 std::cout << " Node x y z V\n";
2804 unsigned int nN = 0;
2805 if (m_elementType == ElementType::Serendipity) {
2806 if (m_degenerate[i]) {
2807 nN = 6;
2808 } else {
2809 nN = 8;
2810 }
2811 } else if (m_elementType == ElementType::CurvedTetrahedron) {
2812 nN = 10;
2813 }
2814 const auto& element = m_elements[i];
2815 for (unsigned int ii = 0; ii < nN; ++ii) {
2816 const Node& node = m_nodes[element.emap[ii]];
2817 const double v = pot[element.emap[ii]];
2818 printf(" %-5d %12g %12g %12g %12g\n", element.emap[ii], node.x, node.y,
2819 node.z, v);
2820 }
2821}

Referenced by Field(), GetMedium(), and Potential().

◆ PrintMaterials()

void Garfield::ComponentFieldMap::PrintMaterials ( )

List all currently defined materials.

Definition at line 511 of file ComponentFieldMap.cc.

511 {
512 // Do not proceed if not properly initialised.
513 if (!m_ready) PrintNotReady("PrintMaterials");
514
515 if (m_materials.empty()) {
516 std::cerr << m_className << "::PrintMaterials:\n"
517 << " No materials are currently defined.\n";
518 return;
519 }
520
521 const size_t nMaterials = m_materials.size();
522 std::cout << m_className << "::PrintMaterials:\n"
523 << " Currently " << nMaterials << " materials are defined.\n"
524 << " Index Permittivity Resistivity Notes\n";
525 for (size_t i = 0; i < nMaterials; ++i) {
526 printf(" %5zu %12g %12g", i, m_materials[i].eps, m_materials[i].ohm);
527 if (m_materials[i].medium) {
528 std::string name = m_materials[i].medium->GetName();
529 std::cout << " " << name;
530 if (m_materials[i].medium->IsDriftable()) std::cout << ", drift medium";
531 if (m_materials[i].medium->IsIonisable()) std::cout << ", ionisable";
532 }
533 if (m_materials[i].driftmedium) {
534 std::cout << " (drift medium)\n";
535 } else {
536 std::cout << "\n";
537 }
538 }
539}

Referenced by Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), and Garfield::ComponentCST::Initialise().

◆ PrintNotReady()

void Garfield::ComponentFieldMap::PrintNotReady ( const std::string & header) const
protected

Definition at line 2780 of file ComponentFieldMap.cc.

2780 {
2781 std::cerr << m_className << "::" << header << ":\n"
2782 << " Field map not yet initialised.\n";
2783}

Referenced by Check(), DriftMedium(), ElectricField(), GetMedium(), NotDriftMedium(), PrintMaterials(), Garfield::ComponentAnsys121::SetWeightingField(), Garfield::ComponentAnsys123::SetWeightingField(), Garfield::ComponentElmer2d::SetWeightingField(), Garfield::ComponentElmer::SetWeightingField(), UpdatePeriodicity2d(), and UpdatePeriodicityCommon().

◆ PrintRange()

void Garfield::ComponentFieldMap::PrintRange ( )

Show x, y, z, V and angular ranges.

Definition at line 2390 of file ComponentFieldMap.cc.

2390 {
2391 std::cout << m_className << "::PrintRange:\n";
2392 std::cout << " Dimensions of the elementary block\n";
2393 printf(" %15g < x < %-15g cm,\n", m_mapmin[0], m_mapmax[0]);
2394 printf(" %15g < y < %-15g cm,\n", m_mapmin[1], m_mapmax[1]);
2395 printf(" %15g < z < %-15g cm,\n", m_mapmin[2], m_mapmax[2]);
2396 printf(" %15g < V < %-15g V.\n", m_mapvmin, m_mapvmax);
2397
2398 std::cout << " Periodicities\n";
2399 const std::array<std::string, 3> axes = {{"x", "y", "z"}};
2400 for (unsigned int i = 0; i < 3; ++i) {
2401 std::cout << " " << axes[i] << ":";
2402 if (m_periodic[i]) {
2403 std::cout << " simple with length " << m_cells[i] << " cm";
2404 }
2405 if (m_mirrorPeriodic[i]) {
2406 std::cout << " mirror with length " << m_cells[i] << " cm";
2407 }
2408 if (m_axiallyPeriodic[i]) {
2409 std::cout << " axial " << int(0.5 + m_mapna[i]) << "-fold repetition";
2410 }
2411 if (m_rotationSymmetric[i]) std::cout << " rotational symmetry";
2412 if (!(m_periodic[i] || m_mirrorPeriodic[i] || m_axiallyPeriodic[i] ||
2413 m_rotationSymmetric[i]))
2414 std::cout << " none";
2415 std::cout << "\n";
2416 }
2417}
Garfield::ComponentFieldMap::m_mapna
std::array< double, 3 > m_mapna
Definition ComponentFieldMap.hh:213
Garfield::ComponentFieldMap::m_cells
std::array< double, 3 > m_cells
Definition ComponentFieldMap.hh:214

Referenced by SetRange(), and UpdatePeriodicityCommon().

◆ PrintWarning()

void Garfield::ComponentFieldMap::PrintWarning ( const std::string & header)
protected

Definition at line 2773 of file ComponentFieldMap.cc.

2773 {
2774 if (!m_warning || m_nWarnings > 10) return;
2775 std::cerr << m_className << "::" << header << ":\n"
2776 << " Warnings have been issued for this field map.\n";
2777 ++m_nWarnings;
2778}

Referenced by DelayedWeightingPotential(), GetMedium(), and Garfield::ComponentCST::WeightingField().

◆ ReadDouble()

double Garfield::ComponentFieldMap::ReadDouble ( char * token,
double def,
bool & error )
staticprotected

Definition at line 2660 of file ComponentFieldMap.cc.

2660 {
2661 if (!token) {
2662 error = true;
2663 return def;
2664 }
2665 return atof(token);
2666}

Referenced by Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentCST::Initialise(), Garfield::ComponentElmer2d::Initialise(), Garfield::ComponentElmer::Initialise(), and Garfield::ComponentCST::SetWeightingField().

◆ ReadInteger()

int Garfield::ComponentFieldMap::ReadInteger ( char * token,
int def,
bool & error )
staticprotected

Definition at line 2651 of file ComponentFieldMap.cc.

2651 {
2652 if (!token) {
2653 error = true;
2654 return def;
2655 }
2656
2657 return atoi(token);
2658}

Referenced by Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentCST::Initialise(), Garfield::ComponentElmer2d::Initialise(), Garfield::ComponentElmer::Initialise(), and Garfield::ComponentCST::SetWeightingField().

◆ Reset()

void Garfield::ComponentFieldMap::Reset ( )
overrideprotectedvirtual

Reset the component.

Implements Garfield::Component.

Definition at line 2099 of file ComponentFieldMap.cc.

2099 {
2100 m_ready = false;
2101
2102 m_elements.clear();
2103 m_degenerate.clear();
2104 m_bbMin.clear();
2105 m_bbMax.clear();
2106 m_w12.clear();
2107 m_nodes.clear();
2108 m_pot.clear();
2109 m_wpot.clear();
2110 m_dwpot.clear();
2111 m_wfieldCopies.clear();
2112 m_materials.clear();
2113 m_hasBoundingBox = false;
2114 m_warning = false;
2115 m_nWarnings = 0;
2116
2117 m_octree.reset(nullptr);
2118 m_cacheElemBoundingBoxes = false;
2119}

Referenced by Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentComsol::Initialise(), Garfield::ComponentCST::Initialise(), Garfield::ComponentElmer2d::Initialise(), and Garfield::ComponentElmer::Initialise().

◆ ScalingFactor()

double Garfield::ComponentFieldMap::ScalingFactor ( std::string unit)
staticprotected

Definition at line 2637 of file ComponentFieldMap.cc.

2637 {
2638 std::transform(unit.begin(), unit.end(), unit.begin(), toupper);
2639 if (unit == "MUM" || unit == "MICRON" || unit == "MICROMETER") {
2640 return 0.0001;
2641 } else if (unit == "MM" || unit == "MILLIMETER") {
2642 return 0.1;
2643 } else if (unit == "CM" || unit == "CENTIMETER") {
2644 return 1.0;
2645 } else if (unit == "M" || unit == "METER") {
2646 return 100.0;
2647 }
2648 return -1.;
2649}

Referenced by Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentComsol::Initialise(), Garfield::ComponentCST::Initialise(), Garfield::ComponentCST::Initialise(), Garfield::ComponentElmer2d::Initialise(), and Garfield::ComponentElmer::Initialise().

◆ SetDefaultDriftMedium()

bool Garfield::ComponentFieldMap::SetDefaultDriftMedium ( )
protected

Find lowest epsilon, check for eps = 0, set default drift media flags.

Definition at line 755 of file ComponentFieldMap.cc.

755 {
756 // Find lowest epsilon and set drift medium flags.
757 const size_t nMaterials = m_materials.size();
758 double epsmin = -1;
759 size_t iepsmin = 0;
760 for (size_t i = 0; i < nMaterials; ++i) {
761 m_materials[i].driftmedium = false;
762 if (m_materials[i].eps < 0) continue;
763 // Check for eps == 0.
764 if (m_materials[i].eps == 0) {
765 std::cerr << m_className << "::SetDefaultDriftMedium:\n"
766 << " Material " << i << " has zero permittivity.\n";
767 m_materials[i].eps = -1.;
768 } else if (epsmin < 0. || epsmin > m_materials[i].eps) {
769 epsmin = m_materials[i].eps;
770 iepsmin = i;
771 }
772 }
773 if (epsmin < 0.) {
774 std::cerr << m_className << "::SetDefaultDriftMedium:\n"
775 << " Found no material with positive permittivity.\n";
776 return false;
777 }
778 m_materials[iepsmin].driftmedium = true;
779 return true;
780}

Referenced by Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentComsol::Initialise(), Garfield::ComponentCST::Initialise(), Garfield::ComponentElmer2d::Initialise(), and Garfield::ComponentElmer::Initialise().

◆ SetGas()

void Garfield::ComponentFieldMap::SetGas ( Medium * medium)

Associate all field map materials with a relative permittivity of unity to a given Medium class.

Definition at line 609 of file ComponentFieldMap.cc.

609 {
610 if (!medium) {
611 std::cerr << m_className << "::SetGas: Null pointer.\n";
612 return;
613 }
614 size_t nMatch = 0;
615 const size_t nMaterials = m_materials.size();
616 for (size_t i = 0; i < nMaterials; ++i) {
617 if (fabs(m_materials[i].eps - 1.) > 1.e-4) continue;
618 m_materials[i].medium = medium;
619 std::cout << m_className << "::SetGas: Associating material " << i
620 << " with " << medium->GetName() << ".\n";
621 ++nMatch;
622 }
623 if (nMatch == 0) {
624 std::cerr << m_className << "::SetGas: Found no material with eps = 1.\n";
625 }
626}

Referenced by main().

◆ SetMedium()

void Garfield::ComponentFieldMap::SetMedium ( const size_t imat,
Medium * medium )

Associate a field map material with a Medium object.

Definition at line 585 of file ComponentFieldMap.cc.

585 {
586 if (imat >= m_materials.size()) {
587 std::cerr << m_className << "::SetMedium: Index out of range.\n";
588 return;
589 }
590 if (!medium) {
591 std::cerr << m_className << "::SetMedium: Null pointer.\n";
592 return;
593 }
594 if (m_debug) {
595 std::cout << m_className << "::SetMedium: Associated material " << imat
596 << " with medium " << medium->GetName() << ".\n";
597 }
598 m_materials[imat].medium = medium;
599}

◆ SetRange()

void Garfield::ComponentFieldMap::SetRange ( )
protectedvirtual

Reimplemented in Garfield::ComponentCST.

Definition at line 2299 of file ComponentFieldMap.cc.

2299 {
2300 // Initial values
2301 m_mapmin.fill(0.);
2302 m_mapmax.fill(0.);
2303 m_mapamin.fill(0.);
2304 m_mapamax.fill(0.);
2305 m_mapvmin = m_mapvmax = 0.;
2306 m_setang.fill(false);
2307
2308 // Make sure the required data is available.
2309 if (!m_ready || m_nodes.empty()) {
2310 std::cerr << m_className << "::SetRange: Field map not yet set.\n";
2311 return;
2312 }
2313
2314 m_mapvmin = *std::min_element(std::begin(m_pot), std::end(m_pot));
2315 m_mapvmax = *std::max_element(std::begin(m_pot), std::end(m_pot));
2316
2317 // Loop over the nodes.
2318 m_mapmin[0] = m_mapmax[0] = m_nodes[0].x;
2319 m_mapmin[1] = m_mapmax[1] = m_nodes[0].y;
2320 m_mapmin[2] = m_mapmax[2] = m_nodes[0].z;
2321
2322 for (const auto& node : m_nodes) {
2323 const std::array<double, 3> pos = {{node.x, node.y, node.z}};
2324 for (unsigned int i = 0; i < 3; ++i) {
2325 m_mapmin[i] = std::min(m_mapmin[i], pos[i]);
2326 m_mapmax[i] = std::max(m_mapmax[i], pos[i]);
2327 }
2328
2329 if (node.y != 0 || node.z != 0) {
2330 const double ang = atan2(node.z, node.y);
2331 if (m_setang[0]) {
2332 m_mapamin[0] = std::min(m_mapamin[0], ang);
2333 m_mapamax[0] = std::max(m_mapamax[0], ang);
2334 } else {
2335 m_mapamin[0] = m_mapamax[0] = ang;
2336 m_setang[0] = true;
2337 }
2338 }
2339
2340 if (node.z != 0 || node.x != 0) {
2341 const double ang = atan2(node.x, node.z);
2342 if (m_setang[1]) {
2343 m_mapamin[1] = std::min(m_mapamin[1], ang);
2344 m_mapamax[1] = std::max(m_mapamax[1], ang);
2345 } else {
2346 m_mapamin[1] = m_mapamax[1] = ang;
2347 m_setang[1] = true;
2348 }
2349 }
2350
2351 if (node.x != 0 || node.y != 0) {
2352 const double ang = atan2(node.y, node.x);
2353 if (m_setang[2]) {
2354 m_mapamin[2] = std::min(m_mapamin[2], ang);
2355 m_mapamax[2] = std::max(m_mapamax[2], ang);
2356 } else {
2357 m_mapamin[2] = m_mapamax[2] = ang;
2358 m_setang[2] = true;
2359 }
2360 }
2361 }
2362
2363 // Fix the angular ranges.
2364 for (unsigned int i = 0; i < 3; ++i) {
2365 if (m_mapamax[i] - m_mapamin[i] > Pi) {
2366 const double aux = m_mapamin[i];
2367 m_mapamin[i] = m_mapamax[i];
2368 m_mapamax[i] = aux + TwoPi;
2369 }
2370 }
2371
2372 // Set provisional cell dimensions.
2373 m_minBoundingBox[0] = m_mapmin[0];
2374 m_maxBoundingBox[0] = m_mapmax[0];
2375 m_minBoundingBox[1] = m_mapmin[1];
2376 m_maxBoundingBox[1] = m_mapmax[1];
2377 if (m_is3d) {
2378 m_minBoundingBox[2] = m_mapmin[2];
2379 m_maxBoundingBox[2] = m_mapmax[2];
2380 } else {
2381 m_mapmin[2] = m_minBoundingBox[2];
2382 m_mapmax[2] = m_maxBoundingBox[2];
2383 }
2384 m_hasBoundingBox = true;
2385
2386 // Display the range if requested.
2387 if (m_debug) PrintRange();
2388}
Garfield::ComponentFieldMap::PrintRange
void PrintRange()
Show x, y, z, V and angular ranges.
Definition ComponentFieldMap.cc:2390
Garfield::ComponentFieldMap::m_setang
std::array< bool, 3 > m_setang
Definition ComponentFieldMap.hh:219

Referenced by Prepare().

◆ TimeInterpolation()

void Garfield::ComponentFieldMap::TimeInterpolation ( const double t,
double & f0,
double & f1,
int & i0,
int & i1 )
protected

Interpolation of potential between two time slices.

Definition at line 2883 of file ComponentFieldMap.cc.

2884 {
2885 const auto it1 = std::upper_bound(m_wdtimes.cbegin(), m_wdtimes.cend(), t);
2886 const auto it0 = std::prev(it1);
2887
2888 const double dt = t - *it0;
2889 i0 = it0 - m_wdtimes.cbegin();
2890 i1 = it1 - m_wdtimes.cbegin();
2891
2892 f1 = dt / (*it1 - *it0);
2893 f0 = 1. - f1;
2894}

Referenced by DelayedWeightingPotential().

◆ UnmapFields()

void Garfield::ComponentFieldMap::UnmapFields ( double & ex,
double & ey,
double & ez,
const double xpos,
const double ypos,
const double zpos,
const bool xmirrored,
const bool ymirrored,
const bool zmirrored,
const double rcoordinate,
const double rotation ) const
protected

Move (ex, ey, ez) to global coordinates.

Definition at line 2563 of file ComponentFieldMap.cc.

2566 {
2567 // Apply mirror imaging.
2568 if (xmirrored) ex = -ex;
2569 if (ymirrored) ey = -ey;
2570 if (zmirrored) ez = -ez;
2571
2572 // Rotate the field.
2573 double er, theta;
2574 if (m_axiallyPeriodic[0]) {
2575 er = sqrt(ey * ey + ez * ez);
2576 theta = atan2(ez, ey);
2577 theta += rotation;
2578 ey = er * cos(theta);
2579 ez = er * sin(theta);
2580 }
2581 if (m_axiallyPeriodic[1]) {
2582 er = sqrt(ez * ez + ex * ex);
2583 theta = atan2(ex, ez);
2584 theta += rotation;
2585 ez = er * cos(theta);
2586 ex = er * sin(theta);
2587 }
2588 if (m_axiallyPeriodic[2]) {
2589 er = sqrt(ex * ex + ey * ey);
2590 theta = atan2(ey, ex);
2591 theta += rotation;
2592 ex = er * cos(theta);
2593 ey = er * sin(theta);
2594 }
2595
2596 // Take care of symmetry.
2597 double eaxis;
2598 er = ex;
2599 eaxis = ey;
2600
2601 // Rotational symmetry
2602 if (m_rotationSymmetric[0]) {
2603 if (rcoordinate <= 0) {
2604 ex = eaxis;
2605 ey = 0;
2606 ez = 0;
2607 } else {
2608 ex = eaxis;
2609 ey = er * ypos / rcoordinate;
2610 ez = er * zpos / rcoordinate;
2611 }
2612 }
2613 if (m_rotationSymmetric[1]) {
2614 if (rcoordinate <= 0) {
2615 ex = 0;
2616 ey = eaxis;
2617 ez = 0;
2618 } else {
2619 ex = er * xpos / rcoordinate;
2620 ey = eaxis;
2621 ez = er * zpos / rcoordinate;
2622 }
2623 }
2624 if (m_rotationSymmetric[2]) {
2625 if (rcoordinate <= 0) {
2626 ex = 0;
2627 ey = 0;
2628 ez = eaxis;
2629 } else {
2630 ex = er * xpos / rcoordinate;
2631 ey = er * ypos / rcoordinate;
2632 ez = eaxis;
2633 }
2634 }
2635}

Referenced by Field().

◆ UpdatePeriodicity()

void Garfield::ComponentFieldMap::UpdatePeriodicity ( )
inlineoverrideprotectedvirtual

Verify periodicities.

Implements Garfield::Component.

Definition at line 244 of file ComponentFieldMap.hh.

244 {
245 if (!m_is3d) UpdatePeriodicity2d();
246 UpdatePeriodicityCommon();
247 }

Referenced by Garfield::ComponentCST::Initialise(), Prepare(), and Garfield::ComponentCST::ShiftComponent().

◆ UpdatePeriodicity2d()

void Garfield::ComponentFieldMap::UpdatePeriodicity2d ( )
protected

Definition at line 2271 of file ComponentFieldMap.cc.

2271 {
2272 // Check the required data is available.
2273 if (!m_ready) {
2274 PrintNotReady("UpdatePeriodicity2d");
2275 return;
2276 }
2277
2278 // No z-periodicity in 2d
2279 if (m_periodic[2] || m_mirrorPeriodic[2]) {
2280 std::cerr << m_className << "::UpdatePeriodicity2d:\n"
2281 << " Simple or mirror periodicity along z\n"
2282 << " requested for a 2d map; reset.\n";
2283 m_periodic[2] = false;
2284 m_mirrorPeriodic[2] = false;
2285 m_warning = true;
2286 }
2287
2288 // Only z-axial periodicity in 2d maps
2289 if (m_axiallyPeriodic[0] || m_axiallyPeriodic[1]) {
2290 std::cerr << m_className << "::UpdatePeriodicity2d:\n"
2291 << " Axial symmetry has been requested \n"
2292 << " around x or y for a 2d map; reset.\n";
2293 m_axiallyPeriodic[0] = false;
2294 m_axiallyPeriodic[1] = false;
2295 m_warning = true;
2296 }
2297}

Referenced by UpdatePeriodicity().

◆ UpdatePeriodicityCommon()

void Garfield::ComponentFieldMap::UpdatePeriodicityCommon ( )
protected

Definition at line 2152 of file ComponentFieldMap.cc.

2152 {
2153 // Check the required data is available.
2154 if (!m_ready) {
2155 PrintNotReady("UpdatePeriodicityCommon");
2156 return;
2157 }
2158
2159 for (size_t i = 0; i < 3; ++i) {
2160 // No regular and mirror periodicity at the same time.
2161 if (m_periodic[i] && m_mirrorPeriodic[i]) {
2162 std::cerr << m_className << "::UpdatePeriodicityCommon:\n"
2163 << " Both simple and mirror periodicity requested. Reset.\n";
2164 m_periodic[i] = false;
2165 m_mirrorPeriodic[i] = false;
2166 m_warning = true;
2167 }
2168 // In case of axial periodicity,
2169 // the range must be an integral part of two pi.
2170 if (m_axiallyPeriodic[i]) {
2171 if (m_mapamin[i] >= m_mapamax[i]) {
2172 m_mapna[i] = 0;
2173 } else {
2174 m_mapna[i] = TwoPi / (m_mapamax[i] - m_mapamin[i]);
2175 }
2176 if (fabs(m_mapna[i] - int(0.5 + m_mapna[i])) > 0.001 ||
2177 m_mapna[i] < 1.5) {
2178 std::cerr << m_className << "::UpdatePeriodicityCommon:\n"
2179 << " Axial symmetry has been requested but map does not\n"
2180 << " cover an integral fraction of 2 pi. Reset.\n";
2181 m_axiallyPeriodic[i] = false;
2182 m_warning = true;
2183 }
2184 }
2185 }
2186
2187 // Not more than 1 rotational symmetry
2188 if ((m_rotationSymmetric[0] && m_rotationSymmetric[1]) ||
2189 (m_rotationSymmetric[0] && m_rotationSymmetric[2]) ||
2190 (m_rotationSymmetric[1] && m_rotationSymmetric[2])) {
2191 std::cerr << m_className << "::UpdatePeriodicityCommon:\n"
2192 << " Only one rotational symmetry allowed; reset.\n";
2193 m_rotationSymmetric.fill(false);
2194 m_warning = true;
2195 }
2196
2197 // No rotational symmetry as well as axial periodicity
2198 if ((m_rotationSymmetric[0] || m_rotationSymmetric[1] ||
2199 m_rotationSymmetric[2]) &&
2200 (m_axiallyPeriodic[0] || m_axiallyPeriodic[1] || m_axiallyPeriodic[2])) {
2201 std::cerr << m_className << "::UpdatePeriodicityCommon:\n"
2202 << " Not allowed to combine rotational symmetry\n"
2203 << " and axial periodicity; reset.\n";
2204 m_axiallyPeriodic.fill(false);
2205 m_rotationSymmetric.fill(false);
2206 m_warning = true;
2207 }
2208
2209 // In case of rotational symmetry, the x-range should not straddle 0.
2210 if (m_rotationSymmetric[0] || m_rotationSymmetric[1] ||
2211 m_rotationSymmetric[2]) {
2212 if (m_mapmin[0] * m_mapmax[0] < 0) {
2213 std::cerr << m_className << "::UpdatePeriodicityCommon:\n"
2214 << " Rotational symmetry requested, \n"
2215 << " but x-range straddles 0; reset.\n";
2216 m_rotationSymmetric.fill(false);
2217 m_warning = true;
2218 }
2219 }
2220
2221 // Recompute the cell ranges.
2222 for (size_t i = 0; i < 3; ++i) {
2223 m_minBoundingBox[i] = m_mapmin[i];
2224 m_maxBoundingBox[i] = m_mapmax[i];
2225 m_cells[i] = fabs(m_mapmax[i] - m_mapmin[i]);
2226 }
2227 for (size_t i = 0; i < 3; ++i) {
2228 if (!m_rotationSymmetric[i]) continue;
2229 const double r = std::max(fabs(m_mapmin[0]), fabs(m_mapmax[0]));
2230 m_minBoundingBox.fill(-r);
2231 m_maxBoundingBox.fill(+r);
2232 m_minBoundingBox[i] = m_mapmin[1];
2233 m_maxBoundingBox[i] = m_mapmax[1];
2234 break;
2235 }
2236
2237 if (m_axiallyPeriodic[0]) {
2238 const double yzmax = std::max({fabs(m_mapmin[1]), fabs(m_mapmax[1]),
2239 fabs(m_mapmin[2]), fabs(m_mapmax[2])});
2240 m_minBoundingBox[1] = -yzmax;
2241 m_maxBoundingBox[1] = +yzmax;
2242 m_minBoundingBox[2] = -yzmax;
2243 m_maxBoundingBox[2] = +yzmax;
2244 } else if (m_axiallyPeriodic[1]) {
2245 const double xzmax = std::max({fabs(m_mapmin[0]), fabs(m_mapmax[0]),
2246 fabs(m_mapmin[2]), fabs(m_mapmax[2])});
2247 m_minBoundingBox[0] = -xzmax;
2248 m_maxBoundingBox[0] = +xzmax;
2249 m_minBoundingBox[2] = -xzmax;
2250 m_maxBoundingBox[2] = +xzmax;
2251 } else if (m_axiallyPeriodic[2]) {
2252 const double xymax = std::max({fabs(m_mapmin[0]), fabs(m_mapmax[0]),
2253 fabs(m_mapmin[1]), fabs(m_mapmax[1])});
2254 m_minBoundingBox[0] = -xymax;
2255 m_maxBoundingBox[0] = +xymax;
2256 m_minBoundingBox[1] = -xymax;
2257 m_maxBoundingBox[1] = +xymax;
2258 }
2259
2260 for (size_t i = 0; i < 3; ++i) {
2261 if (m_periodic[i] || m_mirrorPeriodic[i]) {
2262 m_minBoundingBox[i] = -INFINITY;
2263 m_maxBoundingBox[i] = +INFINITY;
2264 }
2265 }
2266
2267 // Display the range if requested.
2268 if (m_debug) PrintRange();
2269}

Referenced by UpdatePeriodicity().

◆ WeightingField()

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

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

Parameters
x,y,zcoordinates [cm].
wx,wy,wzcomponents of the weighting field [1/cm].
labelname of the electrode

Reimplemented from Garfield::Component.

Definition at line 190 of file ComponentFieldMap.cc.

192 {
193 // Initial values.
194 wx = wy = wz = 0;
195
196 // Do not proceed if not properly initialised.
197 if (!m_ready) return;
198
199 // Do not proceed if the requested weighting field does not exist.
200 if (m_wpot.count(label) == 0) return;
201 if (m_wpot[label].empty()) return;
202
203 int iel = -1;
204 Field(xin, yin, zin, wx, wy, wz, iel, m_wpot[label]);
205}

◆ WeightingPotential()

double Garfield::ComponentFieldMap::WeightingPotential ( const double x,
const double y,
const double z,
const std::string & label )
overridevirtual

Calculate the weighting potential at a given point.

Parameters
x,y,zcoordinates [cm].
labelname of the electrode.
Returns
weighting potential [dimensionless].

Reimplemented from Garfield::Component.

Definition at line 207 of file ComponentFieldMap.cc.

208 {
209 // Do not proceed if not properly initialised.
210 if (!m_ready) return 0.;
211
212 // TODO! From ComponentComsol:
213 // if (!CheckInRange(xin, yin, zin)) return 0.;
214
215 std::string label = label0;
216 if (m_wfieldCopies.count(label0) > 0) {
217 label = m_wfieldCopies[label0].source;
218 TVectorD pos(3);
219 pos(0) = xin;
220 pos(1) = yin;
221 pos(2) = zin;
222 pos = m_wfieldCopies[label0].rot * pos + m_wfieldCopies[label0].trans;
223 xin = pos(0);
224 yin = pos(1);
225 zin = pos(2);
226 }
227 // Do not proceed if the requested weighting field does not exist.
228 if (m_wpot.count(label) == 0) return 0.;
229 if (m_wpot[label].empty()) return 0.;
230
231 return Potential(xin, yin, zin, m_wpot[label]);
232}

Friends And Related Symbol Documentation

◆ ViewFEMesh

friend class ViewFEMesh
friend

Definition at line 135 of file ComponentFieldMap.hh.

Referenced by ViewFEMesh.

Member Data Documentation

◆ m_bbMax

std::vector<std::array<double, 3> > Garfield::ComponentFieldMap::m_bbMax
protected

Definition at line 160 of file ComponentFieldMap.hh.

Referenced by FindElement13(), FindElement5(), and Reset().

◆ m_bbMin

std::vector<std::array<double, 3> > Garfield::ComponentFieldMap::m_bbMin
protected

Definition at line 159 of file ComponentFieldMap.hh.

Referenced by FindElement13(), FindElement5(), and Reset().

◆ m_cells

std::array<double, 3> Garfield::ComponentFieldMap::m_cells = {{0., 0., 0.}}
protected

Definition at line 214 of file ComponentFieldMap.hh.

214{{0., 0., 0.}};

Referenced by PrintRange(), and UpdatePeriodicityCommon().

◆ m_degenerate

std::vector<bool> Garfield::ComponentFieldMap::m_degenerate
protected

Definition at line 157 of file ComponentFieldMap.hh.

Referenced by Check(), DelayedWeightingPotential(), Field(), FindElement5(), GetElement(), Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentComsol::Initialise(), Garfield::ComponentElmer2d::Initialise(), Garfield::ComponentElmer::Initialise(), Potential(), PrintElement(), and Reset().

◆ m_deleteBackground

bool Garfield::ComponentFieldMap::m_deleteBackground = true
protected

Definition at line 222 of file ComponentFieldMap.hh.

Referenced by Garfield::ComponentCST::ComponentCST(), EnableDeleteBackgroundElements(), Garfield::ComponentAnsys121::Initialise(), and Garfield::ComponentAnsys123::Initialise().

◆ m_dwpot

std::map<std::string, std::vector<std::vector<double> > > Garfield::ComponentFieldMap::m_dwpot
protected

Definition at line 176 of file ComponentFieldMap.hh.

Referenced by DelayedWeightingPotential(), Reset(), and Garfield::ComponentComsol::SetDynamicWeightingPotential().

◆ m_elementIndices

std::vector<int> Garfield::ComponentFieldMap::m_elementIndices
protected

Definition at line 155 of file ComponentFieldMap.hh.

Referenced by FindElement13(), FindElement5(), and Prepare().

◆ m_elements

std::vector<Element> Garfield::ComponentFieldMap::m_elements
protected

Definition at line 154 of file ComponentFieldMap.hh.

Referenced by Check(), DelayedWeightingPotential(), ElectricField(), Field(), FindElement13(), FindElement5(), FindElementCube(), GetAspectRatio(), GetElement(), GetElement(), GetElementVolume(), GetMedium(), GetNumberOfElements(), Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentComsol::Initialise(), Garfield::ComponentElmer2d::Initialise(), Garfield::ComponentElmer::Initialise(), Potential(), Prepare(), PrintElement(), and Reset().

◆ m_elementType

ElementType Garfield::ComponentFieldMap::m_elementType = ElementType::CurvedTetrahedron
protected

Definition at line 145 of file ComponentFieldMap.hh.

Referenced by Garfield::ComponentAnsys121::ComponentAnsys121(), Garfield::ComponentElmer2d::ComponentElmer2d(), DelayedWeightingPotential(), Field(), GetAspectRatio(), GetElement(), GetElementVolume(), GetMedium(), Potential(), Prepare(), and PrintElement().

◆ m_hasBoundingBox

bool Garfield::ComponentFieldMap::m_hasBoundingBox = false
protected

Definition at line 204 of file ComponentFieldMap.hh.

Referenced by Reset(), Garfield::ComponentCST::SetRange(), and SetRange().

◆ m_is3d

bool Garfield::ComponentFieldMap::m_is3d = true
protected

Definition at line 138 of file ComponentFieldMap.hh.

Referenced by Check(), Garfield::ComponentAnsys121::ComponentAnsys121(), Garfield::ComponentElmer2d::ComponentElmer2d(), Field(), GetMedium(), Potential(), Garfield::ComponentCST::SetRange(), SetRange(), and UpdatePeriodicity().

◆ m_mapamax

std::array<double, 3> Garfield::ComponentFieldMap::m_mapamax = {{0., 0., 0.}}
protected

Definition at line 212 of file ComponentFieldMap.hh.

212{{0., 0., 0.}};

Referenced by MapCoordinates(), SetRange(), and UpdatePeriodicityCommon().

◆ m_mapamin

std::array<double, 3> Garfield::ComponentFieldMap::m_mapamin = {{0., 0., 0.}}
protected

Definition at line 211 of file ComponentFieldMap.hh.

211{{0., 0., 0.}};

Referenced by MapCoordinates(), SetRange(), and UpdatePeriodicityCommon().

◆ m_mapmax

std::array<double, 3> Garfield::ComponentFieldMap::m_mapmax = {{0., 0., 0.}}
protected

Definition at line 210 of file ComponentFieldMap.hh.

210{{0., 0., 0.}};

Referenced by GetElementaryCell(), MapCoordinates(), PrintRange(), Garfield::ComponentCST::SetRange(), SetRange(), Garfield::ComponentAnsys121::SetRangeZ(), Garfield::ComponentElmer2d::SetRangeZ(), and UpdatePeriodicityCommon().

◆ m_mapmin

std::array<double, 3> Garfield::ComponentFieldMap::m_mapmin = {{0., 0., 0.}}
protected

Definition at line 209 of file ComponentFieldMap.hh.

209{{0., 0., 0.}};

Referenced by GetElementaryCell(), MapCoordinates(), PrintRange(), Garfield::ComponentCST::SetRange(), SetRange(), Garfield::ComponentAnsys121::SetRangeZ(), Garfield::ComponentElmer2d::SetRangeZ(), and UpdatePeriodicityCommon().

◆ m_mapna

std::array<double, 3> Garfield::ComponentFieldMap::m_mapna = {{0., 0., 0.}}
protected

Definition at line 213 of file ComponentFieldMap.hh.

213{{0., 0., 0.}};

Referenced by PrintRange(), and UpdatePeriodicityCommon().

◆ m_mapvmax

double Garfield::ComponentFieldMap::m_mapvmax = 0.
protected

Definition at line 217 of file ComponentFieldMap.hh.

Referenced by GetVoltageRange(), PrintRange(), Garfield::ComponentCST::SetRange(), and SetRange().

◆ m_mapvmin

double Garfield::ComponentFieldMap::m_mapvmin = 0.
protected

Definition at line 216 of file ComponentFieldMap.hh.

Referenced by GetVoltageRange(), PrintRange(), Garfield::ComponentCST::SetRange(), and SetRange().

◆ m_materials

std::vector<Material> Garfield::ComponentFieldMap::m_materials
protected

Definition at line 188 of file ComponentFieldMap.hh.

Referenced by DriftMedium(), ElectricField(), GetConductivity(), Garfield::ComponentCST::GetElement(), GetElement(), Garfield::ComponentCST::GetMedium(), GetMedium(), GetMedium(), GetNumberOfMaterials(), GetPermittivity(), Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentComsol::Initialise(), Garfield::ComponentCST::Initialise(), Garfield::ComponentCST::Initialise(), Garfield::ComponentElmer2d::Initialise(), Garfield::ComponentElmer::Initialise(), NotDriftMedium(), PrintMaterials(), Reset(), SetDefaultDriftMedium(), SetGas(), SetMedium(), and Garfield::ComponentCST::WeightingField().

◆ m_maxBoundingBox

std::array<double, 3> Garfield::ComponentFieldMap::m_maxBoundingBox = {{0., 0., 0.}}
protected

Definition at line 206 of file ComponentFieldMap.hh.

206{{0., 0., 0.}};

Referenced by Garfield::ComponentAnsys121::ComponentAnsys121(), Garfield::ComponentCST::ComponentCST(), Garfield::ComponentElmer2d::ComponentElmer2d(), Field(), GetBoundingBox(), GetMedium(), IsInBoundingBox(), Potential(), Garfield::ComponentCST::SetRange(), SetRange(), Garfield::ComponentAnsys121::SetRangeZ(), Garfield::ComponentCST::SetRangeZ(), Garfield::ComponentElmer2d::SetRangeZ(), and UpdatePeriodicityCommon().

◆ m_minBoundingBox

std::array<double, 3> Garfield::ComponentFieldMap::m_minBoundingBox = {{0., 0., 0.}}
protected

Definition at line 205 of file ComponentFieldMap.hh.

205{{0., 0., 0.}};

Referenced by Garfield::ComponentAnsys121::ComponentAnsys121(), Garfield::ComponentCST::ComponentCST(), Garfield::ComponentElmer2d::ComponentElmer2d(), Field(), GetBoundingBox(), GetMedium(), IsInBoundingBox(), Potential(), Garfield::ComponentCST::SetRange(), SetRange(), Garfield::ComponentAnsys121::SetRangeZ(), Garfield::ComponentCST::SetRangeZ(), Garfield::ComponentElmer2d::SetRangeZ(), and UpdatePeriodicityCommon().

◆ m_nodes

std::vector<Node> Garfield::ComponentFieldMap::m_nodes
protected

Definition at line 169 of file ComponentFieldMap.hh.

Referenced by FindElement13(), FindElement5(), FindElementCube(), GetAspectRatio(), GetElementVolume(), GetNode(), GetNumberOfNodes(), Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentComsol::Initialise(), Garfield::ComponentElmer2d::Initialise(), Garfield::ComponentElmer::Initialise(), Prepare(), PrintElement(), Reset(), Garfield::ComponentComsol::SetDynamicWeightingPotential(), SetRange(), Garfield::ComponentAnsys121::SetWeightingField(), Garfield::ComponentAnsys123::SetWeightingField(), Garfield::ComponentElmer2d::SetWeightingField(), Garfield::ComponentElmer::SetWeightingField(), and Garfield::ComponentComsol::SetWeightingPotential().

◆ m_nWarnings

unsigned int Garfield::ComponentFieldMap::m_nWarnings = 0
protected

Definition at line 226 of file ComponentFieldMap.hh.

Referenced by PrintWarning(), and Reset().

◆ m_pot

std::vector<double> Garfield::ComponentFieldMap::m_pot
protected

Definition at line 172 of file ComponentFieldMap.hh.

Referenced by ElectricField(), ElectricField(), GetMedium(), GetPotential(), Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentComsol::Initialise(), Garfield::ComponentElmer2d::Initialise(), Garfield::ComponentElmer::Initialise(), Reset(), and SetRange().

◆ m_printConvergenceWarnings

bool Garfield::ComponentFieldMap::m_printConvergenceWarnings = true
protected

Definition at line 230 of file ComponentFieldMap.hh.

Referenced by EnableConvergenceWarnings().

◆ m_setang

std::array<bool, 3> Garfield::ComponentFieldMap::m_setang
protected

Definition at line 219 of file ComponentFieldMap.hh.

Referenced by SetRange().

◆ m_w12

std::vector<std::array<std::array<double, 3>, 4> > Garfield::ComponentFieldMap::m_w12
protected

Definition at line 162 of file ComponentFieldMap.hh.

Referenced by FindElement13(), Prepare(), and Reset().

◆ m_warning

bool Garfield::ComponentFieldMap::m_warning = false
protected

Definition at line 225 of file ComponentFieldMap.hh.

Referenced by DelayedWeightingPotential(), GetMedium(), PrintWarning(), Reset(), UpdatePeriodicity2d(), UpdatePeriodicityCommon(), and Garfield::ComponentCST::WeightingField().

◆ m_wdtimes

std::vector<double> Garfield::ComponentFieldMap::m_wdtimes
protected

Definition at line 201 of file ComponentFieldMap.hh.

Referenced by DelayedWeightingPotential(), Garfield::ComponentComsol::SetDynamicWeightingPotential(), Garfield::ComponentComsol::SetTimeInterval(), and TimeInterpolation().

◆ m_wfieldCopies

std::map<std::string, WeightingFieldCopy> Garfield::ComponentFieldMap::m_wfieldCopies
protected

Definition at line 199 of file ComponentFieldMap.hh.

Referenced by CopyWeightingPotential(), DelayedWeightingPotential(), Reset(), and WeightingPotential().

◆ m_wpot

std::map<std::string, std::vector<double> > Garfield::ComponentFieldMap::m_wpot
protected

Definition at line 174 of file ComponentFieldMap.hh.

Referenced by CopyWeightingPotential(), Garfield::ComponentComsol::Initialise(), Reset(), Garfield::ComponentComsol::SetDynamicWeightingPotential(), Garfield::ComponentAnsys121::SetWeightingField(), Garfield::ComponentAnsys123::SetWeightingField(), Garfield::ComponentElmer2d::SetWeightingField(), Garfield::ComponentElmer::SetWeightingField(), Garfield::ComponentComsol::SetWeightingPotential(), WeightingField(), and WeightingPotential().

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