| 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 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 // G4TessellatedSolid implementation << 27 // $Id: G4TessellatedSolid.cc,v 1.4 2006/06/29 18:48:57 gunter Exp $
>> 28 // GEANT4 tag $Name: geant4-08-01 $
28 // 29 //
29 // 31.10.2004, P R Truscott, QinetiQ Ltd, UK - << 30 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
30 // 17.09.2007, P R Truscott, QinetiQ Ltd & Ric << 31 //
31 // Updated extensively prio << 32 // MODULE: G4TessellatedSolid.cc
32 // concaved tessellated sur << 33 //
33 // of Richard Holmberg. Th << 34 // Date: 15/06/2005
34 // to determine with inside << 35 // Author: P R Truscott
35 // random rays from the poi << 36 // Organisation: QinetiQ Ltd, UK
36 // are predefined rather th << 37 // Customer: UK Ministry of Defence : RAO CRP TD Electronic Systems
37 // number generator at run- << 38 // Contract: C/MAT/N03517
38 // 12.10.2012, M Gayer, CERN, complete rewrite << 39 //
39 // requirements more than 5 << 40 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
40 // tens or more depending o << 41 //
41 // to voxelization of surfa << 42 // CHANGE HISTORY
42 // Speedup factor of thousa << 43 // --------------
43 // facets in hundreds of th << 44 //
44 // 23.10.2016, E Tcherniaev, reimplemented Cal << 45 // 22 November 2005, F Lei
45 // use of G4BoundingEnvelop << 46 // - Changed ::DescribeYourselfTo(), line 464
46 // ------------------------------------------- << 47 // - added GetPolyHedron()
>> 48 //
>> 49 // 31 October 2004, P R Truscott, QinetiQ Ltd, UK
>> 50 // - Created.
>> 51 //
>> 52 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
47 53
48 #include "G4TessellatedSolid.hh" 54 #include "G4TessellatedSolid.hh"
>> 55 #include "G4PolyhedronArbitrary.hh"
>> 56 #include "globals.hh"
49 57
50 #if !defined(G4GEOM_USE_UTESSELLATEDSOLID) <<
51 <<
52 #include <algorithm> <<
53 #include <fstream> <<
54 #include <iomanip> <<
55 #include <iostream> 58 #include <iostream>
56 #include <list> <<
57 #include <random> <<
58 #include <stack> <<
59 <<
60 #include "geomdefs.hh" <<
61 #include "Randomize.hh" <<
62 #include "G4SystemOfUnits.hh" <<
63 #include "G4PhysicalConstants.hh" <<
64 #include "G4GeometryTolerance.hh" <<
65 #include "G4VoxelLimits.hh" <<
66 #include "G4AffineTransform.hh" <<
67 #include "G4BoundingEnvelope.hh" <<
68 <<
69 #include "G4VGraphicsScene.hh" <<
70 #include "G4VisExtent.hh" <<
71 <<
72 #include "G4AutoLock.hh" <<
73 <<
74 namespace <<
75 { <<
76 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE <<
77 } <<
78 <<
79 using namespace std; <<
80 59
81 ////////////////////////////////////////////// 60 ///////////////////////////////////////////////////////////////////////////////
82 // 61 //
83 // Standard contructor has blank name and defi << 62 // Standard contructor has blank name and defines no facets.
84 // 63 //
85 G4TessellatedSolid::G4TessellatedSolid () : G4 << 64 G4TessellatedSolid::G4TessellatedSolid ()
>> 65 : G4VSolid("dummy"), fpPolyhedron(0), cubicVolume(0.)
86 { 66 {
87 Initialize(); << 67 dirTolerance = 1.0E-14;
>> 68
>> 69 geometryType = "G4TessellatedSolid";
>> 70 facets.clear();
>> 71 solidClosed = false;
>> 72
>> 73 xMinExtent = kInfinity;
>> 74 xMaxExtent = -kInfinity;
>> 75 yMinExtent = kInfinity;
>> 76 yMaxExtent = -kInfinity;
>> 77 zMinExtent = kInfinity;
>> 78 zMaxExtent = -kInfinity;
88 } 79 }
89 80
90 ////////////////////////////////////////////// 81 ///////////////////////////////////////////////////////////////////////////////
91 // 82 //
92 // Alternative constructor. Simple define name << 83 // Alternative constructor. Simple define name and geometry type - no facets
93 // to detine. 84 // to detine.
94 // 85 //
95 G4TessellatedSolid::G4TessellatedSolid (const << 86 G4TessellatedSolid::G4TessellatedSolid (const G4String &name)
96 : G4VSolid(name) << 87 : G4VSolid(name), fpPolyhedron(0), cubicVolume(0.)
97 { 88 {
98 Initialize(); << 89 dirTolerance = 1.0E-14;
>> 90
>> 91 geometryType = "G4TessellatedSolid";
>> 92 facets.clear();
>> 93 solidClosed = false;
>> 94
>> 95 xMinExtent = kInfinity;
>> 96 xMaxExtent = -kInfinity;
>> 97 yMinExtent = kInfinity;
>> 98 yMaxExtent = -kInfinity;
>> 99 zMinExtent = kInfinity;
>> 100 zMaxExtent = -kInfinity;
99 } 101 }
100 102
101 ////////////////////////////////////////////// 103 ///////////////////////////////////////////////////////////////////////////////
102 // 104 //
103 // Fake default constructor - sets only member << 105 // Destructor.
104 // for usage restri <<
105 // 106 //
106 G4TessellatedSolid::G4TessellatedSolid( __void << 107 G4TessellatedSolid::~G4TessellatedSolid ()
107 { <<
108 Initialize(); <<
109 fMinExtent.set(0,0,0); <<
110 fMaxExtent.set(0,0,0); <<
111 } <<
112 <<
113 <<
114 ////////////////////////////////////////////// <<
115 G4TessellatedSolid::~G4TessellatedSolid() <<
116 { 108 {
117 DeleteObjects(); << 109 DeleteObjects ();
118 } <<
119 <<
120 ////////////////////////////////////////////// <<
121 // <<
122 // Copy constructor. <<
123 // <<
124 G4TessellatedSolid::G4TessellatedSolid (const <<
125 : G4VSolid(ts) <<
126 { <<
127 Initialize(); <<
128 <<
129 CopyObjects(ts); <<
130 } 110 }
131 111
132 ////////////////////////////////////////////// 112 ///////////////////////////////////////////////////////////////////////////////
133 // 113 //
134 // Assignment operator. << 114 // Define assignment operator.
135 // 115 //
136 G4TessellatedSolid& << 116 const G4TessellatedSolid &G4TessellatedSolid::operator=
137 G4TessellatedSolid::operator= (const G4Tessell << 117 (const G4TessellatedSolid &s)
138 { 118 {
139 if (&ts == this) return *this; << 119 if (&s == this) return *this;
140 << 120
141 // Copy base class data <<
142 G4VSolid::operator=(ts); <<
143 <<
144 DeleteObjects (); 121 DeleteObjects ();
145 << 122 CopyObjects (s);
146 Initialize(); << 123
147 <<
148 CopyObjects (ts); <<
149 <<
150 return *this; 124 return *this;
151 } 125 }
152 126
153 ////////////////////////////////////////////// 127 ///////////////////////////////////////////////////////////////////////////////
154 // 128 //
155 void G4TessellatedSolid::Initialize() << 129 void G4TessellatedSolid::DeleteObjects ()
156 { 130 {
157 kCarToleranceHalf = 0.5*kCarTolerance; << 131 for (std::vector<G4VFacet *>::iterator f=facets.begin();
158 << 132 f!=facets.end(); f++) delete *f;
159 fRebuildPolyhedron = false; fpPolyhedron = n << 133 facets.clear();
160 fCubicVolume = 0.; fSurfaceArea = 0.; <<
161 <<
162 fGeometryType = "G4TessellatedSolid"; <<
163 fSolidClosed = false; <<
164 <<
165 fMinExtent.set(kInfinity,kInfinity,kInfinity <<
166 fMaxExtent.set(-kInfinity,-kInfinity,-kInfin <<
167 <<
168 SetRandomVectors(); <<
169 } 134 }
170 135
171 ////////////////////////////////////////////// 136 ///////////////////////////////////////////////////////////////////////////////
172 // 137 //
173 void G4TessellatedSolid::DeleteObjects() << 138 void G4TessellatedSolid::CopyObjects (const G4TessellatedSolid &s)
174 { 139 {
175 std::size_t size = fFacets.size(); << 140 size_t n = s.GetNumberOfFacets();
176 for (std::size_t i = 0; i < size; ++i) { de << 141 for (size_t i=0; i<n; n++)
177 fFacets.clear(); <<
178 delete fpPolyhedron; fpPolyhedron = nullptr; <<
179 } <<
180 <<
181 ////////////////////////////////////////////// <<
182 // <<
183 void G4TessellatedSolid::CopyObjects (const G4 <<
184 { <<
185 G4ThreeVector reductionRatio; <<
186 G4int fmaxVoxels = fVoxels.GetMaxVoxels(redu <<
187 if (fmaxVoxels < 0) <<
188 fVoxels.SetMaxVoxels(reductionRatio); <<
189 else <<
190 fVoxels.SetMaxVoxels(fmaxVoxels); <<
191 <<
192 G4int n = ts.GetNumberOfFacets(); <<
193 for (G4int i = 0; i < n; ++i) <<
194 { 142 {
195 G4VFacet *facetClone = (ts.GetFacet(i))->G << 143 G4VFacet *facetClone = (s.GetFacet(i))->GetClone();
196 AddFacet(facetClone); 144 AddFacet(facetClone);
197 } 145 }
198 if (ts.GetSolidClosed()) SetSolidClosed(true << 146 solidClosed = s.GetSolidClosed();
>> 147
>> 148 // cubicVolume = s.GetCubicVolume();
199 } 149 }
200 150
201 ////////////////////////////////////////////// 151 ///////////////////////////////////////////////////////////////////////////////
202 // 152 //
203 // Add a facet to the facet list. << 153 // Add a facet to the facet list. Note that you can add, but you cannot
204 // Note that you can add, but you cannot delet << 154 // delete.
205 // 155 //
206 G4bool G4TessellatedSolid::AddFacet (G4VFacet* << 156 G4bool G4TessellatedSolid::AddFacet (G4VFacet *aFacet)
207 { 157 {
208 // Add the facet to the vector. 158 // Add the facet to the vector.
209 // << 159
210 if (fSolidClosed) << 160 if (solidClosed)
211 { 161 {
212 G4Exception("G4TessellatedSolid::AddFacet( << 162 G4Exception("G4TessellatedSolid::AddFacet()", "InvalidSetup",
213 JustWarning, "Attempt to add f 163 JustWarning, "Attempt to add facets when solid is closed.");
214 return false; 164 return false;
215 } 165 }
216 else if (aFacet->IsDefined()) 166 else if (aFacet->IsDefined())
217 { 167 {
218 set<G4VertexInfo,G4VertexComparator>::iter << 168 if (facets.size() == 0)
219 = fFacetList.begin(), end = fFacetList.e <<
220 G4ThreeVector p = aFacet->GetCircumcentre( <<
221 G4VertexInfo value; <<
222 value.id = (G4int)fFacetList.size(); <<
223 value.mag2 = p.x() + p.y() + p.z(); <<
224 <<
225 G4bool found = false; <<
226 if (!OutsideOfExtent(p, kCarTolerance)) <<
227 { 169 {
228 G4double kCarTolerance3 = 3 * kCarTolera << 170 facets.push_back(aFacet);
229 pos = fFacetList.lower_bound(value); << 171 }
230 << 172 else
231 it = pos; << 173 {
232 while (!found && it != end) // Loop c << 174 G4bool found = false;
>> 175 FacetI it = facets.begin();
>> 176 do
>> 177 {
>> 178 found = (**it == *aFacet);
>> 179 } while (!found && ++it!=facets.end());
>> 180
>> 181 if (found)
233 { 182 {
234 G4int id = (*it).id; << 183 delete *it;
235 G4VFacet *facet = fFacets[id]; << 184 facets.erase(it);
236 G4ThreeVector q = facet->GetCircumcent <<
237 if ((found = (facet == aFacet))) break <<
238 G4double dif = q.x() + q.y() + q.z() - <<
239 if (dif > kCarTolerance3) break; <<
240 it++; <<
241 } 185 }
242 << 186 else
243 if (fFacets.size() > 1) <<
244 { 187 {
245 it = pos; << 188 facets.push_back(aFacet);
246 while (!found && it != begin) // Lo <<
247 { <<
248 --it; <<
249 G4int id = (*it).id; <<
250 G4VFacet *facet = fFacets[id]; <<
251 G4ThreeVector q = facet->GetCircumce <<
252 found = (facet == aFacet); <<
253 if (found) break; <<
254 G4double dif = value.mag2 - (q.x() + <<
255 if (dif > kCarTolerance3) break; <<
256 } <<
257 } 189 }
258 } 190 }
259 << 191
260 if (!found) <<
261 { <<
262 fFacets.push_back(aFacet); <<
263 fFacetList.insert(value); <<
264 } <<
265 return true; 192 return true;
266 } 193 }
267 else 194 else
268 { 195 {
269 G4Exception("G4TessellatedSolid::AddFacet( << 196 G4Exception("G4TessellatedSolid::AddFacet()", "InvalidSetup",
270 JustWarning, "Attempt to add f 197 JustWarning, "Attempt to add facet not properly defined.");
271 aFacet->StreamInfo(G4cout); << 198 G4cerr << "Facet attributes:" << G4endl;
>> 199 aFacet->StreamInfo(G4cerr);
>> 200 G4cerr << G4endl;
>> 201
272 return false; 202 return false;
273 } 203 }
274 } 204 }
275 205
276 ////////////////////////////////////////////// 206 ///////////////////////////////////////////////////////////////////////////////
277 // 207 //
278 G4int G4TessellatedSolid::SetAllUsingStack(con << 208 void G4TessellatedSolid::SetSolidClosed (const G4bool t)
279 con <<
280 G4b <<
281 { <<
282 vector<G4int> xyz = voxel; <<
283 stack<vector<G4int> > pos; <<
284 pos.push(xyz); <<
285 G4int filled = 0; <<
286 <<
287 vector<G4int> candidates; <<
288 <<
289 while (!pos.empty()) // Loop checking, 13 <<
290 { <<
291 xyz = pos.top(); <<
292 pos.pop(); <<
293 G4int index = fVoxels.GetVoxelsIndex(xyz); <<
294 if (!checked[index]) <<
295 { <<
296 checked.SetBitNumber(index, true); <<
297 <<
298 if (fVoxels.IsEmpty(index)) <<
299 { <<
300 ++filled; <<
301 <<
302 fInsides.SetBitNumber(index, status); <<
303 <<
304 for (auto i = 0; i <= 2; ++i) <<
305 { <<
306 if (xyz[i] < max[i] - 1) <<
307 { <<
308 xyz[i]++; <<
309 pos.push(xyz); <<
310 xyz[i]--; <<
311 } <<
312 <<
313 if (xyz[i] > 0) <<
314 { <<
315 xyz[i]--; <<
316 pos.push(xyz); <<
317 xyz[i]++; <<
318 } <<
319 } <<
320 } <<
321 } <<
322 } <<
323 return filled; <<
324 } <<
325 <<
326 ////////////////////////////////////////////// <<
327 // <<
328 void G4TessellatedSolid::PrecalculateInsides() <<
329 { <<
330 vector<G4int> voxel(3), maxVoxels(3); <<
331 for (auto i = 0; i <= 2; ++i) <<
332 maxVoxels[i] = (G4int)fVoxels.GetBoundary( <<
333 G4int size = maxVoxels[0] * maxVoxels[1] * m <<
334 <<
335 G4SurfBits checked(size-1); <<
336 fInsides.Clear(); <<
337 fInsides.ResetBitNumber(size-1); <<
338 <<
339 G4ThreeVector point; <<
340 for (voxel[2] = 0; voxel[2] < maxVoxels[2] - <<
341 { <<
342 for (voxel[1] = 0; voxel[1] < maxVoxels[1] <<
343 { <<
344 for (voxel[0] = 0; voxel[0] < maxVoxels[ <<
345 { <<
346 G4int index = fVoxels.GetVoxelsIndex(v <<
347 if (!checked[index] && fVoxels.IsEmpty <<
348 { <<
349 for (auto i = 0; i <= 2; ++i) <<
350 { <<
351 point[i] = fVoxels.GetBoundary(i)[ <<
352 } <<
353 auto inside = (G4bool) (InsideNoVoxe <<
354 SetAllUsingStack(voxel, maxVoxels, i <<
355 } <<
356 else checked.SetBitNumber(index); <<
357 } <<
358 } <<
359 } <<
360 } <<
361 <<
362 ////////////////////////////////////////////// <<
363 // <<
364 void G4TessellatedSolid::Voxelize () <<
365 { <<
366 #ifdef G4SPECSDEBUG <<
367 G4cout << "Voxelizing..." << G4endl; <<
368 #endif <<
369 fVoxels.Voxelize(fFacets); <<
370 <<
371 if (fVoxels.Empty().GetNbits() != 0u) <<
372 { <<
373 #ifdef G4SPECSDEBUG <<
374 G4cout << "Precalculating Insides..." << G <<
375 #endif <<
376 PrecalculateInsides(); <<
377 } <<
378 } <<
379 <<
380 ////////////////////////////////////////////// <<
381 // <<
382 // Compute extremeFacets, i.e. find those face <<
383 // planes that bound the volume. <<
384 // Note that this is going to reject concaved <<
385 // note that if the vertex is on the facet, di <<
386 // returns true. So will this work?? Need non <<
387 // "G4bool inside = displacement < 0.0;" <<
388 // or <<
389 // "G4bool inside = displacement <= -0.5*kCarT <<
390 // (Notes from PT 13/08/2007). <<
391 // <<
392 void G4TessellatedSolid::SetExtremeFacets() <<
393 { <<
394 // Copy vertices to local array <<
395 std::size_t vsize = fVertexList.size(); <<
396 std::vector<G4ThreeVector> vertices(vsize); <<
397 for (std::size_t i = 0; i < vsize; ++i) { ve <<
398 <<
399 // Shuffle vertices <<
400 std::mt19937 gen(12345678); <<
401 std::shuffle(vertices.begin(), vertices.end( <<
402 <<
403 // Select six extreme vertices in different <<
404 G4ThreeVector points[6]; <<
405 for (auto & point : points) { point = vertic <<
406 for (std::size_t i=1; i < vsize; ++i) <<
407 { <<
408 if (vertices[i].x() < points[0].x()) point <<
409 if (vertices[i].x() > points[1].x()) point <<
410 if (vertices[i].y() < points[2].y()) point <<
411 if (vertices[i].y() > points[3].y()) point <<
412 if (vertices[i].z() < points[4].z()) point <<
413 if (vertices[i].z() > points[5].z()) point <<
414 } <<
415 <<
416 // Find extreme facets <<
417 std::size_t size = fFacets.size(); <<
418 for (std::size_t j = 0; j < size; ++j) <<
419 { <<
420 G4VFacet &facet = *fFacets[j]; <<
421 <<
422 // Check extreme vertices <<
423 if (!facet.IsInside(points[0])) continue; <<
424 if (!facet.IsInside(points[1])) continue; <<
425 if (!facet.IsInside(points[2])) continue; <<
426 if (!facet.IsInside(points[3])) continue; <<
427 if (!facet.IsInside(points[4])) continue; <<
428 if (!facet.IsInside(points[5])) continue; <<
429 <<
430 // Check vertices <<
431 G4bool isExtreme = true; <<
432 for (std::size_t i=0; i < vsize; ++i) <<
433 { <<
434 if (!facet.IsInside(vertices[i])) <<
435 { <<
436 isExtreme = false; <<
437 break; <<
438 } <<
439 } <<
440 if (isExtreme) fExtremeFacets.insert(&face <<
441 } <<
442 } <<
443 <<
444 ////////////////////////////////////////////// <<
445 // <<
446 void G4TessellatedSolid::CreateVertexList() <<
447 { 209 {
448 // The algorithm: << 210 if (t)
449 // we will have additional vertexListSorted, <<
450 // sorted by magnitude of vertice vector. <<
451 // New candidate for fVertexList - we will d <<
452 // item which would be within its magnitude <<
453 // We will go trough until we will reach > + <<
454 // Comparison (q-p).mag() < 0.5*kCarToleranc <<
455 // They can be just stored in std::vector, w <<
456 // on binary search. <<
457 <<
458 set<G4VertexInfo,G4VertexComparator> vertexL <<
459 set<G4VertexInfo,G4VertexComparator>::iterat <<
460 = vertexListSorted.begin(), end = vertexL <<
461 G4ThreeVector p; <<
462 G4VertexInfo value; <<
463 <<
464 fVertexList.clear(); <<
465 std::size_t size = fFacets.size(); <<
466 <<
467 G4double kCarTolerance24 = kCarTolerance * k <<
468 G4double kCarTolerance3 = 3 * kCarTolerance; <<
469 vector<G4int> newIndex(100); <<
470 <<
471 for (std::size_t k = 0; k < size; ++k) <<
472 { 211 {
473 G4VFacet &facet = *fFacets[k]; << 212 vertexList.clear();
474 G4int max = facet.GetNumberOfVertices(); << 213 for (FacetCI it=facets.begin(); it!=facets.end(); it++)
475 <<
476 for (G4int i = 0; i < max; ++i) <<
477 { 214 {
478 p = facet.GetVertex(i); << 215 size_t m = vertexList.size();
479 value.id = (G4int)fVertexList.size(); << 216 G4ThreeVector p(0.0,0.0,0.0);
480 value.mag2 = p.x() + p.y() + p.z(); << 217 for (size_t i=0; i<(*it)->GetNumberOfVertices(); i++)
481 << 218 {
482 G4bool found = false; << 219 p = (*it)->GetVertex(i);
483 G4int id = 0; << 220 G4bool found = false;
484 if (!OutsideOfExtent(p, kCarTolerance)) << 221 size_t j = 0;
485 { << 222 while (j < m && !found)
486 pos = vertexListSorted.lower_bound(val << 223 {
487 it = pos; << 224 G4ThreeVector q = vertexList[j];
488 while (it != end) // Loop checking, << 225 found = (q-p).mag() < 0.5*kCarTolerance;
489 { << 226 if (!found) j++;
490 id = (*it).id; <<
491 G4ThreeVector q = fVertexList[id]; <<
492 G4double dif = (q-p).mag2(); <<
493 found = (dif < kCarTolerance24); <<
494 if (found) break; <<
495 dif = q.x() + q.y() + q.z() - value. <<
496 if (dif > kCarTolerance3) break; <<
497 ++it; <<
498 } 227 }
499 << 228
500 if (!found && (fVertexList.size() > 1) << 229 if (!found)
501 { 230 {
502 it = pos; << 231 vertexList.push_back(p);
503 while (it != begin) // Loop check << 232 (*it)->SetVertexIndex(i,vertexList.size()-1);
504 { <<
505 --it; <<
506 id = (*it).id; <<
507 G4ThreeVector q = fVertexList[id]; <<
508 G4double dif = (q-p).mag2(); <<
509 found = (dif < kCarTolerance24); <<
510 if (found) break; <<
511 dif = value.mag2 - (q.x() + q.y() <<
512 if (dif > kCarTolerance3) break; <<
513 } <<
514 } 233 }
515 } <<
516 <<
517 if (!found) <<
518 { <<
519 #ifdef G4SPECSDEBUG <<
520 G4cout << p.x() << ":" << p.y() << ":" <<
521 G4cout << "Adding new vertex #" << i < <<
522 << " id " << value.id << G4endl <<
523 G4cout << "===" << G4endl; <<
524 #endif <<
525 fVertexList.push_back(p); <<
526 vertexListSorted.insert(value); <<
527 begin = vertexListSorted.begin(); <<
528 end = vertexListSorted.end(); <<
529 newIndex[i] = value.id; <<
530 // <<
531 // Now update the maximum x, y and z l <<
532 // <<
533 if (value.id == 0) fMinExtent = fMaxEx <<
534 else 234 else
535 { 235 {
536 if (p.x() > fMaxExtent.x()) fMaxExte << 236 (*it)->SetVertexIndex(i,j);
537 else if (p.x() < fMinExtent.x()) fMi <<
538 if (p.y() > fMaxExtent.y()) fMaxExte <<
539 else if (p.y() < fMinExtent.y()) fMi <<
540 if (p.z() > fMaxExtent.z()) fMaxExte <<
541 else if (p.z() < fMinExtent.z()) fMi <<
542 } 237 }
543 } 238 }
544 else <<
545 { <<
546 #ifdef G4SPECSDEBUG <<
547 G4cout << p.x() << ":" << p.y() << ":" <<
548 G4cout << "Vertex #" << i << " of face <<
549 << " found, redirecting to " << <<
550 G4cout << "===" << G4endl; <<
551 #endif <<
552 newIndex[i] = id; <<
553 } <<
554 } 239 }
555 // only now it is possible to change verti <<
556 // 240 //
557 facet.SetVertices(&fVertexList); << 241 // Now update the maximum x, y and z limits of the volume.
558 for (G4int i = 0; i < max; ++i) << 242 //
559 facet.SetVertexIndex(i,newIndex[i]); << 243 for (size_t i=0; i<vertexList.size(); i++)
560 } <<
561 vector<G4ThreeVector>(fVertexList).swap(fVer <<
562 <<
563 #ifdef G4SPECSDEBUG <<
564 G4double previousValue = 0.; <<
565 for (auto res=vertexListSorted.cbegin(); res <<
566 { <<
567 G4int id = (*res).id; <<
568 G4ThreeVector vec = fVertexList[id]; <<
569 G4double mvalue = vec.x() + vec.y() + vec. <<
570 if (previousValue && (previousValue - 1e-9 <<
571 G4cout << "Error in CreateVertexList: pr <<
572 << " is smaller than mvalue " << <<
573 previousValue = mvalue; <<
574 } <<
575 #endif <<
576 } <<
577 <<
578 ////////////////////////////////////////////// <<
579 // <<
580 void G4TessellatedSolid::DisplayAllocatedMemor <<
581 { <<
582 G4int without = AllocatedMemoryWithoutVoxels <<
583 G4int with = AllocatedMemory(); <<
584 G4double ratio = (G4double) with / without; <<
585 G4cout << "G4TessellatedSolid - Allocated me <<
586 << without << "; with " << with << "; <<
587 } <<
588 <<
589 ////////////////////////////////////////////// <<
590 // <<
591 void G4TessellatedSolid::SetSolidClosed (const <<
592 { <<
593 if (t) <<
594 { <<
595 #ifdef G4SPECSDEBUG <<
596 G4cout << "Creating vertex list..." << G4e <<
597 #endif <<
598 CreateVertexList(); <<
599 <<
600 #ifdef G4SPECSDEBUG <<
601 G4cout << "Setting extreme facets..." << G <<
602 #endif <<
603 SetExtremeFacets(); <<
604 <<
605 #ifdef G4SPECSDEBUG <<
606 G4cout << "Voxelizing..." << G4endl; <<
607 #endif <<
608 Voxelize(); <<
609 <<
610 #ifdef G4SPECSDEBUG <<
611 DisplayAllocatedMemory(); <<
612 #endif <<
613 <<
614 #ifdef G4SPECSDEBUG <<
615 G4cout << "Checking Structure..." << G4end <<
616 #endif <<
617 G4int irep = CheckStructure(); <<
618 if (irep != 0) <<
619 { 244 {
620 if ((irep & 1) != 0) << 245 G4ThreeVector p = vertexList[i];
621 { << 246 G4double x = p.x();
622 std::ostringstream message; << 247 G4double y = p.y();
623 message << "Defects in solid: " << Ge << 248 G4double z = p.z();
624 << " - negative cubic volume, << 249
625 G4Exception("G4TessellatedSolid::SetS << 250 if (i > 0)
626 "GeomSolids1001", JustWar << 251 {
627 } << 252 if (x > xMaxExtent) xMaxExtent = x;
628 if ((irep & 2) != 0) << 253 if (x < xMinExtent) xMinExtent = x;
629 { << 254 if (y > yMaxExtent) yMaxExtent = y;
630 std::ostringstream message; << 255 if (y < yMinExtent) yMinExtent = y;
631 message << "Defects in solid: " << Ge << 256 if (z > zMaxExtent) zMaxExtent = z;
632 << " - some facets have wrong << 257 if (z < zMinExtent) zMinExtent = z;
633 G4Exception("G4TessellatedSolid::SetS <<
634 "GeomSolids1001", JustWar <<
635 } 258 }
636 if ((irep & 4) != 0) << 259 else
637 { 260 {
638 std::ostringstream message; << 261 xMaxExtent = x;
639 message << "Defects in solid: " << Ge << 262 xMinExtent = x;
640 << " - there are holes in the << 263 yMaxExtent = y;
641 G4Exception("G4TessellatedSolid::SetS << 264 yMinExtent = y;
642 "GeomSolids1001", JustWar << 265 zMaxExtent = z;
>> 266 zMinExtent = z;
643 } 267 }
644 } 268 }
>> 269 solidClosed = true;
>> 270 }
>> 271 else
>> 272 {
>> 273 solidClosed = false;
645 } 274 }
646 fSolidClosed = t; <<
647 } 275 }
648 276
649 ////////////////////////////////////////////// 277 ///////////////////////////////////////////////////////////////////////////////
650 // 278 //
651 // GetSolidClosed <<
652 // <<
653 // Used to determine whether the solid is clos <<
654 // <<
655 G4bool G4TessellatedSolid::GetSolidClosed () c 279 G4bool G4TessellatedSolid::GetSolidClosed () const
656 { << 280 {return solidClosed;}
657 return fSolidClosed; <<
658 } <<
659 281
660 ////////////////////////////////////////////// 282 ///////////////////////////////////////////////////////////////////////////////
661 // 283 //
662 // CheckStructure << 284 const G4TessellatedSolid &G4TessellatedSolid::operator+=
663 // << 285 (const G4TessellatedSolid &right)
664 // Checks structure of the solid. Return value <<
665 // defect indicators, if any (0 means no defec <<
666 // 1 - cubic volume is negative, wrong orien <<
667 // 2 - some facets have wrong orientation <<
668 // 4 - holes in the surface <<
669 // <<
670 G4int G4TessellatedSolid::CheckStructure() con <<
671 { 286 {
672 G4int nedge = 0; << 287 for (size_t i=0; i<right.GetNumberOfFacets(); i++)
673 std::size_t nface = fFacets.size(); <<
674 <<
675 // Calculate volume <<
676 // <<
677 G4double volume = 0.; <<
678 for (std::size_t i = 0; i < nface; ++i) <<
679 { <<
680 G4VFacet& facet = *fFacets[i]; <<
681 nedge += facet.GetNumberOfVertices(); <<
682 volume += facet.GetArea()*(facet.GetVertex <<
683 } <<
684 auto ivolume = static_cast<G4int>(volume <= <<
685 <<
686 // Create sorted vector of edges <<
687 // <<
688 std::vector<int64_t> iedge(nedge); <<
689 G4int kk = 0; <<
690 for (std::size_t i = 0; i < nface; ++i) <<
691 { <<
692 G4VFacet& facet = *fFacets[i]; <<
693 G4int nnode = facet.GetNumberOfVertices(); <<
694 for (G4int k = 0; k < nnode; ++k) <<
695 { <<
696 int64_t i1 = facet.GetVertexIndex((k == <<
697 int64_t i2 = facet.GetVertexIndex(k); <<
698 auto inverse = static_cast<int64_t>(i2 <<
699 if (inverse != 0) std::swap(i1, i2); <<
700 iedge[kk++] = i1*1000000000 + i2*2 + inv <<
701 } <<
702 } <<
703 std::sort(iedge.begin(), iedge.end()); <<
704 <<
705 // Check edges, correct structure should con <<
706 // with different orientation <<
707 // <<
708 G4int iorder = 0; <<
709 G4int ihole = 0; <<
710 G4int i = 0; <<
711 while (i < nedge - 1) <<
712 { <<
713 if (iedge[i + 1] - iedge[i] == 1) // paire <<
714 { <<
715 i += 2; <<
716 } <<
717 else if (iedge[i + 1] == iedge[i]) // pair <<
718 { <<
719 iorder = 2; <<
720 i += 2; <<
721 } <<
722 else // unpaired edge <<
723 { <<
724 ihole = 4; <<
725 i++; <<
726 } <<
727 } <<
728 return ivolume + iorder + ihole; <<
729 } <<
730 <<
731 ////////////////////////////////////////////// <<
732 // <<
733 // operator+= <<
734 // <<
735 // This operator allows the user to add two te <<
736 // that the solid on the left then includes al <<
737 // on the right. Note that copies of the face <<
738 // using the original facet set of the solid o <<
739 // <<
740 G4TessellatedSolid& <<
741 G4TessellatedSolid::operator+=(const G4Tessell <<
742 { <<
743 G4int size = right.GetNumberOfFacets(); <<
744 for (G4int i = 0; i < size; ++i) <<
745 AddFacet(right.GetFacet(i)->GetClone()); 288 AddFacet(right.GetFacet(i)->GetClone());
746 <<
747 return *this; 289 return *this;
748 } 290 }
749 291
750 ////////////////////////////////////////////// 292 ///////////////////////////////////////////////////////////////////////////////
751 // 293 //
752 // GetNumberOfFacets << 294 G4VFacet *G4TessellatedSolid::GetFacet (size_t i) const
753 // <<
754 G4int G4TessellatedSolid::GetNumberOfFacets() <<
755 { 295 {
756 return (G4int)fFacets.size(); << 296 return facets[i];
757 } 297 }
758 298
759 ////////////////////////////////////////////// 299 ///////////////////////////////////////////////////////////////////////////////
760 // 300 //
761 EInside G4TessellatedSolid::InsideVoxels(const << 301 size_t G4TessellatedSolid::GetNumberOfFacets () const
762 { 302 {
763 // << 303 return facets.size();
764 // First the simple test - check if we're ou <<
765 // of the tessellated solid. <<
766 // <<
767 if (OutsideOfExtent(p, kCarTolerance)) <<
768 return kOutside; <<
769 <<
770 vector<G4int> startingVoxel(3); <<
771 fVoxels.GetVoxel(startingVoxel, p); <<
772 <<
773 const G4double dirTolerance = 1.0E-14; <<
774 <<
775 const vector<G4int> &startingCandidates = <<
776 fVoxels.GetCandidates(startingVoxel); <<
777 std::size_t limit = startingCandidates.size( <<
778 if (limit == 0 && (fInsides.GetNbits() != 0u <<
779 { <<
780 G4int index = fVoxels.GetPointIndex(p); <<
781 EInside location = fInsides[index] ? kInsi <<
782 return location; <<
783 } <<
784 <<
785 G4double minDist = kInfinity; <<
786 <<
787 for(std::size_t i = 0; i < limit; ++i) <<
788 { <<
789 G4int candidate = startingCandidates[i]; <<
790 G4VFacet &facet = *fFacets[candidate]; <<
791 G4double dist = facet.Distance(p,minDist); <<
792 if (dist < minDist) minDist = dist; <<
793 if (dist <= kCarToleranceHalf) <<
794 return kSurface; <<
795 } <<
796 <<
797 // The following is something of an adaptati <<
798 // Rickard Holmberg augmented with informati <<
799 // "Geometric Tools for Computer Graphics," <<
800 // we're trying to determine whether we're i <<
801 // a few rays and determining if the first s <<
802 // vector between 0 to pi/2 (out-going) or p <<
803 // We should also avoid rays which are nearl <<
804 // tessellated surface, and therefore produc <<
805 // For the moment, this is a bit over-engine <<
806 // <<
807 G4double distOut = kInfinity; <<
808 G4double distIn = kInfinity; <<
809 G4double distO = 0.0; <<
810 G4double distI = 0.0; <<
811 G4double distFromSurfaceO = 0.0; <<
812 G4double distFromSurfaceI = 0.0; <<
813 G4ThreeVector normalO, normalI; <<
814 G4bool crossingO = false; <<
815 G4bool crossingI = false; <<
816 EInside location = kOutside; <<
817 G4int sm = 0; <<
818 <<
819 G4bool nearParallel = false; <<
820 do // Loop checking, 13.08.2015, G.Cosmo <<
821 { <<
822 // We loop until we find direction where t <<
823 // to the surface of any facet since this <<
824 // case is that the angles should be suffi <<
825 // are 20 random directions to select from <<
826 // <<
827 distOut = distIn = kInfinity; <<
828 const G4ThreeVector& v = fRandir[sm]; <<
829 ++sm; <<
830 // <<
831 // This code could be voxelized by the sam <<
832 // DistanceToOut(). We will traverse throu <<
833 // intersect only for those, which would b <<
834 // checked before. <<
835 // <<
836 G4ThreeVector currentPoint = p; <<
837 G4ThreeVector direction = v.unit(); <<
838 // G4SurfBits exclusion(fVoxels.GetBitsPer <<
839 vector<G4int> curVoxel(3); <<
840 curVoxel = startingVoxel; <<
841 G4double shiftBonus = kCarTolerance; <<
842 <<
843 G4bool crossed = false; <<
844 G4bool started = true; <<
845 <<
846 do // Loop checking, 13.08.2015, G.Cosm <<
847 { <<
848 const vector<G4int> &candidates = <<
849 started ? startingCandidates : fVoxels <<
850 started = false; <<
851 if (auto candidatesCount = (G4int)candid <<
852 { <<
853 for (G4int i = 0 ; i < candidatesCount <<
854 { <<
855 G4int candidate = candidates[i]; <<
856 // bits.SetBitNumber(candidate); <<
857 G4VFacet& facet = *fFacets[candidate <<
858 <<
859 crossingO = facet.Intersect(p,v,true <<
860 crossingI = facet.Intersect(p,v,fals <<
861 <<
862 if (crossingO || crossingI) <<
863 { <<
864 crossed = true; <<
865 <<
866 nearParallel = (crossingO <<
867 && std::fabs(normalO.dot( <<
868 || (crossingI && std::fab <<
869 if (!nearParallel) <<
870 { <<
871 if (crossingO && distO > 0.0 && <<
872 distOut = distO; <<
873 if (crossingI && distI > 0.0 && <<
874 distIn = distI; <<
875 } <<
876 else break; <<
877 } <<
878 } <<
879 if (nearParallel) break; <<
880 } <<
881 else <<
882 { <<
883 if (!crossed) <<
884 { <<
885 G4int index = fVoxels.GetVoxelsIndex <<
886 G4bool inside = fInsides[index]; <<
887 location = inside ? kInside : kOutsi <<
888 return location; <<
889 } <<
890 } <<
891 <<
892 G4double shift=fVoxels.DistanceToNext(cu <<
893 if (shift == kInfinity) break; <<
894 <<
895 currentPoint += direction * (shift + shi <<
896 } <<
897 while (fVoxels.UpdateCurrentVoxel(currentP <<
898 <<
899 } <<
900 while (nearParallel && sm != fMaxTries); <<
901 // <<
902 // Here we loop through the facets to find o <<
903 // between the ray and that facet. The test <<
904 // the ray is entering the facet or exiting. <<
905 // <<
906 #ifdef G4VERBOSE <<
907 if (sm == fMaxTries) <<
908 { <<
909 // <<
910 // We've run out of random vector directio <<
911 // low (nTries <= 0.5*maxTries) then this <<
912 // something wrong with geometry. <<
913 // <<
914 std::ostringstream message; <<
915 G4long oldprc = message.precision(16); <<
916 message << "Cannot determine whether point <<
917 << G4endl <<
918 << "Solid name = " << GetName() < <<
919 << "Geometry Type = " << fGeometryTyp <<
920 << "Number of facets = " << fFacets.size <<
921 << "Position:" << G4endl << G4endl <<
922 << "p.x() = " << p.x()/mm << " mm" << <<
923 << "p.y() = " << p.y()/mm << " mm" << <<
924 << "p.z() = " << p.z()/mm << " mm"; <<
925 message.precision(oldprc); <<
926 G4Exception("G4TessellatedSolid::Inside()" <<
927 "GeomSolids1002", JustWarning, <<
928 } <<
929 #endif <<
930 <<
931 // In the next if-then-elseif G4String the l <<
932 // (1) You don't hit anything so cannot be i <<
933 // constructed correctly! <<
934 // (2) Distance to inside (ie. nearest facet <<
935 // shorter than distance to outside (nea <<
936 // facet) - on condition of safety dista <<
937 // (3) Distance to outside is shorter than d <<
938 // we're inside. <<
939 // <<
940 if (distIn == kInfinity && distOut == kInfin <<
941 location = kOutside; <<
942 else if (distIn <= distOut - kCarToleranceHa <<
943 location = kOutside; <<
944 else if (distOut <= distIn - kCarToleranceHa <<
945 location = kInside; <<
946 <<
947 return location; <<
948 } 304 }
949 305
950 ////////////////////////////////////////////// 306 ///////////////////////////////////////////////////////////////////////////////
951 // 307 //
952 EInside G4TessellatedSolid::InsideNoVoxels (co << 308 EInside G4TessellatedSolid::Inside (const G4ThreeVector &p) const
953 { <<
954 // <<
955 // First the simple test - check if we're ou <<
956 // of the tessellated solid. <<
957 // <<
958 if (OutsideOfExtent(p, kCarTolerance)) <<
959 return kOutside; <<
960 <<
961 const G4double dirTolerance = 1.0E-14; <<
962 <<
963 G4double minDist = kInfinity; <<
964 // <<
965 // Check if we are close to a surface <<
966 // <<
967 std::size_t size = fFacets.size(); <<
968 for (std::size_t i = 0; i < size; ++i) <<
969 { <<
970 G4VFacet& facet = *fFacets[i]; <<
971 G4double dist = facet.Distance(p,minDist); <<
972 if (dist < minDist) minDist = dist; <<
973 if (dist <= kCarToleranceHalf) <<
974 { <<
975 return kSurface; <<
976 } <<
977 } <<
978 // <<
979 // The following is something of an adaptati <<
980 // Rickard Holmberg augmented with informati <<
981 // "Geometric Tools for Computer Graphics," <<
982 // trying to determine whether we're inside <<
983 // rays and determining if the first surface <<
984 // between 0 to pi/2 (out-going) or pi/2 to <<
985 // avoid rays which are nearly within the pl <<
986 // and therefore produce rays randomly. For <<
987 // over-engineered (belt-braces-and-ducttape <<
988 // <<
989 #if G4SPECSDEBUG <<
990 G4int nTry = 7; <<
991 #else <<
992 G4int nTry = 3; <<
993 #endif <<
994 G4double distOut = kInfinity; <<
995 G4double distIn = kInfinity; <<
996 G4double distO = 0.0; <<
997 G4double distI = 0.0; <<
998 G4double distFromSurfaceO = 0.0; <<
999 G4double distFromSurfaceI = 0.0; <<
1000 G4ThreeVector normalO(0.0,0.0,0.0); <<
1001 G4ThreeVector normalI(0.0,0.0,0.0); <<
1002 G4bool crossingO = false; <<
1003 G4bool crossingI = false; <<
1004 EInside location = kOutside; <<
1005 EInside locationprime = kOutside; <<
1006 G4int sm = 0; <<
1007 <<
1008 for (G4int i=0; i<nTry; ++i) <<
1009 { <<
1010 G4bool nearParallel = false; <<
1011 do // Loop checking, 13.08.2015, G.Cos <<
1012 { <<
1013 // <<
1014 // We loop until we find direction wher <<
1015 // to the surface of any facet since th <<
1016 // case is that the angles should be su <<
1017 // are 20 random directions to select f <<
1018 // <<
1019 distOut = distIn = kInfinity; <<
1020 G4ThreeVector v = fRandir[sm]; <<
1021 sm++; <<
1022 auto f = fFacets.cbegin(); <<
1023 <<
1024 do // Loop checking, 13.08.2015, G.C <<
1025 { <<
1026 // <<
1027 // Here we loop through the facets to <<
1028 // intersection between the ray and t <<
1029 // separately whether the ray is ente <<
1030 // <<
1031 crossingO = ((*f)->Intersect(p,v,true <<
1032 crossingI = ((*f)->Intersect(p,v,fals <<
1033 if (crossingO || crossingI) <<
1034 { <<
1035 nearParallel = (crossingO && std::f <<
1036 || (crossingI && std::f <<
1037 if (!nearParallel) <<
1038 { <<
1039 if (crossingO && distO > 0.0 && d <<
1040 if (crossingI && distI > 0.0 && d <<
1041 } <<
1042 } <<
1043 } while (!nearParallel && ++f != fFacet <<
1044 } while (nearParallel && sm != fMaxTries) <<
1045 <<
1046 #ifdef G4VERBOSE <<
1047 if (sm == fMaxTries) <<
1048 { <<
1049 // <<
1050 // We've run out of random vector direc <<
1051 // sufficiently low (nTries <= 0.5*maxT <<
1052 // that there is something wrong with g <<
1053 // <<
1054 std::ostringstream message; <<
1055 G4long oldprc = message.precision(16); <<
1056 message << "Cannot determine whether po <<
1057 << G4endl <<
1058 << "Solid name = " << GetName() <<
1059 << "Geometry Type = " << fGeometry <<
1060 << "Number of facets = " << fFacets.s <<
1061 << "Position:" << G4endl << G4endl <<
1062 << "p.x() = " << p.x()/mm << " mm" <<
1063 << "p.y() = " << p.y()/mm << " mm" <<
1064 << "p.z() = " << p.z()/mm << " mm"; <<
1065 message.precision(oldprc); <<
1066 G4Exception("G4TessellatedSolid::Inside <<
1067 "GeomSolids1002", JustWarning, messag <<
1068 } <<
1069 #endif <<
1070 // <<
1071 // In the next if-then-elseif G4String th <<
1072 // (1) You don't hit anything so cannot b <<
1073 // constructed correctly! <<
1074 // (2) Distance to inside (ie. nearest fa <<
1075 // shorter than distance to outside ( <<
1076 // facet) - on condition of safety di <<
1077 // (3) Distance to outside is shorter tha <<
1078 // we're inside. <<
1079 // <<
1080 if (distIn == kInfinity && distOut == kIn <<
1081 locationprime = kOutside; <<
1082 else if (distIn <= distOut - kCarToleranc <<
1083 locationprime = kOutside; <<
1084 else if (distOut <= distIn - kCarToleranc <<
1085 locationprime = kInside; <<
1086 <<
1087 if (i == 0) location = locationprime; <<
1088 } <<
1089 <<
1090 return location; <<
1091 } <<
1092 <<
1093 ///////////////////////////////////////////// <<
1094 // <<
1095 // Return index of the facet closest to the p <<
1096 // be located on the surface. Return -1 if no <<
1097 // <<
1098 G4int G4TessellatedSolid::GetFacetIndex (cons <<
1099 { 309 {
1100 G4int index = -1; << 310 G4double minDist = kInfinity;
1101 << 311 G4double dist = 0.0;
1102 if (fVoxels.GetCountOfVoxels() > 1) << 312 typedef std::multimap< G4double, FacetCI, std::less<G4double> > DistMapType;
1103 { << 313 DistMapType distmap;
1104 vector<G4int> curVoxel(3); << 314 size_t purgeIntv = 25;
1105 fVoxels.GetVoxel(curVoxel, p); << 315
1106 const vector<G4int> &candidates = fVoxels << 316 for (FacetCI f=facets.begin(); f!=facets.end(); f++)
1107 if (auto limit = (G4int)candidates.size() <<
1108 { <<
1109 G4double minDist = kInfinity; <<
1110 for(G4int i = 0 ; i < limit ; ++i) <<
1111 { <<
1112 G4int candidate = candidates[i]; <<
1113 G4VFacet& facet = *fFacets[candidate] <<
1114 G4double dist = facet.Distance(p, min <<
1115 if (dist <= kCarToleranceHalf) return <<
1116 if (dist < minDist) <<
1117 { <<
1118 minDist = dist; <<
1119 index = candidate; <<
1120 } <<
1121 } <<
1122 } <<
1123 } <<
1124 else <<
1125 { 317 {
1126 G4double minDist = kInfinity; << 318 dist = (*f)->Distance(p,minDist);
1127 std::size_t size = fFacets.size(); << 319 distmap.insert(DistMapType::value_type(dist,f));
1128 for (std::size_t i = 0; i < size; ++i) << 320 minDist = distmap.begin()->first;
>> 321 if (distmap.size() > purgeIntv)
1129 { 322 {
1130 G4VFacet& facet = *fFacets[i]; << 323 DistMapType::iterator it =
1131 G4double dist = facet.Distance(p, minDi << 324 distmap.lower_bound(minDist + 0.5*kCarTolerance);
1132 if (dist < minDist) << 325 it++;
>> 326 if (it != distmap.end())
1133 { 327 {
1134 minDist = dist; << 328 DistMapType::iterator itend = distmap.end();
1135 index = (G4int)i; << 329 itend--;
>> 330 distmap.erase (it,itend);
1136 } 331 }
>> 332 if (distmap.size() > purgeIntv) purgeIntv = 2*distmap.size();
1137 } 333 }
1138 } 334 }
1139 return index; <<
1140 } <<
1141 335
1142 ///////////////////////////////////////////// << 336 EInside inside = kInside;
1143 // << 337
1144 // Return the outwards pointing unit normal o << 338 if (minDist <= 0.5*kCarTolerance) {inside = kSurface;}
1145 // surface closest to the point at offset p. <<
1146 // <<
1147 G4bool G4TessellatedSolid::Normal (const G4Th <<
1148 G4Th <<
1149 { <<
1150 G4double minDist; <<
1151 G4VFacet* facet = nullptr; <<
1152 <<
1153 if (fVoxels.GetCountOfVoxels() > 1) <<
1154 { <<
1155 vector<G4int> curVoxel(3); <<
1156 fVoxels.GetVoxel(curVoxel, p); <<
1157 const vector<G4int> &candidates = fVoxels <<
1158 // fVoxels.GetCandidatesVoxelArray(p, can <<
1159 <<
1160 if (auto limit = (G4int)candidates.size() <<
1161 { <<
1162 minDist = kInfinity; <<
1163 for(G4int i = 0 ; i < limit ; ++i) <<
1164 { <<
1165 G4int candidate = candidates[i]; <<
1166 G4VFacet &fct = *fFacets[candidate]; <<
1167 G4double dist = fct.Distance(p,minDis <<
1168 if (dist < minDist) minDist = dist; <<
1169 if (dist <= kCarToleranceHalf) <<
1170 { <<
1171 aNormal = fct.GetSurfaceNormal(); <<
1172 return true; <<
1173 } <<
1174 } <<
1175 } <<
1176 minDist = MinDistanceFacet(p, true, facet <<
1177 } <<
1178 else 339 else
1179 { 340 {
1180 minDist = kInfinity; << 341 DistMapType::const_iterator itcut =
1181 std::size_t size = fFacets.size(); << 342 distmap.lower_bound(minDist + 0.5* kCarTolerance);
1182 for (std::size_t i = 0; i < size; ++i) << 343 itcut++;
>> 344 DistMapType::const_iterator it = distmap.begin();
>> 345 do
1183 { 346 {
1184 G4VFacet& f = *fFacets[i]; << 347 if (!((*(it->second))->IsInside(p))) {inside = kOutside;}
1185 G4double dist = f.Distance(p, minDist); << 348 } while (inside == kInside && ++it != itcut);
1186 if (dist < minDist) <<
1187 { <<
1188 minDist = dist; <<
1189 facet = &f; <<
1190 } <<
1191 } <<
1192 } 349 }
1193 350
1194 if (minDist != kInfinity) << 351 return inside;
1195 { <<
1196 if (facet != nullptr) { aNormal = facet- <<
1197 return minDist <= kCarToleranceHalf; <<
1198 } <<
1199 else <<
1200 { <<
1201 #ifdef G4VERBOSE <<
1202 std::ostringstream message; <<
1203 message << "Point p is not on surface !?" <<
1204 << " No facets found for point <<
1205 << " Returning approximated va <<
1206 <<
1207 G4Exception("G4TessellatedSolid::SurfaceN <<
1208 "GeomSolids1002", JustWarning <<
1209 #endif <<
1210 aNormal = (p.z() > 0 ? G4ThreeVector(0,0, <<
1211 return false; <<
1212 } <<
1213 } 352 }
1214 353
1215 ///////////////////////////////////////////// 354 ///////////////////////////////////////////////////////////////////////////////
1216 // 355 //
1217 // G4double DistanceToIn(const G4ThreeVector& << 356 G4ThreeVector G4TessellatedSolid::SurfaceNormal (const G4ThreeVector &p) const
1218 // <<
1219 // Return the distance along the normalised v <<
1220 // from the point at offset p. If there is no <<
1221 // kInfinity. The first intersection resultin <<
1222 // surface/volume is discarded. Hence, this i <<
1223 // surface of shape. <<
1224 // <<
1225 G4double <<
1226 G4TessellatedSolid::DistanceToInNoVoxels (con <<
1227 con <<
1228 <<
1229 { 357 {
1230 G4double minDist = kInfinity; << 358 FacetCI minFacet;
1231 G4double dist = 0.0; << 359 G4double minDist = kInfinity;
1232 G4double distFromSurface = 0.0; << 360 G4double dist = 0.0;
1233 G4ThreeVector normal; << 361
1234 << 362 for (FacetCI f=facets.begin(); f!=facets.end(); f++)
1235 #if G4SPECSDEBUG <<
1236 if (Inside(p) == kInside ) <<
1237 { 363 {
1238 std::ostringstream message; << 364 dist = (*f)->Distance(p,minDist);
1239 G4int oldprc = message.precision(16) ; << 365 if (dist < minDist)
1240 message << "Point p is already inside!?" <<
1241 << "Position:" << G4endl << G4endl <<
1242 << " p.x() = " << p.x()/mm << " mm" <<
1243 << " p.y() = " << p.y()/mm << " mm" <<
1244 << " p.z() = " << p.z()/mm << " mm" <<
1245 << "DistanceToOut(p) == " << DistanceTo <<
1246 message.precision(oldprc) ; <<
1247 G4Exception("G4TriangularFacet::DistanceT <<
1248 "GeomSolids1002", JustWarning <<
1249 } <<
1250 #endif <<
1251 <<
1252 std::size_t size = fFacets.size(); <<
1253 for (std::size_t i = 0; i < size; ++i) <<
1254 { <<
1255 G4VFacet& facet = *fFacets[i]; <<
1256 if (facet.Intersect(p,v,false,dist,distFr <<
1257 { 366 {
1258 // << 367 minDist = dist;
1259 // set minDist to the new distance to c << 368 minFacet = f;
1260 // is in positive direction and point i <<
1261 // within 0.5*kCarTolerance of the surf <<
1262 // zero and leave member function immed <<
1263 // proposed by & credit to Akira Okumur <<
1264 // <<
1265 if (distFromSurface > kCarToleranceHalf <<
1266 { <<
1267 minDist = dist; <<
1268 } <<
1269 else <<
1270 { <<
1271 if (-kCarToleranceHalf <= dist && dis <<
1272 { <<
1273 return 0.0; <<
1274 } <<
1275 else <<
1276 { <<
1277 if (distFromSurface > -kCarToleran <<
1278 && distFromSurface < kCarToleran <<
1279 { <<
1280 minDist = dist; <<
1281 } <<
1282 } <<
1283 } <<
1284 } 369 }
1285 } 370 }
1286 return minDist; << 371
>> 372 return (*minFacet)->GetSurfaceNormal();
1287 } 373 }
1288 374
1289 ///////////////////////////////////////////// 375 ///////////////////////////////////////////////////////////////////////////////
1290 // 376 //
1291 G4double << 377 G4double G4TessellatedSolid::DistanceToIn (const G4ThreeVector &p,
1292 G4TessellatedSolid::DistanceToOutNoVoxels (co << 378 const G4ThreeVector &v) const
1293 co <<
1294 <<
1295 <<
1296 <<
1297 { 379 {
1298 G4double minDist = kInfinity; 380 G4double minDist = kInfinity;
1299 G4double dist = 0.0; 381 G4double dist = 0.0;
1300 G4double distFromSurface = 0.0; 382 G4double distFromSurface = 0.0;
1301 G4ThreeVector normal, minNormal; << 383 G4ThreeVector normal(0.0,0.0,0.0);
1302 << 384
1303 #if G4SPECSDEBUG << 385 for (FacetCI f=facets.begin(); f!=facets.end(); f++)
1304 if ( Inside(p) == kOutside ) <<
1305 { <<
1306 std::ostringstream message; <<
1307 G4int oldprc = message.precision(16) ; <<
1308 message << "Point p is already outside!?" <<
1309 << "Position:" << G4endl << G4endl <<
1310 << " p.x() = " << p.x()/mm << " mm" <<
1311 << " p.y() = " << p.y()/mm << " mm" <<
1312 << " p.z() = " << p.z()/mm << " mm" <<
1313 << "DistanceToIn(p) == " << DistanceToI <<
1314 message.precision(oldprc) ; <<
1315 G4Exception("G4TriangularFacet::DistanceT <<
1316 "GeomSolids1002", JustWarning <<
1317 } <<
1318 #endif <<
1319 <<
1320 G4bool isExtreme = false; <<
1321 std::size_t size = fFacets.size(); <<
1322 for (std::size_t i = 0; i < size; ++i) <<
1323 { 386 {
1324 G4VFacet& facet = *fFacets[i]; << 387 if ((*f)->Intersect(p,v,false,dist,distFromSurface,normal))
1325 if (facet.Intersect(p,v,true,dist,distFro <<
1326 { 388 {
1327 if (distFromSurface > 0.0 && distFromSu << 389 if (dist < minDist) minDist = dist;
1328 && facet.Distance(p,kCarTolerance) <= <<
1329 { <<
1330 // We are on a surface. Return zero. <<
1331 aConvex = (fExtremeFacets.find(&facet <<
1332 // Normal(p, aNormalVector); <<
1333 // aNormalVector = facet.GetSurfaceNo <<
1334 aNormalVector = normal; <<
1335 return 0.0; <<
1336 } <<
1337 if (dist >= 0.0 && dist < minDist) <<
1338 { <<
1339 minDist = dist; <<
1340 minNormal = normal; <<
1341 isExtreme = (fExtremeFacets.find(&fac <<
1342 } <<
1343 } 390 }
1344 } 391 }
1345 if (minDist < kInfinity) <<
1346 { <<
1347 aNormalVector = minNormal; <<
1348 aConvex = isExtreme; <<
1349 return minDist; <<
1350 } <<
1351 else <<
1352 { <<
1353 // No intersection found <<
1354 aConvex = false; <<
1355 Normal(p, aNormalVector); <<
1356 return 0.0; <<
1357 } <<
1358 } <<
1359 <<
1360 ///////////////////////////////////////////// <<
1361 // <<
1362 void G4TessellatedSolid:: <<
1363 DistanceToOutCandidates(const std::vector<G4i <<
1364 const G4ThreeVector& <<
1365 const G4ThreeVector& <<
1366 G4double& minDi <<
1367 G4int& minCandi <<
1368 { <<
1369 auto candidatesCount = (G4int)candidates.si <<
1370 G4double dist = 0.0; <<
1371 G4double distFromSurface = 0.0; <<
1372 G4ThreeVector normal; <<
1373 392
1374 for (G4int i = 0 ; i < candidatesCount; ++i << 393 return minDist;
1375 { <<
1376 G4int candidate = candidates[i]; <<
1377 G4VFacet& facet = *fFacets[candidate]; <<
1378 if (facet.Intersect(aPoint,direction,true <<
1379 { <<
1380 if (distFromSurface > 0.0 && distFromSu <<
1381 && facet.Distance(aPoint,kCarTolerance <<
1382 { <<
1383 // We are on a surface <<
1384 // <<
1385 minDist = 0.0; <<
1386 minNormal = normal; <<
1387 minCandidate = candidate; <<
1388 break; <<
1389 } <<
1390 if (dist >= 0.0 && dist < minDist) <<
1391 { <<
1392 minDist = dist; <<
1393 minNormal = normal; <<
1394 minCandidate = candidate; <<
1395 } <<
1396 } <<
1397 } <<
1398 } 394 }
1399 395
1400 ///////////////////////////////////////////// 396 ///////////////////////////////////////////////////////////////////////////////
1401 // 397 //
1402 G4double << 398 G4double G4TessellatedSolid::DistanceToIn (const G4ThreeVector &p) const
1403 G4TessellatedSolid::DistanceToOutCore(const G <<
1404 const G <<
1405 G <<
1406 G <<
1407 G <<
1408 { 399 {
1409 G4double minDistance; << 400 G4double minDist = kInfinity;
1410 << 401 G4double dist = 0.0;
1411 if (fVoxels.GetCountOfVoxels() > 1) << 402
1412 { << 403 for (FacetCI f=facets.begin(); f!=facets.end(); f++)
1413 minDistance = kInfinity; <<
1414 <<
1415 G4ThreeVector currentPoint = aPoint; <<
1416 G4ThreeVector direction = aDirection.unit <<
1417 G4double totalShift = 0.; <<
1418 vector<G4int> curVoxel(3); <<
1419 if (!fVoxels.Contains(aPoint)) return 0.; <<
1420 <<
1421 fVoxels.GetVoxel(curVoxel, currentPoint); <<
1422 <<
1423 G4double shiftBonus = kCarTolerance; <<
1424 <<
1425 const vector<G4int>* old = nullptr; <<
1426 <<
1427 G4int minCandidate = -1; <<
1428 do // Loop checking, 13.08.2015, G.Cos <<
1429 { <<
1430 const vector<G4int>& candidates = fVoxe <<
1431 if (old == &candidates) <<
1432 ++old; <<
1433 if (old != &candidates && !candidates.e <<
1434 { <<
1435 DistanceToOutCandidates(candidates, a <<
1436 aNormalVector <<
1437 if (minDistance <= totalShift) break; <<
1438 } <<
1439 <<
1440 G4double shift=fVoxels.DistanceToNext(c <<
1441 if (shift == kInfinity) break; <<
1442 <<
1443 totalShift += shift; <<
1444 if (minDistance <= totalShift) break; <<
1445 <<
1446 currentPoint += direction * (shift + sh <<
1447 <<
1448 old = &candidates; <<
1449 } <<
1450 while (fVoxels.UpdateCurrentVoxel(current <<
1451 <<
1452 if (minCandidate < 0) <<
1453 { <<
1454 // No intersection found <<
1455 minDistance = 0.; <<
1456 aConvex = false; <<
1457 Normal(aPoint, aNormalVector); <<
1458 } <<
1459 else <<
1460 { <<
1461 aConvex = (fExtremeFacets.find(fFacets[ <<
1462 != fExtremeFacets.end()); <<
1463 } <<
1464 } <<
1465 else <<
1466 { 404 {
1467 minDistance = DistanceToOutNoVoxels(aPoin << 405 dist = (*f)->Distance(p,minDist,false);
1468 aConv << 406 if (dist < minDist) minDist = dist;
1469 } 407 }
1470 return minDistance; << 408
>> 409 return minDist;
1471 } 410 }
1472 411
1473 ///////////////////////////////////////////// 412 ///////////////////////////////////////////////////////////////////////////////
1474 // 413 //
1475 G4double G4TessellatedSolid:: << 414 G4double G4TessellatedSolid::DistanceToOut (const G4ThreeVector &p,
1476 DistanceToInCandidates(const std::vector<G4in << 415 const G4ThreeVector &v, const G4bool calcNorm,
1477 const G4ThreeVector& a << 416 G4bool *validNorm, G4ThreeVector *n) const
1478 const G4ThreeVector& d <<
1479 { 417 {
1480 auto candidatesCount = (G4int)candidates.si << 418 G4double minDist1 = kInfinity;
>> 419 G4double minDist2 = kInfinity;
1481 G4double dist = 0.0; 420 G4double dist = 0.0;
1482 G4double distFromSurface = 0.0; 421 G4double distFromSurface = 0.0;
1483 G4ThreeVector normal; << 422 G4ThreeVector normal(0.0,0.0,0.0);
1484 << 423 G4ThreeVector minNormal1(0.0,0.0,0.0);
1485 G4double minDistance = kInfinity; << 424 G4ThreeVector minNormal2(0.0,0.0,0.0);
1486 for (G4int i = 0 ; i < candidatesCount; ++i << 425
>> 426 for (FacetCI f=facets.begin(); f!=facets.end(); f++)
1487 { 427 {
1488 G4int candidate = candidates[i]; << 428 if ((*f)->Intersect(p,v,true,dist,distFromSurface,normal))
1489 G4VFacet& facet = *fFacets[candidate]; <<
1490 if (facet.Intersect(aPoint,direction,fals <<
1491 { 429 {
1492 // << 430 if (dist < minDist1)
1493 // Set minDist to the new distance to c <<
1494 // in positive direction and point is n <<
1495 // within 0.5*kCarTolerance of the surf <<
1496 // zero and leave member function immed <<
1497 // proposed by & credit to Akira Okumur <<
1498 // <<
1499 if ( (distFromSurface > kCarToleranceHa <<
1500 && (dist >= 0.0) && (dist < minDistan <<
1501 { 431 {
1502 minDistance = dist; << 432 if (v.dot(normal) > dirTolerance)
1503 } <<
1504 else <<
1505 { <<
1506 if (-kCarToleranceHalf <= dist && dis <<
1507 { 433 {
1508 return 0.0; << 434 minDist1 = dist;
>> 435 minNormal1 = normal;
1509 } 436 }
1510 else if (distFromSurface > -kCarTole << 437 else if (dist < minDist2)
1511 && distFromSurface < kCarTole <<
1512 { 438 {
1513 minDistance = dist; << 439 minDist2 = dist;
>> 440 minNormal2 = normal;
1514 } 441 }
1515 } 442 }
1516 } 443 }
1517 } 444 }
1518 return minDistance; << 445
1519 } << 446 if (minDist1 < kInfinity)
1520 <<
1521 ///////////////////////////////////////////// <<
1522 // <<
1523 G4double <<
1524 G4TessellatedSolid::DistanceToInCore(const G4 <<
1525 const G4 <<
1526 G4 <<
1527 { <<
1528 G4double minDistance; <<
1529 <<
1530 if (fVoxels.GetCountOfVoxels() > 1) <<
1531 { 447 {
1532 minDistance = kInfinity; << 448 if (calcNorm)
1533 G4ThreeVector currentPoint = aPoint; <<
1534 G4ThreeVector direction = aDirection.unit <<
1535 G4double shift = fVoxels.DistanceToFirst( <<
1536 if (shift == kInfinity) return shift; <<
1537 G4double shiftBonus = kCarTolerance; <<
1538 if (shift != 0.0) <<
1539 currentPoint += direction * (shift + sh <<
1540 // if (!fVoxels.Contains(currentPoint)) <<
1541 G4double totalShift = shift; <<
1542 <<
1543 // G4SurfBits exclusion; // (1/*fVoxels.G <<
1544 vector<G4int> curVoxel(3); <<
1545 <<
1546 fVoxels.GetVoxel(curVoxel, currentPoint); <<
1547 do // Loop checking, 13.08.2015, G.Cos <<
1548 { 449 {
1549 const vector<G4int>& candidates = fVoxe << 450 *validNorm = true;
1550 if (!candidates.empty()) << 451 *n = minNormal1;
1551 { <<
1552 G4double distance=DistanceToInCandida <<
1553 if (minDistance > distance) minDistan <<
1554 if (distance < totalShift) break; <<
1555 } <<
1556 <<
1557 shift = fVoxels.DistanceToNext(currentP <<
1558 if (shift == kInfinity /*|| shift == 0* <<
1559 <<
1560 totalShift += shift; <<
1561 if (minDistance < totalShift) break; <<
1562 <<
1563 currentPoint += direction * (shift + sh <<
1564 } 452 }
1565 while (fVoxels.UpdateCurrentVoxel(current << 453 return minDist1;
1566 } 454 }
1567 else 455 else
1568 { 456 {
1569 minDistance = DistanceToInNoVoxels(aPoint << 457 if (calcNorm)
1570 } <<
1571 <<
1572 return minDistance; <<
1573 } <<
1574 <<
1575 ///////////////////////////////////////////// <<
1576 // <<
1577 G4bool <<
1578 G4TessellatedSolid::CompareSortedVoxel(const <<
1579 const <<
1580 { <<
1581 return l.second < r.second; <<
1582 } <<
1583 <<
1584 ///////////////////////////////////////////// <<
1585 // <<
1586 G4double <<
1587 G4TessellatedSolid::MinDistanceFacet(const G4 <<
1588 G4 <<
1589 G4 <<
1590 { <<
1591 G4double minDist = kInfinity; <<
1592 <<
1593 G4int size = fVoxels.GetVoxelBoxesSize(); <<
1594 vector<pair<G4int, G4double> > voxelsSorted <<
1595 <<
1596 pair<G4int, G4double> info; <<
1597 <<
1598 for (G4int i = 0; i < size; ++i) <<
1599 { <<
1600 const G4VoxelBox& voxelBox = fVoxels.GetV <<
1601 <<
1602 G4ThreeVector pointShifted = p - voxelBox <<
1603 G4double safety = fVoxels.MinDistanceToBo <<
1604 info.first = i; <<
1605 info.second = safety; <<
1606 voxelsSorted[i] = info; <<
1607 } <<
1608 <<
1609 std::sort(voxelsSorted.begin(), voxelsSorte <<
1610 &G4TessellatedSolid::CompareSorte <<
1611 <<
1612 for (G4int i = 0; i < size; ++i) <<
1613 { <<
1614 const pair<G4int,G4double>& inf = voxelsS <<
1615 G4double dist = inf.second; <<
1616 if (dist > minDist) break; <<
1617 <<
1618 const vector<G4int>& candidates = fVoxels <<
1619 auto csize = (G4int)candidates.size(); <<
1620 for (G4int j = 0; j < csize; ++j) <<
1621 { 458 {
1622 G4int candidate = candidates[j]; << 459 *validNorm = true;
1623 G4VFacet& facet = *fFacets[candidate]; << 460 *n = minNormal2;
1624 dist = simple ? facet.Distance(p,minDis <<
1625 : facet.Distance(p,minDis <<
1626 if (dist < minDist) <<
1627 { <<
1628 minDist = dist; <<
1629 minFacet = &facet; <<
1630 } <<
1631 } 461 }
>> 462 return minDist2;
1632 } 463 }
1633 return minDist; <<
1634 } 464 }
1635 465
1636 ///////////////////////////////////////////// 466 ///////////////////////////////////////////////////////////////////////////////
1637 // 467 //
1638 G4double G4TessellatedSolid::SafetyFromOutsid << 468 G4double G4TessellatedSolid::DistanceToOut (const G4ThreeVector &p) const
1639 <<
1640 { 469 {
1641 #if G4SPECSDEBUG << 470 G4double minDist = kInfinity;
1642 if ( Inside(p) == kInside ) << 471 G4double dist = 0.0;
1643 { << 472
1644 std::ostringstream message; << 473 for (FacetCI f=facets.begin(); f!=facets.end(); f++)
1645 G4int oldprc = message.precision(16) ; <<
1646 message << "Point p is already inside!?" <<
1647 << "Position:" << G4endl << G4endl <<
1648 << "p.x() = " << p.x()/mm << " mm" << <<
1649 << "p.y() = " << p.y()/mm << " mm" << <<
1650 << "p.z() = " << p.z()/mm << " mm" << <<
1651 << "DistanceToOut(p) == " << DistanceTo <<
1652 message.precision(oldprc) ; <<
1653 G4Exception("G4TriangularFacet::DistanceT <<
1654 "GeomSolids1002", JustWarning <<
1655 } <<
1656 #endif <<
1657 <<
1658 G4double minDist; <<
1659 <<
1660 if (fVoxels.GetCountOfVoxels() > 1) <<
1661 { <<
1662 if (!aAccurate) <<
1663 return fVoxels.DistanceToBoundingBox(p) <<
1664 <<
1665 if (!OutsideOfExtent(p, kCarTolerance)) <<
1666 { <<
1667 vector<G4int> startingVoxel(3); <<
1668 fVoxels.GetVoxel(startingVoxel, p); <<
1669 const vector<G4int> &candidates = fVoxe <<
1670 if (candidates.empty() && (fInsides.Get <<
1671 { <<
1672 G4int index = fVoxels.GetPointIndex(p <<
1673 if (fInsides[index]) return 0.; <<
1674 } <<
1675 } <<
1676 <<
1677 G4VFacet* facet; <<
1678 minDist = MinDistanceFacet(p, true, facet <<
1679 } <<
1680 else <<
1681 { 474 {
1682 minDist = kInfinity; << 475 dist = (*f)->Distance(p,minDist,true);
1683 std::size_t size = fFacets.size(); << 476 if (dist < minDist) minDist = dist;
1684 for (std::size_t i = 0; i < size; ++i) <<
1685 { <<
1686 G4VFacet& facet = *fFacets[i]; <<
1687 G4double dist = facet.Distance(p,minDis <<
1688 if (dist < minDist) minDist = dist; <<
1689 } <<
1690 } 477 }
>> 478
1691 return minDist; 479 return minDist;
1692 } 480 }
1693 481
1694 ///////////////////////////////////////////// 482 ///////////////////////////////////////////////////////////////////////////////
1695 // 483 //
1696 G4double << 484 G4GeometryType G4TessellatedSolid::GetEntityType () const
1697 G4TessellatedSolid::SafetyFromInside (const G <<
1698 { 485 {
1699 #if G4SPECSDEBUG << 486 return geometryType;
1700 if ( Inside(p) == kOutside ) <<
1701 { <<
1702 std::ostringstream message; <<
1703 G4int oldprc = message.precision(16) ; <<
1704 message << "Point p is already outside!?" <<
1705 << "Position:" << G4endl << G4endl <<
1706 << "p.x() = " << p.x()/mm << " mm" << <<
1707 << "p.y() = " << p.y()/mm << " mm" << <<
1708 << "p.z() = " << p.z()/mm << " mm" << <<
1709 << "DistanceToIn(p) == " << DistanceToI <<
1710 message.precision(oldprc) ; <<
1711 G4Exception("G4TriangularFacet::DistanceT <<
1712 "GeomSolids1002", JustWarning <<
1713 } <<
1714 #endif <<
1715 <<
1716 G4double minDist; <<
1717 <<
1718 if (OutsideOfExtent(p, kCarTolerance)) retu <<
1719 <<
1720 if (fVoxels.GetCountOfVoxels() > 1) <<
1721 { <<
1722 G4VFacet* facet; <<
1723 minDist = MinDistanceFacet(p, true, facet <<
1724 } <<
1725 else <<
1726 { <<
1727 minDist = kInfinity; <<
1728 G4double dist = 0.0; <<
1729 std::size_t size = fFacets.size(); <<
1730 for (std::size_t i = 0; i < size; ++i) <<
1731 { <<
1732 G4VFacet& facet = *fFacets[i]; <<
1733 dist = facet.Distance(p,minDist); <<
1734 if (dist < minDist) minDist = dist; <<
1735 } <<
1736 } <<
1737 return minDist; <<
1738 } 487 }
1739 488
1740 ///////////////////////////////////////////// 489 ///////////////////////////////////////////////////////////////////////////////
1741 // 490 //
1742 // G4GeometryType GetEntityType() const; << 491 void G4TessellatedSolid::DescribeYourselfTo (G4VGraphicsScene& scene) const
1743 // <<
1744 // Provide identification of the class of an <<
1745 // <<
1746 G4GeometryType G4TessellatedSolid::GetEntityT <<
1747 { 492 {
1748 return fGeometryType; << 493 scene.AddSolid (*this);
1749 } 494 }
1750 495
1751 ///////////////////////////////////////////// 496 ///////////////////////////////////////////////////////////////////////////////
1752 // 497 //
1753 // IsFaceted << 498 // Dispatch to parameterisation for replication mechanism dimension
1754 // << 499 // computation & modification.
1755 G4bool G4TessellatedSolid::IsFaceted () const << 500 //
1756 { << 501 //void G4TessellatedSolid::ComputeDimensions (G4VPVParameterisation* p,
1757 return true; << 502 // const G4int n, const G4VPhysicalVolume* pRep) const
1758 } << 503 //{
>> 504 // G4VSolid *ptr = 0;
>> 505 // ptr = *this;
>> 506 // p->ComputeDimensions(ptr,n,pRep);
>> 507 //}
1759 508
1760 ///////////////////////////////////////////// 509 ///////////////////////////////////////////////////////////////////////////////
1761 // 510 //
1762 std::ostream &G4TessellatedSolid::StreamInfo( 511 std::ostream &G4TessellatedSolid::StreamInfo(std::ostream &os) const
1763 { 512 {
1764 os << G4endl; 513 os << G4endl;
1765 os << "Solid name = " << GetName() << 514 os << "Geometry Type = " << geometryType << G4endl;
1766 os << "Geometry Type = " << fGeometryTyp << 515 os << "Number of facets = " << facets.size() << G4endl;
1767 os << "Number of facets = " << fFacets.size << 516
1768 << 517 for (FacetCI f = facets.begin(); f != facets.end(); f++)
1769 std::size_t size = fFacets.size(); <<
1770 for (std::size_t i = 0; i < size; ++i) <<
1771 { 518 {
1772 os << "FACET # = " << i + 1 << G << 519 os << "FACET # = " << f-facets.begin()+1 << G4endl;
1773 G4VFacet &facet = *fFacets[i]; << 520 (*f)->StreamInfo(os);
1774 facet.StreamInfo(os); <<
1775 } 521 }
1776 os << G4endl; << 522 os <<G4endl;
1777 << 523
1778 return os; 524 return os;
1779 } 525 }
1780 526
1781 ///////////////////////////////////////////// 527 ///////////////////////////////////////////////////////////////////////////////
1782 // 528 //
1783 // Make a clone of the object << 529 G4Polyhedron *G4TessellatedSolid::CreatePolyhedron () const
1784 // <<
1785 G4VSolid* G4TessellatedSolid::Clone() const <<
1786 { <<
1787 return new G4TessellatedSolid(*this); <<
1788 } <<
1789 <<
1790 ///////////////////////////////////////////// <<
1791 // <<
1792 // EInside G4TessellatedSolid::Inside (const <<
1793 // <<
1794 // This method must return: <<
1795 // * kOutside if the point at offset p is <<
1796 // boundaries plus kCarTolerance/2, <<
1797 // * kSurface if the point is <= kCarToler <<
1798 // * kInside otherwise. <<
1799 // <<
1800 EInside G4TessellatedSolid::Inside (const G4T <<
1801 { <<
1802 EInside location; <<
1803 <<
1804 if (fVoxels.GetCountOfVoxels() > 1) <<
1805 { <<
1806 location = InsideVoxels(aPoint); <<
1807 } <<
1808 else <<
1809 { <<
1810 location = InsideNoVoxels(aPoint); <<
1811 } <<
1812 return location; <<
1813 } <<
1814 <<
1815 ///////////////////////////////////////////// <<
1816 // <<
1817 G4ThreeVector G4TessellatedSolid::SurfaceNorm <<
1818 { <<
1819 G4ThreeVector n; <<
1820 Normal(p, n); <<
1821 return n; <<
1822 } <<
1823 <<
1824 ///////////////////////////////////////////// <<
1825 // <<
1826 // G4double DistanceToIn(const G4ThreeVector& <<
1827 // <<
1828 // Calculate distance to nearest surface of s <<
1829 // distance can be an underestimate. <<
1830 // <<
1831 G4double G4TessellatedSolid::DistanceToIn(con <<
1832 { <<
1833 return SafetyFromOutside(p, false); <<
1834 } <<
1835 <<
1836 ///////////////////////////////////////////// <<
1837 // <<
1838 G4double G4TessellatedSolid::DistanceToIn(con <<
1839 con <<
1840 { 530 {
1841 G4double dist = DistanceToInCore(p,v,kInfin << 531 size_t nVertices = vertexList.size();
1842 #ifdef G4SPECSDEBUG << 532 size_t nFacets = facets.size();
1843 if (dist < kInfinity) << 533 G4PolyhedronArbitrary *polyhedron =
>> 534 new G4PolyhedronArbitrary (nVertices, nFacets);
>> 535 for (G4ThreeVectorList::const_iterator v = vertexList.begin();
>> 536 v!=vertexList.end(); v++) polyhedron->AddVertex(*v);
>> 537
>> 538 for (FacetCI f=facets.begin(); f != facets.end(); f++)
1844 { 539 {
1845 if (Inside(p + dist*v) != kSurface) << 540 size_t v[4];
>> 541 for (size_t j=0; j<4; j++)
1846 { 542 {
1847 std::ostringstream message; << 543 size_t i = (*f)->GetVertexIndex(j);
1848 message << "Invalid response from facet << 544 if (i == 999999999) v[j] = 0;
1849 << G4endl << 545 else v[j] = i+1;
1850 << "at point: " << p << "and d <<
1851 G4Exception("G4TessellatedSolid::Distan <<
1852 "GeomSolids1002", JustWarni <<
1853 } 546 }
1854 } << 547 if ((*f)->GetEntityType() == "G4RectangularFacet")
1855 #endif <<
1856 return dist; <<
1857 } <<
1858 <<
1859 ///////////////////////////////////////////// <<
1860 // <<
1861 // G4double DistanceToOut(const G4ThreeVector <<
1862 // <<
1863 // Calculate distance to nearest surface of s <<
1864 // point. The distance can be an underestimat <<
1865 // <<
1866 G4double G4TessellatedSolid::DistanceToOut(co <<
1867 { <<
1868 return SafetyFromInside(p, false); <<
1869 } <<
1870 <<
1871 ///////////////////////////////////////////// <<
1872 // <<
1873 // G4double DistanceToOut(const G4ThreeVector <<
1874 // const G4bool calcNo <<
1875 // G4bool *validNorm=0 <<
1876 // <<
1877 // Return distance along the normalised vecto <<
1878 // point at an offset p inside or on the surf <<
1879 // shape. Intersections with surfaces, when t <<
1880 // than kCarTolerance/2 from a surface, must <<
1881 // If calcNorm is true, then it must also <<
1882 // * true, if the solid lies entirely beh <<
1883 // surface. Then it must set n to the <<
1884 // (the Magnitude of the vector is not <<
1885 // * false, if the solid does not lie ent <<
1886 // exiting surface. <<
1887 // If calcNorm is false, then validNorm and n <<
1888 // <<
1889 G4double G4TessellatedSolid::DistanceToOut(co <<
1890 co <<
1891 co <<
1892 <<
1893 <<
1894 { <<
1895 G4ThreeVector n; <<
1896 G4bool valid; <<
1897 <<
1898 G4double dist = DistanceToOutCore(p, v, n, <<
1899 if (calcNorm) <<
1900 { <<
1901 *norm = n; <<
1902 *validNorm = valid; <<
1903 } <<
1904 #ifdef G4SPECSDEBUG <<
1905 if (dist < kInfinity) <<
1906 { <<
1907 if (Inside(p + dist*v) != kSurface) <<
1908 { 548 {
1909 std::ostringstream message; << 549 size_t i = v[3];
1910 message << "Invalid response from facet << 550 v[3] = v[2];
1911 << G4endl << 551 v[2] = i;
1912 << "at point: " << p << "and d <<
1913 G4Exception("G4TessellatedSolid::Distan <<
1914 "GeomSolids1002", JustWarni <<
1915 } 552 }
>> 553 polyhedron->AddFacet(v[0],v[1],v[2],v[3]);
1916 } 554 }
1917 #endif << 555
1918 return dist; << 556 return (G4Polyhedron*) polyhedron;
1919 } 557 }
1920 558
1921 ///////////////////////////////////////////// 559 ///////////////////////////////////////////////////////////////////////////////
1922 // 560 //
1923 void G4TessellatedSolid::DescribeYourselfTo ( << 561 G4NURBS *G4TessellatedSolid::CreateNURBS () const
1924 { 562 {
1925 scene.AddSolid (*this); << 563 return 0;
1926 } <<
1927 <<
1928 ///////////////////////////////////////////// <<
1929 // <<
1930 G4Polyhedron* G4TessellatedSolid::CreatePolyh <<
1931 { <<
1932 auto nVertices = (G4int)fVertexList.size(); <<
1933 auto nFacets = (G4int)fFacets.size(); <<
1934 auto polyhedron = new G4Polyhedron(nVertice <<
1935 for (auto i = 0; i < nVertices; ++i) <<
1936 { <<
1937 polyhedron->SetVertex(i+1, fVertexList[i] <<
1938 } <<
1939 <<
1940 for (auto i = 0; i < nFacets; ++i) <<
1941 { <<
1942 G4VFacet* facet = fFacets[i]; <<
1943 G4int v[4] = {0}; <<
1944 G4int n = facet->GetNumberOfVertices(); <<
1945 if (n > 4) n = 4; <<
1946 for (auto j = 0; j < n; ++j) <<
1947 { <<
1948 v[j] = facet->GetVertexIndex(j) + 1; <<
1949 } <<
1950 polyhedron->SetFacet(i+1, v[0], v[1], v[2 <<
1951 } <<
1952 polyhedron->SetReferences(); <<
1953 <<
1954 return polyhedron; <<
1955 } 564 }
1956 565
1957 ///////////////////////////////////////////// 566 ///////////////////////////////////////////////////////////////////////////////
1958 // 567 //
1959 // GetPolyhedron 568 // GetPolyhedron
1960 // 569 //
1961 G4Polyhedron* G4TessellatedSolid::GetPolyhedr << 570 G4Polyhedron* G4TessellatedSolid::GetPolyhedron () const
1962 { 571 {
1963 if (fpPolyhedron == nullptr || << 572 if (!fpPolyhedron ||
1964 fRebuildPolyhedron || <<
1965 fpPolyhedron->GetNumberOfRotationStepsA 573 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
1966 fpPolyhedron->GetNumberOfRotationSteps( 574 fpPolyhedron->GetNumberOfRotationSteps())
1967 { << 575 {
1968 G4AutoLock l(&polyhedronMutex); << 576 delete fpPolyhedron;
1969 delete fpPolyhedron; << 577 fpPolyhedron = CreatePolyhedron();
1970 fpPolyhedron = CreatePolyhedron(); << 578 }
1971 fRebuildPolyhedron = false; <<
1972 l.unlock(); <<
1973 } <<
1974 return fpPolyhedron; 579 return fpPolyhedron;
1975 } 580 }
1976 581
1977 ///////////////////////////////////////////// 582 ///////////////////////////////////////////////////////////////////////////////
1978 // 583 //
1979 // Get bounding box << 584 G4bool G4TessellatedSolid::CalculateExtent(const EAxis pAxis,
1980 // << 585 const G4VoxelLimits& pVoxelLimit, const G4AffineTransform& pTransform,
1981 void G4TessellatedSolid::BoundingLimits(G4Thr << 586 G4double& pMin, G4double& pMax) const
1982 G4Thr << 587 {
1983 { << 588 if (!pTransform.IsRotated())
1984 pMin = fMinExtent; << 589 {
1985 pMax = fMaxExtent; << 590 G4double xoffset,xMin,xMax;
>> 591 G4double yoffset,yMin,yMax;
>> 592 G4double zoffset,zMin,zMax;
>> 593
>> 594 xoffset = pTransform.NetTranslation().x();
>> 595 xMin = xoffset + xMinExtent;
>> 596 xMax = xoffset + xMaxExtent;
>> 597
>> 598 if (pVoxelLimit.IsXLimited())
>> 599 {
>> 600 if ( xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance ||
>> 601 xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance ) return false ;
>> 602 else
>> 603 {
>> 604 if (xMin < pVoxelLimit.GetMinXExtent())
>> 605 {
>> 606 xMin = pVoxelLimit.GetMinXExtent() ;
>> 607 }
>> 608 if (xMax > pVoxelLimit.GetMaxXExtent())
>> 609 {
>> 610 xMax = pVoxelLimit.GetMaxXExtent() ;
>> 611 }
>> 612 }
>> 613 }
>> 614
>> 615 yoffset = pTransform.NetTranslation().y();
>> 616 yMin = yoffset + yMinExtent;
>> 617 yMax = yoffset + yMaxExtent;
>> 618
>> 619 if (pVoxelLimit.IsYLimited())
>> 620 {
>> 621 if ( yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance ||
>> 622 yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance ) return false ;
>> 623 else
>> 624 {
>> 625 if (yMin < pVoxelLimit.GetMinYExtent())
>> 626 {
>> 627 yMin = pVoxelLimit.GetMinYExtent() ;
>> 628 }
>> 629 if (yMax > pVoxelLimit.GetMaxYExtent())
>> 630 {
>> 631 yMax = pVoxelLimit.GetMaxYExtent() ;
>> 632 }
>> 633 }
>> 634 }
1986 635
1987 // Check correctness of the bounding box << 636 zoffset = pTransform.NetTranslation().z();
1988 // << 637 zMin = zoffset + zMinExtent;
1989 if (pMin.x() >= pMax.x() || pMin.y() >= pMa << 638 zMax = zoffset + zMaxExtent;
1990 { << 639
1991 std::ostringstream message; << 640 if (pVoxelLimit.IsZLimited())
1992 message << "Bad bounding box (min >= max) << 641 {
1993 << GetName() << " !" << 642 if ( zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance ||
1994 << "\npMin = " << pMin << 643 zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance ) return false ;
1995 << "\npMax = " << pMax; << 644 else
1996 G4Exception("G4TessellatedSolid::Bounding << 645 {
1997 "GeomMgt0001", JustWarning, m << 646 if (zMin < pVoxelLimit.GetMinZExtent())
1998 DumpInfo(); << 647 {
1999 } << 648 zMin = pVoxelLimit.GetMinZExtent() ;
2000 } << 649 }
>> 650 if (zMax > pVoxelLimit.GetMaxZExtent())
>> 651 {
>> 652 zMax = pVoxelLimit.GetMaxZExtent() ;
>> 653 }
>> 654 }
>> 655 }
2001 656
2002 ///////////////////////////////////////////// << 657 switch (pAxis)
2003 // << 658 {
2004 // Calculate extent under transform and speci << 659 case kXAxis:
2005 // << 660 pMin = xMin ;
2006 G4bool << 661 pMax = xMax ;
2007 G4TessellatedSolid::CalculateExtent(const EAx << 662 break ;
2008 const G4V << 663 case kYAxis:
2009 const G4A << 664 pMin=yMin;
2010 G4d << 665 pMax=yMax;
2011 { << 666 break;
2012 G4ThreeVector bmin, bmax; << 667 case kZAxis:
2013 << 668 pMin=zMin;
2014 // Check bounding box (bbox) << 669 pMax=zMax;
2015 // << 670 break;
2016 BoundingLimits(bmin,bmax); << 671 default:
2017 G4BoundingEnvelope bbox(bmin,bmax); << 672 break;
2018 << 673 }
2019 // Use simple bounding-box to help in the c << 674 pMin -= kCarTolerance ;
2020 // << 675 pMax += kCarTolerance ;
2021 return bbox.CalculateExtent(pAxis,pVoxelLim << 676
2022 << 677 return true;
2023 #if 0 <<
2024 // Precise extent computation (disabled by <<
2025 // <<
2026 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVo <<
2027 { <<
2028 return (pMin < pMax) ? true : false; <<
2029 } 678 }
2030 << 679 else
2031 // The extent is calculated as cumulative e <<
2032 // formed by facets and the center of the b <<
2033 // <<
2034 G4double eminlim = pVoxelLimit.GetMinExtent <<
2035 G4double emaxlim = pVoxelLimit.GetMaxExtent <<
2036 <<
2037 G4ThreeVectorList base; <<
2038 G4ThreeVectorList apex(1); <<
2039 std::vector<const G4ThreeVectorList *> pyra <<
2040 pyramid[0] = &base; <<
2041 pyramid[1] = &apex; <<
2042 apex[0] = (bmin+bmax)*0.5; <<
2043 <<
2044 // main loop along facets <<
2045 pMin = kInfinity; <<
2046 pMax = -kInfinity; <<
2047 for (G4int i=0; i<GetNumberOfFacets(); ++i) <<
2048 { 680 {
2049 G4VFacet* facet = GetFacet(i); <<
2050 if (std::abs((facet->GetSurfaceNormal()). <<
2051 < kCarToleranceHalf) continue; <<
2052 <<
2053 G4int nv = facet->GetNumberOfVertices(); <<
2054 base.resize(nv); <<
2055 for (G4int k=0; k<nv; ++k) { base[k] = fa <<
2056 <<
2057 G4double emin,emax; <<
2058 G4BoundingEnvelope benv(pyramid); <<
2059 if (!benv.CalculateExtent(pAxis,pVoxelLim <<
2060 if (emin < pMin) pMin = emin; <<
2061 if (emax > pMax) pMax = emax; <<
2062 if (eminlim > pMin && emaxlim < pMax) bre <<
2063 } 681 }
2064 return (pMin < pMax); << 682 return false;
2065 #endif <<
2066 } 683 }
2067 684
2068 ///////////////////////////////////////////// 685 ///////////////////////////////////////////////////////////////////////////////
2069 // 686 //
2070 G4double G4TessellatedSolid::GetMinXExtent () 687 G4double G4TessellatedSolid::GetMinXExtent () const
2071 { << 688 {return xMinExtent;}
2072 return fMinExtent.x(); <<
2073 } <<
2074 689
2075 ///////////////////////////////////////////// 690 ///////////////////////////////////////////////////////////////////////////////
2076 // 691 //
2077 G4double G4TessellatedSolid::GetMaxXExtent () 692 G4double G4TessellatedSolid::GetMaxXExtent () const
2078 { << 693 {return xMaxExtent;}
2079 return fMaxExtent.x(); <<
2080 } <<
2081 694
2082 ///////////////////////////////////////////// 695 ///////////////////////////////////////////////////////////////////////////////
2083 // 696 //
2084 G4double G4TessellatedSolid::GetMinYExtent () 697 G4double G4TessellatedSolid::GetMinYExtent () const
2085 { << 698 {return yMinExtent;}
2086 return fMinExtent.y(); <<
2087 } <<
2088 699
2089 ///////////////////////////////////////////// 700 ///////////////////////////////////////////////////////////////////////////////
2090 // 701 //
2091 G4double G4TessellatedSolid::GetMaxYExtent () 702 G4double G4TessellatedSolid::GetMaxYExtent () const
2092 { << 703 {return yMaxExtent;}
2093 return fMaxExtent.y(); <<
2094 } <<
2095 704
2096 ///////////////////////////////////////////// 705 ///////////////////////////////////////////////////////////////////////////////
2097 // 706 //
2098 G4double G4TessellatedSolid::GetMinZExtent () 707 G4double G4TessellatedSolid::GetMinZExtent () const
2099 { << 708 {return zMinExtent;}
2100 return fMinExtent.z(); <<
2101 } <<
2102 709
2103 ///////////////////////////////////////////// 710 ///////////////////////////////////////////////////////////////////////////////
2104 // 711 //
2105 G4double G4TessellatedSolid::GetMaxZExtent () 712 G4double G4TessellatedSolid::GetMaxZExtent () const
2106 { << 713 {return zMaxExtent;}
2107 return fMaxExtent.z(); <<
2108 } <<
2109 714
2110 ///////////////////////////////////////////// 715 ///////////////////////////////////////////////////////////////////////////////
2111 // 716 //
2112 G4VisExtent G4TessellatedSolid::GetExtent () 717 G4VisExtent G4TessellatedSolid::GetExtent () const
2113 { 718 {
2114 return { fMinExtent.x(), fMaxExtent.x(), << 719 return G4VisExtent (xMinExtent, xMaxExtent, yMinExtent, yMaxExtent,
2115 fMinExtent.y(), fMaxExtent.y(), << 720 zMinExtent, zMaxExtent);
2116 fMinExtent.z(), fMaxExtent.z() }; <<
2117 } <<
2118 <<
2119 ///////////////////////////////////////////// <<
2120 // <<
2121 G4double G4TessellatedSolid::GetCubicVolume ( <<
2122 { <<
2123 if (fCubicVolume != 0.) return fCubicVolume <<
2124 <<
2125 // For explanation of the following algorit <<
2126 // https://en.wikipedia.org/wiki/Polyhedron <<
2127 // http://wwwf.imperial.ac.uk/~rn/centroid. <<
2128 <<
2129 std::size_t size = fFacets.size(); <<
2130 for (std::size_t i = 0; i < size; ++i) <<
2131 { <<
2132 G4VFacet &facet = *fFacets[i]; <<
2133 G4double area = facet.GetArea(); <<
2134 G4ThreeVector unit_normal = facet.GetSurf <<
2135 fCubicVolume += area * (facet.GetVertex(0 <<
2136 } <<
2137 fCubicVolume /= 3.; <<
2138 return fCubicVolume; <<
2139 } <<
2140 <<
2141 ///////////////////////////////////////////// <<
2142 // <<
2143 G4double G4TessellatedSolid::GetSurfaceArea ( <<
2144 { <<
2145 if (fSurfaceArea != 0.) return fSurfaceArea <<
2146 <<
2147 std::size_t size = fFacets.size(); <<
2148 for (std::size_t i = 0; i < size; ++i) <<
2149 { <<
2150 G4VFacet &facet = *fFacets[i]; <<
2151 fSurfaceArea += facet.GetArea(); <<
2152 } <<
2153 return fSurfaceArea; <<
2154 } <<
2155 <<
2156 ///////////////////////////////////////////// <<
2157 // <<
2158 G4ThreeVector G4TessellatedSolid::GetPointOnS <<
2159 { <<
2160 // Select randomly a facet and return a ran <<
2161 <<
2162 auto i = (G4int) G4RandFlat::shoot(0., fFac <<
2163 return fFacets[i]->GetPointOnFace(); <<
2164 } <<
2165 <<
2166 ///////////////////////////////////////////// <<
2167 // <<
2168 // SetRandomVectorSet <<
2169 // <<
2170 // This is a set of predefined random vectors <<
2171 // in terms!) used to generate rays from a us <<
2172 // function Inside uses these to determine wh <<
2173 // outside of the tessellated solid. All vec <<
2174 // <<
2175 void G4TessellatedSolid::SetRandomVectors () <<
2176 { <<
2177 fRandir.resize(20); <<
2178 fRandir[0] = <<
2179 G4ThreeVector(-0.9577428892113370, 0.2732 <<
2180 fRandir[1] = <<
2181 G4ThreeVector(-0.8331264504940770,-0.5162 <<
2182 fRandir[2] = <<
2183 G4ThreeVector(-0.1516671651108820, 0.9666 <<
2184 fRandir[3] = <<
2185 G4ThreeVector( 0.6570250350323190,-0.6944 <<
2186 fRandir[4] = <<
2187 G4ThreeVector(-0.4820456281280320,-0.6331 <<
2188 fRandir[5] = <<
2189 G4ThreeVector( 0.7629032554236800 , 0.101 <<
2190 fRandir[6] = <<
2191 G4ThreeVector( 0.7689540409061150, 0.5034 <<
2192 fRandir[7] = <<
2193 G4ThreeVector( 0.5765188359255740, 0.5997 <<
2194 fRandir[8] = <<
2195 G4ThreeVector( 0.6660632777862070,-0.6362 <<
2196 fRandir[9] = <<
2197 G4ThreeVector( 0.3824415020414780, 0.6541 <<
2198 fRandir[10] = <<
2199 G4ThreeVector(-0.5107726564526760, 0.6020 <<
2200 fRandir[11] = <<
2201 G4ThreeVector( 0.7459135439578050, 0.6618 <<
2202 fRandir[12] = <<
2203 G4ThreeVector( 0.1536405855311580, 0.8117 <<
2204 fRandir[13] = <<
2205 G4ThreeVector( 0.0744395301705579,-0.8707 <<
2206 fRandir[14] = <<
2207 G4ThreeVector(-0.1665874645185400, 0.6018 <<
2208 fRandir[15] = <<
2209 G4ThreeVector( 0.7766902003633100, 0.6014 <<
2210 fRandir[16] = <<
2211 G4ThreeVector(-0.8710128685847430,-0.1434 <<
2212 fRandir[17] = <<
2213 G4ThreeVector( 0.8901082092766820,-0.4388 <<
2214 fRandir[18] = <<
2215 G4ThreeVector(-0.6430417431544370,-0.3295 <<
2216 fRandir[19] = <<
2217 G4ThreeVector( 0.6331124368380410, 0.6306 <<
2218 <<
2219 fMaxTries = 20; <<
2220 } 721 }
2221 722
2222 ///////////////////////////////////////////// 723 ///////////////////////////////////////////////////////////////////////////////
2223 // 724 //
2224 G4int G4TessellatedSolid::AllocatedMemoryWith << 725 //G4double G4TessellatedSolid::GetCubicVolume ()
2225 { << 726 // {return cubicVolume;}
2226 G4int base = sizeof(*this); << 727
2227 base += fVertexList.capacity() * sizeof(G4T <<
2228 base += fRandir.capacity() * sizeof(G4Three <<
2229 <<
2230 std::size_t limit = fFacets.size(); <<
2231 for (std::size_t i = 0; i < limit; ++i) <<
2232 { <<
2233 G4VFacet& facet = *fFacets[i]; <<
2234 base += facet.AllocatedMemory(); <<
2235 } <<
2236 <<
2237 for (const auto & fExtremeFacet : fExtremeF <<
2238 { <<
2239 G4VFacet &facet = *fExtremeFacet; <<
2240 base += facet.AllocatedMemory(); <<
2241 } <<
2242 return base; <<
2243 } <<
2244 <<
2245 ///////////////////////////////////////////// <<
2246 // <<
2247 G4int G4TessellatedSolid::AllocatedMemory() <<
2248 { <<
2249 G4int size = AllocatedMemoryWithoutVoxels() <<
2250 G4int sizeInsides = fInsides.GetNbytes(); <<
2251 G4int sizeVoxels = fVoxels.AllocatedMemory( <<
2252 size += sizeInsides + sizeVoxels; <<
2253 return size; <<
2254 } <<
2255 <<
2256 #endif <<
2257 728