| Geant4 Cross Reference (Editor's cut) |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration and of QinetiQ Ltd, *
20 // * subject to DEFCON 705 IPR conditions. *
21 // * By using, copying, modifying or distributing the software (or *
22 // * any work based on the software) you agree to acknowledge its *
23 // * use in resulting scientific publications, and indicate your *
24 // * acceptance of all terms of the Geant4 Software license. *
25 // ********************************************************************
26 //
27 // G4TessellatedSolid implementation
28 //
29 // 31.10.2004, P R Truscott, QinetiQ Ltd, UK - Created.
30 // 17.09.2007, P R Truscott, QinetiQ Ltd & Richard Holmberg
31 // Updated extensively prior to this date to deal with
32 // concaved tessellated surfaces, based on the algorithm
33 // of Richard Holmberg. This had been slightly modified
34 // to determine with inside the geometry by projecting
35 // random rays from the point provided. Now random rays
36 // are predefined rather than making use of random
37 // number generator at run-time.
38 // 12.10.2012, M Gayer, CERN, complete rewrite reducing memory
39 // requirements more than 50% and speedup by a factor of
40 // tens or more depending on the number of facets, thanks
41 // to voxelization of surface and improvements.
42 // Speedup factor of thousands for solids with number of
43 // facets in hundreds of thousands.
44 // 23.10.2016, E Tcherniaev, reimplemented CalculateExtent() to make
45 // use of G4BoundingEnvelope.
46 // --------------------------------------------------------------------
47
48 #include "G4TessellatedSolid.hh"
49
50 #if !defined(G4GEOM_USE_UTESSELLATEDSOLID)
51
52 #include <algorithm>
53 #include <fstream>
54 #include <iomanip>
55 #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_INITIALIZER;
77 }
78
79 using namespace std;
80
81 ///////////////////////////////////////////////////////////////////////////////
82 //
83 // Standard contructor has blank name and defines no fFacets.
84 //
85 G4TessellatedSolid::G4TessellatedSolid () : G4VSolid("dummy")
86 {
87 Initialize();
88 }
89
90 ///////////////////////////////////////////////////////////////////////////////
91 //
92 // Alternative constructor. Simple define name and geometry type - no fFacets
93 // to detine.
94 //
95 G4TessellatedSolid::G4TessellatedSolid (const G4String& name)
96 : G4VSolid(name)
97 {
98 Initialize();
99 }
100
101 ///////////////////////////////////////////////////////////////////////////////
102 //
103 // Fake default constructor - sets only member data and allocates memory
104 // for usage restricted to object persistency.
105 //
106 G4TessellatedSolid::G4TessellatedSolid( __void__& a) : G4VSolid(a)
107 {
108 Initialize();
109 fMinExtent.set(0,0,0);
110 fMaxExtent.set(0,0,0);
111 }
112
113
114 ///////////////////////////////////////////////////////////////////////////////
115 G4TessellatedSolid::~G4TessellatedSolid()
116 {
117 DeleteObjects();
118 }
119
120 ///////////////////////////////////////////////////////////////////////////////
121 //
122 // Copy constructor.
123 //
124 G4TessellatedSolid::G4TessellatedSolid (const G4TessellatedSolid& ts)
125 : G4VSolid(ts)
126 {
127 Initialize();
128
129 CopyObjects(ts);
130 }
131
132 ///////////////////////////////////////////////////////////////////////////////
133 //
134 // Assignment operator.
135 //
136 G4TessellatedSolid&
137 G4TessellatedSolid::operator= (const G4TessellatedSolid &ts)
138 {
139 if (&ts == this) return *this;
140
141 // Copy base class data
142 G4VSolid::operator=(ts);
143
144 DeleteObjects ();
145
146 Initialize();
147
148 CopyObjects (ts);
149
150 return *this;
151 }
152
153 ///////////////////////////////////////////////////////////////////////////////
154 //
155 void G4TessellatedSolid::Initialize()
156 {
157 kCarToleranceHalf = 0.5*kCarTolerance;
158
159 fRebuildPolyhedron = false; fpPolyhedron = nullptr;
160 fCubicVolume = 0.; fSurfaceArea = 0.;
161
162 fGeometryType = "G4TessellatedSolid";
163 fSolidClosed = false;
164
165 fMinExtent.set(kInfinity,kInfinity,kInfinity);
166 fMaxExtent.set(-kInfinity,-kInfinity,-kInfinity);
167
168 SetRandomVectors();
169 }
170
171 ///////////////////////////////////////////////////////////////////////////////
172 //
173 void G4TessellatedSolid::DeleteObjects()
174 {
175 std::size_t size = fFacets.size();
176 for (std::size_t i = 0; i < size; ++i) { delete fFacets[i]; }
177 fFacets.clear();
178 delete fpPolyhedron; fpPolyhedron = nullptr;
179 }
180
181 ///////////////////////////////////////////////////////////////////////////////
182 //
183 void G4TessellatedSolid::CopyObjects (const G4TessellatedSolid &ts)
184 {
185 G4ThreeVector reductionRatio;
186 G4int fmaxVoxels = fVoxels.GetMaxVoxels(reductionRatio);
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 {
195 G4VFacet *facetClone = (ts.GetFacet(i))->GetClone();
196 AddFacet(facetClone);
197 }
198 if (ts.GetSolidClosed()) SetSolidClosed(true);
199 }
200
201 ///////////////////////////////////////////////////////////////////////////////
202 //
203 // Add a facet to the facet list.
204 // Note that you can add, but you cannot delete.
205 //
206 G4bool G4TessellatedSolid::AddFacet (G4VFacet* aFacet)
207 {
208 // Add the facet to the vector.
209 //
210 if (fSolidClosed)
211 {
212 G4Exception("G4TessellatedSolid::AddFacet()", "GeomSolids1002",
213 JustWarning, "Attempt to add facets when solid is closed.");
214 return false;
215 }
216 else if (aFacet->IsDefined())
217 {
218 set<G4VertexInfo,G4VertexComparator>::iterator begin
219 = fFacetList.begin(), end = fFacetList.end(), pos, it;
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 {
228 G4double kCarTolerance3 = 3 * kCarTolerance;
229 pos = fFacetList.lower_bound(value);
230
231 it = pos;
232 while (!found && it != end) // Loop checking, 13.08.2015, G.Cosmo
233 {
234 G4int id = (*it).id;
235 G4VFacet *facet = fFacets[id];
236 G4ThreeVector q = facet->GetCircumcentre();
237 if ((found = (facet == aFacet))) break;
238 G4double dif = q.x() + q.y() + q.z() - value.mag2;
239 if (dif > kCarTolerance3) break;
240 it++;
241 }
242
243 if (fFacets.size() > 1)
244 {
245 it = pos;
246 while (!found && it != begin) // Loop checking, 13.08.2015, G.Cosmo
247 {
248 --it;
249 G4int id = (*it).id;
250 G4VFacet *facet = fFacets[id];
251 G4ThreeVector q = facet->GetCircumcentre();
252 found = (facet == aFacet);
253 if (found) break;
254 G4double dif = value.mag2 - (q.x() + q.y() + q.z());
255 if (dif > kCarTolerance3) break;
256 }
257 }
258 }
259
260 if (!found)
261 {
262 fFacets.push_back(aFacet);
263 fFacetList.insert(value);
264 }
265 return true;
266 }
267 else
268 {
269 G4Exception("G4TessellatedSolid::AddFacet()", "GeomSolids1002",
270 JustWarning, "Attempt to add facet not properly defined.");
271 aFacet->StreamInfo(G4cout);
272 return false;
273 }
274 }
275
276 ///////////////////////////////////////////////////////////////////////////////
277 //
278 G4int G4TessellatedSolid::SetAllUsingStack(const std::vector<G4int>& voxel,
279 const std::vector<G4int>& max,
280 G4bool status, G4SurfBits& checked)
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.08.2015, G.Cosmo
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(i).size();
333 G4int size = maxVoxels[0] * maxVoxels[1] * maxVoxels[2];
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] - 1; ++voxel[2])
341 {
342 for (voxel[1] = 0; voxel[1] < maxVoxels[1] - 1; ++voxel[1])
343 {
344 for (voxel[0] = 0; voxel[0] < maxVoxels[0] - 1; ++voxel[0])
345 {
346 G4int index = fVoxels.GetVoxelsIndex(voxel);
347 if (!checked[index] && fVoxels.IsEmpty(index))
348 {
349 for (auto i = 0; i <= 2; ++i)
350 {
351 point[i] = fVoxels.GetBoundary(i)[voxel[i]];
352 }
353 auto inside = (G4bool) (InsideNoVoxels(point) == kInside);
354 SetAllUsingStack(voxel, maxVoxels, inside, checked);
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..." << G4endl;
375 #endif
376 PrecalculateInsides();
377 }
378 }
379
380 ///////////////////////////////////////////////////////////////////////////////
381 //
382 // Compute extremeFacets, i.e. find those facets that have surface
383 // planes that bound the volume.
384 // Note that this is going to reject concaved surfaces as being extreme. Also
385 // note that if the vertex is on the facet, displacement is zero, so IsInside
386 // returns true. So will this work?? Need non-equality
387 // "G4bool inside = displacement < 0.0;"
388 // or
389 // "G4bool inside = displacement <= -0.5*kCarTolerance"
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) { vertices[i] = fVertexList[i]; }
398
399 // Shuffle vertices
400 std::mt19937 gen(12345678);
401 std::shuffle(vertices.begin(), vertices.end(), gen);
402
403 // Select six extreme vertices in different directions
404 G4ThreeVector points[6];
405 for (auto & point : points) { point = vertices[0]; }
406 for (std::size_t i=1; i < vsize; ++i)
407 {
408 if (vertices[i].x() < points[0].x()) points[0] = vertices[i];
409 if (vertices[i].x() > points[1].x()) points[1] = vertices[i];
410 if (vertices[i].y() < points[2].y()) points[2] = vertices[i];
411 if (vertices[i].y() > points[3].y()) points[3] = vertices[i];
412 if (vertices[i].z() < points[4].z()) points[4] = vertices[i];
413 if (vertices[i].z() > points[5].z()) points[5] = vertices[i];
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(&facet);
441 }
442 }
443
444 ///////////////////////////////////////////////////////////////////////////////
445 //
446 void G4TessellatedSolid::CreateVertexList()
447 {
448 // The algorithm:
449 // we will have additional vertexListSorted, where all the items will be
450 // sorted by magnitude of vertice vector.
451 // New candidate for fVertexList - we will determine the position fo first
452 // item which would be within its magnitude - 0.5*kCarTolerance.
453 // We will go trough until we will reach > +0.5 kCarTolerance.
454 // Comparison (q-p).mag() < 0.5*kCarTolerance will be made.
455 // They can be just stored in std::vector, with custom insertion based
456 // on binary search.
457
458 set<G4VertexInfo,G4VertexComparator> vertexListSorted;
459 set<G4VertexInfo,G4VertexComparator>::iterator begin
460 = vertexListSorted.begin(), end = vertexListSorted.end(), pos, it;
461 G4ThreeVector p;
462 G4VertexInfo value;
463
464 fVertexList.clear();
465 std::size_t size = fFacets.size();
466
467 G4double kCarTolerance24 = kCarTolerance * kCarTolerance / 4.0;
468 G4double kCarTolerance3 = 3 * kCarTolerance;
469 vector<G4int> newIndex(100);
470
471 for (std::size_t k = 0; k < size; ++k)
472 {
473 G4VFacet &facet = *fFacets[k];
474 G4int max = facet.GetNumberOfVertices();
475
476 for (G4int i = 0; i < max; ++i)
477 {
478 p = facet.GetVertex(i);
479 value.id = (G4int)fVertexList.size();
480 value.mag2 = p.x() + p.y() + p.z();
481
482 G4bool found = false;
483 G4int id = 0;
484 if (!OutsideOfExtent(p, kCarTolerance))
485 {
486 pos = vertexListSorted.lower_bound(value);
487 it = pos;
488 while (it != end) // Loop checking, 13.08.2015, G.Cosmo
489 {
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.mag2;
496 if (dif > kCarTolerance3) break;
497 ++it;
498 }
499
500 if (!found && (fVertexList.size() > 1))
501 {
502 it = pos;
503 while (it != begin) // Loop checking, 13.08.2015, G.Cosmo
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() + q.z());
512 if (dif > kCarTolerance3) break;
513 }
514 }
515 }
516
517 if (!found)
518 {
519 #ifdef G4SPECSDEBUG
520 G4cout << p.x() << ":" << p.y() << ":" << p.z() << G4endl;
521 G4cout << "Adding new vertex #" << i << " of facet " << k
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 limits of the volume.
532 //
533 if (value.id == 0) fMinExtent = fMaxExtent = p;
534 else
535 {
536 if (p.x() > fMaxExtent.x()) fMaxExtent.setX(p.x());
537 else if (p.x() < fMinExtent.x()) fMinExtent.setX(p.x());
538 if (p.y() > fMaxExtent.y()) fMaxExtent.setY(p.y());
539 else if (p.y() < fMinExtent.y()) fMinExtent.setY(p.y());
540 if (p.z() > fMaxExtent.z()) fMaxExtent.setZ(p.z());
541 else if (p.z() < fMinExtent.z()) fMinExtent.setZ(p.z());
542 }
543 }
544 else
545 {
546 #ifdef G4SPECSDEBUG
547 G4cout << p.x() << ":" << p.y() << ":" << p.z() << G4endl;
548 G4cout << "Vertex #" << i << " of facet " << k
549 << " found, redirecting to " << id << G4endl;
550 G4cout << "===" << G4endl;
551 #endif
552 newIndex[i] = id;
553 }
554 }
555 // only now it is possible to change vertices pointer
556 //
557 facet.SetVertices(&fVertexList);
558 for (G4int i = 0; i < max; ++i)
559 facet.SetVertexIndex(i,newIndex[i]);
560 }
561 vector<G4ThreeVector>(fVertexList).swap(fVertexList);
562
563 #ifdef G4SPECSDEBUG
564 G4double previousValue = 0.;
565 for (auto res=vertexListSorted.cbegin(); res!=vertexListSorted.cend(); ++res)
566 {
567 G4int id = (*res).id;
568 G4ThreeVector vec = fVertexList[id];
569 G4double mvalue = vec.x() + vec.y() + vec.z();
570 if (previousValue && (previousValue - 1e-9 > mvalue))
571 G4cout << "Error in CreateVertexList: previousValue " << previousValue
572 << " is smaller than mvalue " << mvalue << G4endl;
573 previousValue = mvalue;
574 }
575 #endif
576 }
577
578 ///////////////////////////////////////////////////////////////////////////////
579 //
580 void G4TessellatedSolid::DisplayAllocatedMemory()
581 {
582 G4int without = AllocatedMemoryWithoutVoxels();
583 G4int with = AllocatedMemory();
584 G4double ratio = (G4double) with / without;
585 G4cout << "G4TessellatedSolid - Allocated memory without voxel overhead "
586 << without << "; with " << with << "; ratio: " << ratio << G4endl;
587 }
588
589 ///////////////////////////////////////////////////////////////////////////////
590 //
591 void G4TessellatedSolid::SetSolidClosed (const G4bool t)
592 {
593 if (t)
594 {
595 #ifdef G4SPECSDEBUG
596 G4cout << "Creating vertex list..." << G4endl;
597 #endif
598 CreateVertexList();
599
600 #ifdef G4SPECSDEBUG
601 G4cout << "Setting extreme facets..." << G4endl;
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..." << G4endl;
616 #endif
617 G4int irep = CheckStructure();
618 if (irep != 0)
619 {
620 if ((irep & 1) != 0)
621 {
622 std::ostringstream message;
623 message << "Defects in solid: " << GetName()
624 << " - negative cubic volume, please check orientation of facets!";
625 G4Exception("G4TessellatedSolid::SetSolidClosed()",
626 "GeomSolids1001", JustWarning, message);
627 }
628 if ((irep & 2) != 0)
629 {
630 std::ostringstream message;
631 message << "Defects in solid: " << GetName()
632 << " - some facets have wrong orientation!";
633 G4Exception("G4TessellatedSolid::SetSolidClosed()",
634 "GeomSolids1001", JustWarning, message);
635 }
636 if ((irep & 4) != 0)
637 {
638 std::ostringstream message;
639 message << "Defects in solid: " << GetName()
640 << " - there are holes in the surface!";
641 G4Exception("G4TessellatedSolid::SetSolidClosed()",
642 "GeomSolids1001", JustWarning, message);
643 }
644 }
645 }
646 fSolidClosed = t;
647 }
648
649 ///////////////////////////////////////////////////////////////////////////////
650 //
651 // GetSolidClosed
652 //
653 // Used to determine whether the solid is closed to adding further fFacets.
654 //
655 G4bool G4TessellatedSolid::GetSolidClosed () const
656 {
657 return fSolidClosed;
658 }
659
660 ///////////////////////////////////////////////////////////////////////////////
661 //
662 // CheckStructure
663 //
664 // Checks structure of the solid. Return value is a sum of the following
665 // defect indicators, if any (0 means no defects):
666 // 1 - cubic volume is negative, wrong orientation of facets
667 // 2 - some facets have wrong orientation
668 // 4 - holes in the surface
669 //
670 G4int G4TessellatedSolid::CheckStructure() const
671 {
672 G4int nedge = 0;
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(0).dot(facet.GetSurfaceNormal()));
683 }
684 G4int ivolume = static_cast<G4int>(volume <= 0.);
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 == 0) ? nnode - 1 : k - 1);
697 int64_t i2 = facet.GetVertexIndex(k);
698 int64_t inverse = static_cast<int64_t>(i2 > i1);
699 if (inverse != 0) std::swap(i1, i2);
700 iedge[kk++] = i1*1000000000 + i2*2 + inverse;
701 }
702 }
703 std::sort(iedge.begin(), iedge.end());
704
705 // Check edges, correct structure should consist of paired edges
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) // paired edges with different orientation
714 {
715 i += 2;
716 }
717 else if (iedge[i + 1] == iedge[i]) // paired edges with the same orientation
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 tessellated solids together, so
736 // that the solid on the left then includes all of the facets in the solid
737 // on the right. Note that copies of the facets are generated, rather than
738 // using the original facet set of the solid on the right.
739 //
740 G4TessellatedSolid&
741 G4TessellatedSolid::operator+=(const G4TessellatedSolid& right)
742 {
743 G4int size = right.GetNumberOfFacets();
744 for (G4int i = 0; i < size; ++i)
745 AddFacet(right.GetFacet(i)->GetClone());
746
747 return *this;
748 }
749
750 ///////////////////////////////////////////////////////////////////////////////
751 //
752 // GetNumberOfFacets
753 //
754 G4int G4TessellatedSolid::GetNumberOfFacets() const
755 {
756 return (G4int)fFacets.size();
757 }
758
759 ///////////////////////////////////////////////////////////////////////////////
760 //
761 EInside G4TessellatedSolid::InsideVoxels(const G4ThreeVector& p) const
762 {
763 //
764 // First the simple test - check if we're outside of the X-Y-Z extremes
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] ? kInside : kOutside;
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 adaptation of the method implemented by
798 // Rickard Holmberg augmented with information from Schneider & Eberly,
799 // "Geometric Tools for Computer Graphics," pp700-701, 2003. In essence,
800 // we're trying to determine whether we're inside the volume by projecting
801 // a few rays and determining if the first surface crossed is has a normal
802 // vector between 0 to pi/2 (out-going) or pi/2 to pi (in-going).
803 // We should also avoid rays which are nearly within the plane of the
804 // tessellated surface, and therefore produce rays randomly.
805 // For the moment, this is a bit over-engineered (belt-braces-and-ducttape).
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 the vector is not nearly parallel
823 // to the surface of any facet since this causes ambiguities. The usual
824 // case is that the angles should be sufficiently different, but there
825 // are 20 random directions to select from - hopefully sufficient.
826 //
827 distOut = distIn = kInfinity;
828 const G4ThreeVector& v = fRandir[sm];
829 ++sm;
830 //
831 // This code could be voxelized by the same algorithm, which is used for
832 // DistanceToOut(). We will traverse through fVoxels. we will call
833 // intersect only for those, which would be candidates and was not
834 // checked before.
835 //
836 G4ThreeVector currentPoint = p;
837 G4ThreeVector direction = v.unit();
838 // G4SurfBits exclusion(fVoxels.GetBitsPerSlice());
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.Cosmo
847 {
848 const vector<G4int> &candidates =
849 started ? startingCandidates : fVoxels.GetCandidates(curVoxel);
850 started = false;
851 if (auto candidatesCount = (G4int)candidates.size())
852 {
853 for (G4int i = 0 ; i < candidatesCount; ++i)
854 {
855 G4int candidate = candidates[i];
856 // bits.SetBitNumber(candidate);
857 G4VFacet& facet = *fFacets[candidate];
858
859 crossingO = facet.Intersect(p,v,true,distO,distFromSurfaceO,normalO);
860 crossingI = facet.Intersect(p,v,false,distI,distFromSurfaceI,normalI);
861
862 if (crossingO || crossingI)
863 {
864 crossed = true;
865
866 nearParallel = (crossingO
867 && std::fabs(normalO.dot(v))<dirTolerance)
868 || (crossingI && std::fabs(normalI.dot(v))<dirTolerance);
869 if (!nearParallel)
870 {
871 if (crossingO && distO > 0.0 && distO < distOut)
872 distOut = distO;
873 if (crossingI && distI > 0.0 && distI < distIn)
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(curVoxel);
886 G4bool inside = fInsides[index];
887 location = inside ? kInside : kOutside;
888 return location;
889 }
890 }
891
892 G4double shift=fVoxels.DistanceToNext(currentPoint, direction, curVoxel);
893 if (shift == kInfinity) break;
894
895 currentPoint += direction * (shift + shiftBonus);
896 }
897 while (fVoxels.UpdateCurrentVoxel(currentPoint, direction, curVoxel));
898
899 }
900 while (nearParallel && sm != fMaxTries);
901 //
902 // Here we loop through the facets to find out if there is an intersection
903 // between the ray and that facet. The test if performed separately whether
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 directions. If nTries is set sufficiently
911 // low (nTries <= 0.5*maxTries) then this would indicate that there is
912 // something wrong with geometry.
913 //
914 std::ostringstream message;
915 G4long oldprc = message.precision(16);
916 message << "Cannot determine whether point is inside or outside volume!"
917 << G4endl
918 << "Solid name = " << GetName() << G4endl
919 << "Geometry Type = " << fGeometryType << G4endl
920 << "Number of facets = " << fFacets.size() << G4endl
921 << "Position:" << G4endl << G4endl
922 << "p.x() = " << p.x()/mm << " mm" << G4endl
923 << "p.y() = " << p.y()/mm << " mm" << G4endl
924 << "p.z() = " << p.z()/mm << " mm";
925 message.precision(oldprc);
926 G4Exception("G4TessellatedSolid::Inside()",
927 "GeomSolids1002", JustWarning, message);
928 }
929 #endif
930
931 // In the next if-then-elseif G4String the logic is as follows:
932 // (1) You don't hit anything so cannot be inside volume, provided volume
933 // constructed correctly!
934 // (2) Distance to inside (ie. nearest facet such that you enter facet) is
935 // shorter than distance to outside (nearest facet such that you exit
936 // facet) - on condition of safety distance - therefore we're outside.
937 // (3) Distance to outside is shorter than distance to inside therefore
938 // we're inside.
939 //
940 if (distIn == kInfinity && distOut == kInfinity)
941 location = kOutside;
942 else if (distIn <= distOut - kCarToleranceHalf)
943 location = kOutside;
944 else if (distOut <= distIn - kCarToleranceHalf)
945 location = kInside;
946
947 return location;
948 }
949
950 ///////////////////////////////////////////////////////////////////////////////
951 //
952 EInside G4TessellatedSolid::InsideNoVoxels (const G4ThreeVector &p) const
953 {
954 //
955 // First the simple test - check if we're outside of the X-Y-Z extremes
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 adaptation of the method implemented by
980 // Rickard Holmberg augmented with information from Schneider & Eberly,
981 // "Geometric Tools for Computer Graphics," pp700-701, 2003. In essence, we're
982 // trying to determine whether we're inside the volume by projecting a few
983 // rays and determining if the first surface crossed is has a normal vector
984 // between 0 to pi/2 (out-going) or pi/2 to pi (in-going). We should also
985 // avoid rays which are nearly within the plane of the tessellated surface,
986 // and therefore produce rays randomly. For the moment, this is a bit
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.Cosmo
1012 {
1013 //
1014 // We loop until we find direction where the vector is not nearly parallel
1015 // to the surface of any facet since this causes ambiguities. The usual
1016 // case is that the angles should be sufficiently different, but there
1017 // are 20 random directions to select from - hopefully sufficient.
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.Cosmo
1025 {
1026 //
1027 // Here we loop through the facets to find out if there is an
1028 // intersection between the ray and that facet. The test if performed
1029 // separately whether the ray is entering the facet or exiting.
1030 //
1031 crossingO = ((*f)->Intersect(p,v,true,distO,distFromSurfaceO,normalO));
1032 crossingI = ((*f)->Intersect(p,v,false,distI,distFromSurfaceI,normalI));
1033 if (crossingO || crossingI)
1034 {
1035 nearParallel = (crossingO && std::fabs(normalO.dot(v))<dirTolerance)
1036 || (crossingI && std::fabs(normalI.dot(v))<dirTolerance);
1037 if (!nearParallel)
1038 {
1039 if (crossingO && distO > 0.0 && distO < distOut) distOut = distO;
1040 if (crossingI && distI > 0.0 && distI < distIn) distIn = distI;
1041 }
1042 }
1043 } while (!nearParallel && ++f != fFacets.cend());
1044 } while (nearParallel && sm != fMaxTries);
1045
1046 #ifdef G4VERBOSE
1047 if (sm == fMaxTries)
1048 {
1049 //
1050 // We've run out of random vector directions. If nTries is set
1051 // sufficiently low (nTries <= 0.5*maxTries) then this would indicate
1052 // that there is something wrong with geometry.
1053 //
1054 std::ostringstream message;
1055 G4long oldprc = message.precision(16);
1056 message << "Cannot determine whether point is inside or outside volume!"
1057 << G4endl
1058 << "Solid name = " << GetName() << G4endl
1059 << "Geometry Type = " << fGeometryType << G4endl
1060 << "Number of facets = " << fFacets.size() << G4endl
1061 << "Position:" << G4endl << G4endl
1062 << "p.x() = " << p.x()/mm << " mm" << G4endl
1063 << "p.y() = " << p.y()/mm << " mm" << G4endl
1064 << "p.z() = " << p.z()/mm << " mm";
1065 message.precision(oldprc);
1066 G4Exception("G4TessellatedSolid::Inside()",
1067 "GeomSolids1002", JustWarning, message);
1068 }
1069 #endif
1070 //
1071 // In the next if-then-elseif G4String the logic is as follows:
1072 // (1) You don't hit anything so cannot be inside volume, provided volume
1073 // constructed correctly!
1074 // (2) Distance to inside (ie. nearest facet such that you enter facet) is
1075 // shorter than distance to outside (nearest facet such that you exit
1076 // facet) - on condition of safety distance - therefore we're outside.
1077 // (3) Distance to outside is shorter than distance to inside therefore
1078 // we're inside.
1079 //
1080 if (distIn == kInfinity && distOut == kInfinity)
1081 locationprime = kOutside;
1082 else if (distIn <= distOut - kCarToleranceHalf)
1083 locationprime = kOutside;
1084 else if (distOut <= distIn - kCarToleranceHalf)
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 point p, normally the point should
1096 // be located on the surface. Return -1 if no facet selected.
1097 //
1098 G4int G4TessellatedSolid::GetFacetIndex (const G4ThreeVector& p) const
1099 {
1100 G4int index = -1;
1101
1102 if (fVoxels.GetCountOfVoxels() > 1)
1103 {
1104 vector<G4int> curVoxel(3);
1105 fVoxels.GetVoxel(curVoxel, p);
1106 const vector<G4int> &candidates = fVoxels.GetCandidates(curVoxel);
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, minDist);
1115 if (dist <= kCarToleranceHalf) return index = candidate;
1116 if (dist < minDist)
1117 {
1118 minDist = dist;
1119 index = candidate;
1120 }
1121 }
1122 }
1123 }
1124 else
1125 {
1126 G4double minDist = kInfinity;
1127 std::size_t size = fFacets.size();
1128 for (std::size_t i = 0; i < size; ++i)
1129 {
1130 G4VFacet& facet = *fFacets[i];
1131 G4double dist = facet.Distance(p, minDist);
1132 if (dist < minDist)
1133 {
1134 minDist = dist;
1135 index = (G4int)i;
1136 }
1137 }
1138 }
1139 return index;
1140 }
1141
1142 ///////////////////////////////////////////////////////////////////////////////
1143 //
1144 // Return the outwards pointing unit normal of the shape for the
1145 // surface closest to the point at offset p.
1146 //
1147 G4bool G4TessellatedSolid::Normal (const G4ThreeVector& p,
1148 G4ThreeVector& aNormal) const
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.GetCandidates(curVoxel);
1158 // fVoxels.GetCandidatesVoxelArray(p, candidates, 0);
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,minDist);
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
1179 {
1180 minDist = kInfinity;
1181 std::size_t size = fFacets.size();
1182 for (std::size_t i = 0; i < size; ++i)
1183 {
1184 G4VFacet& f = *fFacets[i];
1185 G4double dist = f.Distance(p, minDist);
1186 if (dist < minDist)
1187 {
1188 minDist = dist;
1189 facet = &f;
1190 }
1191 }
1192 }
1193
1194 if (minDist != kInfinity)
1195 {
1196 if (facet != nullptr) { aNormal = facet->GetSurfaceNormal(); }
1197 return minDist <= kCarToleranceHalf;
1198 }
1199 else
1200 {
1201 #ifdef G4VERBOSE
1202 std::ostringstream message;
1203 message << "Point p is not on surface !?" << G4endl
1204 << " No facets found for point: " << p << " !" << G4endl
1205 << " Returning approximated value for normal.";
1206
1207 G4Exception("G4TessellatedSolid::SurfaceNormal(p)",
1208 "GeomSolids1002", JustWarning, message );
1209 #endif
1210 aNormal = (p.z() > 0 ? G4ThreeVector(0,0,1) : G4ThreeVector(0,0,-1));
1211 return false;
1212 }
1213 }
1214
1215 ///////////////////////////////////////////////////////////////////////////////
1216 //
1217 // G4double DistanceToIn(const G4ThreeVector& p, const G4ThreeVector& v)
1218 //
1219 // Return the distance along the normalised vector v to the shape,
1220 // from the point at offset p. If there is no intersection, return
1221 // kInfinity. The first intersection resulting from 'leaving' a
1222 // surface/volume is discarded. Hence, this is tolerant of points on
1223 // surface of shape.
1224 //
1225 G4double
1226 G4TessellatedSolid::DistanceToInNoVoxels (const G4ThreeVector& p,
1227 const G4ThreeVector& v,
1228 G4double /*aPstep*/) const
1229 {
1230 G4double minDist = kInfinity;
1231 G4double dist = 0.0;
1232 G4double distFromSurface = 0.0;
1233 G4ThreeVector normal;
1234
1235 #if G4SPECSDEBUG
1236 if (Inside(p) == kInside )
1237 {
1238 std::ostringstream message;
1239 G4int oldprc = message.precision(16) ;
1240 message << "Point p is already inside!?" << G4endl
1241 << "Position:" << G4endl << G4endl
1242 << " p.x() = " << p.x()/mm << " mm" << G4endl
1243 << " p.y() = " << p.y()/mm << " mm" << G4endl
1244 << " p.z() = " << p.z()/mm << " mm" << G4endl
1245 << "DistanceToOut(p) == " << DistanceToOut(p);
1246 message.precision(oldprc) ;
1247 G4Exception("G4TriangularFacet::DistanceToIn(p,v)",
1248 "GeomSolids1002", JustWarning, message);
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,distFromSurface,normal))
1257 {
1258 //
1259 // set minDist to the new distance to current facet if distFromSurface
1260 // is in positive direction and point is not at surface. If the point is
1261 // within 0.5*kCarTolerance of the surface, then force distance to be
1262 // zero and leave member function immediately (for efficiency), as
1263 // proposed by & credit to Akira Okumura.
1264 //
1265 if (distFromSurface > kCarToleranceHalf && dist >= 0.0 && dist < minDist)
1266 {
1267 minDist = dist;
1268 }
1269 else
1270 {
1271 if (-kCarToleranceHalf <= dist && dist <= kCarToleranceHalf)
1272 {
1273 return 0.0;
1274 }
1275 else
1276 {
1277 if (distFromSurface > -kCarToleranceHalf
1278 && distFromSurface < kCarToleranceHalf)
1279 {
1280 minDist = dist;
1281 }
1282 }
1283 }
1284 }
1285 }
1286 return minDist;
1287 }
1288
1289 ///////////////////////////////////////////////////////////////////////////////
1290 //
1291 G4double
1292 G4TessellatedSolid::DistanceToOutNoVoxels (const G4ThreeVector& p,
1293 const G4ThreeVector& v,
1294 G4ThreeVector& aNormalVector,
1295 G4bool& aConvex,
1296 G4double /*aPstep*/) const
1297 {
1298 G4double minDist = kInfinity;
1299 G4double dist = 0.0;
1300 G4double distFromSurface = 0.0;
1301 G4ThreeVector normal, minNormal;
1302
1303 #if G4SPECSDEBUG
1304 if ( Inside(p) == kOutside )
1305 {
1306 std::ostringstream message;
1307 G4int oldprc = message.precision(16) ;
1308 message << "Point p is already outside!?" << G4endl
1309 << "Position:" << G4endl << G4endl
1310 << " p.x() = " << p.x()/mm << " mm" << G4endl
1311 << " p.y() = " << p.y()/mm << " mm" << G4endl
1312 << " p.z() = " << p.z()/mm << " mm" << G4endl
1313 << "DistanceToIn(p) == " << DistanceToIn(p);
1314 message.precision(oldprc) ;
1315 G4Exception("G4TriangularFacet::DistanceToOut(p)",
1316 "GeomSolids1002", JustWarning, message);
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 {
1324 G4VFacet& facet = *fFacets[i];
1325 if (facet.Intersect(p,v,true,dist,distFromSurface,normal))
1326 {
1327 if (distFromSurface > 0.0 && distFromSurface <= kCarToleranceHalf
1328 && facet.Distance(p,kCarTolerance) <= kCarToleranceHalf)
1329 {
1330 // We are on a surface. Return zero.
1331 aConvex = (fExtremeFacets.find(&facet) != fExtremeFacets.end());
1332 // Normal(p, aNormalVector);
1333 // aNormalVector = facet.GetSurfaceNormal();
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(&facet) != fExtremeFacets.end());
1342 }
1343 }
1344 }
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<G4int>& candidates,
1364 const G4ThreeVector& aPoint,
1365 const G4ThreeVector& direction,
1366 G4double& minDist, G4ThreeVector& minNormal,
1367 G4int& minCandidate ) const
1368 {
1369 auto candidatesCount = (G4int)candidates.size();
1370 G4double dist = 0.0;
1371 G4double distFromSurface = 0.0;
1372 G4ThreeVector normal;
1373
1374 for (G4int i = 0 ; i < candidatesCount; ++i)
1375 {
1376 G4int candidate = candidates[i];
1377 G4VFacet& facet = *fFacets[candidate];
1378 if (facet.Intersect(aPoint,direction,true,dist,distFromSurface,normal))
1379 {
1380 if (distFromSurface > 0.0 && distFromSurface <= kCarToleranceHalf
1381 && facet.Distance(aPoint,kCarTolerance) <= kCarToleranceHalf)
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 }
1399
1400 ///////////////////////////////////////////////////////////////////////////////
1401 //
1402 G4double
1403 G4TessellatedSolid::DistanceToOutCore(const G4ThreeVector& aPoint,
1404 const G4ThreeVector& aDirection,
1405 G4ThreeVector& aNormalVector,
1406 G4bool &aConvex,
1407 G4double aPstep) const
1408 {
1409 G4double minDistance;
1410
1411 if (fVoxels.GetCountOfVoxels() > 1)
1412 {
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.Cosmo
1429 {
1430 const vector<G4int>& candidates = fVoxels.GetCandidates(curVoxel);
1431 if (old == &candidates)
1432 ++old;
1433 if (old != &candidates && !candidates.empty())
1434 {
1435 DistanceToOutCandidates(candidates, aPoint, direction, minDistance,
1436 aNormalVector, minCandidate);
1437 if (minDistance <= totalShift) break;
1438 }
1439
1440 G4double shift=fVoxels.DistanceToNext(currentPoint, direction, curVoxel);
1441 if (shift == kInfinity) break;
1442
1443 totalShift += shift;
1444 if (minDistance <= totalShift) break;
1445
1446 currentPoint += direction * (shift + shiftBonus);
1447
1448 old = &candidates;
1449 }
1450 while (fVoxels.UpdateCurrentVoxel(currentPoint, direction, curVoxel));
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[minCandidate])
1462 != fExtremeFacets.end());
1463 }
1464 }
1465 else
1466 {
1467 minDistance = DistanceToOutNoVoxels(aPoint, aDirection, aNormalVector,
1468 aConvex, aPstep);
1469 }
1470 return minDistance;
1471 }
1472
1473 ///////////////////////////////////////////////////////////////////////////////
1474 //
1475 G4double G4TessellatedSolid::
1476 DistanceToInCandidates(const std::vector<G4int>& candidates,
1477 const G4ThreeVector& aPoint,
1478 const G4ThreeVector& direction) const
1479 {
1480 auto candidatesCount = (G4int)candidates.size();
1481 G4double dist = 0.0;
1482 G4double distFromSurface = 0.0;
1483 G4ThreeVector normal;
1484
1485 G4double minDistance = kInfinity;
1486 for (G4int i = 0 ; i < candidatesCount; ++i)
1487 {
1488 G4int candidate = candidates[i];
1489 G4VFacet& facet = *fFacets[candidate];
1490 if (facet.Intersect(aPoint,direction,false,dist,distFromSurface,normal))
1491 {
1492 //
1493 // Set minDist to the new distance to current facet if distFromSurface is
1494 // in positive direction and point is not at surface. If the point is
1495 // within 0.5*kCarTolerance of the surface, then force distance to be
1496 // zero and leave member function immediately (for efficiency), as
1497 // proposed by & credit to Akira Okumura.
1498 //
1499 if ( (distFromSurface > kCarToleranceHalf)
1500 && (dist >= 0.0) && (dist < minDistance))
1501 {
1502 minDistance = dist;
1503 }
1504 else
1505 {
1506 if (-kCarToleranceHalf <= dist && dist <= kCarToleranceHalf)
1507 {
1508 return 0.0;
1509 }
1510 else if (distFromSurface > -kCarToleranceHalf
1511 && distFromSurface < kCarToleranceHalf)
1512 {
1513 minDistance = dist;
1514 }
1515 }
1516 }
1517 }
1518 return minDistance;
1519 }
1520
1521 ///////////////////////////////////////////////////////////////////////////////
1522 //
1523 G4double
1524 G4TessellatedSolid::DistanceToInCore(const G4ThreeVector& aPoint,
1525 const G4ThreeVector& aDirection,
1526 G4double aPstep) const
1527 {
1528 G4double minDistance;
1529
1530 if (fVoxels.GetCountOfVoxels() > 1)
1531 {
1532 minDistance = kInfinity;
1533 G4ThreeVector currentPoint = aPoint;
1534 G4ThreeVector direction = aDirection.unit();
1535 G4double shift = fVoxels.DistanceToFirst(currentPoint, direction);
1536 if (shift == kInfinity) return shift;
1537 G4double shiftBonus = kCarTolerance;
1538 if (shift != 0.0)
1539 currentPoint += direction * (shift + shiftBonus);
1540 // if (!fVoxels.Contains(currentPoint)) return minDistance;
1541 G4double totalShift = shift;
1542
1543 // G4SurfBits exclusion; // (1/*fVoxels.GetBitsPerSlice()*/);
1544 vector<G4int> curVoxel(3);
1545
1546 fVoxels.GetVoxel(curVoxel, currentPoint);
1547 do // Loop checking, 13.08.2015, G.Cosmo
1548 {
1549 const vector<G4int>& candidates = fVoxels.GetCandidates(curVoxel);
1550 if (!candidates.empty())
1551 {
1552 G4double distance=DistanceToInCandidates(candidates, aPoint, direction);
1553 if (minDistance > distance) minDistance = distance;
1554 if (distance < totalShift) break;
1555 }
1556
1557 shift = fVoxels.DistanceToNext(currentPoint, direction, curVoxel);
1558 if (shift == kInfinity /*|| shift == 0*/) break;
1559
1560 totalShift += shift;
1561 if (minDistance < totalShift) break;
1562
1563 currentPoint += direction * (shift + shiftBonus);
1564 }
1565 while (fVoxels.UpdateCurrentVoxel(currentPoint, direction, curVoxel));
1566 }
1567 else
1568 {
1569 minDistance = DistanceToInNoVoxels(aPoint, aDirection, aPstep);
1570 }
1571
1572 return minDistance;
1573 }
1574
1575 ///////////////////////////////////////////////////////////////////////////////
1576 //
1577 G4bool
1578 G4TessellatedSolid::CompareSortedVoxel(const std::pair<G4int, G4double>& l,
1579 const std::pair<G4int, G4double>& r)
1580 {
1581 return l.second < r.second;
1582 }
1583
1584 ///////////////////////////////////////////////////////////////////////////////
1585 //
1586 G4double
1587 G4TessellatedSolid::MinDistanceFacet(const G4ThreeVector& p,
1588 G4bool simple,
1589 G4VFacet* &minFacet) const
1590 {
1591 G4double minDist = kInfinity;
1592
1593 G4int size = fVoxels.GetVoxelBoxesSize();
1594 vector<pair<G4int, G4double> > voxelsSorted(size);
1595
1596 pair<G4int, G4double> info;
1597
1598 for (G4int i = 0; i < size; ++i)
1599 {
1600 const G4VoxelBox& voxelBox = fVoxels.GetVoxelBox(i);
1601
1602 G4ThreeVector pointShifted = p - voxelBox.pos;
1603 G4double safety = fVoxels.MinDistanceToBox(pointShifted, voxelBox.hlen);
1604 info.first = i;
1605 info.second = safety;
1606 voxelsSorted[i] = info;
1607 }
1608
1609 std::sort(voxelsSorted.begin(), voxelsSorted.end(),
1610 &G4TessellatedSolid::CompareSortedVoxel);
1611
1612 for (G4int i = 0; i < size; ++i)
1613 {
1614 const pair<G4int,G4double>& inf = voxelsSorted[i];
1615 G4double dist = inf.second;
1616 if (dist > minDist) break;
1617
1618 const vector<G4int>& candidates = fVoxels.GetVoxelBoxCandidates(inf.first);
1619 auto csize = (G4int)candidates.size();
1620 for (G4int j = 0; j < csize; ++j)
1621 {
1622 G4int candidate = candidates[j];
1623 G4VFacet& facet = *fFacets[candidate];
1624 dist = simple ? facet.Distance(p,minDist)
1625 : facet.Distance(p,minDist,false);
1626 if (dist < minDist)
1627 {
1628 minDist = dist;
1629 minFacet = &facet;
1630 }
1631 }
1632 }
1633 return minDist;
1634 }
1635
1636 ///////////////////////////////////////////////////////////////////////////////
1637 //
1638 G4double G4TessellatedSolid::SafetyFromOutside (const G4ThreeVector& p,
1639 G4bool aAccurate) const
1640 {
1641 #if G4SPECSDEBUG
1642 if ( Inside(p) == kInside )
1643 {
1644 std::ostringstream message;
1645 G4int oldprc = message.precision(16) ;
1646 message << "Point p is already inside!?" << G4endl
1647 << "Position:" << G4endl << G4endl
1648 << "p.x() = " << p.x()/mm << " mm" << G4endl
1649 << "p.y() = " << p.y()/mm << " mm" << G4endl
1650 << "p.z() = " << p.z()/mm << " mm" << G4endl
1651 << "DistanceToOut(p) == " << DistanceToOut(p);
1652 message.precision(oldprc) ;
1653 G4Exception("G4TriangularFacet::DistanceToIn(p)",
1654 "GeomSolids1002", JustWarning, message);
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 = fVoxels.GetCandidates(startingVoxel);
1670 if (candidates.empty() && (fInsides.GetNbits() != 0u))
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 {
1682 minDist = kInfinity;
1683 std::size_t size = fFacets.size();
1684 for (std::size_t i = 0; i < size; ++i)
1685 {
1686 G4VFacet& facet = *fFacets[i];
1687 G4double dist = facet.Distance(p,minDist);
1688 if (dist < minDist) minDist = dist;
1689 }
1690 }
1691 return minDist;
1692 }
1693
1694 ///////////////////////////////////////////////////////////////////////////////
1695 //
1696 G4double
1697 G4TessellatedSolid::SafetyFromInside (const G4ThreeVector& p, G4bool) const
1698 {
1699 #if G4SPECSDEBUG
1700 if ( Inside(p) == kOutside )
1701 {
1702 std::ostringstream message;
1703 G4int oldprc = message.precision(16) ;
1704 message << "Point p is already outside!?" << G4endl
1705 << "Position:" << G4endl << G4endl
1706 << "p.x() = " << p.x()/mm << " mm" << G4endl
1707 << "p.y() = " << p.y()/mm << " mm" << G4endl
1708 << "p.z() = " << p.z()/mm << " mm" << G4endl
1709 << "DistanceToIn(p) == " << DistanceToIn(p);
1710 message.precision(oldprc) ;
1711 G4Exception("G4TriangularFacet::DistanceToOut(p)",
1712 "GeomSolids1002", JustWarning, message);
1713 }
1714 #endif
1715
1716 G4double minDist;
1717
1718 if (OutsideOfExtent(p, kCarTolerance)) return 0.0;
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 }
1739
1740 ///////////////////////////////////////////////////////////////////////////////
1741 //
1742 // G4GeometryType GetEntityType() const;
1743 //
1744 // Provide identification of the class of an object
1745 //
1746 G4GeometryType G4TessellatedSolid::GetEntityType () const
1747 {
1748 return fGeometryType;
1749 }
1750
1751 ///////////////////////////////////////////////////////////////////////////////
1752 //
1753 std::ostream &G4TessellatedSolid::StreamInfo(std::ostream &os) const
1754 {
1755 os << G4endl;
1756 os << "Solid name = " << GetName() << G4endl;
1757 os << "Geometry Type = " << fGeometryType << G4endl;
1758 os << "Number of facets = " << fFacets.size() << G4endl;
1759
1760 std::size_t size = fFacets.size();
1761 for (std::size_t i = 0; i < size; ++i)
1762 {
1763 os << "FACET # = " << i + 1 << G4endl;
1764 G4VFacet &facet = *fFacets[i];
1765 facet.StreamInfo(os);
1766 }
1767 os << G4endl;
1768
1769 return os;
1770 }
1771
1772 ///////////////////////////////////////////////////////////////////////////////
1773 //
1774 // Make a clone of the object
1775 //
1776 G4VSolid* G4TessellatedSolid::Clone() const
1777 {
1778 return new G4TessellatedSolid(*this);
1779 }
1780
1781 ///////////////////////////////////////////////////////////////////////////////
1782 //
1783 // EInside G4TessellatedSolid::Inside (const G4ThreeVector &p) const
1784 //
1785 // This method must return:
1786 // * kOutside if the point at offset p is outside the shape
1787 // boundaries plus kCarTolerance/2,
1788 // * kSurface if the point is <= kCarTolerance/2 from a surface, or
1789 // * kInside otherwise.
1790 //
1791 EInside G4TessellatedSolid::Inside (const G4ThreeVector& aPoint) const
1792 {
1793 EInside location;
1794
1795 if (fVoxels.GetCountOfVoxels() > 1)
1796 {
1797 location = InsideVoxels(aPoint);
1798 }
1799 else
1800 {
1801 location = InsideNoVoxels(aPoint);
1802 }
1803 return location;
1804 }
1805
1806 ///////////////////////////////////////////////////////////////////////////////
1807 //
1808 G4ThreeVector G4TessellatedSolid::SurfaceNormal(const G4ThreeVector& p) const
1809 {
1810 G4ThreeVector n;
1811 Normal(p, n);
1812 return n;
1813 }
1814
1815 ///////////////////////////////////////////////////////////////////////////////
1816 //
1817 // G4double DistanceToIn(const G4ThreeVector& p)
1818 //
1819 // Calculate distance to nearest surface of shape from an outside point p. The
1820 // distance can be an underestimate.
1821 //
1822 G4double G4TessellatedSolid::DistanceToIn(const G4ThreeVector& p) const
1823 {
1824 return SafetyFromOutside(p, false);
1825 }
1826
1827 ///////////////////////////////////////////////////////////////////////////////
1828 //
1829 G4double G4TessellatedSolid::DistanceToIn(const G4ThreeVector& p,
1830 const G4ThreeVector& v)const
1831 {
1832 G4double dist = DistanceToInCore(p,v,kInfinity);
1833 #ifdef G4SPECSDEBUG
1834 if (dist < kInfinity)
1835 {
1836 if (Inside(p + dist*v) != kSurface)
1837 {
1838 std::ostringstream message;
1839 message << "Invalid response from facet in solid '" << GetName() << "',"
1840 << G4endl
1841 << "at point: " << p << "and direction: " << v;
1842 G4Exception("G4TessellatedSolid::DistanceToIn(p,v)",
1843 "GeomSolids1002", JustWarning, message);
1844 }
1845 }
1846 #endif
1847 return dist;
1848 }
1849
1850 ///////////////////////////////////////////////////////////////////////////////
1851 //
1852 // G4double DistanceToOut(const G4ThreeVector& p)
1853 //
1854 // Calculate distance to nearest surface of shape from an inside
1855 // point. The distance can be an underestimate.
1856 //
1857 G4double G4TessellatedSolid::DistanceToOut(const G4ThreeVector& p) const
1858 {
1859 return SafetyFromInside(p, false);
1860 }
1861
1862 ///////////////////////////////////////////////////////////////////////////////
1863 //
1864 // G4double DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
1865 // const G4bool calcNorm=false,
1866 // G4bool *validNorm=0, G4ThreeVector *n=0);
1867 //
1868 // Return distance along the normalised vector v to the shape, from a
1869 // point at an offset p inside or on the surface of the
1870 // shape. Intersections with surfaces, when the point is not greater
1871 // than kCarTolerance/2 from a surface, must be ignored.
1872 // If calcNorm is true, then it must also set validNorm to either
1873 // * true, if the solid lies entirely behind or on the exiting
1874 // surface. Then it must set n to the outwards normal vector
1875 // (the Magnitude of the vector is not defined).
1876 // * false, if the solid does not lie entirely behind or on the
1877 // exiting surface.
1878 // If calcNorm is false, then validNorm and n are unused.
1879 //
1880 G4double G4TessellatedSolid::DistanceToOut(const G4ThreeVector& p,
1881 const G4ThreeVector& v,
1882 const G4bool calcNorm,
1883 G4bool* validNorm,
1884 G4ThreeVector* norm) const
1885 {
1886 G4ThreeVector n;
1887 G4bool valid;
1888
1889 G4double dist = DistanceToOutCore(p, v, n, valid);
1890 if (calcNorm)
1891 {
1892 *norm = n;
1893 *validNorm = valid;
1894 }
1895 #ifdef G4SPECSDEBUG
1896 if (dist < kInfinity)
1897 {
1898 if (Inside(p + dist*v) != kSurface)
1899 {
1900 std::ostringstream message;
1901 message << "Invalid response from facet in solid '" << GetName() << "',"
1902 << G4endl
1903 << "at point: " << p << "and direction: " << v;
1904 G4Exception("G4TessellatedSolid::DistanceToOut(p,v,..)",
1905 "GeomSolids1002", JustWarning, message);
1906 }
1907 }
1908 #endif
1909 return dist;
1910 }
1911
1912 ///////////////////////////////////////////////////////////////////////////////
1913 //
1914 void G4TessellatedSolid::DescribeYourselfTo (G4VGraphicsScene& scene) const
1915 {
1916 scene.AddSolid (*this);
1917 }
1918
1919 ///////////////////////////////////////////////////////////////////////////////
1920 //
1921 G4Polyhedron* G4TessellatedSolid::CreatePolyhedron () const
1922 {
1923 auto nVertices = (G4int)fVertexList.size();
1924 auto nFacets = (G4int)fFacets.size();
1925 auto polyhedron = new G4Polyhedron(nVertices, nFacets);
1926 for (auto i = 0; i < nVertices; ++i)
1927 {
1928 polyhedron->SetVertex(i+1, fVertexList[i]);
1929 }
1930
1931 for (auto i = 0; i < nFacets; ++i)
1932 {
1933 G4VFacet* facet = fFacets[i];
1934 G4int v[4] = {0};
1935 G4int n = facet->GetNumberOfVertices();
1936 if (n > 4) n = 4;
1937 for (auto j = 0; j < n; ++j)
1938 {
1939 v[j] = facet->GetVertexIndex(j) + 1;
1940 }
1941 polyhedron->SetFacet(i+1, v[0], v[1], v[2], v[3]);
1942 }
1943 polyhedron->SetReferences();
1944
1945 return polyhedron;
1946 }
1947
1948 ///////////////////////////////////////////////////////////////////////////////
1949 //
1950 // GetPolyhedron
1951 //
1952 G4Polyhedron* G4TessellatedSolid::GetPolyhedron() const
1953 {
1954 if (fpPolyhedron == nullptr ||
1955 fRebuildPolyhedron ||
1956 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
1957 fpPolyhedron->GetNumberOfRotationSteps())
1958 {
1959 G4AutoLock l(&polyhedronMutex);
1960 delete fpPolyhedron;
1961 fpPolyhedron = CreatePolyhedron();
1962 fRebuildPolyhedron = false;
1963 l.unlock();
1964 }
1965 return fpPolyhedron;
1966 }
1967
1968 ///////////////////////////////////////////////////////////////////////////////
1969 //
1970 // Get bounding box
1971 //
1972 void G4TessellatedSolid::BoundingLimits(G4ThreeVector& pMin,
1973 G4ThreeVector& pMax) const
1974 {
1975 pMin = fMinExtent;
1976 pMax = fMaxExtent;
1977
1978 // Check correctness of the bounding box
1979 //
1980 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
1981 {
1982 std::ostringstream message;
1983 message << "Bad bounding box (min >= max) for solid: "
1984 << GetName() << " !"
1985 << "\npMin = " << pMin
1986 << "\npMax = " << pMax;
1987 G4Exception("G4TessellatedSolid::BoundingLimits()",
1988 "GeomMgt0001", JustWarning, message);
1989 DumpInfo();
1990 }
1991 }
1992
1993 ///////////////////////////////////////////////////////////////////////////////
1994 //
1995 // Calculate extent under transform and specified limit
1996 //
1997 G4bool
1998 G4TessellatedSolid::CalculateExtent(const EAxis pAxis,
1999 const G4VoxelLimits& pVoxelLimit,
2000 const G4AffineTransform& pTransform,
2001 G4double& pMin, G4double& pMax) const
2002 {
2003 G4ThreeVector bmin, bmax;
2004
2005 // Check bounding box (bbox)
2006 //
2007 BoundingLimits(bmin,bmax);
2008 G4BoundingEnvelope bbox(bmin,bmax);
2009
2010 // Use simple bounding-box to help in the case of complex meshes
2011 //
2012 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
2013
2014 #if 0
2015 // Precise extent computation (disabled by default for this shape)
2016 //
2017 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
2018 {
2019 return (pMin < pMax) ? true : false;
2020 }
2021
2022 // The extent is calculated as cumulative extent of the pyramids
2023 // formed by facets and the center of the bounding box.
2024 //
2025 G4double eminlim = pVoxelLimit.GetMinExtent(pAxis);
2026 G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis);
2027
2028 G4ThreeVectorList base;
2029 G4ThreeVectorList apex(1);
2030 std::vector<const G4ThreeVectorList *> pyramid(2);
2031 pyramid[0] = &base;
2032 pyramid[1] = &apex;
2033 apex[0] = (bmin+bmax)*0.5;
2034
2035 // main loop along facets
2036 pMin = kInfinity;
2037 pMax = -kInfinity;
2038 for (G4int i=0; i<GetNumberOfFacets(); ++i)
2039 {
2040 G4VFacet* facet = GetFacet(i);
2041 if (std::abs((facet->GetSurfaceNormal()).dot(facet->GetVertex(0)-apex[0]))
2042 < kCarToleranceHalf) continue;
2043
2044 G4int nv = facet->GetNumberOfVertices();
2045 base.resize(nv);
2046 for (G4int k=0; k<nv; ++k) { base[k] = facet->GetVertex(k); }
2047
2048 G4double emin,emax;
2049 G4BoundingEnvelope benv(pyramid);
2050 if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue;
2051 if (emin < pMin) pMin = emin;
2052 if (emax > pMax) pMax = emax;
2053 if (eminlim > pMin && emaxlim < pMax) break; // max possible extent
2054 }
2055 return (pMin < pMax);
2056 #endif
2057 }
2058
2059 ///////////////////////////////////////////////////////////////////////////////
2060 //
2061 G4double G4TessellatedSolid::GetMinXExtent () const
2062 {
2063 return fMinExtent.x();
2064 }
2065
2066 ///////////////////////////////////////////////////////////////////////////////
2067 //
2068 G4double G4TessellatedSolid::GetMaxXExtent () const
2069 {
2070 return fMaxExtent.x();
2071 }
2072
2073 ///////////////////////////////////////////////////////////////////////////////
2074 //
2075 G4double G4TessellatedSolid::GetMinYExtent () const
2076 {
2077 return fMinExtent.y();
2078 }
2079
2080 ///////////////////////////////////////////////////////////////////////////////
2081 //
2082 G4double G4TessellatedSolid::GetMaxYExtent () const
2083 {
2084 return fMaxExtent.y();
2085 }
2086
2087 ///////////////////////////////////////////////////////////////////////////////
2088 //
2089 G4double G4TessellatedSolid::GetMinZExtent () const
2090 {
2091 return fMinExtent.z();
2092 }
2093
2094 ///////////////////////////////////////////////////////////////////////////////
2095 //
2096 G4double G4TessellatedSolid::GetMaxZExtent () const
2097 {
2098 return fMaxExtent.z();
2099 }
2100
2101 ///////////////////////////////////////////////////////////////////////////////
2102 //
2103 G4VisExtent G4TessellatedSolid::GetExtent () const
2104 {
2105 return { fMinExtent.x(), fMaxExtent.x(),
2106 fMinExtent.y(), fMaxExtent.y(),
2107 fMinExtent.z(), fMaxExtent.z() };
2108 }
2109
2110 ///////////////////////////////////////////////////////////////////////////////
2111 //
2112 G4double G4TessellatedSolid::GetCubicVolume ()
2113 {
2114 if (fCubicVolume != 0.) return fCubicVolume;
2115
2116 // For explanation of the following algorithm see:
2117 // https://en.wikipedia.org/wiki/Polyhedron#Volume
2118 // http://wwwf.imperial.ac.uk/~rn/centroid.pdf
2119
2120 std::size_t size = fFacets.size();
2121 for (std::size_t i = 0; i < size; ++i)
2122 {
2123 G4VFacet &facet = *fFacets[i];
2124 G4double area = facet.GetArea();
2125 G4ThreeVector unit_normal = facet.GetSurfaceNormal();
2126 fCubicVolume += area * (facet.GetVertex(0).dot(unit_normal));
2127 }
2128 fCubicVolume /= 3.;
2129 return fCubicVolume;
2130 }
2131
2132 ///////////////////////////////////////////////////////////////////////////////
2133 //
2134 G4double G4TessellatedSolid::GetSurfaceArea ()
2135 {
2136 if (fSurfaceArea != 0.) return fSurfaceArea;
2137
2138 std::size_t size = fFacets.size();
2139 for (std::size_t i = 0; i < size; ++i)
2140 {
2141 G4VFacet &facet = *fFacets[i];
2142 fSurfaceArea += facet.GetArea();
2143 }
2144 return fSurfaceArea;
2145 }
2146
2147 ///////////////////////////////////////////////////////////////////////////////
2148 //
2149 G4ThreeVector G4TessellatedSolid::GetPointOnSurface() const
2150 {
2151 // Select randomly a facet and return a random point on it
2152
2153 auto i = (G4int) G4RandFlat::shoot(0., fFacets.size());
2154 return fFacets[i]->GetPointOnFace();
2155 }
2156
2157 ///////////////////////////////////////////////////////////////////////////////
2158 //
2159 // SetRandomVectorSet
2160 //
2161 // This is a set of predefined random vectors (if that isn't a contradition
2162 // in terms!) used to generate rays from a user-defined point. The member
2163 // function Inside uses these to determine whether the point is inside or
2164 // outside of the tessellated solid. All vectors should be unit vectors.
2165 //
2166 void G4TessellatedSolid::SetRandomVectors ()
2167 {
2168 fRandir.resize(20);
2169 fRandir[0] =
2170 G4ThreeVector(-0.9577428892113370, 0.2732676269591740, 0.0897405271949221);
2171 fRandir[1] =
2172 G4ThreeVector(-0.8331264504940770,-0.5162067214954600,-0.1985722492445700);
2173 fRandir[2] =
2174 G4ThreeVector(-0.1516671651108820, 0.9666292616127460, 0.2064580868390110);
2175 fRandir[3] =
2176 G4ThreeVector( 0.6570250350323190,-0.6944539025883300, 0.2933460081893360);
2177 fRandir[4] =
2178 G4ThreeVector(-0.4820456281280320,-0.6331060000098690,-0.6056474264406270);
2179 fRandir[5] =
2180 G4ThreeVector( 0.7629032554236800 , 0.1016854697539910,-0.6384658864065180);
2181 fRandir[6] =
2182 G4ThreeVector( 0.7689540409061150, 0.5034929891988220, 0.3939600142169160);
2183 fRandir[7] =
2184 G4ThreeVector( 0.5765188359255740, 0.5997271636278330,-0.5549354566343150);
2185 fRandir[8] =
2186 G4ThreeVector( 0.6660632777862070,-0.6362809868288380, 0.3892379937580790);
2187 fRandir[9] =
2188 G4ThreeVector( 0.3824415020414780, 0.6541792713761380,-0.6525243125110690);
2189 fRandir[10] =
2190 G4ThreeVector(-0.5107726564526760, 0.6020905056811610, 0.6136760679616570);
2191 fRandir[11] =
2192 G4ThreeVector( 0.7459135439578050, 0.6618796061649330, 0.0743530220183488);
2193 fRandir[12] =
2194 G4ThreeVector( 0.1536405855311580, 0.8117477913978260,-0.5634359711967240);
2195 fRandir[13] =
2196 G4ThreeVector( 0.0744395301705579,-0.8707110101772920,-0.4861286795736560);
2197 fRandir[14] =
2198 G4ThreeVector(-0.1665874645185400, 0.6018553940549240,-0.7810369397872780);
2199 fRandir[15] =
2200 G4ThreeVector( 0.7766902003633100, 0.6014617505959970,-0.1870724331097450);
2201 fRandir[16] =
2202 G4ThreeVector(-0.8710128685847430,-0.1434320216603030,-0.4698551243971010);
2203 fRandir[17] =
2204 G4ThreeVector( 0.8901082092766820,-0.4388411398893870, 0.1229871120030100);
2205 fRandir[18] =
2206 G4ThreeVector(-0.6430417431544370,-0.3295938228697690, 0.6912779675984150);
2207 fRandir[19] =
2208 G4ThreeVector( 0.6331124368380410, 0.6306211461665000, 0.4488714875425340);
2209
2210 fMaxTries = 20;
2211 }
2212
2213 ///////////////////////////////////////////////////////////////////////////////
2214 //
2215 G4int G4TessellatedSolid::AllocatedMemoryWithoutVoxels()
2216 {
2217 G4int base = sizeof(*this);
2218 base += fVertexList.capacity() * sizeof(G4ThreeVector);
2219 base += fRandir.capacity() * sizeof(G4ThreeVector);
2220
2221 std::size_t limit = fFacets.size();
2222 for (std::size_t i = 0; i < limit; ++i)
2223 {
2224 G4VFacet& facet = *fFacets[i];
2225 base += facet.AllocatedMemory();
2226 }
2227
2228 for (const auto & fExtremeFacet : fExtremeFacets)
2229 {
2230 G4VFacet &facet = *fExtremeFacet;
2231 base += facet.AllocatedMemory();
2232 }
2233 return base;
2234 }
2235
2236 ///////////////////////////////////////////////////////////////////////////////
2237 //
2238 G4int G4TessellatedSolid::AllocatedMemory()
2239 {
2240 G4int size = AllocatedMemoryWithoutVoxels();
2241 G4int sizeInsides = fInsides.GetNbytes();
2242 G4int sizeVoxels = fVoxels.AllocatedMemory();
2243 size += sizeInsides + sizeVoxels;
2244 return size;
2245 }
2246
2247 #endif
2248