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