46 const scalar nearestDistSqr,
54 for (
direction dir = 0; dir < vector::nComponents; dir++)
56 scalar d0 = p0[dir] - sample[dir];
57 scalar d1 = p1[dir] - sample[dir];
59 if ((d0 > 0) != (d1 > 0))
64 else if (
mag(d0) <
mag(d1))
73 if (distSqr > nearestDistSqr)
90 const scalar nearestDistSqr,
109 if (octant & treeBoundBox::RIGHTHALF)
118 if (octant & treeBoundBox::TOPHALF)
127 if (octant & treeBoundBox::FRONTHALF)
136 const point mid(0.5*(min+max));
138 return overlaps(mid, other, nearestDistSqr, sample);
156 for (
direction octant = 0; octant < 8; octant++)
169 for (
direction octant = 0; octant < 8; octant++)
171 subBbs[octant] = bb.
subBbox(octant);
176 label shapeI = indices()[i];
178 for (
direction octant = 0; octant < 8; octant++)
180 if (shapes_.overlaps(shapeI, subBbs[octant]))
182 result[octant]().append(shapeI);
195 const label contentI,
196 const label parentNodeIndex,
197 const label octantToBeDivided
206 bb.
min()[0] >= bb.
max()[0]
207 || bb.
min()[1] >= bb.
max()[1]
208 || bb.
min()[2] >= bb.
max()[2]
212 <<
"Badly formed bounding box:" << bb
220 divide(indices, bb, dividedIndices);
225 bool replaced =
false;
227 for (
direction octant = 0; octant < dividedIndices.
size(); octant++)
231 if (subIndices().size())
235 contents_[contentI]().transfer(subIndices());
236 nod.
subNodes_[octant] = contentPlusOctant(contentI, octant);
244 label sz = contents_.size();
254 contents_[sz]().transfer(subIndices());
256 nod.
subNodes_[octant] = contentPlusOctant(sz, octant);
262 nod.
subNodes_[octant] = emptyPlusOctant(octant);
267 if (parentNodeIndex != -1)
271 label sz = nodes_.size();
275 nodes_[parentNodeIndex].subNodes_[octantToBeDivided]
276 = nodePlusOctant(sz, octantToBeDivided);
287 const label contentI,
288 const label parentIndex,
295 contents_[contentI]().size() > minSize_
296 && nLevels < maxLevels_
300 node nod =
divide(subBb, contentI, parentIndex, octant);
307 for (
direction subOct = 0; subOct < 8; subOct++)
311 if (isContent(subNodeLabel))
315 const label subContentI = getContent(subNodeLabel);
317 const label parentNodeIndex = nodes_.size() - 1;
342 const node& nod = nodes_[nodeI];
355 subType = calcVolumeType(getNode(index));
357 else if (isContent(index))
361 subType = volumeType::MIXED;
371 shapes_.getVolumeType(*
this, subBb.
midpoint())
376 nodeTypes_.set((nodeI<<3)+octant, subType);
380 if (myType == volumeType::UNKNOWN)
384 else if (subType != myType)
386 myType = volumeType::MIXED;
400 const node& nod = nodes_[nodeI];
406 if (octantType == volumeType::INSIDE)
410 else if (octantType == volumeType::OUTSIDE)
414 else if (octantType == volumeType::UNKNOWN)
419 else if (octantType == volumeType::MIXED)
426 volumeType subType = getVolumeType(getNode(index), sample);
430 else if (isContent(index))
433 return volumeType(shapes_.getVolumeType(*
this, sample));
441 "dynamicIndexedOctree<Type>::getVolumeType" 442 "(const label, const point&)" 443 ) <<
"Sample:" << sample <<
" node:" << nodeI
444 <<
" with bb:" << nodes_[nodeI].bb_ <<
nl 445 <<
"Empty subnode has invalid volume type MIXED." 448 return volumeType::UNKNOWN;
455 "dynamicIndexedOctree<Type>::getVolumeType" 456 "(const label, const point&)" 457 ) <<
"Sample:" << sample <<
" at node:" << nodeI
458 <<
" octant:" << octant
460 <<
"Node has invalid volume type " << octantType
463 return volumeType::UNKNOWN;
471 const vector& outsideNormal,
475 if ((outsideNormal&vec) >= 0)
477 return volumeType::OUTSIDE;
481 return volumeType::INSIDE;
498 scalar& nearestDistSqr,
499 label& nearestShapeI,
503 const node& nod = nodes_[nodeI];
518 label subNodeI = getNode(index);
522 if (overlaps(subBb.
min(), subBb.
max(), nearestDistSqr, sample))
535 else if (isContent(index))
550 contents_[getContent(index)],
571 label& nearestShapeI,
576 const node& nod = nodes_[nodeI];
608 else if (isContent(index))
616 contents_[getContent(index)],
633 const label parentNodeI,
638 const node& nod = nodes_[parentNodeI];
644 return nodes_[getNode(index)].bb_;
661 const bool pushInside
665 const vector perturbVec = perturbTol_*bb.
span();
667 point perturbedPt(pt);
674 for (
direction dir = 0; dir < vector::nComponents; dir++)
676 if (
mag(pt[dir]-bb.
min()[dir]) <
mag(perturbVec[dir]))
679 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
680 perturbedPt[dir] = bb.
min()[dir] + perturbDist;
682 else if (
mag(pt[dir]-bb.
max()[dir]) <
mag(perturbVec[dir]))
685 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
686 perturbedPt[dir] = bb.
max()[dir] - perturbDist;
692 for (
direction dir = 0; dir < vector::nComponents; dir++)
694 if (
mag(pt[dir]-bb.
min()[dir]) <
mag(perturbVec[dir]))
696 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
697 perturbedPt[dir] = bb.
min()[dir] - perturbDist;
699 else if (
mag(pt[dir]-bb.
max()[dir]) <
mag(perturbVec[dir]))
701 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
702 perturbedPt[dir] = bb.
max()[dir] + perturbDist;
709 if (pushInside != bb.
contains(perturbedPt))
711 FatalErrorIn(
"dynamicIndexedOctree<Type>::pushPoint(..)")
712 <<
"pushed point:" << pt
713 <<
" to:" << perturbedPt
714 <<
" wanted side:" << pushInside
715 <<
" obtained side:" << bb.
contains(perturbedPt)
733 const bool pushInside
737 const vector perturbVec = perturbTol_*bb.
span();
739 point perturbedPt(pt);
746 FatalErrorIn(
"dynamicIndexedOctree<Type>::pushPoint(..)")
750 if (faceID & treeBoundBox::LEFTBIT)
754 perturbedPt[0] = bb.
min()[0] + (perturbVec[0] + ROOTVSMALL);
758 perturbedPt[0] = bb.
min()[0] - (perturbVec[0] + ROOTVSMALL);
761 else if (faceID & treeBoundBox::RIGHTBIT)
765 perturbedPt[0] = bb.
max()[0] - (perturbVec[0] + ROOTVSMALL);
769 perturbedPt[0] = bb.
max()[0] + (perturbVec[0] + ROOTVSMALL);
773 if (faceID & treeBoundBox::BOTTOMBIT)
777 perturbedPt[1] = bb.
min()[1] + (perturbVec[1] + ROOTVSMALL);
781 perturbedPt[1] = bb.
min()[1] - (perturbVec[1] + ROOTVSMALL);
784 else if (faceID & treeBoundBox::TOPBIT)
788 perturbedPt[1] = bb.
max()[1] - (perturbVec[1] + ROOTVSMALL);
792 perturbedPt[1] = bb.
max()[1] + (perturbVec[1] + ROOTVSMALL);
796 if (faceID & treeBoundBox::BACKBIT)
800 perturbedPt[2] = bb.
min()[2] + (perturbVec[2] + ROOTVSMALL);
804 perturbedPt[2] = bb.
min()[2] - (perturbVec[2] + ROOTVSMALL);
807 else if (faceID & treeBoundBox::FRONTBIT)
811 perturbedPt[2] = bb.
max()[2] - (perturbVec[2] + ROOTVSMALL);
815 perturbedPt[2] = bb.
max()[2] + (perturbVec[2] + ROOTVSMALL);
821 if (pushInside != bb.
contains(perturbedPt))
823 FatalErrorIn(
"dynamicIndexedOctree<Type>::pushPoint(..)")
824 <<
"pushed point:" << pt <<
" on face:" << faceString(faceID)
825 <<
" to:" << perturbedPt
826 <<
" wanted side:" << pushInside
827 <<
" obtained side:" << bb.
contains(perturbedPt)
852 FatalErrorIn(
"dynamicIndexedOctree<Type>::pushPointIntoFace(..)")
853 <<
" bb:" << bb <<
endl 868 if (ptFaceID & treeBoundBox::LEFTBIT)
870 faceIndices[nFaces++] = treeBoundBox::LEFT;
872 else if (ptFaceID & treeBoundBox::RIGHTBIT)
874 faceIndices[nFaces++] = treeBoundBox::RIGHT;
877 if (ptFaceID & treeBoundBox::BOTTOMBIT)
879 faceIndices[nFaces++] = treeBoundBox::BOTTOM;
881 else if (ptFaceID & treeBoundBox::TOPBIT)
883 faceIndices[nFaces++] = treeBoundBox::TOP;
886 if (ptFaceID & treeBoundBox::BACKBIT)
888 faceIndices[nFaces++] = treeBoundBox::BACK;
890 else if (ptFaceID & treeBoundBox::FRONTBIT)
892 faceIndices[nFaces++] = treeBoundBox::FRONT;
905 else if (nFaces == 1)
908 keepFaceID = faceIndices[0];
915 keepFaceID = faceIndices[0];
916 scalar maxInproduct =
mag(treeBoundBox::faceNormals[keepFaceID] & dir);
921 scalar
s =
mag(treeBoundBox::faceNormals[face] & dir);
922 if (s > maxInproduct)
938 if (keepFaceID == treeBoundBox::LEFT)
940 facePoint.
x() = bb.
min().
x();
941 faceID = treeBoundBox::LEFTBIT;
943 else if (keepFaceID == treeBoundBox::RIGHT)
945 facePoint.
x() = bb.
max().
x();
946 faceID = treeBoundBox::RIGHTBIT;
948 else if (keepFaceID == treeBoundBox::BOTTOM)
950 facePoint.
y() = bb.
min().
y();
951 faceID = treeBoundBox::BOTTOMBIT;
953 else if (keepFaceID == treeBoundBox::TOP)
955 facePoint.
y() = bb.
max().
y();
956 faceID = treeBoundBox::TOPBIT;
958 else if (keepFaceID == treeBoundBox::BACK)
960 facePoint.
z() = bb.
min().
z();
961 faceID = treeBoundBox::BACKBIT;
963 else if (keepFaceID == treeBoundBox::FRONT)
965 facePoint.
z() = bb.
max().
z();
966 faceID = treeBoundBox::FRONTBIT;
972 if (faceID != bb.
faceBits(facePoint))
974 FatalErrorIn(
"dynamicIndexedOctree<Type>::pushPointIntoFace(..)")
975 <<
"Pushed point from " << pt
976 <<
" on face:" << ptFaceID <<
" of bb:" << bb <<
endl 977 <<
"onto " << facePoint
978 <<
" on face:" << faceID
979 <<
" which is not consistent with geometric face " 983 if (bb.
posBits(facePoint) != 0)
985 FatalErrorIn(
"dynamicIndexedOctree<Type>::pushPointIntoFace(..)")
986 <<
" bb:" << bb <<
endl 987 <<
"does not contain perturbed point " 1196 template<
class Type>
1206 parentNodeI = nodes_[nodeI].parent_;
1208 if (parentNodeI == -1)
1214 const node& parentNode = nodes_[parentNodeI];
1223 if (isNode(index) && getNode(index) == nodeI)
1230 if (parentOctant == 255)
1233 <<
"Problem: no parent found for octant:" << octant
1234 <<
" node:" << nodeI
1246 template<
class Type>
1249 const point& facePoint,
1255 label oldNodeI = nodeI;
1263 const direction X = treeBoundBox::RIGHTHALF;
1265 const direction Z = treeBoundBox::FRONTHALF;
1270 if ((faceID & treeBoundBox::LEFTBIT) != 0)
1276 else if ((faceID & treeBoundBox::RIGHTBIT) != 0)
1281 if ((faceID & treeBoundBox::BOTTOMBIT) != 0)
1287 else if ((faceID & treeBoundBox::TOPBIT) != 0)
1293 if ((faceID & treeBoundBox::BACKBIT) != 0)
1298 else if ((faceID & treeBoundBox::FRONTBIT) != 0)
1330 while (wantedValue != (octant & octantMask))
1336 if (wantedValue & X)
1356 if (wantedValue & Y)
1373 if (wantedValue & Z)
1393 walkToParent(nodeI, octant, parentNodeI, parentOctant);
1395 if (parentNodeI == -1)
1410 nodeI = parentNodeI;
1411 octant = parentOctant;
1417 octant ^= octantMask;
1429 FatalErrorIn(
"dynamicIndexedOctree<Type>::walkToNeighbour(..)")
1430 <<
"When searching for " << facePoint
1431 <<
" ended up in node:" << nodeI
1432 <<
" octant:" << octant
1433 <<
" with bb:" << subBb
1441 labelBits index = nodes_[nodeI].subNodes_[octant];
1447 nodeI = getNode(node);
1448 octant = getOctant(node);
1456 if (nodeI == oldNodeI && octant == oldOctant)
1458 FatalErrorIn(
"dynamicIndexedOctree<Type>::walkToNeighbour(..)")
1459 <<
"Did not go to neighbour when searching for " << facePoint
1461 <<
" starting from face:" << faceString(faceID)
1462 <<
" node:" << nodeI
1463 <<
" octant:" << octant
1470 FatalErrorIn(
"dynamicIndexedOctree<Type>::walkToNeighbour(..)")
1471 <<
"When searching for " << facePoint
1472 <<
" ended up in node:" << nodeI
1473 <<
" octant:" << octant
1484 template<
class Type>
1496 if (faceID & treeBoundBox::LEFTBIT)
1498 if (!desc.empty()) desc +=
"+";
1501 if (faceID & treeBoundBox::RIGHTBIT)
1503 if (!desc.empty()) desc +=
"+";
1506 if (faceID & treeBoundBox::BOTTOMBIT)
1508 if (!desc.empty()) desc +=
"+";
1511 if (faceID & treeBoundBox::TOPBIT)
1513 if (!desc.empty()) desc +=
"+";
1516 if (faceID & treeBoundBox::BACKBIT)
1518 if (!desc.empty()) desc +=
"+";
1521 if (faceID & treeBoundBox::FRONTBIT)
1523 if (!desc.empty()) desc +=
"+";
1536 template<
class Type>
1540 const point& treeStart,
1556 if (octantBb.
posBits(start) != 0)
1558 FatalErrorIn(
"dynamicIndexedOctree<Type>::traverseNode(..)")
1559 <<
"Node:" << nodeI <<
" octant:" << octant
1560 <<
" bb:" << octantBb <<
endl 1566 const node& nod = nodes_[nodeI];
1570 if (isContent(index))
1572 const labelList& indices = contents_[getContent(index)];
1582 label shapeI = indices[elemI];
1585 bool hit = shapes_.intersects(shapeI, start, end, pt);
1606 point nearestPoint(end);
1610 label shapeI = indices[elemI];
1613 bool hit = shapes_.intersects
1679 point perturbedEnd(pushPoint(octantBb, end,
false));
1699 template<
class Type>
1703 const point& treeStart,
1704 const point& treeEnd,
1705 const label startNodeI,
1710 const vector treeVec(treeEnd - treeStart);
1713 label nodeI = startNodeI;
1718 Pout<<
"findLine : treeStart:" << treeStart
1719 <<
" treeEnd:" << treeEnd <<
endl 1721 <<
" octant:" << octant
1722 <<
" bb:" << subBbox(nodeI, octant) <<
endl;
1730 for (; i < 100000; i++)
1751 <<
" at current:" << hitInfo.
rawPoint()
1752 <<
" (perturbed:" << startPoint <<
")" <<
endl 1753 <<
" node:" << nodeI
1754 <<
" octant:" << octant
1755 <<
" bb:" << subBbox(nodeI, octant) <<
endl;
1786 if (hitFaceID == 0 || hitInfo.
rawPoint() == treeEnd)
1793 point perturbedPoint
1806 Pout<<
" iter:" << i
1807 <<
" hit face:" << faceString(hitFaceID)
1809 <<
" node:" << nodeI
1810 <<
" octant:" << octant
1811 <<
" bb:" << subBbox(nodeI, octant) <<
nl 1812 <<
" walking to neighbour containing:" << perturbedPoint
1821 bool ok = walkToNeighbour
1840 <<
" to neighbour node:" << nodeI
1841 <<
" octant:" << octant
1843 <<
" of octantBb:" << octantBb <<
endl 1866 FatalErrorIn(
"dynamicIndexedOctree<Type>::findLine(..)")
1867 <<
"Got stuck in loop raytracing from:" << treeStart
1868 <<
" to:" << treeEnd <<
endl 1869 <<
"inside top box:" << subBbox(startNodeI, startOctant)
1874 WarningIn(
"dynamicIndexedOctree<Type>::findLine(..)")
1875 <<
"Got stuck in loop raytracing from:" << treeStart
1876 <<
" to:" << treeEnd <<
endl 1877 <<
"inside top box:" << subBbox(startNodeI, startOctant)
1887 template<
class Type>
1907 if ((startBit & endBit) != 0)
1916 point trackStart(start);
1917 point trackEnd(end);
1922 if (!treeBb.
intersects(start, end, trackStart))
1931 if (!treeBb.
intersects(end, trackStart, trackEnd))
1939 labelBits index = findNode(0, trackStart);
1941 label parentNodeI = getNode(index);
1958 template<
class Type>
1966 const node& nod = nodes_[nodeI];
1975 const treeBoundBox& subBb = nodes_[getNode(index)].bb_;
1979 findBox(getNode(index), searchBox, elements);
1982 else if (isContent(index))
1988 const labelList& indices = contents_[getContent(index)];
1992 label shapeI = indices[i];
1994 if (shapes_.overlaps(shapeI, searchBox))
2005 template<
class Type>
2009 const point& centre,
2010 const scalar radiusSqr,
2014 const node& nod = nodes_[nodeI];
2023 const treeBoundBox& subBb = nodes_[getNode(index)].bb_;
2025 if (subBb.
overlaps(centre, radiusSqr))
2027 findSphere(getNode(index), centre, radiusSqr, elements);
2030 else if (isContent(index))
2034 if (subBb.
overlaps(centre, radiusSqr))
2036 const labelList& indices = contents_[getContent(index)];
2040 label shapeI = indices[i];
2042 if (shapes_.overlaps(shapeI, centre, radiusSqr))
2053 template<
class Type>
2054 template<
class CompareOp>
2057 const scalar nearDist,
2068 const vector nearSpan(nearDist, nearDist, nearDist);
2070 if (tree1.
isNode(index1))
2079 if (tree2.
isNode(index2))
2146 if (tree2.
isNode(index2))
2189 label shape1 = indices1[i];
2193 label shape2 = indices2[j];
2195 if ((&tree1 != &tree2) || (shape1 != shape2))
2228 template<
class Type>
2239 label nodeI = getNode(index);
2241 const node& nod = nodes_[nodeI];
2245 nElems += countElements(nod.
subNodes_[octant]);
2248 else if (isContent(index))
2250 nElems += contents_[getContent(index)]().size();
2261 template<
class Type>
2273 labelBits index = nodes_[nodeI].subNodes_[octant];
2279 subBb = nodes_[getNode(index)].bb_;
2281 else if (isContent(index) || isEmpty(index))
2283 subBb = nodes_[nodeI].bb_.
subBbox(octant);
2286 Pout<<
"dumpContentNode : writing node:" << nodeI <<
" octant:" << octant
2294 const point& pt = bbPoints[i];
2296 str<<
"v " << pt.
x() <<
' ' << pt.
y() <<
' ' << pt.
z() <<
endl;
2299 forAll(treeBoundBox::edges, i)
2301 const edge&
e = treeBoundBox::edges[i];
2303 str<<
"l " << e[0] + 1 <<
' ' << e[1] + 1 <<
nl;
2310 template<
class Type>
2315 const label maxLevels,
2316 const scalar maxLeafRatio,
2317 const scalar maxDuplicity
2322 maxLevels_(maxLevels),
2324 maxLeafRatio_(maxLeafRatio),
2325 minSize_(
label(maxLeafRatio)),
2326 maxDuplicity_(maxDuplicity),
2327 nodes_(
label(shapes.size() / maxLeafRatio_)),
2328 contents_(
label(shapes.size() / maxLeafRatio_)),
2331 if (shapes_.size() == 0)
2336 insert(0, shapes_.size());
2347 template<
class Type>
2354 template<
class Type>
2357 const point& sample,
2358 const scalar startDistSqr
2361 scalar nearestDistSqr = startDistSqr;
2362 label nearestShapeI = -1;
2363 point nearestPoint = vector::zero;
2378 return pointIndexHit(nearestShapeI != -1, nearestPoint, nearestShapeI);
2382 template<
class Type>
2390 label nearestShapeI = -1;
2408 nearestPoint = vector::zero;
2411 return pointIndexHit(nearestShapeI != -1, nearestPoint, nearestShapeI);
2416 template<
class Type>
2423 return findLine(
false, start, end);
2428 template<
class Type>
2435 return findLine(
true, start, end);
2439 template<
class Type>
2450 findBox(0, searchBox, elements);
2453 return elements.
toc();
2457 template<
class Type>
2460 const point& centre,
2461 const scalar radiusSqr
2469 findSphere(0, centre, radiusSqr, elements);
2472 return elements.
toc();
2477 template<
class Type>
2487 return nodePlusOctant(nodeI, 0);
2490 const node& nod = nodes_[nodeI];
2497 <<
"Cannot find " << sample <<
" in node " << nodeI
2509 return findNode(getNode(index), sample);
2511 else if (isContent(index))
2514 return nodePlusOctant(nodeI, octant);
2519 return nodePlusOctant(nodeI, octant);
2524 template<
class Type>
2532 const node& nod = nodes_[getNode(index)];
2537 if (isContent(contentIndex))
2539 labelList indices = contents_[getContent(contentIndex)];
2543 label shapeI = indices[elemI];
2545 if (shapes_.contains(shapeI, sample))
2556 template<
class Type>
2564 const node& nod = nodes_[getNode(index)];
2569 if (isContent(contentIndex))
2571 return contents_[getContent(contentIndex)];
2575 return emptyList<label>();
2581 template<
class Type>
2589 return volumeType::UNKNOWN;
2592 if (nodeTypes_.size() != 8*nodes_.size())
2596 nodeTypes_.setSize(8*nodes_.size());
2597 nodeTypes_ = volumeType::UNKNOWN;
2612 if (type == volumeType::UNKNOWN)
2616 else if (type == volumeType::MIXED)
2620 else if (type == volumeType::INSIDE)
2624 else if (type == volumeType::OUTSIDE)
2634 Pout<<
"dynamicIndexedOctree<Type>::getVolumeType : " 2636 <<
" nodes_:" << nodes_.size()
2637 <<
" nodeTypes_:" << nodeTypes_.size()
2638 <<
" nUNKNOWN:" << nUNKNOWN
2639 <<
" nMIXED:" << nMIXED
2640 <<
" nINSIDE:" << nINSIDE
2641 <<
" nOUTSIDE:" << nOUTSIDE
2646 return getVolumeType(0, sample);
2650 template<
class Type>
2651 template<
class CompareOp>
2654 const scalar nearDist,
2664 nodePlusOctant(0, 0),
2667 nodePlusOctant(0, 0),
2674 template<
class Type>
2677 if (startIndex == endIndex)
2692 contents_[0]().append(0);
2697 nodes_.append(topNode);
2704 for (
label pI = startIndex; pI < endIndex; ++pI)
2708 if (!insertIndex(0, pI, nLevels))
2713 nLevelsMax_ =
max(nLevels, nLevelsMax_);
2720 template<
class Type>
2723 const label nodIndex,
2728 bool shapeInserted =
false;
2730 for (
direction octant = 0; octant < 8; octant++)
2732 const labelBits& subNodeLabel = nodes_[nodIndex].subNodes_[octant];
2734 if (isNode(subNodeLabel))
2736 const treeBoundBox& subBb = nodes_[getNode(subNodeLabel)].bb_;
2738 if (shapes().overlaps(index, subBb))
2742 if (insertIndex(getNode(subNodeLabel), index, nLevels))
2744 shapeInserted =
true;
2748 else if (isContent(subNodeLabel))
2752 if (shapes().overlaps(index, subBb))
2754 const label contentI = getContent(subNodeLabel);
2756 contents_[contentI]().append(index);
2758 recursiveSubDivision
2767 shapeInserted =
true;
2774 if (shapes().overlaps(index, subBb))
2776 label sz = contents_.size();
2783 contents_[sz]().append(index);
2785 nodes_[nodIndex].subNodes_[octant]
2786 = contentPlusOctant(sz, octant);
2789 shapeInserted =
true;
2793 return shapeInserted;
2797 template<
class Type>
2805 removeIndex(0, index);
2811 template<
class Type>
2814 const label nodIndex,
2818 label totalContents = 0;
2820 for (
direction octant = 0; octant < 8; octant++)
2822 const labelBits& subNodeLabel = nodes_[nodIndex].subNodes_[octant];
2824 if (isNode(subNodeLabel))
2826 const treeBoundBox& subBb = nodes_[getNode(subNodeLabel)].bb_;
2828 if (shapes().overlaps(index, subBb))
2831 label childContentsSize
2832 = removeIndex(getNode(subNodeLabel), index);
2834 totalContents += childContentsSize;
2836 if (childContentsSize == 0)
2838 nodes_[nodIndex].subNodes_[octant]
2839 = emptyPlusOctant(octant);
2848 else if (isContent(subNodeLabel))
2852 const label contentI = getContent(subNodeLabel);
2854 if (shapes().overlaps(index, subBb))
2862 const label oldIndex = contentList[pI];
2864 if (oldIndex != index)
2866 newContent.
append(oldIndex);
2870 newContent.shrink();
2872 if (newContent.size() == 0)
2875 nodes_[nodIndex].subNodes_[octant]
2876 = emptyPlusOctant(octant);
2882 totalContents += contents_[contentI]().size();
2890 return totalContents;
2895 template<
class Type>
2899 const bool printContents,
2903 const node& nod = nodes_[nodeI];
2906 os <<
"nodeI:" << nodeI <<
" bb:" << bb <<
nl 2908 <<
"n:" << countElements(nodePlusOctant(nodeI, 0)) <<
nl;
2919 label subNodeI = getNode(index);
2921 os <<
"octant:" << octant
2922 <<
" node: n:" << countElements(index)
2923 <<
" bb:" << subBb <<
endl;
2925 string oldPrefix = os.
prefix();
2926 os.
prefix() =
" " + oldPrefix;
2928 print(os, printContents, subNodeI);
2932 else if (isContent(index))
2934 const labelList& indices = contents_[getContent(index)];
2941 os <<
"octant:" << octant
2942 <<
" content: n:" << indices.
size()
2950 os <<
' ' << indices[j];
2962 os <<
"octant:" << octant <<
" empty:" << subBb <<
endl;
2968 template<
class Type>
2974 nEntries += contents_[i]().size();
2977 Pout<<
"indexedOctree<Type>::indexedOctree" 2978 <<
" : finished construction of tree of:" << shapes().typeName
2980 <<
" bounding box: " << this->bb() <<
nl 2981 <<
" shapes: " << shapes().size() <<
nl 2982 <<
" treeNodes: " << nodes_.size() <<
nl 2983 <<
" nEntries: " << nEntries <<
nl 2984 <<
" levels/maxLevels: " << nLevelsMax_ <<
"/" << maxLevels_ <<
nl 2985 <<
" minSize: " << minSize_ <<
nl 2986 <<
" per treeLeaf: " 2987 << scalar(nEntries)/contents_.size() <<
nl 2988 <<
" per shape (duplicity):" 2989 << scalar(nEntries)/shapes().size() <<
nl 2995 template<
class Type>
3004 template<
class Type>
3006 Foam::operator<<(Ostream& os, const dynamicIndexedOctree<Type>& t)
3008 os << t.bb() << token::SPACE << t.nodes() <<
endl;
3012 os << t.contents()[cI]() <<
endl;
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,.
const string & prefix() const
Return the prefix of the stream.
const List< node > & nodes() const
List of all nodes.
labelBits findNode(const label nodeI, const point &) const
Find deepest node (as parent+octant) containing point. Starts.
void transfer(List< T > &)
Transfer contents of the argument List into this.
const point & min() const
Minimum describing the bounding box.
bool insert(label startIndex, label endIndex)
Insert a new object into the tree.
bool remove(const label index)
Remove an object from the tree.
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().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 label getContent(const labelBits i)
dimensioned< scalar > mag(const dimensioned< Type > &)
dynamicIndexedOctree(const Type &shapes, const treeBoundBox &bb, const label maxLevels, const scalar maxLeafRatio, const scalar maxDuplicity )
Construct from shapes.
word name(const complex &)
Return a string representation of a complex.
bool contains(const vector &dir, const point &) const
Contains point (inside or on edge) and moving in direction.
pointIndexHit findLineAny(const point &start, const point &end) const
Find any intersection of line between start and end.
const Point & rawPoint() const
Return point with no checking.
bool write(Ostream &os) const
static bool isNode(const labelBits i)
vector span() const
The bounding box span (from minimum to maximum)
const Type & shapes() const
Reference to shape.
direction subOctant(const point &pt) const
Returns octant number given point and the calculated midpoint.
A class for handling words, derived from string.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Standard boundBox + extra functionality for use in octree.
label removeIndex(const label nodIndex, const label index)
const labelList & findIndices(const point &) const
Find the shape indices that occupy the result of findNode.
A 29bits label and 3bits direction packed into single label.
void size(const label)
Override size to be inconsistent with allocated storage.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Various functions to operate on Lists.
PointIndexHit< point > pointIndexHit
bool hit() const
Is there a hit.
const contentListList & contents() const
List of all contents (referenced by those nodes that are.
static scalar & perturbTol()
Get the perturbation tolerance.
A 1D vector of objects of type <T> with a fixed size <Size>.
tmp< pointField > points() const
Vertex coordinates. In octant coding.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
timeIndices insert(timeIndex, timeDirs[timeI].value())
label facePoint(const int facei, const block &block, const label i, const label j)
treeBoundBox bb_
Bounding box of this node.
static volumeType getSide(const vector &outsideNormal, const vector &vec)
Helper function to return the side. Returns outside if.
A face is a list of labels corresponding to mesh vertices.
static label getNode(const labelBits i)
Ostream & endl(Ostream &os)
Add newline and flush stream.
#define WarningIn(functionName)
Report a warning using Foam::Warning.
label parent_
Parent node (index into nodes_ of tree)
const dimensionedScalar e
Elementary charge.
void searchOrder(const point &pt, FixedList< direction, 8 > &octantOrder) const
Calculates optimal order to look for nearest to point.
PtrList< volScalarField > & Y
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
const fileName & name() const
Return the name of the stream.
FixedList< labelBits, 8 > subNodes_
IDs of the 8 nodes on all sides of the mid point.
bool good() const
Return true if next operation might succeed.
direction faceBits(const point &) const
Code position of point on bounding box faces.
const point & max() const
Maximum describing the bounding box.
List< Key > toc() const
Return the table of contents.
errorManip< error > abort(error &err)
void setIndex(const label index)
direction posBits(const point &) const
Position of point relative to bounding box.
fileName::Type type(const fileName &)
Return the file type: DIRECTORY or FILE.
static bool isContent(const labelBits i)
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
point midpoint() const
The midpoint of the bounding box.
void writeTreeInfo() const
bool insertIndex(const label nodIndex, const label index, label &nLevels)
void setPoint(const Point &p)
Point centre() const
Return centre (centroid)
An Ostream is an abstract base class for all output systems (streams, files, token lists...
treeBoundBox subBbox(const direction) const
Sub box given by octant number. Midpoint calculated.
Version of OSstream which prints a prefix on each line.
Non-pointer based hierarchical recursive searching. Storage is dynamic, so elements can be deleted...
const treeBoundBox & bb() const
Top bounding box.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
label findInside(const point &) const
Find shape containing point. Only implemented for certain.
void print(prefixOSstream &, const bool printContents, const label) const
Print tree. Either print all indices (printContent = true) or.
bool overlaps(const boundBox &) const
Overlaps/touches boundingBox?
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Tree node. Has up pointer and down pointers.
void divide(FieldField< Field, Type > &f, const FieldField< Field, Type > &f1, const FieldField< Field, scalar > &f2)
prefixOSstream Pout(cout,"Pout")
bool insert(const Key &key)
Insert a new entry.