| 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 result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio << 19 // * technical work of the GEANT4 collaboration and of QinetiQ Ltd, *
>> 20 // * subject DEFCON 705 IPR conditions. *
20 // * By using, copying, modifying or distri 21 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 22 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 23 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 24 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 25 // ********************************************************************
25 // 26 //
26 // G4TriangularFacet implementation << 27 // $Id: G4TriangularFacet.cc,v 1.7 2007/02/15 17:03:49 gcosmo Exp $
>> 28 // GEANT4 tag $Name: geant4-08-03 $
27 // 29 //
28 // 31.10.2004, P R Truscott, QinetiQ Ltd, UK - << 30 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
29 // 01.08.2007, P R Truscott, QinetiQ Ltd, UK << 31 //
30 // Significant modification t << 32 // MODULE: G4TriangularFacet.cc
31 // based on patches/observati << 33 //
32 // Holmberg. << 34 // Date: 15/06/2005
33 // 12.10.2012, M Gayer, CERN << 35 // Author: P R Truscott
34 // New implementation reducin << 36 // Organisation: QinetiQ Ltd, UK
35 // and considerable CPU speed << 37 // Customer: UK Ministry of Defence : RAO CRP TD Electronic Systems
36 // implementation of G4Tessel << 38 // Contract: C/MAT/N03517
37 // 23.02.2016, E Tcherniaev, CERN << 39 //
38 // Improved test to detect de << 40 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
39 // too narrow) triangles. << 41 //
40 // ------------------------------------------- << 42 // CHANGE HISTORY
>> 43 // --------------
>> 44 //
>> 45 // 31 October 2004, P R Truscott, QinetiQ Ltd, UK - Created.
>> 46 //
>> 47 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
41 48
42 #include "G4TriangularFacet.hh" 49 #include "G4TriangularFacet.hh"
43 << 50 #include "globals.hh"
44 #include "Randomize.hh" 51 #include "Randomize.hh"
45 #include "G4TessellatedGeometryAlgorithms.hh" <<
46 <<
47 using namespace std; <<
48 52
49 ////////////////////////////////////////////// 53 ///////////////////////////////////////////////////////////////////////////////
50 // 54 //
51 // Definition of triangular facet using absolu << 55 // Definition of triangular facet using absolute vectors to vertices.
52 // From this for first vector is retained to d 56 // From this for first vector is retained to define the facet location and
53 // two relative vectors (E0 and E1) define the 57 // two relative vectors (E0 and E1) define the sides and orientation of
54 // the outward surface normal. 58 // the outward surface normal.
55 // 59 //
56 G4TriangularFacet::G4TriangularFacet (const G4 << 60 G4TriangularFacet::G4TriangularFacet (const G4ThreeVector Pt0,
57 const G4 << 61 const G4ThreeVector vt1, const G4ThreeVector vt2,
58 const G4 << 62 G4FacetVertexType vertexType)
59 G4 << 63 : G4VFacet()
60 { 64 {
61 fVertices = new vector<G4ThreeVector>(3); << 65 P0 = Pt0;
62 << 66 nVertices = 3;
63 SetVertex(0, vt0); <<
64 if (vertexType == ABSOLUTE) 67 if (vertexType == ABSOLUTE)
65 { 68 {
66 SetVertex(1, vt1); << 69 P.push_back(vt1);
67 SetVertex(2, vt2); << 70 P.push_back(vt2);
68 fE1 = vt1 - vt0; << 71
69 fE2 = vt2 - vt0; << 72 E.push_back(vt1 - P0);
>> 73 E.push_back(vt2 - P0);
70 } 74 }
71 else 75 else
72 { 76 {
73 SetVertex(1, vt0 + vt1); << 77 P.push_back(P0 + vt1);
74 SetVertex(2, vt0 + vt2); << 78 P.push_back(P0 + vt2);
75 fE1 = vt1; << 79
76 fE2 = vt2; << 80 E.push_back(vt1);
77 } << 81 E.push_back(vt2);
78 << 82 }
79 G4ThreeVector E1xE2 = fE1.cross(fE2); << 83
80 fArea = 0.5 * E1xE2.mag(); << 84 G4double Emag1 = E[0].mag();
81 for (G4int i = 0; i < 3; ++i) fIndices[i] = << 85 G4double Emag2 = E[1].mag();
82 << 86 G4double Emag3 = (E[1]-E[0]).mag();
83 fIsDefined = true; << 87
84 G4double delta = kCarTolerance; // Set toler << 88 if (Emag1 <= kCarTolerance || Emag2 <= kCarTolerance ||
85 << 89 Emag3 <= kCarTolerance)
86 // Check length of edges << 90 {
87 // << 91 G4Exception("G4TriangularFacet::G4TriangularFacet()", "InvalidSetup",
88 G4double leng1 = fE1.mag(); << 92 JustWarning, "Length of sides of facet are too small.");
89 G4double leng2 = (fE2-fE1).mag(); << 93 G4cerr << G4endl;
90 G4double leng3 = fE2.mag(); << 94 G4cerr << "P0 = " << P0 << G4endl;
91 if (leng1 <= delta || leng2 <= delta || leng << 95 G4cerr << "P1 = " << P[0] << G4endl;
92 { << 96 G4cerr << "P2 = " << P[1] << G4endl;
93 fIsDefined = false; << 97 G4cerr << "Side lengths = P0->P1" << Emag1 << G4endl;
94 } << 98 G4cerr << "Side lengths = P0->P2" << Emag2 << G4endl;
95 << 99 G4cerr << "Side lengths = P1->P2" << Emag3 << G4endl;
96 // Check min height of triangle << 100 G4cerr << G4endl;
97 // << 101
98 if (fIsDefined) << 102 isDefined = false;
99 { << 103 geometryType = "G4TriangularFacet";
100 if (2.*fArea/std::max(std::max(leng1,leng2 << 104 surfaceNormal = G4ThreeVector(0.0,0.0,0.0);
101 { << 105 a = 0.0;
102 fIsDefined = false; << 106 b = 0.0;
103 } << 107 c = 0.0;
104 } << 108 det = 0.0;
105 <<
106 // Define facet <<
107 // <<
108 if (!fIsDefined) <<
109 { <<
110 ostringstream message; <<
111 message << "Facet is too small or too narr <<
112 << "Triangle area = " << fArea << <<
113 << "P0 = " << GetVertex(0) << G4en <<
114 << "P1 = " << GetVertex(1) << G4en <<
115 << "P2 = " << GetVertex(2) << G4en <<
116 << "Side1 length (P0->P1) = " << l <<
117 << "Side2 length (P1->P2) = " << l <<
118 << "Side3 length (P2->P0) = " << l <<
119 G4Exception("G4TriangularFacet::G4Triangul <<
120 "GeomSolids1001", JustWarning, message); <<
121 fSurfaceNormal.set(0,0,0); <<
122 fA = fB = fC = 0.0; <<
123 fDet = 0.0; <<
124 fCircumcentre = vt0 + 0.5*fE1 + 0.5*fE2; <<
125 fArea = fRadius = 0.0; <<
126 } 109 }
127 else 110 else
128 { <<
129 fSurfaceNormal = E1xE2.unit(); <<
130 fA = fE1.mag2(); <<
131 fB = fE1.dot(fE2); <<
132 fC = fE2.mag2(); <<
133 fDet = std::fabs(fA*fC - fB*fB); <<
134 <<
135 fCircumcentre = <<
136 vt0 + (E1xE2.cross(fE1)*fC + fE2.cross(E <<
137 fRadius = (fCircumcentre - vt0).mag(); <<
138 } <<
139 } <<
140 <<
141 ////////////////////////////////////////////// <<
142 // <<
143 G4TriangularFacet::G4TriangularFacet () <<
144 { <<
145 fVertices = new vector<G4ThreeVector>(3); <<
146 G4ThreeVector zero(0,0,0); <<
147 SetVertex(0, zero); <<
148 SetVertex(1, zero); <<
149 SetVertex(2, zero); <<
150 for (G4int i = 0; i < 3; ++i) fIndices[i] = <<
151 fIsDefined = false; <<
152 fSurfaceNormal.set(0,0,0); <<
153 fA = fB = fC = 0; <<
154 fE1 = zero; <<
155 fE2 = zero; <<
156 fDet = 0.0; <<
157 fArea = fRadius = 0.0; <<
158 } <<
159 <<
160 ////////////////////////////////////////////// <<
161 // <<
162 G4TriangularFacet::~G4TriangularFacet () <<
163 { <<
164 SetVertices(nullptr); <<
165 } <<
166 <<
167 ////////////////////////////////////////////// <<
168 // <<
169 void G4TriangularFacet::CopyFrom (const G4Tria <<
170 { <<
171 auto p = (char *) &rhs; <<
172 copy(p, p + sizeof(*this), (char *)this); <<
173 <<
174 if (fIndices[0] < 0 && fVertices == nullptr) <<
175 { <<
176 fVertices = new vector<G4ThreeVector>(3); <<
177 for (G4int i = 0; i < 3; ++i) (*fVertices) <<
178 } <<
179 } <<
180 <<
181 ////////////////////////////////////////////// <<
182 // <<
183 void G4TriangularFacet::MoveFrom (G4Triangular <<
184 { <<
185 fSurfaceNormal = std::move(rhs.fSurfaceNorma <<
186 fArea = rhs.fArea; <<
187 fCircumcentre = std::move(rhs.fCircumcentre) <<
188 fRadius = rhs.fRadius; <<
189 fIndices = rhs.fIndices; <<
190 fA = rhs.fA; fB = rhs.fB; fC = rhs.fC; <<
191 fDet = rhs.fDet; <<
192 fSqrDist = rhs.fSqrDist; <<
193 fE1 = std::move(rhs.fE1); fE2 = std::move(rh <<
194 fIsDefined = rhs.fIsDefined; <<
195 fVertices = rhs.fVertices; <<
196 rhs.fVertices = nullptr; <<
197 } <<
198 <<
199 ////////////////////////////////////////////// <<
200 // <<
201 G4TriangularFacet::G4TriangularFacet (const G4 <<
202 : G4VFacet(rhs) <<
203 { <<
204 CopyFrom(rhs); <<
205 } <<
206 <<
207 ////////////////////////////////////////////// <<
208 // <<
209 G4TriangularFacet::G4TriangularFacet (G4Triang <<
210 : G4VFacet(rhs) <<
211 { <<
212 MoveFrom(rhs); <<
213 } <<
214 <<
215 ////////////////////////////////////////////// <<
216 // <<
217 G4TriangularFacet& <<
218 G4TriangularFacet::operator=(const G4Triangula <<
219 { <<
220 SetVertices(nullptr); <<
221 <<
222 if (this != &rhs) <<
223 { 111 {
224 delete fVertices; << 112 isDefined = true;
225 CopyFrom(rhs); << 113 geometryType = "G4TriangularFacet";
>> 114 surfaceNormal = E[0].cross(E[1]).unit();
>> 115 a = E[0].mag2();
>> 116 b = E[0].dot(E[1]);
>> 117 c = E[1].mag2();
>> 118 det = std::abs(a*c - b*b);
>> 119
>> 120 sMin = -0.5*kCarTolerance/std::sqrt(a);
>> 121 sMax = 1.0 - sMin;
>> 122 tMin = -0.5*kCarTolerance/std::sqrt(c);
>> 123
>> 124 G4ThreeVector vtmp = 0.25 * (E[0] + E[1]);
>> 125 centroid = P0 + vtmp;
>> 126 radiusSqr = vtmp.mag2();
>> 127 radius = std::sqrt(radiusSqr);
>> 128
>> 129 for (size_t i=0; i<3; i++) I.push_back(0);
226 } 130 }
227 <<
228 return *this; <<
229 } 131 }
230 132
231 ////////////////////////////////////////////// 133 ///////////////////////////////////////////////////////////////////////////////
232 // 134 //
233 G4TriangularFacet& << 135 G4TriangularFacet::~G4TriangularFacet ()
234 G4TriangularFacet::operator=(G4TriangularFacet <<
235 { 136 {
236 SetVertices(nullptr); << 137 P.clear();
237 << 138 E.clear();
238 if (this != &rhs) << 139 I.clear();
239 { <<
240 delete fVertices; <<
241 MoveFrom(rhs); <<
242 } <<
243 <<
244 return *this; <<
245 } 140 }
246 141
247 ////////////////////////////////////////////// 142 ///////////////////////////////////////////////////////////////////////////////
248 // 143 //
249 // GetClone << 144 G4VFacet *G4TriangularFacet::GetClone ()
250 // <<
251 // Simple member function to generate fA dupli <<
252 // <<
253 G4VFacet* G4TriangularFacet::GetClone () <<
254 { 145 {
255 auto fc = new G4TriangularFacet (GetVertex(0 << 146 G4TriangularFacet *fc = new G4TriangularFacet (P0, P[0], P[1], ABSOLUTE);
256 GetVertex(2 << 147 G4VFacet *cc = 0;
257 return fc; << 148 cc = fc;
>> 149 return cc;
258 } 150 }
259 151
260 ////////////////////////////////////////////// 152 ///////////////////////////////////////////////////////////////////////////////
261 // 153 //
262 // GetFlippedFacet << 154 G4TriangularFacet *G4TriangularFacet::GetFlippedFacet ()
263 // <<
264 // Member function to generate an identical fa <<
265 // pointing at 180 degrees. <<
266 // <<
267 G4TriangularFacet* G4TriangularFacet::GetFlipp <<
268 { 155 {
269 auto flipped = new G4TriangularFacet (GetVer << 156 G4TriangularFacet *flipped = new G4TriangularFacet (P0, P[1], P[0], ABSOLUTE);
270 GetVer <<
271 return flipped; 157 return flipped;
272 } 158 }
273 159
274 ////////////////////////////////////////////// 160 ///////////////////////////////////////////////////////////////////////////////
275 // 161 //
276 // Distance (G4ThreeVector) << 162 // Determine the closest distance from the facet to the point p. If the
277 // << 163 // direction of the vector to the closest point is outward-going and outgoing
278 // Determines the vector between p and the clo << 164 // is true or the vector is in-going and outgoing is false then the distance
279 // This is based on the algorithm published in << 165 // is returned. Otherwise kInfinity is returned.
280 // Graphics," Philip J Scheider and David H Eb << 166 //
281 // 2003. at the time of writing, the algorith << 167 G4ThreeVector G4TriangularFacet::Distance (const G4ThreeVector &p)
282 // technical note "Distance between point and << 168 {
283 // at http://www.geometrictools.com/Documentat << 169 G4ThreeVector D = P0 - p;
284 // << 170 G4double d = E[0].dot(D);
285 // The by-product is the square-distance fSqrD << 171 G4double e = E[1].dot(D);
286 // in case needed by the other "Distance" memb << 172 G4double f = D.mag2();
287 // << 173 G4double s = b*e - c*d;
288 G4ThreeVector G4TriangularFacet::Distance (con << 174 G4double t = b*d - a*e;
289 { << 175 G4double sqrDist = 0.0;
290 G4ThreeVector D = GetVertex(0) - p; << 176
291 G4double d = fE1.dot(D); << 177 if (s+t <= det)
292 G4double e = fE2.dot(D); <<
293 G4double f = D.mag2(); <<
294 G4double q = fB*e - fC*d; <<
295 G4double t = fB*d - fA*e; <<
296 fSqrDist = 0.; <<
297 <<
298 if (q+t <= fDet) <<
299 { 178 {
300 if (q < 0.0) << 179 if (s < 0.0)
301 { 180 {
302 if (t < 0.0) 181 if (t < 0.0)
303 { 182 {
304 // << 183 //
305 // We are in region 4. << 184 // We are in region 4.
306 // << 185 //
307 if (d < 0.0) 186 if (d < 0.0)
308 { 187 {
309 t = 0.0; 188 t = 0.0;
310 if (-d >= fA) {q = 1.0; fSqrDist = f << 189 if (-d >= a) {s = 1.0; sqrDist = a + 2.0*d + f;}
311 else {q = -d/fA; fSqrDist = << 190 else {s = -d/a; sqrDist = d*s + f;}
312 } 191 }
313 else 192 else
314 { 193 {
315 q = 0.0; << 194 s = 0.0;
316 if (e >= 0.0) {t = 0.0; fSqrDi << 195 if (e >= 0.0) {t = 0.0; sqrDist = f;}
317 else if (-e >= fC) {t = 1.0; fSqrD << 196 else if (-e >= c) {t = 1.0; sqrDist = c + 2.0*e + f;}
318 else {t = -e/fC; fSqr << 197 else {t = -e/c; sqrDist = e*t + f;}
319 } 198 }
320 } 199 }
321 else 200 else
322 { 201 {
323 // << 202 //
324 // We are in region 3. << 203 // We are in region 3.
325 // << 204 //
326 q = 0.0; << 205 s = 0.0;
327 if (e >= 0.0) {t = 0.0; fSqrDist << 206 if (e >= 0.0) {t = 0.0; sqrDist = f;}
328 else if (-e >= fC) {t = 1.0; fSqrDist << 207 else if (-e >= c) {t = 0.0; sqrDist = c + 2.0*e + f;}
329 else {t = -e/fC; fSqrDis << 208 else {t = -e/c; sqrDist = e*t + f;}
330 } 209 }
331 } 210 }
332 else if (t < 0.0) 211 else if (t < 0.0)
333 { 212 {
334 // << 213 //
335 // We are in region 5. << 214 // We are in region 5.
336 // << 215 //
337 t = 0.0; 216 t = 0.0;
338 if (d >= 0.0) {q = 0.0; fSqrDist = << 217 if (d >= 0.0) {s = 0.0; sqrDist = f;}
339 else if (-d >= fA) {q = 1.0; fSqrDist = << 218 else if (-d >= a) {s = 1.0; sqrDist = a + 2.0*d + f;}
340 else {q = -d/fA; fSqrDist << 219 else {s = -d/a; sqrDist = d*s + f;}
341 } 220 }
342 else 221 else
343 { 222 {
344 // << 223 //
345 // We are in region 0. << 224 // We are in region 0.
346 // << 225 //
347 G4double dist = fSurfaceNormal.dot(D); << 226 G4double invDet = 1.0 / det;
348 fSqrDist = dist*dist; << 227 s *= invDet;
349 return fSurfaceNormal*dist; << 228 t *= invDet;
>> 229 sqrDist = s*(a*s + b*t + 2.0*d) + t*(b*s + c*t + 2.0*e) + f;
350 } 230 }
351 } 231 }
352 else 232 else
353 { 233 {
354 if (q < 0.0) << 234 G4double tmp0 = 0.0;
355 { << 235 G4double tmp1 = 0.0;
356 // << 236 G4double numer = 0.0;
357 // We are in region 2. << 237 G4double denom = 0.0;
358 // << 238 if (s < 0.0)
359 G4double tmp0 = fB + d; << 239 {
360 G4double tmp1 = fC + e; << 240 //
>> 241 // We are in region 2.
>> 242 //
>> 243 tmp0 = b + d;
>> 244 tmp1 = c + e;
361 if (tmp1 > tmp0) 245 if (tmp1 > tmp0)
362 { 246 {
363 G4double numer = tmp1 - tmp0; << 247 numer = tmp1 - tmp0;
364 G4double denom = fA - 2.0*fB + fC; << 248 denom = a - 2.0*b*c;
365 if (numer >= denom) {q = 1.0; t = 0.0; << 249 if (numer >= denom) {s = 1.0; t = 0.0; sqrDist = a + 2.0*d + f;}
366 else 250 else
367 { 251 {
368 q = numer/denom; << 252 s = numer/denom;
369 t = 1.0 - q; << 253 t = 1.0 - s;
370 fSqrDist = q*(fA*q + fB*t +2.0*d) + << 254 sqrDist = s*(a*s + b*t +2.0*d) + t*(b*s + c*t + 2.0*e) + f;
371 } 255 }
372 } 256 }
373 else 257 else
374 { 258 {
375 q = 0.0; << 259 s = 0.0;
376 if (tmp1 <= 0.0) {t = 1.0; fSqrDi << 260 if (tmp1 <= 0.0) {t = 1.0; sqrDist = c + 2.0*e + f;}
377 else if (e >= 0.0) {t = 0.0; fSqrDi << 261 else if (e >= 0.0) {t = 0.0; sqrDist = f;}
378 else {t = -e/fC; fSqr << 262 else {t = -e/c; sqrDist = e*t + f;}
379 } 263 }
380 } 264 }
381 else if (t < 0.0) 265 else if (t < 0.0)
382 { 266 {
383 // << 267 //
384 // We are in region 6. << 268 // We are in region 6.
385 // << 269 //
386 G4double tmp0 = fB + e; << 270 tmp0 = b + e;
387 G4double tmp1 = fA + d; << 271 tmp1 = a + d;
388 if (tmp1 > tmp0) 272 if (tmp1 > tmp0)
389 { 273 {
390 G4double numer = tmp1 - tmp0; << 274 numer = tmp1 - tmp0;
391 G4double denom = fA - 2.0*fB + fC; << 275 denom = a - 2.0*b*c;
392 if (numer >= denom) {t = 1.0; q = 0.0; << 276 if (numer >= denom) {t = 1.0; s = 0.0; sqrDist = c + 2.0*e + f;}
393 else 277 else
394 { 278 {
395 t = numer/denom; 279 t = numer/denom;
396 q = 1.0 - t; << 280 s = 1.0 - t;
397 fSqrDist = q*(fA*q + fB*t +2.0*d) + << 281 sqrDist = s*(a*s + b*t +2.0*d) + t*(b*s + c*t + 2.0*e) + f;
398 } 282 }
399 } 283 }
400 else 284 else
401 { 285 {
402 t = 0.0; 286 t = 0.0;
403 if (tmp1 <= 0.0) {q = 1.0; fSqrDi << 287 if (tmp1 <= 0.0) {s = 1.0; sqrDist = a + 2.0*d + f;}
404 else if (d >= 0.0) {q = 0.0; fSqrDi << 288 else if (d >= 0.0) {s = 0.0; sqrDist = f;}
405 else {q = -d/fA; fSqr << 289 else {s = -d/a; sqrDist = d*s + f;}
406 } 290 }
407 } 291 }
408 else 292 else
409 // << 293 //
410 // We are in region 1. << 294 // We are in region 1.
411 // << 295 //
412 { 296 {
413 G4double numer = fC + e - fB - d; << 297 numer = c + f - b - d;
414 if (numer <= 0.0) 298 if (numer <= 0.0)
415 { 299 {
416 q = 0.0; << 300 s = 0.0;
417 t = 1.0; 301 t = 1.0;
418 fSqrDist = fC + 2.0*e + f; << 302 sqrDist = c + 2.0*e*f;
419 } 303 }
420 else 304 else
421 { 305 {
422 G4double denom = fA - 2.0*fB + fC; << 306 denom = a - 2.0*b*c;
423 if (numer >= denom) {q = 1.0; t = 0.0; << 307 if (numer >= denom) {s = 1.0; t = 0.0; sqrDist = a + 2.0*d + f;}
424 else 308 else
425 { 309 {
426 q = numer/denom; << 310 s = numer/denom;
427 t = 1.0 - q; << 311 t = 1.0 - s;
428 fSqrDist = q*(fA*q + fB*t + 2.0*d) + << 312 sqrDist = s*(a*s + b*t + 2.0*d) + t*(b*s + c*t + 2.0*e) + f;
429 } 313 }
430 } 314 }
431 } 315 }
432 } << 316 }
433 // << 317 return D + s*E[0] + t*E[1];
434 // <<
435 // Do fA check for rounding errors in the di <<
436 // the conventional methods for calculating <<
437 // near to or at the surface (as required by <<
438 // We'll therefore also use the magnitude-sq <<
439 // (Note that I've also tried to get around <<
440 // existing equations for <<
441 // <<
442 // fSqrDist = function(fA,fB,fC,d,q,t) <<
443 // <<
444 // and use fA more accurate addition process <<
445 // breakdown of cummutitivity [where (A+B)+C <<
446 // doesn't work. <<
447 // Calculation from u = D + q*fE1 + t*fE2 is <<
448 // more robust. <<
449 // <<
450 if (fSqrDist < 0.0) fSqrDist = 0.; <<
451 G4ThreeVector u = D + q*fE1 + t*fE2; <<
452 G4double u2 = u.mag2(); <<
453 // <<
454 // The following (part of the roundoff error <<
455 // updates. <<
456 // <<
457 if (fSqrDist > u2) fSqrDist = u2; <<
458 <<
459 return u; <<
460 } 318 }
461 319
462 ////////////////////////////////////////////// 320 ///////////////////////////////////////////////////////////////////////////////
463 // 321 //
464 // Distance (G4ThreeVector, G4double) << 322 G4double G4TriangularFacet::Distance (const G4ThreeVector &p,
465 // << 323 const G4double minDist)
466 // Determines the closest distance between poi <<
467 // use of G4ThreeVector G4TriangularFacet::Dis <<
468 // square of the distance in variable fSqrDist <<
469 // the distance is to be greater than minDist, <<
470 // computation and return fA very large number <<
471 // <<
472 G4double G4TriangularFacet::Distance (const G4 <<
473 G4 <<
474 { 324 {
475 // << 325 /*G4ThreeVector D = P0 - p;
476 // Start with quicky test to determine if th << 326 G4double d = E[0].dot(D);
477 // the triangle is any closer to p than minD << 327 G4double e = E[1].dot(D);
478 // about more accurate test. << 328 G4double s = b*e - c*d;
479 // << 329 G4double t = b*d - a*e;*/
480 G4double dist = kInfinity; 330 G4double dist = kInfinity;
481 if ((p-fCircumcentre).mag()-fRadius < minDis << 331
>> 332 /*if (s+t > 1.0 || s < 0.0 || t < 0.0)
482 { 333 {
483 // << 334 G4ThreeVector D0 = P0 - p;
484 // It's possible that the triangle is clos << 335 G4ThreeVector D1 = P[0] - p;
485 // so do more accurate assessment. << 336 G4ThreeVector D2 = P[1] - p;
486 // << 337
487 dist = Distance(p).mag(); << 338 G4double d0 = D0.mag();
488 } << 339 G4double d1 = D1.mag();
>> 340 G4double d2 = D2.mag();
>> 341
>> 342 dist = min(d0, min(d1, d2));
>> 343 if (dist > minDist) return kInfinity;
>> 344 }*/
>> 345
>> 346 dist = Distance(p).mag();
>> 347 if (dist > minDist) return kInfinity;
>> 348
489 return dist; 349 return dist;
490 } 350 }
491 351
492 ////////////////////////////////////////////// 352 ///////////////////////////////////////////////////////////////////////////////
493 // 353 //
494 // Distance (G4ThreeVector, G4double, G4bool) << 354 // Determine the distance to point p bearing in mind that if the distance is
495 // << 355 // likely to be longer than minDist, forget doing further calculation and
496 // Determine the distance to point p. kInfini << 356 // return kInfinity.
497 // (1) outgoing is TRUE and the dot product of << 357 //
498 // and the displacement vector from p to t << 358 G4double G4TriangularFacet::Distance (const G4ThreeVector &p,
499 // (2) outgoing is FALSE and the dot product o << 359 const G4double, const G4bool outgoing)
500 // and the displacement vector from p to t << 360 {
501 // If approximate methods show the distance is << 361 /*G4ThreeVector D = P0 - p;
502 // forget about further computation and return << 362 G4double d = E[0].dot(D);
503 // << 363 G4double e = E[1].dot(D);
504 // This method has been heavily modified thank << 364 G4double s = b*e - c*d;
505 // corrections of Rickard Holmberg. << 365 G4double t = b*d - a*e;*/
506 // <<
507 G4double G4TriangularFacet::Distance (const G4 <<
508 G4 <<
509 const G4 <<
510 { <<
511 // <<
512 // Start with quicky test to determine if th <<
513 // the triangle is any closer to p than minD <<
514 // about more accurate test. <<
515 // <<
516 G4double dist = kInfinity; 366 G4double dist = kInfinity;
517 if ((p-fCircumcentre).mag()-fRadius < minDis << 367
>> 368 /*if (s+t > 1.0 || s < 0.0 || t < 0.0)
518 { 369 {
519 // << 370 G4ThreeVector D0 = P0 - p;
520 // It's possible that the triangle is clos << 371 G4ThreeVector D1 = P[0] - p;
521 // so do more accurate assessment. << 372 G4ThreeVector D2 = P[1] - p;
522 // << 373
523 G4ThreeVector v = Distance(p); << 374 G4double d0 = D0.mag();
524 G4double dist1 = sqrt(fSqrDist); << 375 G4double d1 = D1.mag();
525 G4double dir = v.dot(fSurfaceNormal); << 376 G4double d2 = D2.mag();
526 G4bool wrongSide = (dir > 0.0 && !outgoing << 377
527 if (dist1 <= kCarTolerance) << 378 dist = min(d0, min(d1, d2));
528 { << 379 if (dist > minDist ||
529 // << 380 (D0.dot(surfaceNormal) > 0.0 && !outgoing) ||
530 // Point p is very close to triangle. C << 381 (D0.dot(surfaceNormal) < 0.0 && outgoing)) return kInfinity;
531 // in which case return distance of 0.0 << 382 }*/
532 // << 383
533 if (wrongSide) dist = 0.0; << 384 G4ThreeVector v = Distance(p);
534 else dist = dist1; << 385 G4double dir = v.dot(surfaceNormal);
535 } << 386 if ((dir > dirTolerance && !outgoing) ||
536 else if (!wrongSide) dist = dist1; << 387 (dir <-dirTolerance && outgoing)) dist = kInfinity;
537 } << 388 else dist = v.mag();
>> 389
538 return dist; 390 return dist;
539 } 391 }
540 392
541 ////////////////////////////////////////////// 393 ///////////////////////////////////////////////////////////////////////////////
542 // 394 //
543 // Extent <<
544 // <<
545 // Calculates the furthest the triangle extend <<
546 // defined by the vector axis. <<
547 // <<
548 G4double G4TriangularFacet::Extent (const G4Th 395 G4double G4TriangularFacet::Extent (const G4ThreeVector axis)
549 { 396 {
550 G4double ss = GetVertex(0).dot(axis); << 397 G4double s = P0.dot(axis);
551 G4double sp = GetVertex(1).dot(axis); << 398 G4double sp = P[0].dot(axis);
552 if (sp > ss) ss = sp; << 399 if (sp > s) s = sp;
553 sp = GetVertex(2).dot(axis); << 400 sp = P[1].dot(axis);
554 if (sp > ss) ss = sp; << 401 if (sp > s) s = sp;
555 return ss; << 402
>> 403 return s;
556 } 404 }
557 405
558 ////////////////////////////////////////////// 406 ///////////////////////////////////////////////////////////////////////////////
559 // 407 //
560 // Intersect << 408 G4bool G4TriangularFacet::Intersect (const G4ThreeVector &p,
561 // << 409 const G4ThreeVector &v, G4bool outgoing, G4double &distance,
562 // Member function to find the next intersecti << 410 G4double &distFromSurface, G4ThreeVector &normal)
563 // direction of v. If: << 411 {
564 // (1) "outgoing" is TRUE, only consider the f << 412 G4ThreeVector D = P0 - p;
565 // the face. << 413 G4double d = E[0].dot(D);
566 // (2) "outgoing" is FALSE, only consider the << 414 G4double e = E[1].dot(D);
567 // the face. << 415 G4double g = E[0].dot(v);
568 // Member functions returns TRUE if there is a << 416 G4double h = E[1].dot(v);
569 // Sets the distance (distance along w), distF << 417 G4double q = D.dot(v);
570 // and normal. << 418
571 // << 419 G4double A00 = a - g*g;
572 // Also considers intersections that happen wi << 420 G4double A11 = c - h*h;
573 // distances of distFromSurface = 0.5*kCarTole << 421 G4double A01 = b - g*h;
574 // This is to detect kSurface without doing fA << 422 G4double det2 = A00*A11 - A01*A01;
575 // G4TessellatedSolid::Distance(p,v) calculati << 423
576 // << 424 G4double s = kInfinity;
577 // This member function is thanks the valuable << 425 G4double t = kInfinity;
578 // However, "gotos" are the Work of the Devil << 426
579 // extreme prejudice!! << 427 G4double dist = kInfinity;
580 // << 428 G4bool intersect = false;
581 // IMPORTANT NOTE: These calculations are pre << 429 G4double normalComp = 0.0;
582 // vector. If G4TessellatedSolid or other cla << 430
583 // with |v| != 1 then there will be errors. << 431
584 // << 432 if (det2 != 0.0)
585 G4bool G4TriangularFacet::Intersect (const G4T << 433 {
586 const G4T << 434 G4double B0 = q*g - d;
587 G4b << 435 G4double B1 = q*h - e;
588 G4d << 436 s = (A11*B0 - A01*B1)/det2;
589 G4d << 437 if ((s >= sMin) && (s <= sMax))
590 G4T << 438 {
591 { << 439 t = (A00*B1 - A01*B0)/det2;
592 // << 440 if ((t >= tMin) && (t < 1.0 - s + std::fabs(sMin)))
593 // Check whether the direction of the facet << 441 { //THIS IS A FUDGE FOR THE MOMENT
594 // and the need to be outgoing or ingoing. << 442 dist = q + g*s + h*t;
595 // return false. << 443 normalComp = v.dot(surfaceNormal);
596 // << 444 // intersect = (dist >= 0.0 &&
597 G4double w = v.dot(fSurfaceNormal); << 445 // ((outgoing && normalComp > 0.0) || (!outgoing && normalComp < 0.0)));
598 if ((outgoing && w < -dirTolerance) || (!out << 446 intersect = (dist >= -kCarTolerance*0.5 &&
599 { << 447 ((outgoing && normalComp > dirTolerance) ||
600 distance = kInfinity; << 448 (!outgoing && normalComp <-dirTolerance))); //FUDGE FOR THE MOMENT
601 distFromSurface = kInfinity; <<
602 normal.set(0,0,0); <<
603 return false; <<
604 } <<
605 // <<
606 // Calculate the orthogonal distance from p <<
607 // triangle. Then determine if we're on the <<
608 // surface (at fA distance greater than kCar <<
609 // "outgoing". <<
610 // <<
611 const G4ThreeVector& p0 = GetVertex(0); <<
612 G4ThreeVector D = p0 - p; <<
613 distFromSurface = D.dot(fSurfaceNormal); <<
614 G4bool wrongSide = (outgoing && distFromSurf <<
615 (!outgoing && distFromSurface > 0.5*kCarT <<
616 <<
617 if (wrongSide) <<
618 { <<
619 distance = kInfinity; <<
620 distFromSurface = kInfinity; <<
621 normal.set(0,0,0); <<
622 return false; <<
623 } <<
624 <<
625 wrongSide = (outgoing && distFromSurface < 0 <<
626 || (!outgoing && distFromSurface > <<
627 if (wrongSide) <<
628 { <<
629 // <<
630 // We're slightly on the wrong side of the <<
631 // enough using fA precise distance calcul <<
632 // <<
633 G4ThreeVector u = Distance(p); <<
634 if (fSqrDist <= kCarTolerance*kCarToleranc <<
635 { <<
636 // <<
637 // We're very close. Therefore return f <<
638 // to pretend we intersect. <<
639 // <<
640 // distance = -0.5*kCarTolerance <<
641 distance = 0.0; <<
642 normal = fSurfaceNormal; <<
643 return true; <<
644 } <<
645 else <<
646 { <<
647 // <<
648 // We're close to the surface containing <<
649 // far from the triangle, and on the wro <<
650 // of the surface normal and v. There i <<
651 // <<
652 distance = kInfinity; <<
653 distFromSurface = kInfinity; <<
654 normal.set(0,0,0); <<
655 return false; <<
656 } <<
657 } <<
658 if (w < dirTolerance && w > -dirTolerance) <<
659 { <<
660 // <<
661 // The ray is within the plane of the tria <<
662 // in the plane of the triangle. First try <<
663 // mu and nu, where mu is fE1/|fE1|. This <<
664 // the original algorithm due to Rickard H <<
665 // mathematical justification than the ori <<
666 // beware Rickard's was less time-consumin <<
667 // <<
668 // Note that vprime is not fA unit vector. <<
669 // since the values of distance along vpri <<
670 // with the triangle will be used to deter <<
671 // same time. <<
672 // <<
673 G4ThreeVector mu = fE1.unit(); <<
674 G4ThreeVector nu = fSurfaceNormal.cross(mu <<
675 G4TwoVector pprime(p.dot(mu), p.dot(nu)); <<
676 G4TwoVector vprime(v.dot(mu), v.dot(nu)); <<
677 G4TwoVector P0prime(p0.dot(mu), p0.dot(nu) <<
678 G4TwoVector E0prime(fE1.mag(), 0.0); <<
679 G4TwoVector E1prime(fE2.dot(mu), fE2.dot(n <<
680 G4TwoVector loc[2]; <<
681 if (G4TessellatedGeometryAlgorithms::Inter <<
682 vprime, P0 <<
683 { <<
684 // <<
685 // There is an intersection between the <<
686 // Now check which part of the line inte <<
687 // containing the triangle in 3D. <<
688 // <<
689 G4double vprimemag = vprime.mag(); <<
690 G4double s0 = (loc[0] - pprime).m <<
691 G4double s1 = (loc[1] - pprime).m <<
692 G4double normDist0 = fSurfaceNormal.dot( <<
693 G4double normDist1 = fSurfaceNormal.dot( <<
694 <<
695 if ((normDist0 < 0.0 && normDist1 < 0.0) <<
696 || (normDist0 > 0.0 && normDist1 > 0.0) <<
697 || (normDist0 == 0.0 && normDist1 == 0. <<
698 { <<
699 distance = kInfinity; <<
700 distFromSurface = kInfinity; <<
701 normal.set(0,0,0); <<
702 return false; <<
703 } <<
704 else <<
705 { <<
706 G4double dnormDist = normDist1 - normD <<
707 if (fabs(dnormDist) < DBL_EPSILON) <<
708 { <<
709 distance = s0; <<
710 normal = fSurfaceNormal; <<
711 if (!outgoing) distFromSurface = -di <<
712 return true; <<
713 } <<
714 else <<
715 { <<
716 distance = s0 - normDist0*(s1-s0)/dn <<
717 normal = fSurfaceNormal; <<
718 if (!outgoing) distFromSurface = -di <<
719 return true; <<
720 } <<
721 } 449 }
722 } 450 }
723 else <<
724 { <<
725 distance = kInfinity; <<
726 distFromSurface = kInfinity; <<
727 normal.set(0,0,0); <<
728 return false; <<
729 } <<
730 } 451 }
731 // << 452
732 // << 453 if (intersect)
733 // Use conventional algorithm to determine t << 454 {
734 // intersection. This involves determining << 455 if (dist < kCarTolerance * 0.5) { dist = 0.0; }
735 // line with the plane containing the triang << 456 distance = dist;
736 // point is within the triangle. << 457 distFromSurface = dist * normalComp;
737 // << 458 normal = surfaceNormal;
738 distance = distFromSurface / w; <<
739 G4ThreeVector pp = p + v*distance; <<
740 G4ThreeVector DD = p0 - pp; <<
741 G4double d = fE1.dot(DD); <<
742 G4double e = fE2.dot(DD); <<
743 G4double ss = fB*e - fC*d; <<
744 G4double t = fB*d - fA*e; <<
745 <<
746 G4double sTolerance = (fabs(fB)+ fabs(fC) + <<
747 G4double tTolerance = (fabs(fA)+ fabs(fB) + <<
748 G4double detTolerance = (fabs(fA)+ fabs(fC) <<
749 <<
750 //if (ss < 0.0 || t < 0.0 || ss+t > fDet) <<
751 if (ss < -sTolerance || t < -tTolerance || ( <<
752 { <<
753 // <<
754 // The intersection is outside of the tria <<
755 // <<
756 distance = distFromSurface = kInfinity; <<
757 normal.set(0,0,0); <<
758 return false; <<
759 } 459 }
760 else 460 else
761 { 461 {
762 // << 462 distance = kInfinity;
763 // There is an intersection. Now we only << 463 distFromSurface = kInfinity;
764 // << 464 normal = G4ThreeVector(0.0,0.0,0.0);
765 normal = fSurfaceNormal; <<
766 if (!outgoing) distFromSurface = -distFrom <<
767 return true; <<
768 } 465 }
>> 466
>> 467 return intersect;
769 } 468 }
770 469
771 ////////////////////////////////////////////// 470 ////////////////////////////////////////////////////////////////////////
772 // 471 //
773 // GetPointOnFace 472 // GetPointOnFace
774 // 473 //
775 // Auxiliary method, returns a uniform random << 474 // Auxiliary method for get a random point on surface
776 // << 475
777 G4ThreeVector G4TriangularFacet::GetPointOnFac 476 G4ThreeVector G4TriangularFacet::GetPointOnFace() const
778 { 477 {
779 G4double u = G4UniformRand(); << 478 G4double lambda1,lambda2;
780 G4double v = G4UniformRand(); << 479 G4ThreeVector v, w;
781 if (u+v > 1.) { u = 1. - u; v = 1. - v; } << 480
782 return GetVertex(0) + u*fE1 + v*fE2; << 481 v = P[1] - P[0];
>> 482 w = P[0] - P0;
>> 483
>> 484 lambda1 = CLHEP::RandFlat::shoot(0.,1.);
>> 485 lambda2 = CLHEP::RandFlat::shoot(0.,lambda1);
>> 486
>> 487 return (P0 + lambda1*w + lambda2*v);
783 } 488 }
784 489
785 ////////////////////////////////////////////// 490 ////////////////////////////////////////////////////////////////////////
786 // 491 //
787 // GetArea 492 // GetArea
788 // 493 //
789 // Auxiliary method for returning the surface << 494 // Auxiliary method for returning the surface area
790 // <<
791 G4double G4TriangularFacet::GetArea() const <<
792 { <<
793 return fArea; <<
794 } <<
795 495
796 ////////////////////////////////////////////// << 496 G4double G4TriangularFacet::GetArea()
797 // <<
798 G4GeometryType G4TriangularFacet::GetEntityTyp <<
799 { 497 {
800 return "G4TriangularFacet"; << 498 if (area) { return area; }
801 } <<
802 499
803 ////////////////////////////////////////////// << 500 G4ThreeVector v, w;
804 // <<
805 G4ThreeVector G4TriangularFacet::GetSurfaceNor <<
806 { <<
807 return fSurfaceNormal; <<
808 } <<
809 501
810 ////////////////////////////////////////////// << 502 v = P[1] - P[0];
811 // << 503 w = P[0] - P0;
812 void G4TriangularFacet::SetSurfaceNormal (cons << 504 area = 0.5*(v.cross(w)).mag();
813 { << 505
814 fSurfaceNormal = normal; << 506 return area;
815 } 507 }
816 508