| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Implementation for G4UTet wrapper class
27 //
28 // 1.11.13 G.Cosmo, CERN
29 // -------------------------------------------
30
31 #include "G4Tet.hh"
32 #include "G4UTet.hh"
33
34 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G
35
36 #include "G4AffineTransform.hh"
37 #include "G4VPVParameterisation.hh"
38 #include "G4BoundingEnvelope.hh"
39
40 using namespace CLHEP;
41
42 //////////////////////////////////////////////
43 //
44 // Constructor - create a tetrahedron
45 // This class is implemented separately from g
46 // because the simplex geometry can be compute
47 // which may become important in situations im
48 // in which a very large number of G4Tets are
49 // A Tet has all of its geometrical informatio
50 //
51 G4UTet::G4UTet(const G4String& pName,
52 const G4ThreeVector& anchor,
53 const G4ThreeVector& p1,
54 const G4ThreeVector& p2,
55 const G4ThreeVector& p3, G4bool
56 : Base_t(pName, U3Vector(anchor.x(),anchor.y
57 U3Vector(p1.x(), p1.y(), p1.
58 U3Vector(p2.x(), p2.y(), p2.
59 U3Vector(p3.x(), p3.y(), p3.
60 {
61 // Check for degeneracy
62 G4bool degenerate = CheckDegeneracy(anchor,
63 if(degeneracyFlag != nullptr) *degeneracyFla
64 else if (degenerate)
65 {
66 G4Exception("G4UTet::G4UTet()", "GeomSolid
67 "Degenerate tetrahedron not al
68 }
69
70 // Set bounding box
71 for (G4int i = 0; i < 3; ++i)
72 {
73 fBmin[i] = std::min(std::min(std::min(anch
74 fBmax[i] = std::max(std::max(std::max(anch
75 }
76 }
77
78 //////////////////////////////////////////////
79 //
80 // Fake default constructor - sets only member
81 // for usage restri
82 //
83 G4UTet::G4UTet( __void__& a )
84 : Base_t(a)
85 {
86 }
87
88 //////////////////////////////////////////////
89 //
90 // Destructor
91 //
92 G4UTet::~G4UTet() = default;
93
94 //////////////////////////////////////////////
95 //
96 // Copy constructor
97 //
98 G4UTet::G4UTet(const G4UTet& rhs)
99 : Base_t(rhs)
100 {
101 fBmin = rhs.fBmin;
102 fBmax = rhs.fBmax;
103 }
104
105
106 //////////////////////////////////////////////
107 //
108 // Assignment operator
109 //
110 G4UTet& G4UTet::operator = (const G4UTet& rhs)
111 {
112 // Check assignment to self
113 if (this == &rhs) { return *this; }
114
115 // Copy base class data
116 Base_t::operator=(rhs);
117
118 // Copy bounding box
119 fBmin = rhs.fBmin;
120 fBmax = rhs.fBmax;
121
122 return *this;
123 }
124
125 //////////////////////////////////////////////
126 //
127 // Return true if tetrahedron is degenerate
128 // Tetrahedron is concidered as degenerate in
129 // height is less than the degeneracy toleranc
130 //
131 G4bool G4UTet::CheckDegeneracy(const G4ThreeVe
132 const G4ThreeVe
133 const G4ThreeVe
134 const G4ThreeVe
135 {
136 G4double hmin = 4. * kCarTolerance; // degen
137
138 // Calculate volume
139 G4double vol = std::abs((p1 - p0).cross(p2 -
140
141 // Calculate face areas squared
142 G4double ss[4];
143 ss[0] = ((p1 - p0).cross(p2 - p0)).mag2();
144 ss[1] = ((p2 - p0).cross(p3 - p0)).mag2();
145 ss[2] = ((p3 - p0).cross(p1 - p0)).mag2();
146 ss[3] = ((p2 - p1).cross(p3 - p1)).mag2();
147
148 // Find face with max area
149 G4int k = 0;
150 for (G4int i = 1; i < 4; ++i) { if (ss[i] >
151
152 // Check: vol^2 / s^2 <= hmin^2
153 return (vol*vol <= ss[k]*hmin*hmin);
154 }
155
156 //////////////////////////////////////////////
157 //
158 // Dispatch to parameterisation for replicatio
159 // computation & modification.
160 //
161 void G4UTet::ComputeDimensions(G4VPVParameteri
162 const G4int,
163 const G4VPhysic
164 {
165 }
166
167 //////////////////////////////////////////////
168 //
169 // Make a clone of the object
170 //
171 G4VSolid* G4UTet::Clone() const
172 {
173 return new G4UTet(*this);
174 }
175
176 //////////////////////////////////////////////
177 //
178 // Modifier
179 //
180 void G4UTet::SetVertices(const G4ThreeVector&
181 const G4ThreeVector&
182 const G4ThreeVector&
183 const G4ThreeVector&
184 G4bool* degeneracyFla
185 {
186 // Check for degeneracy
187 G4bool degenerate = CheckDegeneracy(anchor,
188 if(degeneracyFlag != nullptr) *degeneracyFla
189 else if (degenerate)
190 {
191 G4Exception("G4UTet::SetVertices()", "Geom
192 "Degenerate tetrahedron not al
193 }
194
195 // Change tetrahedron
196 *this = G4UTet(GetName(), anchor, p1, p2, p3
197 }
198
199 //////////////////////////////////////////////
200 //
201 // Accessors
202 //
203 void G4UTet::GetVertices(G4ThreeVector& anchor
204 G4ThreeVector& p1,
205 G4ThreeVector& p2,
206 G4ThreeVector& p3) co
207 {
208 std::vector<U3Vector> vec(4);
209 Base_t::GetVertices(vec[0], vec[1], vec[2],
210 anchor = G4ThreeVector(vec[0].x(), vec[0].y(
211 p1 = G4ThreeVector(vec[1].x(), vec[1].y(), v
212 p2 = G4ThreeVector(vec[2].x(), vec[2].y(), v
213 p3 = G4ThreeVector(vec[3].x(), vec[3].y(), v
214 }
215
216 std::vector<G4ThreeVector> G4UTet::GetVertices
217 {
218 std::vector<U3Vector> vec(4);
219 Base_t::GetVertices(vec[0], vec[1], vec[2],
220 std::vector<G4ThreeVector> vertices;
221 for (unsigned int i=0; i<4; ++i)
222 {
223 G4ThreeVector v(vec[i].x(), vec[i].y(), ve
224 vertices.push_back(v);
225 }
226 return vertices;
227 }
228
229 //////////////////////////////////////////////
230 //
231 // Set bounding box
232 //
233 void G4UTet::SetBoundingLimits(const G4ThreeVe
234 const G4ThreeVe
235 {
236 G4ThreeVector fVertex[4];
237 GetVertices(fVertex[0], fVertex[1], fVertex[
238
239 G4int iout[4] = { 0, 0, 0, 0 };
240 for (G4int i = 0; i < 4; ++i)
241 {
242 iout[i] = (G4int)(fVertex[i].x() < pMin.x(
243 fVertex[i].y() < pMin.y(
244 fVertex[i].z() < pMin.z(
245 fVertex[i].x() > pMax.x(
246 fVertex[i].y() > pMax.y(
247 fVertex[i].z() > pMax.z(
248 }
249 if (iout[0] + iout[1] + iout[2] + iout[3] !=
250 {
251 std::ostringstream message;
252 message << "Attempt to set bounding box th
253 << GetName() << " !\n"
254 << " Specified bounding box limit
255 << " pmin: " << pMin << "\n"
256 << " pmax: " << pMax << "\n"
257 << " Tetrahedron vertices:\n"
258 << " anchor " << fVertex[0] <<
259 << " p1 " << fVertex[1] <<
260 << " p2 " << fVertex[2] <<
261 << " p3 " << fVertex[3] <<
262 G4Exception("G4UTet::SetBoundingLimits()",
263 FatalException, message);
264 }
265 fBmin = pMin;
266 fBmax = pMax;
267 }
268
269 //////////////////////////////////////////////
270 //
271 // Get bounding box
272
273 void G4UTet::BoundingLimits(G4ThreeVector& pMi
274 {
275 pMin = fBmin;
276 pMax = fBmax;
277 }
278
279 //////////////////////////////////////////////
280 //
281 // Calculate extent under transform and specif
282
283 G4bool
284 G4UTet::CalculateExtent(const EAxis pAxis,
285 const G4VoxelLimits& p
286 const G4AffineTransfor
287 G4double& pMin,
288 {
289 G4ThreeVector bmin, bmax;
290
291 // Check bounding box (bbox)
292 //
293 BoundingLimits(bmin,bmax);
294 G4BoundingEnvelope bbox(bmin,bmax);
295
296 // Use simple bounding-box to help in the ca
297 //
298 return bbox.CalculateExtent(pAxis,pVoxelLimi
299
300 #if 0
301 // Precise extent computation (disabled by d
302 //
303 G4bool exist;
304 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox
305 {
306 return exist = (pMin < pMax) ? true : fals
307 }
308
309 // Set bounding envelope (benv) and calculat
310 //
311 std::vector<G4ThreeVector> vec = GetVertices
312
313 G4ThreeVectorList anchor(1);
314 anchor[0] = vec[0];
315
316 G4ThreeVectorList base(3);
317 base[0] = vec[1];
318 base[1] = vec[2];
319 base[2] = vec[3];
320
321 std::vector<const G4ThreeVectorList *> polyg
322 polygons[0] = &anchor;
323 polygons[1] = &base;
324
325 G4BoundingEnvelope benv(bmin,bmax,polygons);
326 return exists = benv.CalculateExtent(pAxis,p
327 #endif
328 }
329
330 //////////////////////////////////////////////
331 //
332 // CreatePolyhedron
333 //
334 G4Polyhedron* G4UTet::CreatePolyhedron() const
335 {
336 std::vector<U3Vector> vec(4);
337 Base_t::GetVertices(vec[0], vec[1], vec[2],
338
339 G4double xyz[4][3];
340 const G4int faces[4][4] = {{1,3,2,0},{1,4,3,
341 for (unsigned int i=0; i<4; ++i)
342 {
343 xyz[i][0] = vec[i].x();
344 xyz[i][1] = vec[i].y();
345 xyz[i][2] = vec[i].z();
346 }
347
348 auto ph = new G4Polyhedron;
349 ph->createPolyhedron(4,4,xyz,faces);
350 return ph;
351 }
352
353 #endif // G4GEOM_USE_USOLIDS
354