155 if (&ts ==
this)
return *
this;
207 for (
G4int i = 0; i <
n; ++i)
226 G4Exception(
"G4TessellatedSolid::AddFacet()",
"GeomSolids1002",
227 JustWarning,
"Attempt to add facets when solid is closed.");
232 set<G4VertexInfo,G4VertexComparator>::iterator begin
237 value.
mag2 = p.x() + p.y() + p.z();
246 while (!found && it != end)
251 if ((found = (facet == aFacet)))
break;
253 if (dif > kCarTolerance3)
break;
260 while (!found && it != begin)
266 found = (facet == aFacet);
269 if (dif > kCarTolerance3)
break;
284 G4Exception(
"G4TessellatedSolid::AddFacet()",
"GeomSolids1002",
285 JustWarning,
"Attempt to add facet not properly defined.");
294 const std::vector<G4int> &
max,
297 vector<G4int> xyz = voxel;
298 stack<vector<G4int> >
pos;
301 G4int cc = 0, nz = 0;
303 vector<G4int> candidates;
321 for (
G4int i = 0; i <= 2; ++i)
323 if (xyz[i] < max[i] - 1)
351 vector<G4int> voxel(3), maxVoxels(3);
353 G4int size = maxVoxels[0] * maxVoxels[1] * maxVoxels[2];
360 for (voxel[2] = 0; voxel[2] < maxVoxels[2] - 1; ++voxel[2])
362 for (voxel[1] = 0; voxel[1] < maxVoxels[1] - 1; ++voxel[1])
364 for (voxel[0] = 0; voxel[0] < maxVoxels[0] - 1; ++voxel[0])
412 for (
G4int j = 0; j < size; ++j)
418 for (
G4int i=0; i < vsize; ++i)
444 set<G4VertexInfo,G4VertexComparator> vertexListSorted;
445 set<G4VertexInfo,G4VertexComparator>::iterator begin
446 = vertexListSorted.begin(), end = vertexListSorted.end(),
pos, it;
455 vector<G4int> newIndex(100);
457 for (
G4int k = 0; k < size; ++k)
466 value.
mag2 = p.x() + p.y() + p.z();
472 pos = vertexListSorted.lower_bound(value);
479 found = (dif < kCarTolerance24);
481 dif = q.x() + q.y() + q.z() - value.
mag2;
482 if (dif > kCarTolerance3)
break;
495 found = (dif < kCarTolerance24);
497 dif = value.
mag2 - (q.x() + q.y() + q.z());
498 if (dif > kCarTolerance3)
break;
506 G4cout << p.x() <<
":" << p.y() <<
":" << p.z() <<
G4endl;
507 G4cout <<
"Adding new vertex #" << i <<
" of facet " << k
512 vertexListSorted.insert(value);
513 begin = vertexListSorted.begin();
514 end = vertexListSorted.end();
515 newIndex[i] = value.
id;
533 G4cout << p.x() <<
":" << p.y() <<
":" << p.z() <<
G4endl;
534 G4cout <<
"Vertex #" << i <<
" of facet " << k
535 <<
" found, redirecting to " <<
id <<
G4endl;
551 for (set<G4VertexInfo,G4VertexComparator>::iterator res=
552 vertexListSorted.begin(); res!=vertexListSorted.end(); ++res)
554 G4int id = (*res).id;
556 G4double mvalue = vec.x() + vec.y() + vec.z();
557 if (previousValue && (previousValue - 1e-9 > mvalue))
558 G4cout <<
"Error in CreateVertexList: previousValue " << previousValue
559 <<
" is smaller than mvalue " << mvalue <<
G4endl;
560 previousValue = mvalue;
572 G4cout <<
"G4TessellatedSolid - Allocated memory without voxel overhead "
573 << without <<
"; with " << with <<
"; ratio: " << ratio <<
G4endl;
629 for (
G4int i = 0; i < size; ++i)
655 vector<G4int> startingVoxel(3);
658 const G4double dirTolerance = 1.0E-14;
660 const vector<G4int> &startingCandidates =
662 G4int limit = startingCandidates.size();
672 for(
G4int i = 0; i < limit; ++i)
674 G4int candidate = startingCandidates[i];
677 if (dist < minDist) minDist = dist;
704 G4bool nearParallel =
false;
724 vector<G4int> curVoxel(3);
725 curVoxel = startingVoxel;
733 const vector<G4int> &candidates =
736 if (
G4int candidatesCount = candidates.size())
738 for (
G4int i = 0 ; i < candidatesCount; ++i)
740 G4int candidate = candidates[i];
744 crossingO = facet.
Intersect(p,v,
true,distO,distFromSurfaceO,normalO);
745 crossingI = facet.
Intersect(p,v,
false,distI,distFromSurfaceI,normalI);
747 if (crossingO || crossingI)
751 nearParallel = (crossingO
752 && std::fabs(normalO.dot(v))<dirTolerance)
753 || (crossingI && std::fabs(normalI.dot(v))<dirTolerance);
756 if (crossingO && distO > 0.0 && distO < distOut)
758 if (crossingI && distI > 0.0 && distI < distIn)
764 if (nearParallel)
break;
780 currentPoint += direction * (shift + shiftBonus);
799 std::ostringstream message;
800 G4int oldprc = message.precision(16);
801 message <<
"Cannot determine whether point is inside or outside volume!"
807 <<
"p.x() = " << p.x()/
mm <<
" mm" <<
G4endl
808 <<
"p.y() = " << p.y()/
mm <<
" mm" <<
G4endl
809 <<
"p.z() = " << p.z()/
mm <<
" mm";
810 message.precision(oldprc);
846 const G4double dirTolerance = 1.0E-14;
853 for (
G4int i = 0; i < size; ++i)
857 if (dist < minDist) minDist = dist;
893 for (
G4int i=0; i<nTry; ++i)
895 G4bool nearParallel =
false;
907 vector<G4VFacet*>::const_iterator f =
fFacets.begin();
916 crossingO = ((*f)->Intersect(p,v,
true,distO,distFromSurfaceO,normalO));
917 crossingI = ((*f)->Intersect(p,v,
false,distI,distFromSurfaceI,normalI));
918 if (crossingO || crossingI)
920 nearParallel = (crossingO && std::fabs(normalO.dot(v))<dirTolerance)
921 || (crossingI && std::fabs(normalI.dot(v))<dirTolerance);
924 if (crossingO && distO > 0.0 && distO < distOut) distOut = distO;
925 if (crossingI && distI > 0.0 && distI < distIn) distIn = distI;
928 }
while (!nearParallel && ++f!=
fFacets.end());
939 std::ostringstream message;
940 G4int oldprc = message.precision(16);
941 message <<
"Cannot determine whether point is inside or outside volume!"
947 <<
"p.x() = " << p.x()/
mm <<
" mm" <<
G4endl
948 <<
"p.y() = " << p.y()/
mm <<
" mm" <<
G4endl
949 <<
"p.z() = " << p.z()/
mm <<
" mm";
950 message.precision(oldprc);
972 if (i == 0) location = locationprime;
991 vector<G4int> curVoxel(3);
996 if (
G4int limit = candidates.size())
999 for(
G4int i = 0 ; i < limit ; ++i)
1001 G4int candidate = candidates[i];
1004 if (dist < minDist) minDist = dist;
1018 for (
G4int i = 0; i < size; ++i)
1038 std::ostringstream message;
1039 message <<
"Point p is not on surface !?" <<
G4endl
1040 <<
" No facets found for point: " << p <<
" !" <<
G4endl
1041 <<
" Returning approximated value for normal.";
1043 G4Exception(
"G4TessellatedSolid::SurfaceNormal(p)",
1074 std::ostringstream message;
1075 G4int oldprc = message.precision(16) ;
1076 message <<
"Point p is already inside!?" <<
G4endl
1078 <<
" p.x() = " << p.x()/
mm <<
" mm" <<
G4endl
1079 <<
" p.y() = " << p.y()/
mm <<
" mm" <<
G4endl
1080 <<
" p.z() = " << p.z()/
mm <<
" mm" <<
G4endl
1082 message.precision(oldprc) ;
1083 G4Exception(
"G4TriangularFacet::DistanceToIn(p,v)",
1089 for (
G4int i = 0; i < size; ++i)
1092 if (facet.
Intersect(p,v,
false,dist,distFromSurface,normal))
1142 std::ostringstream message;
1143 G4int oldprc = message.precision(16) ;
1144 message <<
"Point p is already outside!?" <<
G4endl
1146 <<
" p.x() = " << p.x()/
mm <<
" mm" <<
G4endl
1147 <<
" p.y() = " << p.y()/
mm <<
" mm" <<
G4endl
1148 <<
" p.z() = " << p.z()/
mm <<
" mm" <<
G4endl
1150 message.precision(oldprc) ;
1151 G4Exception(
"G4TriangularFacet::DistanceToOut(p)",
1156 G4bool isExtreme =
false;
1158 for (
G4int i = 0; i < size; ++i)
1161 if (facet.
Intersect(p,v,
true,dist,distFromSurface,normal))
1173 if (dist >= 0.0 && dist < minDist)
1183 aNormalVector = minNormal;
1184 aConvex = isExtreme;
1191 Normal(p, aNormalVector);
1203 G4int &minCandidate )
const
1205 G4int candidatesCount = candidates.size();
1210 for (
G4int i = 0 ; i < candidatesCount; ++i)
1212 G4int candidate = candidates[i];
1214 if (facet.
Intersect(aPoint,direction,
true,dist,distFromSurface,normal))
1223 minCandidate = candidate;
1226 if (dist >= 0.0 && dist < minDist)
1230 minCandidate = candidate;
1254 vector<G4int> curVoxel(3);
1261 const vector<G4int> *old = 0;
1263 G4int minCandidate = -1;
1267 if (old == &candidates)
1269 if (old != &candidates && candidates.size())
1272 aNormalVector, minCandidate);
1273 if (minDistance <= totalShift)
break;
1279 totalShift += shift;
1280 if (minDistance <= totalShift)
break;
1282 currentPoint += direction * (shift + shiftBonus);
1288 if (minCandidate < 0)
1293 Normal(aPoint, aNormalVector);
1316 G4int candidatesCount = candidates.size();
1322 for (
G4int i = 0 ; i < candidatesCount; ++i)
1324 G4int candidate = candidates[i];
1326 if (facet.
Intersect(aPoint,direction,
false,dist,distFromSurface,normal))
1336 && (dist >= 0.0) && (dist < minDistance))
1375 currentPoint += direction * (shift + shiftBonus);
1380 vector<G4int> curVoxel(3);
1386 if (candidates.size())
1389 if (minDistance > distance) minDistance = distance;
1390 if (distance < totalShift)
break;
1396 totalShift += shift;
1397 if (minDistance < totalShift)
break;
1399 currentPoint += direction * (shift + shiftBonus);
1415 const std::pair<G4int, G4double> &r)
1417 return l.second < r.second;
1430 vector<pair<G4int, G4double> > voxelsSorted(size);
1432 pair<G4int, G4double> info;
1434 for (
G4int i = 0; i < size; ++i)
1441 info.second = safety;
1442 voxelsSorted[i] = info;
1445 std::sort(voxelsSorted.begin(), voxelsSorted.end(),
1448 for (
G4int i = 0; i < size; ++i)
1450 const pair<G4int,G4double> &inf = voxelsSorted[i];
1452 if (dist > minDist)
break;
1455 G4int csize = candidates.size();
1456 for (
G4int j = 0; j < csize; ++j)
1458 G4int candidate = candidates[j];
1460 dist = simple ? facet.
Distance(p,minDist)
1480 std::ostringstream message;
1481 G4int oldprc = message.precision(16) ;
1482 message <<
"Point p is already inside!?" <<
G4endl
1484 <<
"p.x() = " << p.x()/
mm <<
" mm" <<
G4endl
1485 <<
"p.y() = " << p.y()/
mm <<
" mm" <<
G4endl
1486 <<
"p.z() = " << p.z()/
mm <<
" mm" <<
G4endl
1488 message.precision(oldprc) ;
1503 vector<G4int> startingVoxel(3);
1520 for (
G4int i = 0; i < size; ++i)
1524 if (dist < minDist) minDist = dist;
1538 std::ostringstream message;
1539 G4int oldprc = message.precision(16) ;
1540 message <<
"Point p is already outside!?" <<
G4endl
1542 <<
"p.x() = " << p.x()/
mm <<
" mm" <<
G4endl
1543 <<
"p.y() = " << p.y()/
mm <<
" mm" <<
G4endl
1544 <<
"p.z() = " << p.z()/
mm <<
" mm" <<
G4endl
1546 message.precision(oldprc) ;
1547 G4Exception(
"G4TriangularFacet::DistanceToOut(p)",
1566 for (
G4int i = 0; i < size; ++i)
1570 if (dist < minDist) minDist = dist;
1596 for (
G4int i = 0; i < size; ++i)
1598 os <<
"FACET # = " << i + 1 <<
G4endl;
1733 for (G4ThreeVectorList::const_iterator v=
fVertexList.begin();
1740 for (
G4int i = 0; i < size; ++i)
1746 else if (n == 3) v[3] = 0;
1747 for (
G4int j=0; j<
n; ++j)
1752 polyhedron->
AddFacet(v[0],v[1],v[2],v[3]);
1791 for (
G4int i=0; i < size; ++i)
1800 size = transVertexList.size();
1801 for (
G4int i=0; i< size; ++i)
1803 for (
G4int axis=G4ThreeVector::X; axis < G4ThreeVector::SIZE; ++axis)
1805 G4double coordinate = transVertexList[i][axis];
1806 if (coordinate < minExtent[axis])
1807 { minExtent[axis] = coordinate; }
1808 if (coordinate > maxExtent[axis])
1809 { maxExtent[axis] = coordinate; }
1814 for (
G4int axis=G4ThreeVector::X; axis < G4ThreeVector::SIZE; ++axis)
1819 case G4ThreeVector::X: geomAxis =
kXAxis;
break;
1820 case G4ThreeVector::Y: geomAxis =
kYAxis;
break;
1821 case G4ThreeVector::Z: geomAxis =
kZAxis;
break;
1836 if (minExtent[axis] < voxelMinExtent)
1838 minExtent[axis] = voxelMinExtent ;
1840 if (maxExtent[axis] > voxelMaxExtent)
1842 maxExtent[axis] = voxelMaxExtent;
1852 case kXAxis: vecAxis = G4ThreeVector::X;
break;
1853 case kYAxis: vecAxis = G4ThreeVector::Y;
break;
1854 case kZAxis: vecAxis = G4ThreeVector::Z;
break;
1929 for (
G4int i = 0; i < size; ++i)
1944 return fFacets[i]->GetPointOnFace();
1960 G4ThreeVector(-0.9577428892113370, 0.2732676269591740, 0.0897405271949221);
1962 G4ThreeVector(-0.8331264504940770,-0.5162067214954600,-0.1985722492445700);
1964 G4ThreeVector(-0.1516671651108820, 0.9666292616127460, 0.2064580868390110);
1966 G4ThreeVector( 0.6570250350323190,-0.6944539025883300, 0.2933460081893360);
1968 G4ThreeVector(-0.4820456281280320,-0.6331060000098690,-0.6056474264406270);
1970 G4ThreeVector( 0.7629032554236800 , 0.1016854697539910,-0.6384658864065180);
1972 G4ThreeVector( 0.7689540409061150, 0.5034929891988220, 0.3939600142169160);
1974 G4ThreeVector( 0.5765188359255740, 0.5997271636278330,-0.5549354566343150);
1976 G4ThreeVector( 0.6660632777862070,-0.6362809868288380, 0.3892379937580790);
1978 G4ThreeVector( 0.3824415020414780, 0.6541792713761380,-0.6525243125110690);
1980 G4ThreeVector(-0.5107726564526760, 0.6020905056811610, 0.6136760679616570);
1982 G4ThreeVector( 0.7459135439578050, 0.6618796061649330, 0.0743530220183488);
1984 G4ThreeVector( 0.1536405855311580, 0.8117477913978260,-0.5634359711967240);
1986 G4ThreeVector( 0.0744395301705579,-0.8707110101772920,-0.4861286795736560);
1988 G4ThreeVector(-0.1665874645185400, 0.6018553940549240,-0.7810369397872780);
1990 G4ThreeVector( 0.7766902003633100, 0.6014617505959970,-0.1870724331097450);
1992 G4ThreeVector(-0.8710128685847430,-0.1434320216603030,-0.4698551243971010);
1994 G4ThreeVector( 0.8901082092766820,-0.4388411398893870, 0.1229871120030100);
1996 G4ThreeVector(-0.6430417431544370,-0.3295938228697690, 0.6912779675984150);
1998 G4ThreeVector( 0.6331124368380410, 0.6306211461665000, 0.4488714875425340);
2007 G4int base =
sizeof(*this);
2012 for (
G4int i = 0; i < limit; i++)
2018 std::set<G4VFacet *>::const_iterator beg, end, it;
2021 for (it = beg; it != end; it++)
2036 size += sizeInsides + sizeVoxels;
G4bool Contains(const G4ThreeVector &point) const
void ResetBitNumber(unsigned int bitnumber)
unsigned int GetNbits() const
ThreeVector shoot(const G4int Ap, const G4int Af)
void SetSolidClosed(const G4bool t)
EInside InsideVoxels(const G4ThreeVector &aPoint) const
G4double DistanceToInCore(const G4ThreeVector &p, const G4ThreeVector &v, G4double aPstep=kInfinity) const
static const G4double kInfinity
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 PrecalculateInsides()
void CopyObjects(const G4TessellatedSolid &s)
std::vector< G4ThreeVector > fVertexList
G4int GetCandidates(std::vector< G4int > &curVoxel, std::vector< G4int > *&candidates, std::vector< G4int > &space) const
virtual G4bool Intersect(const G4ThreeVector &, const G4ThreeVector &, const G4bool, G4double &, G4double &, G4ThreeVector &)=0
virtual G4VisExtent GetExtent() const
G4double DistanceToFirst(const G4ThreeVector &point, const G4ThreeVector &direction) const
G4double DistanceToBoundingBox(const G4ThreeVector &point) const
virtual G4double GetCubicVolume()
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
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
G4int GetPointIndex(const G4ThreeVector &p) const
void SetBitNumber(unsigned int bitnumber, G4bool value=true)
G4VFacet * GetFacet(G4int i) const
static double normal(HepRandomEngine *eptr)
G4double kCarToleranceHalf
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
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
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
void SetMaxVoxels(G4int max)
EInside InsideNoVoxels(const G4ThreeVector &p) const
G4TessellatedSolid & operator+=(const G4TessellatedSolid &right)
G4SurfaceVoxelizer fVoxels
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)
std::set< G4VertexInfo, G4VertexComparator > fFacetList
G4int AllocatedMemoryWithoutVoxels()
virtual G4Polyhedron * GetPolyhedron() const
G4int GetNumberOfFacets() const
G4bool IsInside(const G4ThreeVector &p) const
virtual std::ostream & StreamInfo(std::ostream &os) const
virtual G4int GetVertexIndex(G4int i) const =0
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double MinDistanceFacet(const G4ThreeVector &p, G4bool simple, G4VFacet *&facet) const
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
void AddFacet(const G4int iv1, const G4int iv2, const G4int iv3, const G4int iv4=0)
virtual ~G4TessellatedSolid()
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
G4VSolid & operator=(const G4VSolid &rhs)
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
const G4SurfBits & Empty() const
G4int GetMaxVoxels(G4ThreeVector &ratioOfReduction)
const std::vector< G4int > & GetVoxelBoxCandidates(G4int i) const
void DisplayAllocatedMemory()
virtual G4VFacet * GetClone()=0
virtual G4bool Normal(const G4ThreeVector &p, G4ThreeVector &n) const
G4int GetVoxelBoxesSize() const
G4double GetMaxExtent(const EAxis pAxis) const
virtual G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
virtual G4ThreeVector GetVertex(G4int i) const =0
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
G4TessellatedSolid & operator=(const G4TessellatedSolid &right)
G4bool GetSolidClosed() const
virtual G4bool IsDefined() const =0