62 : fBoundingBox(
"VoxBBox", 1, 1, 1)
85 const std::vector<G4int> empty(0);
88 unsigned int size =
max[0] *
max[1] *
max[2];
94 for (xyz[2] = 0; xyz[2] < max[2]; ++xyz[2])
96 for (xyz[1] = 0; xyz[1] < max[1]; ++xyz[1])
98 for (xyz[0] = 0; xyz[0] < max[0]; ++xyz[0])
109 std::vector<G4int> &c = (
fCandidates[index] = empty);
110 c.reserve(candidates.size());
111 c.assign(candidates.begin(), candidates.end());
123 std::vector<G4Transform3D>& transforms)
134 if (
G4int numNodes = solids.size())
144 for (
G4int i = 0; i < numNodes; ++i)
156 orbToleranceVector.
set(tolerance,tolerance,tolerance);
157 min -= orbToleranceVector;
158 max += orbToleranceVector;
162 min -= toleranceVector;
163 max += toleranceVector;
182 if (
G4int numNodes = facets.size())
189 for (
G4int i = 0; i < numNodes; ++i)
197 min -= toleranceVector;
198 max += toleranceVector;
201 fBoxes[i].pos = min + hlen;
214 for(
G4int i = 0; i < numNodes; ++i)
216 G4cout << setw(10) << setiosflags(ios::fixed) <<
217 " -> Node " << i+1 <<
":\n" <<
218 "\t * [x,y,z] = " <<
fBoxes[i].hlen <<
219 "\t * [x,y,z] = " <<
fBoxes[i].pos <<
"\n";
221 G4cout.precision(oldprec);
236 for(
G4int i = 0 ; i < numNodes; ++i)
248 boundary[2*i] = p -
d;
249 boundary[2*i+1] = p +
d;
251 std::sort(boundary.begin(), boundary.end());
272 std::vector<G4double> sortedBoundary(2*numNodes);
276 for (
G4int j = 0; j <= 2; ++j)
284 for(
G4int i = 0 ; i < 2*numNodes; ++i)
286 G4double newBoundary = sortedBoundary[i];
288 if (j == 0)
G4cout <<
"Examining " << newBoundary <<
"..." <<
G4endl;
290 G4int size = boundary.size();
291 if(!size || std::abs(boundary[size-1] - newBoundary) > tolerance)
296 if (j == 0)
G4cout <<
"Adding boundary " << newBoundary <<
"..."
299 boundary.push_back(newBoundary);
307 G4int n = boundary.size();
311 G4int skip = n / (max /2);
313 std::vector<G4double> reduced;
314 for (
G4int i = 0; i <
n; ++i)
317 G4int size = boundary.size();
318 if (i % skip == 0 || i == 0 || i == size - 1)
325 reduced.push_back(boundary[i]);
337 char axis[3] = {
'X',
'Y',
'Z'};
338 for (
G4int i = 0; i <= 2; ++i)
340 G4cout <<
" * " << axis[i] <<
" axis:" <<
G4endl <<
" | ";
350 G4int count = boundaries.size();
352 for(
G4int i = 0; i < count; ++i)
354 G4cout << setw(10) << setiosflags(ios::fixed) << boundaries[i];
355 if(i != count-1)
G4cout <<
"-> ";
358 G4cout.precision(oldprec);
371 for (
G4int k = 0; k < 3; ++k)
374 std::vector<G4double> &boundary = boundaries[k];
375 G4int voxelsCount = boundary.size() - 1;
384 bitmask.
SetBitNumber(voxelsCount*bitsPerSlice-1,
false);
389 candidatesCount.resize(voxelsCount);
391 for(
G4int i = 0 ; i < voxelsCount; ++i) { candidatesCount[i] = 0; }
395 for(
G4int j = 0 ; j < numNodes; ++j)
406 if (i < 0) { i = 0; }
414 candidatesCount[i]++;
418 while (max > boundary[i] && i < voxelsCount);
434 for(
G4int i=0; i<numNodes; ++i)
446 char axis[3] = {
'X',
'Y',
'Z'};
450 for (
G4int j = 0; j <= 2; ++j)
454 for(
G4int i=0; i < count-1; ++i)
483 for (
G4int i = 0; i <= 2; ++i)
516 if (maxVoxels > 0 && maxVoxels < maxTotal)
519 ratio = std::pow(ratio, 1./3.);
520 if (ratio > 1) { ratio = 1; }
521 reductionRatio.
set(ratio,ratio,ratio);
529 for (
G4int k = 0; k <= 2; ++k)
533 std::vector<G4VoxelInfo> voxels(max);
535 std::set<G4int, G4VoxelComparator> voxelSet(comp);
536 std::vector<G4int> mergings;
541 voxel.
count = candidatesCount[j];
547 for (
G4int j = 0; j < max - 1; ++j) { voxelSet.insert(j); }
550 G4double reduction = reductionRatio[k];
553 G4int count = 0, currentCount;
554 while ((currentCount = voxelSet.size()) > 2)
557 if ((currentRatio <= reduction) && (currentCount <= 1000))
559 const G4int pos = *voxelSet.begin();
560 mergings.push_back(pos + 1);
565 if (voxelSet.erase(pos) != 1)
569 if (voxel.
next != max - 1)
570 if (voxelSet.erase(voxel.
next) != 1)
575 if (voxelSet.erase(voxel.
previous) != 1)
583 if (voxel.
next != max - 1)
584 voxelSet.insert(voxel.
next);
597 std::sort(mergings.begin(), mergings.end());
599 const std::vector<G4double>& boundary = boundaries[k];
600 int mergingsSize = mergings.size();
601 vector<G4double> reducedBoundary;
602 G4int skip = mergings[0], i = 0;
603 max = boundary.size();
608 reducedBoundary.push_back(boundary[j]);
610 else if (++i < mergingsSize)
615 boundaries[k] = reducedBoundary;
684 for (
G4int k = 0; k <= 2; ++k)
689 for (
G4int i = 0; i <
max; ++i) total += candidatesCount[i];
691 G4double reduction = reductionRatio[k];
695 G4int destination = (
G4int) (reduction * max) + 1;
696 if (destination > 1000) destination = 1000;
697 if (destination < 2) destination = 2;
700 std::vector<G4int> mergings;
702 std::vector<G4double> &boundary = boundaries[k];
703 std::vector<G4double> reducedBoundary(destination);
708 sum += candidatesCount[i];
709 if (sum > average * (cur + 1) || i == 0)
712 reducedBoundary[cur] = val;
714 if (cur == destination)
718 reducedBoundary[destination-1] = boundary[
max];
719 boundaries[k] = reducedBoundary;
725 std::vector<G4Transform3D>& transforms)
735 for (
G4int i = 0; i < 3; ++i)
745 std::vector<G4int> voxel(3), maxVoxels(3);
746 for (
G4int i = 0; i <= 2; ++i) maxVoxels[i] = boundaries[i].size();
749 for (voxel[2] = 0; voxel[2] < maxVoxels[2] - 1; ++voxel[2])
751 for (voxel[1] = 0; voxel[1] < maxVoxels[1] - 1; ++voxel[1])
753 for (voxel[0] = 0; voxel[0] < maxVoxels[0] - 1; ++voxel[0])
755 std::vector<G4int> candidates;
760 for (
G4int i = 0; i <= 2; ++i)
762 G4int index = voxel[i];
763 const std::vector<G4double> &boundary = boundaries[i];
764 G4double hlen = 0.5 * (boundary[index+1] - boundary[index]);
766 box.
pos[i] = boundary[index] + hlen;
769 std::vector<G4int>(candidates).swap(candidates);
783 G4int size = facets.size();
786 for (
G4int i = 0; i < (
G4int) facets.size(); ++i)
788 if (facets[i]->GetNumberOfVertices() > 3) size++;
792 if ((size >= 10 || maxVoxels > 0) && maxVoxels != 0 && maxVoxels != 1)
831 G4cout <<
"Total number of voxels after reduction: "
847 std::vector<G4double> miniBoundaries[3];
864 G4cout <<
"Total number of mini voxels: " << total <<
G4endl;
892 G4cout <<
"Deallocating unnecessary fields during runtime..." <<
G4endl;
897 for (
G4int i = 0; i < 3; ++i)
912 G4cout <<
" Candidates in voxel [" << voxels[0] <<
" ; " << voxels[1]
913 <<
" ; " << voxels[2] <<
"]: ";
914 std::vector<G4int> candidates;
917 for (
G4int i = 0; i < count; ++i)
G4cout << candidates[i];
923 std::vector<G4int> &list,
G4int i)
925 for (
G4int byte = 0; byte < (
G4int) (
sizeof(
unsigned int)); byte++)
927 if (
G4int maskByte = mask & 0xFF)
929 for (
G4int bit = 0; bit < 8; bit++)
932 { list.push_back(8*(
sizeof(
unsigned int)*i+ byte) + bit); }
933 if (!(maskByte >>= 1))
break;
965 for (
G4int i = 0 ; i < limit; ++i)
972 if (current.
x() > max.
x()) max.
setX(current.
x());
973 if (current.
x() < min.
x()) min.
setX(current.
x());
975 if (current.
y() > max.
y()) max.
setY(current.
y());
976 if (current.
y() < min.
y()) min.
setY(current.
y());
978 if (current.
z() > max.
z()) max.
setZ(current.
z());
979 if (current.
z() < min.
z()) min.
setZ(current.
z());
985 std::vector<G4int> &list,
G4SurfBits *crossed)
const
991 for (
G4int i = 0; i <= 2; ++i)
1006 unsigned int mask = 0xFFffFFff;
1011 if (!(mask = ((
unsigned int*)
fBitmasks[0].fAllBits)[slice]))
1017 if (!(mask &= ((
unsigned int*)
fBitmasks[1].fAllBits)[slice]))
1023 if (!(mask &= ((
unsigned int*)
fBitmasks[2].fAllBits)[slice]))
1026 if (crossed && (!(mask &= ~((
unsigned int*)crossed->
fAllBits)[0])))
1033 unsigned int* masks[3],
mask;
1034 for (
G4int i = 0; i <= 2; ++i)
1037 masks[i] = ((
unsigned int*)
fBitmasks[i].fAllBits)
1040 unsigned int* maskCrossed = crossed
1041 ? (
unsigned int*)crossed->
fAllBits : 0;
1048 if (!(mask = masks[0][i]))
continue;
1049 if (!(mask &= masks[1][i]))
continue;
1050 if (!(mask &= masks[2][i]))
continue;
1051 if (maskCrossed && !(mask &= ~maskCrossed[i]))
continue;
1107 std::vector<G4int> &list,
1122 if (!(mask = ((
unsigned int *) bitmasks[0].fAllBits)[voxels[0]]))
1124 if (!(mask &= ((
unsigned int *) bitmasks[1].fAllBits)[voxels[1]]))
1126 if (!(mask &= ((
unsigned int *) bitmasks[2].fAllBits)[voxels[2]]))
1128 if (crossed && (!(mask &= ~((
unsigned int *)crossed->
fAllBits)[0])))
1135 unsigned int *masks[3],
mask;
1136 for (
G4int i = 0; i <= 2; ++i)
1138 masks[i] = ((
unsigned int *) bitmasks[i].fAllBits)
1141 unsigned int *maskCrossed = crossed
1142 ? (
unsigned int *)crossed->
fAllBits : 0;
1149 if (!(mask = masks[0][i]))
continue;
1150 if (!(mask &= masks[1][i]))
continue;
1151 if (!(mask &= masks[2][i]))
continue;
1152 if (maskCrossed && !(mask &= ~maskCrossed[i]))
continue;
1164 std::vector<G4int> &list,
G4SurfBits *crossed)
const
1174 for (
G4int i = 0; i < 3; ++i)
1211 safe = safx = -f.
x() + std::abs(aPoint.
x());
1212 safy = -f.
y() + std::abs(aPoint.
y());
1213 if ( safy > safe ) safe = safy;
1214 safz = -f.
z() + std::abs(aPoint.
z());
1215 if ( safz > safe ) safe = safz;
1216 if (safe < 0.0)
return 0.0;
1220 if ( safx > 0 ) { safsq += safx*safx; count++; }
1221 if ( safy > 0 ) { safsq += safy*safy; count++; }
1222 if ( safz > 0 ) { safsq += safz*safz; count++; }
1223 if (count == 1)
return safe;
1224 return std::sqrt(safsq);
1231 std::vector<G4int> &curVoxel)
const
1236 for (
int i = 0; i <= 2; ++i)
1240 const std::vector<G4double>& boundary =
fBoundaries[i];
1241 G4int index = curVoxel[i];
1242 if (direction[i] >= 1
e-10)
1248 if (direction[i] > -1
e-10)
1251 G4double dif = boundary[index] - point[i];
1252 G4double distance = dif / direction[i];
1254 if (shift > distance)
1266 if (direction[cur] > 0)
1273 if (--curVoxel[cur] < 0)
1316 std::vector<G4int> &curVoxel)
const
1318 for (
G4int i = 0; i <= 2; ++i)
1320 G4int index = curVoxel[i];
1321 const std::vector<G4double> &boundary =
fBoundaries[i];
1323 if (direction[i] > 0)
1325 if (point[i] >= boundary[++index])
1326 if (++curVoxel[i] >= (
G4int) boundary.size() - 1)
1331 if (point[i] < boundary[index])
1332 if (--curVoxel[i] < 0)
1337 if (curVoxel[i] != indexOK)
1338 curVoxel[i] = indexOK;
1383 for (
G4int i = 0; i < csize; ++i)
void SetXHalfLength(G4double dx)
void set(double x, double y, double z)
G4double total(Particle const *const p1, Particle const *const p2)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
CLHEP::Hep3Vector G4ThreeVector
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4ThreeVector fBoundingBoxSize
std::vector< G4double > fBoundaries[3]
void ResetBitNumber(unsigned int bitnumber)
void Voxelize(std::vector< G4VSolid * > &solids, std::vector< G4Transform3D > &transforms)
G4bool Contains(const G4ThreeVector &point) const
static const G4double pos
static const G4double kInfinity
void BuildReduceVoxels(std::vector< G4double > fBoundaries[], G4ThreeVector reductionRatio)
G4ThreeVector GetGlobalPoint(const G4Transform3D &trans, const G4ThreeVector &lpoint) const
G4ThreeVector fReductionRatio
G4int GetCandidatesVoxelArray(const G4ThreeVector &point, std::vector< G4int > &list, G4SurfBits *crossed=0) const
virtual void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
std::map< G4int, std::vector< G4int > > fCandidates
unsigned int GetNbytes() const
void BuildVoxelLimits(std::vector< G4VSolid * > &solids, std::vector< G4Transform3D > &transforms)
void SetZHalfLength(G4double dz)
void DisplayVoxelLimits() const
void BuildBitmasks(std::vector< G4double > fBoundaries[], G4SurfBits bitmasks[], G4bool countsOnly=false)
std::vector< G4int > fCandidatesCounts[3]
void SetMaxVoxels(G4int max)
void DisplayListNodes() const
void SetReductionRatio(G4int maxVoxels, G4ThreeVector &reductionRatio)
G4double DistanceToFirst(const G4ThreeVector &point, const G4ThreeVector &direction) const
std::vector< std::vector< G4int > > fVoxelBoxesCandidates
G4String GetCandidatesAsString(const G4SurfBits &bits) const
std::vector< G4VoxelBox > fVoxelBoxes
void GetCandidatesVoxel(std::vector< G4int > &voxels)
void ResetAllBits(G4bool value=false)
static G4double MinDistanceToBox(const G4ThreeVector &aPoint, const G4ThreeVector &f)
static void FindComponentsFastest(unsigned int mask, std::vector< G4int > &list, G4int i)
void set(unsigned int nbits, const char *array)
void CreateSortedBoundary(std::vector< G4double > &boundaryRaw, G4int axis)
G4bool TestBitNumber(unsigned int bitnumber) const
G4double G4ParticleHPJENDLHEData::G4double result
long long CountVoxels(std::vector< G4double > boundaries[]) const
void SetYHalfLength(G4double dy)
static G4SolidStore * GetInstance()
G4double DistanceToNext(const G4ThreeVector &point, const G4ThreeVector &direction, std::vector< G4int > &curVoxel) const
std::vector< G4VoxelBox > fBoxes
G4GLOB_DLL std::ostream G4cout
G4int GetBitsPerSlice() const
G4ThreeVector fBoundingBoxCenter
G4bool UpdateCurrentVoxel(const G4ThreeVector &point, const G4ThreeVector &direction, std::vector< G4int > &curVoxel) const
static G4ThreadLocal G4int fDefaultVoxelsCount
static G4GeometryTolerance * GetInstance()
void SetBitNumber(unsigned int bitnumber, G4bool value=true)
G4double DistanceToBoundingBox(const G4ThreeVector &point) const
void BuildReduceVoxels2(std::vector< G4double > fBoundaries[], G4ThreeVector reductionRatio)
static void DeRegister(G4VSolid *pSolid)
void TransformLimits(G4ThreeVector &min, G4ThreeVector &max, const G4Transform3D &transformation) const
virtual G4GeometryType GetEntityType() const =0
G4int GetVoxelsIndex(G4int x, G4int y, G4int z) const
void CreateMiniVoxels(std::vector< G4double > fBoundaries[], G4SurfBits bitmasks[])
static G4int BinarySearch(const std::vector< T > &vec, T value)
virtual G4double Extent(const G4ThreeVector)=0
G4double GetSurfaceTolerance() const
static G4int GetDefaultVoxelsCount()
static void SetDefaultVoxelsCount(G4int count)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double GetRadialTolerance() const