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

Additional Inherited Members

- Protected Types inherited from Garfield::ComponentFieldMap
enum class  ElementType { Unknown = 0 , Serendipity = 5 , CurvedTetrahedron = 13 }
 
- Protected Member Functions inherited from Garfield::ComponentFieldMap
void Reset () override
 Reset the component.
 
void Prepare ()
 
virtual void SetRange ()
 
void UpdatePeriodicity () override
 Verify periodicities.
 
void UpdatePeriodicity2d ()
 
void UpdatePeriodicityCommon ()
 
bool SetDefaultDriftMedium ()
 Find lowest epsilon, check for eps = 0, set default drift media flags.
 
int Field (const double x, const double y, const double z, double &fx, double &fy, double &fz, int &iel, const std::vector< double > &potentials) const
 Compute the electric/weighting field.
 
double Potential (const double x, const double y, const double z, const std::vector< double > &potentials) const
 Compute the electrostatic/weighting potential.
 
int FindElement5 (const double x, const double y, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det) const
 Find the element for a point in curved quadratic quadrilaterals.
 
int FindElement13 (const double x, const double y, const double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det) const
 Find the element for a point in curved quadratic tetrahedra.
 
int FindElementCube (const double x, const double y, const double z, double &t1, double &t2, double &t3, TMatrixD *&jac, std::vector< TMatrixD * > &dN) const
 Find the element for a point in a cube.
 
void MapCoordinates (double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
 Move (xpos, ypos, zpos) to field map coordinates.
 
void UnmapFields (double &ex, double &ey, double &ez, const double xpos, const double ypos, const double zpos, const bool xmirrored, const bool ymirrored, const bool zmirrored, const double rcoordinate, const double rotation) const
 Move (ex, ey, ez) to global coordinates.
 
virtual double GetElementVolume (const size_t i) const
 
virtual void GetAspectRatio (const size_t i, double &dmin, double &dmax) const
 
void PrintWarning (const std::string &header)
 
void PrintNotReady (const std::string &header) const
 
void PrintCouldNotOpen (const std::string &header, const std::string &filename) const
 
void PrintElement (const std::string &header, const double x, const double y, const double z, const double t1, const double t2, const double t3, const double t4, const size_t i, const std::vector< double > &potential) const
 
void TimeInterpolation (const double t, double &f0, double &f1, int &i0, int &i1)
 Interpolation of potential between two time slices.
 
- Static Protected Member Functions inherited from Garfield::ComponentFieldMap
static double ScalingFactor (std::string unit)
 
static double Potential3 (const std::array< double, 6 > &v, const std::array< double, 3 > &t)
 Interpolate the potential in a triangle.
 
static void Field3 (const std::array< double, 6 > &v, const std::array< double, 3 > &t, double jac[4][4], const double det, double &ex, double &ey)
 Interpolate the field in a triangle.
 
static double Potential5 (const std::array< double, 8 > &v, const std::array< double, 2 > &t)
 Interpolate the potential in a curved quadrilateral.
 
static void Field5 (const std::array< double, 8 > &v, const std::array< double, 2 > &t, double jac[4][4], const double det, double &ex, double &ey)
 Interpolate the field in a curved quadrilateral.
 
static double Potential13 (const std::array< double, 10 > &v, const std::array< double, 4 > &t)
 Interpolate the potential in a curved quadratic tetrahedron.
 
static void Field13 (const std::array< double, 10 > &v, const std::array< double, 4 > &t, double jac[4][4], const double det, double &ex, double &ey, double &ez)
 Interpolate the field in a curved quadratic tetrahedron.
 
static int ReadInteger (char *token, int def, bool &error)
 
static double ReadDouble (char *token, double def, bool &error)
 
- Protected Attributes inherited from Garfield::ComponentFieldMap
bool m_is3d = true
 
ElementType m_elementType = ElementType::CurvedTetrahedron
 
std::vector< Element > m_elements
 
std::vector< int > m_elementIndices
 
std::vector< bool > m_degenerate
 
std::vector< std::array< double, 3 > > m_bbMin
 
std::vector< std::array< double, 3 > > m_bbMax
 
std::vector< std::array< std::array< double, 3 >, 4 > > m_w12
 
std::vector< Node > m_nodes
 
std::vector< double > m_pot
 
std::map< std::string, std::vector< double > > m_wpot
 
std::map< std::string, std::vector< std::vector< double > > > m_dwpot
 
std::vector< Material > m_materials
 
std::map< std::string, WeightingFieldCopy > m_wfieldCopies
 
std::vector< double > m_wdtimes
 
bool m_hasBoundingBox = false
 
std::array< double, 3 > m_minBoundingBox = {{0., 0., 0.}}
 
std::array< double, 3 > m_maxBoundingBox = {{0., 0., 0.}}
 
std::array< double, 3 > m_mapmin = {{0., 0., 0.}}
 
std::array< double, 3 > m_mapmax = {{0., 0., 0.}}
 
std::array< double, 3 > m_mapamin = {{0., 0., 0.}}
 
std::array< double, 3 > m_mapamax = {{0., 0., 0.}}
 
std::array< double, 3 > m_mapna = {{0., 0., 0.}}
 
std::array< double, 3 > m_cells = {{0., 0., 0.}}
 
double m_mapvmin = 0.
 
double m_mapvmax = 0.
 
std::array< bool, 3 > m_setang
 
bool m_deleteBackground = true
 
bool m_warning = false
 
unsigned int m_nWarnings = 0
 
bool m_printConvergenceWarnings = true
 
- Protected Attributes inherited from Garfield::Component
std::string m_className = "Component"
 Class name.
 
Geometry * m_geometry = nullptr
 Pointer to the geometry.
 
std::array< double, 3 > m_b0 = {{0., 0., 0.}}
 Constant magnetic field.
 
bool m_ready = false
 Ready for use?
 
bool m_debug = false
 Switch on/off debugging messages.
 
std::array< bool, 3 > m_periodic = {{false, false, false}}
 Simple periodicity in x, y, z.
 
std::array< bool, 3 > m_mirrorPeriodic = {{false, false, false}}
 Mirror periodicity in x, y, z.
 
std::array< bool, 3 > m_axiallyPeriodic = {{false, false, false}}
 Axial periodicity in x, y, z.
 
std::array< bool, 3 > m_rotationSymmetric = {{false, false, false}}
 Rotation symmetry around x-axis, y-axis, z-axis.
 

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

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

Referenced by ComponentComsol().

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

52 : ComponentComsol() {
53 Initialise(mesh, mplist, field, unit);
54}
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:56
Garfield::ComponentComsol::ComponentComsol
ComponentComsol()
Default constructor.
Definition ComponentComsol.cc:46

◆ ~ComponentComsol()

Garfield::ComponentComsol::~ComponentComsol ( )
inline

Destructor.

Definition at line 17 of file ComponentComsol.hh.

17{}

Member Function Documentation

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

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

Referenced by ComponentComsol().

◆ SetDynamicWeightingPotential()

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

Import the time-dependent weighting field maps.

Definition at line 432 of file ComponentComsol.cc.

433 {
434 if (!m_ready) {
435 std::cerr << m_className << "::SetDynamicWeightingPotential:\n"
436 << " No valid field map is present.\n"
437 << " Weighting fields cannot be added.\n";
438 return false;
439 }
440
441 if (!m_timeset && !GetTimeInterval(field)) return false;
442
443 if (!m_timeset) {
444 std::cerr << m_className << "::SetDynamicWeightingPotential:\n"
445 << " No valid times slices of potential set.\n"
446 << " Please add the time slices.\n";
447 return false;
448 }
449
450 const int T = m_wdtimes.size();
451
452 // Open the voltage list.
453 std::ifstream ffield(field);
454 if (!ffield) {
455 PrintCouldNotOpen("SetDynamicWeightingPotential", field);
456 return false;
457 }
458
459 // Check if a weighting field with the same label already exists.
460 if (m_wpot.count(label) > 0) {
461 std::cout << m_className << "::SetDynamicWeightingPotential:\n"
462 << " Replacing existing weighting field " << label << ".\n";
463 m_wpot[label].clear();
464 }
465 if (m_dwpot.count(label) > 0) {
466 m_dwpot[label].clear();
467 }
468
469 std::vector<double> pot(m_nodes.size(), 0.);
470 std::vector<std::vector<double> > dpot(m_nodes.size());
471
472 // Build a k-d tree from the node coordinates.
473 std::vector<std::vector<double> > points;
474 for (const auto &node : m_nodes) {
475 std::vector<double> point = {node.x, node.y, node.z};
476 points.push_back(std::move(point));
477 }
478 KDTree kdtree(points);
479
480 std::string line;
481 int nLines = 1;
482 const int nNodes = m_nodes.size();
483 const unsigned int nPrint =
484 std::pow(10, static_cast<unsigned int>(
485 std::max(std::floor(std::log10(nNodes)) - 1, 1.)));
486 std::cout << m_className << "::SetDynamicWeightingPotential:\n"
487 << " Reading weighting potentials for " << label << ".\n";
488 PrintProgress(0.);
489
490 while (std::getline(ffield, line)) {
491 // Skip empty lines.
492 if (line.empty()) continue;
493 // Skip comments.
494 if (isComment(line)) continue;
495
496 std::istringstream data(line);
497 double x = 0., y = 0., z = 0.;
498 data >> x >> y >> z;
499 x *= m_unit;
500 y *= m_unit;
501 z *= m_unit;
502 if (!CheckInRange(x, y, z)) continue;
503 std::vector<KDTreeResult> res;
504 kdtree.n_nearest({x, y, z}, 1, res);
505 if (res.empty()) {
506 std::cerr << m_className << "::SetDynamicWeightingPotential:\n"
507 << " Could not find a matching mesh node for point ("
508 << x << ", " << y << ", " << z << ")\n.";
509 ffield.close();
510 return false;
511 }
512 if (res[0].dis > MaxNodeDistance) continue;
513
514 double p = 0.;
515 double p0 = 0.;
516 std::vector<double> pvect;
517 for (int i = 0; i < T; i++) {
518 data >> p;
519 if (i == 0) p0 = p;
520 pvect.push_back(p - p0);
521 }
522 const size_t k = res[0].idx;
523 dpot[k] = pvect;
524 pot[k] = p0;
525
526 if ((nLines + 1) % nPrint == 0) {
527 PrintProgress(double(nLines + 1) / nNodes);
528 }
529 nLines++;
530 }
531
532 PrintProgress(1.);
533 std::cout << std::endl
534 << m_className << "::SetDynamicWeightingPotential: Done.\n";
535 ffield.close();
536 m_wpot[label] = std::move(pot);
537 m_dwpot[label] = std::move(dpot);
538 return true;
539}
Garfield::ComponentFieldMap::m_wdtimes
std::vector< double > m_wdtimes
Definition ComponentFieldMap.hh:201
Garfield::ComponentFieldMap::m_dwpot
std::map< std::string, std::vector< std::vector< double > > > m_dwpot
Definition ComponentFieldMap.hh:176

◆ SetImportRange()

void Garfield::ComponentComsol::SetImportRange ( const double xmin,
const double xmax,
const double ymin,
const double ymax,
const double zmin,
const double zmax )
inline

Definition at line 19 of file ComponentComsol.hh.

20 {
21 // TODO: Must happen before initialise function
22 m_range.set = true;
23 m_range.xmin = xmin;
24 m_range.ymin = ymin;
25 m_range.zmin = zmin;
26
27 m_range.xmax = xmax;
28 m_range.ymax = ymax;
29 m_range.zmax = zmax;
30 }

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

542 {
543 std::cout << std::endl
544 << m_className
545 << "::SetTimeInterval: Overwriting time interval of weighting "
546 "potential.\n";
547
548 if (m_wdtimes.empty()) {
549 double t = mint;
550 while (t <= maxt) {
551 m_wdtimes.push_back(t);
552 t += stept;
553 }
554 }
555 m_timeset = true;
556
557 std::cout << std::endl
558 << m_className
559 << "::SetTimeInterval: Time of weighting potential set for t in ["
560 << mint << "," << maxt << "].\n";
561}

◆ SetWeightingField()

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

Definition at line 45 of file ComponentComsol.hh.

45 {
46 return SetWeightingPotential(file, label);
47 }
Garfield::ComponentComsol::SetWeightingPotential
bool SetWeightingPotential(const std::string &file, const std::string &label)
Import the weighting potential maps.
Definition ComponentComsol.cc:341

◆ SetWeightingPotential()

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

Import the weighting potential maps.

Definition at line 341 of file ComponentComsol.cc.

342 {
343 if (!m_ready) {
344 std::cerr << m_className << "::SetWeightingPotential:\n"
345 << " No valid field map is present.\n"
346 << " Weighting fields cannot be added.\n";
347 return false;
348 }
349
350 std::cout << m_className << "::SetWeightingPotential:\n"
351 << " Reading field map for electrode " << label << ".\n";
352
353 // Check if a weighting field with the same label already exists.
354 if (m_wpot.count(label) > 0) {
355 std::cout << " Replacing existing weighting field.\n";
356 m_wpot[label].clear();
357 }
358
359 std::vector<double> pot(m_nodes.size(), 0.);
360 if (!LoadPotentials(field, pot)) return false;
361 m_wpot[label] = pot;
362 return true;
363}

Referenced by SetWeightingField().

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