Garfield++ 3.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 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:

Classes

struct  nodeCmp
 

Public Member Functions

 ComponentComsol ()
 Default constructor.
 
 ComponentComsol (std::string mesh, std::string mplist, std::string field)
 
 ~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
 
Medium * GetMedium (const double x, const double y, const double z) override
 Get the medium at a given location (x, y, z).
 
bool Initialise (std::string header="mesh.mphtxt", std::string mplist="dielectrics.dat", std::string field="field.txt")
 
bool SetWeightingField (std::string file, std::string label)
 
- Public Member Functions inherited from Garfield::ComponentFieldMap
 ComponentFieldMap ()
 Constructor.
 
virtual ~ComponentFieldMap ()
 Destructor.
 
virtual void SetRange ()
 Calculate x, y, z, V and angular ranges.
 
void PrintRange ()
 Show x, y, z, V and angular ranges.
 
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 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.
 
unsigned int GetNumberOfMaterials () const
 Return the number of materials in the field map.
 
double GetPermittivity (const unsigned int imat) const
 Return the permittivity of a field map material.
 
double GetConductivity (const unsigned int imat) const
 Return the conductivity of a field map material.
 
void SetMedium (const unsigned int imat, Medium *medium)
 Associate a field map material with a Medium class.
 
Medium * GetMedium (const unsigned int i) const
 Return the Medium associated to a field map material.
 
unsigned int GetNumberOfMedia () const
 
int GetNumberOfElements () const
 Return the number of mesh elements.
 
bool GetElement (const unsigned int i, double &vol, double &dmin, double &dmax)
 Return the volume and aspect ratio of a mesh element.
 
void EnableCheckMapIndices ()
 
void DisableCheckMapIndices ()
 
void EnableDeleteBackgroundElements (const bool on=true)
 Option to eliminate mesh elements in conductors (default: on).
 
void EnableTetrahedralTreeForElementSearch (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::ComponentBase
 ComponentBase ()
 Constructor.
 
virtual ~ComponentBase ()
 Destructor.
 
virtual void SetGeometry (GeometryBase *geo)
 Define the geometry.
 
virtual void Clear ()
 Reset.
 
virtual Medium * GetMedium (const double x, const double y, const double z)
 Get the medium at a given location (x, y, z).
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status)=0
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status)=0
 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 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.
 
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 IntegrateFlux (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)
 
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 DisablePeriodicityX ()
 
void EnablePeriodicityY (const bool on=true)
 Enable simple periodicity in the $y$ direction.
 
void DisablePeriodicityY ()
 
void EnablePeriodicityZ (const bool on=true)
 Enable simple periodicity in the $z$ direction.
 
void DisablePeriodicityZ ()
 
void EnableMirrorPeriodicityX (const bool on=true)
 Enable mirror periodicity in the $x$ direction.
 
void DisableMirrorPeriodicityX ()
 
void EnableMirrorPeriodicityY (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void DisableMirrorPeriodicityY ()
 
void EnableMirrorPeriodicityZ (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void DisableMirrorPeriodicityZ ()
 
void EnableAxialPeriodicityX (const bool on=true)
 Enable axial periodicity in the $x$ direction.
 
void DisableAxialPeriodicityX ()
 
void EnableAxialPeriodicityY (const bool on=true)
 Enable axial periodicity in the $y$ direction.
 
void DisableAxialPeriodicityY ()
 
void EnableAxialPeriodicityZ (const bool on=true)
 Enable axial periodicity in the $z$ direction.
 
void DisableAxialPeriodicityZ ()
 
void EnableRotationSymmetryX (const bool on=true)
 Enable rotation symmetry around the $x$ axis.
 
void DisableRotationSymmetryX ()
 
void EnableRotationSymmetryY (const bool on=true)
 Enable rotation symmetry around the $y$ axis.
 
void DisableRotationSymmetryY ()
 
void EnableRotationSymmetryZ (const bool on=true)
 Enable rotation symmetry around the $z$ axis.
 
void DisableRotationSymmetryZ ()
 
void EnableDebugging ()
 Switch on debugging messages.
 
void DisableDebugging ()
 Switch off debugging messages.
 
void ActivateTraps ()
 Request trapping to be taken care of by the component (for TCAD).
 
void DeactivateTraps ()
 
bool IsTrapActive ()
 
void ActivateVelocityMap ()
 Request velocity to be taken care of by the component (for TCAD).
 
void DectivateVelocityMap ()
 
bool IsVelocityActive ()
 
virtual bool ElectronAttachment (const double, const double, const double, double &eta)
 Get the electron attachment coefficient.
 
virtual bool HoleAttachment (const double, const double, const double, double &eta)
 Get the hole attachment coefficient.
 
virtual void ElectronVelocity (const double, const double, const double, double &vx, double &vy, double &vz, Medium *&, int &status)
 Get the electron drift velocity.
 
virtual void HoleVelocity (const double, const double, const double, double &vx, double &vy, double &vz, Medium *&, int &status)
 Get the hole drift velocity.
 
virtual bool GetElectronLifetime (const double, const double, const double, double &etau)
 
virtual bool GetHoleLifetime (const double, const double, const double, double &htau)
 

Protected Member Functions

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

- Protected Attributes inherited from Garfield::ComponentFieldMap
bool m_is3d = true
 
int nElements = -1
 
std::vector< Element > elements
 
int nNodes = -1
 
std::vector< Node > nodes
 
unsigned int m_nMaterials = 0
 
std::vector< Material > materials
 
int nWeightingFields = 0
 
std::vector< std::string > wfields
 
std::vector< bool > wfieldsOk
 
bool 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
 
double m_mapvmax
 
std::array< bool, 3 > m_setang
 
bool m_deleteBackground = true
 
bool m_warning = false
 
unsigned int m_nWarnings = 0
 
- Protected Attributes inherited from Garfield::ComponentBase
std::string m_className = "ComponentBase"
 Class name.
 
GeometryBase * m_geometry = nullptr
 Pointer to the geometry.
 
bool m_ready = false
 Ready for use?
 
bool m_activeTraps = false
 Does the component have traps?
 
bool m_hasVelocityMap = false
 Does the component have velocity maps?
 
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.
 
double m_bx0 = 0.
 
double m_by0 = 0.
 
double m_bz0 = 0.
 
bool m_debug = false
 Switch on/off debugging messages.
 

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 14 of file ComponentComsol.cc.

14 : ComponentFieldMap() {
15 m_className = "ComponentComsol";
16}
Garfield::ComponentBase::m_className
std::string m_className
Class name.
Definition: ComponentBase.hh:285

◆ ComponentComsol() [2/2]

Garfield::ComponentComsol::ComponentComsol ( std::string  mesh,
std::string  mplist,
std::string  field 
)

Definition at line 18 of file ComponentComsol.cc.

20 : ComponentFieldMap() {
21 m_className = "ComponentComsol";
22 Initialise(mesh, mplist, field);
23}
Garfield::ComponentComsol::Initialise
bool Initialise(std::string header="mesh.mphtxt", std::string mplist="dielectrics.dat", std::string field="field.txt")
Definition: ComponentComsol.cc:36

◆ ~ComponentComsol()

Garfield::ComponentComsol::~ComponentComsol ( )
inline

Destructor.

Definition at line 15 of file ComponentComsol.hh.

15{}

Member Function Documentation

◆ 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::ComponentBase.

Definition at line 340 of file ComponentComsol.cc.

343 {
344 // Copy the coordinates
345 double x = xin, y = yin, z = zin;
346
347 // Map the coordinates onto field map coordinates
348 bool xmirr, ymirr, zmirr;
349 double rcoordinate, rotation;
350 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
351
352 // Initial values
353 ex = ey = ez = volt = 0.;
354 status = 0;
355 m = NULL;
356
357 // Do not proceed if not properly initialised.
358 if (!m_ready) {
359 status = -10;
360 PrintNotReady("ElectricField");
361 return;
362 }
363
364 if (m_warning) PrintWarning("ElectricField");
365
366 // Find the element that contains this point
367 double t1, t2, t3, t4, jac[4][4], det;
368 const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
369 if (imap < 0) {
370 if (m_debug) {
371 std::cout << m_className << "::ElectricField:\n";
372 std::cout << " Point (" << x << ", " << y << ", " << z
373 << " not in the mesh.\n";
374 }
375 status = -6;
376 return;
377 }
378
379 const Element& element = elements[imap];
380 if (m_debug) {
381 PrintElement("ElectricField", x, y, z, t1, t2, t3, t4, element, 10);
382 }
383 const Node& n0 = nodes[element.emap[0]];
384 const Node& n1 = nodes[element.emap[1]];
385 const Node& n2 = nodes[element.emap[2]];
386 const Node& n3 = nodes[element.emap[3]];
387 const Node& n4 = nodes[element.emap[4]];
388 const Node& n5 = nodes[element.emap[5]];
389 const Node& n6 = nodes[element.emap[6]];
390 const Node& n7 = nodes[element.emap[7]];
391 const Node& n8 = nodes[element.emap[8]];
392 const Node& n9 = nodes[element.emap[9]];
393 // Tetrahedral field
394 volt = n0.v * t1 * (2 * t1 - 1) + n1.v * t2 * (2 * t2 - 1) +
395 n2.v * t3 * (2 * t3 - 1) + n3.v * t4 * (2 * t4 - 1) +
396 4 * n4.v * t1 * t2 + 4 * n5.v * t1 * t3 + 4 * n6.v * t1 * t4 +
397 4 * n7.v * t2 * t3 + 4 * n8.v * t2 * t4 + 4 * n9.v * t3 * t4;
398 ex = -(n0.v * (4 * t1 - 1) * jac[0][1] + n1.v * (4 * t2 - 1) * jac[1][1] +
399 n2.v * (4 * t3 - 1) * jac[2][1] + n3.v * (4 * t4 - 1) * jac[3][1] +
400 n4.v * (4 * t2 * jac[0][1] + 4 * t1 * jac[1][1]) +
401 n5.v * (4 * t3 * jac[0][1] + 4 * t1 * jac[2][1]) +
402 n6.v * (4 * t4 * jac[0][1] + 4 * t1 * jac[3][1]) +
403 n7.v * (4 * t3 * jac[1][1] + 4 * t2 * jac[2][1]) +
404 n8.v * (4 * t4 * jac[1][1] + 4 * t2 * jac[3][1]) +
405 n9.v * (4 * t4 * jac[2][1] + 4 * t3 * jac[3][1])) /
406 det;
407 ey = -(n0.v * (4 * t1 - 1) * jac[0][2] + n1.v * (4 * t2 - 1) * jac[1][2] +
408 n2.v * (4 * t3 - 1) * jac[2][2] + n3.v * (4 * t4 - 1) * jac[3][2] +
409 n4.v * (4 * t2 * jac[0][2] + 4 * t1 * jac[1][2]) +
410 n5.v * (4 * t3 * jac[0][2] + 4 * t1 * jac[2][2]) +
411 n6.v * (4 * t4 * jac[0][2] + 4 * t1 * jac[3][2]) +
412 n7.v * (4 * t3 * jac[1][2] + 4 * t2 * jac[2][2]) +
413 n8.v * (4 * t4 * jac[1][2] + 4 * t2 * jac[3][2]) +
414 n9.v * (4 * t4 * jac[2][2] + 4 * t3 * jac[3][2])) /
415 det;
416 ez = -(n0.v * (4 * t1 - 1) * jac[0][3] + n1.v * (4 * t2 - 1) * jac[1][3] +
417 n2.v * (4 * t3 - 1) * jac[2][3] + n3.v * (4 * t4 - 1) * jac[3][3] +
418 n4.v * (4 * t2 * jac[0][3] + 4 * t1 * jac[1][3]) +
419 n5.v * (4 * t3 * jac[0][3] + 4 * t1 * jac[2][3]) +
420 n6.v * (4 * t4 * jac[0][3] + 4 * t1 * jac[3][3]) +
421 n7.v * (4 * t3 * jac[1][3] + 4 * t2 * jac[2][3]) +
422 n8.v * (4 * t4 * jac[1][3] + 4 * t2 * jac[3][3]) +
423 n9.v * (4 * t4 * jac[2][3] + 4 * t3 * jac[3][3])) /
424 det;
425
426 // Transform field to global coordinates
427 UnmapFields(ex, ey, ez, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
428 // std::cout << "ef @(" << xin << ", " << yin << ", " << zin << ") = " <<
429 // volt << "\n";
430
431 // Drift medium?
432 if (m_debug) {
433 std::cout << m_className << "::ElectricField:\n";
434 std::cout << " Material " << element.matmap << ", drift flag "
435 << materials[element.matmap].driftmedium << "\n";
436 }
437 m = materials[element.matmap].medium;
438 status = -5;
439 if (materials[element.matmap].driftmedium) {
440 if (m && m->IsDriftable()) status = 0;
441 }
442}
Garfield::ComponentBase::m_ready
bool m_ready
Ready for use?
Definition: ComponentBase.hh:291
Garfield::ComponentBase::m_debug
bool m_debug
Switch on/off debugging messages.
Definition: ComponentBase.hh:310
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:2074
Garfield::ComponentFieldMap::PrintWarning
void PrintWarning(const std::string &header)
Definition: ComponentFieldMap.hh:193
Garfield::ComponentFieldMap::materials
std::vector< Material > materials
Definition: ComponentFieldMap.hh:125
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:298
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:2396
Garfield::ComponentFieldMap::elements
std::vector< Element > elements
Definition: ComponentFieldMap.hh:100
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:2198
Garfield::ComponentFieldMap::PrintNotReady
void PrintNotReady(const std::string &header) const
Definition: ComponentFieldMap.hh:199

◆ 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::ComponentBase.

Definition at line 333 of file ComponentComsol.cc.

335 {
336 double v = 0.;
337 ElectricField(x, y, z, ex, ey, ez, v, m, status);
338}
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:333

Referenced by ElectricField().

◆ GetAspectRatio()

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

Implements Garfield::ComponentFieldMap.

Definition at line 666 of file ComponentComsol.cc.

667 {
668 if (i >= elements.size()) {
669 dmin = dmax = 0.;
670 return;
671 }
672
673 const Element& element = elements[i];
674 const int np = 4;
675 // Loop over all pairs of vertices.
676 for (int j = 0; j < np - 1; ++j) {
677 const Node& nj = nodes[element.emap[j]];
678 for (int k = j + 1; k < np; ++k) {
679 const Node& nk = nodes[element.emap[k]];
680 // Compute distance.
681 const double dx = nj.x - nk.x;
682 const double dy = nj.y - nk.y;
683 const double dz = nj.z - nk.z;
684 const double dist = sqrt(dx * dx + dy * dy + dz * dz);
685 if (k == 1) {
686 dmin = dmax = dist;
687 } else {
688 if (dist < dmin) dmin = dist;
689 if (dist > dmax) dmax = dist;
690 }
691 }
692 }
693}
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 645 of file ComponentComsol.cc.

645 {
646 if (i >= elements.size()) return 0.;
647 const Element& element = elements[i];
648 const Node& n0 = nodes[element.emap[0]];
649 const Node& n1 = nodes[element.emap[1]];
650 const Node& n2 = nodes[element.emap[2]];
651 const Node& n3 = nodes[element.emap[3]];
652
653 // Uses formula V = |a (dot) b x c|/6
654 // with a => "3", b => "1", c => "2" and origin "0"
655 const double vol =
656 fabs((n3.x - n0.x) *
657 ((n1.y - n0.y) * (n2.z - n0.z) - (n2.y - n0.y) * (n1.z - n0.z)) +
658 (n3.y - n0.y) *
659 ((n1.z - n0.z) * (n2.x - n0.x) - (n2.z - n0.z) * (n1.x - n0.x)) +
660 (n3.z - n0.z) * ((n1.x - n0.x) * (n2.y - n0.y) -
661 (n3.x - n0.x) * (n1.y - n0.y))) /
662 6.;
663 return vol;
664}
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::ComponentBase.

Definition at line 600 of file ComponentComsol.cc.

601 {
602 // Copy the coordinates
603 double x = xin, y = yin, z = zin;
604
605 // Map the coordinates onto field map coordinates
606 bool xmirr, ymirr, zmirr;
607 double rcoordinate, rotation;
608 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
609
610 // Do not proceed if not properly initialised.
611 if (!m_ready) {
612 PrintNotReady("GetMedium");
613 return nullptr;
614 }
615 if (m_warning) PrintWarning("GetMedium");
616
617 // Find the element that contains this point
618 double t1, t2, t3, t4, jac[4][4], det;
619 const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
620 if (imap < 0) {
621 if (m_debug) {
622 std::cout << m_className << "::GetMedium:\n";
623 std::cout << " Point (" << x << ", " << y << ", " << z
624 << ") not in the mesh.\n";
625 }
626 return nullptr;
627 }
628 const Element& element = elements[imap];
629 if (element.matmap >= m_nMaterials) {
630 if (m_debug) {
631 std::cerr << m_className << "::GetMedium:\n";
632 std::cerr << " Point (" << x << ", " << y
633 << ") has out of range material number " << imap << ".\n";
634 }
635 return nullptr;
636 }
637
638 if (m_debug) {
639 PrintElement("GetMedium", x, y, z, t1, t2, t3, t4, element, 10);
640 }
641
642 return materials[element.matmap].medium;
643}

◆ Initialise()

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

Definition at line 36 of file ComponentComsol.cc.

37 {
38 m_ready = false;
39 m_warning = false;
40 m_nWarnings = 0;
41
42 double unit = 100.0; // m
43
44 std::string line;
45
46 // Open the materials file.
47 materials.clear();
48 std::ifstream fmplist;
49 fmplist.open(mplist.c_str(), std::ios::in);
50 if (fmplist.fail()) {
51 std::cerr << m_className << "::Initialise:\n";
52 std::cerr << " Could not open result file " << mplist
53 << " for reading.\n";
54 return false;
55 }
56 fmplist >> m_nMaterials;
57 for (unsigned int i = 0; i < m_nMaterials; ++i) {
58 Material newMaterial;
59 newMaterial.driftmedium = true;
60 newMaterial.medium = nullptr;
61 newMaterial.ohm = -1;
62 fmplist >> newMaterial.eps;
63 materials.push_back(newMaterial);
64 }
65 {
66 // add default material
67 Material newMaterial;
68 newMaterial.driftmedium = false;
69 newMaterial.medium = nullptr;
70 newMaterial.eps = newMaterial.ohm = -1;
71 materials.push_back(newMaterial);
72 m_nMaterials++;
73 }
74 std::map<int, int> domain2material;
75 int d2msize;
76 fmplist >> d2msize;
77 for (int i = 0; i < d2msize; ++i) {
78 int domain;
79 fmplist >> domain;
80 fmplist >> domain2material[domain];
81 }
82 fmplist.close();
83
84 nodes.clear();
85 std::ifstream fmesh;
86 fmesh.open(mesh.c_str(), std::ios::in);
87 if (fmesh.fail()) {
88 std::cerr << m_className << "::Initialise:\n";
89 std::cerr << " Could not open nodes file " << mesh << " for reading.\n";
90 return false;
91 }
92
93 do {
94 std::getline(fmesh, line);
95 } while (!ends_with(line, "# number of mesh points"));
96 nNodes = readInt(line);
97
98 std::cout << m_className << "::Initialise:\n";
99 std::cout << " Read " << nNodes << " nodes from file " << mesh << ".\n";
100 do {
101 std::getline(fmesh, line);
102 } while (line != "# Mesh point coordinates");
103 double minx = 1e100, miny = 1e100, minz = 1e100, maxx = -1e100, maxy = -1e100,
104 maxz = -1e100;
105 for (int i = 0; i < nNodes; ++i) {
106 Node newNode;
107 fmesh >> newNode.x >> newNode.y >> newNode.z;
108 newNode.x *= unit;
109 newNode.y *= unit;
110 newNode.z *= unit;
111 nodes.push_back(newNode);
112 minx = std::min(minx, newNode.x);
113 maxx = std::max(maxx, newNode.x);
114 miny = std::min(miny, newNode.y);
115 maxy = std::max(maxy, newNode.y);
116 minz = std::min(minz, newNode.z);
117 maxz = std::max(maxz, newNode.z);
118 }
119 std::cout << minx << " < x < " << maxx << "\n";
120 std::cout << miny << " < y < " << maxy << "\n";
121 std::cout << minz << " < z < " << maxz << "\n";
122
123 do {
124 std::getline(fmesh, line);
125 } while (line != "4 tet2 # type name");
126 do {
127 std::getline(fmesh, line);
128 } while (!ends_with(line, "# number of elements"));
129 nElements = readInt(line);
130 elements.clear();
131 std::cout << m_className << "::Initialise:\n";
132 std::cout << " Read " << nElements << " elements from file " << mesh
133 << ".\n";
134 std::getline(fmesh, line);
135 // elements 6 & 7 are swapped due to differences in COMSOL and ANSYS
136 // representation
137 int perm[10] = {0, 1, 2, 3, 4, 5, 7, 6, 8, 9};
138 for (int i = 0; i < nElements; ++i) {
139 Element newElement;
140 newElement.degenerate = false;
141 for (int j = 0; j < 10; ++j) {
142 fmesh >> newElement.emap[perm[j]];
143 }
144 elements.push_back(newElement);
145 }
146
147 do {
148 std::getline(fmesh, line);
149 } while (line != "# Geometric entity indices");
150 for (int i = 0; i < nElements; ++i) {
151 int domain;
152 fmesh >> domain;
153 elements[i].matmap = domain2material.count(domain) ? domain2material[domain]
154 : m_nMaterials - 1;
155 }
156 fmesh.close();
157
158 std::map<Node, std::vector<int>, nodeCmp> nodeIdx;
159 for (int i = 0; i < nNodes; ++i) {
160 nodeIdx[nodes[i]].push_back(i);
161 }
162 std::cout << "Map size: " << nodeIdx.size() << std::endl;
163
164 std::ifstream ffield;
165 ffield.open(field.c_str(), std::ios::in);
166 if (ffield.fail()) {
167 std::cerr << m_className << "::Initialise:\n";
168 std::cerr << " Could not open field potentials file " << field
169 << " for reading.\n";
170 return false;
171 }
172 do {
173 std::getline(ffield, line);
174 } while (line.substr(0, 81) !=
175 "% x y z "
176 " V (V)");
177 {
178 std::istringstream sline(line);
179 std::string token;
180 sline >> token; // %
181 sline >> token; // x
182 sline >> token; // y
183 sline >> token; // z
184 sline >> token; // V
185 sline >> token; // (V)
186 while (sline >> token) {
187 std::cout << m_className << "::Initialise:\n";
188 std::cout << " Reading data for weighting field " << token << ".\n";
189 nWeightingFields++;
190 wfields.push_back(token);
191 wfieldsOk.push_back(true);
192 sline >> token; // (V)
193 }
194 }
195 for (int i = 0; i < nNodes; ++i) {
196 Node tmp;
197 ffield >> tmp.x >> tmp.y >> tmp.z >> tmp.v;
198 tmp.x *= unit;
199 tmp.y *= unit;
200 tmp.z *= unit;
201 for (int j = 0; j < nWeightingFields; ++j) {
202 double w;
203 ffield >> w;
204 tmp.w.push_back(w);
205 }
206 int closest = -1;
207 double closestDist = 1;
208 const unsigned int nIdx = nodeIdx[tmp].size();
209 // for (int j : nodeIdx[tmp]) {
210 for (unsigned int k = 0; k < nIdx; ++k) {
211 int j = nodeIdx[tmp][k];
212 double dist = (tmp.x - nodes[j].x) * (tmp.x - nodes[j].x) +
213 (tmp.y - nodes[j].y) * (tmp.y - nodes[j].y) +
214 (tmp.z - nodes[j].z) * (tmp.z - nodes[j].z);
215 if (dist < closestDist) {
216 closestDist = dist;
217 closest = j;
218 }
219 }
220 if (closest == -1) {
221 std::cerr << m_className << "::Initialise:\n";
222 std::cerr << " Could not match the node from field potentials file: "
223 << tmp.x << " " << tmp.y << " " << tmp.z << "\n.";
224 return false;
225 }
226 nodes[closest].v = tmp.v;
227 nodes[closest].w = tmp.w;
228 }
229
230 m_ready = true;
231
232 // for (int i = 0; i < nNodes; ++i) {
233 // double ex, ey, ez, v;
234 // Medium* m;
235 // int status;
236 // ElectricField(nodes[i].x, nodes[i].y, nodes[i].z, ex, ey, ez, v, m,
237 // status);
238 // std::cout << "Field at " << nodes[i].x << " " << nodes[i].y << " " <<
239 // nodes[i].z << ": " << ex << " " << ey << " " << ez << " " << v << "\n";
240 // }
241
242 // Establish the ranges.
243 SetRange();
244 UpdatePeriodicity();
245 return true;
246}
Garfield::ComponentComsol::UpdatePeriodicity
void UpdatePeriodicity() override
Verify periodicities.
Definition: ComponentComsol.hh:39
Garfield::ComponentFieldMap::SetRange
virtual void SetRange()
Calculate x, y, z, V and angular ranges.
Definition: ComponentFieldMap.cc:1934
Garfield::ComponentFieldMap::wfields
std::vector< std::string > wfields
Definition: ComponentFieldMap.hh:128
Garfield::readInt
int readInt(std::string s)
Definition: ComponentComsol.cc:29
Garfield::ends_with
bool ends_with(std::string s, std::string t)
Definition: ComponentComsol.cc:25

Referenced by ComponentComsol().

◆ SetWeightingField()

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

Definition at line 248 of file ComponentComsol.cc.

248 {
249 double unit = 100.0; // m;
250
251 if (!m_ready) {
252 std::cerr << m_className << "::SetWeightingField:\n";
253 std::cerr << " No valid field map is present.\n";
254 std::cerr << " Weighting field cannot be added.\n";
255 return false;
256 }
257
258 // Open the voltage list.
259 std::ifstream ffield;
260 ffield.open(field.c_str(), std::ios::in);
261 if (ffield.fail()) {
262 std::cerr << m_className << "::Initialise:\n";
263 std::cerr << " Could not open field potentials file " << field
264 << " for reading.\n";
265 return false;
266 }
267
268 // Check if a weighting field with the same label alm_ready exists.
269 int iw = nWeightingFields;
270 for (int i = nWeightingFields; i--;) {
271 if (wfields[i] == label) {
272 iw = i;
273 break;
274 }
275 }
276 if (iw == nWeightingFields) {
277 ++nWeightingFields;
278 wfields.resize(nWeightingFields);
279 wfieldsOk.resize(nWeightingFields);
280 for (int j = 0; j < nNodes; ++j) {
281 nodes[j].w.resize(nWeightingFields);
282 }
283 } else {
284 std::cout << m_className << "::SetWeightingField:\n";
285 std::cout << " Replacing existing weighting field " << label << ".\n";
286 }
287 wfields[iw] = label;
288 wfieldsOk[iw] = false;
289 std::map<Node, std::vector<int>, nodeCmp> nodeIdx;
290 for (int i = 0; i < nNodes; ++i) {
291 nodeIdx[nodes[i]].push_back(i);
292 }
293 std::cout << "Map size: " << nodeIdx.size() << std::endl;
294
295 std::string line;
296 do {
297 std::getline(ffield, line);
298 } while (line !=
299 "% x y z "
300 " V (V)");
301 for (int i = 0; i < nNodes; ++i) {
302 Node tmp;
303 ffield >> tmp.x >> tmp.y >> tmp.z >> tmp.v;
304 tmp.x *= unit;
305 tmp.y *= unit;
306 tmp.z *= unit;
307 int closest = -1;
308 double closestDist = 1;
309 const unsigned int nIdx = nodeIdx[tmp].size();
310 // for (int j : nodeIdx[tmp]) {
311 for (unsigned int k = 0; k < nIdx; ++k) {
312 int j = nodeIdx[tmp][k];
313 double dist = (tmp.x - nodes[j].x) * (tmp.x - nodes[j].x) +
314 (tmp.y - nodes[j].y) * (tmp.y - nodes[j].y) +
315 (tmp.z - nodes[j].z) * (tmp.z - nodes[j].z);
316 if (dist < closestDist) {
317 closestDist = dist;
318 closest = j;
319 }
320 }
321 if (closest == -1) {
322 std::cerr << m_className << "::Initialise:\n";
323 std::cerr << " Could not match the node from field potentials file: "
324 << tmp.x << " " << tmp.y << " " << tmp.z << "\n.";
325 return false;
326 }
327 nodes[closest].w[iw] = tmp.v;
328 }
329
330 return true;
331}

◆ UpdatePeriodicity()

void Garfield::ComponentComsol::UpdatePeriodicity ( )
inlineoverrideprotectedvirtual

Verify periodicities.

Implements Garfield::ComponentBase.

Definition at line 39 of file ComponentComsol.hh.

39{ UpdatePeriodicityCommon(); }

Referenced by Initialise().

◆ 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::ComponentBase.

Definition at line 444 of file ComponentComsol.cc.

446 {
447 // Initial values
448 wx = wy = wz = 0;
449
450 // Do not proceed if not properly initialised.
451 if (!m_ready) return;
452
453 // Look for the label.
454 int iw = 0;
455 bool found = false;
456 for (int i = nWeightingFields; i--;) {
457 if (wfields[i] == label) {
458 iw = i;
459 found = true;
460 break;
461 }
462 }
463
464 // Do not proceed if the requested weighting field does not exist.
465 if (!found) return;
466 // Check if the weighting field is properly initialised.
467 if (!wfieldsOk[iw]) return;
468
469 // Copy the coordinates.
470 double x = xin, y = yin, z = zin;
471
472 // Map the coordinates onto field map coordinates
473 bool xmirr, ymirr, zmirr;
474 double rcoordinate, rotation;
475 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
476
477 if (m_warning) PrintWarning("WeightingField");
478
479 // Find the element that contains this point.
480 double t1, t2, t3, t4, jac[4][4], det;
481 const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
482 // Check if the point is in the mesh.
483 if (imap < 0) return;
484
485 const Element& element = elements[imap];
486 if (m_debug) {
487 PrintElement("WeightingField", x, y, z, t1, t2, t3, t4, element, 10, iw);
488 }
489 const Node& n0 = nodes[element.emap[0]];
490 const Node& n1 = nodes[element.emap[1]];
491 const Node& n2 = nodes[element.emap[2]];
492 const Node& n3 = nodes[element.emap[3]];
493 const Node& n4 = nodes[element.emap[4]];
494 const Node& n5 = nodes[element.emap[5]];
495 const Node& n6 = nodes[element.emap[6]];
496 const Node& n7 = nodes[element.emap[7]];
497 const Node& n8 = nodes[element.emap[8]];
498 const Node& n9 = nodes[element.emap[9]];
499 // Tetrahedral field
500 wx = -(n0.w[iw] * (4 * t1 - 1) * jac[0][1] +
501 n1.w[iw] * (4 * t2 - 1) * jac[1][1] +
502 n2.w[iw] * (4 * t3 - 1) * jac[2][1] +
503 n3.w[iw] * (4 * t4 - 1) * jac[3][1] +
504 n4.w[iw] * (4 * t2 * jac[0][1] + 4 * t1 * jac[1][1]) +
505 n5.w[iw] * (4 * t3 * jac[0][1] + 4 * t1 * jac[2][1]) +
506 n6.w[iw] * (4 * t4 * jac[0][1] + 4 * t1 * jac[3][1]) +
507 n7.w[iw] * (4 * t3 * jac[1][1] + 4 * t2 * jac[2][1]) +
508 n8.w[iw] * (4 * t4 * jac[1][1] + 4 * t2 * jac[3][1]) +
509 n9.w[iw] * (4 * t4 * jac[2][1] + 4 * t3 * jac[3][1])) /
510 det;
511
512 wy = -(n0.w[iw] * (4 * t1 - 1) * jac[0][2] +
513 n1.w[iw] * (4 * t2 - 1) * jac[1][2] +
514 n2.w[iw] * (4 * t3 - 1) * jac[2][2] +
515 n3.w[iw] * (4 * t4 - 1) * jac[3][2] +
516 n4.w[iw] * (4 * t2 * jac[0][2] + 4 * t1 * jac[1][2]) +
517 n5.w[iw] * (4 * t3 * jac[0][2] + 4 * t1 * jac[2][2]) +
518 n6.w[iw] * (4 * t4 * jac[0][2] + 4 * t1 * jac[3][2]) +
519 n7.w[iw] * (4 * t3 * jac[1][2] + 4 * t2 * jac[2][2]) +
520 n8.w[iw] * (4 * t4 * jac[1][2] + 4 * t2 * jac[3][2]) +
521 n9.w[iw] * (4 * t4 * jac[2][2] + 4 * t3 * jac[3][2])) /
522 det;
523
524 wz = -(n0.w[iw] * (4 * t1 - 1) * jac[0][3] +
525 n1.w[iw] * (4 * t2 - 1) * jac[1][3] +
526 n2.w[iw] * (4 * t3 - 1) * jac[2][3] +
527 n3.w[iw] * (4 * t4 - 1) * jac[3][3] +
528 n4.w[iw] * (4 * t2 * jac[0][3] + 4 * t1 * jac[1][3]) +
529 n5.w[iw] * (4 * t3 * jac[0][3] + 4 * t1 * jac[2][3]) +
530 n6.w[iw] * (4 * t4 * jac[0][3] + 4 * t1 * jac[3][3]) +
531 n7.w[iw] * (4 * t3 * jac[1][3] + 4 * t2 * jac[2][3]) +
532 n8.w[iw] * (4 * t4 * jac[1][3] + 4 * t2 * jac[3][3]) +
533 n9.w[iw] * (4 * t4 * jac[2][3] + 4 * t3 * jac[3][3])) /
534 det;
535
536 // Transform field to global coordinates
537 UnmapFields(wx, wy, wz, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
538}

◆ 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::ComponentBase.

Definition at line 540 of file ComponentComsol.cc.

542 {
543 // Do not proceed if not properly initialised.
544 if (!m_ready) return 0.;
545
546 // Look for the label.
547 int iw = 0;
548 bool found = false;
549 for (int i = nWeightingFields; i--;) {
550 if (wfields[i] == label) {
551 iw = i;
552 found = true;
553 break;
554 }
555 }
556
557 // Do not proceed if the requested weighting field does not exist.
558 if (!found) return 0.;
559 // Check if the weighting field is properly initialised.
560 if (!wfieldsOk[iw]) return 0.;
561
562 // Copy the coordinates.
563 double x = xin, y = yin, z = zin;
564
565 // Map the coordinates onto field map coordinates.
566 bool xmirr, ymirr, zmirr;
567 double rcoordinate, rotation;
568 MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
569
570 if (m_warning) PrintWarning("WeightingPotential");
571
572 // Find the element that contains this point.
573 double t1, t2, t3, t4, jac[4][4], det;
574 const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
575 if (imap < 0) return 0.;
576
577 const Element& element = elements[imap];
578 if (m_debug) {
579 PrintElement("WeightingPotential", x, y, z, t1, t2, t3, t4, element, 10,
580 iw);
581 }
582 const Node& n0 = nodes[element.emap[0]];
583 const Node& n1 = nodes[element.emap[1]];
584 const Node& n2 = nodes[element.emap[2]];
585 const Node& n3 = nodes[element.emap[3]];
586 const Node& n4 = nodes[element.emap[4]];
587 const Node& n5 = nodes[element.emap[5]];
588 const Node& n6 = nodes[element.emap[6]];
589 const Node& n7 = nodes[element.emap[7]];
590 const Node& n8 = nodes[element.emap[8]];
591 const Node& n9 = nodes[element.emap[9]];
592 // Tetrahedral field
593 return n0.w[iw] * t1 * (2 * t1 - 1) + n1.w[iw] * t2 * (2 * t2 - 1) +
594 n2.w[iw] * t3 * (2 * t3 - 1) + n3.w[iw] * t4 * (2 * t4 - 1) +
595 4 * n4.w[iw] * t1 * t2 + 4 * n5.w[iw] * t1 * t3 +
596 4 * n6.w[iw] * t1 * t4 + 4 * n7.w[iw] * t2 * t3 +
597 4 * n8.w[iw] * t2 * t4 + 4 * n9.w[iw] * t3 * t4;
598}
The documentation for this class was generated from the following files: