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