Geant4  10.02.p01
G4TessellatedSolid.cc
Go to the documentation of this file.
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 // $Id: G4TessellatedSolid.cc 95424 2016-02-10 14:38:11Z gcosmo $
28 //
29 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
30 //
31 // CHANGE HISTORY
32 // --------------
33 // 12 October 2012, M Gayer, CERN, complete rewrite reducing memory
34 // requirements more than 50% and speedup by a factor of
35 // tens or more depending on the number of facets, thanks
36 // to voxelization of surface and improvements.
37 // Speedup factor of thousands for solids with number of
38 // facets in hundreds of thousands.
39 //
40 // 22 August 2011, I Hrivnacova, Orsay, fix in DistanceToOut(p) and
41 // DistanceToIn(p) to exactly compute distance from facet
42 // avoiding use of 'outgoing' flag shortcut variant.
43 //
44 // 04 August 2011, T Nikitina, CERN, added SetReferences() to
45 // CreatePolyhedron() for Visualization of Boolean Operations
46 //
47 // 12 April 2010, P R Truscott, QinetiQ, bug fixes to treat optical
48 // photon transport, in particular internal reflection
49 // at surface.
50 //
51 // 14 November 2007, P R Truscott, QinetiQ & Stan Seibert, U Texas
52 // Bug fixes to CalculateExtent
53 //
54 // 17 September 2007, P R Truscott, QinetiQ Ltd & Richard Holmberg
55 // Updated extensively prior to this date to deal with
56 // concaved tessellated surfaces, based on the algorithm
57 // of Richard Holmberg. This had been slightly modified
58 // to determine with inside the geometry by projecting
59 // random rays from the point provided. Now random rays
60 // are predefined rather than making use of random
61 // number generator at run-time.
62 //
63 // 22 November 2005, F Lei
64 // - Changed ::DescribeYourselfTo(), line 464
65 // - added GetPolyHedron()
66 //
67 // 31 October 2004, P R Truscott, QinetiQ Ltd, UK
68 // - Created.
69 //
70 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
71 
72 #include <iostream>
73 #include <stack>
74 #include <iostream>
75 #include <iomanip>
76 #include <fstream>
77 #include <algorithm>
78 #include <list>
79 
80 #include "G4TessellatedSolid.hh"
81 
82 #include "geomdefs.hh"
83 #include "Randomize.hh"
84 #include "G4SystemOfUnits.hh"
85 #include "G4PhysicalConstants.hh"
86 #include "G4GeometryTolerance.hh"
87 #include "G4VFacet.hh"
88 #include "G4VoxelLimits.hh"
89 #include "G4AffineTransform.hh"
90 
91 #include "G4PolyhedronArbitrary.hh"
92 #include "G4VGraphicsScene.hh"
93 #include "G4VisExtent.hh"
94 
95 #include "G4AutoLock.hh"
96 
97 namespace
98 {
99  G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
100 }
101 
102 using namespace std;
103 
105 //
106 // Standard contructor has blank name and defines no fFacets.
107 //
109 {
110  Initialize();
111 }
112 
114 //
115 // Alternative constructor. Simple define name and geometry type - no fFacets
116 // to detine.
117 //
119  : G4VSolid(name)
120 {
121  Initialize();
122 }
123 
125 //
126 // Fake default constructor - sets only member data and allocates memory
127 // for usage restricted to object persistency.
128 //
130 {
131  Initialize();
132  fMinExtent.set(0,0,0);
133  fMaxExtent.set(0,0,0);
134 }
135 
136 
139 {
140  DeleteObjects ();
141 }
142 
144 //
145 // Copy constructor.
146 //
148  : G4VSolid(ts), fpPolyhedron(0)
149 {
150  Initialize();
151 
152  CopyObjects(ts);
153 }
154 
156 //
157 // Assignment operator.
158 //
161 {
162  if (&ts == this) return *this;
163 
164  // Copy base class data
166 
167  DeleteObjects ();
168 
169  Initialize();
170 
171  CopyObjects (ts);
172 
173  return *this;
174 }
175 
177 //
179 {
181 
182  fRebuildPolyhedron = false; fpPolyhedron = 0;
183  fCubicVolume = 0.; fSurfaceArea = 0.;
184 
185  fGeometryType = "G4TessellatedSolid";
186  fSolidClosed = false;
187 
190 
192 }
193 
195 //
197 {
198  G4int size = fFacets.size();
199  for (G4int i = 0; i < size; ++i) { delete fFacets[i]; }
200  fFacets.clear();
201  delete fpPolyhedron; fpPolyhedron = 0;
202 }
203 
205 //
207 {
208  G4ThreeVector reductionRatio;
209  G4int fmaxVoxels = fVoxels.GetMaxVoxels(reductionRatio);
210  if (fmaxVoxels < 0)
211  fVoxels.SetMaxVoxels(reductionRatio);
212  else
213  fVoxels.SetMaxVoxels(fmaxVoxels);
214 
215  G4int n = ts.GetNumberOfFacets();
216  for (G4int i = 0; i < n; ++i)
217  {
218  G4VFacet *facetClone = (ts.GetFacet(i))->GetClone();
219  AddFacet(facetClone);
220  }
221  if (ts.GetSolidClosed()) SetSolidClosed(true);
222 }
223 
225 //
226 // Add a facet to the facet list.
227 // Note that you can add, but you cannot delete.
228 //
230 {
231  // Add the facet to the vector.
232  //
233  if (fSolidClosed)
234  {
235  G4Exception("G4TessellatedSolid::AddFacet()", "GeomSolids1002",
236  JustWarning, "Attempt to add facets when solid is closed.");
237  return false;
238  }
239  else if (aFacet->IsDefined())
240  {
241  set<G4VertexInfo,G4VertexComparator>::iterator begin
242  = fFacetList.begin(), end = fFacetList.end(), pos, it;
243  G4ThreeVector p = aFacet->GetCircumcentre();
244  G4VertexInfo value;
245  value.id = fFacetList.size();
246  value.mag2 = p.x() + p.y() + p.z();
247 
248  G4bool found = false;
250  {
251  G4double kCarTolerance3 = 3 * kCarTolerance;
252  pos = fFacetList.lower_bound(value);
253 
254  it = pos;
255  while (!found && it != end) // Loop checking, 13.08.2015, G.Cosmo
256  {
257  G4int id = (*it).id;
258  G4VFacet *facet = fFacets[id];
259  G4ThreeVector q = facet->GetCircumcentre();
260  if ((found = (facet == aFacet))) break;
261  G4double dif = q.x() + q.y() + q.z() - value.mag2;
262  if (dif > kCarTolerance3) break;
263  it++;
264  }
265 
266  if (fFacets.size() > 1)
267  {
268  it = pos;
269  while (!found && it != begin) // Loop checking, 13.08.2015, G.Cosmo
270  {
271  --it;
272  G4int id = (*it).id;
273  G4VFacet *facet = fFacets[id];
274  G4ThreeVector q = facet->GetCircumcentre();
275  found = (facet == aFacet);
276  if (found) break;
277  G4double dif = value.mag2 - (q.x() + q.y() + q.z());
278  if (dif > kCarTolerance3) break;
279  }
280  }
281  }
282 
283  if (!found)
284  {
285  fFacets.push_back(aFacet);
286  fFacetList.insert(value);
287  }
288 
289  return true;
290  }
291  else
292  {
293  G4Exception("G4TessellatedSolid::AddFacet()", "GeomSolids1002",
294  JustWarning, "Attempt to add facet not properly defined.");
295  aFacet->StreamInfo(G4cout);
296  return false;
297  }
298 }
299 
301 //
302 G4int G4TessellatedSolid::SetAllUsingStack(const std::vector<G4int> &voxel,
303  const std::vector<G4int> &max,
304  G4bool status, G4SurfBits &checked)
305 {
306  vector<G4int> xyz = voxel;
307  stack<vector<G4int> > pos;
308  pos.push(xyz);
309  G4int filled = 0;
310  G4int cc = 0, nz = 0;
311 
312  vector<G4int> candidates;
313 
314  while (!pos.empty()) // Loop checking, 13.08.2015, G.Cosmo
315  {
316  xyz = pos.top();
317  pos.pop();
318  G4int index = fVoxels.GetVoxelsIndex(xyz);
319  if (!checked[index])
320  {
321  checked.SetBitNumber(index, true);
322  cc++;
323 
324  if (fVoxels.IsEmpty(index))
325  {
326  filled++;
327 
328  fInsides.SetBitNumber(index, status);
329 
330  for (G4int i = 0; i <= 2; ++i)
331  {
332  if (xyz[i] < max[i] - 1)
333  {
334  xyz[i]++;
335  pos.push(xyz);
336  xyz[i]--;
337  }
338 
339  if (xyz[i] > 0)
340  {
341  xyz[i]--;
342  pos.push(xyz);
343  xyz[i]++;
344  }
345  }
346  }
347  else
348  {
349  nz++;
350  }
351  }
352  }
353  return filled;
354 }
355 
357 //
359 {
360  vector<G4int> voxel(3), maxVoxels(3);
361  for (G4int i = 0; i <= 2; ++i) maxVoxels[i] = fVoxels.GetBoundary(i).size();
362  G4int size = maxVoxels[0] * maxVoxels[1] * maxVoxels[2];
363 
364  G4SurfBits checked(size-1);
365  fInsides.Clear();
366  fInsides.ResetBitNumber(size-1);
367 
368  G4ThreeVector point;
369  for (voxel[2] = 0; voxel[2] < maxVoxels[2] - 1; ++voxel[2])
370  {
371  for (voxel[1] = 0; voxel[1] < maxVoxels[1] - 1; ++voxel[1])
372  {
373  for (voxel[0] = 0; voxel[0] < maxVoxels[0] - 1; ++voxel[0])
374  {
375  G4int index = fVoxels.GetVoxelsIndex(voxel);
376  if (!checked[index] && fVoxels.IsEmpty(index))
377  {
378  for (G4int i = 0; i <= 2; ++i) point[i] = fVoxels.GetBoundary(i)[voxel[i]];
379  G4bool inside = (G4bool) (InsideNoVoxels(point) == kInside);
380  SetAllUsingStack(voxel, maxVoxels, inside, checked);
381  }
382  else checked.SetBitNumber(index);
383  }
384  }
385  }
386 }
387 
389 //
391 {
392 #ifdef G4SPECSDEBUG
393  G4cout << "Voxelizing..." << G4endl;
394 #endif
396 
397  if (fVoxels.Empty().GetNbits())
398  {
399 #ifdef G4SPECSDEBUG
400  G4cout << "Precalculating Insides..." << G4endl;
401 #endif
403  }
404 }
405 
407 //
408 // Compute extremeFacets, i.e. find those facets that have surface
409 // planes that bound the volume.
410 // Note that this is going to reject concaved surfaces as being extreme. Also
411 // note that if the vertex is on the facet, displacement is zero, so IsInside
412 // returns true. So will this work?? Need non-equality
413 // "G4bool inside = displacement < 0.0;"
414 // or
415 // "G4bool inside = displacement <= -0.5*kCarTolerance"
416 // (Notes from PT 13/08/2007).
417 //
419 {
420  G4int size = fFacets.size();
421  for (G4int j = 0; j < size; ++j)
422  {
423  G4VFacet &facet = *fFacets[j];
424 
425  G4bool isExtreme = true;
426  G4int vsize = fVertexList.size();
427  for (G4int i=0; i < vsize; ++i)
428  {
429  if (!facet.IsInside(fVertexList[i]))
430  {
431  isExtreme = false;
432  break;
433  }
434  }
435  if (isExtreme) fExtremeFacets.insert(&facet);
436  }
437 }
438 
440 //
442 {
443  // The algorithm:
444  // we will have additional vertexListSorted, where all the items will be
445  // sorted by magnitude of vertice vector.
446  // New candidate for fVertexList - we will determine the position fo first
447  // item which would be within its magnitude - 0.5*kCarTolerance.
448  // We will go trough until we will reach > +0.5 kCarTolerance.
449  // Comparison (q-p).mag() < 0.5*kCarTolerance will be made.
450  // They can be just stored in std::vector, with custom insertion based
451  // on binary search.
452 
453  set<G4VertexInfo,G4VertexComparator> vertexListSorted;
454  set<G4VertexInfo,G4VertexComparator>::iterator begin
455  = vertexListSorted.begin(), end = vertexListSorted.end(), pos, it;
456  G4ThreeVector p;
457  G4VertexInfo value;
458 
459  fVertexList.clear();
460  G4int size = fFacets.size();
461 
462  G4double kCarTolerance24 = kCarTolerance * kCarTolerance / 4.0;
463  G4double kCarTolerance3 = 3 * kCarTolerance;
464  vector<G4int> newIndex(100);
465 
466  for (G4int k = 0; k < size; ++k)
467  {
468  G4VFacet &facet = *fFacets[k];
469  G4int max = facet.GetNumberOfVertices();
470 
471  for (G4int i = 0; i < max; ++i)
472  {
473  p = facet.GetVertex(i);
474  value.id = fVertexList.size();
475  value.mag2 = p.x() + p.y() + p.z();
476 
477  G4bool found = false;
478  G4int id = 0;
480  {
481  pos = vertexListSorted.lower_bound(value);
482  it = pos;
483  while (it != end) // Loop checking, 13.08.2015, G.Cosmo
484  {
485  id = (*it).id;
486  G4ThreeVector q = fVertexList[id];
487  G4double dif = (q-p).mag2();
488  found = (dif < kCarTolerance24);
489  if (found) break;
490  dif = q.x() + q.y() + q.z() - value.mag2;
491  if (dif > kCarTolerance3) break;
492  it++;
493  }
494 
495  if (!found && (fVertexList.size() > 1))
496  {
497  it = pos;
498  while (it != begin) // Loop checking, 13.08.2015, G.Cosmo
499  {
500  --it;
501  id = (*it).id;
502  G4ThreeVector q = fVertexList[id];
503  G4double dif = (q-p).mag2();
504  found = (dif < kCarTolerance24);
505  if (found) break;
506  dif = value.mag2 - (q.x() + q.y() + q.z());
507  if (dif > kCarTolerance3) break;
508  }
509  }
510  }
511 
512  if (!found)
513  {
514 #ifdef G4SPECSDEBUG
515  G4cout << p.x() << ":" << p.y() << ":" << p.z() << G4endl;
516  G4cout << "Adding new vertex #" << i << " of facet " << k
517  << " id " << value.id << G4endl;
518  G4cout << "===" << G4endl;
519 #endif
520  fVertexList.push_back(p);
521  vertexListSorted.insert(value);
522  begin = vertexListSorted.begin();
523  end = vertexListSorted.end();
524  newIndex[i] = value.id;
525  //
526  // Now update the maximum x, y and z limits of the volume.
527  //
528  if (value.id == 0) fMinExtent = fMaxExtent = p;
529  else
530  {
531  if (p.x() > fMaxExtent.x()) fMaxExtent.setX(p.x());
532  else if (p.x() < fMinExtent.x()) fMinExtent.setX(p.x());
533  if (p.y() > fMaxExtent.y()) fMaxExtent.setY(p.y());
534  else if (p.y() < fMinExtent.y()) fMinExtent.setY(p.y());
535  if (p.z() > fMaxExtent.z()) fMaxExtent.setZ(p.z());
536  else if (p.z() < fMinExtent.z()) fMinExtent.setZ(p.z());
537  }
538  }
539  else
540  {
541 #ifdef G4SPECSDEBUG
542  G4cout << p.x() << ":" << p.y() << ":" << p.z() << G4endl;
543  G4cout << "Vertex #" << i << " of facet " << k
544  << " found, redirecting to " << id << G4endl;
545  G4cout << "===" << G4endl;
546 #endif
547  newIndex[i] = id;
548  }
549  }
550  // only now it is possible to change vertices pointer
551  //
552  facet.SetVertices(&fVertexList);
553  for (G4int i = 0; i < max; i++)
554  facet.SetVertexIndex(i,newIndex[i]);
555  }
556  vector<G4ThreeVector>(fVertexList).swap(fVertexList);
557 
558 #ifdef G4SPECSDEBUG
559  G4double previousValue = 0;
560  for (set<G4VertexInfo,G4VertexComparator>::iterator res=
561  vertexListSorted.begin(); res!=vertexListSorted.end(); ++res)
562  {
563  G4int id = (*res).id;
564  G4ThreeVector vec = fVertexList[id];
565  G4double mvalue = vec.x() + vec.y() + vec.z();
566  if (previousValue && (previousValue - 1e-9 > mvalue))
567  G4cout << "Error in CreateVertexList: previousValue " << previousValue
568  << " is smaller than mvalue " << mvalue << G4endl;
569  previousValue = mvalue;
570  }
571 #endif
572 }
573 
575 //
577 {
579  G4int with = AllocatedMemory();
580  G4double ratio = (G4double) with / without;
581  G4cout << "G4TessellatedSolid - Allocated memory without voxel overhead "
582  << without << "; with " << with << "; ratio: " << ratio << G4endl;
583 }
584 
586 //
588 {
589  if (t)
590  {
591 #ifdef G4SPECSDEBUG
592  G4cout << "Creating vertex list..." << G4endl;
593 #endif
595 
596 #ifdef G4SPECSDEBUG
597  G4cout << "Setting extreme facets..." << G4endl;
598 #endif
600 
601 #ifdef G4SPECSDEBUG
602  G4cout << "Voxelizing..." << G4endl;
603 #endif
604  Voxelize();
605 
606 #ifdef G4SPECSDEBUG
608 #endif
609 
610  }
611  fSolidClosed = t;
612 }
613 
615 //
616 // GetSolidClosed
617 //
618 // Used to determine whether the solid is closed to adding further fFacets.
619 //
621 {
622  return fSolidClosed;
623 }
624 
626 //
627 // operator+=
628 //
629 // This operator allows the user to add two tessellated solids together, so
630 // that the solid on the left then includes all of the facets in the solid
631 // on the right. Note that copies of the facets are generated, rather than
632 // using the original facet set of the solid on the right.
633 //
636 {
637  G4int size = right.GetNumberOfFacets();
638  for (G4int i = 0; i < size; ++i)
639  AddFacet(right.GetFacet(i)->GetClone());
640 
641  return *this;
642 }
643 
645 //
646 // GetNumberOfFacets
647 //
649 {
650  return fFacets.size();
651 }
652 
654 //
656 {
657  //
658  // First the simple test - check if we're outside of the X-Y-Z extremes
659  // of the tessellated solid.
660  //
662  return kOutside;
663 
664  vector<G4int> startingVoxel(3);
665  fVoxels.GetVoxel(startingVoxel, p);
666 
667  const G4double dirTolerance = 1.0E-14;
668 
669  const vector<G4int> &startingCandidates =
670  fVoxels.GetCandidates(startingVoxel);
671  G4int limit = startingCandidates.size();
672  if (limit == 0 && fInsides.GetNbits())
673  {
674  G4int index = fVoxels.GetPointIndex(p);
675  EInside location = fInsides[index] ? kInside : kOutside;
676  return location;
677  }
678 
679  G4double minDist = kInfinity;
680 
681  for(G4int i = 0; i < limit; ++i)
682  {
683  G4int candidate = startingCandidates[i];
684  G4VFacet &facet = *fFacets[candidate];
685  G4double dist = facet.Distance(p,minDist);
686  if (dist < minDist) minDist = dist;
687  if (dist <= kCarToleranceHalf)
688  return kSurface;
689  }
690 
691  // The following is something of an adaptation of the method implemented by
692  // Rickard Holmberg augmented with information from Schneider & Eberly,
693  // "Geometric Tools for Computer Graphics," pp700-701, 2003. In essence,
694  // we're trying to determine whether we're inside the volume by projecting
695  // a few rays and determining if the first surface crossed is has a normal
696  // vector between 0 to pi/2 (out-going) or pi/2 to pi (in-going).
697  // We should also avoid rays which are nearly within the plane of the
698  // tessellated surface, and therefore produce rays randomly.
699  // For the moment, this is a bit over-engineered (belt-braces-and-ducttape).
700  //
701  G4double distOut = kInfinity;
702  G4double distIn = kInfinity;
703  G4double distO = 0.0;
704  G4double distI = 0.0;
705  G4double distFromSurfaceO = 0.0;
706  G4double distFromSurfaceI = 0.0;
707  G4ThreeVector normalO, normalI;
708  G4bool crossingO = false;
709  G4bool crossingI = false;
710  EInside location = kOutside;
711  G4int sm = 0;
712 
713  G4bool nearParallel = false;
714  do // Loop checking, 13.08.2015, G.Cosmo
715  {
716  // We loop until we find direction where the vector is not nearly parallel
717  // to the surface of any facet since this causes ambiguities. The usual
718  // case is that the angles should be sufficiently different, but there
719  // are 20 random directions to select from - hopefully sufficient.
720  //
721  distOut = distIn = kInfinity;
722  const G4ThreeVector &v = fRandir[sm];
723  sm++;
724  //
725  // This code could be voxelized by the same algorithm, which is used for
726  // DistanceToOut(). We will traverse through fVoxels. we will call
727  // intersect only for those, which would be candidates and was not
728  // checked before.
729  //
730  G4ThreeVector currentPoint = p;
731  G4ThreeVector direction = v.unit();
732  // G4SurfBits exclusion(fVoxels.GetBitsPerSlice());
733  vector<G4int> curVoxel(3);
734  curVoxel = startingVoxel;
735  G4double shiftBonus = kCarTolerance;
736 
737  G4bool crossed = false;
738  G4bool started = true;
739 
740  do // Loop checking, 13.08.2015, G.Cosmo
741  {
742  const vector<G4int> &candidates =
743  started ? startingCandidates : fVoxels.GetCandidates(curVoxel);
744  started = false;
745  if (G4int candidatesCount = candidates.size())
746  {
747  for (G4int i = 0 ; i < candidatesCount; ++i)
748  {
749  G4int candidate = candidates[i];
750  // bits.SetBitNumber(candidate);
751  G4VFacet &facet = *fFacets[candidate];
752 
753  crossingO = facet.Intersect(p,v,true,distO,distFromSurfaceO,normalO);
754  crossingI = facet.Intersect(p,v,false,distI,distFromSurfaceI,normalI);
755 
756  if (crossingO || crossingI)
757  {
758  crossed = true;
759 
760  nearParallel = (crossingO
761  && std::fabs(normalO.dot(v))<dirTolerance)
762  || (crossingI && std::fabs(normalI.dot(v))<dirTolerance);
763  if (!nearParallel)
764  {
765  if (crossingO && distO > 0.0 && distO < distOut)
766  distOut = distO;
767  if (crossingI && distI > 0.0 && distI < distIn)
768  distIn = distI;
769  }
770  else break;
771  }
772  }
773  if (nearParallel) break;
774  }
775  else
776  {
777  if (!crossed)
778  {
779  G4int index = fVoxels.GetVoxelsIndex(curVoxel);
780  G4bool inside = fInsides[index];
781  location = inside ? kInside : kOutside;
782  return location;
783  }
784  }
785 
786  G4double shift=fVoxels.DistanceToNext(currentPoint, direction, curVoxel);
787  if (shift == kInfinity) break;
788 
789  currentPoint += direction * (shift + shiftBonus);
790  }
791  while (fVoxels.UpdateCurrentVoxel(currentPoint, direction, curVoxel));
792 
793  }
794  while (nearParallel && sm!=fMaxTries);
795  //
796  // Here we loop through the facets to find out if there is an intersection
797  // between the ray and that facet. The test if performed separately whether
798  // the ray is entering the facet or exiting.
799  //
800 #ifdef G4VERBOSE
801  if (sm == fMaxTries)
802  {
803  //
804  // We've run out of random vector directions. If nTries is set sufficiently
805  // low (nTries <= 0.5*maxTries) then this would indicate that there is
806  // something wrong with geometry.
807  //
808  std::ostringstream message;
809  G4int oldprc = message.precision(16);
810  message << "Cannot determine whether point is inside or outside volume!"
811  << G4endl
812  << "Solid name = " << GetName() << G4endl
813  << "Geometry Type = " << fGeometryType << G4endl
814  << "Number of facets = " << fFacets.size() << G4endl
815  << "Position:" << G4endl << G4endl
816  << "p.x() = " << p.x()/mm << " mm" << G4endl
817  << "p.y() = " << p.y()/mm << " mm" << G4endl
818  << "p.z() = " << p.z()/mm << " mm";
819  message.precision(oldprc);
820  G4Exception("G4TessellatedSolid::Inside()",
821  "GeomSolids1002", JustWarning, message);
822  }
823 #endif
824 
825  // In the next if-then-elseif G4String the logic is as follows:
826  // (1) You don't hit anything so cannot be inside volume, provided volume
827  // constructed correctly!
828  // (2) Distance to inside (ie. nearest facet such that you enter facet) is
829  // shorter than distance to outside (nearest facet such that you exit
830  // facet) - on condition of safety distance - therefore we're outside.
831  // (3) Distance to outside is shorter than distance to inside therefore
832  // we're inside.
833  //
834  if (distIn == kInfinity && distOut == kInfinity)
835  location = kOutside;
836  else if (distIn <= distOut - kCarToleranceHalf)
837  location = kOutside;
838  else if (distOut <= distIn - kCarToleranceHalf)
839  location = kInside;
840 
841  return location;
842 }
843 
845 //
847 {
848  //
849  // First the simple test - check if we're outside of the X-Y-Z extremes
850  // of the tessellated solid.
851  //
853  return kOutside;
854 
855  const G4double dirTolerance = 1.0E-14;
856 
857  G4double minDist = kInfinity;
858  //
859  // Check if we are close to a surface
860  //
861  G4int size = fFacets.size();
862  for (G4int i = 0; i < size; ++i)
863  {
864  G4VFacet &facet = *fFacets[i];
865  G4double dist = facet.Distance(p,minDist);
866  if (dist < minDist) minDist = dist;
867  if (dist <= kCarToleranceHalf)
868  {
869  return kSurface;
870  }
871  }
872  //
873  // The following is something of an adaptation of the method implemented by
874  // Rickard Holmberg augmented with information from Schneider & Eberly,
875  // "Geometric Tools for Computer Graphics," pp700-701, 2003. In essence, we're
876  // trying to determine whether we're inside the volume by projecting a few
877  // rays and determining if the first surface crossed is has a normal vector
878  // between 0 to pi/2 (out-going) or pi/2 to pi (in-going). We should also
879  // avoid rays which are nearly within the plane of the tessellated surface,
880  // and therefore produce rays randomly. For the moment, this is a bit
881  // over-engineered (belt-braces-and-ducttape).
882  //
883 #if G4SPECSDEBUG
884  G4int nTry = 7;
885 #else
886  G4int nTry = 3;
887 #endif
888  G4double distOut = kInfinity;
889  G4double distIn = kInfinity;
890  G4double distO = 0.0;
891  G4double distI = 0.0;
892  G4double distFromSurfaceO = 0.0;
893  G4double distFromSurfaceI = 0.0;
894  G4ThreeVector normalO(0.0,0.0,0.0);
895  G4ThreeVector normalI(0.0,0.0,0.0);
896  G4bool crossingO = false;
897  G4bool crossingI = false;
898  EInside location = kOutside;
899  EInside locationprime = kOutside;
900  G4int sm = 0;
901 
902  for (G4int i=0; i<nTry; ++i)
903  {
904  G4bool nearParallel = false;
905  do // Loop checking, 13.08.2015, G.Cosmo
906  {
907  //
908  // We loop until we find direction where the vector is not nearly parallel
909  // to the surface of any facet since this causes ambiguities. The usual
910  // case is that the angles should be sufficiently different, but there
911  // are 20 random directions to select from - hopefully sufficient.
912  //
913  distOut = distIn = kInfinity;
914  G4ThreeVector v = fRandir[sm];
915  sm++;
916  vector<G4VFacet*>::const_iterator f = fFacets.begin();
917 
918  do // Loop checking, 13.08.2015, G.Cosmo
919  {
920  //
921  // Here we loop through the facets to find out if there is an
922  // intersection between the ray and that facet. The test if performed
923  // separately whether the ray is entering the facet or exiting.
924  //
925  crossingO = ((*f)->Intersect(p,v,true,distO,distFromSurfaceO,normalO));
926  crossingI = ((*f)->Intersect(p,v,false,distI,distFromSurfaceI,normalI));
927  if (crossingO || crossingI)
928  {
929  nearParallel = (crossingO && std::fabs(normalO.dot(v))<dirTolerance)
930  || (crossingI && std::fabs(normalI.dot(v))<dirTolerance);
931  if (!nearParallel)
932  {
933  if (crossingO && distO > 0.0 && distO < distOut) distOut = distO;
934  if (crossingI && distI > 0.0 && distI < distIn) distIn = distI;
935  }
936  }
937  } while (!nearParallel && ++f!=fFacets.end());
938  } while (nearParallel && sm!=fMaxTries);
939 
940 #ifdef G4VERBOSE
941  if (sm == fMaxTries)
942  {
943  //
944  // We've run out of random vector directions. If nTries is set
945  // sufficiently low (nTries <= 0.5*maxTries) then this would indicate
946  // that there is something wrong with geometry.
947  //
948  std::ostringstream message;
949  G4int oldprc = message.precision(16);
950  message << "Cannot determine whether point is inside or outside volume!"
951  << G4endl
952  << "Solid name = " << GetName() << G4endl
953  << "Geometry Type = " << fGeometryType << G4endl
954  << "Number of facets = " << fFacets.size() << G4endl
955  << "Position:" << G4endl << G4endl
956  << "p.x() = " << p.x()/mm << " mm" << G4endl
957  << "p.y() = " << p.y()/mm << " mm" << G4endl
958  << "p.z() = " << p.z()/mm << " mm";
959  message.precision(oldprc);
960  G4Exception("G4TessellatedSolid::Inside()",
961  "GeomSolids1002", JustWarning, message);
962  }
963 #endif
964  //
965  // In the next if-then-elseif G4String the logic is as follows:
966  // (1) You don't hit anything so cannot be inside volume, provided volume
967  // constructed correctly!
968  // (2) Distance to inside (ie. nearest facet such that you enter facet) is
969  // shorter than distance to outside (nearest facet such that you exit
970  // facet) - on condition of safety distance - therefore we're outside.
971  // (3) Distance to outside is shorter than distance to inside therefore
972  // we're inside.
973  //
974  if (distIn == kInfinity && distOut == kInfinity)
975  locationprime = kOutside;
976  else if (distIn <= distOut - kCarToleranceHalf)
977  locationprime = kOutside;
978  else if (distOut <= distIn - kCarToleranceHalf)
979  locationprime = kInside;
980 
981  if (i == 0) location = locationprime;
982  }
983 
984  return location;
985 }
986 
988 //
989 // Return the outwards pointing unit normal of the shape for the
990 // surface closest to the point at offset p.
991 //
993  G4ThreeVector &aNormal) const
994 {
995  G4double minDist;
996  G4VFacet *facet = 0;
997 
998  if (fVoxels.GetCountOfVoxels() > 1)
999  {
1000  vector<G4int> curVoxel(3);
1001  fVoxels.GetVoxel(curVoxel, p);
1002  const vector<G4int> &candidates = fVoxels.GetCandidates(curVoxel);
1003  // fVoxels.GetCandidatesVoxelArray(p, candidates, 0);
1004 
1005  if (G4int limit = candidates.size())
1006  {
1007  minDist = kInfinity;
1008  for(G4int i = 0 ; i < limit ; ++i)
1009  {
1010  G4int candidate = candidates[i];
1011  G4VFacet &fct = *fFacets[candidate];
1012  G4double dist = fct.Distance(p,minDist);
1013  if (dist < minDist) minDist = dist;
1014  if (dist <= kCarToleranceHalf)
1015  {
1016  aNormal = fct.GetSurfaceNormal();
1017  return true;
1018  }
1019  }
1020  }
1021  minDist = MinDistanceFacet(p, true, facet);
1022  }
1023  else
1024  {
1025  minDist = kInfinity;
1026  G4int size = fFacets.size();
1027  for (G4int i = 0; i < size; ++i)
1028  {
1029  G4VFacet &f = *fFacets[i];
1030  G4double dist = f.Distance(p, minDist);
1031  if (dist < minDist)
1032  {
1033  minDist = dist;
1034  facet = &f;
1035  }
1036  }
1037  }
1038 
1039  if (minDist != kInfinity)
1040  {
1041  if (facet) { aNormal = facet->GetSurfaceNormal(); }
1042  return minDist <= kCarToleranceHalf;
1043  }
1044  else
1045  {
1046 #ifdef G4VERBOSE
1047  std::ostringstream message;
1048  message << "Point p is not on surface !?" << G4endl
1049  << " No facets found for point: " << p << " !" << G4endl
1050  << " Returning approximated value for normal.";
1051 
1052  G4Exception("G4TessellatedSolid::SurfaceNormal(p)",
1053  "GeomSolids1002", JustWarning, message );
1054 #endif
1055  aNormal = (p.z() > 0 ? G4ThreeVector(0,0,1) : G4ThreeVector(0,0,-1));
1056  return false;
1057  }
1058 }
1059 
1061 //
1062 // G4double DistanceToIn(const G4ThreeVector& p, const G4ThreeVector& v)
1063 //
1064 // Return the distance along the normalised vector v to the shape,
1065 // from the point at offset p. If there is no intersection, return
1066 // kInfinity. The first intersection resulting from 'leaving' a
1067 // surface/volume is discarded. Hence, this is tolerant of points on
1068 // surface of shape.
1069 //
1070 G4double
1072  const G4ThreeVector &v,
1073  G4double /*aPstep*/) const
1074 {
1075  G4double minDist = kInfinity;
1076  G4double dist = 0.0;
1077  G4double distFromSurface = 0.0;
1079 
1080 #if G4SPECSDEBUG
1081  if (Inside(p) == kInside )
1082  {
1083  std::ostringstream message;
1084  G4int oldprc = message.precision(16) ;
1085  message << "Point p is already inside!?" << G4endl
1086  << "Position:" << G4endl << G4endl
1087  << " p.x() = " << p.x()/mm << " mm" << G4endl
1088  << " p.y() = " << p.y()/mm << " mm" << G4endl
1089  << " p.z() = " << p.z()/mm << " mm" << G4endl
1090  << "DistanceToOut(p) == " << DistanceToOut(p);
1091  message.precision(oldprc) ;
1092  G4Exception("G4TriangularFacet::DistanceToIn(p,v)",
1093  "GeomSolids1002", JustWarning, message);
1094  }
1095 #endif
1096 
1097  G4int size = fFacets.size();
1098  for (G4int i = 0; i < size; ++i)
1099  {
1100  G4VFacet &facet = *fFacets[i];
1101  if (facet.Intersect(p,v,false,dist,distFromSurface,normal))
1102  {
1103  //
1104  // set minDist to the new distance to current facet if distFromSurface
1105  // is in positive direction and point is not at surface. If the point is
1106  // within 0.5*kCarTolerance of the surface, then force distance to be
1107  // zero and leave member function immediately (for efficiency), as
1108  // proposed by & credit to Akira Okumura.
1109  //
1110  if (distFromSurface > kCarToleranceHalf && dist >= 0.0 && dist < minDist)
1111  {
1112  minDist = dist;
1113  }
1114  else
1115  {
1116  if (-kCarToleranceHalf <= dist && dist <= kCarToleranceHalf)
1117  {
1118  return 0.0;
1119  }
1120  else
1121  {
1122  if (distFromSurface > -kCarToleranceHalf
1123  && distFromSurface < kCarToleranceHalf)
1124  {
1125  minDist = dist;
1126  }
1127  }
1128  }
1129  }
1130  }
1131  return minDist;
1132 }
1133 
1135 //
1136 G4double
1138  const G4ThreeVector &v,
1139  G4ThreeVector &aNormalVector,
1140  G4bool &aConvex,
1141  G4double /*aPstep*/) const
1142 {
1143  G4double minDist = kInfinity;
1144  G4double dist = 0.0;
1145  G4double distFromSurface = 0.0;
1146  G4ThreeVector normal, minNormal;
1147 
1148 #if G4SPECSDEBUG
1149  if ( Inside(p) == kOutside )
1150  {
1151  std::ostringstream message;
1152  G4int oldprc = message.precision(16) ;
1153  message << "Point p is already outside!?" << G4endl
1154  << "Position:" << G4endl << G4endl
1155  << " p.x() = " << p.x()/mm << " mm" << G4endl
1156  << " p.y() = " << p.y()/mm << " mm" << G4endl
1157  << " p.z() = " << p.z()/mm << " mm" << G4endl
1158  << "DistanceToIn(p) == " << DistanceToIn(p);
1159  message.precision(oldprc) ;
1160  G4Exception("G4TriangularFacet::DistanceToOut(p)",
1161  "GeomSolids1002", JustWarning, message);
1162  }
1163 #endif
1164 
1165  G4bool isExtreme = false;
1166  G4int size = fFacets.size();
1167  for (G4int i = 0; i < size; ++i)
1168  {
1169  G4VFacet &facet = *fFacets[i];
1170  if (facet.Intersect(p,v,true,dist,distFromSurface,normal))
1171  {
1172  if (distFromSurface > 0.0 && distFromSurface <= kCarToleranceHalf &&
1174  {
1175  // We are on a surface. Return zero.
1176  aConvex = (fExtremeFacets.find(&facet) != fExtremeFacets.end());
1177  // Normal(p, aNormalVector);
1178  // aNormalVector = facet.GetSurfaceNormal();
1179  aNormalVector = normal;
1180  return 0.0;
1181  }
1182  if (dist >= 0.0 && dist < minDist)
1183  {
1184  minDist = dist;
1185  minNormal = normal;
1186  isExtreme = (fExtremeFacets.find(&facet) != fExtremeFacets.end());
1187  }
1188  }
1189  }
1190  if (minDist < kInfinity)
1191  {
1192  aNormalVector = minNormal;
1193  aConvex = isExtreme;
1194  return minDist;
1195  }
1196  else
1197  {
1198  // No intersection found
1199  aConvex = false;
1200  Normal(p, aNormalVector);
1201  return 0.0;
1202  }
1203 }
1204 
1206 //
1208 DistanceToOutCandidates(const std::vector<G4int> &candidates,
1209  const G4ThreeVector &aPoint,
1210  const G4ThreeVector &direction,
1211  G4double &minDist, G4ThreeVector &minNormal,
1212  G4int &minCandidate ) const
1213 {
1214  G4int candidatesCount = candidates.size();
1215  G4double dist = 0.0;
1216  G4double distFromSurface = 0.0;
1218 
1219  for (G4int i = 0 ; i < candidatesCount; ++i)
1220  {
1221  G4int candidate = candidates[i];
1222  G4VFacet &facet = *fFacets[candidate];
1223  if (facet.Intersect(aPoint,direction,true,dist,distFromSurface,normal))
1224  {
1225  if (distFromSurface > 0.0 && distFromSurface <= kCarToleranceHalf
1226  && facet.Distance(aPoint,kCarTolerance) <= kCarToleranceHalf)
1227  {
1228  // We are on a surface
1229  //
1230  minDist = 0.0;
1231  minNormal = normal;
1232  minCandidate = candidate;
1233  break;
1234  }
1235  if (dist >= 0.0 && dist < minDist)
1236  {
1237  minDist = dist;
1238  minNormal = normal;
1239  minCandidate = candidate;
1240  }
1241  }
1242  }
1243 }
1244 
1246 //
1247 G4double
1249  const G4ThreeVector &aDirection,
1250  G4ThreeVector &aNormalVector,
1251  G4bool &aConvex,
1252  G4double aPstep) const
1253 {
1254  G4double minDistance;
1255 
1256  if (fVoxels.GetCountOfVoxels() > 1)
1257  {
1258  minDistance = kInfinity;
1259 
1260  G4ThreeVector currentPoint = aPoint;
1261  G4ThreeVector direction = aDirection.unit();
1262  G4double totalShift = 0;
1263  vector<G4int> curVoxel(3);
1264  if (!fVoxels.Contains(aPoint)) return 0;
1265 
1266  fVoxels.GetVoxel(curVoxel, currentPoint);
1267 
1268  G4double shiftBonus = kCarTolerance;
1269 
1270  const vector<G4int> *old = 0;
1271 
1272  G4int minCandidate = -1;
1273  do // Loop checking, 13.08.2015, G.Cosmo
1274  {
1275  const vector<G4int> &candidates = fVoxels.GetCandidates(curVoxel);
1276  if (old == &candidates)
1277  old++;
1278  if (old != &candidates && candidates.size())
1279  {
1280  DistanceToOutCandidates(candidates, aPoint, direction, minDistance,
1281  aNormalVector, minCandidate);
1282  if (minDistance <= totalShift) break;
1283  }
1284 
1285  G4double shift=fVoxels.DistanceToNext(currentPoint, direction, curVoxel);
1286  if (shift == kInfinity) break;
1287 
1288  totalShift += shift;
1289  if (minDistance <= totalShift) break;
1290 
1291  currentPoint += direction * (shift + shiftBonus);
1292 
1293  old = &candidates;
1294  }
1295  while (fVoxels.UpdateCurrentVoxel(currentPoint, direction, curVoxel));
1296 
1297  if (minCandidate < 0)
1298  {
1299  // No intersection found
1300  minDistance = 0;
1301  aConvex = false;
1302  Normal(aPoint, aNormalVector);
1303  }
1304  else
1305  {
1306  aConvex = (fExtremeFacets.find(fFacets[minCandidate])
1307  != fExtremeFacets.end());
1308  }
1309  }
1310  else
1311  {
1312  minDistance = DistanceToOutNoVoxels(aPoint, aDirection, aNormalVector,
1313  aConvex, aPstep);
1314  }
1315  return minDistance;
1316 }
1317 
1319 //
1321 DistanceToInCandidates(const std::vector<G4int> &candidates,
1322  const G4ThreeVector &aPoint,
1323  const G4ThreeVector &direction) const
1324 {
1325  G4int candidatesCount = candidates.size();
1326  G4double dist = 0.0;
1327  G4double distFromSurface = 0.0;
1329 
1330  G4double minDistance = kInfinity;
1331  for (G4int i = 0 ; i < candidatesCount; ++i)
1332  {
1333  G4int candidate = candidates[i];
1334  G4VFacet &facet = *fFacets[candidate];
1335  if (facet.Intersect(aPoint,direction,false,dist,distFromSurface,normal))
1336  {
1337  //
1338  // Set minDist to the new distance to current facet if distFromSurface is
1339  // in positive direction and point is not at surface. If the point is
1340  // within 0.5*kCarTolerance of the surface, then force distance to be
1341  // zero and leave member function immediately (for efficiency), as
1342  // proposed by & credit to Akira Okumura.
1343  //
1344  if ( (distFromSurface > kCarToleranceHalf)
1345  && (dist >= 0.0) && (dist < minDistance))
1346  {
1347  minDistance = dist;
1348  }
1349  else
1350  {
1351  if (-kCarToleranceHalf <= dist && dist <= kCarToleranceHalf)
1352  {
1353  return 0.0;
1354  }
1355  else if (distFromSurface > -kCarToleranceHalf
1356  && distFromSurface < kCarToleranceHalf)
1357  {
1358  minDistance = dist;
1359  }
1360  }
1361  }
1362  }
1363  return minDistance;
1364 }
1365 
1367 //
1368 G4double
1370  const G4ThreeVector &aDirection,
1371  G4double aPstep) const
1372 {
1373  G4double minDistance;
1374 
1375  if (fVoxels.GetCountOfVoxels() > 1)
1376  {
1377  minDistance = kInfinity;
1378  G4ThreeVector currentPoint = aPoint;
1379  G4ThreeVector direction = aDirection.unit();
1380  G4double shift = fVoxels.DistanceToFirst(currentPoint, direction);
1381  if (shift == kInfinity) return shift;
1382  G4double shiftBonus = kCarTolerance;
1383  if (shift)
1384  currentPoint += direction * (shift + shiftBonus);
1385  // if (!fVoxels.Contains(currentPoint)) return minDistance;
1386  G4double totalShift = shift;
1387 
1388  // G4SurfBits exclusion; // (1/*fVoxels.GetBitsPerSlice()*/);
1389  vector<G4int> curVoxel(3);
1390 
1391  fVoxels.GetVoxel(curVoxel, currentPoint);
1392  do // Loop checking, 13.08.2015, G.Cosmo
1393  {
1394  const vector<G4int> &candidates = fVoxels.GetCandidates(curVoxel);
1395  if (candidates.size())
1396  {
1397  G4double distance=DistanceToInCandidates(candidates, aPoint, direction);
1398  if (minDistance > distance) minDistance = distance;
1399  if (distance < totalShift) break;
1400  }
1401 
1402  shift = fVoxels.DistanceToNext(currentPoint, direction, curVoxel);
1403  if (shift == kInfinity /*|| shift == 0*/) break;
1404 
1405  totalShift += shift;
1406  if (minDistance < totalShift) break;
1407 
1408  currentPoint += direction * (shift + shiftBonus);
1409  }
1410  while (fVoxels.UpdateCurrentVoxel(currentPoint, direction, curVoxel));
1411  }
1412  else
1413  {
1414  minDistance = DistanceToInNoVoxels(aPoint, aDirection, aPstep);
1415  }
1416 
1417  return minDistance;
1418 }
1419 
1421 //
1422 G4bool
1423 G4TessellatedSolid::CompareSortedVoxel(const std::pair<G4int, G4double> &l,
1424  const std::pair<G4int, G4double> &r)
1425 {
1426  return l.second < r.second;
1427 }
1428 
1430 //
1431 G4double
1433  G4bool simple,
1434  G4VFacet * &minFacet) const
1435 {
1436  G4double minDist = kInfinity;
1437 
1438  G4int size = fVoxels.GetVoxelBoxesSize();
1439  vector<pair<G4int, G4double> > voxelsSorted(size);
1440 
1441  pair<G4int, G4double> info;
1442 
1443  for (G4int i = 0; i < size; ++i)
1444  {
1445  const G4VoxelBox &voxelBox = fVoxels.GetVoxelBox(i);
1446 
1447  G4ThreeVector pointShifted = p - voxelBox.pos;
1448  G4double safety = fVoxels.MinDistanceToBox(pointShifted, voxelBox.hlen);
1449  info.first = i;
1450  info.second = safety;
1451  voxelsSorted[i] = info;
1452  }
1453 
1454  std::sort(voxelsSorted.begin(), voxelsSorted.end(),
1456 
1457  for (G4int i = 0; i < size; ++i)
1458  {
1459  const pair<G4int,G4double> &inf = voxelsSorted[i];
1460  G4double dist = inf.second;
1461  if (dist > minDist) break;
1462 
1463  const vector<G4int> &candidates = fVoxels.GetVoxelBoxCandidates(inf.first);
1464  G4int csize = candidates.size();
1465  for (G4int j = 0; j < csize; ++j)
1466  {
1467  G4int candidate = candidates[j];
1468  G4VFacet &facet = *fFacets[candidate];
1469  dist = simple ? facet.Distance(p,minDist)
1470  : facet.Distance(p,minDist,false);
1471  if (dist < minDist)
1472  {
1473  minDist = dist;
1474  minFacet = &facet;
1475  }
1476  }
1477  }
1478  return minDist;
1479 }
1480 
1482 //
1484  G4bool aAccurate) const
1485 {
1486 #if G4SPECSDEBUG
1487  if ( Inside(p) == kInside )
1488  {
1489  std::ostringstream message;
1490  G4int oldprc = message.precision(16) ;
1491  message << "Point p is already inside!?" << G4endl
1492  << "Position:" << G4endl << G4endl
1493  << "p.x() = " << p.x()/mm << " mm" << G4endl
1494  << "p.y() = " << p.y()/mm << " mm" << G4endl
1495  << "p.z() = " << p.z()/mm << " mm" << G4endl
1496  << "DistanceToOut(p) == " << DistanceToOut(p);
1497  message.precision(oldprc) ;
1498  G4Exception("G4TriangularFacet::DistanceToIn(p)",
1499  "GeomSolids1002", JustWarning, message);
1500  }
1501 #endif
1502 
1503  G4double minDist;
1504 
1505  if (fVoxels.GetCountOfVoxels() > 1)
1506  {
1507  if (!aAccurate)
1508  return fVoxels.DistanceToBoundingBox(p);
1509 
1510  if (!OutsideOfExtent(p, kCarTolerance))
1511  {
1512  vector<G4int> startingVoxel(3);
1513  fVoxels.GetVoxel(startingVoxel, p);
1514  const vector<G4int> &candidates = fVoxels.GetCandidates(startingVoxel);
1515  if (candidates.size() == 0 && fInsides.GetNbits())
1516  {
1517  G4int index = fVoxels.GetPointIndex(p);
1518  if (fInsides[index]) return 0.;
1519  }
1520  }
1521 
1522  G4VFacet *facet;
1523  minDist = MinDistanceFacet(p, true, facet);
1524  }
1525  else
1526  {
1527  minDist = kInfinity;
1528  G4int size = fFacets.size();
1529  for (G4int i = 0; i < size; ++i)
1530  {
1531  G4VFacet &facet = *fFacets[i];
1532  G4double dist = facet.Distance(p,minDist);
1533  if (dist < minDist) minDist = dist;
1534  }
1535  }
1536  return minDist;
1537 }
1538 
1540 //
1541 G4double
1543 {
1544 #if G4SPECSDEBUG
1545  if ( Inside(p) == kOutside )
1546  {
1547  std::ostringstream message;
1548  G4int oldprc = message.precision(16) ;
1549  message << "Point p is already outside!?" << G4endl
1550  << "Position:" << G4endl << G4endl
1551  << "p.x() = " << p.x()/mm << " mm" << G4endl
1552  << "p.y() = " << p.y()/mm << " mm" << G4endl
1553  << "p.z() = " << p.z()/mm << " mm" << G4endl
1554  << "DistanceToIn(p) == " << DistanceToIn(p);
1555  message.precision(oldprc) ;
1556  G4Exception("G4TriangularFacet::DistanceToOut(p)",
1557  "GeomSolids1002", JustWarning, message);
1558  }
1559 #endif
1560 
1561  G4double minDist;
1562 
1563  if (OutsideOfExtent(p, kCarTolerance)) return 0.0;
1564 
1565  if (fVoxels.GetCountOfVoxels() > 1)
1566  {
1567  G4VFacet *facet;
1568  minDist = MinDistanceFacet(p, true, facet);
1569  }
1570  else
1571  {
1572  minDist = kInfinity;
1573  G4double dist = 0.0;
1574  G4int size = fFacets.size();
1575  for (G4int i = 0; i < size; ++i)
1576  {
1577  G4VFacet &facet = *fFacets[i];
1578  dist = facet.Distance(p,minDist);
1579  if (dist < minDist) minDist = dist;
1580  }
1581  }
1582  return minDist;
1583 }
1584 
1586 //
1587 // G4GeometryType GetEntityType() const;
1588 //
1589 // Provide identification of the class of an object
1590 //
1592 {
1593  return fGeometryType;
1594 }
1595 
1597 //
1598 std::ostream &G4TessellatedSolid::StreamInfo(std::ostream &os) const
1599 {
1600  os << G4endl;
1601  os << "Geometry Type = " << fGeometryType << G4endl;
1602  os << "Number of facets = " << fFacets.size() << G4endl;
1603 
1604  G4int size = fFacets.size();
1605  for (G4int i = 0; i < size; ++i)
1606  {
1607  os << "FACET # = " << i + 1 << G4endl;
1608  G4VFacet &facet = *fFacets[i];
1609  facet.StreamInfo(os);
1610  }
1611  os << G4endl;
1612 
1613  return os;
1614 }
1615 
1617 //
1618 // Make a clone of the object
1619 //
1621 {
1622  return new G4TessellatedSolid(*this);
1623 }
1624 
1626 //
1627 // EInside G4TessellatedSolid::Inside (const G4ThreeVector &p) const
1628 //
1629 // This method must return:
1630 // * kOutside if the point at offset p is outside the shape
1631 // boundaries plus kCarTolerance/2,
1632 // * kSurface if the point is <= kCarTolerance/2 from a surface, or
1633 // * kInside otherwise.
1634 //
1636 {
1637  EInside location;
1638 
1639  if (fVoxels.GetCountOfVoxels() > 1)
1640  {
1641  location = InsideVoxels(aPoint);
1642  }
1643  else
1644  {
1645  location = InsideNoVoxels(aPoint);
1646  }
1647  return location;
1648 }
1649 
1651 //
1653 {
1654  G4ThreeVector n;
1655  Normal(p, n);
1656  return n;
1657 }
1658 
1660 //
1661 // G4double DistanceToIn(const G4ThreeVector& p)
1662 //
1663 // Calculate distance to nearest surface of shape from an outside point p. The
1664 // distance can be an underestimate.
1665 //
1667 {
1668  return SafetyFromOutside(p,false);
1669 }
1670 
1672 //
1674  const G4ThreeVector& v)const
1675 {
1676  G4double dist = DistanceToInCore(p,v,kInfinity);
1677 #ifdef G4SPECSDEBUG
1678  if (dist < kInfinity)
1679  {
1680  if (Inside(p + dist*v) != kSurface)
1681  {
1682  std::ostringstream message;
1683  message << "Invalid response from facet in solid '" << GetName() << "',"
1684  << G4endl
1685  << "at point: " << p << "and direction: " << v;
1686  G4Exception("G4TessellatedSolid::DistanceToIn(p,v)",
1687  "GeomSolids1002", JustWarning, message);
1688  }
1689  }
1690 #endif
1691  return dist;
1692 }
1693 
1695 //
1696 // G4double DistanceToOut(const G4ThreeVector& p)
1697 //
1698 // Calculate distance to nearest surface of shape from an inside
1699 // point. The distance can be an underestimate.
1700 //
1702 {
1703  return SafetyFromInside(p,false);
1704 }
1705 
1707 //
1708 // G4double DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
1709 // const G4bool calcNorm=false,
1710 // G4bool *validNorm=0, G4ThreeVector *n=0);
1711 //
1712 // Return distance along the normalised vector v to the shape, from a
1713 // point at an offset p inside or on the surface of the
1714 // shape. Intersections with surfaces, when the point is not greater
1715 // than kCarTolerance/2 from a surface, must be ignored.
1716 // If calcNorm is true, then it must also set validNorm to either
1717 // * true, if the solid lies entirely behind or on the exiting
1718 // surface. Then it must set n to the outwards normal vector
1719 // (the Magnitude of the vector is not defined).
1720 // * false, if the solid does not lie entirely behind or on the
1721 // exiting surface.
1722 // If calcNorm is false, then validNorm and n are unused.
1723 //
1725  const G4ThreeVector& v,
1726  const G4bool calcNorm,
1727  G4bool *validNorm,
1728  G4ThreeVector *norm) const
1729 {
1730  G4ThreeVector n;
1731  G4bool valid;
1732 
1733  G4double dist = DistanceToOutCore(p, v, n, valid);
1734  if (calcNorm)
1735  {
1736  *norm = n;
1737  *validNorm = valid;
1738  }
1739 #ifdef G4SPECSDEBUG
1740  if (dist < kInfinity)
1741  {
1742  if (Inside(p + dist*v) != kSurface)
1743  {
1744  std::ostringstream message;
1745  message << "Invalid response from facet in solid '" << GetName() << "',"
1746  << G4endl
1747  << "at point: " << p << "and direction: " << v;
1748  G4Exception("G4TessellatedSolid::DistanceToOut(p,v,..)",
1749  "GeomSolids1002", JustWarning, message);
1750  }
1751  }
1752 #endif
1753  return dist;
1754 }
1755 
1757 //
1759 {
1760  scene.AddSolid (*this);
1761 }
1762 
1764 //
1766 {
1767  G4int nVertices = fVertexList.size();
1768  G4int nFacets = fFacets.size();
1769  G4PolyhedronArbitrary *polyhedron =
1770  new G4PolyhedronArbitrary (nVertices, nFacets);
1771  for (G4ThreeVectorList::const_iterator v= fVertexList.begin();
1772  v!=fVertexList.end(); ++v)
1773  {
1774  polyhedron->AddVertex(*v);
1775  }
1776 
1777  G4int size = fFacets.size();
1778  for (G4int i = 0; i < size; ++i)
1779  {
1780  G4VFacet &facet = *fFacets[i];
1781  G4int v[4];
1782  G4int n = facet.GetNumberOfVertices();
1783  if (n > 4) n = 4;
1784  else if (n == 3) v[3] = 0;
1785  for (G4int j=0; j<n; ++j)
1786  {
1787  G4int k = facet.GetVertexIndex(j);
1788  v[j] = k+1;
1789  }
1790  polyhedron->AddFacet(v[0],v[1],v[2],v[3]);
1791  }
1792  polyhedron->SetReferences();
1793 
1794  return (G4Polyhedron*) polyhedron;
1795 }
1796 
1798 //
1799 // GetPolyhedron
1800 //
1802 {
1803  if (!fpPolyhedron ||
1806  fpPolyhedron->GetNumberOfRotationSteps())
1807  {
1808  G4AutoLock l(&polyhedronMutex);
1809  delete fpPolyhedron;
1811  fRebuildPolyhedron = false;
1812  l.unlock();
1813  }
1814  return fpPolyhedron;
1815 }
1816 
1818 //
1819 // CalculateExtent
1820 //
1821 // Based on correction provided by Stan Seibert, University of Texas.
1822 //
1823 G4bool
1825  const G4VoxelLimits& pVoxelLimit,
1826  const G4AffineTransform& pTransform,
1827  G4double& pMin, G4double& pMax) const
1828 {
1829  G4ThreeVectorList transVertexList(fVertexList);
1830  G4int size = fVertexList.size();
1831 
1832  // Put solid into transformed frame
1833  for (G4int i=0; i < size; ++i)
1834  {
1835  pTransform.ApplyPointTransform(transVertexList[i]);
1836  }
1837 
1838  // Find min and max extent in each dimension
1840  G4ThreeVector maxExtent(-kInfinity, -kInfinity, -kInfinity);
1841 
1842  size = transVertexList.size();
1843  for (G4int i=0; i< size; ++i)
1844  {
1845  for (G4int axis=G4ThreeVector::X; axis < G4ThreeVector::SIZE; ++axis)
1846  {
1847  G4double coordinate = transVertexList[i][axis];
1848  if (coordinate < minExtent[axis])
1849  { minExtent[axis] = coordinate; }
1850  if (coordinate > maxExtent[axis])
1851  { maxExtent[axis] = coordinate; }
1852  }
1853  }
1854 
1855  // Check for containment and clamp to voxel boundaries
1856  for (G4int axis=G4ThreeVector::X; axis < G4ThreeVector::SIZE; ++axis)
1857  {
1858  EAxis geomAxis = kXAxis; // U geom classes use different index type
1859  switch(axis)
1860  {
1861  case G4ThreeVector::X: geomAxis = kXAxis; break;
1862  case G4ThreeVector::Y: geomAxis = kYAxis; break;
1863  case G4ThreeVector::Z: geomAxis = kZAxis; break;
1864  }
1865  G4bool isLimited = pVoxelLimit.IsLimited(geomAxis);
1866  G4double voxelMinExtent = pVoxelLimit.GetMinExtent(geomAxis);
1867  G4double voxelMaxExtent = pVoxelLimit.GetMaxExtent(geomAxis);
1868 
1869  if (isLimited)
1870  {
1871  if ( minExtent[axis] > voxelMaxExtent+kCarTolerance ||
1872  maxExtent[axis] < voxelMinExtent-kCarTolerance )
1873  {
1874  return false ;
1875  }
1876  else
1877  {
1878  if (minExtent[axis] < voxelMinExtent)
1879  {
1880  minExtent[axis] = voxelMinExtent ;
1881  }
1882  if (maxExtent[axis] > voxelMaxExtent)
1883  {
1884  maxExtent[axis] = voxelMaxExtent;
1885  }
1886  }
1887  }
1888  }
1889 
1890  // Convert pAxis into G4ThreeVector index
1891  G4int vecAxis=0;
1892  switch(pAxis)
1893  {
1894  case kXAxis: vecAxis = G4ThreeVector::X; break;
1895  case kYAxis: vecAxis = G4ThreeVector::Y; break;
1896  case kZAxis: vecAxis = G4ThreeVector::Z; break;
1897  default: break;
1898  }
1899 
1900  pMin = minExtent[vecAxis] - kCarTolerance;
1901  pMax = maxExtent[vecAxis] + kCarTolerance;
1902 
1903  return true;
1904 }
1905 
1907 //
1909 {
1910  return fMinExtent.x();
1911 }
1912 
1914 //
1916 {
1917  return fMaxExtent.x();
1918 }
1919 
1921 //
1923 {
1924  return fMinExtent.y();
1925 }
1926 
1928 //
1930 {
1931  return fMaxExtent.y();
1932 }
1933 
1935 //
1937 {
1938  return fMinExtent.z();
1939 }
1940 
1942 //
1944 {
1945  return fMaxExtent.z();
1946 }
1947 
1949 //
1951 {
1952  return G4VisExtent (fMinExtent.x(), fMaxExtent.x(), fMinExtent.y(), fMaxExtent.y(), fMinExtent.z(), fMaxExtent.z());
1953 }
1954 
1956 //
1958 {
1959  if (fCubicVolume != 0.) return fCubicVolume;
1960 
1961  // For explanation of the following algorithm see:
1962  // https://en.wikipedia.org/wiki/Polyhedron#Volume
1963  // http://wwwf.imperial.ac.uk/~rn/centroid.pdf
1964 
1965  G4int size = fFacets.size();
1966  for (G4int i = 0; i < size; ++i)
1967  {
1968  G4VFacet &facet = *fFacets[i];
1969  G4double area = facet.GetArea();
1970  G4ThreeVector unit_normal = facet.GetSurfaceNormal();
1971  fCubicVolume += area * (facet.GetVertex(0).dot(unit_normal));
1972  }
1973  fCubicVolume /= 3.;
1974  return fCubicVolume;
1975 }
1976 
1978 //
1980 {
1981  if (fSurfaceArea != 0.) return fSurfaceArea;
1982 
1983  G4int size = fFacets.size();
1984  for (G4int i = 0; i < size; ++i)
1985  {
1986  G4VFacet &facet = *fFacets[i];
1987  fSurfaceArea += facet.GetArea();
1988  }
1989  return fSurfaceArea;
1990 }
1991 
1993 //
1995 {
1996  // Select randomly a facet and return a random point on it
1997 
1998  G4int i = (G4int) G4RandFlat::shoot(0., fFacets.size());
1999  return fFacets[i]->GetPointOnFace();
2000 }
2001 
2003 //
2004 // SetRandomVectorSet
2005 //
2006 // This is a set of predefined random vectors (if that isn't a contradition
2007 // in terms!) used to generate rays from a user-defined point. The member
2008 // function Inside uses these to determine whether the point is inside or
2009 // outside of the tessellated solid. All vectors should be unit vectors.
2010 //
2012 {
2013  fRandir.resize(20);
2014  fRandir[0] =
2015  G4ThreeVector(-0.9577428892113370, 0.2732676269591740, 0.0897405271949221);
2016  fRandir[1] =
2017  G4ThreeVector(-0.8331264504940770,-0.5162067214954600,-0.1985722492445700);
2018  fRandir[2] =
2019  G4ThreeVector(-0.1516671651108820, 0.9666292616127460, 0.2064580868390110);
2020  fRandir[3] =
2021  G4ThreeVector( 0.6570250350323190,-0.6944539025883300, 0.2933460081893360);
2022  fRandir[4] =
2023  G4ThreeVector(-0.4820456281280320,-0.6331060000098690,-0.6056474264406270);
2024  fRandir[5] =
2025  G4ThreeVector( 0.7629032554236800 , 0.1016854697539910,-0.6384658864065180);
2026  fRandir[6] =
2027  G4ThreeVector( 0.7689540409061150, 0.5034929891988220, 0.3939600142169160);
2028  fRandir[7] =
2029  G4ThreeVector( 0.5765188359255740, 0.5997271636278330,-0.5549354566343150);
2030  fRandir[8] =
2031  G4ThreeVector( 0.6660632777862070,-0.6362809868288380, 0.3892379937580790);
2032  fRandir[9] =
2033  G4ThreeVector( 0.3824415020414780, 0.6541792713761380,-0.6525243125110690);
2034  fRandir[10] =
2035  G4ThreeVector(-0.5107726564526760, 0.6020905056811610, 0.6136760679616570);
2036  fRandir[11] =
2037  G4ThreeVector( 0.7459135439578050, 0.6618796061649330, 0.0743530220183488);
2038  fRandir[12] =
2039  G4ThreeVector( 0.1536405855311580, 0.8117477913978260,-0.5634359711967240);
2040  fRandir[13] =
2041  G4ThreeVector( 0.0744395301705579,-0.8707110101772920,-0.4861286795736560);
2042  fRandir[14] =
2043  G4ThreeVector(-0.1665874645185400, 0.6018553940549240,-0.7810369397872780);
2044  fRandir[15] =
2045  G4ThreeVector( 0.7766902003633100, 0.6014617505959970,-0.1870724331097450);
2046  fRandir[16] =
2047  G4ThreeVector(-0.8710128685847430,-0.1434320216603030,-0.4698551243971010);
2048  fRandir[17] =
2049  G4ThreeVector( 0.8901082092766820,-0.4388411398893870, 0.1229871120030100);
2050  fRandir[18] =
2051  G4ThreeVector(-0.6430417431544370,-0.3295938228697690, 0.6912779675984150);
2052  fRandir[19] =
2053  G4ThreeVector( 0.6331124368380410, 0.6306211461665000, 0.4488714875425340);
2054 
2055  fMaxTries = 20;
2056 }
2057 
2059 //
2061 {
2062  G4int base = sizeof(*this);
2063  base += fVertexList.capacity() * sizeof(G4ThreeVector);
2064  base += fRandir.capacity() * sizeof(G4ThreeVector);
2065 
2066  G4int limit = fFacets.size();
2067  for (G4int i = 0; i < limit; i++)
2068  {
2069  G4VFacet &facet = *fFacets[i];
2070  base += facet.AllocatedMemory();
2071  }
2072 
2073  std::set<G4VFacet *>::const_iterator beg, end, it;
2074  beg = fExtremeFacets.begin();
2075  end = fExtremeFacets.end();
2076  for (it = beg; it != end; it++)
2077  {
2078  G4VFacet &facet = *(*it);
2079  base += facet.AllocatedMemory();
2080  }
2081  return base;
2082 }
2083 
2085 //
2087 {
2089  G4int sizeInsides = fInsides.GetNbytes();
2090  G4int sizeVoxels = fVoxels.AllocatedMemory();
2091  size += sizeInsides + sizeVoxels;
2092  return size;
2093 }
G4bool Contains(const G4ThreeVector &point) const
void ResetBitNumber(unsigned int bitnumber)
Definition: G4SurfBits.hh:161
G4String GetName() const
unsigned int GetNbits() const
Definition: G4SurfBits.hh:101
ThreeVector shoot(const G4int Ap, const G4int Af)
void Clear()
Definition: G4SurfBits.cc:92
void SetSolidClosed(const G4bool t)
double Y(double density)
EInside InsideVoxels(const G4ThreeVector &aPoint) const
G4double DistanceToInCore(const G4ThreeVector &p, const G4ThreeVector &v, G4double aPstep=kInfinity) const
static const G4double kInfinity
Definition: geomdefs.hh:42
G4double DistanceToInNoVoxels(const G4ThreeVector &p, const G4ThreeVector &v, G4double aPstep=kInfinity) const
void Voxelize(std::vector< G4VFacet * > &facets)
virtual G4double GetArea()=0
CLHEP::Hep3Vector G4ThreeVector
void DistanceToOutCandidates(const std::vector< G4int > &candidates, const G4ThreeVector &aPoint, const G4ThreeVector &direction, G4double &minDist, G4ThreeVector &minNormal, G4int &minCandidate) const
virtual void DescribeYourselfTo(G4VGraphicsScene &scene) const
void CopyObjects(const G4TessellatedSolid &s)
std::vector< G4ThreeVector > fVertexList
virtual G4bool Intersect(const G4ThreeVector &, const G4ThreeVector &, const G4bool, G4double &, G4double &, G4ThreeVector &)=0
virtual G4VisExtent GetExtent() const
G4String name
Definition: TRTMaterials.hh:40
G4double DistanceToFirst(const G4ThreeVector &point, const G4ThreeVector &direction) const
G4double DistanceToBoundingBox(const G4ThreeVector &point) const
G4double GetMaxXExtent() const
virtual G4int GetNumberOfVertices() const =0
virtual G4ThreeVector GetCircumcentre() const =0
virtual G4double GetCubicVolume()
virtual G4double Distance(const G4ThreeVector &, G4double)=0
G4int GetVoxelsIndex(G4int x, G4int y, G4int z) const
G4double GetMinXExtent() const
G4double a
Definition: TRTMaterials.hh:39
std::vector< G4ThreeVector > fRandir
virtual G4double GetSurfaceArea()
static G4double MinDistanceToBox(const G4ThreeVector &aPoint, const G4ThreeVector &f)
virtual void SetVertexIndex(G4int i, G4int j)=0
G4bool UpdateCurrentVoxel(const G4ThreeVector &point, const G4ThreeVector &direction, std::vector< G4int > &curVoxel) const
G4double GetMaxZExtent() const
static G4bool CompareSortedVoxel(const std::pair< G4int, G4double > &l, const std::pair< G4int, G4double > &r)
virtual void AddSolid(const G4Box &)=0
int G4int
Definition: G4Types.hh:78
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:175
G4int GetPointIndex(const G4ThreeVector &p) const
void SetBitNumber(unsigned int bitnumber, G4bool value=true)
Definition: G4SurfBits.hh:123
G4VFacet * GetFacet(G4int i) const
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
G4bool IsEmpty(G4int index) const
virtual G4ThreeVector GetSurfaceNormal() const =0
G4double DistanceToInCandidates(const std::vector< G4int > &candidates, const G4ThreeVector &aPoint, const G4ThreeVector &aDirection) const
G4double GetMinZExtent() const
G4GLOB_DLL std::ostream G4cout
G4bool IsLimited() const
virtual G4double DistanceToOut(const G4ThreeVector &p) const
virtual G4Polyhedron * CreatePolyhedron() const
void GetVoxel(std::vector< G4int > &curVoxel, const G4ThreeVector &point) const
virtual void SetVertices(std::vector< G4ThreeVector > *vertices)=0
G4int SetAllUsingStack(const std::vector< G4int > &voxel, const std::vector< G4int > &max, G4bool status, G4SurfBits &checked)
virtual EInside Inside(const G4ThreeVector &p) const
bool G4bool
Definition: G4Types.hh:79
G4bool OutsideOfExtent(const G4ThreeVector &p, G4double tolerance=0) const
G4bool AddFacet(G4VFacet *aFacet)
G4double DistanceToNext(const G4ThreeVector &point, const G4ThreeVector &direction, const std::vector< G4int > &curVoxel) const
std::vector< G4ThreeVector > G4ThreeVectorList
Definition: G4VSolid.hh:79
void SetMaxVoxels(G4int max)
EInside InsideNoVoxels(const G4ThreeVector &p) const
G4TessellatedSolid & operator+=(const G4TessellatedSolid &right)
G4SurfaceVoxelizer fVoxels
const G4int n
const std::vector< G4int > & GetCandidates(std::vector< G4int > &curVoxel) const
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
void AddVertex(const G4ThreeVector &v)
virtual G4VSolid * Clone() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
std::set< G4VertexInfo, G4VertexComparator > fFacetList
virtual G4Polyhedron * GetPolyhedron() const
G4int G4Mutex
Definition: G4Threading.hh:173
G4int GetNumberOfFacets() const
G4bool IsInside(const G4ThreeVector &p) const
Definition: G4VFacet.cc:114
virtual std::ostream & StreamInfo(std::ostream &os) const
virtual G4int GetVertexIndex(G4int i) const =0
G4ThreeVector pos
T max(const T t1, const T t2)
brief Return the largest of the two arguments
EInside
Definition: geomdefs.hh:58
G4double MinDistanceFacet(const G4ThreeVector &p, G4bool simple, G4VFacet *&facet) const
EAxis
Definition: geomdefs.hh:54
G4Polyhedron * fpPolyhedron
G4double DistanceToOutNoVoxels(const G4ThreeVector &p, const G4ThreeVector &v, G4ThreeVector &aNormalVector, G4bool &aConvex, G4double aPstep=kInfinity) const
virtual G4double SafetyFromInside(const G4ThreeVector &p, G4bool aAccurate=false) const
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4VFacet.cc:98
void AddFacet(const G4int iv1, const G4int iv2, const G4int iv3, const G4int iv4=0)
G4GeometryType fGeometryType
std::set< G4VFacet * > fExtremeFacets
std::vector< G4VFacet * > fFacets
long long GetCountOfVoxels() const
virtual G4GeometryType GetEntityType() const
const std::vector< G4double > & GetBoundary(G4int index) const
#define G4endl
Definition: G4ios.hh:61
G4VSolid & operator=(const G4VSolid &rhs)
Definition: G4VSolid.cc:110
G4double kCarTolerance
Definition: G4VSolid.hh:304
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
const G4SurfBits & Empty() const
G4int GetMaxVoxels(G4ThreeVector &ratioOfReduction)
const std::vector< G4int > & GetVoxelBoxCandidates(G4int i) const
double G4double
Definition: G4Types.hh:76
virtual G4VFacet * GetClone()=0
virtual G4bool Normal(const G4ThreeVector &p, G4ThreeVector &n) const
G4int GetVoxelBoxesSize() const
G4double GetMaxExtent(const EAxis pAxis) const
G4ThreeVector hlen
virtual G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
virtual G4ThreeVector GetVertex(G4int i) const =0
static const double mm
Definition: G4SIunits.hh:114
const G4VoxelBox & GetVoxelBox(G4int i) const
G4double GetMaxYExtent() const
virtual G4ThreeVector GetPointOnSurface() const
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
G4double DistanceToOutCore(const G4ThreeVector &p, const G4ThreeVector &v, G4ThreeVector &aNormalVector, G4bool &aConvex, G4double aPstep=kInfinity) const
static const G4double pos
G4double GetMinExtent(const EAxis pAxis) const
virtual G4int AllocatedMemory()=0
virtual G4double SafetyFromOutside(const G4ThreeVector &p, G4bool aAccurate=false) const
G4double GetMinYExtent() const
unsigned int GetNbytes() const
Definition: G4SurfBits.hh:102
G4TessellatedSolid & operator=(const G4TessellatedSolid &right)
G4bool GetSolidClosed() const
void ApplyPointTransform(G4ThreeVector &vec) const
virtual G4bool IsDefined() const =0