40 const treeBoundBox& bb,
41 FixedList<DynamicList<label>, 8>& dividedIndices
44 const label presize = (indices.size()/8);
46 for (
direction octant = 0; octant < 8; ++octant)
48 const treeBoundBox subBbs(bb.subBbox(octant));
50 auto& contains = dividedIndices[octant];
53 contains.reserve_nocopy(presize);
55 for (
const label index : indices)
57 if (shapes_.overlaps(index, subBbs))
59 contains.push_back(index);
70 const treeBoundBox& bb,
72 const label parentNodeIndex,
73 const label octantToBeDivided
78 bb.min().x() >= bb.max().x()
79 || bb.min().y() >= bb.max().y()
80 || bb.min().z() >= bb.max().z()
84 <<
"Badly formed bounding box:" << bb
93 const DynamicList<label>& indices = contents_[contentIndex];
95 FixedList<DynamicList<label>, 8> dividedIndices;
96 divide(indices, bb, dividedIndices);
102 bool replaceNode =
true;
104 for (
direction octant = 0; octant < 8; ++octant)
106 auto& subIndices = dividedIndices[octant];
108 if (subIndices.size())
113 contents_[contentIndex] = std::move(subIndices);
119 contentIndex = contents_.size();
120 contents_.push_back(std::move(subIndices));
123 nod.subNodes_[octant] = contentPlusOctant(contentIndex, octant);
128 nod.subNodes_[octant] = emptyPlusOctant(octant);
133 if (parentNodeIndex != -1)
135 nod.parent_ = parentNodeIndex;
137 const label newNodeId = nodes_.size();
139 nodes_.push_back(nod);
141 nodes_[parentNodeIndex].subNodes_[octantToBeDivided]
142 = nodePlusOctant(newNodeId, octantToBeDivided);
152 const treeBoundBox& subBb,
153 const label contentI,
154 const label parentIndex,
161 contents_[contentI].size() > minSize_
162 && nLevels < maxLevels_
166 node nod =
divide(subBb, contentI, parentIndex, octant);
173 for (
direction subOct = 0; subOct < node::nChildren; ++subOct)
175 const labelBits subNodeLabel = nod.subNodes_[subOct];
177 if (isContent(subNodeLabel))
179 const treeBoundBox subBb = nod.bb_.subBbox(subOct);
181 const label subContentI = getContent(subNodeLabel);
183 const label parentNodeIndex = nodes_.size() - 1;
209 const node& nod = nodes_[nodeI];
211 volumeType myType = volumeType::UNKNOWN;
213 for (
direction octant = 0; octant < node::nChildren; ++octant)
217 labelBits index = nod.subNodes_[octant];
222 subType = calcVolumeType(getNode(index));
224 else if (isContent(index))
234 const treeBoundBox subBb = nod.bb_.subBbox(octant);
238 shapes_.getVolumeType(*
this, subBb.centre())
243 nodeTypes_.set(labelBits::pack(nodeI, octant), subType);
247 if (myType == volumeType::UNKNOWN)
251 else if (subType != myType)
267 const node& nod = nodes_[nodeI];
269 direction octant = nod.bb_.subOctant(sample);
271 volumeType octantType =
274 nodeTypes_.get(labelBits::pack(nodeI, octant))
277 if (octantType == volumeType::INSIDE)
281 else if (octantType == volumeType::OUTSIDE)
285 else if (octantType == volumeType::UNKNOWN)
292 labelBits index = nod.subNodes_[octant];
297 volumeType subType = getVolumeType(getNode(index), sample);
301 else if (isContent(index))
304 return volumeType(shapes_.getVolumeType(*
this, sample));
310 <<
"Sample:" << sample <<
" node:" << nodeI
311 <<
" with bb:" << nodes_[nodeI].bb_ <<
nl 312 <<
"Empty subnode has invalid volume type MIXED." 315 return volumeType::UNKNOWN;
319 <<
"Sample:" << sample <<
" at node:" << nodeI
320 <<
" octant:" << octant
321 <<
" with bb:" << nod.bb_.subBbox(octant) <<
nl 322 <<
"Node has invalid volume type " << octantType
325 return volumeType::UNKNOWN;
332 const vector& outsideNormal,
336 if ((outsideNormal&vec) >= 0)
338 return volumeType::OUTSIDE;
342 return volumeType::INSIDE;
353 scalar& nearestDistSqr,
354 label& nearestShapeI,
358 const node& nod = nodes_[nodeI];
363 for (
const direction octant : nod.bb_.searchOrder(sample))
365 labelBits index = nod.subNodes_[octant];
369 label subNodeI = getNode(index);
371 const treeBoundBox& subBb = nodes_[subNodeI].bb_;
373 if (subBb.overlaps(sample, nearestDistSqr))
386 else if (isContent(index))
388 if (nod.bb_.subOverlaps(octant, sample, nearestDistSqr))
392 contents_[getContent(index)],
411 treeBoundBox& tightest,
412 label& nearestShapeI,
417 const node& nod = nodes_[nodeI];
418 const treeBoundBox& nodeBb = nod.bb_;
423 for (
const direction octant : nod.bb_.searchOrder(
ln.centre()))
425 labelBits index = nod.subNodes_[octant];
429 const treeBoundBox& subBb = nodes_[getNode(index)].bb_;
431 if (subBb.overlaps(tightest))
445 else if (isContent(index))
447 if (nodeBb.subOverlaps(octant, tightest))
451 contents_[getContent(index)],
468 const label parentNodeI,
473 const node& nod = nodes_[parentNodeI];
474 labelBits index = nod.subNodes_[octant];
479 return nodes_[getNode(index)].bb_;
483 return nod.bb_.
subBbox(octant);
490 const treeBoundBox& bb,
492 const bool pushInside
499 const vector perturbVec = perturbTol_*bb.span();
501 point perturbedPt(pt);
510 if (
mag(pt[dir]-bb.min()[dir]) <
mag(perturbVec[dir]))
513 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
514 perturbedPt[dir] = bb.min()[dir] + perturbDist;
516 else if (
mag(pt[dir]-bb.max()[dir]) <
mag(perturbVec[dir]))
519 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
520 perturbedPt[dir] = bb.max()[dir] - perturbDist;
528 if (
mag(pt[dir]-bb.min()[dir]) <
mag(perturbVec[dir]))
530 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
531 perturbedPt[dir] = bb.min()[dir] - perturbDist;
533 else if (
mag(pt[dir]-bb.max()[dir]) <
mag(perturbVec[dir]))
535 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
536 perturbedPt[dir] = bb.max()[dir] + perturbDist;
543 if (pushInside != bb.contains(perturbedPt))
546 <<
"pushed point:" << pt
547 <<
" to:" << perturbedPt
548 <<
" wanted side:" << pushInside
549 <<
" obtained side:" << bb.contains(perturbedPt)
550 <<
" of bb:" << bb <<
nl;
566 const treeBoundBox& bb,
569 const bool pushInside
576 const vector perturbVec = perturbTol_*bb.span();
578 point perturbedPt(pt);
592 if (faceID & treeBoundBox::LEFTBIT)
597 ? (bb.min()[dir] + (perturbVec[dir] + ROOTVSMALL))
598 : (bb.
min()[dir] - (perturbVec[dir] + ROOTVSMALL))
601 else if (faceID & treeBoundBox::RIGHTBIT)
606 ? (bb.max()[dir] - (perturbVec[dir] + ROOTVSMALL))
607 : (bb.
max()[dir] + (perturbVec[dir] + ROOTVSMALL))
615 if (faceID & treeBoundBox::BOTTOMBIT)
620 ? (bb.min()[dir] + (perturbVec[dir] + ROOTVSMALL))
621 : (bb.
min()[dir] - (perturbVec[dir] + ROOTVSMALL))
624 else if (faceID & treeBoundBox::TOPBIT)
629 ? (bb.max()[dir] - (perturbVec[dir] + ROOTVSMALL))
630 : (bb.
max()[dir] + (perturbVec[dir] + ROOTVSMALL))
638 if (faceID & treeBoundBox::BACKBIT)
643 ? (bb.min()[dir] + (perturbVec[dir] + ROOTVSMALL))
644 : (bb.
min()[dir] - (perturbVec[dir] + ROOTVSMALL))
647 else if (faceID & treeBoundBox::FRONTBIT)
652 ? (bb.max()[dir] - (perturbVec[dir] + ROOTVSMALL))
653 : (bb.
max()[dir] + (perturbVec[dir] + ROOTVSMALL))
660 if (pushInside != bb.contains(perturbedPt))
663 <<
"pushed point:" << pt <<
" on face:" << faceString(faceID)
664 <<
" to:" << perturbedPt
665 <<
" wanted side:" << pushInside
666 <<
" obtained side:" << bb.contains(perturbedPt)
667 <<
" of bb:" << bb <<
nl;
683 const treeBoundBox& bb,
690 if (bb.posBits(pt) != 0)
693 <<
" bb:" << bb <<
endl 694 <<
"does not contain point " << pt <<
nl;
711 FixedList<direction, 3> faceIndices;
713 if (ptFaceID & treeBoundBox::LEFTBIT)
715 faceIndices[nFaces++] = treeBoundBox::LEFT;
717 else if (ptFaceID & treeBoundBox::RIGHTBIT)
719 faceIndices[nFaces++] = treeBoundBox::RIGHT;
722 if (ptFaceID & treeBoundBox::BOTTOMBIT)
724 faceIndices[nFaces++] = treeBoundBox::BOTTOM;
726 else if (ptFaceID & treeBoundBox::TOPBIT)
728 faceIndices[nFaces++] = treeBoundBox::TOP;
731 if (ptFaceID & treeBoundBox::BACKBIT)
733 faceIndices[nFaces++] = treeBoundBox::BACK;
735 else if (ptFaceID & treeBoundBox::FRONTBIT)
737 faceIndices[nFaces++] = treeBoundBox::FRONT;
750 else if (nFaces == 1)
753 keepFaceID = faceIndices[0];
760 keepFaceID = faceIndices[0];
761 scalar maxInproduct =
mag(treeBoundBox::faceNormals[keepFaceID] & dir);
766 scalar
s =
mag(treeBoundBox::faceNormals[face] & dir);
767 if (
s > maxInproduct)
783 if (keepFaceID == treeBoundBox::LEFT)
786 faceID = treeBoundBox::LEFTBIT;
788 else if (keepFaceID == treeBoundBox::RIGHT)
791 faceID = treeBoundBox::RIGHTBIT;
793 else if (keepFaceID == treeBoundBox::BOTTOM)
796 faceID = treeBoundBox::BOTTOMBIT;
798 else if (keepFaceID == treeBoundBox::TOP)
801 faceID = treeBoundBox::TOPBIT;
803 else if (keepFaceID == treeBoundBox::BACK)
806 faceID = treeBoundBox::BACKBIT;
808 else if (keepFaceID == treeBoundBox::FRONT)
811 faceID = treeBoundBox::FRONTBIT;
820 <<
"Pushed point from " << pt
821 <<
" on face:" << ptFaceID <<
" of bb:" << bb <<
nl 823 <<
" on face:" << faceID
824 <<
" which is not consistent with geometric face " 835 <<
" bb:" << bb <<
nl 836 <<
"does not contain perturbed point " 861 parentNodeI = nodes_[nodeI].parent_;
863 if (parentNodeI == -1)
869 const node& parentNode = nodes_[parentNodeI];
874 for (
direction i = 0; i < node::nChildren; ++i)
876 labelBits index = parentNode.subNodes_[i];
878 if (isNode(index) && getNode(index) == nodeI)
885 if (parentOctant == 255)
888 <<
"Problem: no parent found for octant:" << octant
912 label oldNodeI = nodeI;
920 const direction X = treeBoundBox::RIGHTHALF;
922 const direction Z = treeBoundBox::FRONTHALF;
927 if ((faceID & treeBoundBox::LEFTBIT) != 0)
933 else if ((faceID & treeBoundBox::RIGHTBIT) != 0)
938 if ((faceID & treeBoundBox::BOTTOMBIT) != 0)
944 else if ((faceID & treeBoundBox::TOPBIT) != 0)
950 if ((faceID & treeBoundBox::BACKBIT) != 0)
955 else if ((faceID & treeBoundBox::FRONTBIT) != 0)
987 while (wantedValue != (octant & octantMask))
1013 if (wantedValue &
Y)
1030 if (wantedValue & Z)
1050 walkToParent(nodeI, octant, parentNodeI, parentOctant);
1052 if (parentNodeI == -1)
1067 nodeI = parentNodeI;
1068 octant = parentOctant;
1074 octant ^= octantMask;
1082 const treeBoundBox subBb(subBbox(nodeI, octant));
1088 <<
" ended up in node:" << nodeI
1089 <<
" octant:" << octant
1090 <<
" with bb:" << subBb <<
nl;
1102 labelBits index = nodes_[nodeI].subNodes_[octant];
1106 labelBits node = findNode(getNode(index),
facePoint);
1108 nodeI = getNode(node);
1115 const treeBoundBox subBb(subBbox(nodeI, octant));
1117 if (nodeI == oldNodeI && octant == oldOctant)
1120 <<
"Did not go to neighbour when searching for " <<
facePoint 1122 <<
" starting from face:" << faceString(faceID)
1123 <<
" node:" << nodeI
1124 <<
" octant:" << octant
1125 <<
" bb:" << subBb <<
nl;
1137 <<
" ended up in node:" << nodeI
1138 <<
" octant:" << octant
1139 <<
" bb:" << subBb <<
nl;
1153 template<
class Type>
1165 if (faceID & treeBoundBox::LEFTBIT)
1167 if (!desc.empty()) desc +=
"+";
1170 if (faceID & treeBoundBox::RIGHTBIT)
1172 if (!desc.empty()) desc +=
"+";
1175 if (faceID & treeBoundBox::BOTTOMBIT)
1177 if (!desc.empty()) desc +=
"+";
1180 if (faceID & treeBoundBox::TOPBIT)
1182 if (!desc.empty()) desc +=
"+";
1185 if (faceID & treeBoundBox::BACKBIT)
1187 if (!desc.empty()) desc +=
"+";
1190 if (faceID & treeBoundBox::FRONTBIT)
1192 if (!desc.empty()) desc +=
"+";
1199 template<
class Type>
1203 const point& treeStart,
1217 const treeBoundBox octantBb(subBbox(nodeI, octant));
1219 if (octantBb.posBits(start) != 0)
1222 <<
"Node:" << nodeI <<
" octant:" << octant
1223 <<
" bb:" << octantBb <<
nl 1224 <<
"does not contain point " << start <<
nl;
1234 const node& nod = nodes_[nodeI];
1236 labelBits index = nod.subNodes_[octant];
1238 if (isContent(index))
1240 const labelList& indices = contents_[getContent(index)];
1250 label shapeI = indices[elemI];
1253 bool hit = shapes_.intersects(shapeI, start,
end, pt);
1261 hitInfo.hitPoint(pt);
1262 hitInfo.setIndex(shapeI);
1271 const treeBoundBox octantBb(subBbox(nodeI, octant));
1277 label shapeI = indices[elemI];
1280 bool hit = shapes_.intersects
1293 if (hit && octantBb.contains(pt))
1298 hitInfo.hitPoint(pt);
1299 hitInfo.setIndex(shapeI);
1317 const treeBoundBox octantBb(subBbox(nodeI, octant));
1320 bool intersected = octantBb.intersects
1336 hitInfo.setPoint(pt);
1345 point perturbedEnd(pushPoint(octantBb,
end,
false));
1364 template<
class Type>
1368 const point& treeStart,
1369 const point& treeEnd,
1370 const label startNodeI,
1375 const vector treeVec(treeEnd - treeStart);
1378 label nodeI = startNodeI;
1383 Pout<<
"findLine : treeStart:" << treeStart
1384 <<
" treeEnd:" << treeEnd <<
endl 1386 <<
" octant:" << octant
1387 <<
" bb:" << subBbox(nodeI, octant) <<
endl;
1395 for (; i < 100000; i++)
1400 const treeBoundBox octantBb(subBbox(nodeI, octant));
1416 <<
" at current:" << hitInfo.point()
1417 <<
" (perturbed:" << startPoint <<
")" <<
endl 1418 <<
" node:" << nodeI
1419 <<
" octant:" << octant
1420 <<
" bb:" << subBbox(nodeI, octant) <<
endl;
1451 if (hitFaceID == 0 || hitInfo.point() == treeEnd)
1458 point perturbedPoint
1471 Pout<<
" iter:" << i
1472 <<
" hit face:" << faceString(hitFaceID)
1473 <<
" at:" << hitInfo.point() <<
nl 1474 <<
" node:" << nodeI
1475 <<
" octant:" << octant
1476 <<
" bb:" << subBbox(nodeI, octant) <<
nl 1477 <<
" walking to neighbour containing:" << perturbedPoint
1486 bool ok = walkToNeighbour
1503 const treeBoundBox octantBb(subBbox(nodeI, octant));
1504 Pout<<
" walked for point:" << hitInfo.point() <<
endl 1505 <<
" to neighbour node:" << nodeI
1506 <<
" octant:" << octant
1507 <<
" face:" << faceString(octantBb.faceBits(hitInfo.point()))
1508 <<
" of octantBb:" << octantBb <<
endl 1532 <<
"Got stuck in loop raytracing from:" << treeStart
1533 <<
" to:" << treeEnd <<
endl 1534 <<
"inside top box:" << subBbox(startNodeI, startOctant)
1540 <<
"Got stuck in loop raytracing from:" << treeStart
1541 <<
" to:" << treeEnd <<
endl 1542 <<
"inside top box:" << subBbox(startNodeI, startOctant)
1551 template<
class Type>
1563 const treeBoundBox& treeBb = nodes_[0].bb_;
1568 direction startBit = treeBb.posBits(start);
1571 if ((startBit & endBit) != 0)
1580 point trackStart(start);
1586 if (!treeBb.intersects(start,
end, trackStart))
1595 if (!treeBb.intersects(
end, trackStart, trackEnd))
1603 labelBits index = findNode(0, trackStart);
1605 label parentNodeI = getNode(index);
1622 template<
class Type>
1626 const treeBoundBox& searchBox,
1630 const node& nod = nodes_[nodeI];
1631 const treeBoundBox& nodeBb = nod.bb_;
1633 bool foundAny =
false;
1635 for (
direction octant = 0; octant < node::nChildren; ++octant)
1637 labelBits index = nod.subNodes_[octant];
1641 const treeBoundBox& subBb = nodes_[getNode(index)].bb_;
1643 if (subBb.overlaps(searchBox))
1645 if (findBox(getNode(index), searchBox, elements))
1648 if (!elements)
return true;
1654 else if (isContent(index))
1656 if (nodeBb.subOverlaps(octant, searchBox))
1658 const labelList& indices = contents_[getContent(index)];
1660 for (
const label index : indices)
1662 if (shapes_.overlaps(index, searchBox))
1665 if (!elements)
return true;
1668 elements->insert(index);
1679 template<
class Type>
1683 const point& centre,
1684 const scalar radiusSqr,
1688 const node& nod = nodes_[nodeI];
1689 const treeBoundBox& nodeBb = nod.bb_;
1691 bool foundAny =
false;
1693 for (
direction octant = 0; octant < node::nChildren; ++octant)
1695 labelBits index = nod.subNodes_[octant];
1699 const treeBoundBox& subBb = nodes_[getNode(index)].bb_;
1701 if (subBb.overlaps(centre, radiusSqr))
1703 if (findSphere(getNode(index), centre, radiusSqr, elements))
1706 if (!elements)
return true;
1712 else if (isContent(index))
1714 if (nodeBb.subOverlaps(octant, centre, radiusSqr))
1716 const labelList& indices = contents_[getContent(index)];
1718 for (
const label index : indices)
1720 if (shapes_.overlaps(index, centre, radiusSqr))
1723 if (!elements)
return true;
1726 elements->insert(index);
1737 template<
class Type>
1738 template<
class CompareOp>
1741 const scalar nearDist,
1743 const dynamicIndexedOctree<Type>& tree1,
1744 const labelBits index1,
1745 const treeBoundBox& bb1,
1746 const dynamicIndexedOctree<Type>& tree2,
1747 const labelBits index2,
1748 const treeBoundBox& bb2,
1752 const vector nearSpan(nearDist, nearDist, nearDist);
1754 if (tree1.isNode(index1))
1756 const node& nod1 = tree1.nodes()[tree1.getNode(index1)];
1757 const treeBoundBox searchBox
1763 if (tree2.isNode(index2))
1765 if (bb2.overlaps(searchBox))
1767 const node& nod2 = tree2.nodes()[tree2.getNode(index2)];
1769 for (
direction i2 = 0; i2 < node::nChildren; ++i2)
1771 labelBits subIndex2 = nod2.subNodes_[i2];
1772 const treeBoundBox subBb2
1774 tree2.isNode(subIndex2)
1775 ? tree2.nodes()[tree2.getNode(subIndex2)].bb_
1794 else if (tree2.isContent(index2))
1797 for (
direction i1 = 0; i1 < node::nChildren; ++i1)
1799 labelBits subIndex1 = nod1.subNodes_[i1];
1800 const treeBoundBox subBb1
1802 tree1.isNode(subIndex1)
1803 ? tree1.nodes()[tree1.getNode(subIndex1)].bb_
1822 else if (tree1.isContent(index1))
1824 const treeBoundBox searchBox
1830 if (tree2.isNode(index2))
1833 tree2.nodes()[tree2.getNode(index2)];
1835 if (bb2.overlaps(searchBox))
1837 for (
direction i2 = 0; i2 < node::nChildren; ++i2)
1839 labelBits subIndex2 = nod2.subNodes_[i2];
1840 const treeBoundBox subBb2
1842 tree2.isNode(subIndex2)
1843 ? tree2.nodes()[tree2.getNode(subIndex2)].bb_
1862 else if (tree2.isContent(index2))
1867 tree1.contents()[tree1.getContent(index1)];
1869 tree2.contents()[tree2.getContent(index2)];
1873 label shape1 = indices1[i];
1877 label shape2 = indices2[j];
1879 if ((&tree1 != &tree2) || (shape1 != shape2))
1911 template<
class Type>
1914 const labelBits index
1922 label nodeI = getNode(index);
1924 const node& nod = nodes_[nodeI];
1926 for (
direction octant = 0; octant < node::nChildren; ++octant)
1928 nElems += countElements(nod.subNodes_[octant]);
1931 else if (isContent(index))
1933 nElems += contents_[getContent(index)].size();
1944 template<
class Type>
1950 const bool leavesOnly,
1951 const bool writeLinesOnly
1954 const node& nod = nodes_[nodeI];
1955 const treeBoundBox& bb = nod.bb_;
1957 for (
direction octant = 0; octant < node::nChildren; ++octant)
1959 const treeBoundBox subBb(bb.subBbox(octant));
1961 labelBits index = nod.subNodes_[octant];
1965 label subNodeI = getNode(index);
1967 writeOBJ(subNodeI,
os, vertIndex, leavesOnly, writeLinesOnly);
1969 else if (isContent(index))
1973 else if (isEmpty(index) && !leavesOnly)
1981 template<
class Type>
1988 labelBits index = nodes_[nodeI].subNodes_[octant];
1994 subBb = nodes_[getNode(index)].bb_;
1996 else if (isContent(index) || isEmpty(index))
1998 subBb = nodes_[nodeI].bb_.
subBbox(octant);
2006 Pout<<
"dumpContentNode : writing node:" << nodeI <<
" octant:" << octant
2007 <<
" to " <<
os.name() <<
endl;
2009 bool writeLinesOnly(
false);
2015 template<
class Type>
2018 if (!nodes_.empty())
2029 template<
class Type>
2034 const label maxLevels,
2035 const scalar maxLeafRatio,
2036 const scalar maxDuplicity
2041 maxLevels_(maxLevels),
2043 maxLeafRatio_(maxLeafRatio),
2044 minSize_(label(maxLeafRatio)),
2045 maxDuplicity_(maxDuplicity),
2046 nodes_(label(shapes.size() / maxLeafRatio_)),
2047 contents_(label(shapes.size() / maxLeafRatio_)),
2050 if (shapes_.size() == 0)
2055 insert(0, shapes_.size());
2066 template<
class Type>
2069 const point& sample,
2070 const scalar startDistSqr
2073 scalar nearestDistSqr = startDistSqr;
2074 label nearestShapeI = -1;
2090 return pointIndexHit(nearestShapeI != -1, nearestPoint, nearestShapeI);
2094 template<
class Type>
2102 label nearestShapeI = -1;
2119 return pointIndexHit(nearestShapeI != -1, nearestPoint, nearestShapeI);
2123 template<
class Type>
2130 return findLine(
false, start,
end);
2134 template<
class Type>
2141 return findLine(
true, start,
end);
2145 template<
class Type>
2152 return !nodes_.empty() && findBox(0, searchBox,
nullptr);
2156 template<
class Type>
2165 if (!nodes_.empty())
2168 label estimatedCount(
max(128, (shapes_.size() / 100)));
2169 elements.
reserve(estimatedCount);
2172 findBox(0, searchBox, &elements);
2175 return elements.
size();
2179 template<
class Type>
2192 findBox(searchBox, elements);
2195 return elements.
toc();
2199 template<
class Type>
2202 const point& centre,
2203 const scalar radiusSqr
2207 return !nodes_.empty() && findSphere(0, centre, radiusSqr,
nullptr);
2211 template<
class Type>
2214 const point& centre,
2215 const scalar radiusSqr,
2221 if (!nodes_.empty())
2224 label estimatedCount(
max(128, (shapes_.size() / 100)));
2225 elements.
reserve(estimatedCount);
2228 findSphere(0, centre, radiusSqr, &elements);
2231 return elements.
size();
2235 template<
class Type>
2238 const point& centre,
2239 const scalar radiusSqr
2249 findSphere(centre, radiusSqr, elements);
2252 return elements.
toc();
2256 template<
class Type>
2266 return nodePlusOctant(nodeI, 0);
2269 const node& nod = nodes_[nodeI];
2273 if (!nod.bb_.contains(sample))
2276 <<
"Cannot find " << sample <<
" in node " << nodeI
2281 direction octant = nod.bb_.subOctant(sample);
2283 labelBits index = nod.subNodes_[octant];
2288 return findNode(getNode(index), sample);
2290 else if (isContent(index))
2293 return nodePlusOctant(nodeI, octant);
2298 return nodePlusOctant(nodeI, octant);
2303 template<
class Type>
2309 labelBits index = findNode(0, sample);
2311 const node& nod = nodes_[getNode(index)];
2313 labelBits contentIndex = nod.subNodes_[
getOctant(index)];
2316 if (isContent(contentIndex))
2318 const labelList& indices = contents_[getContent(contentIndex)];
2322 label shapeI = indices[elemI];
2324 if (shapes_.contains(shapeI, sample))
2335 template<
class Type>
2343 const node& nod = nodes_[getNode(index)];
2348 if (isContent(contentIndex))
2350 return contents_[getContent(contentIndex)];
2353 return labelList::null();
2357 template<
class Type>
2365 return volumeType::UNKNOWN;
2368 if (nodeTypes_.size() != 8*nodes_.size())
2372 nodeTypes_.setSize(8*nodes_.size());
2373 nodeTypes_ = volumeType::UNKNOWN;
2388 if (
type == volumeType::UNKNOWN)
2396 else if (
type == volumeType::INSIDE)
2400 else if (
type == volumeType::OUTSIDE)
2410 Pout<<
"dynamicIndexedOctree::getVolumeType : " 2412 <<
" nodes_:" << nodes_.size()
2413 <<
" nodeTypes_:" << nodeTypes_.size()
2414 <<
" nUNKNOWN:" << nUNKNOWN
2415 <<
" nMIXED:" << nMIXED
2416 <<
" nINSIDE:" << nINSIDE
2417 <<
" nOUTSIDE:" << nOUTSIDE
2422 return getVolumeType(0, sample);
2426 template<
class Type>
2427 template<
class CompareOp>
2430 const scalar nearDist,
2440 nodePlusOctant(0, 0),
2443 nodePlusOctant(0, 0),
2450 template<
class Type>
2453 if (startIndex == endIndex)
2460 contents_.emplace_back(1).push_back(0);
2463 node topNode =
divide(bb_, 0, -1, 0);
2465 nodes_.push_back(topNode);
2472 for (label pI = startIndex; pI < endIndex; ++pI)
2476 if (!insertIndex(0, pI, nLevels))
2481 nLevelsMax_ =
max(nLevels, nLevelsMax_);
2488 template<
class Type>
2491 const label nodIndex,
2496 bool shapeInserted =
false;
2498 for (
direction octant = 0; octant < node::nChildren; ++octant)
2500 const labelBits subNodeLabel = nodes_[nodIndex].subNodes_[octant];
2502 if (isNode(subNodeLabel))
2504 const treeBoundBox& subBb = nodes_[getNode(subNodeLabel)].bb_;
2506 if (shapes().overlaps(index, subBb))
2510 if (insertIndex(getNode(subNodeLabel), index, nLevels))
2512 shapeInserted =
true;
2516 else if (isContent(subNodeLabel))
2520 if (shapes().overlaps(index, subBb))
2522 const label contentI = getContent(subNodeLabel);
2524 contents_[contentI].push_back(index);
2526 recursiveSubDivision
2535 shapeInserted =
true;
2542 if (shapes().overlaps(index, subBb))
2544 const label sz = contents_.size();
2546 contents_.emplace_back(1).push_back(index);
2548 nodes_[nodIndex].subNodes_[octant]
2549 = contentPlusOctant(sz, octant);
2552 shapeInserted =
true;
2556 return shapeInserted;
2560 template<
class Type>
2568 removeIndex(0, index);
2574 template<
class Type>
2577 const label nodIndex,
2581 label totalContents = 0;
2583 for (
direction octant = 0; octant < node::nChildren; ++octant)
2585 const labelBits subNodeLabel = nodes_[nodIndex].subNodes_[octant];
2587 if (isNode(subNodeLabel))
2589 const treeBoundBox& subBb = nodes_[getNode(subNodeLabel)].bb_;
2591 if (shapes().overlaps(index, subBb))
2594 label childContentsSize
2595 = removeIndex(getNode(subNodeLabel), index);
2597 totalContents += childContentsSize;
2599 if (childContentsSize == 0)
2601 nodes_[nodIndex].subNodes_[octant]
2602 = emptyPlusOctant(octant);
2611 else if (isContent(subNodeLabel))
2615 const label contentI = getContent(subNodeLabel);
2617 if (shapes().overlaps(index, subBb))
2625 const label oldIndex = contentList[pI];
2627 if (oldIndex != index)
2633 newContent.shrink();
2635 if (newContent.empty())
2638 nodes_[nodIndex].subNodes_[octant]
2639 = emptyPlusOctant(octant);
2645 totalContents += contents_[contentI].size();
2653 return totalContents;
2657 template<
class Type>
2661 const bool printContents,
2665 const node& nod = nodes_[nodeI];
2668 os <<
"nodeI:" << nodeI <<
" bb:" << bb <<
nl 2669 <<
"parent:" << nod.parent_ <<
nl 2670 <<
"n:" << countElements(nodePlusOctant(nodeI, 0)) <<
nl;
2672 for (
direction octant = 0; octant < node::nChildren; ++octant)
2676 labelBits index = nod.subNodes_[octant];
2681 label subNodeI = getNode(index);
2683 os <<
"octant:" << octant
2684 <<
" node: n:" << countElements(index)
2685 <<
" bb:" << subBb <<
endl;
2687 string oldPrefix =
os.prefix();
2688 os.prefix() =
" " + oldPrefix;
2690 print(
os, printContents, subNodeI);
2692 os.prefix() = oldPrefix;
2694 else if (isContent(index))
2696 const labelList& indices = contents_[getContent(index)];
2703 os <<
"octant:" << octant
2704 <<
" content: n:" << indices.
size()
2712 os <<
' ' << indices[j];
2724 os <<
"octant:" << octant <<
" empty:" << subBb <<
endl;
2730 template<
class Type>
2734 for (
const auto& subContents : contents_)
2736 nEntries += subContents.size();
2739 Pout<<
"indexedOctree::indexedOctree" 2740 <<
" : finished construction of tree of:" << shapes().typeName
2742 <<
" bounding box: " << this->bb() <<
nl 2743 <<
" shapes: " << shapes().size() <<
nl 2744 <<
" treeNodes: " << nodes_.size() <<
nl 2745 <<
" nEntries: " << nEntries <<
nl 2746 <<
" levels/maxLevels: " << nLevelsMax_ <<
"/" << maxLevels_ <<
nl 2747 <<
" minSize: " << minSize_ <<
nl 2748 <<
" per treeLeaf: " 2749 << scalar(nEntries)/contents_.size() <<
nl 2750 <<
" per shape (duplicity):" 2751 << scalar(nEntries)/shapes().size() <<
nl 2756 template<
class Type>
2765 template<
class Type>
2767 Foam::operator<<(Ostream& os, const dynamicIndexedOctree<Type>& t)
2769 os << t.bb() << token::SPACE << t.nodes() <<
endl;
2771 for (
const auto& subContents : t.contents())
2773 os << subContents <<
endl;
const labelList & findIndices(const point &) const
Find the shape indices that occupy the result of findNode.
void size(const label n)
Older name for setAddressableSize.
void divide(DimensionedField< Type, GeoMesh > &result, const DimensionedField< Type, GeoMesh > &f1, const DimensionedField< scalar, GeoMesh > &f2)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Mixed uniform/non-uniform (eg, after reduction)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Ostream & print(Ostream &os, UIntType value, char off='0', char on='1')
Print 0/1 bits in the (unsigned) integral type.
srcOptions insert("case", fileName(rootDirSource/caseDirSource))
void print(prefixOSstream &, const bool printContents, const label) const
Print tree. Either print all indices (printContent = true) or.
constexpr char nl
The newline '\n' character (0x0a)
Ostream & endl(Ostream &os)
Add newline and flush stream.
PointIndexHit< point > pointIndexHit
A PointIndexHit with a 3D point.
bool insert(label startIndex, label endIndex)
Insert a new object into the tree.
This class describes the interaction of an object (often a face) and a point. It carries the info of ...
::Foam::direction nComponents(const expressions::valueTypeCode) noexcept
The number of components associated with given valueTypeCode.
bool remove(const label index)
Remove an object from the tree.
An enumeration wrapper for classification of a location as being inside/outside of a volume...
UList< label > labelUList
A UList of labels.
Various functions to operate on Lists.
#define forAll(list, i)
Loop across all elements in list.
label removeIndex(const label nodIndex, const label index)
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
label size() const noexcept
The number of elements in table.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
line< point, const point & > linePointRef
A line using referred points.
Version of OSstream that prints a prefix on each line.
A class for handling words, derived from Foam::string.
void clear()
Remove all entries from table.
void writeTreeInfo() const
Non-pointer based hierarchical recursive searching. Storage is dynamic, so elements can be deleted...
treeBoundBox subBbox(const direction octant) const
Sub-box of given octant. Midpoint calculated.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
bool ln(const fileName &src, const fileName &dst)
Create a softlink. dst should not exist. Returns true if successful.
An Ostream is an abstract base class for all output systems (streams, files, token lists...
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
bool insertIndex(const label nodIndex, const label index, label &nLevels)
int debug
Static debugging option.
label facePoint(const int facei, const block &block, const label i, const label j)
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
OBJstream os(runTime.globalPath()/outputName)
void push_back(const T &val)
Copy append an element to the end of this list.
dynamicIndexedOctree(const Type &shapes, const treeBoundBox &bb, const label maxLevels, const scalar maxLeafRatio, const scalar maxDuplicity)
Construct from shapes.
static volumeType getSide(const vector &outsideNormal, const vector &vec)
Helper function to return the side. Returns outside if.
bool overlaps(const treeBoundBox &bb) const
True if any shapes overlap the bounding box.
vector point
Point is a vector.
#define WarningInFunction
Report a warning using Foam::Warning.
pointIndexHit findLineAny(const point &start, const point &end) const
Find any intersection of line between start and end.
PtrList< volScalarField > & Y
void reserve(label numEntries)
Reserve space for at least the specified number of elements (not the number of buckets) and regenerat...
void transfer(List< T > &list)
Transfer contents of the argument List into this.
Standard boundBox with extra functionality for use in octree.
A 29bits (or 61bits) integer label with 3bits direction (eg, octant) packed into single label...
labelBits findNode(const label nodeI, const point &) const
Find deepest node (as parent+octant) containing point. Starts.
List< Key > toc() const
The table of contents (the keys) in unsorted order.
List< label > labelList
A List of labels.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
static unsigned getOctant(const point &p)
label findInside(const point &) const
Find shape containing point. Only implemented for certain.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
bool write(Ostream &os) const
Tree node. Has up pointer and down pointers.
static constexpr const zero Zero
Global zero (0)
const treeBoundBox & bb() const
Top bounding box.