Garfield++ v1r0
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 | List of all members
Garfield::ComponentTcad3d Class Reference

#include <ComponentTcad3d.hh>

+ Inheritance diagram for Garfield::ComponentTcad3d:

Public Member Functions

 ComponentTcad3d ()
 
 ~ComponentTcad3d ()
 
void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status)
 
void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status)
 
Medium * GetMedium (const double &x, const double &y, const double &z)
 
bool GetVoltageRange (double &vmin, double &vmax)
 
bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
 
bool Initialise (const std::string gridfilename, const std::string datafilename)
 
void PrintRegions ()
 
int GetNumberOfRegions () const
 
void GetRegion (const int ireg, std::string &name, bool &active)
 
void SetDriftRegion (const int ireg)
 
void UnsetDriftRegion (const int ireg)
 
void SetMedium (const int ireg, Medium *m)
 
bool GetMedium (const int ireg, Medium *&m) const
 
int GetNumberOfElements () const
 
bool GetElement (const int i, double &vol, double &dmin, double &dmax, int &type)
 
bool GetElement (const int i, double &vol, double &dmin, double &dmax, int &type, int &node1, int &node2, int &node3, int &node4, int &node5, int &node6, int &node7, int &reg)
 
int GetNumberOfNodes () const
 
bool GetNode (const int i, double &x, double &y, double &z, double &v, double &ex, double &ey, double &ez)
 
- Public Member Functions inherited from Garfield::ComponentBase
 ComponentBase ()
 
virtual ~ComponentBase ()
 
virtual void SetGeometry (GeometryBase *geo)
 
virtual void Clear ()
 
virtual Medium * GetMedium (const double &x, const double &y, const double &z)
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status)=0
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status)=0
 
virtual bool GetVoltageRange (double &vmin, double &vmax)=0
 
virtual void WeightingField (const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string label)
 
virtual double WeightingPotential (const double x, const double y, const double z, const std::string label)
 
virtual void MagneticField (const double x, const double y, const double z, double &bx, double &by, double &bz, int &status)
 
void SetMagneticField (const double bx, const double by, const double bz)
 
virtual bool IsReady ()
 
virtual bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
 
virtual bool IsWireCrossed (const double x0, const double y0, const double z0, const double x1, const double y1, const double z1, double &xc, double &yc, double &zc)
 
virtual bool IsInTrapRadius (double x0, double y0, double z0, double &xw, double &yw, double &rw)
 
void EnablePeriodicityX ()
 
void DisablePeriodicityX ()
 
void EnablePeriodicityY ()
 
void DisablePeriodicityY ()
 
void EnablePeriodicityZ ()
 
void DisablePeriodicityZ ()
 
void EnableMirrorPeriodicityX ()
 
void DisableMirrorPeriodicityX ()
 
void EnableMirrorPeriodicityY ()
 
void DisableMirrorPeriodicityY ()
 
void EnableMirrorPeriodicityZ ()
 
void DisableMirrorPeriodicityZ ()
 
void EnableAxialPeriodicityX ()
 
void DisableAxialPeriodicityX ()
 
void EnableAxialPeriodicityY ()
 
void DisableAxialPeriodicityY ()
 
void EnableAxialPeriodicityZ ()
 
void DisableAxialPeriodicityZ ()
 
void EnableRotationSymmetryX ()
 
void DisableRotationSymmetryX ()
 
void EnableRotationSymmetryY ()
 
void DisableRotationSymmetryY ()
 
void EnableRotationSymmetryZ ()
 
void DisableRotationSymmetryZ ()
 
void EnableDebugging ()
 
void DisableDebugging ()
 

Additional Inherited Members

virtual void Reset ()=0
 
virtual void UpdatePeriodicity ()=0
 
- Protected Attributes inherited from Garfield::ComponentBase
std::string m_className
 
GeometryBase * theGeometry
 
bool ready
 
bool xPeriodic
 
bool yPeriodic
 
bool zPeriodic
 
bool xMirrorPeriodic
 
bool yMirrorPeriodic
 
bool zMirrorPeriodic
 
bool xAxiallyPeriodic
 
bool yAxiallyPeriodic
 
bool zAxiallyPeriodic
 
bool xRotationSymmetry
 
bool yRotationSymmetry
 
bool zRotationSymmetry
 
double bx0
 
double by0
 
double bz0
 
bool debug
 

Detailed Description

Definition at line 10 of file ComponentTcad3d.hh.

Constructor & Destructor Documentation

◆ ComponentTcad3d()

Garfield::ComponentTcad3d::ComponentTcad3d ( )

Definition at line 13 of file ComponentTcad3d.cc.

14 : ComponentBase(),
15 m_nRegions(0),
16 m_nVertices(0),
17 m_nElements(0),
18 m_lastElement(0) {
19
20 m_className = "ComponentTcad3d";
21
22 m_regions.clear();
23 m_vertices.clear();
24 m_elements.clear();
25
26 for (int i = nMaxVertices; i--;) m_w[i] = 0.;
27}

◆ ~ComponentTcad3d()

Garfield::ComponentTcad3d::~ComponentTcad3d ( )
inline

Definition at line 16 of file ComponentTcad3d.hh.

16{}

Member Function Documentation

◆ ElectricField() [1/2]

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

Implements Garfield::ComponentBase.

Definition at line 29 of file ComponentTcad3d.cc.

32 {
33
34 m = 0;
35 // Make sure the field map has been loaded.
36 if (!ready) {
37 std::cerr << m_className << "::ElectricField:\n";
38 std::cerr << " Field map is not available for interpolation.\n";
39 status = -10;
40 return;
41 }
42
43 // Initialise the electric field and potential.
44 ex = ey = ez = p = 0.;
45
46 double x = xin, y = yin, z = zin;
47 // In case of periodicity, reduce to the cell volume.
48 bool xMirrored = false;
49 const double cellsx = m_xMaxBoundingBox - m_xMinBoundingBox;
50 if (xPeriodic) {
51 x = m_xMinBoundingBox + fmod(x - m_xMinBoundingBox, cellsx);
52 if (x < m_xMinBoundingBox) x += cellsx;
53 } else if (xMirrorPeriodic) {
54 double xNew = m_xMinBoundingBox + fmod(x - m_xMinBoundingBox, cellsx);
55 if (xNew < m_xMinBoundingBox) xNew += cellsx;
56 int nx = int(floor(0.5 + (xNew - x) / cellsx));
57 if (nx != 2 * (nx / 2)) {
58 xNew = m_xMinBoundingBox + m_xMaxBoundingBox - xNew;
59 xMirrored = true;
60 }
61 x = xNew;
62 }
63 bool yMirrored = false;
64 const double cellsy = m_yMaxBoundingBox - m_yMinBoundingBox;
65 if (yPeriodic) {
66 y = m_yMinBoundingBox + fmod(y - m_yMinBoundingBox, cellsy);
67 if (y < m_yMinBoundingBox) y += cellsy;
68 } else if (yMirrorPeriodic) {
69 double yNew = m_yMinBoundingBox + fmod(y - m_yMinBoundingBox, cellsy);
70 if (yNew < m_yMinBoundingBox) yNew += cellsy;
71 int ny = int(floor(0.5 + (yNew - y) / cellsy));
72 if (ny != 2 * (ny / 2)) {
73 yNew = m_yMinBoundingBox + m_yMaxBoundingBox - yNew;
74 yMirrored = true;
75 }
76 y = yNew;
77 }
78 bool zMirrored = false;
79 const double cellsz = m_zMaxBoundingBox - m_zMinBoundingBox;
80 if (zPeriodic) {
81 z = m_zMinBoundingBox + fmod(z - m_zMinBoundingBox, cellsz);
82 if (z < m_zMinBoundingBox) z += cellsz;
83 } else if (zMirrorPeriodic) {
84 double zNew = m_zMinBoundingBox + fmod(z - m_zMinBoundingBox, cellsz);
85 if (zNew < m_zMinBoundingBox) zNew += cellsz;
86 int nz = int(floor(0.5 + (zNew - z) / cellsz));
87 if (nz != 2 * (nz / 2)) {
88 zNew = m_zMinBoundingBox + m_zMaxBoundingBox - zNew;
89 zMirrored = true;
90 }
91 z = zNew;
92 }
93
94 // Check if the point is inside the bounding box.
95 if (x < m_xMinBoundingBox || x > m_xMaxBoundingBox || y < m_yMinBoundingBox ||
96 y > m_yMaxBoundingBox || z < m_zMinBoundingBox || z > m_zMaxBoundingBox) {
97 if (debug) {
98 std::cerr << m_className << "::ElectricField:\n";
99 std::cerr << " Point (" << x << ", " << y << ", " << z
100 << ") is outside the bounding box.\n";
101 }
102 status = -11;
103 return;
104 }
105
106 // Assume this will work.
107 status = 0;
108 // Check if the point is still located in the previously found element.
109 int i = m_lastElement;
110 switch (m_elements[i].type) {
111 case 2:
112 if (CheckTriangle(x, y, z, i)) {
113 ex = m_w[0] * m_vertices[m_elements[i].vertex[0]].ex +
114 m_w[1] * m_vertices[m_elements[i].vertex[1]].ex +
115 m_w[2] * m_vertices[m_elements[i].vertex[2]].ex;
116 ey = m_w[0] * m_vertices[m_elements[i].vertex[0]].ey +
117 m_w[1] * m_vertices[m_elements[i].vertex[1]].ey +
118 m_w[2] * m_vertices[m_elements[i].vertex[2]].ey;
119 ez = m_w[0] * m_vertices[m_elements[i].vertex[0]].ez +
120 m_w[1] * m_vertices[m_elements[i].vertex[1]].ez +
121 m_w[2] * m_vertices[m_elements[i].vertex[2]].ez;
122 p = m_w[0] * m_vertices[m_elements[i].vertex[0]].p +
123 m_w[1] * m_vertices[m_elements[i].vertex[1]].p +
124 m_w[2] * m_vertices[m_elements[i].vertex[2]].p;
125 if (xMirrored) ex = -ex;
126 if (yMirrored) ey = -ey;
127 if (zMirrored) ez = -ez;
128 m = m_regions[m_elements[i].region].medium;
129 if (!m_regions[m_elements[i].region].drift || m == 0) status = -5;
130 return;
131 }
132 break;
133 case 5:
134 if (CheckTetrahedron(x, y, z, i)) {
135 ex = m_w[0] * m_vertices[m_elements[i].vertex[0]].ex +
136 m_w[1] * m_vertices[m_elements[i].vertex[1]].ex +
137 m_w[2] * m_vertices[m_elements[i].vertex[2]].ex +
138 m_w[3] * m_vertices[m_elements[i].vertex[3]].ex;
139 ey = m_w[0] * m_vertices[m_elements[i].vertex[0]].ey +
140 m_w[1] * m_vertices[m_elements[i].vertex[1]].ey +
141 m_w[2] * m_vertices[m_elements[i].vertex[2]].ey +
142 m_w[3] * m_vertices[m_elements[i].vertex[3]].ey;
143 ez = m_w[0] * m_vertices[m_elements[i].vertex[0]].ez +
144 m_w[1] * m_vertices[m_elements[i].vertex[1]].ez +
145 m_w[2] * m_vertices[m_elements[i].vertex[2]].ez +
146 m_w[3] * m_vertices[m_elements[i].vertex[3]].ez;
147 p = m_w[0] * m_vertices[m_elements[i].vertex[0]].p +
148 m_w[1] * m_vertices[m_elements[i].vertex[1]].p +
149 m_w[2] * m_vertices[m_elements[i].vertex[2]].p +
150 m_w[3] * m_vertices[m_elements[i].vertex[3]].p;
151 if (xMirrored) ex = -ex;
152 if (yMirrored) ey = -ey;
153 if (zMirrored) ez = -ez;
154 m = m_regions[m_elements[i].region].medium;
155 if (!m_regions[m_elements[i].region].drift || m == 0) status = -5;
156 return;
157 }
158 break;
159 default:
160 std::cerr << m_className << "::ElectricField:\n";
161 std::cerr << " Unknown element type (" << m_elements[i].type << ").\n";
162 status = -11;
163 return;
164 break;
165 }
166
167 // The point is not in the previous element.
168 // We have to loop over all elements.
169 for (i = m_nElements; i--;) {
170 switch (m_elements[i].type) {
171 case 2:
172 if (CheckTriangle(x, y, z, i)) {
173 ex = m_w[0] * m_vertices[m_elements[i].vertex[0]].ex +
174 m_w[1] * m_vertices[m_elements[i].vertex[1]].ex +
175 m_w[2] * m_vertices[m_elements[i].vertex[2]].ex;
176 ey = m_w[0] * m_vertices[m_elements[i].vertex[0]].ey +
177 m_w[1] * m_vertices[m_elements[i].vertex[1]].ey +
178 m_w[2] * m_vertices[m_elements[i].vertex[2]].ey;
179 ez = m_w[0] * m_vertices[m_elements[i].vertex[0]].ez +
180 m_w[1] * m_vertices[m_elements[i].vertex[1]].ez +
181 m_w[2] * m_vertices[m_elements[i].vertex[2]].ez;
182 p = m_w[0] * m_vertices[m_elements[i].vertex[0]].p +
183 m_w[1] * m_vertices[m_elements[i].vertex[1]].p +
184 m_w[2] * m_vertices[m_elements[i].vertex[2]].p;
185 if (xMirrored) ex = -ex;
186 if (yMirrored) ey = -ey;
187 if (zMirrored) ez = -ez;
188 m_lastElement = i;
189 m = m_regions[m_elements[i].region].medium;
190 if (!m_regions[m_elements[i].region].drift || m == 0) status = -5;
191 return;
192 }
193 break;
194 case 5:
195 if (CheckTetrahedron(x, y, z, i)) {
196 ex = m_w[0] * m_vertices[m_elements[i].vertex[0]].ex +
197 m_w[1] * m_vertices[m_elements[i].vertex[1]].ex +
198 m_w[2] * m_vertices[m_elements[i].vertex[2]].ex +
199 m_w[3] * m_vertices[m_elements[i].vertex[3]].ex;
200 ey = m_w[0] * m_vertices[m_elements[i].vertex[0]].ey +
201 m_w[1] * m_vertices[m_elements[i].vertex[1]].ey +
202 m_w[2] * m_vertices[m_elements[i].vertex[2]].ey +
203 m_w[3] * m_vertices[m_elements[i].vertex[3]].ey;
204 ez = m_w[0] * m_vertices[m_elements[i].vertex[0]].ez +
205 m_w[1] * m_vertices[m_elements[i].vertex[1]].ez +
206 m_w[2] * m_vertices[m_elements[i].vertex[2]].ez +
207 m_w[3] * m_vertices[m_elements[i].vertex[3]].ez;
208 p = m_w[0] * m_vertices[m_elements[i].vertex[0]].p +
209 m_w[1] * m_vertices[m_elements[i].vertex[1]].p +
210 m_w[2] * m_vertices[m_elements[i].vertex[2]].p +
211 m_w[3] * m_vertices[m_elements[i].vertex[3]].p;
212 if (xMirrored) ex = -ex;
213 if (yMirrored) ey = -ey;
214 if (zMirrored) ez = -ez;
215 m_lastElement = i;
216 m = m_regions[m_elements[i].region].medium;
217 if (!m_regions[m_elements[i].region].drift || m == 0) status = -5;
218 return;
219 }
220 break;
221 default:
222 std::cerr << m_className << "::ElectricField:\n";
223 std::cerr << " Invalid element type (" << m_elements[i].type
224 << ").\n";
225 status = -11;
226 return;
227 break;
228 }
229 }
230 // Point is outside the mesh.
231 if (debug) {
232 std::cerr << m_className << "::ElectricField:\n";
233 std::cerr << " Point (" << x << ", " << y << ", " << z
234 << ") is outside the mesh.\n";
235 }
236 status = -6;
237 return;
238}

Referenced by ElectricField().

◆ ElectricField() [2/2]

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

Implements Garfield::ComponentBase.

Definition at line 240 of file ComponentTcad3d.cc.

242 {
243
244 double v = 0.;
245 ElectricField(x, y, z, ex, ey, ez, v, m, status);
246}
Garfield::ComponentTcad3d::ElectricField
void ElectricField(const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status)
Definition: ComponentTcad3d.cc:29

◆ GetBoundingBox()

bool Garfield::ComponentTcad3d::GetBoundingBox ( double &  xmin,
double &  ymin,
double &  zmin,
double &  xmax,
double &  ymax,
double &  zmax 
)
virtual

Reimplemented from Garfield::ComponentBase.

Definition at line 545 of file ComponentTcad3d.cc.

546 {
547
548 if (!ready) return false;
549 xmin = m_xMinBoundingBox;
550 ymin = m_yMinBoundingBox;
551 zmin = m_zMinBoundingBox;
552 xmax = m_xMaxBoundingBox;
553 ymax = m_yMaxBoundingBox;
554 zmax = m_zMaxBoundingBox;
555 if (xPeriodic || xMirrorPeriodic) {
556 xmin = -INFINITY;
557 xmax = +INFINITY;
558 }
559 if (yPeriodic || yMirrorPeriodic) {
560 ymin = -INFINITY;
561 ymax = +INFINITY;
562 }
563 if (zPeriodic || zMirrorPeriodic) {
564 zmin = -INFINITY;
565 zmax = +INFINITY;
566 }
567 return true;
568}

◆ GetElement() [1/2]

bool Garfield::ComponentTcad3d::GetElement ( const int  i,
double &  vol,
double &  dmin,
double &  dmax,
int &  type 
)

Definition at line 672 of file ComponentTcad3d.cc.

673 {
674
675 if (i < 0 || i >= m_nElements) {
676 std::cerr << m_className << "::GetElement:\n";
677 std::cerr << " Element index (" << i << ") out of range.\n";
678 return false;
679 }
680
681 type = m_elements[i].type;
682 if (m_elements[i].type == 2) {
683 // Triangle
684 const double vx = (m_vertices[m_elements[i].vertex[1]].y -
685 m_vertices[m_elements[i].vertex[0]].y) *
686 (m_vertices[m_elements[i].vertex[2]].z -
687 m_vertices[m_elements[i].vertex[0]].z) -
688 (m_vertices[m_elements[i].vertex[1]].z -
689 m_vertices[m_elements[i].vertex[0]].z) *
690 (m_vertices[m_elements[i].vertex[2]].y -
691 m_vertices[m_elements[i].vertex[0]].y);
692 const double vy = (m_vertices[m_elements[i].vertex[1]].z -
693 m_vertices[m_elements[i].vertex[0]].z) *
694 (m_vertices[m_elements[i].vertex[2]].x -
695 m_vertices[m_elements[i].vertex[0]].x) -
696 (m_vertices[m_elements[i].vertex[1]].x -
697 m_vertices[m_elements[i].vertex[0]].x) *
698 (m_vertices[m_elements[i].vertex[2]].z -
699 m_vertices[m_elements[i].vertex[0]].z);
700 const double vz = (m_vertices[m_elements[i].vertex[1]].x -
701 m_vertices[m_elements[i].vertex[0]].x) *
702 (m_vertices[m_elements[i].vertex[2]].y -
703 m_vertices[m_elements[i].vertex[0]].y) -
704 (m_vertices[m_elements[i].vertex[1]].y -
705 m_vertices[m_elements[i].vertex[0]].y) *
706 (m_vertices[m_elements[i].vertex[2]].x -
707 m_vertices[m_elements[i].vertex[0]].x);
708 vol = sqrt(vx * vx + vy * vy + vz * vz);
709 const double a = sqrt(pow(m_vertices[m_elements[i].vertex[1]].x -
710 m_vertices[m_elements[i].vertex[0]].x,
711 2) +
712 pow(m_vertices[m_elements[i].vertex[1]].y -
713 m_vertices[m_elements[i].vertex[0]].y,
714 2) +
715 pow(m_vertices[m_elements[i].vertex[1]].z -
716 m_vertices[m_elements[i].vertex[0]].z,
717 2));
718 const double b = sqrt(pow(m_vertices[m_elements[i].vertex[2]].x -
719 m_vertices[m_elements[i].vertex[0]].x,
720 2) +
721 pow(m_vertices[m_elements[i].vertex[2]].y -
722 m_vertices[m_elements[i].vertex[0]].y,
723 2) +
724 pow(m_vertices[m_elements[i].vertex[2]].z -
725 m_vertices[m_elements[i].vertex[0]].z,
726 2));
727 const double c = sqrt(pow(m_vertices[m_elements[i].vertex[1]].x -
728 m_vertices[m_elements[i].vertex[2]].x,
729 2) +
730 pow(m_vertices[m_elements[i].vertex[1]].y -
731 m_vertices[m_elements[i].vertex[2]].y,
732 2) +
733 pow(m_vertices[m_elements[i].vertex[1]].z -
734 m_vertices[m_elements[i].vertex[2]].z,
735 2));
736 dmin = dmax = a;
737 if (b < dmin) dmin = b;
738 if (c < dmin) dmin = c;
739 if (b > dmax) dmax = b;
740 if (c > dmax) dmax = c;
741 } else if (m_elements[i].type == 5) {
742 // Tetrahedron
743 vol = fabs((m_vertices[m_elements[i].vertex[3]].x -
744 m_vertices[m_elements[i].vertex[0]].x) *
745 ((m_vertices[m_elements[i].vertex[1]].y -
746 m_vertices[m_elements[i].vertex[0]].y) *
747 (m_vertices[m_elements[i].vertex[2]].z -
748 m_vertices[m_elements[i].vertex[0]].z) -
749 (m_vertices[m_elements[i].vertex[2]].y -
750 m_vertices[m_elements[i].vertex[0]].y) *
751 (m_vertices[m_elements[i].vertex[1]].z -
752 m_vertices[m_elements[i].vertex[0]].z)) +
753 (m_vertices[m_elements[i].vertex[3]].y -
754 m_vertices[m_elements[i].vertex[0]].y) *
755 ((m_vertices[m_elements[i].vertex[1]].z -
756 m_vertices[m_elements[i].vertex[0]].z) *
757 (m_vertices[m_elements[i].vertex[2]].x -
758 m_vertices[m_elements[i].vertex[0]].x) -
759 (m_vertices[m_elements[i].vertex[2]].z -
760 m_vertices[m_elements[i].vertex[0]].z) *
761 (m_vertices[m_elements[i].vertex[1]].x -
762 m_vertices[m_elements[i].vertex[0]].x)) +
763 (m_vertices[m_elements[i].vertex[3]].z -
764 m_vertices[m_elements[i].vertex[0]].z) *
765 ((m_vertices[m_elements[i].vertex[1]].x -
766 m_vertices[m_elements[i].vertex[0]].x) *
767 (m_vertices[m_elements[i].vertex[2]].y -
768 m_vertices[m_elements[i].vertex[0]].y) -
769 (m_vertices[m_elements[i].vertex[3]].x -
770 m_vertices[m_elements[i].vertex[0]].x) *
771 (m_vertices[m_elements[i].vertex[1]].y -
772 m_vertices[m_elements[i].vertex[0]].y))) /
773 6.;
774 // Loop over all pairs of m_vertices.
775 for (int j = 0; j < nMaxVertices - 1; ++j) {
776 for (int k = j + 1; k < nMaxVertices; ++k) {
777 // Compute distance.
778 const double dist = sqrt(pow(m_vertices[m_elements[i].vertex[j]].x -
779 m_vertices[m_elements[i].vertex[k]].x,
780 2) +
781 pow(m_vertices[m_elements[i].vertex[j]].y -
782 m_vertices[m_elements[i].vertex[k]].y,
783 2) +
784 pow(m_vertices[m_elements[i].vertex[j]].z -
785 m_vertices[m_elements[i].vertex[k]].z,
786 2));
787 if (k == 1) {
788 dmin = dist;
789 dmax = dist;
790 } else {
791 if (dist < dmin) dmin = dist;
792 if (dist > dmax) dmax = dist;
793 }
794 }
795 }
796 } else {
797 std::cerr << m_className << "::GetElement:\n";
798 std::cerr << " Unexpected element type (" << type << ").\n";
799 return false;
800 }
801 return true;
802}
pow
DoubleAc pow(const DoubleAc &f, double p)
Definition: DoubleAc.cpp:336
sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:313
fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:616

Referenced by GetElement().

◆ GetElement() [2/2]

bool Garfield::ComponentTcad3d::GetElement ( const int  i,
double &  vol,
double &  dmin,
double &  dmax,
int &  type,
int &  node1,
int &  node2,
int &  node3,
int &  node4,
int &  node5,
int &  node6,
int &  node7,
int &  reg 
)

Definition at line 804 of file ComponentTcad3d.cc.

807 {
808
809 if (!GetElement(i, vol, dmin, dmax, type)) return false;
810 node1 = m_elements[i].vertex[0];
811 node2 = m_elements[i].vertex[1];
812 node3 = m_elements[i].vertex[2];
813 node4 = m_elements[i].vertex[3];
814 node5 = m_elements[i].vertex[4];
815 node6 = m_elements[i].vertex[5];
816 node7 = m_elements[i].vertex[6];
817 reg = m_elements[i].region;
818 return true;
819}
Garfield::ComponentTcad3d::GetElement
bool GetElement(const int i, double &vol, double &dmin, double &dmax, int &type)
Definition: ComponentTcad3d.cc:672

◆ GetMedium() [1/2]

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

Reimplemented from Garfield::ComponentBase.

Definition at line 248 of file ComponentTcad3d.cc.

249 {
250
251 // Make sure the field map has been loaded.
252 if (!ready) {
253 std::cerr << m_className << "::GetMedium:\n";
254 std::cerr << " Field map not available for interpolation.\n";
255 return NULL;
256 }
257
258 double x = xin, y = yin, z = zin;
259 // In case of periodicity, reduce to the cell volume.
260 const double cellsx = m_xMaxBoundingBox - m_xMinBoundingBox;
261 if (xPeriodic) {
262 x = m_xMinBoundingBox + fmod(x - m_xMinBoundingBox, cellsx);
263 if (x < m_xMinBoundingBox) x += cellsx;
264 } else if (xMirrorPeriodic) {
265 double xNew = m_xMinBoundingBox + fmod(x - m_xMinBoundingBox, cellsx);
266 if (xNew < m_xMinBoundingBox) xNew += cellsx;
267 int nx = int(floor(0.5 + (xNew - x) / cellsx));
268 if (nx != 2 * (nx / 2)) {
269 xNew = m_xMinBoundingBox + m_xMaxBoundingBox - xNew;
270 }
271 x = xNew;
272 }
273 const double cellsy = m_yMaxBoundingBox - m_yMinBoundingBox;
274 if (yPeriodic) {
275 y = m_yMinBoundingBox + fmod(y - m_yMinBoundingBox, cellsy);
276 if (y < m_yMinBoundingBox) y += cellsy;
277 } else if (yMirrorPeriodic) {
278 double yNew = m_yMinBoundingBox + fmod(y - m_yMinBoundingBox, cellsy);
279 if (yNew < m_yMinBoundingBox) yNew += cellsy;
280 int ny = int(floor(0.5 + (yNew - y) / cellsy));
281 if (ny != 2 * (ny / 2)) {
282 yNew = m_yMinBoundingBox + m_yMaxBoundingBox - yNew;
283 }
284 y = yNew;
285 }
286 const double cellsz = m_zMaxBoundingBox - m_zMinBoundingBox;
287 if (zPeriodic) {
288 z = m_zMinBoundingBox + fmod(z - m_zMinBoundingBox, cellsz);
289 if (z < m_zMinBoundingBox) z += cellsz;
290 } else if (zMirrorPeriodic) {
291 double zNew = m_zMinBoundingBox + fmod(z - m_zMinBoundingBox, cellsz);
292 if (zNew < m_zMinBoundingBox) zNew += cellsz;
293 int nz = int(floor(0.5 + (zNew - z) / cellsz));
294 if (nz != 2 * (nz / 2)) {
295 zNew = m_zMinBoundingBox + m_zMaxBoundingBox - zNew;
296 }
297 z = zNew;
298 }
299
300 // Check if the point is inside the bounding box.
301 if (x < m_xMinBoundingBox || x > m_xMaxBoundingBox || y < m_yMinBoundingBox ||
302 y > m_yMaxBoundingBox || z < m_zMinBoundingBox || z > m_zMaxBoundingBox) {
303 return NULL;
304 }
305
306 // Check if the point is still located in the previous element.
307 int i = m_lastElement;
308 switch (m_elements[i].type) {
309 case 2:
310 if (CheckTriangle(x, y, z, i)) {
311 return m_regions[m_elements[i].region].medium;
312 }
313 break;
314 case 5:
315 if (CheckTetrahedron(x, y, z, i)) {
316 return m_regions[m_elements[i].region].medium;
317 }
318 break;
319 default:
320 std::cerr << m_className << "::GetMedium:\n";
321 std::cerr << " Invalid element type (" << m_elements[i].type << ").\n";
322 return NULL;
323 break;
324 }
325
326 // The point is not in the previous element.
327 // We have to loop over all elements.
328 for (i = m_nElements; i--;) {
329 switch (m_elements[i].type) {
330 case 2:
331 if (CheckTriangle(x, y, z, i)) {
332 m_lastElement = i;
333 return m_regions[m_elements[i].region].medium;
334 }
335 break;
336 case 5:
337 if (CheckTetrahedron(x, y, z, i)) {
338 m_lastElement = i;
339 return m_regions[m_elements[i].region].medium;
340 }
341 break;
342 default:
343 std::cerr << m_className << "::GetMedium:\n";
344 std::cerr << " Invalid element type (" << m_elements[i].type
345 << ").\n";
346 return NULL;
347 break;
348 }
349 }
350 // The point is outside the mesh.
351 return NULL;
352}

◆ GetMedium() [2/2]

bool Garfield::ComponentTcad3d::GetMedium ( const int  ireg,
Medium *&  m 
) const

Definition at line 659 of file ComponentTcad3d.cc.

659 {
660
661 if (i < 0 || i >= m_nRegions) {
662 std::cerr << m_className << "::GetMedium:\n";
663 std::cerr << " Region " << i << " does not exist.\n";
664 return false;
665 }
666
667 m = m_regions[i].medium;
668 if (m == 0) return false;
669 return true;
670}

◆ GetNode()

bool Garfield::ComponentTcad3d::GetNode ( const int  i,
double &  x,
double &  y,
double &  z,
double &  v,
double &  ex,
double &  ey,
double &  ez 
)

Definition at line 821 of file ComponentTcad3d.cc.

822 {
823
824 if (i < 0 || i >= m_nVertices) {
825 std::cerr << m_className << "::GetNode:\n";
826 std::cerr << " Node index (" << i << ") out of range.\n";
827 return false;
828 }
829
830 x = m_vertices[i].x;
831 y = m_vertices[i].y;
832 z = m_vertices[i].z;
833 v = m_vertices[i].p;
834 ex = m_vertices[i].ex;
835 ey = m_vertices[i].ey;
836 ez = m_vertices[i].ez;
837 return true;
838}

◆ GetNumberOfElements()

int Garfield::ComponentTcad3d::GetNumberOfElements ( ) const
inline

Definition at line 45 of file ComponentTcad3d.hh.

45{ return m_nElements; }

◆ GetNumberOfNodes()

int Garfield::ComponentTcad3d::GetNumberOfNodes ( ) const
inline

Definition at line 51 of file ComponentTcad3d.hh.

51{ return m_nVertices; }

◆ GetNumberOfRegions()

int Garfield::ComponentTcad3d::GetNumberOfRegions ( ) const
inline

Definition at line 37 of file ComponentTcad3d.hh.

37{ return m_nRegions; }

◆ GetRegion()

void Garfield::ComponentTcad3d::GetRegion ( const int  ireg,
std::string &  name,
bool &  active 
)

Definition at line 611 of file ComponentTcad3d.cc.

611 {
612
613 if (i < 0 || i >= m_nRegions) {
614 std::cerr << m_className << "::GetRegion:\n";
615 std::cerr << " Region " << i << " does not exist.\n";
616 return;
617 }
618 name = m_regions[i].name;
619 active = m_regions[i].drift;
620}

◆ GetVoltageRange()

bool Garfield::ComponentTcad3d::GetVoltageRange ( double &  vmin,
double &  vmax 
)
virtual

Implements Garfield::ComponentBase.

Definition at line 570 of file ComponentTcad3d.cc.

570 {
571
572 if (!ready) return false;
573 vmin = m_pMin;
574 vmax = m_pMax;
575 return true;
576}

◆ Initialise()

bool Garfield::ComponentTcad3d::Initialise ( const std::string  gridfilename,
const std::string  datafilename 
)

Definition at line 354 of file ComponentTcad3d.cc.

355 {
356
357 ready = false;
358 // Import mesh data from .grd file.
359 if (!LoadGrid(gridfilename)) {
360 std::cerr << m_className << "::Initialise:\n";
361 std::cerr << " Importing mesh data failed.\n";
362 return false;
363 }
364
365 // Import electric field and potential from .dat file.
366 if (!LoadData(datafilename)) {
367 std::cerr << m_className << "::Initialise:\n";
368 std::cerr << " Importing electric field and potential failed.\n";
369 return false;
370 }
371
372 // Find min./max. coordinates and potentials.
373 m_xMaxBoundingBox = m_vertices[m_elements[0].vertex[0]].x;
374 m_yMaxBoundingBox = m_vertices[m_elements[0].vertex[0]].y;
375 m_zMaxBoundingBox = m_vertices[m_elements[0].vertex[0]].z;
376 m_xMinBoundingBox = m_xMaxBoundingBox;
377 m_yMinBoundingBox = m_yMaxBoundingBox;
378 m_zMinBoundingBox = m_zMaxBoundingBox;
379 m_pMax = m_pMin = m_vertices[m_elements[0].vertex[0]].p;
380 for (int i = m_nElements; i--;) {
381 for (int j = 0; j <= m_elements[i].type; ++j) {
382 if (m_vertices[m_elements[i].vertex[j]].x < m_xMinBoundingBox) {
383 m_xMinBoundingBox = m_vertices[m_elements[i].vertex[j]].x;
384 } else if (m_vertices[m_elements[i].vertex[j]].x > m_xMaxBoundingBox) {
385 m_xMaxBoundingBox = m_vertices[m_elements[i].vertex[j]].x;
386 }
387 if (m_vertices[m_elements[i].vertex[j]].y < m_yMinBoundingBox) {
388 m_yMinBoundingBox = m_vertices[m_elements[i].vertex[j]].y;
389 } else if (m_vertices[m_elements[i].vertex[j]].y > m_yMaxBoundingBox) {
390 m_yMaxBoundingBox = m_vertices[m_elements[i].vertex[j]].y;
391 }
392 if (m_vertices[m_elements[i].vertex[j]].z < m_zMinBoundingBox) {
393 m_zMinBoundingBox = m_vertices[m_elements[i].vertex[j]].z;
394 } else if (m_vertices[m_elements[i].vertex[j]].z > m_zMaxBoundingBox) {
395 m_zMaxBoundingBox = m_vertices[m_elements[i].vertex[j]].z;
396 }
397 if (m_vertices[m_elements[i].vertex[j]].p < m_pMin) {
398 m_pMin = m_vertices[m_elements[i].vertex[j]].p;
399 } else if (m_vertices[m_elements[i].vertex[j]].p > m_pMax) {
400 m_pMax = m_vertices[m_elements[i].vertex[j]].p;
401 }
402 }
403 }
404
405 std::cout << m_className << "::Initialise:\n";
406 std::cout << " Bounding box:\n";
407 std::cout << " " << m_xMinBoundingBox << " < x [cm] < "
408 << m_xMaxBoundingBox << "\n";
409 std::cout << " " << m_yMinBoundingBox << " < y [cm] < "
410 << m_yMaxBoundingBox << "\n";
411 std::cout << " " << m_zMinBoundingBox << " < z [cm] < "
412 << m_zMaxBoundingBox << "\n";
413 std::cout << " Voltage range:\n";
414 std::cout << " " << m_pMin << " < V < " << m_pMax << "\n";
415
416 bool ok = true;
417
418 // Count the number of elements belonging to a region.
419 std::vector<int> nElementsRegion;
420 nElementsRegion.resize(m_nRegions);
421 for (int i = m_nRegions; i--;) nElementsRegion[i] = 0;
422
423 // Count the different element shapes.
424 int nTriangles = 0;
425 int nTetrahedra = 0;
426 int nOtherShapes = 0;
427
428 // Check if there are elements which are not part of any region.
429 int nLoose = 0;
430 std::vector<int> looseElements;
431 looseElements.clear();
432
433 // Check if there are degenerate elements.
434 int nDegenerate = 0;
435 std::vector<int> degenerateElements;
436 degenerateElements.clear();
437
438 for (int i = m_nElements; i--;) {
439 if (m_elements[i].type == 2) {
440 ++nTriangles;
441 if (m_elements[i].vertex[0] == m_elements[i].vertex[1] ||
442 m_elements[i].vertex[1] == m_elements[i].vertex[2] ||
443 m_elements[i].vertex[2] == m_elements[i].vertex[0]) {
444 degenerateElements.push_back(i);
445 ++nDegenerate;
446 }
447 } else if (m_elements[i].type == 5) {
448 if (m_elements[i].vertex[0] == m_elements[i].vertex[1] ||
449 m_elements[i].vertex[0] == m_elements[i].vertex[2] ||
450 m_elements[i].vertex[0] == m_elements[i].vertex[3] ||
451 m_elements[i].vertex[1] == m_elements[i].vertex[2] ||
452 m_elements[i].vertex[1] == m_elements[i].vertex[3] ||
453 m_elements[i].vertex[2] == m_elements[i].vertex[3]) {
454 degenerateElements.push_back(i);
455 ++nDegenerate;
456 }
457 ++nTetrahedra;
458 } else {
459 // Other shapes should not occur, since they were excluded in LoadGrid.
460 ++nOtherShapes;
461 }
462 if (m_elements[i].region >= 0 && m_elements[i].region < m_nRegions) {
463 ++nElementsRegion[m_elements[i].region];
464 } else {
465 looseElements.push_back(i);
466 ++nLoose;
467 }
468 }
469
470 if (nDegenerate > 0) {
471 std::cerr << m_className << "::Initialise:\n";
472 std::cerr << " The following elements are degenerate:\n";
473 for (int i = nDegenerate; i--;) {
474 std::cerr << " " << degenerateElements[i] << "\n";
475 }
476 ok = false;
477 }
478
479 if (nLoose > 0) {
480 std::cerr << m_className << "::Initialise:\n";
481 std::cerr << " The following elements are not part of any region:\n";
482 for (int i = nLoose; i--;) {
483 std::cerr << " " << looseElements[i] << "\n";
484 }
485 ok = false;
486 }
487
488 std::cout << m_className << "::Initialise:\n";
489 std::cout << " Number of regions: " << m_nRegions << "\n";
490 for (int i = 0; i < m_nRegions; ++i) {
491 std::cout << " " << i << ": " << m_regions[i].name << ", "
492 << nElementsRegion[i] << " elements\n";
493 }
494
495 std::cout << " Number of elements: " << m_nElements << "\n";
496 if (nTriangles > 0) {
497 std::cout << " " << nTriangles << " triangles\n";
498 }
499 if (nTetrahedra > 0) {
500 std::cout << " " << nTetrahedra << " tetrahedra\n";
501 }
502 if (nOtherShapes > 0) {
503 std::cout << " " << nOtherShapes << " elements of unknown type\n";
504 std::cout << " Program bug!\n";
505 ready = false;
506 Cleanup();
507 return false;
508 }
509 if (debug) {
510 // For each element, print the indices of the constituting vertices.
511 for (int i = 0; i < m_nElements; ++i) {
512 if (m_elements[i].type == 2) {
513 std::cout << " " << i << ": " << m_elements[i].vertex[0] << " "
514 << m_elements[i].vertex[1] << " " << m_elements[i].vertex[2]
515 << " (triangle, region " << m_elements[i].region << ")\n";
516 } else if (m_elements[i].type == 5) {
517 std::cout << " " << i << ": " << m_elements[i].vertex[0] << " "
518 << m_elements[i].vertex[1] << " " << m_elements[i].vertex[2]
519 << " " << m_elements[i].vertex[3] << " (tetrahedron, region "
520 << m_elements[i].region << ")\n";
521 }
522 }
523 }
524
525 std::cout << " Number of vertices: " << m_nVertices << "\n";
526 if (debug) {
527 for (int i = 0; i < m_nVertices; ++i) {
528 std::cout << " " << i << ": (x, y, z) = (" << m_vertices[i].x << ", "
529 << m_vertices[i].y << ", " << m_vertices[i].z
530 << "), V = " << m_vertices[i].p << "\n";
531 }
532 }
533
534 if (!ok) {
535 ready = false;
536 Cleanup();
537 return false;
538 }
539
540 ready = true;
541 UpdatePeriodicity();
542 return true;
543}

◆ PrintRegions()

void Garfield::ComponentTcad3d::PrintRegions ( )

Definition at line 578 of file ComponentTcad3d.cc.

578 {
579
580 // Do not proceed if not properly initialised.
581 if (!ready) {
582 std::cerr << m_className << "::PrintRegions:\n";
583 std::cerr << " Field map not yet initialised.\n";
584 return;
585 }
586
587 if (m_nRegions < 0) {
588 std::cerr << m_className << "::PrintRegions:\n";
589 std::cerr << " No regions are currently defined.\n";
590 return;
591 }
592
593 std::cout << m_className << "::PrintRegions:\n";
594 std::cout << " Currently " << m_nRegions << " regions are defined.\n";
595 std::cout << " Index Name Medium\n";
596 for (int i = 0; i < m_nRegions; ++i) {
597 std::cout << " " << i << " " << m_regions[i].name;
598 if (m_regions[i].medium == 0) {
599 std::cout << " none ";
600 } else {
601 std::cout << " " << m_regions[i].medium->GetName();
602 }
603 if (m_regions[i].drift) {
604 std::cout << " (active region)\n";
605 } else {
606 std::cout << "\n";
607 }
608 }
609}

◆ SetDriftRegion()

void Garfield::ComponentTcad3d::SetDriftRegion ( const int  ireg)

Definition at line 622 of file ComponentTcad3d.cc.

622 {
623
624 if (i < 0 || i >= m_nRegions) {
625 std::cerr << m_className << "::SetDriftRegion:\n";
626 std::cerr << " Region " << i << " does not exist.\n";
627 return;
628 }
629 m_regions[i].drift = true;
630}

◆ SetMedium()

void Garfield::ComponentTcad3d::SetMedium ( const int  ireg,
Medium *  m 
)

Definition at line 642 of file ComponentTcad3d.cc.

642 {
643
644 if (i < 0 || i >= m_nRegions) {
645 std::cerr << m_className << "::SetMedium:\n";
646 std::cerr << " Region " << i << " does not exist.\n";
647 return;
648 }
649
650 if (medium == 0) {
651 std::cerr << m_className << "::SetMedium:\n";
652 std::cerr << " Medium pointer is null.\n";
653 return;
654 }
655
656 m_regions[i].medium = medium;
657}

◆ UnsetDriftRegion()

void Garfield::ComponentTcad3d::UnsetDriftRegion ( const int  ireg)

Definition at line 632 of file ComponentTcad3d.cc.

632 {
633
634 if (i < 0 || i >= m_nRegions) {
635 std::cerr << m_className << "::UnsetDriftRegion:\n";
636 std::cerr << " Region " << i << " does not exist.\n";
637 return;
638 }
639 m_regions[i].drift = false;
640}
The documentation for this class was generated from the following files: