43 const scalar nearestDistSqr,
51 for (
direction dir = 0; dir < vector::nComponents; dir++)
53 scalar d0 = p0[dir] - sample[dir];
54 scalar d1 = p1[dir] - sample[dir];
56 if ((d0 > 0) != (d1 > 0))
61 else if (
mag(d0) <
mag(d1))
70 if (distSqr > nearestDistSqr)
85 const scalar nearestDistSqr,
104 if (octant & treeBoundBox::RIGHTHALF)
113 if (octant & treeBoundBox::TOPHALF)
122 if (octant & treeBoundBox::FRONTHALF)
131 const point mid(0.5*(min+max));
133 return overlaps(mid, other, nearestDistSqr, sample);
145 for (
direction octant = 0; octant < 8; octant++)
158 for (
direction octant = 0; octant < 8; octant++)
160 subBbs[octant] = bb.
subBbox(octant);
165 label shapeI = indices()[i];
167 for (
direction octant = 0; octant < 8; octant++)
169 if (shapes_.overlaps(shapeI, subBbs[octant]))
171 result[octant]().append(shapeI);
183 const label contentI,
184 const label parentNodeIndex,
185 const label octantToBeDivided
194 bb.
min()[0] >= bb.
max()[0]
195 || bb.
min()[1] >= bb.
max()[1]
196 || bb.
min()[2] >= bb.
max()[2]
200 <<
"Badly formed bounding box:" << bb
208 divide(indices, bb, dividedIndices);
213 bool replaced =
false;
215 for (
direction octant = 0; octant < dividedIndices.
size(); octant++)
219 if (subIndices().size())
223 contents_[contentI]().transfer(subIndices());
224 nod.
subNodes_[octant] = contentPlusOctant(contentI, octant);
232 label sz = contents_.size();
242 contents_[sz]().transfer(subIndices());
244 nod.
subNodes_[octant] = contentPlusOctant(sz, octant);
250 nod.
subNodes_[octant] = emptyPlusOctant(octant);
255 if (parentNodeIndex != -1)
259 label sz = nodes_.size();
263 nodes_[parentNodeIndex].subNodes_[octantToBeDivided]
264 = nodePlusOctant(sz, octantToBeDivided);
275 const label contentI,
276 const label parentIndex,
283 contents_[contentI]().size() > minSize_
284 && nLevels < maxLevels_
288 node nod =
divide(subBb, contentI, parentIndex, octant);
295 for (
direction subOct = 0; subOct < 8; subOct++)
299 if (isContent(subNodeLabel))
303 const label subContentI = getContent(subNodeLabel);
305 const label parentNodeIndex = nodes_.size() - 1;
331 const node& nod = nodes_[nodeI];
344 subType = calcVolumeType(getNode(index));
346 else if (isContent(index))
350 subType = volumeType::MIXED;
360 shapes_.getVolumeType(*
this, subBb.
midpoint())
365 nodeTypes_.set((nodeI<<3)+octant, subType);
369 if (myType == volumeType::UNKNOWN)
373 else if (subType != myType)
375 myType = volumeType::MIXED;
389 const node& nod = nodes_[nodeI];
395 if (octantType == volumeType::INSIDE)
399 else if (octantType == volumeType::OUTSIDE)
403 else if (octantType == volumeType::UNKNOWN)
408 else if (octantType == volumeType::MIXED)
415 volumeType subType = getVolumeType(getNode(index), sample);
419 else if (isContent(index))
422 return volumeType(shapes_.getVolumeType(*
this, sample));
429 <<
"Sample:" << sample <<
" node:" << nodeI
430 <<
" with bb:" << nodes_[nodeI].bb_ <<
nl 431 <<
"Empty subnode has invalid volume type MIXED." 434 return volumeType::UNKNOWN;
440 <<
"Sample:" << sample <<
" at node:" << nodeI
441 <<
" octant:" << octant
443 <<
"Node has invalid volume type " << octantType
446 return volumeType::UNKNOWN;
454 const vector& outsideNormal,
458 if ((outsideNormal&vec) >= 0)
460 return volumeType::OUTSIDE;
464 return volumeType::INSIDE;
475 scalar& nearestDistSqr,
476 label& nearestShapeI,
480 const node& nod = nodes_[nodeI];
495 label subNodeI = getNode(index);
499 if (overlaps(subBb.
min(), subBb.
max(), nearestDistSqr, sample))
512 else if (isContent(index))
527 contents_[getContent(index)],
547 label& nearestShapeI,
552 const node& nod = nodes_[nodeI];
584 else if (isContent(index))
592 contents_[getContent(index)],
609 const label parentNodeI,
614 const node& nod = nodes_[parentNodeI];
620 return nodes_[getNode(index)].bb_;
635 const bool pushInside
642 const vector perturbVec = perturbTol_*bb.
span();
644 point perturbedPt(pt);
651 for (
direction dir = 0; dir < vector::nComponents; dir++)
653 if (
mag(pt[dir]-bb.
min()[dir]) <
mag(perturbVec[dir]))
656 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
657 perturbedPt[dir] = bb.
min()[dir] + perturbDist;
659 else if (
mag(pt[dir]-bb.
max()[dir]) <
mag(perturbVec[dir]))
662 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
663 perturbedPt[dir] = bb.
max()[dir] - perturbDist;
669 for (
direction dir = 0; dir < vector::nComponents; dir++)
671 if (
mag(pt[dir]-bb.
min()[dir]) <
mag(perturbVec[dir]))
673 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
674 perturbedPt[dir] = bb.
min()[dir] - perturbDist;
676 else if (
mag(pt[dir]-bb.
max()[dir]) <
mag(perturbVec[dir]))
678 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
679 perturbedPt[dir] = bb.
max()[dir] + perturbDist;
686 if (pushInside != bb.
contains(perturbedPt))
689 <<
"pushed point:" << pt
690 <<
" to:" << perturbedPt
691 <<
" wanted side:" << pushInside
692 <<
" obtained side:" << bb.
contains(perturbedPt)
708 const bool pushInside
715 const vector perturbVec = perturbTol_*bb.
span();
717 point perturbedPt(pt);
728 if (faceID & treeBoundBox::LEFTBIT)
732 perturbedPt[0] = bb.
min()[0] + (perturbVec[0] + ROOTVSMALL);
736 perturbedPt[0] = bb.
min()[0] - (perturbVec[0] + ROOTVSMALL);
739 else if (faceID & treeBoundBox::RIGHTBIT)
743 perturbedPt[0] = bb.
max()[0] - (perturbVec[0] + ROOTVSMALL);
747 perturbedPt[0] = bb.
max()[0] + (perturbVec[0] + ROOTVSMALL);
751 if (faceID & treeBoundBox::BOTTOMBIT)
755 perturbedPt[1] = bb.
min()[1] + (perturbVec[1] + ROOTVSMALL);
759 perturbedPt[1] = bb.
min()[1] - (perturbVec[1] + ROOTVSMALL);
762 else if (faceID & treeBoundBox::TOPBIT)
766 perturbedPt[1] = bb.
max()[1] - (perturbVec[1] + ROOTVSMALL);
770 perturbedPt[1] = bb.
max()[1] + (perturbVec[1] + ROOTVSMALL);
774 if (faceID & treeBoundBox::BACKBIT)
778 perturbedPt[2] = bb.
min()[2] + (perturbVec[2] + ROOTVSMALL);
782 perturbedPt[2] = bb.
min()[2] - (perturbVec[2] + ROOTVSMALL);
785 else if (faceID & treeBoundBox::FRONTBIT)
789 perturbedPt[2] = bb.
max()[2] - (perturbVec[2] + ROOTVSMALL);
793 perturbedPt[2] = bb.
max()[2] + (perturbVec[2] + ROOTVSMALL);
799 if (pushInside != bb.
contains(perturbedPt))
802 <<
"pushed point:" << pt <<
" on face:" << faceString(faceID)
803 <<
" to:" << perturbedPt
804 <<
" wanted side:" << pushInside
805 <<
" obtained side:" << bb.
contains(perturbedPt)
828 <<
" bb:" << bb <<
endl 843 if (ptFaceID & treeBoundBox::LEFTBIT)
845 faceIndices[nFaces++] = treeBoundBox::LEFT;
847 else if (ptFaceID & treeBoundBox::RIGHTBIT)
849 faceIndices[nFaces++] = treeBoundBox::RIGHT;
852 if (ptFaceID & treeBoundBox::BOTTOMBIT)
854 faceIndices[nFaces++] = treeBoundBox::BOTTOM;
856 else if (ptFaceID & treeBoundBox::TOPBIT)
858 faceIndices[nFaces++] = treeBoundBox::TOP;
861 if (ptFaceID & treeBoundBox::BACKBIT)
863 faceIndices[nFaces++] = treeBoundBox::BACK;
865 else if (ptFaceID & treeBoundBox::FRONTBIT)
867 faceIndices[nFaces++] = treeBoundBox::FRONT;
880 else if (nFaces == 1)
883 keepFaceID = faceIndices[0];
890 keepFaceID = faceIndices[0];
891 scalar maxInproduct =
mag(treeBoundBox::faceNormals[keepFaceID] & dir);
896 scalar
s =
mag(treeBoundBox::faceNormals[face] & dir);
897 if (s > maxInproduct)
913 if (keepFaceID == treeBoundBox::LEFT)
915 facePoint.
x() = bb.
min().
x();
916 faceID = treeBoundBox::LEFTBIT;
918 else if (keepFaceID == treeBoundBox::RIGHT)
920 facePoint.
x() = bb.
max().
x();
921 faceID = treeBoundBox::RIGHTBIT;
923 else if (keepFaceID == treeBoundBox::BOTTOM)
925 facePoint.
y() = bb.
min().
y();
926 faceID = treeBoundBox::BOTTOMBIT;
928 else if (keepFaceID == treeBoundBox::TOP)
930 facePoint.
y() = bb.
max().
y();
931 faceID = treeBoundBox::TOPBIT;
933 else if (keepFaceID == treeBoundBox::BACK)
935 facePoint.
z() = bb.
min().
z();
936 faceID = treeBoundBox::BACKBIT;
938 else if (keepFaceID == treeBoundBox::FRONT)
940 facePoint.
z() = bb.
max().
z();
941 faceID = treeBoundBox::FRONTBIT;
947 if (faceID != bb.
faceBits(facePoint))
950 <<
"Pushed point from " << pt
951 <<
" on face:" << ptFaceID <<
" of bb:" << bb <<
endl 952 <<
"onto " << facePoint
953 <<
" on face:" << faceID
954 <<
" which is not consistent with geometric face " 958 if (bb.
posBits(facePoint) != 0)
961 <<
" bb:" << bb <<
endl 962 <<
"does not contain perturbed point " 982 parentNodeI = nodes_[nodeI].parent_;
984 if (parentNodeI == -1)
990 const node& parentNode = nodes_[parentNodeI];
999 if (isNode(index) && getNode(index) == nodeI)
1006 if (parentOctant == 255)
1009 <<
"Problem: no parent found for octant:" << octant
1010 <<
" node:" << nodeI
1019 template<
class Type>
1022 const point& facePoint,
1033 label oldNodeI = nodeI;
1041 const direction X = treeBoundBox::RIGHTHALF;
1043 const direction Z = treeBoundBox::FRONTHALF;
1048 if ((faceID & treeBoundBox::LEFTBIT) != 0)
1054 else if ((faceID & treeBoundBox::RIGHTBIT) != 0)
1059 if ((faceID & treeBoundBox::BOTTOMBIT) != 0)
1065 else if ((faceID & treeBoundBox::TOPBIT) != 0)
1071 if ((faceID & treeBoundBox::BACKBIT) != 0)
1076 else if ((faceID & treeBoundBox::FRONTBIT) != 0)
1108 while (wantedValue != (octant & octantMask))
1114 if (wantedValue & X)
1134 if (wantedValue & Y)
1151 if (wantedValue & Z)
1171 walkToParent(nodeI, octant, parentNodeI, parentOctant);
1173 if (parentNodeI == -1)
1188 nodeI = parentNodeI;
1189 octant = parentOctant;
1195 octant ^= octantMask;
1208 <<
"When searching for " << facePoint
1209 <<
" ended up in node:" << nodeI
1210 <<
" octant:" << octant
1211 <<
" with bb:" << subBb
1219 labelBits index = nodes_[nodeI].subNodes_[octant];
1225 nodeI = getNode(node);
1226 octant = getOctant(node);
1234 if (nodeI == oldNodeI && octant == oldOctant)
1237 <<
"Did not go to neighbour when searching for " << facePoint
1239 <<
" starting from face:" << faceString(faceID)
1240 <<
" node:" << nodeI
1241 <<
" octant:" << octant
1249 <<
"When searching for " << facePoint
1250 <<
" ended up in node:" << nodeI
1251 <<
" octant:" << octant
1262 template<
class Type>
1274 if (faceID & treeBoundBox::LEFTBIT)
1276 if (!desc.empty()) desc +=
"+";
1279 if (faceID & treeBoundBox::RIGHTBIT)
1281 if (!desc.empty()) desc +=
"+";
1284 if (faceID & treeBoundBox::BOTTOMBIT)
1286 if (!desc.empty()) desc +=
"+";
1289 if (faceID & treeBoundBox::TOPBIT)
1291 if (!desc.empty()) desc +=
"+";
1294 if (faceID & treeBoundBox::BACKBIT)
1296 if (!desc.empty()) desc +=
"+";
1299 if (faceID & treeBoundBox::FRONTBIT)
1301 if (!desc.empty()) desc +=
"+";
1308 template<
class Type>
1312 const point& treeStart,
1328 if (octantBb.
posBits(start) != 0)
1331 <<
"Node:" << nodeI <<
" octant:" << octant
1332 <<
" bb:" << octantBb <<
endl 1338 const node& nod = nodes_[nodeI];
1342 if (isContent(index))
1344 const labelList& indices = contents_[getContent(index)];
1354 label shapeI = indices[elemI];
1357 bool hit = shapes_.intersects(shapeI, start, end, pt);
1378 point nearestPoint(end);
1382 label shapeI = indices[elemI];
1385 bool hit = shapes_.intersects
1451 point perturbedEnd(pushPoint(octantBb, end,
false));
1470 template<
class Type>
1474 const point& treeStart,
1475 const point& treeEnd,
1476 const label startNodeI,
1481 const vector treeVec(treeEnd - treeStart);
1484 label nodeI = startNodeI;
1489 Pout<<
"findLine : treeStart:" << treeStart
1490 <<
" treeEnd:" << treeEnd <<
endl 1492 <<
" octant:" << octant
1493 <<
" bb:" << subBbox(nodeI, octant) <<
endl;
1501 for (; i < 100000; i++)
1522 <<
" at current:" << hitInfo.
rawPoint()
1523 <<
" (perturbed:" << startPoint <<
")" <<
endl 1524 <<
" node:" << nodeI
1525 <<
" octant:" << octant
1526 <<
" bb:" << subBbox(nodeI, octant) <<
endl;
1557 if (hitFaceID == 0 || hitInfo.
rawPoint() == treeEnd)
1564 point perturbedPoint
1577 Pout<<
" iter:" << i
1578 <<
" hit face:" << faceString(hitFaceID)
1580 <<
" node:" << nodeI
1581 <<
" octant:" << octant
1582 <<
" bb:" << subBbox(nodeI, octant) <<
nl 1583 <<
" walking to neighbour containing:" << perturbedPoint
1592 bool ok = walkToNeighbour
1611 <<
" to neighbour node:" << nodeI
1612 <<
" octant:" << octant
1614 <<
" of octantBb:" << octantBb <<
endl 1638 <<
"Got stuck in loop raytracing from:" << treeStart
1639 <<
" to:" << treeEnd <<
endl 1640 <<
"inside top box:" << subBbox(startNodeI, startOctant)
1646 <<
"Got stuck in loop raytracing from:" << treeStart
1647 <<
" to:" << treeEnd <<
endl 1648 <<
"inside top box:" << subBbox(startNodeI, startOctant)
1657 template<
class Type>
1677 if ((startBit & endBit) != 0)
1686 point trackStart(start);
1687 point trackEnd(end);
1692 if (!treeBb.
intersects(start, end, trackStart))
1701 if (!treeBb.
intersects(end, trackStart, trackEnd))
1709 labelBits index = findNode(0, trackStart);
1711 label parentNodeI = getNode(index);
1728 template<
class Type>
1736 const node& nod = nodes_[nodeI];
1745 const treeBoundBox& subBb = nodes_[getNode(index)].bb_;
1749 findBox(getNode(index), searchBox, elements);
1752 else if (isContent(index))
1758 const labelList& indices = contents_[getContent(index)];
1762 label shapeI = indices[i];
1764 if (shapes_.overlaps(shapeI, searchBox))
1775 template<
class Type>
1779 const point& centre,
1780 const scalar radiusSqr,
1784 const node& nod = nodes_[nodeI];
1793 const treeBoundBox& subBb = nodes_[getNode(index)].bb_;
1795 if (subBb.
overlaps(centre, radiusSqr))
1797 findSphere(getNode(index), centre, radiusSqr, elements);
1800 else if (isContent(index))
1804 if (subBb.
overlaps(centre, radiusSqr))
1806 const labelList& indices = contents_[getContent(index)];
1810 label shapeI = indices[i];
1812 if (shapes_.overlaps(shapeI, centre, radiusSqr))
1823 template<
class Type>
1824 template<
class CompareOp>
1827 const scalar nearDist,
1838 const vector nearSpan(nearDist, nearDist, nearDist);
1840 if (tree1.
isNode(index1))
1849 if (tree2.
isNode(index2))
1916 if (tree2.
isNode(index2))
1959 label shape1 = indices1[i];
1963 label shape2 = indices2[j];
1965 if ((&tree1 != &tree2) || (shape1 != shape2))
1997 template<
class Type>
2008 label nodeI = getNode(index);
2010 const node& nod = nodes_[nodeI];
2014 nElems += countElements(nod.
subNodes_[octant]);
2017 else if (isContent(index))
2019 nElems += contents_[getContent(index)]().size();
2030 template<
class Type>
2042 labelBits index = nodes_[nodeI].subNodes_[octant];
2048 subBb = nodes_[getNode(index)].bb_;
2050 else if (isContent(index) || isEmpty(index))
2052 subBb = nodes_[nodeI].bb_.
subBbox(octant);
2055 Pout<<
"dumpContentNode : writing node:" << nodeI <<
" octant:" << octant
2063 const point& pt = bbPoints[i];
2065 str<<
"v " << pt.
x() <<
' ' << pt.
y() <<
' ' << pt.
z() <<
endl;
2068 forAll(treeBoundBox::edges, i)
2070 const edge&
e = treeBoundBox::edges[i];
2072 str<<
"l " << e[0] + 1 <<
' ' << e[1] + 1 <<
nl;
2079 template<
class Type>
2084 const label maxLevels,
2085 const scalar maxLeafRatio,
2086 const scalar maxDuplicity
2091 maxLevels_(maxLevels),
2093 maxLeafRatio_(maxLeafRatio),
2094 minSize_(
label(maxLeafRatio)),
2095 maxDuplicity_(maxDuplicity),
2096 nodes_(
label(shapes.size() / maxLeafRatio_)),
2097 contents_(
label(shapes.size() / maxLeafRatio_)),
2100 if (shapes_.size() == 0)
2105 insert(0, shapes_.size());
2116 template<
class Type>
2123 template<
class Type>
2126 const point& sample,
2127 const scalar startDistSqr
2130 scalar nearestDistSqr = startDistSqr;
2131 label nearestShapeI = -1;
2147 return pointIndexHit(nearestShapeI != -1, nearestPoint, nearestShapeI);
2151 template<
class Type>
2159 label nearestShapeI = -1;
2177 nearestPoint =
Zero;
2180 return pointIndexHit(nearestShapeI != -1, nearestPoint, nearestShapeI);
2184 template<
class Type>
2191 return findLine(
false, start, end);
2195 template<
class Type>
2202 return findLine(
true, start, end);
2206 template<
class Type>
2217 findBox(0, searchBox, elements);
2220 return elements.
toc();
2224 template<
class Type>
2227 const point& centre,
2228 const scalar radiusSqr
2236 findSphere(0, centre, radiusSqr, elements);
2239 return elements.
toc();
2243 template<
class Type>
2253 return nodePlusOctant(nodeI, 0);
2256 const node& nod = nodes_[nodeI];
2263 <<
"Cannot find " << sample <<
" in node " << nodeI
2275 return findNode(getNode(index), sample);
2277 else if (isContent(index))
2280 return nodePlusOctant(nodeI, octant);
2285 return nodePlusOctant(nodeI, octant);
2290 template<
class Type>
2298 const node& nod = nodes_[getNode(index)];
2303 if (isContent(contentIndex))
2305 labelList indices = contents_[getContent(contentIndex)];
2309 label shapeI = indices[elemI];
2311 if (shapes_.contains(shapeI, sample))
2322 template<
class Type>
2330 const node& nod = nodes_[getNode(index)];
2335 if (isContent(contentIndex))
2337 return contents_[getContent(contentIndex)];
2341 return emptyList<label>();
2346 template<
class Type>
2354 return volumeType::UNKNOWN;
2357 if (nodeTypes_.size() != 8*nodes_.size())
2361 nodeTypes_.setSize(8*nodes_.size());
2362 nodeTypes_ = volumeType::UNKNOWN;
2377 if (type == volumeType::UNKNOWN)
2381 else if (type == volumeType::MIXED)
2385 else if (type == volumeType::INSIDE)
2389 else if (type == volumeType::OUTSIDE)
2399 Pout<<
"dynamicIndexedOctree<Type>::getVolumeType : " 2401 <<
" nodes_:" << nodes_.size()
2402 <<
" nodeTypes_:" << nodeTypes_.size()
2403 <<
" nUNKNOWN:" << nUNKNOWN
2404 <<
" nMIXED:" << nMIXED
2405 <<
" nINSIDE:" << nINSIDE
2406 <<
" nOUTSIDE:" << nOUTSIDE
2411 return getVolumeType(0, sample);
2415 template<
class Type>
2416 template<
class CompareOp>
2419 const scalar nearDist,
2429 nodePlusOctant(0, 0),
2432 nodePlusOctant(0, 0),
2439 template<
class Type>
2442 if (startIndex == endIndex)
2457 contents_[0]().append(0);
2462 nodes_.append(topNode);
2469 for (
label pI = startIndex; pI < endIndex; ++pI)
2473 if (!insertIndex(0, pI, nLevels))
2478 nLevelsMax_ =
max(nLevels, nLevelsMax_);
2485 template<
class Type>
2488 const label nodIndex,
2493 bool shapeInserted =
false;
2495 for (
direction octant = 0; octant < 8; octant++)
2497 const labelBits& subNodeLabel = nodes_[nodIndex].subNodes_[octant];
2499 if (isNode(subNodeLabel))
2501 const treeBoundBox& subBb = nodes_[getNode(subNodeLabel)].bb_;
2503 if (shapes().overlaps(index, subBb))
2507 if (insertIndex(getNode(subNodeLabel), index, nLevels))
2509 shapeInserted =
true;
2513 else if (isContent(subNodeLabel))
2517 if (shapes().overlaps(index, subBb))
2519 const label contentI = getContent(subNodeLabel);
2521 contents_[contentI]().append(index);
2523 recursiveSubDivision
2532 shapeInserted =
true;
2539 if (shapes().overlaps(index, subBb))
2541 label sz = contents_.size();
2548 contents_[sz]().append(index);
2550 nodes_[nodIndex].subNodes_[octant]
2551 = contentPlusOctant(sz, octant);
2554 shapeInserted =
true;
2558 return shapeInserted;
2562 template<
class Type>
2570 removeIndex(0, index);
2576 template<
class Type>
2579 const label nodIndex,
2583 label totalContents = 0;
2585 for (
direction octant = 0; octant < 8; octant++)
2587 const labelBits& subNodeLabel = nodes_[nodIndex].subNodes_[octant];
2589 if (isNode(subNodeLabel))
2591 const treeBoundBox& subBb = nodes_[getNode(subNodeLabel)].bb_;
2593 if (shapes().overlaps(index, subBb))
2596 label childContentsSize
2597 = removeIndex(getNode(subNodeLabel), index);
2599 totalContents += childContentsSize;
2601 if (childContentsSize == 0)
2603 nodes_[nodIndex].subNodes_[octant]
2604 = emptyPlusOctant(octant);
2613 else if (isContent(subNodeLabel))
2617 const label contentI = getContent(subNodeLabel);
2619 if (shapes().overlaps(index, subBb))
2627 const label oldIndex = contentList[pI];
2629 if (oldIndex != index)
2631 newContent.
append(oldIndex);
2635 newContent.shrink();
2637 if (newContent.size() == 0)
2640 nodes_[nodIndex].subNodes_[octant]
2641 = emptyPlusOctant(octant);
2647 totalContents += contents_[contentI]().size();
2655 return totalContents;
2659 template<
class Type>
2663 const bool printContents,
2667 const node& nod = nodes_[nodeI];
2670 os <<
"nodeI:" << nodeI <<
" bb:" << bb <<
nl 2672 <<
"n:" << countElements(nodePlusOctant(nodeI, 0)) <<
nl;
2683 label subNodeI = getNode(index);
2685 os <<
"octant:" << octant
2686 <<
" node: n:" << countElements(index)
2687 <<
" bb:" << subBb <<
endl;
2689 string oldPrefix = os.
prefix();
2690 os.
prefix() =
" " + oldPrefix;
2692 print(os, printContents, subNodeI);
2696 else if (isContent(index))
2698 const labelList& indices = contents_[getContent(index)];
2705 os <<
"octant:" << octant
2706 <<
" content: n:" << indices.
size()
2714 os <<
' ' << indices[j];
2726 os <<
"octant:" << octant <<
" empty:" << subBb <<
endl;
2732 template<
class Type>
2738 nEntries += contents_[i]().size();
2741 Pout<<
"indexedOctree<Type>::indexedOctree" 2742 <<
" : finished construction of tree of:" << shapes().typeName
2744 <<
" bounding box: " << this->bb() <<
nl 2745 <<
" shapes: " << shapes().size() <<
nl 2746 <<
" treeNodes: " << nodes_.size() <<
nl 2747 <<
" nEntries: " << nEntries <<
nl 2748 <<
" levels/maxLevels: " << nLevelsMax_ <<
"/" << maxLevels_ <<
nl 2749 <<
" minSize: " << minSize_ <<
nl 2750 <<
" per treeLeaf: " 2751 << scalar(nEntries)/contents_.size() <<
nl 2752 <<
" per shape (duplicity):" 2753 << scalar(nEntries)/shapes().size() <<
nl 2758 template<
class Type>
2767 template<
class Type>
2769 Foam::operator<<(Ostream& os, const dynamicIndexedOctree<Type>& t)
2771 os << t.bb() << token::SPACE << t.nodes() <<
endl;
2775 os << t.contents()[cI]() <<
endl;
void divide(FieldField< Field, Type > &f, const FieldField< Field, Type > &f1, const FieldField< Field, scalar > &f2)
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Tree node. Has up pointer and down pointers.
A face is a list of labels corresponding to mesh vertices.
A 1D vector of objects of type <T> with a fixed size <Size>.
const point & min() const
Minimum describing the bounding box.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
bool write(Ostream &os) const
bool hit() const
Is there a hit.
label parent_
Parent node (index into nodes_ of tree)
void size(const label)
Override size to be inconsistent with allocated storage.
void writeTreeInfo() const
void setIndex(const label index)
bool good() const
Return true if next operation might succeed.
Ostream & endl(Ostream &os)
Add newline and flush stream.
bool insert(label startIndex, label endIndex)
Insert a new object into the tree.
const labelList & findIndices(const point &) const
Find the shape indices that occupy the result of findNode.
const point & max() const
Maximum describing the bounding box.
point midpoint() const
The midpoint of the bounding box.
const List< node > & nodes() const
List of all nodes.
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
PointIndexHit< point > pointIndexHit
bool remove(const label index)
Remove an object from the tree.
bool insert(const Key &key)
Insert a new entry.
void print(prefixOSstream &, const bool printContents, const label) const
Print tree. Either print all indices (printContent = true) or.
static bool isContent(const labelBits i)
void setPoint(const Point &p)
FixedList< labelBits, 8 > subNodes_
IDs of the 8 nodes on all sides of the mid point.
Various functions to operate on Lists.
direction faceBits(const point &) const
Code position of point on bounding box faces.
static label getContent(const labelBits i)
label removeIndex(const label nodIndex, const label index)
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,.
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))
direction posBits(const point &) const
Position of point relative to bounding box.
const dimensionedScalar e
Elementary charge.
void searchOrder(const point &pt, FixedList< direction, 8 > &octantOrder) const
Calculates optimal order to look for nearest to point.
Version of OSstream which prints a prefix on each line.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
const Type & shapes() const
Reference to shape.
A class for handling words, derived from string.
direction subOctant(const point &pt) const
Returns octant number given point and the calculated midpoint.
timeIndices insert(timeIndex, timeDirs[timeI].value())
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Non-pointer based hierarchical recursive searching. Storage is dynamic, so elements can be deleted...
label findInside(const point &) const
Find shape containing point. Only implemented for certain.
treeBoundBox bb_
Bounding box of this node.
vector span() const
The bounding box span (from minimum to maximum)
const Point & rawPoint() const
Return point with no checking.
errorManip< error > abort(error &err)
List< Key > toc() const
Return the table of contents.
An Ostream is an abstract base class for all output systems (streams, files, token lists...
prefixOSstream Pout(cout,"Pout")
bool insertIndex(const label nodIndex, const label index, label &nLevels)
label facePoint(const int facei, const block &block, const label i, const label j)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
word name(const complex &)
Return a string representation of a complex.
dynamicIndexedOctree(const Type &shapes, const treeBoundBox &bb, const label maxLevels, const scalar maxLeafRatio, const scalar maxDuplicity)
Construct from shapes.
bool contains(const vector &dir, const point &) const
Contains point (inside or on edge) and moving in direction.
static scalar & perturbTol()
Get the perturbation tolerance.
static volumeType getSide(const vector &outsideNormal, const vector &vec)
Helper function to return the side. Returns outside if.
fileName::Type type(const fileName &)
Return the file type: DIRECTORY or FILE.
#define WarningInFunction
Report a warning using Foam::Warning.
tmp< pointField > points() const
Vertex coordinates. In octant coding.
PtrList< volScalarField > & Y
Standard boundBox + extra functionality for use in octree.
A 29bits label and 3bits direction packed into single label.
dimensioned< scalar > mag(const dimensioned< Type > &)
static bool isNode(const labelBits i)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
labelBits findNode(const label nodeI, const point &) const
Find deepest node (as parent+octant) containing point. Starts.
const treeBoundBox & bb() const
Top bounding box.
treeBoundBox subBbox(const direction) const
Sub box given by octant number. Midpoint calculated.
const string & prefix() const
Return the prefix of the stream.
void transfer(List< T > &)
Transfer contents of the argument List into this.
static label getNode(const labelBits i)
const contentListList & contents() const
List of all contents (referenced by those nodes that are.
pointIndexHit findLineAny(const point &start, const point &end) const
Find any intersection of line between start and end.
bool overlaps(const boundBox &) const
Overlaps/touches boundingBox?
const fileName & name() const
Return the name of the stream.
Point centre() const
Return centre (centroid)