| 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 <<
27 // 26 //
28 // 1.11.13 G.Cosmo, CERN << 27 //
>> 28 // Implementation for G4UTet wrapper class
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 G4ThreeVector anchor,
53 const G4ThreeVector& p1, << 53 G4ThreeVector p2,
54 const G4ThreeVector& p2, << 54 G4ThreeVector p3,
55 const G4ThreeVector& p3, G4bool << 55 G4ThreeVector p4, 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. <<
58 U3Vector(p2.x(), p2.y(), p2. 57 U3Vector(p2.x(), p2.y(), p2.z()),
59 U3Vector(p3.x(), p3.y(), p3. << 58 U3Vector(p3.x(), p3.y(), p3.z()),
>> 59 U3Vector(p4.x(), p4.y(), p4.z()))
60 { 60 {
61 // Check for degeneracy << 61 G4double fXMin=std::min(std::min(std::min(anchor.x(), p2.x()),p3.x()),p4.x());
62 G4bool degenerate = CheckDegeneracy(anchor, << 62 G4double fXMax=std::max(std::max(std::max(anchor.x(), p2.x()),p3.x()),p4.x());
63 if(degeneracyFlag != nullptr) *degeneracyFla << 63 G4double fYMin=std::min(std::min(std::min(anchor.y(), p2.y()),p3.y()),p4.y());
>> 64 G4double fYMax=std::max(std::max(std::max(anchor.y(), p2.y()),p3.y()),p4.y());
>> 65 G4double fZMin=std::min(std::min(std::min(anchor.z(), p2.z()),p3.z()),p4.z());
>> 66 G4double fZMax=std::max(std::max(std::max(anchor.z(), p2.z()),p3.z()),p4.z());
>> 67
>> 68 G4ThreeVector fMiddle=G4ThreeVector(fXMax+fXMin,fYMax+fYMin,fZMax+fZMin)*0.5;
>> 69 G4double fMaxSize=std::max(std::max(std::max((anchor-fMiddle).mag(),
>> 70 (p2-fMiddle).mag()),
>> 71 (p3-fMiddle).mag()),
>> 72 (p4-fMiddle).mag());
>> 73 // fV<x><y> is vector from vertex <y> to vertex <x>
>> 74 //
>> 75 G4ThreeVector fV21=p2-anchor;
>> 76 G4ThreeVector fV31=p3-anchor;
>> 77 G4ThreeVector fV41=p4-anchor;
>> 78
>> 79 // make sure this is a correctly oriented set of points for the tetrahedron
>> 80 //
>> 81 G4double signed_vol=fV21.cross(fV31).dot(fV41);
>> 82 G4bool degenerate=std::fabs(signed_vol) < 1e-9*fMaxSize*fMaxSize*fMaxSize;
>> 83
>> 84 if(degeneracyFlag) *degeneracyFlag=degenerate;
64 else if (degenerate) 85 else if (degenerate)
65 { 86 {
66 G4Exception("G4UTet::G4UTet()", "GeomSolid 87 G4Exception("G4UTet::G4UTet()", "GeomSolids0002", FatalException,
67 "Degenerate tetrahedron not al 88 "Degenerate tetrahedron not allowed.");
68 } 89 }
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 } 90 }
77 91
78 ////////////////////////////////////////////// 92 //////////////////////////////////////////////////////////////////////////
79 // 93 //
80 // Fake default constructor - sets only member 94 // Fake default constructor - sets only member data and allocates memory
81 // for usage restri 95 // for usage restricted to object persistency.
82 // 96 //
83 G4UTet::G4UTet( __void__& a ) 97 G4UTet::G4UTet( __void__& a )
84 : Base_t(a) 98 : Base_t(a)
85 { 99 {
86 } 100 }
87 101
88 ////////////////////////////////////////////// 102 //////////////////////////////////////////////////////////////////////////
89 // 103 //
90 // Destructor 104 // Destructor
91 // 105 //
92 G4UTet::~G4UTet() = default; << 106 G4UTet::~G4UTet()
>> 107 {
>> 108 }
93 109
94 ////////////////////////////////////////////// 110 ///////////////////////////////////////////////////////////////////////////////
95 // 111 //
96 // Copy constructor 112 // Copy constructor
97 // 113 //
98 G4UTet::G4UTet(const G4UTet& rhs) 114 G4UTet::G4UTet(const G4UTet& rhs)
99 : Base_t(rhs) 115 : Base_t(rhs)
100 { 116 {
101 fBmin = rhs.fBmin; <<
102 fBmax = rhs.fBmax; <<
103 } 117 }
104 118
105 119
106 ////////////////////////////////////////////// 120 ///////////////////////////////////////////////////////////////////////////////
107 // 121 //
108 // Assignment operator 122 // Assignment operator
109 // 123 //
110 G4UTet& G4UTet::operator = (const G4UTet& rhs) << 124 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 { 125 {
136 G4double hmin = 4. * kCarTolerance; // degen << 126 // Check assignment to self
137 << 127 //
138 // Calculate volume << 128 if (this == &rhs) { return *this; }
139 G4double vol = std::abs((p1 - p0).cross(p2 - <<
140 129
141 // Calculate face areas squared << 130 // Copy base class data
142 G4double ss[4]; << 131 //
143 ss[0] = ((p1 - p0).cross(p2 - p0)).mag2(); << 132 Base_t::operator=(rhs);
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 133
152 // Check: vol^2 / s^2 <= hmin^2 << 134 return *this;
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 } 135 }
198 136
199 ////////////////////////////////////////////// 137 ///////////////////////////////////////////////////////////////////////////////
200 // 138 //
201 // Accessors 139 // Accessors
202 // 140 //
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 141 std::vector<G4ThreeVector> G4UTet::GetVertices() const
217 { 142 {
218 std::vector<U3Vector> vec(4); 143 std::vector<U3Vector> vec(4);
219 Base_t::GetVertices(vec[0], vec[1], vec[2], 144 Base_t::GetVertices(vec[0], vec[1], vec[2], vec[3]);
220 std::vector<G4ThreeVector> vertices; 145 std::vector<G4ThreeVector> vertices;
221 for (unsigned int i=0; i<4; ++i) 146 for (unsigned int i=0; i<4; ++i)
222 { 147 {
223 G4ThreeVector v(vec[i].x(), vec[i].y(), ve 148 G4ThreeVector v(vec[i].x(), vec[i].y(), vec[i].z());
224 vertices.push_back(v); 149 vertices.push_back(v);
225 } 150 }
226 return vertices; 151 return vertices;
227 } 152 }
228 153
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 ////////////////////////////////////////////// 154 //////////////////////////////////////////////////////////////////////////
270 // 155 //
271 // Get bounding box 156 // Get bounding box
272 157
273 void G4UTet::BoundingLimits(G4ThreeVector& pMi 158 void G4UTet::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
274 { 159 {
275 pMin = fBmin; << 160 U3Vector vmin, vmax;
276 pMax = fBmax; << 161 Base_t::Extent(vmin,vmax);
>> 162 pMin.set(vmin.x(),vmin.y(),vmin.z());
>> 163 pMax.set(vmax.x(),vmax.y(),vmax.z());
>> 164
>> 165 // Check correctness of the bounding box
>> 166 //
>> 167 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
>> 168 {
>> 169 std::ostringstream message;
>> 170 message << "Bad bounding box (min >= max) for solid: "
>> 171 << GetName() << " !"
>> 172 << "\npMin = " << pMin
>> 173 << "\npMax = " << pMax;
>> 174 G4Exception("G4UTet::BoundingLimits()", "GeomMgt0001",
>> 175 JustWarning, message);
>> 176 StreamInfo(G4cout);
>> 177 }
277 } 178 }
278 179
279 ////////////////////////////////////////////// 180 //////////////////////////////////////////////////////////////////////////
280 // 181 //
281 // Calculate extent under transform and specif 182 // Calculate extent under transform and specified limit
282 183
283 G4bool 184 G4bool
284 G4UTet::CalculateExtent(const EAxis pAxis, 185 G4UTet::CalculateExtent(const EAxis pAxis,
285 const G4VoxelLimits& p 186 const G4VoxelLimits& pVoxelLimit,
286 const G4AffineTransfor 187 const G4AffineTransform& pTransform,
287 G4double& pMin, 188 G4double& pMin, G4double& pMax) const
288 { 189 {
289 G4ThreeVector bmin, bmax; 190 G4ThreeVector bmin, bmax;
>> 191 G4bool exist;
290 192
291 // Check bounding box (bbox) 193 // Check bounding box (bbox)
292 // 194 //
293 BoundingLimits(bmin,bmax); 195 BoundingLimits(bmin,bmax);
294 G4BoundingEnvelope bbox(bmin,bmax); 196 G4BoundingEnvelope bbox(bmin,bmax);
295 << 197 #ifdef G4BBOX_EXTENT
296 // Use simple bounding-box to help in the ca << 198 if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
297 // << 199 #endif
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 200 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
305 { 201 {
306 return exist = (pMin < pMax) ? true : fals 202 return exist = (pMin < pMax) ? true : false;
307 } 203 }
308 204
309 // Set bounding envelope (benv) and calculat 205 // Set bounding envelope (benv) and calculate extent
310 // 206 //
311 std::vector<G4ThreeVector> vec = GetVertices 207 std::vector<G4ThreeVector> vec = GetVertices();
312 208
313 G4ThreeVectorList anchor(1); 209 G4ThreeVectorList anchor(1);
314 anchor[0] = vec[0]; 210 anchor[0] = vec[0];
315 211
316 G4ThreeVectorList base(3); 212 G4ThreeVectorList base(3);
317 base[0] = vec[1]; 213 base[0] = vec[1];
318 base[1] = vec[2]; 214 base[1] = vec[2];
319 base[2] = vec[3]; 215 base[2] = vec[3];
320 216
321 std::vector<const G4ThreeVectorList *> polyg 217 std::vector<const G4ThreeVectorList *> polygons(2);
322 polygons[0] = &anchor; 218 polygons[0] = &anchor;
323 polygons[1] = &base; 219 polygons[1] = &base;
324 220
325 G4BoundingEnvelope benv(bmin,bmax,polygons); 221 G4BoundingEnvelope benv(bmin,bmax,polygons);
326 return exists = benv.CalculateExtent(pAxis,p << 222 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
327 #endif << 223 return exist;
328 } 224 }
329 225
330 ////////////////////////////////////////////// 226 ////////////////////////////////////////////////////////////////////////
331 // 227 //
332 // CreatePolyhedron 228 // CreatePolyhedron
333 // 229 //
334 G4Polyhedron* G4UTet::CreatePolyhedron() const 230 G4Polyhedron* G4UTet::CreatePolyhedron() const
335 { 231 {
336 std::vector<U3Vector> vec(4); << 232 G4int index = 0;
337 Base_t::GetVertices(vec[0], vec[1], vec[2], << 233 G4double array[12];
>> 234 Base_t::GetParametersList(index, array);
338 235
>> 236 G4Polyhedron *ph=new G4Polyhedron;
339 G4double xyz[4][3]; 237 G4double xyz[4][3];
340 const G4int faces[4][4] = {{1,3,2,0},{1,4,3, << 238 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) << 239 xyz[0][0]=array[0]; xyz[0][1]=array[1]; xyz[0][2]=array[2]; // fAnchor
342 { << 240 xyz[1][0]=array[3]; xyz[1][1]=array[4]; xyz[1][2]=array[5]; // fP2
343 xyz[i][0] = vec[i].x(); << 241 xyz[2][0]=array[6]; xyz[2][1]=array[7]; xyz[2][2]=array[8]; // fP3
344 xyz[i][1] = vec[i].y(); << 242 xyz[3][0]=array[9]; xyz[3][1]=array[10]; xyz[3][2]=array[11]; // fP4
345 xyz[i][2] = vec[i].z(); <<
346 } <<
347 243
348 auto ph = new G4Polyhedron; <<
349 ph->createPolyhedron(4,4,xyz,faces); 244 ph->createPolyhedron(4,4,xyz,faces);
>> 245
350 return ph; 246 return ph;
351 } 247 }
352 248
353 #endif // G4GEOM_USE_USOLIDS 249 #endif // G4GEOM_USE_USOLIDS
354 250