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