41 const treeBoundBox& bb,
42 FixedList<labelList, 8>& dividedIndices
46 DynamicList<label> contains(indices.size());
48 for (
direction octant = 0; octant < 8; ++octant)
50 const treeBoundBox subBbs(bb.subBbox(octant));
54 for (
const label index : indices)
56 if (shapes_.overlaps(index, subBbs))
58 contains.push_back(index);
63 dividedIndices[octant] = contains;
72 const treeBoundBox& bb,
73 DynamicList<labelList>& contents,
79 bb.min().x() >= bb.max().x()
80 || bb.min().y() >= bb.max().y()
81 || bb.min().z() >= bb.max().z()
85 <<
"Badly formed bounding box:" << bb
93 const labelList& indices = contents[contentIndex];
95 FixedList<labelList, 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);
140 DynamicList<indexedOctreeBase::node>& nodes,
141 DynamicList<labelList>& contents
144 label currentSize = nodes.size();
149 for (label nodeI = 0; nodeI < currentSize; nodeI++)
151 for (
direction octant = 0; octant < node::nChildren; ++octant)
153 labelBits index = nodes[nodeI].subNodes_[octant];
159 else if (isContent(index))
161 label contentI = getContent(index);
163 if (contents[contentI].size() > minSize)
168 const node& nod = nodes[nodeI];
169 const treeBoundBox bb(nod.bb_.subBbox(octant));
171 node subNode(
divide(bb, contents, contentI));
172 subNode.parent_ = nodeI;
173 label sz = nodes.size();
174 nodes.push_back(subNode);
175 nodes[nodeI].subNodes_[octant] = nodePlusOctant(sz, octant);
186 DynamicList<node>& nodes,
187 DynamicList<labelList>& contents,
188 const label compactLevel,
192 List<labelList>& compactedContents,
196 const node& nod = nodes[nodeI];
200 if (level < compactLevel)
202 for (
direction octant = 0; octant < node::nChildren; ++octant)
204 labelBits index = nod.subNodes_[octant];
208 nNodes += compactContents
221 else if (level == compactLevel)
224 for (
direction octant = 0; octant < node::nChildren; ++octant)
226 labelBits index = nod.subNodes_[octant];
228 if (isContent(index))
230 label contentI = getContent(index);
232 compactedContents[compactI].transfer(contents[contentI]);
235 nodes[nodeI].subNodes_[octant] =
236 contentPlusOctant(compactI, octant);
240 else if (isNode(index))
260 const node& nod = nodes_[nodeI];
262 volumeType myType = volumeType::UNKNOWN;
264 for (
direction octant = 0; octant < node::nChildren; ++octant)
268 labelBits index = nod.subNodes_[octant];
273 subType = calcVolumeType(getNode(index));
275 else if (isContent(index))
279 subType = volumeType::MIXED;
285 const treeBoundBox subBb = nod.bb_.subBbox(octant);
287 subType = shapes_.getVolumeType(*
this, subBb.centre());
291 nodeTypes_.set(labelBits::pack(nodeI, octant), subType);
295 if (myType == volumeType::UNKNOWN)
299 else if (subType != myType)
301 myType = volumeType::MIXED;
315 const node& nod = nodes_[nodeI];
317 direction octant = nod.bb_.subOctant(sample);
321 nodeTypes_.get(labelBits::pack(nodeI, octant))
324 if (octantType == volumeType::INSIDE)
328 else if (octantType == volumeType::OUTSIDE)
332 else if (octantType == volumeType::UNKNOWN)
337 else if (octantType == volumeType::MIXED)
339 labelBits index = nod.subNodes_[octant];
344 volumeType subType = getVolumeType(getNode(index), sample);
348 else if (isContent(index))
351 return volumeType(shapes_.getVolumeType(*
this, sample));
357 <<
"Sample:" << sample <<
" node:" << nodeI
358 <<
" with bb:" << nodes_[nodeI].bb_ <<
nl 359 <<
"Empty subnode has invalid volume type MIXED." 362 return volumeType::UNKNOWN;
366 <<
"Sample:" << sample <<
" at node:" << nodeI
367 <<
" octant:" << octant
368 <<
" with bb:" << nod.bb_.subBbox(octant) <<
nl 369 <<
"Node has invalid volume type " << octantType
372 return volumeType::UNKNOWN;
379 const vector& outsideNormal,
383 if ((outsideNormal&vec) >= 0)
385 return volumeType::OUTSIDE;
389 return volumeType::INSIDE;
395 template<
class FindNearestOp>
401 scalar& nearestDistSqr,
402 label& nearestShapeI,
405 const FindNearestOp& fnOp
408 const node& nod = nodes_[nodeI];
413 for (
const direction octant : nod.bb_.searchOrder(sample))
415 labelBits index = nod.subNodes_[octant];
419 label subNodeI = getNode(index);
421 const treeBoundBox& subBb = nodes_[subNodeI].bb_;
423 if (subBb.overlaps(sample, nearestDistSqr))
438 else if (isContent(index))
440 if (nod.bb_.subOverlaps(octant, sample, nearestDistSqr))
444 contents_[getContent(index)],
458 template<
class FindNearestOp>
464 treeBoundBox& tightest,
465 label& nearestShapeI,
469 const FindNearestOp& fnOp
472 const node& nod = nodes_[nodeI];
473 const treeBoundBox& nodeBb = nod.bb_;
478 for (
const direction octant : nod.bb_.searchOrder(
ln.centre()))
480 labelBits index = nod.subNodes_[octant];
484 const treeBoundBox& subBb = nodes_[getNode(index)].bb_;
486 if (subBb.overlaps(tightest))
502 else if (isContent(index))
504 if (nodeBb.subOverlaps(octant, tightest))
508 contents_[getContent(index)],
525 const label parentNodeI,
530 const node& nod = nodes_[parentNodeI];
531 labelBits index = nod.subNodes_[octant];
536 return nodes_[getNode(index)].bb_;
540 return nod.bb_.subBbox(octant);
547 const treeBoundBox& bb,
549 const bool pushInside
553 const vector perturbVec = perturbTol_*bb.span();
555 point perturbedPt(pt);
562 for (
direction dir = 0; dir < vector::nComponents; dir++)
564 if (
mag(pt[dir]-bb.min()[dir]) <
mag(perturbVec[dir]))
567 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
568 perturbedPt[dir] = bb.min()[dir] + perturbDist;
570 else if (
mag(pt[dir]-bb.max()[dir]) <
mag(perturbVec[dir]))
573 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
574 perturbedPt[dir] = bb.max()[dir] - perturbDist;
580 for (
direction dir = 0; dir < vector::nComponents; dir++)
582 if (
mag(pt[dir]-bb.min()[dir]) <
mag(perturbVec[dir]))
584 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
585 perturbedPt[dir] = bb.min()[dir] - perturbDist;
587 else if (
mag(pt[dir]-bb.max()[dir]) <
mag(perturbVec[dir]))
589 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
590 perturbedPt[dir] = bb.max()[dir] + perturbDist;
597 if (pushInside != bb.contains(perturbedPt))
600 <<
"pushed point:" << pt
601 <<
" to:" << perturbedPt
602 <<
" wanted side:" << pushInside
603 <<
" obtained side:" << bb.contains(perturbedPt)
604 <<
" of bb:" << bb <<
nl;
620 const treeBoundBox& bb,
623 const bool pushInside
627 const vector perturbVec = perturbTol_*bb.span();
629 point perturbedPt(pt);
643 if (faceID & treeBoundBox::LEFTBIT)
648 ? (bb.min()[dir] + (perturbVec[dir] + ROOTVSMALL))
649 : (bb.
min()[dir] - (perturbVec[dir] + ROOTVSMALL))
652 else if (faceID & treeBoundBox::RIGHTBIT)
657 ? (bb.max()[dir] - (perturbVec[dir] + ROOTVSMALL))
658 : (bb.
max()[dir] + (perturbVec[dir] + ROOTVSMALL))
666 if (faceID & treeBoundBox::BOTTOMBIT)
671 ? (bb.min()[dir] + (perturbVec[dir] + ROOTVSMALL))
672 : (bb.
min()[dir] - (perturbVec[dir] + ROOTVSMALL))
675 else if (faceID & treeBoundBox::TOPBIT)
680 ? (bb.max()[dir] - (perturbVec[dir] + ROOTVSMALL))
681 : (bb.
max()[dir] + (perturbVec[dir] + ROOTVSMALL))
689 if (faceID & treeBoundBox::BACKBIT)
694 ? (bb.min()[dir] + (perturbVec[dir] + ROOTVSMALL))
695 : (bb.
min()[dir] - (perturbVec[dir] + ROOTVSMALL))
698 else if (faceID & treeBoundBox::FRONTBIT)
703 ? (bb.max()[dir] - (perturbVec[dir] + ROOTVSMALL))
704 : (bb.
max()[dir] + (perturbVec[dir] + ROOTVSMALL))
711 if (pushInside != bb.contains(perturbedPt))
714 <<
"pushed point:" << pt <<
" on face:" << faceString(faceID)
715 <<
" to:" << perturbedPt
716 <<
" wanted side:" << pushInside
717 <<
" obtained side:" << bb.contains(perturbedPt)
718 <<
" of bb:" << bb <<
nl;
734 const treeBoundBox& bb,
741 if (bb.posBits(pt) != 0)
744 <<
" bb:" << bb <<
endl 745 <<
"does not contain point " << pt <<
nl;
762 FixedList<direction, 3> faceIndices;
764 if (ptFaceID & treeBoundBox::LEFTBIT)
766 faceIndices[nFaces++] = treeBoundBox::LEFT;
768 else if (ptFaceID & treeBoundBox::RIGHTBIT)
770 faceIndices[nFaces++] = treeBoundBox::RIGHT;
773 if (ptFaceID & treeBoundBox::BOTTOMBIT)
775 faceIndices[nFaces++] = treeBoundBox::BOTTOM;
777 else if (ptFaceID & treeBoundBox::TOPBIT)
779 faceIndices[nFaces++] = treeBoundBox::TOP;
782 if (ptFaceID & treeBoundBox::BACKBIT)
784 faceIndices[nFaces++] = treeBoundBox::BACK;
786 else if (ptFaceID & treeBoundBox::FRONTBIT)
788 faceIndices[nFaces++] = treeBoundBox::FRONT;
801 else if (nFaces == 1)
804 keepFaceID = faceIndices[0];
811 keepFaceID = faceIndices[0];
812 scalar maxInproduct =
mag(treeBoundBox::faceNormals[keepFaceID] & dir);
817 scalar
s =
mag(treeBoundBox::faceNormals[face] & dir);
818 if (
s > maxInproduct)
834 if (keepFaceID == treeBoundBox::LEFT)
837 faceID = treeBoundBox::LEFTBIT;
839 else if (keepFaceID == treeBoundBox::RIGHT)
842 faceID = treeBoundBox::RIGHTBIT;
844 else if (keepFaceID == treeBoundBox::BOTTOM)
847 faceID = treeBoundBox::BOTTOMBIT;
849 else if (keepFaceID == treeBoundBox::TOP)
852 faceID = treeBoundBox::TOPBIT;
854 else if (keepFaceID == treeBoundBox::BACK)
857 faceID = treeBoundBox::BACKBIT;
859 else if (keepFaceID == treeBoundBox::FRONT)
862 faceID = treeBoundBox::FRONTBIT;
871 <<
"Pushed point from " << pt
872 <<
" on face:" << ptFaceID <<
" of bb:" << bb <<
nl 874 <<
" on face:" << faceID
875 <<
" which is not consistent with geometric face " 886 <<
" bb:" << bb <<
nl 887 <<
"does not contain perturbed point " 912 parentNodeI = nodes_[nodeI].parent_;
914 if (parentNodeI == -1)
920 const node& parentNode = nodes_[parentNodeI];
925 for (
direction i = 0; i < node::nChildren; ++i)
927 labelBits index = parentNode.subNodes_[i];
929 if (isNode(index) && getNode(index) == nodeI)
936 if (parentOctant == 255)
939 <<
"Problem: no parent found for octant:" << octant
961 label oldNodeI = nodeI;
969 const direction X = treeBoundBox::RIGHTHALF;
971 const direction Z = treeBoundBox::FRONTHALF;
976 if ((faceID & treeBoundBox::LEFTBIT) != 0)
982 else if ((faceID & treeBoundBox::RIGHTBIT) != 0)
987 if ((faceID & treeBoundBox::BOTTOMBIT) != 0)
993 else if ((faceID & treeBoundBox::TOPBIT) != 0)
999 if ((faceID & treeBoundBox::BACKBIT) != 0)
1004 else if ((faceID & treeBoundBox::FRONTBIT) != 0)
1036 while (wantedValue != (octant & octantMask))
1042 if (wantedValue & X)
1062 if (wantedValue &
Y)
1079 if (wantedValue & Z)
1099 walkToParent(nodeI, octant, parentNodeI, parentOctant);
1101 if (parentNodeI == -1)
1116 nodeI = parentNodeI;
1117 octant = parentOctant;
1123 octant ^= octantMask;
1131 const treeBoundBox subBb(subBbox(nodeI, octant));
1137 <<
" ended up in node:" << nodeI
1138 <<
" octant:" << octant
1139 <<
" with bb:" << subBb <<
nl;
1151 labelBits index = nodes_[nodeI].subNodes_[octant];
1155 labelBits node = findNode(getNode(index),
facePoint);
1157 nodeI = getNode(node);
1164 const treeBoundBox subBb(subBbox(nodeI, octant));
1166 if (nodeI == oldNodeI && octant == oldOctant)
1169 <<
"Did not go to neighbour when searching for " <<
facePoint 1171 <<
" starting from face:" << faceString(faceID)
1172 <<
" node:" << nodeI
1173 <<
" octant:" << octant
1174 <<
" bb:" << subBb <<
nl;
1186 <<
" ended up in node:" << nodeI
1187 <<
" octant:" << octant
1188 <<
" bb:" << subBb <<
nl;
1202 template<
class Type>
1214 if (faceID & treeBoundBox::LEFTBIT)
1216 if (!desc.empty()) desc +=
"+";
1219 if (faceID & treeBoundBox::RIGHTBIT)
1221 if (!desc.empty()) desc +=
"+";
1224 if (faceID & treeBoundBox::BOTTOMBIT)
1226 if (!desc.empty()) desc +=
"+";
1229 if (faceID & treeBoundBox::TOPBIT)
1231 if (!desc.empty()) desc +=
"+";
1234 if (faceID & treeBoundBox::BACKBIT)
1236 if (!desc.empty()) desc +=
"+";
1239 if (faceID & treeBoundBox::FRONTBIT)
1241 if (!desc.empty()) desc +=
"+";
1248 template<
class Type>
1249 template<
class FindIntersectOp>
1253 const point& treeStart,
1264 const FindIntersectOp& fiOp
1277 const treeBoundBox octantBb(subBbox(nodeI, octant));
1279 if (octantBb.posBits(start) != 0)
1282 <<
"Node:" << nodeI <<
" octant:" << octant
1283 <<
" bb:" << octantBb <<
nl 1284 <<
"does not contain point " << start <<
nl;
1293 const node& nod = nodes_[nodeI];
1295 labelBits index = nod.subNodes_[octant];
1297 if (isContent(index))
1299 const labelList& indices = contents_[getContent(index)];
1309 label shapeI = indices[elemI];
1312 bool hit = fiOp(shapeI, start,
end, pt);
1320 hitInfo.hitPoint(pt);
1321 hitInfo.setIndex(shapeI);
1330 const treeBoundBox octantBb(subBbox(nodeI, octant));
1336 label shapeI = indices[elemI];
1339 bool hit = fiOp(shapeI, start, nearestPoint, pt);
1346 if (hit && octantBb.contains(pt))
1351 hitInfo.hitPoint(pt);
1352 hitInfo.setIndex(shapeI);
1370 const treeBoundBox octantBb(subBbox(nodeI, octant));
1373 bool intersected = octantBb.intersects
1389 hitInfo.setPoint(pt);
1398 point perturbedEnd(pushPoint(octantBb,
end,
false));
1419 template<
class Type>
1420 template<
class FindIntersectOp>
1424 const point& treeStart,
1425 const point& treeEnd,
1426 const label startNodeI,
1428 const FindIntersectOp& fiOp,
1432 const vector treeVec(treeEnd - treeStart);
1435 label nodeI = startNodeI;
1440 Pout<<
"findLine : treeStart:" << treeStart
1441 <<
" treeEnd:" << treeEnd <<
endl 1443 <<
" octant:" << octant
1444 <<
" bb:" << subBbox(nodeI, octant) <<
endl;
1452 for (; i < 100000; i++)
1457 const treeBoundBox octantBb(subBbox(nodeI, octant));
1473 <<
" at current:" << hitInfo.point()
1474 <<
" (perturbed:" << startPoint <<
")" <<
endl 1475 <<
" node:" << nodeI
1476 <<
" octant:" << octant
1477 <<
" bb:" << subBbox(nodeI, octant) <<
endl;
1510 if (hitFaceID == 0 || hitInfo.point() == treeEnd)
1517 point perturbedPoint
1530 Pout<<
" iter:" << i
1531 <<
" hit face:" << faceString(hitFaceID)
1532 <<
" at:" << hitInfo.point() <<
nl 1533 <<
" node:" << nodeI
1534 <<
" octant:" << octant
1535 <<
" bb:" << subBbox(nodeI, octant) <<
nl 1536 <<
" walking to neighbour containing:" << perturbedPoint
1545 bool ok = walkToNeighbour
1562 const treeBoundBox octantBb(subBbox(nodeI, octant));
1563 Pout<<
" walked for point:" << hitInfo.point() <<
endl 1564 <<
" to neighbour node:" << nodeI
1565 <<
" octant:" << octant
1566 <<
" face:" << faceString(octantBb.faceBits(hitInfo.point()))
1567 <<
" of octantBb:" << octantBb <<
endl 1592 <<
"Got stuck in loop raytracing from:" << treeStart
1593 <<
" to:" << treeEnd <<
endl 1594 <<
"inside top box:" << subBbox(startNodeI, startOctant)
1600 <<
"Got stuck in loop raytracing from:" << treeStart
1601 <<
" to:" << treeEnd <<
endl 1602 <<
"inside top box:" << subBbox(startNodeI, startOctant)
1611 template<
class Type>
1612 template<
class FindIntersectOp>
1618 const FindIntersectOp& fiOp
1633 if ((startBit & endBit) != 0)
1642 point trackStart(start);
1665 labelBits index = findNode(0, trackStart);
1667 label parentNodeI = getNode(index);
1685 template<
class Type>
1689 const treeBoundBox& searchBox,
1693 const node& nod = nodes_[nodeI];
1694 const treeBoundBox& nodeBb = nod.bb_;
1696 bool foundAny =
false;
1698 for (
direction octant = 0; octant < node::nChildren; ++octant)
1700 labelBits index = nod.subNodes_[octant];
1704 const treeBoundBox& subBb = nodes_[getNode(index)].bb_;
1706 if (subBb.overlaps(searchBox))
1708 if (findBox(getNode(index), searchBox, elements))
1711 if (!elements)
return true;
1717 else if (isContent(index))
1719 if (nodeBb.subOverlaps(octant, searchBox))
1721 const labelList& indices = contents_[getContent(index)];
1723 for (
const label index : indices)
1725 if (shapes_.overlaps(index, searchBox))
1728 if (!elements)
return true;
1731 elements->insert(index);
1742 template<
class Type>
1746 const point& centre,
1747 const scalar radiusSqr,
1751 const node& nod = nodes_[nodeI];
1752 const treeBoundBox& nodeBb = nod.bb_;
1754 bool foundAny =
false;
1756 for (
direction octant = 0; octant < node::nChildren; ++octant)
1758 labelBits index = nod.subNodes_[octant];
1762 const treeBoundBox& subBb = nodes_[getNode(index)].bb_;
1764 if (subBb.overlaps(centre, radiusSqr))
1766 if (findSphere(getNode(index), centre, radiusSqr, elements))
1769 if (!elements)
return true;
1775 else if (isContent(index))
1777 if (nodeBb.subOverlaps(octant, centre, radiusSqr))
1779 const labelList& indices = contents_[getContent(index)];
1781 for (
const label index : indices)
1783 if (shapes_.overlaps(index, centre, radiusSqr))
1786 if (!elements)
return true;
1789 elements->insert(index);
1800 template<
class Type>
1801 template<
class CompareOp>
1804 const scalar nearDist,
1806 const indexedOctree<Type>& tree1,
1807 const labelBits index1,
1808 const treeBoundBox& bb1,
1809 const indexedOctree<Type>& tree2,
1810 const labelBits index2,
1811 const treeBoundBox& bb2,
1815 const vector nearSpan(nearDist, nearDist, nearDist);
1817 if (tree1.isNode(index1))
1819 const node& nod1 = tree1.nodes()[tree1.getNode(index1)];
1820 const treeBoundBox searchBox
1826 if (tree2.isNode(index2))
1828 if (bb2.overlaps(searchBox))
1830 const node& nod2 = tree2.nodes()[tree2.getNode(index2)];
1832 for (
direction i2 = 0; i2 < node::nChildren; ++i2)
1834 labelBits subIndex2 = nod2.subNodes_[i2];
1835 const treeBoundBox subBb2
1837 tree2.isNode(subIndex2)
1838 ? tree2.nodes()[tree2.getNode(subIndex2)].bb_
1857 else if (tree2.isContent(index2))
1860 for (
direction i1 = 0; i1 < node::nChildren; ++i1)
1862 labelBits subIndex1 = nod1.subNodes_[i1];
1863 const treeBoundBox subBb1
1865 tree1.isNode(subIndex1)
1866 ? tree1.nodes()[tree1.getNode(subIndex1)].bb_
1885 else if (tree1.isContent(index1))
1887 const treeBoundBox searchBox
1893 if (tree2.isNode(index2))
1896 tree2.nodes()[tree2.getNode(index2)];
1898 if (bb2.overlaps(searchBox))
1900 for (
direction i2 = 0; i2 < node::nChildren; ++i2)
1902 labelBits subIndex2 = nod2.subNodes_[i2];
1903 const treeBoundBox subBb2
1905 tree2.isNode(subIndex2)
1906 ? tree2.nodes()[tree2.getNode(subIndex2)].bb_
1925 else if (tree2.isContent(index2))
1930 tree1.contents()[tree1.getContent(index1)];
1932 tree2.contents()[tree2.getContent(index2)];
1936 label shape1 = indices1[i];
1940 label shape2 = indices2[j];
1942 if ((&tree1 != &tree2) || (shape1 != shape2))
1974 template<
class Type>
1979 const node& nod = nodes_[nodeI];
1981 for (
direction octant = 0; octant < node::nChildren; ++octant)
1983 labelBits index = nod.subNodes_[octant];
1987 total += countLeafs(getNode(index));
1989 else if (isContent(index))
1999 template<
class Type>
2002 const labelBits index
2010 label nodeI = getNode(index);
2012 const node& nod = nodes_[nodeI];
2014 for (
direction octant = 0; octant < node::nChildren; ++octant)
2016 nElems += countElements(nod.subNodes_[octant]);
2019 else if (isContent(index))
2021 nElems += contents_[getContent(index)].size();
2032 template<
class Type>
2038 const bool leavesOnly,
2039 const bool writeLinesOnly
2042 const node& nod = nodes_[nodeI];
2043 const treeBoundBox& bb = nod.bb_;
2045 for (
direction octant = 0; octant < node::nChildren; ++octant)
2047 const treeBoundBox subBb(bb.subBbox(octant));
2049 labelBits index = nod.subNodes_[octant];
2053 label subNodeI = getNode(index);
2055 writeOBJ(subNodeI,
os, vertIndex, leavesOnly, writeLinesOnly);
2057 else if (isContent(index))
2061 else if (isEmpty(index) && !leavesOnly)
2069 template<
class Type>
2076 labelBits index = nodes_[nodeI].subNodes_[octant];
2082 subBb = nodes_[getNode(index)].bb_;
2084 else if (isContent(index) || isEmpty(index))
2086 subBb = nodes_[nodeI].bb_.subBbox(octant);
2094 Pout<<
"dumpContentNode : writing node:" << nodeI <<
" octant:" << octant
2095 <<
" to " <<
os.name() <<
endl;
2097 bool writeLinesOnly(
false);
2103 template<
class Type>
2106 if (!nodes_.empty())
2117 template<
class Type>
2127 template<
class Type>
2137 contents_(contents),
2142 template<
class Type>
2147 const label maxLevels,
2148 const scalar maxLeafRatio,
2149 const scalar maxDuplicity
2160 Pout<<
"indexedOctree::indexedOctree:" <<
nl 2161 <<
" shapes:" << shapes.size() <<
nl 2162 <<
" bb:" << bb <<
nl 2167 if (shapes.size() == 0)
2173 DynamicList<node> nodes(label(shapes.size() / maxLeafRatio));
2174 DynamicList<labelList> contents(label(shapes.size() / maxLeafRatio));
2175 contents.push_back(
identity(shapes.size()));
2178 node topNode(
divide(bb, contents, 0));
2179 nodes.push_back(topNode);
2187 for (; nLevels < maxLevels; nLevels++)
2191 label minEntries = shapes.size();
2192 label maxEntries = 0;
2193 for (
const auto& subContents : contents)
2195 const label num = subContents.size();
2198 minEntries =
min(minEntries, num);
2199 maxEntries =
max(maxEntries, num);
2204 Pout<<
"indexedOctree::indexedOctree:" <<
nl 2205 <<
" nLevels:" << nLevels <<
nl 2206 <<
" nEntries:" << nEntries <<
nl 2207 <<
" - min per leaf: " << minEntries <<
nl 2208 <<
" - max per leaf: " << maxEntries <<
nl 2209 <<
" - avg per leaf: " 2210 << scalar(nEntries)/contents.size() <<
nl 2211 <<
" - per shape (duplicity): " 2212 << scalar(nEntries)/shapes_.size() <<
nl 2220 nEntries > maxDuplicity*shapes.size()
2228 label nOldNodes = nodes.size();
2231 label(maxLeafRatio),
2236 if (nOldNodes == nodes.size())
2250 contents_.setSize(contents.size());
2257 label nNodes = compactContents
2268 if (compactI == 0 && nNodes == 0)
2274 if (compactI == contents_.size())
2282 nodes_.transfer(nodes);
2288 label minEntries = shapes.size();
2289 label maxEntries = 0;
2290 for (
const auto& subContents : contents_)
2292 const label num = subContents.size();
2295 minEntries =
min(minEntries, num);
2296 maxEntries =
max(maxEntries, num);
2299 label memSize = memInfo().size();
2301 Pout<<
"indexedOctree::indexedOctree" 2302 <<
" : finished construction of tree of:" << shapes.typeName
2304 <<
" bb:" << this->bb() <<
nl 2305 <<
" shapes:" << shapes.size() <<
nl 2306 <<
" nLevels:" << nLevels <<
nl 2307 <<
" treeNodes:" << nodes_.size() <<
nl 2308 <<
" nEntries:" << nEntries <<
nl 2309 <<
" - min per leaf:" << minEntries <<
nl 2310 <<
" - max per leaf:" << maxEntries <<
nl 2311 <<
" - avg per leaf:" 2312 << scalar(nEntries)/contents.size() <<
nl 2313 <<
" - per shape (duplicity):" 2314 << scalar(nEntries)/shapes.size() <<
nl 2315 <<
" total memory:" << memSize-oldMemSize <<
nl 2321 template<
class Type>
2337 template<
class Type>
2340 const point& sample,
2341 const scalar startDistSqr
2348 typename Type::findNearestOp(*
this)
2353 template<
class Type>
2354 template<
class FindNearestOp>
2357 const point& sample,
2358 const scalar startDistSqr,
2360 const FindNearestOp& fnOp
2363 scalar nearestDistSqr = startDistSqr;
2364 label nearestShapeI = -1;
2382 return pointIndexHit(nearestShapeI != -1, nearestPoint, nearestShapeI);
2386 template<
class Type>
2399 typename Type::findNearestOp(*
this)
2404 template<
class Type>
2405 template<
class FindNearestOp>
2412 const FindNearestOp& fnOp
2415 label nearestShapeI = -1;
2434 return pointIndexHit(nearestShapeI != -1, nearestPoint, nearestShapeI);
2438 template<
class Type>
2450 typename Type::findIntersectOp(*
this)
2455 template<
class Type>
2467 typename Type::findIntersectOp(*
this)
2472 template<
class Type>
2473 template<
class FindIntersectOp>
2478 const FindIntersectOp& fiOp
2481 return findLine(
false, start,
end, fiOp);
2485 template<
class Type>
2486 template<
class FindIntersectOp>
2491 const FindIntersectOp& fiOp
2494 return findLine(
true, start,
end, fiOp);
2498 template<
class Type>
2505 return !nodes_.empty() && findBox(0, searchBox,
nullptr);
2509 template<
class Type>
2518 if (!nodes_.empty())
2523 label estimatedCapacity(
max(256, 2*(shapes_.size() / 100)));
2524 elements.
resize(estimatedCapacity);
2528 findBox(0, searchBox, &elements);
2531 return elements.
size();
2535 template<
class Type>
2548 findBox(searchBox, elements);
2551 return elements.
toc();
2555 template<
class Type>
2558 const point& centre,
2559 const scalar radiusSqr
2563 return !nodes_.empty() && findSphere(0, centre, radiusSqr,
nullptr);
2567 template<
class Type>
2570 const point& centre,
2571 const scalar radiusSqr,
2577 if (!nodes_.empty())
2582 label estimatedCapacity(
max(256, 2*(shapes_.size() / 100)));
2583 elements.
resize(estimatedCapacity);
2587 findSphere(0, centre, radiusSqr, &elements);
2595 template<
class Type>
2598 const point& centre,
2599 const scalar radiusSqr
2609 findSphere(centre, radiusSqr, elements);
2612 return elements.
toc();
2616 template<
class Type>
2626 return nodePlusOctant(nodeI, 0);
2629 const node& nod = nodes_[nodeI];
2631 direction octant = nod.bb_.subOctant(sample);
2633 labelBits index = nod.subNodes_[octant];
2638 return findNode(getNode(index), sample);
2640 else if (isContent(index))
2643 return nodePlusOctant(nodeI, octant);
2648 return nodePlusOctant(nodeI, octant);
2653 template<
class Type>
2661 labelBits index = findNode(0, sample);
2663 const node& nod = nodes_[getNode(index)];
2665 labelBits contentIndex = nod.subNodes_[
getOctant(index)];
2668 if (isContent(contentIndex))
2670 labelList indices = contents_[getContent(contentIndex)];
2674 label shapeI = indices[elemI];
2676 if (shapes_.contains(shapeI, sample))
2687 template<
class Type>
2695 return labelList::null();
2698 labelBits index = findNode(0, sample);
2700 const node& nod = nodes_[getNode(index)];
2702 labelBits contentIndex = nod.subNodes_[
getOctant(index)];
2705 if (isContent(contentIndex))
2707 return contents_[getContent(contentIndex)];
2710 return labelList::null();
2714 template<
class Type>
2722 return volumeType::UNKNOWN;
2725 if (nodeTypes_.size() != 8*nodes_.size())
2729 nodeTypes_.setSize(8*nodes_.size());
2730 nodeTypes_ = volumeType::UNKNOWN;
2745 if (
type == volumeType::UNKNOWN)
2749 else if (
type == volumeType::MIXED)
2753 else if (
type == volumeType::INSIDE)
2757 else if (
type == volumeType::OUTSIDE)
2767 Pout<<
"indexedOctree::getVolumeType : " 2769 <<
" nodes_:" << nodes_.size()
2770 <<
" nodeTypes_:" << nodeTypes_.size()
2771 <<
" nUNKNOWN:" << nUNKNOWN
2772 <<
" nMIXED:" << nMIXED
2773 <<
" nINSIDE:" << nINSIDE
2774 <<
" nOUTSIDE:" << nOUTSIDE
2779 return getVolumeType(0, sample);
2783 template<
class Type>
2784 template<
class CompareOp>
2787 const scalar nearDist,
2792 if (!nodes_.empty())
2799 nodePlusOctant(0, 0),
2802 nodePlusOctant(0, 0),
2810 template<
class Type>
2814 const bool printContents,
2823 const node& nod = nodes_[nodeI];
2824 const treeBoundBox& bb = nod.bb_;
2826 os <<
"nodeI:" << nodeI <<
" bb:" << bb <<
nl 2827 <<
"parent:" << nod.parent_ <<
nl 2828 <<
"n:" << countElements(nodePlusOctant(nodeI, 0)) <<
nl;
2830 for (
direction octant = 0; octant < node::nChildren; ++octant)
2832 const treeBoundBox subBb(bb.subBbox(octant));
2834 labelBits index = nod.subNodes_[octant];
2839 label subNodeI = getNode(index);
2841 os <<
"octant:" << octant
2842 <<
" node: n:" << countElements(index)
2843 <<
" bb:" << subBb <<
endl;
2845 string oldPrefix =
os.prefix();
2846 os.prefix() =
" " + oldPrefix;
2848 print(
os, printContents, subNodeI);
2850 os.prefix() = oldPrefix;
2852 else if (isContent(index))
2854 const labelList& indices = contents_[getContent(index)];
2861 os <<
"octant:" << octant
2862 <<
" content: n:" << indices.size()
2870 os <<
' ' << indices[j];
2882 os <<
"octant:" << octant <<
" empty:" << subBb <<
endl;
2888 template<
class Type>
2891 if (nodes_.size() < 2)
2894 return nodes_.size();
2897 return countLeafs(0);
2901 template<
class Type>
2910 template<
class Type>
2911 Foam::Ostream& Foam::operator<<(Ostream& os, const indexedOctree<Type>& t)
2914 os << t.bb() << token::SPACE << t.nodes()
2915 << token::SPACE << t.contents();
void divide(DimensionedField< Type, GeoMesh > &result, const DimensionedField< Type, GeoMesh > &f1, const DimensionedField< scalar, GeoMesh > &f2)
bool overlaps(const treeBoundBox &bb) const
True if any shapes overlap the bounding box.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
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.
pointIndexHit findLineAny(const point &start, const point &end) const
Find any intersection of line between start and end.
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.
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.
int size() const noexcept
Memory size (VmSize in /proc/PID/status) at last update()
This class describes the interaction of an object (often a face) and a point. It carries the info of ...
static volumeType getSide(const vector &outsideNormal, const vector &vec)
Helper function to return the side. Returns outside if.
An enumeration wrapper for classification of a location as being inside/outside of a volume...
labelBits findNode(const label nodeI, const point &) const
Find deepest node (as parent+octant) containing point. Starts.
UList< label > labelUList
A UList of labels.
void resize(const label sz)
Resize the hash table for efficiency.
Various functions to operate on Lists.
#define forAll(list, i)
Loop across all elements in list.
direction posBits(const point &pt) const
Position of point relative to bounding box.
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
label size() const noexcept
The number of elements in table.
bool intersects(const point &overallStart, const vector &overallVec, const point &start, const point &end, point &pt, direction &ptBits) const
Intersects segment; set point to intersection position and face,.
label capacity() const noexcept
The size of the underlying table.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
line< point, const point & > linePointRef
A line using referred points.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
A class for handling words, derived from Foam::string.
void clear()
Clear all entries from table.
const treeBoundBox & bb() const
Top bounding box.
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.
indexedOctree(const Type &shapes)
Construct null.
An Ostream is an abstract base class for all output systems (streams, files, token lists...
void print(prefixOSstream &, const bool printContents, const label) const
Print tree. Either print all indices (printContent = true) or.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
int debug
Static debugging option.
label findInside(const point &) const
Find shape containing point. Only implemented for certain.
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)
Memory usage information for the current process, and the system memory that is free.
const labelList & findIndices(const point &) const
Find the shape indices that occupy the result of findNode.
vector point
Point is a vector.
Non-pointer based hierarchical recursive searching.
#define WarningInFunction
Report a warning using Foam::Warning.
PtrList< volScalarField > & Y
Standard boundBox with extra functionality for use in octree.
A 29bits (or 61bits) integer label with 3bits direction (eg, octant) packed into single label...
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 nLeafs() const
Return the number of leaf nodes.
bool write(Ostream &os) const
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Tree node. Has up pointer and down pointers.
static constexpr const zero Zero
Global zero (0)