Garfield++ 4.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Public Member Functions | Protected Member Functions | List of all members
Garfield::ComponentComsol Class Reference

Component for importing and interpolating Comsol field maps. More...

#include <ComponentComsol.hh>

+ Inheritance diagram for Garfield::ComponentComsol:

Public Member Functions

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

Protected Member Functions

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

Additional Inherited Members

- Static Protected Member Functions inherited from Garfield::ComponentFieldMap
static double ScalingFactor (std::string unit)
 
- Protected Attributes inherited from Garfield::ComponentFieldMap
bool m_is3d = true
 
std::vector< Element > m_elements
 
std::vector< Node > m_nodes
 
std::vector< Material > m_materials
 
std::vector< std::string > m_wfields
 
std::vector< bool > m_wfieldsOk
 
std::vector< bool > m_dwfieldsOk
 
std::vector< double > m_wdtimes
 
bool m_hasBoundingBox = false
 
std::array< double, 3 > m_minBoundingBox
 
std::array< double, 3 > m_maxBoundingBox
 
std::array< double, 3 > m_mapmin
 
std::array< double, 3 > m_mapmax
 
std::array< double, 3 > m_mapamin
 
std::array< double, 3 > m_mapamax
 
std::array< double, 3 > m_mapna
 
std::array< double, 3 > m_cells
 
double m_mapvmin = 0.
 
double m_mapvmax = 0.
 
std::array< bool, 3 > m_setang
 
bool m_deleteBackground = true
 
bool m_warning = false
 
unsigned int m_nWarnings = 0
 
bool m_printConvergenceWarnings = true
 
- Protected Attributes inherited from Garfield::Component
std::string m_className = "Component"
 Class name.
 
Geometry * m_geometry = nullptr
 Pointer to the geometry.
 
std::array< double, 3 > m_b0 = {{0., 0., 0.}}
 Constant magnetic field.
 
bool m_ready = false
 Ready for use?
 
bool m_debug = false
 Switch on/off debugging messages.
 
std::array< bool, 3 > m_periodic = {{false, false, false}}
 Simple periodicity in x, y, z.
 
std::array< bool, 3 > m_mirrorPeriodic = {{false, false, false}}
 Mirror periodicity in x, y, z.
 
std::array< bool, 3 > m_axiallyPeriodic = {{false, false, false}}
 Axial periodicity in x, y, z.
 
std::array< bool, 3 > m_rotationSymmetric = {{false, false, false}}
 Rotation symmetry around x-axis, y-axis, z-axis.
 

Detailed Description

Component for importing and interpolating Comsol field maps.

Definition at line 9 of file ComponentComsol.hh.

Constructor & Destructor Documentation

◆ ComponentComsol() [1/2]

Garfield::ComponentComsol::ComponentComsol ( )

Default constructor.

Definition at line 52 of file ComponentComsol.cc.

52: ComponentFieldMap("Comsol") {}
Garfield::ComponentFieldMap::ComponentFieldMap
ComponentFieldMap()=delete
Default constructor.

◆ ComponentComsol() [2/2]

Garfield::ComponentComsol::ComponentComsol ( const std::string &  mesh,
const std::string &  mplist,
const std::string &  field,
const std::string &  unit = "m" 
)

Constructor from file names.

Definition at line 54 of file ComponentComsol.cc.

58 : ComponentComsol() {
59 Initialise(mesh, mplist, field, unit);
60}
Garfield::ComponentComsol::Initialise
bool Initialise(const std::string &header="mesh.mphtxt", const std::string &mplist="dielectrics.dat", const std::string &field="field.txt", const std::string &unit="m")
Definition: ComponentComsol.cc:62
Garfield::ComponentComsol::ComponentComsol
ComponentComsol()
Default constructor.
Definition: ComponentComsol.cc:52

◆ ~ComponentComsol()

Garfield::ComponentComsol::~ComponentComsol ( )
inline

Destructor.

Definition at line 17 of file ComponentComsol.hh.

17{}

Member Function Documentation

◆ DelayedWeightingPotential()

double Garfield::ComponentComsol::DelayedWeightingPotential ( const double  x,
const double  y,
const 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 818 of file ComponentComsol.cc.

822 {
823 if (m_wdtimes.empty()) return 0.;
824 // Assume no weighting field for times outside the range of available maps.
825 if (t < m_wdtimes.front() || t > m_wdtimes.back()) return 0.;
826
827 // Do not proceed if not properly initialised.
828 if (!m_ready) return 0.;
829 // Look for the label.
830 const size_t iw = GetWeightingFieldIndex(label);
831 // Do not proceed if the requested weighting field does not exist.
832 if (iw == m_wfields.size()) return 0.;
833 // Check if the weighting field is properly initialised.
834
835 // Copy the coordinates.
836 double x = xin, y = yin, z = zin;
837
838 // Map the coordinates onto field map coordinates.
839 bool xmirr, ymirr, zmirr;
840 double rcoordinate, rotation;
841 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
842
843 if (m_warning) PrintWarning("WeightingPotential");
844
845 // Find the element that contains this point.
846 double t1, t2, t3, t4, jac[4][4], det;
847 const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
848 if (imap < 0) return 0.;
849 const Element& element = m_elements[imap];
850 if (m_debug) {
851 PrintElement("WeightingPotential", x, y, z, t1, t2, t3, t4, element, 10,
852 iw);
853 }
854 const Node& n0 = m_nodes[element.emap[0]];
855 const Node& n1 = m_nodes[element.emap[1]];
856 const Node& n2 = m_nodes[element.emap[2]];
857 const Node& n3 = m_nodes[element.emap[3]];
858 const Node& n4 = m_nodes[element.emap[4]];
859 const Node& n5 = m_nodes[element.emap[5]];
860 const Node& n6 = m_nodes[element.emap[6]];
861 const Node& n7 = m_nodes[element.emap[7]];
862 const Node& n8 = m_nodes[element.emap[8]];
863 const Node& n9 = m_nodes[element.emap[9]];
864 // Tetrahedral field
865
866 const auto it1 = std::upper_bound(m_wdtimes.cbegin(), m_wdtimes.cend(), t);
867 const auto it0 = std::prev(it1);
868
869 const double dt = t - *it0;
870 const unsigned int i0 = it0 - m_wdtimes.cbegin();
871 double dp0 =
872 n0.dw[iw][i0] * t1 * (2 * t1 - 1) + n1.dw[iw][i0] * t2 * (2 * t2 - 1) +
873 n2.dw[iw][i0] * t3 * (2 * t3 - 1) + n3.dw[iw][i0] * t4 * (2 * t4 - 1) +
874 4 * n4.dw[iw][i0] * t1 * t2 + 4 * n5.dw[iw][i0] * t1 * t3 +
875 4 * n6.dw[iw][i0] * t1 * t4 + 4 * n7.dw[iw][i0] * t2 * t3 +
876 4 * n8.dw[iw][i0] * t2 * t4 + 4 * n9.dw[iw][i0] * t3 * t4;
877
878 const unsigned int i1 = it1 - m_wdtimes.cbegin();
879
880 double dp1 =
881 n0.dw[iw][i1] * t1 * (2 * t1 - 1) + n1.dw[iw][i1] * t2 * (2 * t2 - 1) +
882 n2.dw[iw][i1] * t3 * (2 * t3 - 1) + n3.dw[iw][i1] * t4 * (2 * t4 - 1) +
883 4 * n4.dw[iw][i1] * t1 * t2 + 4 * n5.dw[iw][i1] * t1 * t3 +
884 4 * n6.dw[iw][i1] * t1 * t4 + 4 * n7.dw[iw][i1] * t2 * t3 +
885 4 * n8.dw[iw][i1] * t2 * t4 + 4 * n9.dw[iw][i1] * t3 * t4;
886
887 const double f1 = dt / (*it1 - *it0);
888 const double f0 = 1. - f1;
889
890 return f0 * dp0 + f1 * dp1;
891}
Garfield::ComponentFieldMap::MapCoordinates
void MapCoordinates(double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
Move (xpos, ypos, zpos) to field map coordinates.
Definition: ComponentFieldMap.cc:2021
Garfield::ComponentFieldMap::PrintWarning
void PrintWarning(const std::string &header)
Definition: ComponentFieldMap.cc:2375
Garfield::ComponentFieldMap::GetWeightingFieldIndex
size_t GetWeightingFieldIndex(const std::string &label) const
Definition: ComponentFieldMap.cc:2348
Garfield::ComponentFieldMap::m_wdtimes
std::vector< double > m_wdtimes
Definition: ComponentFieldMap.hh:144
Garfield::ComponentFieldMap::m_wfields
std::vector< std::string > m_wfields
Definition: ComponentFieldMap.hh:140
Garfield::ComponentFieldMap::FindElement13
int FindElement13(const double x, const double y, const double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det)
Find the element for a point in curved quadratic tetrahedra.
Definition: ComponentFieldMap.cc:337
Garfield::ComponentFieldMap::m_elements
std::vector< Element > m_elements
Definition: ComponentFieldMap.hh:113
Garfield::ComponentFieldMap::PrintElement
void PrintElement(const std::string &header, const double x, const double y, const double z, const double t1, const double t2, const double t3, const double t4, const Element &element, const unsigned int n, const int iw=-1) const
Definition: ComponentFieldMap.cc:2395
Garfield::Component::m_debug
bool m_debug
Switch on/off debugging messages.
Definition: Component.hh:341
Garfield::Component::m_ready
bool m_ready
Ready for use?
Definition: Component.hh:338
Element
Definition: neBEM.h:155

◆ ElectricField() [1/2]

void Garfield::ComponentComsol::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 378 of file ComponentComsol.cc.

381 {
382 // Copy the coordinates
383 double x = xin, y = yin, z = zin;
384
385 // Map the coordinates onto field map coordinates
386 bool xmirr, ymirr, zmirr;
387 double rcoordinate, rotation;
388 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
389
390 // Initial values
391 ex = ey = ez = volt = 0.;
392 status = 0;
393 m = nullptr;
394
395 // Do not proceed if not properly initialised.
396 if (!m_ready) {
397 status = -10;
398 PrintNotReady("ElectricField");
399 return;
400 }
401
402 if (m_warning) PrintWarning("ElectricField");
403
404 // Find the element that contains this point
405 double t1, t2, t3, t4, jac[4][4], det;
406 const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
407 if (imap < 0) {
408 if (m_debug) {
409 std::cout << m_className << "::ElectricField:\n"
410 << " Point (" << x << ", " << y << ", " << z
411 << " not in the mesh.\n";
412 }
413 status = -6;
414 return;
415 }
416
417 const Element& element = m_elements[imap];
418 if (m_debug) {
419 PrintElement("ElectricField", x, y, z, t1, t2, t3, t4, element, 10);
420 }
421 const Node& n0 = m_nodes[element.emap[0]];
422 const Node& n1 = m_nodes[element.emap[1]];
423 const Node& n2 = m_nodes[element.emap[2]];
424 const Node& n3 = m_nodes[element.emap[3]];
425 const Node& n4 = m_nodes[element.emap[4]];
426 const Node& n5 = m_nodes[element.emap[5]];
427 const Node& n6 = m_nodes[element.emap[6]];
428 const Node& n7 = m_nodes[element.emap[7]];
429 const Node& n8 = m_nodes[element.emap[8]];
430 const Node& n9 = m_nodes[element.emap[9]];
431 // Tetrahedral field
432 volt = n0.v * t1 * (2 * t1 - 1) + n1.v * t2 * (2 * t2 - 1) +
433 n2.v * t3 * (2 * t3 - 1) + n3.v * t4 * (2 * t4 - 1) +
434 4 * n4.v * t1 * t2 + 4 * n5.v * t1 * t3 + 4 * n6.v * t1 * t4 +
435 4 * n7.v * t2 * t3 + 4 * n8.v * t2 * t4 + 4 * n9.v * t3 * t4;
436 ex = -(n0.v * (4 * t1 - 1) * jac[0][1] + n1.v * (4 * t2 - 1) * jac[1][1] +
437 n2.v * (4 * t3 - 1) * jac[2][1] + n3.v * (4 * t4 - 1) * jac[3][1] +
438 n4.v * (4 * t2 * jac[0][1] + 4 * t1 * jac[1][1]) +
439 n5.v * (4 * t3 * jac[0][1] + 4 * t1 * jac[2][1]) +
440 n6.v * (4 * t4 * jac[0][1] + 4 * t1 * jac[3][1]) +
441 n7.v * (4 * t3 * jac[1][1] + 4 * t2 * jac[2][1]) +
442 n8.v * (4 * t4 * jac[1][1] + 4 * t2 * jac[3][1]) +
443 n9.v * (4 * t4 * jac[2][1] + 4 * t3 * jac[3][1])) /
444 det;
445 ey = -(n0.v * (4 * t1 - 1) * jac[0][2] + n1.v * (4 * t2 - 1) * jac[1][2] +
446 n2.v * (4 * t3 - 1) * jac[2][2] + n3.v * (4 * t4 - 1) * jac[3][2] +
447 n4.v * (4 * t2 * jac[0][2] + 4 * t1 * jac[1][2]) +
448 n5.v * (4 * t3 * jac[0][2] + 4 * t1 * jac[2][2]) +
449 n6.v * (4 * t4 * jac[0][2] + 4 * t1 * jac[3][2]) +
450 n7.v * (4 * t3 * jac[1][2] + 4 * t2 * jac[2][2]) +
451 n8.v * (4 * t4 * jac[1][2] + 4 * t2 * jac[3][2]) +
452 n9.v * (4 * t4 * jac[2][2] + 4 * t3 * jac[3][2])) /
453 det;
454 ez = -(n0.v * (4 * t1 - 1) * jac[0][3] + n1.v * (4 * t2 - 1) * jac[1][3] +
455 n2.v * (4 * t3 - 1) * jac[2][3] + n3.v * (4 * t4 - 1) * jac[3][3] +
456 n4.v * (4 * t2 * jac[0][3] + 4 * t1 * jac[1][3]) +
457 n5.v * (4 * t3 * jac[0][3] + 4 * t1 * jac[2][3]) +
458 n6.v * (4 * t4 * jac[0][3] + 4 * t1 * jac[3][3]) +
459 n7.v * (4 * t3 * jac[1][3] + 4 * t2 * jac[2][3]) +
460 n8.v * (4 * t4 * jac[1][3] + 4 * t2 * jac[3][3]) +
461 n9.v * (4 * t4 * jac[2][3] + 4 * t3 * jac[3][3])) /
462 det;
463
464 // Transform field to global coordinates
465 UnmapFields(ex, ey, ez, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
466
467 // Drift medium?
468 if (m_debug) {
469 std::cout << m_className << "::ElectricField:\n"
470 << " Material " << element.matmap << ", drift flag "
471 << m_materials[element.matmap].driftmedium << "\n";
472 }
473 m = m_materials[element.matmap].medium;
474 status = -5;
475 if (m_materials[element.matmap].driftmedium) {
476 if (m && m->IsDriftable()) status = 0;
477 }
478}
Garfield::ComponentFieldMap::m_materials
std::vector< Material > m_materials
Definition: ComponentFieldMap.hh:138
Garfield::ComponentFieldMap::UnmapFields
void UnmapFields(double &ex, double &ey, double &ez, double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
Move (ex, ey, ez) to global coordinates.
Definition: ComponentFieldMap.cc:2138
Garfield::ComponentFieldMap::PrintNotReady
void PrintNotReady(const std::string &header) const
Definition: ComponentFieldMap.cc:2383
Garfield::Component::m_className
std::string m_className
Class name.
Definition: Component.hh:329

◆ ElectricField() [2/2]

void Garfield::ComponentComsol::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 371 of file ComponentComsol.cc.

373 {
374 double v = 0.;
375 ElectricField(x, y, z, ex, ey, ez, v, m, status);
376}
Garfield::ComponentComsol::ElectricField
void ElectricField(const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status) override
Definition: ComponentComsol.cc:371

Referenced by ElectricField().

◆ GetAspectRatio()

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

Implements Garfield::ComponentFieldMap.

Definition at line 684 of file ComponentComsol.cc.

685 {
686 if (i >= m_elements.size()) {
687 dmin = dmax = 0.;
688 return;
689 }
690
691 const Element& element = m_elements[i];
692 const int np = 4;
693 // Loop over all pairs of vertices.
694 for (int j = 0; j < np - 1; ++j) {
695 const Node& nj = m_nodes[element.emap[j]];
696 for (int k = j + 1; k < np; ++k) {
697 const Node& nk = m_nodes[element.emap[k]];
698 // Compute distance.
699 const double dx = nj.x - nk.x;
700 const double dy = nj.y - nk.y;
701 const double dz = nj.z - nk.z;
702 const double dist = sqrt(dx * dx + dy * dy + dz * dz);
703 if (k == 1) {
704 dmin = dmax = dist;
705 } else {
706 if (dist < dmin) dmin = dist;
707 if (dist > dmax) dmax = dist;
708 }
709 }
710 }
711}
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314

◆ GetElementVolume()

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

Implements Garfield::ComponentFieldMap.

Definition at line 663 of file ComponentComsol.cc.

663 {
664 if (i >= m_elements.size()) return 0.;
665 const Element& element = m_elements[i];
666 const Node& n0 = m_nodes[element.emap[0]];
667 const Node& n1 = m_nodes[element.emap[1]];
668 const Node& n2 = m_nodes[element.emap[2]];
669 const Node& n3 = m_nodes[element.emap[3]];
670
671 // Uses formula V = |a (dot) b x c|/6
672 // with a => "3", b => "1", c => "2" and origin "0"
673 const double vol =
674 fabs((n3.x - n0.x) *
675 ((n1.y - n0.y) * (n2.z - n0.z) - (n2.y - n0.y) * (n1.z - n0.z)) +
676 (n3.y - n0.y) *
677 ((n1.z - n0.z) * (n2.x - n0.x) - (n2.z - n0.z) * (n1.x - n0.x)) +
678 (n3.z - n0.z) * ((n1.x - n0.x) * (n2.y - n0.y) -
679 (n3.x - n0.x) * (n1.y - n0.y))) /
680 6.;
681 return vol;
682}
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615

◆ GetMedium()

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

619 {
620 // Copy the coordinates
621 double x = xin, y = yin, z = zin;
622
623 // Map the coordinates onto field map coordinates
624 bool xmirr, ymirr, zmirr;
625 double rcoordinate, rotation;
626 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
627
628 // Do not proceed if not properly initialised.
629 if (!m_ready) {
630 PrintNotReady("GetMedium");
631 return nullptr;
632 }
633 if (m_warning) PrintWarning("GetMedium");
634
635 // Find the element that contains this point
636 double t1, t2, t3, t4, jac[4][4], det;
637 const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
638 if (imap < 0) {
639 if (m_debug) {
640 std::cout << m_className << "::GetMedium:\n"
641 << " Point (" << x << ", " << y << ", " << z
642 << ") not in the mesh.\n";
643 }
644 return nullptr;
645 }
646 const Element& element = m_elements[imap];
647 if (element.matmap >= m_materials.size()) {
648 if (m_debug) {
649 std::cerr << m_className << "::GetMedium:\n"
650 << " Point (" << x << ", " << y << ", " << z
651 << ") has out of range material number " << imap << ".\n";
652 }
653 return nullptr;
654 }
655
656 if (m_debug) {
657 PrintElement("GetMedium", x, y, z, t1, t2, t3, t4, element, 10);
658 }
659
660 return m_materials[element.matmap].medium;
661}

◆ Initialise()

bool Garfield::ComponentComsol::Initialise ( const std::string &  header = "mesh.mphtxt",
const std::string &  mplist = "dielectrics.dat",
const std::string &  field = "field.txt",
const std::string &  unit = "m" 
)

Import a field map.

Parameters
headername of the file containing the list of nodes
mplistname of the file containing the material properties
fieldname of the file containing the potentials at the nodes
unitlength unit

Definition at line 62 of file ComponentComsol.cc.

65 {
66 Reset();
67
68 // Get the conversion factor to be applied to the coordinates.
69 m_unit = ScalingFactor(unit);
70 if (m_unit <= 0.) {
71 std::cerr << m_className << "::Initialise:\n Unknown length unit "
72 << unit << ". Will use default (m).\n";
73 m_unit = 100.;
74 }
75 // Open the materials file.
76 m_materials.clear();
77 std::ifstream fmplist;
78 fmplist.open(mplist.c_str(), std::ios::in);
79 if (fmplist.fail()) {
80 PrintCouldNotOpen("Initialise", mplist);
81 return false;
82 }
83 unsigned int nMaterials;
84 fmplist >> nMaterials;
85 for (unsigned int i = 0; i < nMaterials; ++i) {
86 Material material;
87 material.driftmedium = false;
88 material.medium = nullptr;
89 material.ohm = -1;
90 fmplist >> material.eps;
91 m_materials.push_back(std::move(material));
92 }
93 if (m_materials.empty()) {
94 // Add default material
95 Material material;
96 material.driftmedium = false;
97 material.medium = nullptr;
98 material.eps = material.ohm = -1;
99 m_materials.push_back(std::move(material));
100 nMaterials = 1;
101 }
102 std::map<int, int> domain2material;
103 int d2msize;
104 fmplist >> d2msize;
105 for (int i = 0; i < d2msize; ++i) {
106 int domain;
107 fmplist >> domain;
108 fmplist >> domain2material[domain];
109 }
110 fmplist.close();
111
112 // Find lowest epsilon, check for eps = 0, set default drift medium.
113 if (!SetDefaultDriftMedium()) return false;
114
115 m_nodes.clear();
116 std::ifstream fmesh;
117 fmesh.open(mesh.c_str(), std::ios::in);
118 if (fmesh.fail()) {
119 PrintCouldNotOpen("Initialise", mesh);
120 return false;
121 }
122
123 std::string line;
124 do {
125 if (!std::getline(fmesh, line)) {
126 std::cerr << m_className << "::Initialise:\n"
127 << " Could not read number of nodes from " << mesh << ".\n";
128 fmesh.close();
129 return false;
130 }
131 } while (!ends_with(line, "# number of mesh points") &&
132 !ends_with(line, "# number of mesh vertices"));
133 const int nNodes = readInt(line);
134
135 std::cout << m_className << "::Initialise: " << nNodes << " nodes.\n";
136 do {
137 if (!std::getline(fmesh, line)) {
138 std::cerr << m_className << "::Initialise:\n"
139 << " Error parsing " << mesh << ".\n";
140 fmesh.close();
141 return false;
142 }
143 } while (line.find("# Mesh point coordinates") == std::string::npos &&
144 line.find("# Mesh vertex coordinates") == std::string::npos);
145 for (int i = 0; i < nNodes; ++i) {
146 Node newNode;
147 fmesh >> newNode.x >> newNode.y >> newNode.z;
148 newNode.x *= m_unit;
149 newNode.y *= m_unit;
150 newNode.z *= m_unit;
151 m_nodes.push_back(std::move(newNode));
152 }
153
154 do {
155 if (!std::getline(fmesh, line)) {
156 std::cerr << m_className << "::Initialise:\n"
157 << " Error parsing " << mesh << ".\n";
158 fmesh.close();
159 return false;
160 }
161 } while (line.find("4 tet2 # type name") == std::string::npos);
162 do {
163 if (!std::getline(fmesh, line)) {
164 std::cerr << m_className << "::Initialise:\n"
165 << " Error parsing " << mesh << ".\n";
166 fmesh.close();
167 return false;
168 }
169 } while (!ends_with(line, "# number of elements"));
170 const int nElements = readInt(line);
171 m_elements.clear();
172 std::cout << m_className << "::Initialise: " << nElements << " elements.\n";
173 std::getline(fmesh, line);
174 // Elements 6 & 7 are swapped due to differences in COMSOL and ANSYS
175 // representation
176 int perm[10] = {0, 1, 2, 3, 4, 5, 7, 6, 8, 9};
177 for (int i = 0; i < nElements; ++i) {
178 Element newElement;
179 newElement.degenerate = false;
180 for (int j = 0; j < 10; ++j) {
181 fmesh >> newElement.emap[perm[j]];
182 }
183 m_elements.push_back(std::move(newElement));
184 }
185
186 do {
187 if (!std::getline(fmesh, line)) {
188 std::cerr << m_className << "::Initialise:\n"
189 << " Error parsing " << mesh << ".\n";
190 fmesh.close();
191 return false;
192 }
193 } while (line.find("# Geometric entity indices") == std::string::npos);
194 for (auto& element : m_elements) {
195 int domain;
196 fmesh >> domain;
197 element.matmap = domain2material.count(domain) ? domain2material[domain]
198 : nMaterials - 1;
199 }
200 fmesh.close();
201
202 std::ifstream ffield;
203 ffield.open(field.c_str(), std::ios::in);
204 if (ffield.fail()) {
205 PrintCouldNotOpen("Initialise", field);
206 return false;
207 }
208
209 const std::string hdr1 =
210 "% x y z "
211 " V (V)";
212
213 const std::string hdr2 =
214 "% x y z V (V)";
215 do {
216 if (!std::getline(ffield, line)) {
217 std::cerr << m_className << "::Initialise:\n"
218 << " Error parsing " << field << ".\n";
219 ffield.close();
220 return false;
221 }
222 } while (line.find(hdr1) == std::string::npos &&
223 line.find(hdr2) == std::string::npos);
224 std::istringstream sline(line);
225 std::string token;
226 sline >> token; // %
227 sline >> token; // x
228 sline >> token; // y
229 sline >> token; // z
230 sline >> token; // V
231 sline >> token; // (V)
232 m_wfields.clear();
233 m_wfieldsOk.clear();
234 while (sline >> token) {
235 std::cout << m_className << "::Initialise:\n"
236 << " Reading data for weighting field " << token << ".\n";
237 m_wfields.push_back(token);
238 m_wfieldsOk.push_back(true);
239 sline >> token; // (V)
240 }
241 const size_t nWeightingFields = m_wfields.size();
242
243 const unsigned int nPrint =
244 std::pow(10, static_cast<unsigned int>(
245 std::max(std::floor(std::log10(nNodes)) - 1, 1.)));
246 std::cout << m_className << "::Initialise: Reading potentials.\n";
247 PrintProgress(0.);
248 // Build a k-d tree from the node coordinates.
249 std::vector<std::vector<double> > points;
250 for (const auto& node : m_nodes) {
251 std::vector<double> point = {node.x, node.y, node.z};
252 points.push_back(std::move(point));
253 }
254 KDTree kdtree(points);
255 std::vector<bool> used(nNodes, false);
256 for (int i = 0; i < nNodes; ++i) {
257 double x, y, z, v;
258 ffield >> x >> y >> z >> v;
259 x *= m_unit;
260 y *= m_unit;
261 z *= m_unit;
262 std::vector<double> w;
263 for (size_t j = 0; j < nWeightingFields; ++j) {
264 double p;
265 ffield >> p;
266 w.push_back(p);
267 }
268 const std::vector<double> pt = {x, y, z};
269 std::vector<KDTreeResult> res;
270 kdtree.n_nearest(pt, 1, res);
271 if (res.empty()) {
272 std::cerr << std::endl
273 << m_className << "::Initialise:\n"
274 << " Could not find a matching mesh node for point (" << x
275 << ", " << y << ", " << z << ")\n.";
276 ffield.close();
277 return false;
278 }
279 const size_t k = res[0].idx;
280 used[k] = true;
281 m_nodes[k].v = v;
282 m_nodes[k].w = w;
283 if ((i + 1) % nPrint == 0) PrintProgress(double(i + 1) / nNodes);
284 }
285 PrintProgress(1.);
286 ffield.close();
287 auto nMissing = std::count(used.begin(), used.end(), false);
288 if (nMissing > 0) {
289 std::cerr << std::endl
290 << m_className << "::Initialise:\n"
291 << " Missing potentials for " << nMissing << " nodes.\n";
292 return false;
293 }
294
295 m_ready = true;
296 Prepare();
297 std::cout << std::endl << m_className << "::Initialise: Done.\n";
298 return true;
299}
Garfield::ComponentFieldMap::SetDefaultDriftMedium
bool SetDefaultDriftMedium()
Find lowest epsilon, check for eps = 0, set default drift media flags.
Definition: ComponentFieldMap.cc:162
Garfield::ComponentFieldMap::ScalingFactor
static double ScalingFactor(std::string unit)
Definition: ComponentFieldMap.cc:2213
Garfield::ComponentFieldMap::m_wfieldsOk
std::vector< bool > m_wfieldsOk
Definition: ComponentFieldMap.hh:141
Garfield::ComponentFieldMap::PrintCouldNotOpen
void PrintCouldNotOpen(const std::string &header, const std::string &filename) const
Definition: ComponentFieldMap.cc:2388
Garfield::ComponentFieldMap::Reset
void Reset() override
Reset the component.
Definition: ComponentFieldMap.cc:1695

Referenced by ComponentComsol().

◆ SetDelayedWeightingPotential()

bool Garfield::ComponentComsol::SetDelayedWeightingPotential ( const std::string &  file,
const std::string &  label 
)

Import the time-dependent weighting field maps.

Definition at line 713 of file ComponentComsol.cc.

714 {
715 if (!m_ready) {
716 std::cerr << m_className << "::SetDelayedWeightingField:\n"
717 << " No valid field map is present.\n"
718 << " Weighting fields cannot be added.\n";
719 return false;
720 }
721
722 if (!GetTimeInterval(field)) return false;
723
724 if (!m_timeset) {
725 std::cerr << m_className << "::SetDelayedWeightingField:\n"
726 << " No valid times slices of potential set.\n"
727 << " Please add the time slices.\n";
728 return false;
729 }
730
731 const int T = m_wdtimes.size();
732
733 double x, y, z;
734
735 // Open the voltage list.
736 std::ifstream ffield;
737 ffield.open(field.c_str(), std::ios::in);
738
739 // Check if a weighting field with the same label already exists.
740 const size_t iw = GetOrCreateWeightingFieldIndex(label);
741
742 if (iw + 1 != m_wfields.size()) {
743 std::cout << m_className << "::SetDelayedWeightingField:\n"
744 << " Replacing existing weighting field " << label << ".\n";
745 }
746
747 m_wfieldsOk[iw] = false; //????
748
749 // Build a k-d tree from the node coordinates.
750
751 std::vector<std::vector<double> > points;
752 for (const auto& node : m_nodes) {
753 std::vector<double> point = {node.x, node.y, node.z};
754 points.push_back(std::move(point));
755 }
756
757 KDTree kdtree(points);
758
759 std::string line;
760 int Linecount = 1;
761 const int nNodes = m_nodes.size();
762 const unsigned int nPrint =
763 std::pow(10, static_cast<unsigned int>(
764 std::max(std::floor(std::log10(nNodes)) - 1, 1.)));
765 std::cout << m_className << "::SetDelayedWeightingField:\n"
766 << " Reading weighting potentials.\n";
767 PrintProgress(0.);
768
769 while (std::getline(ffield, line)) {
770 // Skip empty lines.
771 if (line.empty()) continue;
772 // Skip lines that are not comments.
773 if (isComment(line)) continue;
774
775 std::vector<double> pvect;
776
777 std::istringstream data;
778 data.str(line);
779 data >> x >> y >> z;
780 x *= m_unit;
781 y *= m_unit;
782 z *= m_unit;
783 const std::vector<double> pt = {x, y, z};
784 std::vector<KDTreeResult> res;
785 kdtree.n_nearest(pt, 1, res);
786
787 if (res.empty()) {
788 std::cerr << m_className << "::SetDelayedWeightingField:\n"
789 << " Could not find a matching mesh node for point (" << x
790 << ", " << y << ", " << z << ")\n.";
791 ffield.close();
792 return false;
793 }
794
795 double pholder = 0.;
796 double pholder0 = 0.;
797 for (int i = 0; i < T; i++) {
798 data >> pholder;
799 if (i == 0) pholder0 = pholder;
800 pvect.push_back(pholder - pholder0);
801 }
802 const size_t k = res[0].idx;
803 m_nodes[k].dw[iw] = pvect;
804 m_nodes[k].w[iw] = pholder0;
805
806 if ((Linecount + 1) % nPrint == 0)
807 PrintProgress(double(Linecount + 1) / nNodes);
808 Linecount++;
809 }
810
811 PrintProgress(1.);
812 std::cout << std::endl
813 << m_className << "::SetDelayedWeightingField: Done.\n";
814 ffield.close();
815 return true;
816}
Garfield::ComponentFieldMap::GetOrCreateWeightingFieldIndex
size_t GetOrCreateWeightingFieldIndex(const std::string &label)
Definition: ComponentFieldMap.cc:2357

◆ SetTimeInterval()

void Garfield::ComponentComsol::SetTimeInterval ( const double  mint,
const double  maxt,
const double  stept 
)

Set the time interval of the time-dependent weighting field.

Definition at line 893 of file ComponentComsol.cc.

894 {
895 std::cout << std::endl
896 << m_className
897 << "::SetTimeInterval: Overwriting time interval of weighting "
898 "potential.\n";
899
900 if (m_wdtimes.empty()) {
901 double t = mint;
902 while (t <= maxt) {
903 m_wdtimes.push_back(t);
904 t += stept;
905 }
906 }
907 m_timeset = true;
908
909 std::cout << std::endl
910 << m_className
911 << "::SetTimeInterval: Time of weighting potential set for t in ["
912 << mint << "," << maxt << "].\n";
913}

◆ SetWeightingField()

bool Garfield::ComponentComsol::SetWeightingField ( const std::string &  file,
const std::string &  label 
)

Import the time-independent weighting field maps.

Definition at line 301 of file ComponentComsol.cc.

302 {
303 if (!m_ready) {
304 std::cerr << m_className << "::SetWeightingField:\n"
305 << " No valid field map is present.\n"
306 << " Weighting fields cannot be added.\n";
307 return false;
308 }
309
310 // Open the voltage list.
311 std::ifstream ffield;
312 ffield.open(field.c_str(), std::ios::in);
313 if (ffield.fail()) {
314 PrintCouldNotOpen("SetWeightingField", field);
315 return false;
316 }
317
318 // Check if a weighting field with the same label already exists.
319 const size_t iw = GetOrCreateWeightingFieldIndex(label);
320 if (iw + 1 != m_wfields.size()) {
321 std::cout << m_className << "::SetWeightingField:\n"
322 << " Replacing existing weighting field " << label << ".\n";
323 }
324 m_wfieldsOk[iw] = false;
325
326 // Build a k-d tree from the node coordinates.
327 std::vector<std::vector<double> > points;
328 for (const auto& node : m_nodes) {
329 std::vector<double> point = {node.x, node.y, node.z};
330 points.push_back(std::move(point));
331 }
332 KDTree kdtree(points);
333
334 const std::string hdr =
335 "% x y z "
336 " V (V)";
337 std::string line;
338 do {
339 if (!std::getline(ffield, line)) {
340 std::cerr << m_className << "::SetWeightingField:\n"
341 << " Error parsing " << field << ".\n";
342 ffield.close();
343 return false;
344 }
345 } while (line.find(hdr) == std::string::npos);
346 const int nNodes = m_nodes.size();
347 for (int i = 0; i < nNodes; ++i) {
348 double x, y, z, v;
349 ffield >> x >> y >> z >> v;
350 x *= m_unit;
351 y *= m_unit;
352 z *= m_unit;
353 // Find the closest mesh node.
354 const std::vector<double> pt = {x, y, z};
355 std::vector<KDTreeResult> res;
356 kdtree.n_nearest(pt, 1, res);
357 if (res.empty()) {
358 std::cerr << m_className << "::SetWeightingField:\n"
359 << " Could not find a matching mesh node for point (" << x
360 << ", " << y << ", " << z << ")\n.";
361 ffield.close();
362 return false;
363 }
364 const size_t k = res[0].idx;
365 m_nodes[k].w[iw] = v;
366 }
367 ffield.close();
368 return true;
369}

◆ UpdatePeriodicity()

void Garfield::ComponentComsol::UpdatePeriodicity ( )
inlineoverrideprotectedvirtual

Verify periodicities.

Implements Garfield::Component.

Definition at line 57 of file ComponentComsol.hh.

57{ UpdatePeriodicityCommon(); }

◆ WeightingField()

void Garfield::ComponentComsol::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 480 of file ComponentComsol.cc.

482 {
483 // Initial values
484 wx = wy = wz = 0;
485
486 // Do not proceed if not properly initialised.
487 if (!m_ready) return;
488
489 // Look for the label.
490 const size_t iw = GetWeightingFieldIndex(label);
491 // Do not proceed if the requested weighting field does not exist.
492 if (iw == m_wfields.size()) return;
493 // Check if the weighting field is properly initialised.
494 if (!m_wfieldsOk[iw]) return;
495
496 // Copy the coordinates.
497 double x = xin, y = yin, z = zin;
498
499 // Map the coordinates onto field map coordinates
500 bool xmirr, ymirr, zmirr;
501 double rcoordinate, rotation;
502 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
503
504 if (m_warning) PrintWarning("WeightingField");
505
506 // Find the element that contains this point.
507 double t1, t2, t3, t4, jac[4][4], det;
508 const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
509 // Check if the point is in the mesh.
510 if (imap < 0) return;
511
512 const Element& element = m_elements[imap];
513 if (m_debug) {
514 PrintElement("WeightingField", x, y, z, t1, t2, t3, t4, element, 10, iw);
515 }
516 const Node& n0 = m_nodes[element.emap[0]];
517 const Node& n1 = m_nodes[element.emap[1]];
518 const Node& n2 = m_nodes[element.emap[2]];
519 const Node& n3 = m_nodes[element.emap[3]];
520 const Node& n4 = m_nodes[element.emap[4]];
521 const Node& n5 = m_nodes[element.emap[5]];
522 const Node& n6 = m_nodes[element.emap[6]];
523 const Node& n7 = m_nodes[element.emap[7]];
524 const Node& n8 = m_nodes[element.emap[8]];
525 const Node& n9 = m_nodes[element.emap[9]];
526 // Tetrahedral field
527 wx = -(n0.w[iw] * (4 * t1 - 1) * jac[0][1] +
528 n1.w[iw] * (4 * t2 - 1) * jac[1][1] +
529 n2.w[iw] * (4 * t3 - 1) * jac[2][1] +
530 n3.w[iw] * (4 * t4 - 1) * jac[3][1] +
531 n4.w[iw] * (4 * t2 * jac[0][1] + 4 * t1 * jac[1][1]) +
532 n5.w[iw] * (4 * t3 * jac[0][1] + 4 * t1 * jac[2][1]) +
533 n6.w[iw] * (4 * t4 * jac[0][1] + 4 * t1 * jac[3][1]) +
534 n7.w[iw] * (4 * t3 * jac[1][1] + 4 * t2 * jac[2][1]) +
535 n8.w[iw] * (4 * t4 * jac[1][1] + 4 * t2 * jac[3][1]) +
536 n9.w[iw] * (4 * t4 * jac[2][1] + 4 * t3 * jac[3][1])) /
537 det;
538
539 wy = -(n0.w[iw] * (4 * t1 - 1) * jac[0][2] +
540 n1.w[iw] * (4 * t2 - 1) * jac[1][2] +
541 n2.w[iw] * (4 * t3 - 1) * jac[2][2] +
542 n3.w[iw] * (4 * t4 - 1) * jac[3][2] +
543 n4.w[iw] * (4 * t2 * jac[0][2] + 4 * t1 * jac[1][2]) +
544 n5.w[iw] * (4 * t3 * jac[0][2] + 4 * t1 * jac[2][2]) +
545 n6.w[iw] * (4 * t4 * jac[0][2] + 4 * t1 * jac[3][2]) +
546 n7.w[iw] * (4 * t3 * jac[1][2] + 4 * t2 * jac[2][2]) +
547 n8.w[iw] * (4 * t4 * jac[1][2] + 4 * t2 * jac[3][2]) +
548 n9.w[iw] * (4 * t4 * jac[2][2] + 4 * t3 * jac[3][2])) /
549 det;
550
551 wz = -(n0.w[iw] * (4 * t1 - 1) * jac[0][3] +
552 n1.w[iw] * (4 * t2 - 1) * jac[1][3] +
553 n2.w[iw] * (4 * t3 - 1) * jac[2][3] +
554 n3.w[iw] * (4 * t4 - 1) * jac[3][3] +
555 n4.w[iw] * (4 * t2 * jac[0][3] + 4 * t1 * jac[1][3]) +
556 n5.w[iw] * (4 * t3 * jac[0][3] + 4 * t1 * jac[2][3]) +
557 n6.w[iw] * (4 * t4 * jac[0][3] + 4 * t1 * jac[3][3]) +
558 n7.w[iw] * (4 * t3 * jac[1][3] + 4 * t2 * jac[2][3]) +
559 n8.w[iw] * (4 * t4 * jac[1][3] + 4 * t2 * jac[3][3]) +
560 n9.w[iw] * (4 * t4 * jac[2][3] + 4 * t3 * jac[3][3])) /
561 det;
562
563 // Transform field to global coordinates
564 UnmapFields(wx, wy, wz, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
565}

◆ WeightingPotential()

double Garfield::ComponentComsol::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 567 of file ComponentComsol.cc.

569 {
570 // Do not proceed if not properly initialised.
571 if (!m_ready) return 0.;
572
573 // Look for the label.
574 const size_t iw = GetWeightingFieldIndex(label);
575 // Do not proceed if the requested weighting field does not exist.
576 if (iw == m_wfields.size()) return 0.;
577 // Check if the weighting field is properly initialised.
578 // if (!m_wfieldsOk[iw]) return 0.;
579
580 // Copy the coordinates.
581 double x = xin, y = yin, z = zin;
582
583 // Map the coordinates onto field map coordinates.
584 bool xmirr, ymirr, zmirr;
585 double rcoordinate, rotation;
586 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
587
588 if (m_warning) PrintWarning("WeightingPotential");
589
590 // Find the element that contains this point.
591 double t1, t2, t3, t4, jac[4][4], det;
592 const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
593 if (imap < 0) return 0.;
594
595 const Element& element = m_elements[imap];
596 if (m_debug) {
597 PrintElement("WeightingPotential", x, y, z, t1, t2, t3, t4, element, 10,
598 iw);
599 }
600 const Node& n0 = m_nodes[element.emap[0]];
601 const Node& n1 = m_nodes[element.emap[1]];
602 const Node& n2 = m_nodes[element.emap[2]];
603 const Node& n3 = m_nodes[element.emap[3]];
604 const Node& n4 = m_nodes[element.emap[4]];
605 const Node& n5 = m_nodes[element.emap[5]];
606 const Node& n6 = m_nodes[element.emap[6]];
607 const Node& n7 = m_nodes[element.emap[7]];
608 const Node& n8 = m_nodes[element.emap[8]];
609 const Node& n9 = m_nodes[element.emap[9]];
610 // Tetrahedral field
611 return n0.w[iw] * t1 * (2 * t1 - 1) + n1.w[iw] * t2 * (2 * t2 - 1) +
612 n2.w[iw] * t3 * (2 * t3 - 1) + n3.w[iw] * t4 * (2 * t4 - 1) +
613 4 * n4.w[iw] * t1 * t2 + 4 * n5.w[iw] * t1 * t3 +
614 4 * n6.w[iw] * t1 * t4 + 4 * n7.w[iw] * t2 * t3 +
615 4 * n8.w[iw] * t2 * t4 + 4 * n9.w[iw] * t3 * t4;
616}
The documentation for this class was generated from the following files: