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