| 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 to DEFCON 705 IPR conditions. 20 // * subject to DEFCON 705 IPR conditions. *
21 // * By using, copying, modifying or distri 21 // * By using, copying, modifying or distributing the software (or *
22 // * any work based on the software) you ag 22 // * any work based on the software) you agree to acknowledge its *
23 // * use in resulting scientific publicati 23 // * use in resulting scientific publications, and indicate your *
24 // * acceptance of all terms of the Geant4 Sof 24 // * acceptance of all terms of the Geant4 Software license. *
25 // ******************************************* 25 // ********************************************************************
26 // 26 //
27 // G4QuadrangularFacet class implementation. <<
28 // 27 //
29 // 31 October 2004 P R Truscott, QinetiQ Ltd, << 28 // $Id: G4QuadrangularFacet.cc,v 1.7 2008/12/18 12:57:34 gunter Exp $
30 // 12 October 2012 M Gayer, CERN << 29 // GEANT4 tag $Name: geant4-09-02 $
31 // New implementation reducing << 30 //
32 // and considerable CPU speedu << 31 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33 // implementation of G4Tessell << 32 //
34 // 29 February 2016 E Tcherniaev, CERN << 33 // MODULE: G4QuadrangularFacet.cc
35 // Added exhaustive tests to c << 34 //
36 // quadrangular facet: colline << 35 // Date: 15/06/2005
37 // degenerate, concave or self << 36 // Author: P R Truscott
38 // ------------------------------------------- << 37 // Organisation: QinetiQ Ltd, UK
>> 38 // Customer: UK Ministry of Defence : RAO CRP TD Electronic Systems
>> 39 // Contract: C/MAT/N03517
>> 40 //
>> 41 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>> 42 //
>> 43 // CHANGE HISTORY
>> 44 // --------------
>> 45 //
>> 46 // 31 October 2004, P R Truscott, QinetiQ Ltd, UK - Created.
>> 47 //
>> 48 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
39 49
40 #include "G4QuadrangularFacet.hh" 50 #include "G4QuadrangularFacet.hh"
41 #include "geomdefs.hh" << 51 #include "globals.hh"
42 #include "Randomize.hh" 52 #include "Randomize.hh"
43 <<
44 using namespace std; <<
45 53
46 ////////////////////////////////////////////// 54 ///////////////////////////////////////////////////////////////////////////////
47 // 55 //
48 // Constructing two adjacent G4TriangularFacet << 56 // !!!THIS IS A FUDGE!!! IT'S TWO ADJACENT G4TRIANGULARFACETS
49 // Not efficient, but practical... << 57 // --- NOT EFFICIENT BUT PRACTICAL.
50 // 58 //
51 G4QuadrangularFacet::G4QuadrangularFacet (cons << 59 G4QuadrangularFacet::G4QuadrangularFacet (const G4ThreeVector Pt0,
52 cons << 60 const G4ThreeVector vt1, const G4ThreeVector vt2,
53 cons << 61 const G4ThreeVector vt3, G4FacetVertexType vertexType)
54 cons << 62 : G4VFacet()
55 <<
56 { 63 {
57 G4double delta = 1.0 * kCarTolerance; // << 64 P0 = Pt0;
58 G4double epsilon = 0.01 * kCarTolerance; // << 65 nVertices = 4;
59 <<
60 G4ThreeVector e1, e2, e3; <<
61 SetVertex(0, vt0); <<
62 if (vertexType == ABSOLUTE) 66 if (vertexType == ABSOLUTE)
63 { 67 {
64 SetVertex(1, vt1); << 68 P.push_back(vt1);
65 SetVertex(2, vt2); << 69 P.push_back(vt2);
66 SetVertex(3, vt3); << 70 P.push_back(vt3);
67 << 71
68 e1 = vt1 - vt0; << 72 E.push_back(vt1 - P0);
69 e2 = vt2 - vt0; << 73 E.push_back(vt2 - P0);
70 e3 = vt3 - vt0; << 74 E.push_back(vt3 - P0);
71 } 75 }
72 else 76 else
73 { 77 {
74 SetVertex(1, vt0 + vt1); << 78 P.push_back(P0 + vt1);
75 SetVertex(2, vt0 + vt2); << 79 P.push_back(P0 + vt2);
76 SetVertex(3, vt0 + vt3); << 80 P.push_back(P0 + vt3);
77 << 81
78 e1 = vt1; << 82 E.push_back(vt1);
79 e2 = vt2; << 83 E.push_back(vt2);
80 e3 = vt3; << 84 E.push_back(vt3);
81 } <<
82 <<
83 // Check length of sides and diagonals <<
84 // <<
85 G4double leng1 = e1.mag(); <<
86 G4double leng2 = (e2-e1).mag(); <<
87 G4double leng3 = (e3-e2).mag(); <<
88 G4double leng4 = e3.mag(); <<
89 <<
90 G4double diag1 = e2.mag(); <<
91 G4double diag2 = (e3-e1).mag(); <<
92 <<
93 if (leng1 <= delta || leng2 <= delta || leng <<
94 diag1 <= delta || diag2 <= delta) <<
95 { <<
96 ostringstream message; <<
97 message << "Sides/diagonals of facet are t <<
98 << "P0 = " << GetVertex(0) << G4en <<
99 << "P1 = " << GetVertex(1) << G4en <<
100 << "P2 = " << GetVertex(2) << G4en <<
101 << "P3 = " << GetVertex(3) << G4en <<
102 << "Side1 length (P0->P1) = " << l <<
103 << "Side2 length (P1->P2) = " << l <<
104 << "Side3 length (P2->P3) = " << l <<
105 << "Side4 length (P3->P0) = " << l <<
106 << "Diagonal1 length (P0->P2) = " <<
107 << "Diagonal2 length (P1->P3) = " <<
108 G4Exception("G4QuadrangularFacet::G4Quadra <<
109 "GeomSolids1001", JustWarning, <<
110 return; <<
111 } <<
112 <<
113 // Check that vertices are not collinear <<
114 // <<
115 G4double s1 = (e1.cross(e2)).mag()*0.5; <<
116 G4double s2 = ((e2-e1).cross(e3-e2)).mag()*0 <<
117 G4double s3 = (e2.cross(e3)).mag()*0.5; <<
118 G4double s4 = (e1.cross(e3)).mag()*0.5; <<
119 <<
120 G4double h1 = 2.*s1 / std::max(std::max(leng <<
121 G4double h2 = 2.*s2 / std::max(std::max(leng <<
122 G4double h3 = 2.*s3 / std::max(std::max(leng <<
123 G4double h4 = 2.*s4 / std::max(std::max(leng <<
124 <<
125 if (h1 <= delta || h2 <= delta || h3 <= delt <<
126 { <<
127 ostringstream message; <<
128 message << "Facet has three or more collin <<
129 << "P0 = " << GetVertex(0) << G4en <<
130 << "P1 = " << GetVertex(1) << G4en <<
131 << "P2 = " << GetVertex(2) << G4en <<
132 << "P3 = " << GetVertex(3) << G4en <<
133 << "Smallest heights:" << G4endl <<
134 << " in triangle P0-P1-P2 = " << <<
135 << " in triangle P1-P2-P3 = " << <<
136 << " in triangle P2-P3-P0 = " << <<
137 << " in triangle P3-P0-P1 = " << <<
138 G4Exception("G4QuadrangularFacet::G4Quadra <<
139 "GeomSolids1001", JustWarning, messa <<
140 return; <<
141 } 85 }
142 86
143 // Check that vertices are coplanar by compu << 87 G4double length1 = E[0].mag();
144 // height of tetrahedron comprising of verti << 88 G4double length2 = (P[1]-P[0]).mag();
145 // << 89 G4double length3 = (P[2]-P[1]).mag();
146 G4double smax = std::max( std::max(s1,s2), s << 90 G4double length4 = E[2].mag();
147 G4double hmin = 0.5 * std::fabs( e1.dot(e2.c << 91
148 if (hmin >= epsilon) << 92 G4ThreeVector normal1 = E[0].cross(E[1]).unit();
149 { << 93 G4ThreeVector normal2 = E[1].cross(E[2]).unit();
150 ostringstream message; <<
151 message << "Facet is not planar." << G4end <<
152 << "Disrepancy = " << hmin << G4en <<
153 << "P0 = " << GetVertex(0) << G4en <<
154 << "P1 = " << GetVertex(1) << G4en <<
155 << "P2 = " << GetVertex(2) << G4en <<
156 << "P3 = " << GetVertex(3); <<
157 G4Exception("G4QuadrangularFacet::G4Quadra <<
158 "GeomSolids1001", JustWarning, messa <<
159 return; <<
160 } <<
161 94
162 // Check that facet is convex by computing c << 95 if (length1 <= kCarTolerance || length2 <= kCarTolerance ||
163 // of diagonals << 96 length3 <= kCarTolerance || length4 <= kCarTolerance ||
164 // << 97 normal1.dot(normal2) < 0.9999999999)
165 G4ThreeVector normal = e2.cross(e3-e1); <<
166 G4double s = kInfinity, t = kInfinity, magni <<
167 if (magnitude2 > delta*delta) // check: magn <<
168 { 98 {
169 s = normal.dot(e1.cross(e3-e1)) / magnitud << 99 G4Exception("G4QuadrangularFacet::G4QuadrangularFacet()",
170 t = normal.dot(e1.cross(e2)) / magnitud << 100 "InvalidSetup", JustWarning,
>> 101 "Length of sides of facet are too small or sides not planar.");
>> 102 G4cerr << G4endl;
>> 103 G4cerr << "P0 = " << P0 << G4endl;
>> 104 G4cerr << "P1 = " << P[0] << G4endl;
>> 105 G4cerr << "P2 = " << P[1] << G4endl;
>> 106 G4cerr << "P3 = " << P[2] << G4endl;
>> 107 G4cerr << "Side lengths = P0->P1" << length1 << G4endl;
>> 108 G4cerr << "Side lengths = P1->P2" << length2 << G4endl;
>> 109 G4cerr << "Side lengths = P2->P3" << length3 << G4endl;
>> 110 G4cerr << "Side lengths = P3->P0" << length4 << G4endl;
>> 111 G4cerr << G4endl;
>> 112
>> 113 isDefined = false;
>> 114 geometryType = "G4QuadragularFacet";
>> 115 surfaceNormal = G4ThreeVector(0.0,0.0,0.0);
171 } 116 }
172 if (s <= 0. || s >= 1. || t <= 0. || t >= 1. << 117 else
173 { 118 {
174 ostringstream message; << 119 isDefined = true;
175 message << "Facet is not convex." << G4end << 120 geometryType = "G4QuadrangularFacet";
176 << "Parameters of crosspoint of di << 121
177 << s << " and " << t << G4endl << 122 facet1 = new G4TriangularFacet(P0,P[0],P[1],ABSOLUTE);
178 << "should both be within (0,1) ra << 123 facet2 = new G4TriangularFacet(P0,P[1],P[2],ABSOLUTE);
179 << "P0 = " << GetVertex(0) << G4endl << 124 surfaceNormal = normal1;
180 << "P1 = " << GetVertex(1) << G4endl << 125
181 << "P2 = " << GetVertex(2) << G4endl << 126 G4ThreeVector vtmp = 0.5 * (E[0] + E[1]);
182 << "P3 = " << GetVertex(3); << 127 circumcentre = P0 + vtmp;
183 G4Exception("G4QuadrangularFacet::G4Quadra << 128 radiusSqr = vtmp.mag2();
184 "GeomSolids1001", JustWarning, messa << 129 radius = std::sqrt(radiusSqr);
185 return; << 130
>> 131 for (size_t i=0; i<4; i++) I.push_back(0);
186 } 132 }
187 <<
188 // Define facet <<
189 // <<
190 fFacet1 = G4TriangularFacet(GetVertex(0),Get <<
191 fFacet2 = G4TriangularFacet(GetVertex(0),Get <<
192 <<
193 normal = normal.unit(); <<
194 fFacet1.SetSurfaceNormal(normal); <<
195 fFacet2.SetSurfaceNormal(normal); <<
196 <<
197 G4ThreeVector vtmp = 0.5 * (e1 + e2); <<
198 fCircumcentre = GetVertex(0) + vtmp; <<
199 G4double radiusSqr = vtmp.mag2(); <<
200 fRadius = std::sqrt(radiusSqr); <<
201 // 29.02.2016 Remark by E.Tcherniaev: comput <<
202 // of fCircumcenter and fRadius is wrong, ho <<
203 // it did not create any problem till now. <<
204 // Bizarre! Need to investigate! <<
205 } 133 }
206 134
207 ////////////////////////////////////////////// 135 ///////////////////////////////////////////////////////////////////////////////
208 // 136 //
209 G4QuadrangularFacet::~G4QuadrangularFacet () = << 137 G4QuadrangularFacet::~G4QuadrangularFacet ()
210 <<
211 ////////////////////////////////////////////// <<
212 // <<
213 G4QuadrangularFacet::G4QuadrangularFacet (cons <<
214 : G4VFacet(rhs) <<
215 { 138 {
216 fFacet1 = rhs.fFacet1; << 139 if (facet1 != 0) delete facet1;
217 fFacet2 = rhs.fFacet2; << 140 if (facet2 != 0) delete facet2;
218 fRadius = 0.0; << 141
219 } << 142 P.clear();
220 << 143 E.clear();
221 ////////////////////////////////////////////// << 144 I.clear();
222 // <<
223 G4QuadrangularFacet & <<
224 G4QuadrangularFacet::operator=(const G4Quadran <<
225 { <<
226 if (this == &rhs) return *this; <<
227 <<
228 fFacet1 = rhs.fFacet1; <<
229 fFacet2 = rhs.fFacet2; <<
230 fRadius = 0.0; <<
231 <<
232 return *this; <<
233 } 145 }
234 146
235 ////////////////////////////////////////////// 147 ///////////////////////////////////////////////////////////////////////////////
236 // 148 //
237 G4VFacet* G4QuadrangularFacet::GetClone () << 149 G4VFacet *G4QuadrangularFacet::GetClone ()
238 { 150 {
239 auto c = new G4QuadrangularFacet (GetVertex( << 151 G4QuadrangularFacet *c =
240 GetVertex( << 152 new G4QuadrangularFacet (P0, P[0], P[1], P[2], ABSOLUTE);
241 return c; << 153 G4VFacet *cc = 0;
>> 154 cc = c;
>> 155 return cc;
242 } 156 }
243 157
244 ////////////////////////////////////////////// 158 ///////////////////////////////////////////////////////////////////////////////
245 // 159 //
246 G4ThreeVector G4QuadrangularFacet::Distance (c << 160 G4ThreeVector G4QuadrangularFacet::Distance (const G4ThreeVector &p)
247 { 161 {
248 G4ThreeVector v1 = fFacet1.Distance(p); << 162 G4ThreeVector v1 = facet1->Distance(p);
249 G4ThreeVector v2 = fFacet2.Distance(p); << 163 G4ThreeVector v2 = facet2->Distance(p);
250 << 164
251 if (v1.mag2() < v2.mag2()) return v1; 165 if (v1.mag2() < v2.mag2()) return v1;
252 else return v2; 166 else return v2;
253 } 167 }
254 168
255 ////////////////////////////////////////////// 169 ///////////////////////////////////////////////////////////////////////////////
256 // 170 //
257 G4double G4QuadrangularFacet::Distance (const << 171 G4double G4QuadrangularFacet::Distance (const G4ThreeVector &p,
258 << 172 const G4double)
259 { << 173 {
260 G4double dist = Distance(p).mag(); << 174 /*G4ThreeVector D = P0 - p;
>> 175 G4double d = E[0].dot(D);
>> 176 G4double e = E[1].dot(D);
>> 177 G4double s = b*e - c*d;
>> 178 G4double t = b*d - a*e;*/
>> 179 G4double dist = kInfinity;
>> 180
>> 181 /*if (s+t > 1.0 || s < 0.0 || t < 0.0)
>> 182 {
>> 183 G4ThreeVector D0 = P0 - p;
>> 184 G4ThreeVector D1 = P[0] - p;
>> 185 G4ThreeVector D2 = P[1] - p;
>> 186
>> 187 G4double d0 = D0.mag();
>> 188 G4double d1 = D1.mag();
>> 189 G4double d2 = D2.mag();
>> 190
>> 191 dist = min(d0, min(d1, d2));
>> 192 if (dist > minDist) return kInfinity;
>> 193 }*/
>> 194
>> 195 dist = Distance(p).mag();
>> 196
261 return dist; 197 return dist;
262 } 198 }
263 199
264 ////////////////////////////////////////////// 200 ///////////////////////////////////////////////////////////////////////////////
265 // 201 //
266 G4double G4QuadrangularFacet::Distance (const << 202 G4double G4QuadrangularFacet::Distance (const G4ThreeVector &p,
267 const << 203 const G4double, const G4bool outgoing)
268 { 204 {
269 G4double dist; << 205 /*G4ThreeVector D = P0 - p;
270 << 206 G4double d = E[0].dot(D);
>> 207 G4double e = E[1].dot(D);
>> 208 G4double s = b*e - c*d;
>> 209 G4double t = b*d - a*e;*/
>> 210 G4double dist = kInfinity;
>> 211
>> 212 /*if (s+t > 1.0 || s < 0.0 || t < 0.0)
>> 213 {
>> 214 G4ThreeVector D0 = P0 - p;
>> 215 G4ThreeVector D1 = P[0] - p;
>> 216 G4ThreeVector D2 = P[1] - p;
>> 217
>> 218 G4double d0 = D0.mag();
>> 219 G4double d1 = D1.mag();
>> 220 G4double d2 = D2.mag();
>> 221
>> 222 dist = min(d0, min(d1, d2));
>> 223 if (dist > minDist ||
>> 224 (D0.dot(surfaceNormal) > 0.0 && !outgoing) ||
>> 225 (D0.dot(surfaceNormal) < 0.0 && outgoing)) return kInfinity;
>> 226 }*/
>> 227
271 G4ThreeVector v = Distance(p); 228 G4ThreeVector v = Distance(p);
272 G4double dir = v.dot(GetSurfaceNormal()); << 229 G4double dir = v.dot(surfaceNormal);
273 if ( ((dir > dirTolerance) && (!outgoing)) << 230 if ((dir > dirTolerance && !outgoing) ||
274 || ((dir < -dirTolerance) && outgoing)) << 231 (dir <-dirTolerance && outgoing)) dist = kInfinity;
275 dist = kInfinity; << 232 else dist = v.mag();
276 else << 233
277 dist = v.mag(); <<
278 return dist; 234 return dist;
279 } 235 }
280 236
281 ////////////////////////////////////////////// 237 ///////////////////////////////////////////////////////////////////////////////
282 // 238 //
283 G4double G4QuadrangularFacet::Extent (const G4 239 G4double G4QuadrangularFacet::Extent (const G4ThreeVector axis)
284 { 240 {
285 G4double ss = 0; << 241 G4double s = P0.dot(axis);
286 << 242 for (G4ThreeVectorList::iterator it=P.begin(); it!=P.end(); it++)
287 for (G4int i = 0; i <= 3; ++i) <<
288 { 243 {
289 G4double sp = GetVertex(i).dot(axis); << 244 G4double sp = it->dot(axis);
290 if (sp > ss) ss = sp; << 245 if (sp > s) s = sp;
291 } 246 }
292 return ss; << 247
>> 248 return s;
293 } 249 }
294 250
295 ////////////////////////////////////////////// 251 ///////////////////////////////////////////////////////////////////////////////
296 // 252 //
297 G4bool G4QuadrangularFacet::Intersect (const G << 253 G4bool G4QuadrangularFacet::Intersect (const G4ThreeVector &p,
298 const G << 254 const G4ThreeVector &v, G4bool outgoing, G4double &distance,
299 G << 255 G4double &distFromSurface, G4ThreeVector &normal)
300 G <<
301 G <<
302 G <<
303 { 256 {
304 G4bool intersect = 257 G4bool intersect =
305 fFacet1.Intersect(p,v,outgoing,distance,di << 258 facet1->Intersect(p,v,outgoing,distance,distFromSurface,normal);
306 if (!intersect) intersect = << 259 if (!intersect)
307 fFacet2.Intersect(p,v,outgoing,distance,di << 260 {
>> 261 intersect = facet2->Intersect(p,v,outgoing,distance,distFromSurface,normal);
>> 262 }
>> 263
308 if (!intersect) 264 if (!intersect)
309 { 265 {
310 distance = distFromSurface = kInfinity; << 266 distance = kInfinity;
311 normal.set(0,0,0); << 267 distFromSurface = kInfinity;
>> 268 normal = G4ThreeVector(0.0,0.0,0.0);
312 } 269 }
>> 270
313 return intersect; 271 return intersect;
314 } 272 }
315 273
316 ////////////////////////////////////////////// << 274 ////////////////////////////////////////////////////////////////////////
317 // 275 //
318 // Auxiliary method to get a uniform random po << 276 // GetPointOnFace
319 // 277 //
>> 278 // Auxiliary method for get a random point on surface
>> 279
320 G4ThreeVector G4QuadrangularFacet::GetPointOnF 280 G4ThreeVector G4QuadrangularFacet::GetPointOnFace() const
321 { 281 {
322 G4double s1 = fFacet1.GetArea(); << 282 G4ThreeVector pr;
323 G4double s2 = fFacet2.GetArea(); << 283
324 return ((s1+s2)*G4UniformRand() < s1) ? << 284 if ( G4UniformRand() < 0.5 )
325 fFacet1.GetPointOnFace() : fFacet2.GetPoin << 285 {
>> 286 pr = facet1->GetPointOnFace();
>> 287 }
>> 288 else
>> 289 {
>> 290 pr = facet2->GetPointOnFace();
>> 291 }
>> 292
>> 293 return pr;
326 } 294 }
327 295
328 ////////////////////////////////////////////// << 296 ////////////////////////////////////////////////////////////////////////
329 // 297 //
330 // Auxiliary method for returning the surface << 298 // GetArea
331 // 299 //
332 G4double G4QuadrangularFacet::GetArea() const << 300 // Auxiliary method for returning the surface area
333 { <<
334 G4double area = fFacet1.GetArea() + fFacet2. <<
335 return area; <<
336 } <<
337 301
338 ////////////////////////////////////////////// << 302 G4double G4QuadrangularFacet::GetArea()
339 // <<
340 G4String G4QuadrangularFacet::GetEntityType () <<
341 { 303 {
342 return "G4QuadrangularFacet"; << 304 if (!area) { area = facet1->GetArea() + facet2->GetArea(); }
343 } <<
344 305
345 ////////////////////////////////////////////// << 306 return area;
346 // <<
347 G4ThreeVector G4QuadrangularFacet::GetSurfaceN <<
348 { <<
349 return fFacet1.GetSurfaceNormal(); <<
350 } 307 }
351 308