61 : fBoundingBox(
"TessBBox", 1, 1, 1)
84 const vector<G4int> empty(0);
87 unsigned int size =
max[0] *
max[1] *
max[2];
93 for (xyz[2] = 0; xyz[2] < max[2]; ++xyz[2])
95 for (xyz[1] = 0; xyz[1] < max[1]; ++xyz[1])
97 for (xyz[0] = 0; xyz[0] < max[0]; ++xyz[0])
109 c.reserve(candidates.size());
110 c.assign(candidates.begin(), candidates.end());
128 if (
G4int numNodes = facets.size())
136 for (
G4int i = 0; i < numNodes; ++i)
144 min -= toleranceVector; max += toleranceVector;
147 fBoxes[i].pos = min + hlen;
160 for(
G4int i = 0; i < numNodes; ++i)
162 G4cout << setw(10) << setiosflags(ios::fixed) <<
163 " -> Node " << i+1 <<
":\n" <<
164 "\t * [x,y,z] = " <<
fBoxes[i].hlen <<
165 "\t * [x,y,z] = " <<
fBoxes[i].pos <<
"\n";
167 G4cout.precision(oldprec);
182 for(
G4int i = 0 ; i < numNodes; ++i)
192 G4cout <<
"Boundary " << p - d <<
" - " << p + d <<
G4endl;
194 boundary[2*i] = p - d;
195 boundary[2*i+1] = p + d;
197 sort(boundary.begin(), boundary.end());
218 vector<G4double> sortedBoundary(2*numNodes);
222 for (
G4int j = 0; j <= 2; ++j)
230 for(
G4int i = 0 ; i < 2*numNodes; ++i)
232 G4double newBoundary = sortedBoundary[i];
234 if (j == 0)
G4cout <<
"Examining " << newBoundary <<
"..." <<
G4endl;
236 G4int size = boundary.size();
237 if(!size || abs(boundary[size-1] - newBoundary) > tolerance)
242 if (j == 0)
G4cout <<
"Adding boundary " << newBoundary <<
"..."
245 boundary.push_back(newBoundary);
253 G4int n = boundary.size();
257 G4int skip = n / (max /2);
259 vector<G4double> reduced;
260 for (
G4int i = 0; i <
n; i++)
263 G4int size = boundary.size();
264 if (i % skip == 0 || i == 0 || i == size - 1)
271 reduced.push_back(boundary[i]);
283 char axis[3] = {
'X',
'Y',
'Z'};
284 for (
G4int i = 0; i <= 2; ++i)
286 G4cout <<
" * " << axis[i] <<
" axis:" <<
G4endl <<
" | ";
296 G4int count = boundaries.size();
298 for(
G4int i = 0; i < count; ++i)
300 G4cout << setw(10) << setiosflags(ios::fixed) << boundaries[i];
301 if(i != count-1)
G4cout <<
"-> ";
304 G4cout.precision(oldprec);
317 for (
G4int k = 0; k < 3; k++)
320 vector<G4double> &boundary = boundaries[k];
321 G4int voxelsCount = boundary.size() - 1;
330 bitmask.
SetBitNumber(voxelsCount*bitsPerSlice-1,
false);
335 candidatesCount.resize(voxelsCount);
337 for(
G4int i = 0 ; i < voxelsCount; ++i) { candidatesCount[i] = 0; }
341 for(
G4int j = 0 ; j < numNodes; j++)
352 if (i < 0) { i = 0; }
358 G4int index = i*bitsPerSlice+j;
361 candidatesCount[i]++;
365 while (max > boundary[i] && i < voxelsCount);
381 for(
G4int i=0; i<numNodes; ++i)
393 char axis[3] = {
'X',
'Y',
'Z'};
397 for (
G4int j = 0; j <= 2; ++j)
401 for(
G4int i=0; i < count-1; ++i)
418 for (
G4int i = 0; i <= 2; ++i)
451 if (maxVoxels > 0 && maxVoxels < maxTotal)
454 ratio = pow(ratio, 1./3.);
455 if (ratio > 1) { ratio = 1; }
456 reductionRatio.set(ratio,ratio,ratio);
464 for (
G4int k = 0; k <= 2; ++k)
468 vector<G4VoxelInfo> voxels(max);
470 set<G4int, G4VoxelComparator> voxelSet(comp);
471 vector<G4int> mergings;
476 voxel.
count = candidatesCount[j];
482 for (
G4int j = 0; j < max - 1; ++j) { voxelSet.insert(j); }
488 G4double reduction = reductionRatio[k];
491 G4int currentCount = voxelSet.size();
492 if (currentCount <= 2)
495 if (currentRatio <= reduction && currentCount <= 1000)
497 const G4int pos = *voxelSet.begin();
498 mergings.push_back(pos);
504 if (voxel.
next != max - 1) { voxelSet.erase(voxel.
next); }
511 if (voxel.
next != max - 1)
512 voxelSet.insert(voxel.
next);
524 sort(mergings.begin(), mergings.end());
526 vector<G4double> &boundary = boundaries[k];
527 vector<G4double> reducedBoundary(boundary.size() - mergings.size());
528 G4int skip = mergings[0] + 1, cur = 0, i = 0;
529 max = boundary.size();
534 reducedBoundary[cur++] = boundary[j];
538 if (++i < (
G4int)mergings.size()) { skip = mergings[i] + 1; }
541 boundaries[k] = reducedBoundary;
550 for (
G4int k = 0; k <= 2; ++k)
555 for (
G4int i = 0; i <
max; i++) total += candidatesCount[i];
557 G4double reduction = reductionRatio[k];
561 G4int destination = (
G4int) (reduction * max) + 1;
562 if (destination > 1000) destination = 1000;
563 if (destination < 2) destination = 2;
566 vector<G4int> mergings;
568 vector<G4double> &boundary = boundaries[k];
569 vector<G4double> reducedBoundary(destination);
571 G4int sum = 0, cur = 0;
574 sum += candidatesCount[i];
575 if (sum > average * (cur + 1) || i == 0)
578 reducedBoundary[cur] = val;
580 if (cur == destination)
584 reducedBoundary[destination-1] = boundary[
max];
585 boundaries[k] = reducedBoundary;
593 vector<G4int> voxel(3), maxVoxels(3);
594 for (
G4int i = 0; i <= 2; ++i) maxVoxels[i] = boundaries[i].size();
597 for (voxel[2] = 0; voxel[2] < maxVoxels[2] - 1; ++voxel[2])
599 for (voxel[1] = 0; voxel[1] < maxVoxels[1] - 1; ++voxel[1])
601 for (voxel[0] = 0; voxel[0] < maxVoxels[0] - 1; ++voxel[0])
603 vector<G4int> candidates;
608 for (
G4int i = 0; i <= 2; ++i)
610 G4int index = voxel[i];
611 const vector<G4double> &boundary = boundaries[i];
612 G4double hlen = 0.5 * (boundary[index+1] - boundary[index]);
614 box.
pos[i] = boundary[index] + hlen;
617 vector<G4int>(candidates).swap(candidates);
631 G4int size = facets.size();
634 for (
G4int i = 0; i < (
G4int) facets.size(); i++)
636 if (facets[i]->GetNumberOfVertices() > 3) size++;
640 if ((size >= 10 || maxVoxels > 0) && maxVoxels != 0 && maxVoxels != 1)
679 G4cout <<
"Total number of voxels after reduction: "
695 vector<G4double> miniBoundaries[3];
715 G4cout <<
"Total number of mini voxels: " << total <<
G4endl;
743 G4cout <<
"Deallocating unnecessary fields during runtime..." <<
G4endl;
748 for (
G4int i = 0; i < 3; i++)
763 G4cout <<
" Candidates in voxel [" << voxels[0] <<
" ; " << voxels[1]
764 <<
" ; " << voxels[2] <<
"]: ";
765 vector<G4int> candidates;
768 for (
G4int i = 0; i < count; ++i)
G4cout << candidates[i];
774 std::vector<G4int> &list,
G4int i)
776 for (
G4int byte = 0; byte < (
G4int) (
sizeof(
unsigned int)); byte++)
778 if (
G4int maskByte = mask & 0xFF)
780 for (
G4int bit = 0; bit < 8; bit++)
783 { list.push_back(8*(
sizeof(
unsigned int)*i+ byte) + bit); }
784 if (!(maskByte >>= 1))
break;
793 std::vector<G4int> &list,
G4SurfBits *crossed)
const
799 for (
G4int i = 0; i <= 2; ++i)
816 if (!(mask = ((
unsigned int *)
fBitmasks[0].fAllBits)[slice]
819 if (!(mask &= ((
unsigned int *)
fBitmasks[1].fAllBits)[slice]
822 if (!(mask &= ((
unsigned int *)
fBitmasks[2].fAllBits)[slice]
824 if (crossed && (!(mask &= ~((
unsigned int *)crossed->
fAllBits)[0])))
831 unsigned int *masks[3], mask;
832 for (
G4int i = 0; i <= 2; ++i)
837 unsigned int *maskCrossed =
838 crossed ? (
unsigned int *)crossed->
fAllBits : 0;
845 if (!(mask = masks[0][i]))
continue;
846 if (!(mask &= masks[1][i]))
continue;
847 if (!(mask &= masks[2][i]))
continue;
848 if (maskCrossed && !(mask &= ~maskCrossed[i]))
continue;
861 std::vector<G4int> &list,
876 if (!(mask = ((
unsigned int *) bitmasks[0].fAllBits)[voxels[0]]
878 if (!(mask &= ((
unsigned int *) bitmasks[1].fAllBits)[voxels[1]]
880 if (!(mask &= ((
unsigned int *) bitmasks[2].fAllBits)[voxels[2]]
882 if (crossed && (!(mask &= ~((
unsigned int *)crossed->
fAllBits)[0])))
889 unsigned int *masks[3], mask;
890 for (
G4int i = 0; i <= 2; ++i)
892 masks[i] = ((
unsigned int *) bitmasks[i].fAllBits)
895 unsigned int *maskCrossed =
896 crossed ? (
unsigned int *)crossed->
fAllBits : 0;
903 if (!(mask = masks[0][i]))
continue;
904 if (!(mask &= masks[1][i]))
continue;
905 if (!(mask &= masks[2][i]))
continue;
906 if (maskCrossed && !(mask &= ~maskCrossed[i]))
continue;
918 std::vector<G4int> &list,
G4SurfBits *crossed)
const
928 for (
G4int i = 0; i < 3; ++i)
965 safe = safx = -f.x() + std::abs(aPoint.x());
966 safy = -f.y() + std::abs(aPoint.y());
967 if ( safy > safe ) safe = safy;
968 safz = -f.z() + std::abs(aPoint.z());
969 if ( safz > safe ) safe = safz;
970 if (safe < 0.0)
return 0.0;
974 if ( safx > 0 ) { safsq += safx*safx; count++; }
975 if ( safy > 0 ) { safsq += safy*safy; count++; }
976 if ( safz > 0 ) { safsq += safz*safz; count++; }
977 if (count == 1)
return safe;
978 return std::sqrt(safsq);
985 const std::vector<G4int> &curVoxel)
const
989 for (
G4int i = 0; i <= 2; ++i)
994 G4int cur = curVoxel[i];
995 if(direction[i] >= 1e-10)
998 if (++cur >= (
G4int) boundary.size())
1003 if(direction[i] <= -1e-10)
1012 G4double dif = boundary[cur] - point[i];
1013 G4double distance = dif / direction[i];
1015 if (shift > distance)
1026 std::vector<G4int> &curVoxel)
const
1028 for (
G4int i = 0; i <= 2; ++i)
1030 G4int index = curVoxel[i];
1031 const vector<G4double> &boundary =
fBoundaries[i];
1033 if (direction[i] > 0)
1035 if (point[i] >= boundary[++index])
1036 if (++curVoxel[i] >= (
G4int) boundary.size() - 1)
1041 if (point[i] < boundary[index])
1042 if (--curVoxel[i] < 0)
1047 if (curVoxel[i] != indexOK)
1048 curVoxel[i] = indexOK;
1095 for (
G4int i = 0; i < csize; i++)
G4bool Contains(const G4ThreeVector &point) const
void ResetBitNumber(unsigned int bitnumber)
G4ThreeVector fBoundingBoxSize
void BuildReduceVoxels(std::vector< G4double > fBoundaries[], G4ThreeVector reductionRatio)
void DisplayListNodes() const
static const G4double kInfinity
G4int GetCandidatesVoxelArray(const G4ThreeVector &point, std::vector< G4int > &list, G4SurfBits *crossed=0) const
std::vector< G4VoxelBox > fVoxelBoxes
void Voxelize(std::vector< G4VFacet * > &facets)
CLHEP::Hep3Vector G4ThreeVector
static void DeRegister(G4VSolid *pSolid)
G4ThreeVector fBoundingBoxCenter
G4double DistanceToFirst(const G4ThreeVector &point, const G4ThreeVector &direction) const
G4double GetSurfaceTolerance() const
G4double DistanceToBoundingBox(const G4ThreeVector &point) const
std::vector< G4int > fCandidatesCounts[3]
void BuildVoxelLimits(std::vector< G4VFacet * > &facets)
G4int GetVoxelsIndex(G4int x, G4int y, G4int z) const
G4String GetCandidatesAsString(const G4SurfBits &bits) const
void DisplayVoxelLimits()
static G4double MinDistanceToBox(const G4ThreeVector &aPoint, const G4ThreeVector &f)
G4bool UpdateCurrentVoxel(const G4ThreeVector &point, const G4ThreeVector &direction, std::vector< G4int > &curVoxel) const
G4int GetBitsPerSlice() const
static G4ThreadLocal G4int fDefaultVoxelsCount
void SetBitNumber(unsigned int bitnumber, G4bool value=true)
std::map< G4int, std::vector< G4int > > fCandidates
static G4int BinarySearch(const std::vector< T > &vec, T value)
G4bool TestBitNumber(unsigned int bitnumber) const
void SetReductionRatio(G4int maxVoxels, G4ThreeVector &reductionRatio)
G4GLOB_DLL std::ostream G4cout
std::vector< std::vector< G4int > > fVoxelBoxesCandidates
void CreateMiniVoxels(std::vector< G4double > fBoundaries[], G4SurfBits bitmasks[])
void BuildBitmasks(std::vector< G4double > fBoundaries[], G4SurfBits bitmasks[], G4bool countsOnly=false)
std::vector< G4VoxelBox > fBoxes
G4double DistanceToNext(const G4ThreeVector &point, const G4ThreeVector &direction, const std::vector< G4int > &curVoxel) const
static G4SolidStore * GetInstance()
void SetMaxVoxels(G4int max)
std::vector< G4double > fBoundaries[3]
void ResetAllBits(G4bool value=false)
void set(unsigned int nbits, const char *array)
static G4int SetDefaultVoxelsCount(G4int count)
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4double total(Particle const *const p1, Particle const *const p2)
virtual G4double Extent(const G4ThreeVector)=0
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4ThreeVector fReductionRatio
static void FindComponentsFastest(unsigned int mask, std::vector< G4int > &list, G4int i)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void BuildReduceVoxels2(std::vector< G4double > fBoundaries[], G4ThreeVector reductionRatio)
void GetCandidatesVoxel(std::vector< G4int > &voxels)
void CreateSortedBoundary(std::vector< G4double > &boundaryRaw, G4int axis)
long long CountVoxels(std::vector< G4double > boundaries[]) const
static const G4double pos
static G4GeometryTolerance * GetInstance()
static G4int GetDefaultVoxelsCount()
unsigned int GetNbytes() const