| 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 intell 18 // * This code implementation is the intellectual property of the *
19 // * Vanderbilt University Free Electron Laser 19 // * Vanderbilt University Free Electron Laser Center *
20 // * Vanderbilt University, Nashville, TN, USA 20 // * Vanderbilt University, Nashville, TN, USA *
21 // * Development supported by: 21 // * Development supported by: *
22 // * United States MFEL program under grant F 22 // * United States MFEL program under grant FA9550-04-1-0045 *
23 // * and NASA under contract number NNG04CT05P 23 // * and NASA under contract number NNG04CT05P *
24 // * Written by Marcus H. Mendenhall and Rober 24 // * Written by Marcus H. Mendenhall and Robert A. Weller. *
25 // * 25 // * *
26 // * Contributed to the Geant4 Core, January, 26 // * Contributed to the Geant4 Core, January, 2005. *
27 // * 27 // * *
28 // ******************************************* 28 // ********************************************************************
29 // 29 //
30 // Implementation for G4Tet class 30 // Implementation for G4Tet class
31 // 31 //
32 // 03.09.2004 - Marcus Mendenhall, created 32 // 03.09.2004 - Marcus Mendenhall, created
33 // 08.01.2020 - Evgueni Tcherniaev, complete r <<
34 // ------------------------------------------- 33 // --------------------------------------------------------------------
35 34
36 #include "G4Tet.hh" 35 #include "G4Tet.hh"
37 36
38 #if !defined(G4GEOM_USE_UTET) 37 #if !defined(G4GEOM_USE_UTET)
39 38
40 #include "G4VoxelLimits.hh" 39 #include "G4VoxelLimits.hh"
41 #include "G4AffineTransform.hh" 40 #include "G4AffineTransform.hh"
42 #include "G4BoundingEnvelope.hh" 41 #include "G4BoundingEnvelope.hh"
43 42
44 #include "G4VPVParameterisation.hh" 43 #include "G4VPVParameterisation.hh"
45 44
46 #include "G4QuickRand.hh" << 45 #include "Randomize.hh"
47 46
48 #include "G4VGraphicsScene.hh" 47 #include "G4VGraphicsScene.hh"
49 #include "G4Polyhedron.hh" 48 #include "G4Polyhedron.hh"
50 #include "G4VisExtent.hh" 49 #include "G4VisExtent.hh"
51 50
>> 51 #include "G4ThreeVector.hh"
>> 52
>> 53 #include <cmath>
>> 54
52 #include "G4AutoLock.hh" 55 #include "G4AutoLock.hh"
53 56
54 namespace 57 namespace
55 { 58 {
56 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE 59 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
57 } 60 }
58 61
59 using namespace CLHEP; 62 using namespace CLHEP;
60 63
61 ////////////////////////////////////////////// 64 ////////////////////////////////////////////////////////////////////////
62 // 65 //
63 // Constructor - create a tetrahedron 66 // Constructor - create a tetrahedron
>> 67 // This class is implemented separately from general polyhedra,
>> 68 // because the simplex geometry can be computed very quickly,
>> 69 // which may become important in situations imported from mesh generators,
>> 70 // in which a very large number of G4Tets are created.
64 // A Tet has all of its geometrical informatio 71 // A Tet has all of its geometrical information precomputed
65 // 72 //
66 G4Tet::G4Tet(const G4String& pName, 73 G4Tet::G4Tet(const G4String& pName,
67 const G4ThreeVector& p0, << 74 G4ThreeVector anchor,
68 const G4ThreeVector& p1, << 75 G4ThreeVector p2,
69 const G4ThreeVector& p2, << 76 G4ThreeVector p3,
70 const G4ThreeVector& p3, G4bool* << 77 G4ThreeVector p4, G4bool* degeneracyFlag)
71 : G4VSolid(pName) 78 : G4VSolid(pName)
72 { 79 {
73 // Check for degeneracy << 80 // fV<x><y> is vector from vertex <y> to vertex <x>
74 G4bool degenerate = CheckDegeneracy(p0, p1, << 81 //
75 if (degeneracyFlag != nullptr) << 82 G4ThreeVector fV21=p2-anchor;
>> 83 G4ThreeVector fV31=p3-anchor;
>> 84 G4ThreeVector fV41=p4-anchor;
>> 85
>> 86 // make sure this is a correctly oriented set of points for the tetrahedron
>> 87 //
>> 88 G4double signed_vol=fV21.cross(fV31).dot(fV41);
>> 89 if(signed_vol<0.0)
76 { 90 {
77 *degeneracyFlag = degenerate; << 91 G4ThreeVector temp(p4);
78 } << 92 p4=p3;
>> 93 p3=temp;
>> 94 temp=fV41;
>> 95 fV41=fV31;
>> 96 fV31=temp;
>> 97 }
>> 98 fCubicVolume = std::fabs(signed_vol) / 6.;
>> 99
>> 100 G4ThreeVector fV24=p2-p4;
>> 101 G4ThreeVector fV43=p4-p3;
>> 102 G4ThreeVector fV32=p3-p2;
>> 103
>> 104 fXMin=std::min(std::min(std::min(anchor.x(), p2.x()),p3.x()),p4.x());
>> 105 fXMax=std::max(std::max(std::max(anchor.x(), p2.x()),p3.x()),p4.x());
>> 106 fYMin=std::min(std::min(std::min(anchor.y(), p2.y()),p3.y()),p4.y());
>> 107 fYMax=std::max(std::max(std::max(anchor.y(), p2.y()),p3.y()),p4.y());
>> 108 fZMin=std::min(std::min(std::min(anchor.z(), p2.z()),p3.z()),p4.z());
>> 109 fZMax=std::max(std::max(std::max(anchor.z(), p2.z()),p3.z()),p4.z());
>> 110
>> 111 fDx=(fXMax-fXMin)*0.5; fDy=(fYMax-fYMin)*0.5; fDz=(fZMax-fZMin)*0.5;
>> 112
>> 113 fMiddle=G4ThreeVector(fXMax+fXMin, fYMax+fYMin, fZMax+fZMin)*0.5;
>> 114 fMaxSize=std::max(std::max(std::max((anchor-fMiddle).mag(),
>> 115 (p2-fMiddle).mag()),
>> 116 (p3-fMiddle).mag()),
>> 117 (p4-fMiddle).mag());
>> 118
>> 119 G4bool degenerate=std::fabs(signed_vol) < 1e-9*fMaxSize*fMaxSize*fMaxSize;
>> 120
>> 121 if(degeneracyFlag) *degeneracyFlag=degenerate;
79 else if (degenerate) 122 else if (degenerate)
80 { 123 {
81 std::ostringstream message; << 124 G4Exception("G4Tet::G4Tet()", "GeomSolids0002", FatalException,
82 message << "Degenerate tetrahedron: " << G << 125 "Degenerate tetrahedron not allowed.");
83 << " anchor: " << p0 << "\n" <<
84 << " p1 : " << p1 << "\n" <<
85 << " p2 : " << p2 << "\n" <<
86 << " p3 : " << p3 << "\n" <<
87 << " volume: " <<
88 << std::abs((p1 - p0).cross(p2 - p <<
89 G4Exception("G4Tet::G4Tet()", "GeomSolids0 <<
90 } 126 }
91 127
92 // Define surface thickness << 128 fTol=1e-9*(std::fabs(fXMin)+std::fabs(fXMax)+std::fabs(fYMin)
93 halfTolerance = 0.5 * kCarTolerance; << 129 +std::fabs(fYMax)+std::fabs(fZMin)+std::fabs(fZMax));
94 << 130 // fTol=kCarTolerance;
95 // Set data members << 131
96 Initialize(p0, p1, p2, p3); << 132 fAnchor=anchor;
>> 133 fP2=p2;
>> 134 fP3=p3;
>> 135 fP4=p4;
>> 136
>> 137 G4ThreeVector fCenter123=(anchor+p2+p3)*(1.0/3.0); // face center
>> 138 G4ThreeVector fCenter134=(anchor+p4+p3)*(1.0/3.0);
>> 139 G4ThreeVector fCenter142=(anchor+p4+p2)*(1.0/3.0);
>> 140 G4ThreeVector fCenter234=(p2+p3+p4)*(1.0/3.0);
>> 141
>> 142 // compute area of each triangular face by cross product
>> 143 // and sum for total surface area
>> 144
>> 145 G4ThreeVector normal123=fV31.cross(fV21);
>> 146 G4ThreeVector normal134=fV41.cross(fV31);
>> 147 G4ThreeVector normal142=fV21.cross(fV41);
>> 148 G4ThreeVector normal234=fV32.cross(fV43);
>> 149
>> 150 fSurfaceArea=(
>> 151 normal123.mag()+
>> 152 normal134.mag()+
>> 153 normal142.mag()+
>> 154 normal234.mag()
>> 155 )/2.0;
>> 156
>> 157 fNormal123=normal123.unit();
>> 158 fNormal134=normal134.unit();
>> 159 fNormal142=normal142.unit();
>> 160 fNormal234=normal234.unit();
>> 161
>> 162 fCdotN123=fCenter123.dot(fNormal123);
>> 163 fCdotN134=fCenter134.dot(fNormal134);
>> 164 fCdotN142=fCenter142.dot(fNormal142);
>> 165 fCdotN234=fCenter234.dot(fNormal234);
97 } 166 }
98 167
99 ////////////////////////////////////////////// << 168 //////////////////////////////////////////////////////////////////////////
100 // 169 //
101 // Fake default constructor - sets only member 170 // Fake default constructor - sets only member data and allocates memory
102 // for usage restri 171 // for usage restricted to object persistency.
103 // 172 //
104 G4Tet::G4Tet( __void__& a ) 173 G4Tet::G4Tet( __void__& a )
105 : G4VSolid(a) << 174 : G4VSolid(a),
>> 175 fAnchor(0,0,0), fP2(0,0,0), fP3(0,0,0), fP4(0,0,0), fMiddle(0,0,0),
>> 176 fNormal123(0,0,0), fNormal142(0,0,0), fNormal134(0,0,0),
>> 177 fNormal234(0,0,0),
>> 178 fCdotN123(0.), fCdotN142(0.), fCdotN134(0.), fCdotN234(0.),
>> 179 fXMin(0.), fXMax(0.), fYMin(0.), fYMax(0.), fZMin(0.), fZMax(0.),
>> 180 fDx(0.), fDy(0.), fDz(0.), fTol(0.), fMaxSize(0.)
106 { 181 {
107 } 182 }
108 183
109 ////////////////////////////////////////////// << 184 //////////////////////////////////////////////////////////////////////////
110 // 185 //
111 // Destructor 186 // Destructor
112 // 187 //
113 G4Tet::~G4Tet() 188 G4Tet::~G4Tet()
114 { 189 {
115 delete fpPolyhedron; fpPolyhedron = nullptr; << 190 delete fpPolyhedron; fpPolyhedron = nullptr;
116 } 191 }
117 192
118 ////////////////////////////////////////////// << 193 ///////////////////////////////////////////////////////////////////////////////
119 // 194 //
120 // Copy constructor 195 // Copy constructor
121 // 196 //
122 G4Tet::G4Tet(const G4Tet& rhs) 197 G4Tet::G4Tet(const G4Tet& rhs)
123 : G4VSolid(rhs) << 198 : G4VSolid(rhs),
>> 199 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
>> 200 fAnchor(rhs.fAnchor),
>> 201 fP2(rhs.fP2), fP3(rhs.fP3), fP4(rhs.fP4), fMiddle(rhs.fMiddle),
>> 202 fNormal123(rhs.fNormal123), fNormal142(rhs.fNormal142),
>> 203 fNormal134(rhs.fNormal134), fNormal234(rhs.fNormal234),
>> 204 warningFlag(rhs.warningFlag), fCdotN123(rhs.fCdotN123),
>> 205 fCdotN142(rhs.fCdotN142), fCdotN134(rhs.fCdotN134),
>> 206 fCdotN234(rhs.fCdotN234), fXMin(rhs.fXMin), fXMax(rhs.fXMax),
>> 207 fYMin(rhs.fYMin), fYMax(rhs.fYMax), fZMin(rhs.fZMin), fZMax(rhs.fZMax),
>> 208 fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz), fTol(rhs.fTol),
>> 209 fMaxSize(rhs.fMaxSize)
124 { 210 {
125 halfTolerance = rhs.halfTolerance; <<
126 fCubicVolume = rhs.fCubicVolume; <<
127 fSurfaceArea = rhs.fSurfaceArea; <<
128 for (G4int i = 0; i < 4; ++i) { fVertex[i] <<
129 for (G4int i = 0; i < 4; ++i) { fNormal[i] <<
130 for (G4int i = 0; i < 4; ++i) { fDist[i] = <<
131 for (G4int i = 0; i < 4; ++i) { fArea[i] = <<
132 fBmin = rhs.fBmin; <<
133 fBmax = rhs.fBmax; <<
134 } 211 }
135 212
136 ////////////////////////////////////////////// << 213
>> 214 ///////////////////////////////////////////////////////////////////////////////
137 // 215 //
138 // Assignment operator 216 // Assignment operator
139 // 217 //
140 G4Tet& G4Tet::operator = (const G4Tet& rhs) << 218 G4Tet& G4Tet::operator = (const G4Tet& rhs)
141 { 219 {
142 // Check assignment to self 220 // Check assignment to self
143 // 221 //
144 if (this == &rhs) { return *this; } 222 if (this == &rhs) { return *this; }
145 223
146 // Copy base class data 224 // Copy base class data
147 // 225 //
148 G4VSolid::operator=(rhs); 226 G4VSolid::operator=(rhs);
149 227
150 // Copy data 228 // Copy data
151 // 229 //
152 halfTolerance = rhs.halfTolerance; << 230 fCubicVolume = rhs.fCubicVolume; fSurfaceArea = rhs.fSurfaceArea;
153 fCubicVolume = rhs.fCubicVolume; << 231 fAnchor = rhs.fAnchor;
154 fSurfaceArea = rhs.fSurfaceArea; << 232 fP2 = rhs.fP2; fP3 = rhs.fP3; fP4 = rhs.fP4; fMiddle = rhs.fMiddle;
155 for (G4int i = 0; i < 4; ++i) { fVertex[i] << 233 fNormal123 = rhs.fNormal123; fNormal142 = rhs.fNormal142;
156 for (G4int i = 0; i < 4; ++i) { fNormal[i] << 234 fNormal134 = rhs.fNormal134; fNormal234 = rhs.fNormal234;
157 for (G4int i = 0; i < 4; ++i) { fDist[i] = << 235 warningFlag = rhs.warningFlag; fCdotN123 = rhs.fCdotN123;
158 for (G4int i = 0; i < 4; ++i) { fArea[i] = << 236 fCdotN142 = rhs.fCdotN142; fCdotN134 = rhs.fCdotN134;
159 fBmin = rhs.fBmin; << 237 fCdotN234 = rhs.fCdotN234; fXMin = rhs.fXMin; fXMax = rhs.fXMax;
160 fBmax = rhs.fBmax; << 238 fYMin = rhs.fYMin; fYMax = rhs.fYMax; fZMin = rhs.fZMin; fZMax = rhs.fZMax;
>> 239 fDx = rhs.fDx; fDy = rhs.fDy; fDz = rhs.fDz; fTol = rhs.fTol;
>> 240 fMaxSize = rhs.fMaxSize;
161 fRebuildPolyhedron = false; 241 fRebuildPolyhedron = false;
162 delete fpPolyhedron; fpPolyhedron = nullptr 242 delete fpPolyhedron; fpPolyhedron = nullptr;
163 243
164 return *this; 244 return *this;
165 } 245 }
166 246
167 ////////////////////////////////////////////// << 247 //////////////////////////////////////////////////////////////////////////
168 // 248 //
169 // Return true if tetrahedron is degenerate << 249 // CheckDegeneracy
170 // Tetrahedron is concidered as degenerate in <<
171 // height is less than degeneracy tolerance <<
172 // <<
173 G4bool G4Tet::CheckDegeneracy(const G4ThreeVec <<
174 const G4ThreeVec <<
175 const G4ThreeVec <<
176 const G4ThreeVec <<
177 { <<
178 G4double hmin = 4. * kCarTolerance; // degen <<
179 <<
180 // Calculate volume <<
181 G4double vol = std::abs((p1 - p0).cross(p2 - <<
182 <<
183 // Calculate face areas squared <<
184 G4double ss[4]; <<
185 ss[0] = ((p1 - p0).cross(p2 - p0)).mag2(); <<
186 ss[1] = ((p2 - p0).cross(p3 - p0)).mag2(); <<
187 ss[2] = ((p3 - p0).cross(p1 - p0)).mag2(); <<
188 ss[3] = ((p2 - p1).cross(p3 - p1)).mag2(); <<
189 <<
190 // Find face with max area <<
191 G4int k = 0; <<
192 for (G4int i = 1; i < 4; ++i) { if (ss[i] > <<
193 <<
194 // Check: vol^2 / s^2 <= hmin^2 <<
195 return (vol*vol <= ss[k]*hmin*hmin); <<
196 } <<
197 <<
198 ////////////////////////////////////////////// <<
199 // <<
200 // Set data members <<
201 // 250 //
202 void G4Tet::Initialize(const G4ThreeVector& p0 << 251 G4bool G4Tet::CheckDegeneracy( G4ThreeVector anchor,
203 const G4ThreeVector& p1 << 252 G4ThreeVector p2,
204 const G4ThreeVector& p2 << 253 G4ThreeVector p3,
205 const G4ThreeVector& p3 << 254 G4ThreeVector p4 )
206 { 255 {
207 // Set vertices << 256 G4bool result;
208 fVertex[0] = p0; << 257 G4Tet* object=new G4Tet("temp",anchor,p2,p3,p4,&result);
209 fVertex[1] = p1; << 258 delete object;
210 fVertex[2] = p2; << 259 return result;
211 fVertex[3] = p3; <<
212 <<
213 G4ThreeVector norm[4]; <<
214 norm[0] = (p2 - p0).cross(p1 - p0); <<
215 norm[1] = (p3 - p0).cross(p2 - p0); <<
216 norm[2] = (p1 - p0).cross(p3 - p0); <<
217 norm[3] = (p2 - p1).cross(p3 - p1); <<
218 G4double volume = norm[0].dot(p3 - p0); <<
219 if (volume > 0.) <<
220 { <<
221 for (auto & i : norm) { i = -i; } <<
222 } <<
223 <<
224 // Set normals to face planes <<
225 for (G4int i = 0; i < 4; ++i) { fNormal[i] = <<
226 <<
227 // Set distances to planes <<
228 for (G4int i = 0; i < 3; ++i) { fDist[i] = f <<
229 fDist[3] = fNormal[3].dot(p1); <<
230 <<
231 // Set face areas <<
232 for (G4int i = 0; i < 4; ++i) { fArea[i] = 0 <<
233 <<
234 // Set bounding box <<
235 for (G4int i = 0; i < 3; ++i) <<
236 { <<
237 fBmin[i] = std::min(std::min(std::min(p0[i <<
238 fBmax[i] = std::max(std::max(std::max(p0[i <<
239 } <<
240 <<
241 // Set volume and surface area <<
242 fCubicVolume = std::abs(volume)/6.; <<
243 fSurfaceArea = fArea[0] + fArea[1] + fArea[2 <<
244 } <<
245 <<
246 ////////////////////////////////////////////// <<
247 // <<
248 // Set vertices <<
249 // <<
250 void G4Tet::SetVertices(const G4ThreeVector& p <<
251 const G4ThreeVector& p <<
252 const G4ThreeVector& p <<
253 const G4ThreeVector& p <<
254 { <<
255 // Check for degeneracy <<
256 G4bool degenerate = CheckDegeneracy(p0, p1, <<
257 if (degeneracyFlag != nullptr) <<
258 { <<
259 *degeneracyFlag = degenerate; <<
260 } <<
261 else if (degenerate) <<
262 { <<
263 std::ostringstream message; <<
264 message << "Degenerate tetrahedron is not <<
265 << " anchor: " << p0 << "\n" <<
266 << " p1 : " << p1 << "\n" <<
267 << " p2 : " << p2 << "\n" <<
268 << " p3 : " << p3 << "\n" <<
269 << " volume: " <<
270 << std::abs((p1 - p0).cross(p2 - p <<
271 G4Exception("G4Tet::SetVertices()", "GeomS <<
272 FatalException, message); <<
273 } <<
274 <<
275 // Set data members <<
276 Initialize(p0, p1, p2, p3); <<
277 <<
278 // Set flag to rebuild polyhedron <<
279 fRebuildPolyhedron = true; <<
280 } 260 }
281 261
282 ////////////////////////////////////////////// << 262 //////////////////////////////////////////////////////////////////////////
283 // <<
284 // Return four vertices <<
285 // <<
286 void G4Tet::GetVertices(G4ThreeVector& p0, <<
287 G4ThreeVector& p1, <<
288 G4ThreeVector& p2, <<
289 G4ThreeVector& p3) con <<
290 { <<
291 p0 = fVertex[0]; <<
292 p1 = fVertex[1]; <<
293 p2 = fVertex[2]; <<
294 p3 = fVertex[3]; <<
295 } <<
296 <<
297 ////////////////////////////////////////////// <<
298 // <<
299 // Return std::vector of vertices <<
300 // <<
301 std::vector<G4ThreeVector> G4Tet::GetVertices( <<
302 { <<
303 std::vector<G4ThreeVector> vertices(4); <<
304 for (G4int i = 0; i < 4; ++i) { vertices[i] <<
305 return vertices; <<
306 } <<
307 <<
308 ////////////////////////////////////////////// <<
309 // 263 //
310 // Dispatch to parameterisation for replicatio 264 // Dispatch to parameterisation for replication mechanism dimension
311 // computation & modification. 265 // computation & modification.
312 // 266 //
313 void G4Tet::ComputeDimensions(G4VPVParameteris 267 void G4Tet::ComputeDimensions(G4VPVParameterisation* ,
314 const G4int , 268 const G4int ,
315 const G4VPhysica 269 const G4VPhysicalVolume* )
316 { 270 {
>> 271 G4Exception("G4Tet::ComputeDimensions()",
>> 272 "GeomSolids0001", FatalException,
>> 273 "G4Tet does not support Parameterisation.");
317 } 274 }
318 275
319 ////////////////////////////////////////////// << 276 //////////////////////////////////////////////////////////////////////////
320 // 277 //
321 // Set bounding box << 278 // Get bounding box
322 // 279 //
323 void G4Tet::SetBoundingLimits(const G4ThreeVec << 280 void G4Tet::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
324 const G4ThreeVec <<
325 { 281 {
326 G4int iout[4] = { 0, 0, 0, 0 }; << 282 pMin.set(fXMin,fYMin,fZMin);
327 for (G4int i = 0; i < 4; ++i) << 283 pMax.set(fXMax,fYMax,fZMax);
328 { << 284
329 iout[i] = (G4int)(fVertex[i].x() < pMin.x( << 285 // Check correctness of the bounding box
330 fVertex[i].y() < pMin.y( << 286 //
331 fVertex[i].z() < pMin.z( << 287 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
332 fVertex[i].x() > pMax.x( <<
333 fVertex[i].y() > pMax.y( <<
334 fVertex[i].z() > pMax.z( <<
335 } <<
336 if (iout[0] + iout[1] + iout[2] + iout[3] != <<
337 { 288 {
338 std::ostringstream message; 289 std::ostringstream message;
339 message << "Attempt to set bounding box th << 290 message << "Bad bounding box (min >= max) for solid: "
340 << GetName() << " !\n" << 291 << GetName() << " !"
341 << " Specified bounding box limit << 292 << "\npMin = " << pMin
342 << " pmin: " << pMin << "\n" << 293 << "\npMax = " << pMax;
343 << " pmax: " << pMax << "\n" << 294 G4Exception("G4Tet::BoundingLimits()", "GeomMgt0001", JustWarning, message);
344 << " Tetrahedron vertices:\n" << 295 DumpInfo();
345 << " anchor " << fVertex[0] << <<
346 << " p1 " << fVertex[1] << <<
347 << " p2 " << fVertex[2] << <<
348 << " p3 " << fVertex[3] << <<
349 G4Exception("G4Tet::SetBoundingLimits()", <<
350 FatalException, message); <<
351 } 296 }
352 fBmin = pMin; <<
353 fBmax = pMax; <<
354 } <<
355 <<
356 ////////////////////////////////////////////// <<
357 // <<
358 // Return bounding box <<
359 // <<
360 void G4Tet::BoundingLimits(G4ThreeVector& pMin <<
361 { <<
362 pMin = fBmin; <<
363 pMax = fBmax; <<
364 } 297 }
365 298
366 ////////////////////////////////////////////// << 299 //////////////////////////////////////////////////////////////////////////
367 // 300 //
368 // Calculate extent under transform and specif 301 // Calculate extent under transform and specified limit
369 // 302 //
370 G4bool G4Tet::CalculateExtent(const EAxis pAxi 303 G4bool G4Tet::CalculateExtent(const EAxis pAxis,
371 const G4VoxelLim 304 const G4VoxelLimits& pVoxelLimit,
372 const G4AffineTr 305 const G4AffineTransform& pTransform,
373 G4double& 306 G4double& pMin, G4double& pMax) const
374 { 307 {
375 G4ThreeVector bmin, bmax; 308 G4ThreeVector bmin, bmax;
376 309
377 // Check bounding box (bbox) 310 // Check bounding box (bbox)
378 // 311 //
379 BoundingLimits(bmin,bmax); 312 BoundingLimits(bmin,bmax);
380 G4BoundingEnvelope bbox(bmin,bmax); 313 G4BoundingEnvelope bbox(bmin,bmax);
381 314
382 // Use simple bounding-box to help in the ca 315 // Use simple bounding-box to help in the case of complex 3D meshes
383 // 316 //
384 return bbox.CalculateExtent(pAxis,pVoxelLimi 317 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
385 318
386 #if 0 319 #if 0
387 // Precise extent computation (disabled by d 320 // Precise extent computation (disabled by default for this shape)
388 // 321 //
389 G4bool exist; 322 G4bool exist;
390 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox 323 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
391 { 324 {
392 return exist = (pMin < pMax) ? true : fals 325 return exist = (pMin < pMax) ? true : false;
393 } 326 }
394 327
395 // Set bounding envelope (benv) and calculat 328 // Set bounding envelope (benv) and calculate extent
396 // 329 //
397 std::vector<G4ThreeVector> vec = GetVertices 330 std::vector<G4ThreeVector> vec = GetVertices();
398 331
399 G4ThreeVectorList anchor(1); 332 G4ThreeVectorList anchor(1);
400 anchor[0].set(vec[0].x(),vec[0].y(),vec[0].z 333 anchor[0].set(vec[0].x(),vec[0].y(),vec[0].z());
401 334
402 G4ThreeVectorList base(3); 335 G4ThreeVectorList base(3);
403 base[0].set(vec[1].x(),vec[1].y(),vec[1].z() 336 base[0].set(vec[1].x(),vec[1].y(),vec[1].z());
404 base[1].set(vec[2].x(),vec[2].y(),vec[2].z() 337 base[1].set(vec[2].x(),vec[2].y(),vec[2].z());
405 base[2].set(vec[3].x(),vec[3].y(),vec[3].z() 338 base[2].set(vec[3].x(),vec[3].y(),vec[3].z());
406 339
407 std::vector<const G4ThreeVectorList *> polyg 340 std::vector<const G4ThreeVectorList *> polygons(2);
408 polygons[0] = &anchor; 341 polygons[0] = &anchor;
409 polygons[1] = &base; 342 polygons[1] = &base;
410 343
411 G4BoundingEnvelope benv(bmin,bmax,polygons); 344 G4BoundingEnvelope benv(bmin,bmax,polygons);
412 return exists = benv.CalculateExtent(pAxis,p 345 return exists = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
413 #endif 346 #endif
414 } 347 }
415 348
416 ////////////////////////////////////////////// << 349 /////////////////////////////////////////////////////////////////////////
417 // 350 //
418 // Return whether point inside/outside/on surf << 351 // Return whether point inside/outside/on surface, using tolerance
419 // 352 //
420 EInside G4Tet::Inside(const G4ThreeVector& p) 353 EInside G4Tet::Inside(const G4ThreeVector& p) const
421 { 354 {
422 G4double dd[4]; << 355 G4double r123, r134, r142, r234;
423 for (G4int i = 0; i < 4; ++i) { dd[i] = fNor << 356
>> 357 // this is written to allow if-statement truncation so the outside test
>> 358 // (where most of the world is) can fail very quickly and efficiently
424 359
425 G4double dist = std::max(std::max(std::max(d << 360 if ( (r123=p.dot(fNormal123)-fCdotN123) > fTol ||
426 return (dist <= -halfTolerance) ? << 361 (r134=p.dot(fNormal134)-fCdotN134) > fTol ||
427 kInside : ((dist <= halfTolerance) ? kSurf << 362 (r142=p.dot(fNormal142)-fCdotN142) > fTol ||
>> 363 (r234=p.dot(fNormal234)-fCdotN234) > fTol )
>> 364 {
>> 365 return kOutside; // at least one is out!
>> 366 }
>> 367 else if( (r123 < -fTol)&&(r134 < -fTol)&&(r142 < -fTol)&&(r234 < -fTol) )
>> 368 {
>> 369 return kInside; // all are definitively inside
>> 370 }
>> 371 else
>> 372 {
>> 373 return kSurface; // too close to tell
>> 374 }
428 } 375 }
429 376
430 ////////////////////////////////////////////// << 377 ///////////////////////////////////////////////////////////////////////
431 // 378 //
432 // Return unit normal to surface at p << 379 // Calculate side nearest to p, and return normal
>> 380 // If two sides are equidistant, normal of first side (x/y/z)
>> 381 // encountered returned.
>> 382 // This assumes that we are looking from the inside!
433 // 383 //
434 G4ThreeVector G4Tet::SurfaceNormal( const G4Th 384 G4ThreeVector G4Tet::SurfaceNormal( const G4ThreeVector& p) const
435 { 385 {
436 G4double k[4]; << 386 G4double r123=std::fabs(p.dot(fNormal123)-fCdotN123);
437 for (G4int i = 0; i < 4; ++i) << 387 G4double r134=std::fabs(p.dot(fNormal134)-fCdotN134);
>> 388 G4double r142=std::fabs(p.dot(fNormal142)-fCdotN142);
>> 389 G4double r234=std::fabs(p.dot(fNormal234)-fCdotN234);
>> 390
>> 391 const G4double delta = 0.5*kCarTolerance;
>> 392 G4ThreeVector sumnorm(0., 0., 0.);
>> 393 G4int noSurfaces=0;
>> 394
>> 395 if (r123 <= delta)
438 { 396 {
439 k[i] = (G4double)(std::abs(fNormal[i].dot( << 397 ++noSurfaces;
>> 398 sumnorm= fNormal123;
440 } 399 }
441 G4double nsurf = k[0] + k[1] + k[2] + k[3]; <<
442 G4ThreeVector norm = <<
443 k[0]*fNormal[0] + k[1]*fNormal[1] + k[2]*f <<
444 400
445 if (nsurf == 1.) return norm; << 401 if (r134 <= delta)
446 else if (nsurf > 1.) return norm.unit(); // <<
447 { 402 {
448 #ifdef G4SPECSDEBUG << 403 ++noSurfaces;
449 std::ostringstream message; << 404 sumnorm += fNormal134;
450 G4long oldprc = message.precision(16); <<
451 message << "Point p is not on surface (!?) <<
452 << GetName() << "\n"; <<
453 message << "Position:\n"; <<
454 message << " p.x() = " << p.x()/mm << " <<
455 message << " p.y() = " << p.y()/mm << " <<
456 message << " p.z() = " << p.z()/mm << " <<
457 G4cout.precision(oldprc); <<
458 G4Exception("G4Tet::SurfaceNormal(p)", "Ge <<
459 JustWarning, message ); <<
460 DumpInfo(); <<
461 #endif <<
462 return ApproxSurfaceNormal(p); <<
463 } 405 }
464 } << 406
465 << 407 if (r142 <= delta)
466 ////////////////////////////////////////////// <<
467 // <<
468 // Find surface nearest to point and return co <<
469 // This method normally should not be called <<
470 // <<
471 G4ThreeVector G4Tet::ApproxSurfaceNormal(const <<
472 { <<
473 G4double dist = -DBL_MAX; <<
474 G4int iside = 0; <<
475 for (G4int i = 0; i < 4; ++i) <<
476 { 408 {
477 G4double d = fNormal[i].dot(p) - fDist[i]; << 409 ++noSurfaces;
478 if (d > dist) { dist = d; iside = i; } << 410 sumnorm += fNormal142;
479 } 411 }
480 return fNormal[iside]; << 412 if (r234 <= delta)
481 } << 413 {
>> 414 ++noSurfaces;
>> 415 sumnorm += fNormal234;
>> 416 }
>> 417
>> 418 if( noSurfaces > 0 )
>> 419 {
>> 420 if( noSurfaces == 1 )
>> 421 {
>> 422 return sumnorm;
>> 423 }
>> 424 else
>> 425 {
>> 426 return sumnorm.unit();
>> 427 }
>> 428 }
>> 429 else // Approximative Surface Normal
>> 430 {
482 431
483 ////////////////////////////////////////////// << 432 if( (r123<=r134) && (r123<=r142) && (r123<=r234) ) { return fNormal123; }
>> 433 else if ( (r134<=r142) && (r134<=r234) ) { return fNormal134; }
>> 434 else if (r142 <= r234) { return fNormal142; }
>> 435 return fNormal234;
>> 436 }
>> 437 }
>> 438 ///////////////////////////////////////////////////////////////////////////
484 // 439 //
485 // Calculate distance to surface from outside, << 440 // Calculate distance to box from an outside point
486 // return kInfinity if no intersection << 441 // - return kInfinity if no intersection.
>> 442 // All this is very unrolled, for speed.
487 // 443 //
488 G4double G4Tet::DistanceToIn(const G4ThreeVect 444 G4double G4Tet::DistanceToIn(const G4ThreeVector& p,
489 const G4ThreeVect 445 const G4ThreeVector& v) const
490 { 446 {
491 G4double tin = -DBL_MAX, tout = DBL_MAX; << 447 G4ThreeVector vu(v.unit()), hp;
492 for (G4int i = 0; i < 4; ++i) << 448 G4double vdotn, t, tmin=kInfinity;
493 { << 449
494 G4double cosa = fNormal[i].dot(v); << 450 G4double extraDistance=10.0*fTol; // a little ways into the solid
495 G4double dist = fNormal[i].dot(p) - fDist[ << 451
496 if (dist >= -halfTolerance) << 452 vdotn=-vu.dot(fNormal123);
497 { << 453 if(vdotn > 1e-12)
498 if (cosa >= 0.) { return kInfinity; } << 454 { // this is a candidate face, since it is pointing at us
499 tin = std::max(tin, -dist/cosa); << 455 t=(p.dot(fNormal123)-fCdotN123)/vdotn; // # distance to intersection
>> 456 if( (t>=-fTol) && (t<tmin) )
>> 457 { // if not true, we're going away from this face or it's not close
>> 458 hp=p+vu*(t+extraDistance); // a little beyond point of intersection
>> 459 if ( ( hp.dot(fNormal134)-fCdotN134 < 0.0 ) &&
>> 460 ( hp.dot(fNormal142)-fCdotN142 < 0.0 ) &&
>> 461 ( hp.dot(fNormal234)-fCdotN234 < 0.0 ) )
>> 462 {
>> 463 tmin=t;
>> 464 }
>> 465 }
500 } 466 }
501 else if (cosa > 0.) << 467
502 { << 468 vdotn=-vu.dot(fNormal134);
503 tout = std::min(tout, -dist/cosa); << 469 if(vdotn > 1e-12)
>> 470 { // # this is a candidate face, since it is pointing at us
>> 471 t=(p.dot(fNormal134)-fCdotN134)/vdotn; // # distance to intersection
>> 472 if( (t>=-fTol) && (t<tmin) )
>> 473 { // if not true, we're going away from this face
>> 474 hp=p+vu*(t+extraDistance); // a little beyond point of intersection
>> 475 if ( ( hp.dot(fNormal123)-fCdotN123 < 0.0 ) &&
>> 476 ( hp.dot(fNormal142)-fCdotN142 < 0.0 ) &&
>> 477 ( hp.dot(fNormal234)-fCdotN234 < 0.0 ) )
>> 478 {
>> 479 tmin=t;
>> 480 }
>> 481 }
>> 482 }
>> 483
>> 484 vdotn=-vu.dot(fNormal142);
>> 485 if(vdotn > 1e-12)
>> 486 { // # this is a candidate face, since it is pointing at us
>> 487 t=(p.dot(fNormal142)-fCdotN142)/vdotn; // # distance to intersection
>> 488 if( (t>=-fTol) && (t<tmin) )
>> 489 { // if not true, we're going away from this face
>> 490 hp=p+vu*(t+extraDistance); // a little beyond point of intersection
>> 491 if ( ( hp.dot(fNormal123)-fCdotN123 < 0.0 ) &&
>> 492 ( hp.dot(fNormal134)-fCdotN134 < 0.0 ) &&
>> 493 ( hp.dot(fNormal234)-fCdotN234 < 0.0 ) )
>> 494 {
>> 495 tmin=t;
>> 496 }
>> 497 }
>> 498 }
>> 499
>> 500 vdotn=-vu.dot(fNormal234);
>> 501 if(vdotn > 1e-12)
>> 502 { // # this is a candidate face, since it is pointing at us
>> 503 t=(p.dot(fNormal234)-fCdotN234)/vdotn; // # distance to intersection
>> 504 if( (t>=-fTol) && (t<tmin) )
>> 505 { // if not true, we're going away from this face
>> 506 hp=p+vu*(t+extraDistance); // a little beyond point of intersection
>> 507 if ( ( hp.dot(fNormal123)-fCdotN123 < 0.0 ) &&
>> 508 ( hp.dot(fNormal134)-fCdotN134 < 0.0 ) &&
>> 509 ( hp.dot(fNormal142)-fCdotN142 < 0.0 ) )
>> 510 {
>> 511 tmin=t;
>> 512 }
>> 513 }
504 } 514 }
505 } <<
506 515
507 return (tout - tin <= halfTolerance) ? << 516 return std::max(0.0,tmin);
508 kInfinity : ((tin < halfTolerance) ? 0. : <<
509 } 517 }
510 518
511 ////////////////////////////////////////////// << 519 //////////////////////////////////////////////////////////////////////////
512 // << 520 //
513 // Estimate safety distance to surface from ou << 521 // Approximate distance to tet.
>> 522 // returns distance to sphere centered on bounding box
>> 523 // - If inside return 0
514 // 524 //
515 G4double G4Tet::DistanceToIn(const G4ThreeVect 525 G4double G4Tet::DistanceToIn(const G4ThreeVector& p) const
516 { 526 {
517 G4double dd[4]; << 527 G4double dd=(p-fMiddle).mag() - fMaxSize - fTol;
518 for (G4int i = 0; i < 4; ++i) { dd[i] = fNor << 528 return std::max(0.0, dd);
519 <<
520 G4double dist = std::max(std::max(std::max(d <<
521 return (dist > 0.) ? dist : 0.; <<
522 } 529 }
523 530
524 ////////////////////////////////////////////// << 531 /////////////////////////////////////////////////////////////////////////
525 // 532 //
526 // Calcluate distance to surface from inside << 533 // Calcluate distance to surface of box from inside
>> 534 // by calculating distances to box's x/y/z planes.
>> 535 // Smallest distance is exact distance to exiting.
527 // 536 //
528 G4double G4Tet::DistanceToOut(const G4ThreeVec << 537 G4double G4Tet::DistanceToOut( const G4ThreeVector& p,const G4ThreeVector& v,
529 const G4ThreeVec << 538 const G4bool calcNorm,
530 const G4bool cal << 539 G4bool* validNorm, G4ThreeVector* n) const
531 G4bool* va << 540 {
532 G4ThreeVec << 541 G4ThreeVector vu(v.unit());
533 { << 542 G4double t1=kInfinity,t2=kInfinity,t3=kInfinity,t4=kInfinity, vdotn, tt;
534 // Calculate distances and cosines <<
535 G4double cosa[4], dist[4]; <<
536 G4int ind[4] = {0}, nside = 0; <<
537 for (G4int i = 0; i < 4; ++i) <<
538 { <<
539 G4double tmp = fNormal[i].dot(v); <<
540 cosa[i] = tmp; <<
541 ind[nside] = (G4int)(tmp > 0) * i; <<
542 nside += (G4int)(tmp > 0); <<
543 dist[i] = fNormal[i].dot(p) - fDist[i]; <<
544 } <<
545 <<
546 // Find intersection (in most of cases nside <<
547 G4double tout = DBL_MAX; <<
548 G4int iside = 0; <<
549 for (G4int i = 0; i < nside; ++i) <<
550 { <<
551 G4int k = ind[i]; <<
552 // Check: leaving the surface <<
553 if (dist[k] >= -halfTolerance) { tout = 0. <<
554 // Compute distance to intersection <<
555 G4double tmp = -dist[k]/cosa[k]; <<
556 if (tmp < tout) { tout = tmp; iside = k; } <<
557 } <<
558 543
559 // Set normal, if required, and return dista << 544 vdotn=vu.dot(fNormal123);
560 if (calcNorm) << 545 if(vdotn > 1e-12) // #we're heading towards this face, so it is a candidate
561 { << 546 {
562 *validNorm = true; << 547 t1=(fCdotN123-p.dot(fNormal123))/vdotn; // # distance to intersection
563 *n = fNormal[iside]; << 548 }
564 } << 549
565 return tout; << 550 vdotn=vu.dot(fNormal134);
>> 551 if(vdotn > 1e-12) // #we're heading towards this face, so it is a candidate
>> 552 {
>> 553 t2=(fCdotN134-p.dot(fNormal134))/vdotn; // # distance to intersection
>> 554 }
>> 555
>> 556 vdotn=vu.dot(fNormal142);
>> 557 if(vdotn > 1e-12) // #we're heading towards this face, so it is a candidate
>> 558 {
>> 559 t3=(fCdotN142-p.dot(fNormal142))/vdotn; // # distance to intersection
>> 560 }
>> 561
>> 562 vdotn=vu.dot(fNormal234);
>> 563 if(vdotn > 1e-12) // #we're heading towards this face, so it is a candidate
>> 564 {
>> 565 t4=(fCdotN234-p.dot(fNormal234))/vdotn; // # distance to intersection
>> 566 }
>> 567
>> 568 tt=std::min(std::min(std::min(t1,t2),t3),t4);
>> 569
>> 570 if (warningFlag && (tt == kInfinity || tt < -fTol))
>> 571 {
>> 572 DumpInfo();
>> 573 std::ostringstream message;
>> 574 message << "No good intersection found or already outside!?" << G4endl
>> 575 << "p = " << p / mm << "mm" << G4endl
>> 576 << "v = " << v << G4endl
>> 577 << "t1, t2, t3, t4 (mm) "
>> 578 << t1/mm << ", " << t2/mm << ", " << t3/mm << ", " << t4/mm;
>> 579 G4Exception("G4Tet::DistanceToOut(p,v,...)", "GeomSolids1002",
>> 580 JustWarning, message);
>> 581 if(validNorm)
>> 582 {
>> 583 *validNorm=false; // flag normal as meaningless
>> 584 }
>> 585 }
>> 586 else if(calcNorm && n != nullptr)
>> 587 {
>> 588 G4ThreeVector normal;
>> 589 if(tt==t1) { normal=fNormal123; }
>> 590 else if (tt==t2) { normal=fNormal134; }
>> 591 else if (tt==t3) { normal=fNormal142; }
>> 592 else if (tt==t4) { normal=fNormal234; }
>> 593 *n=normal;
>> 594 if(validNorm) { *validNorm=true; }
>> 595 }
>> 596
>> 597 return std::max(tt,0.0); // avoid tt<0.0 by a tiny bit
>> 598 // if we are right on a face
566 } 599 }
567 600
568 ////////////////////////////////////////////// << 601 ////////////////////////////////////////////////////////////////////////////
569 // 602 //
570 // Calculate safety distance to surface from i << 603 // Calculate exact shortest distance to any boundary from inside
>> 604 // - If outside return zero
571 // 605 //
572 G4double G4Tet::DistanceToOut(const G4ThreeVec 606 G4double G4Tet::DistanceToOut(const G4ThreeVector& p) const
573 { 607 {
574 G4double dd[4]; << 608 G4double t1,t2,t3,t4;
575 for (G4int i = 0; i < 4; ++i) { dd[i] = fDis << 609 t1=fCdotN123-p.dot(fNormal123); // distance to plane, positive if inside
>> 610 t2=fCdotN134-p.dot(fNormal134); // distance to plane
>> 611 t3=fCdotN142-p.dot(fNormal142); // distance to plane
>> 612 t4=fCdotN234-p.dot(fNormal234); // distance to plane
>> 613
>> 614 // if any one of these is negative, we are outside,
>> 615 // so return zero in that case
576 616
577 G4double dist = std::min(std::min(std::min(d << 617 G4double tmin=std::min(std::min(std::min(t1,t2),t3),t4);
578 return (dist > 0.) ? dist : 0.; << 618 return (tmin < fTol)? 0:tmin;
579 } 619 }
580 620
581 ////////////////////////////////////////////// << 621 //////////////////////////////////////////////////////////////////////////
582 // 622 //
583 // GetEntityType 623 // GetEntityType
584 // 624 //
585 G4GeometryType G4Tet::GetEntityType() const 625 G4GeometryType G4Tet::GetEntityType() const
586 { 626 {
587 return {"G4Tet"}; << 627 return G4String("G4Tet");
588 } <<
589 <<
590 ////////////////////////////////////////////// <<
591 // <<
592 // IsFaceted <<
593 // <<
594 G4bool G4Tet::IsFaceted() const <<
595 { <<
596 return true; <<
597 } 628 }
598 629
599 ////////////////////////////////////////////// << 630 //////////////////////////////////////////////////////////////////////////
600 // 631 //
601 // Make a clone of the object 632 // Make a clone of the object
602 // 633 //
603 G4VSolid* G4Tet::Clone() const 634 G4VSolid* G4Tet::Clone() const
604 { 635 {
605 return new G4Tet(*this); 636 return new G4Tet(*this);
606 } 637 }
607 638
608 ////////////////////////////////////////////// << 639 //////////////////////////////////////////////////////////////////////////
609 // 640 //
610 // Stream object contents to an output stream 641 // Stream object contents to an output stream
611 // 642 //
612 std::ostream& G4Tet::StreamInfo(std::ostream& 643 std::ostream& G4Tet::StreamInfo(std::ostream& os) const
613 { 644 {
614 G4long oldprc = os.precision(16); << 645 G4int oldprc = os.precision(16);
615 os << "------------------------------------- 646 os << "-----------------------------------------------------------\n"
616 << " *** Dump for solid - " << GetName << 647 << " *** Dump for solid - " << GetName() << " ***\n"
617 << " ================================= << 648 << " ===================================================\n"
618 << " Solid type: " << GetEntityType() << << 649 << " Solid type: G4Tet\n"
619 << " Parameters: \n" << 650 << " Parameters: \n"
620 << " anchor: " << fVertex[0]/mm << " m << 651 << " anchor: " << fAnchor/mm << " mm \n"
621 << " p1 : " << fVertex[1]/mm << " m << 652 << " p2: " << fP2/mm << " mm \n"
622 << " p2 : " << fVertex[2]/mm << " m << 653 << " p3: " << fP3/mm << " mm \n"
623 << " p3 : " << fVertex[3]/mm << " m << 654 << " p4: " << fP4/mm << " mm \n"
624 << "------------------------------------- << 655 << " normal123: " << fNormal123 << " \n"
>> 656 << " normal134: " << fNormal134 << " \n"
>> 657 << " normal142: " << fNormal142 << " \n"
>> 658 << " normal234: " << fNormal234 << " \n"
>> 659 << "-----------------------------------------------------------\n";
625 os.precision(oldprc); 660 os.precision(oldprc);
>> 661
626 return os; 662 return os;
627 } 663 }
628 664
629 ////////////////////////////////////////////// 665 ////////////////////////////////////////////////////////////////////////
630 // 666 //
631 // Return random point on the surface << 667 // GetPointOnFace
>> 668 //
>> 669 // Auxiliary method for get point on surface
>> 670 //
>> 671 G4ThreeVector G4Tet::GetPointOnFace(G4ThreeVector p1, G4ThreeVector p2,
>> 672 G4ThreeVector p3, G4double& area) const
>> 673 {
>> 674 G4double lambda1,lambda2;
>> 675 G4ThreeVector v, w;
>> 676
>> 677 v = p3 - p1;
>> 678 w = p1 - p2;
>> 679
>> 680 lambda1 = G4RandFlat::shoot(0.,1.);
>> 681 lambda2 = G4RandFlat::shoot(0.,lambda1);
>> 682
>> 683 area = 0.5*(v.cross(w)).mag();
>> 684
>> 685 return (p2 + lambda1*w + lambda2*v);
>> 686 }
>> 687
>> 688 ////////////////////////////////////////////////////////////////////////////
>> 689 //
>> 690 // GetPointOnSurface
632 // 691 //
633 G4ThreeVector G4Tet::GetPointOnSurface() const 692 G4ThreeVector G4Tet::GetPointOnSurface() const
634 { 693 {
635 constexpr G4int iface[4][3] = { {0,1,2}, {0, << 694 G4double chose,aOne,aTwo,aThree,aFour;
>> 695 G4ThreeVector p1, p2, p3, p4;
>> 696
>> 697 p1 = GetPointOnFace(fAnchor,fP2,fP3,aOne);
>> 698 p2 = GetPointOnFace(fAnchor,fP4,fP3,aTwo);
>> 699 p3 = GetPointOnFace(fAnchor,fP4,fP2,aThree);
>> 700 p4 = GetPointOnFace(fP4,fP3,fP2,aFour);
>> 701
>> 702 chose = G4RandFlat::shoot(0.,aOne+aTwo+aThree+aFour);
>> 703 if( (chose>=0.) && (chose <aOne) ) {return p1;}
>> 704 else if( (chose>=aOne) && (chose < aOne+aTwo) ) {return p2;}
>> 705 else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThree) ) {return p3;}
>> 706 return p4;
>> 707 }
>> 708
>> 709 ////////////////////////////////////////////////////////////////////////
>> 710 //
>> 711 // GetVertices
>> 712 //
>> 713 std::vector<G4ThreeVector> G4Tet::GetVertices() const
>> 714 {
>> 715 std::vector<G4ThreeVector> vertices(4);
>> 716 vertices[0] = fAnchor;
>> 717 vertices[1] = fP2;
>> 718 vertices[2] = fP3;
>> 719 vertices[3] = fP4;
636 720
637 // Select face << 721 return vertices;
638 G4double select = fSurfaceArea*G4QuickRand() <<
639 G4int i = 0; <<
640 i += (G4int)(select > fArea[0]); <<
641 i += (G4int)(select > fArea[0] + fArea[1]); <<
642 i += (G4int)(select > fArea[0] + fArea[1] + <<
643 <<
644 // Set selected triangle <<
645 G4ThreeVector p0 = fVertex[iface[i][0]]; <<
646 G4ThreeVector e1 = fVertex[iface[i][1]] - p0 <<
647 G4ThreeVector e2 = fVertex[iface[i][2]] - p0 <<
648 <<
649 // Return random point <<
650 G4double r1 = G4QuickRand(); <<
651 G4double r2 = G4QuickRand(); <<
652 return (r1 + r2 > 1.) ? <<
653 p0 + e1*(1. - r1) + e2*(1. - r2) : p0 + e1 <<
654 } 722 }
655 723
656 ////////////////////////////////////////////// 724 ////////////////////////////////////////////////////////////////////////
657 // 725 //
658 // Return volume of the tetrahedron << 726 // GetCubicVolume
659 // 727 //
660 G4double G4Tet::GetCubicVolume() 728 G4double G4Tet::GetCubicVolume()
661 { 729 {
662 return fCubicVolume; 730 return fCubicVolume;
663 } 731 }
664 732
665 ////////////////////////////////////////////// 733 ////////////////////////////////////////////////////////////////////////
666 // 734 //
667 // Return surface area of the tetrahedron << 735 // GetSurfaceArea
668 // 736 //
669 G4double G4Tet::GetSurfaceArea() 737 G4double G4Tet::GetSurfaceArea()
670 { 738 {
671 return fSurfaceArea; 739 return fSurfaceArea;
672 } 740 }
673 741
>> 742 // Methods for visualisation
>> 743
674 ////////////////////////////////////////////// 744 ////////////////////////////////////////////////////////////////////////
675 // 745 //
676 // Methods for visualisation << 746 // DescribeYourselfTo
677 // 747 //
678 void G4Tet::DescribeYourselfTo (G4VGraphicsSce << 748 void G4Tet::DescribeYourselfTo (G4VGraphicsScene& scene) const
679 { 749 {
680 scene.AddSolid (*this); 750 scene.AddSolid (*this);
681 } 751 }
682 752
683 ////////////////////////////////////////////// 753 ////////////////////////////////////////////////////////////////////////
684 // 754 //
685 // Return VisExtent << 755 // GetExtent
686 // 756 //
687 G4VisExtent G4Tet::GetExtent() const << 757 G4VisExtent G4Tet::GetExtent() const
688 { 758 {
689 return { fBmin.x(), fBmax.x(), << 759 return G4VisExtent (fXMin, fXMax, fYMin, fYMax, fZMin, fZMax);
690 fBmin.y(), fBmax.y(), <<
691 fBmin.z(), fBmax.z() }; <<
692 } 760 }
693 761
694 ////////////////////////////////////////////// 762 ////////////////////////////////////////////////////////////////////////
695 // 763 //
696 // CreatePolyhedron 764 // CreatePolyhedron
697 // 765 //
698 G4Polyhedron* G4Tet::CreatePolyhedron() const << 766 G4Polyhedron* G4Tet::CreatePolyhedron() const
699 { 767 {
700 // Check orientation of vertices << 768 G4Polyhedron* ph = new G4Polyhedron;
701 G4ThreeVector v1 = fVertex[1] - fVertex[0]; <<
702 G4ThreeVector v2 = fVertex[2] - fVertex[0]; <<
703 G4ThreeVector v3 = fVertex[3] - fVertex[0]; <<
704 G4bool invert = v1.cross(v2).dot(v3) < 0.; <<
705 G4int k2 = (invert) ? 3 : 2; <<
706 G4int k3 = (invert) ? 2 : 3; <<
707 <<
708 // Set coordinates of vertices <<
709 G4double xyz[4][3]; 769 G4double xyz[4][3];
710 for (G4int i = 0; i < 3; ++i) << 770 const G4int faces[4][4]={{1,3,2,0},{1,4,3,0},{1,2,4,0},{2,3,4,0}};
711 { << 771 xyz[0][0]=fAnchor.x(); xyz[0][1]=fAnchor.y(); xyz[0][2]=fAnchor.z();
712 xyz[0][i] = fVertex[0][i]; << 772 xyz[1][0]=fP2.x(); xyz[1][1]=fP2.y(); xyz[1][2]=fP2.z();
713 xyz[1][i] = fVertex[1][i]; << 773 xyz[2][0]=fP3.x(); xyz[2][1]=fP3.y(); xyz[2][2]=fP3.z();
714 xyz[2][i] = fVertex[k2][i]; << 774 xyz[3][0]=fP4.x(); xyz[3][1]=fP4.y(); xyz[3][2]=fP4.z();
715 xyz[3][i] = fVertex[k3][i]; <<
716 } <<
717 775
718 // Create polyhedron <<
719 G4int faces[4][4] = { {1,3,2,0}, {1,4,3,0}, <<
720 auto ph = new G4Polyhedron; <<
721 ph->createPolyhedron(4,4,xyz,faces); 776 ph->createPolyhedron(4,4,xyz,faces);
722 777
723 return ph; 778 return ph;
724 } 779 }
725 780
726 ////////////////////////////////////////////// 781 ////////////////////////////////////////////////////////////////////////
727 // 782 //
728 // GetPolyhedron 783 // GetPolyhedron
729 // 784 //
730 G4Polyhedron* G4Tet::GetPolyhedron() const 785 G4Polyhedron* G4Tet::GetPolyhedron() const
731 { 786 {
732 if (fpPolyhedron == nullptr || 787 if (fpPolyhedron == nullptr ||
733 fRebuildPolyhedron || 788 fRebuildPolyhedron ||
734 fpPolyhedron->GetNumberOfRotationStepsAt 789 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
735 fpPolyhedron->GetNumberOfRotationSteps() 790 fpPolyhedron->GetNumberOfRotationSteps())
736 { << 791 {
737 G4AutoLock l(&polyhedronMutex); << 792 G4AutoLock l(&polyhedronMutex);
738 delete fpPolyhedron; << 793 delete fpPolyhedron;
739 fpPolyhedron = CreatePolyhedron(); << 794 fpPolyhedron = CreatePolyhedron();
740 fRebuildPolyhedron = false; << 795 fRebuildPolyhedron = false;
741 l.unlock(); << 796 l.unlock();
742 } << 797 }
743 return fpPolyhedron; 798 return fpPolyhedron;
744 } 799 }
745 800
746 #endif 801 #endif
747 802