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