dynamicIndexedOctree.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9 Copyright (C) 2019 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
30#include "linePointRef.H"
31#include "OFstream.H"
32#include "ListOps.H"
33
34// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35
36template<class Type>
38
39// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
40
41template<class Type>
43(
44 const point& p0,
45 const point& p1,
46 const scalar nearestDistSqr,
47 const point& sample
48)
49{
50 // Find out where sample is in relation to bb.
51 // Find nearest point on bb.
52 scalar distSqr = 0;
53
54 for (direction dir = 0; dir < vector::nComponents; dir++)
55 {
56 scalar d0 = p0[dir] - sample[dir];
57 scalar d1 = p1[dir] - sample[dir];
58
59 if ((d0 > 0) != (d1 > 0))
60 {
61 // sample inside both extrema. This component does not add any
62 // distance.
63 }
64 else if (mag(d0) < mag(d1))
65 {
66 distSqr += d0*d0;
67 }
68 else
69 {
70 distSqr += d1*d1;
71 }
72
73 if (distSqr > nearestDistSqr)
74 {
75 return false;
76 }
77 }
78
79 return true;
80}
81
82
83template<class Type>
85(
86 const treeBoundBox& parentBb,
87 const direction octant,
88 const scalar nearestDistSqr,
89 const point& sample
90)
91{
92 //- Accelerated version of
93 // treeBoundBox subBb(parentBb.subBbox(mid, octant))
94 // overlaps
95 // (
96 // subBb.min(),
97 // subBb.max(),
98 // nearestDistSqr,
99 // sample
100 // )
101
102 const point& min = parentBb.min();
103 const point& max = parentBb.max();
104
105 point other;
106
107 if (octant & treeBoundBox::RIGHTHALF)
108 {
109 other.x() = max.x();
110 }
111 else
112 {
113 other.x() = min.x();
114 }
115
116 if (octant & treeBoundBox::TOPHALF)
117 {
118 other.y() = max.y();
119 }
120 else
121 {
122 other.y() = min.y();
123 }
124
125 if (octant & treeBoundBox::FRONTHALF)
126 {
127 other.z() = max.z();
128 }
129 else
130 {
131 other.z() = min.z();
132 }
133
134 const point mid(0.5*(min+max));
135
136 return overlaps(mid, other, nearestDistSqr, sample);
137}
138
139
140template<class Type>
142(
143 const autoPtr<DynamicList<label>>& indices,
144 const treeBoundBox& bb,
145 contentListList& result
146) const
147{
148 for (direction octant = 0; octant < 8; octant++)
149 {
150 result.append
151 (
152 autoPtr<DynamicList<label>>
153 (
154 new DynamicList<label>(indices().size()/8)
155 )
156 );
157 }
158
159 // Precalculate bounding boxes.
160 FixedList<treeBoundBox, 8> subBbs;
161 for (direction octant = 0; octant < 8; octant++)
162 {
163 subBbs[octant] = bb.subBbox(octant);
164 }
165
166 forAll(indices(), i)
167 {
168 label shapeI = indices()[i];
169
170 for (direction octant = 0; octant < 8; octant++)
171 {
172 if (shapes_.overlaps(shapeI, subBbs[octant]))
173 {
174 result[octant]().append(shapeI);
175 }
176 }
177 }
178}
179
180
181template<class Type>
184(
185 const treeBoundBox& bb,
186 const label contentI,
187 const label parentNodeIndex,
188 const label octantToBeDivided
189)
190{
191 const autoPtr<DynamicList<label>>& indices = contents_[contentI];
192
193 node nod;
194
195 if
196 (
197 bb.min()[0] >= bb.max()[0]
198 || bb.min()[1] >= bb.max()[1]
199 || bb.min()[2] >= bb.max()[2]
200 )
201 {
203 << "Badly formed bounding box:" << bb
204 << abort(FatalError);
205 }
206
207 nod.bb_ = bb;
208 nod.parent_ = -1;
209
210 contentListList dividedIndices(8);
211 divide(indices, bb, dividedIndices);
212
213 // Have now divided the indices into 8 (possibly empty) subsets.
214 // Replace current contentI with the first (non-empty) subset.
215 // Append the rest.
216 bool replaced = false;
217
218 for (direction octant = 0; octant < dividedIndices.size(); octant++)
219 {
220 autoPtr<DynamicList<label>>& subIndices = dividedIndices[octant];
221
222 if (subIndices().size())
223 {
224 if (!replaced)
225 {
226 contents_[contentI]->transfer(subIndices());
227 nod.subNodes_[octant] = contentPlusOctant(contentI, octant);
228
229 replaced = true;
230 }
231 else
232 {
233 // Store at end of contents.
234 // note dummy append + transfer trick
235 label sz = contents_.size();
236
237 contents_.append
238 (
239 autoPtr<DynamicList<label>>
240 (
241 new DynamicList<label>()
242 )
243 );
244
245 contents_[sz]->transfer(subIndices());
246
247 nod.subNodes_[octant] = contentPlusOctant(sz, octant);
248 }
249 }
250 else
251 {
252 // Mark node as empty
253 nod.subNodes_[octant] = emptyPlusOctant(octant);
254 }
255 }
256
257 // Don't update the parent node if it is the first node.
258 if (parentNodeIndex != -1)
259 {
260 nod.parent_ = parentNodeIndex;
261
262 label sz = nodes_.size();
263
264 nodes_.append(nod);
265
266 nodes_[parentNodeIndex].subNodes_[octantToBeDivided]
267 = nodePlusOctant(sz, octantToBeDivided);
268 }
269
270 return nod;
271}
272
273
274template<class Type>
276(
277 const treeBoundBox& subBb,
278 const label contentI,
279 const label parentIndex,
280 const label octant,
281 label& nLevels
282)
283{
284 if
285 (
286 contents_[contentI]->size() > minSize_
287 && nLevels < maxLevels_
288 )
289 {
290 // Create node for content
291 node nod = divide(subBb, contentI, parentIndex, octant);
292
293 // Increment the number of levels in the tree
294 nLevels++;
295
296 // Recursively divide the contents until maxLevels_ is
297 // reached or the content sizes are less than minSize_
298 for (direction subOct = 0; subOct < 8; subOct++)
299 {
300 const labelBits& subNodeLabel = nod.subNodes_[subOct];
301
302 if (isContent(subNodeLabel))
303 {
304 const treeBoundBox subBb = nod.bb_.subBbox(subOct);
305
306 const label subContentI = getContent(subNodeLabel);
307
308 const label parentNodeIndex = nodes_.size() - 1;
309
310 recursiveSubDivision
311 (
312 subBb,
313 subContentI,
314 parentNodeIndex,
315 subOct,
316 nLevels
317 );
318 }
319 }
320 }
321}
322
323
324template<class Type>
326(
327 const label nodeI
328) const
329{
330 // Pre-calculates wherever possible the volume status per node/subnode.
331 // Recurses to determine status of lowest level boxes. Level above is
332 // combination of octants below.
333
334 const node& nod = nodes_[nodeI];
335
336 volumeType myType = volumeType::UNKNOWN;
337
338 for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
339 {
340 volumeType subType;
341
342 labelBits index = nod.subNodes_[octant];
343
344 if (isNode(index))
345 {
346 // tree node. Recurse.
347 subType = calcVolumeType(getNode(index));
348 }
349 else if (isContent(index))
350 {
351 // Contents. Depending on position in box might be on either
352 // side.
353 subType = volumeType::MIXED;
354 }
355 else
356 {
357 // No data in this octant. Set type for octant acc. to the mid
358 // of its bounding box.
359 const treeBoundBox subBb = nod.bb_.subBbox(octant);
360
361 subType = volumeType
362 (
363 shapes_.getVolumeType(*this, subBb.centre())
364 );
365 }
366
367 // Store octant type
368 nodeTypes_.set((nodeI<<3)+octant, subType);
369
370 // Combine sub node types into type for treeNode. Result is 'mixed' if
371 // types differ among subnodes.
372 if (myType == volumeType::UNKNOWN)
373 {
374 myType = subType;
375 }
376 else if (subType != myType)
377 {
378 myType = volumeType::MIXED;
379 }
380 }
381 return myType;
382}
383
384
385template<class Type>
387(
388 const label nodeI,
389 const point& sample
390) const
391{
392 const node& nod = nodes_[nodeI];
393
394 direction octant = nod.bb_.subOctant(sample);
395
396 volumeType octantType = volumeType::type(nodeTypes_.get((nodeI<<3)+octant));
397
398 if (octantType == volumeType::INSIDE)
399 {
400 return octantType;
401 }
402 else if (octantType == volumeType::OUTSIDE)
404 return octantType;
405 }
406 else if (octantType == volumeType::UNKNOWN)
407 {
408 // Can happen for e.g. non-manifold surfaces.
409 return octantType;
410 }
411 else if (octantType == volumeType::MIXED)
412 {
413 labelBits index = nod.subNodes_[octant];
414
415 if (isNode(index))
416 {
417 // Recurse
418 volumeType subType = getVolumeType(getNode(index), sample);
419
420 return subType;
421 }
422 else if (isContent(index))
423 {
424 // Content. Defer to shapes.
425 return volumeType(shapes_.getVolumeType(*this, sample));
426 }
427
428 // Empty node. Cannot have 'mixed' as its type since not divided
429 // up and has no items inside it.
431 << "Sample:" << sample << " node:" << nodeI
432 << " with bb:" << nodes_[nodeI].bb_ << nl
433 << "Empty subnode has invalid volume type MIXED."
434 << abort(FatalError);
435
436 return volumeType::UNKNOWN;
437 }
438
440 << "Sample:" << sample << " at node:" << nodeI
441 << " octant:" << octant
442 << " with bb:" << nod.bb_.subBbox(octant) << nl
443 << "Node has invalid volume type " << octantType
444 << abort(FatalError);
445
446 return volumeType::UNKNOWN;
447}
448
449
450template<class Type>
452(
453 const vector& outsideNormal,
454 const vector& vec
455)
456{
457 if ((outsideNormal&vec) >= 0)
458 {
459 return volumeType::OUTSIDE;
460 }
461 else
462 {
463 return volumeType::INSIDE;
464 }
465}
466
467
468template<class Type>
470(
471 const label nodeI,
472 const point& sample,
473
474 scalar& nearestDistSqr,
475 label& nearestShapeI,
476 point& nearestPoint
477) const
478{
479 const node& nod = nodes_[nodeI];
480
481 // Determine order to walk through octants
482 FixedList<direction, 8> octantOrder;
483 nod.bb_.searchOrder(sample, octantOrder);
484
485 // Go into all suboctants (one containing sample first) and update nearest.
486 for (direction i = 0; i < 8; i++)
487 {
488 direction octant = octantOrder[i];
489
490 labelBits index = nod.subNodes_[octant];
491
492 if (isNode(index))
493 {
494 label subNodeI = getNode(index);
495
496 const treeBoundBox& subBb = nodes_[subNodeI].bb_;
497
498 if (overlaps(subBb.min(), subBb.max(), nearestDistSqr, sample))
499 {
500 findNearest
501 (
502 subNodeI,
503 sample,
504
505 nearestDistSqr,
506 nearestShapeI,
507 nearestPoint
508 );
509 }
510 }
511 else if (isContent(index))
512 {
513 if
514 (
515 overlaps
516 (
517 nod.bb_,
518 octant,
519 nearestDistSqr,
520 sample
521 )
522 )
523 {
524 shapes_.findNearest
525 (
526 *(contents_[getContent(index)]),
527 sample,
528
529 nearestDistSqr,
530 nearestShapeI,
531 nearestPoint
532 );
533 }
534 }
536}
537
538
539template<class Type>
541(
542 const label nodeI,
543 const linePointRef& ln,
544
545 treeBoundBox& tightest,
546 label& nearestShapeI,
547 point& linePoint,
548 point& nearestPoint
549) const
551 const node& nod = nodes_[nodeI];
552 const treeBoundBox& nodeBb = nod.bb_;
553
554 // Determine order to walk through octants
555 FixedList<direction, 8> octantOrder;
556 nod.bb_.searchOrder(ln.centre(), octantOrder);
557
558 // Go into all suboctants (one containing sample first) and update nearest.
559 for (direction i = 0; i < 8; i++)
560 {
561 direction octant = octantOrder[i];
562
563 labelBits index = nod.subNodes_[octant];
564
565 if (isNode(index))
566 {
567 const treeBoundBox& subBb = nodes_[getNode(index)].bb_;
568
569 if (subBb.overlaps(tightest))
570 {
571 findNearest
573 getNode(index),
574 ln,
575
576 tightest,
577 nearestShapeI,
578 linePoint,
579 nearestPoint
580 );
581 }
582 }
583 else if (isContent(index))
584 {
585 const treeBoundBox subBb(nodeBb.subBbox(octant));
586
587 if (subBb.overlaps(tightest))
588 {
589 shapes_.findNearest
590 (
591 *(contents_[getContent(index)]),
592 ln,
593
594 tightest,
595 nearestShapeI,
596 linePoint,
597 nearestPoint
598 );
599 }
600 }
601 }
602}
603
604
605template<class Type>
607(
608 const label parentNodeI,
609 const direction octant
610) const
611{
612 // Get type of node at octant
613 const node& nod = nodes_[parentNodeI];
614 labelBits index = nod.subNodes_[octant];
615
616 if (isNode(index))
618 // Use stored bb
619 return nodes_[getNode(index)].bb_;
620 }
621 else
622 {
623 // Calculate subBb
624 return nod.bb_.subBbox(octant);
625 }
626}
628
629template<class Type>
631(
632 const treeBoundBox& bb,
633 const point& pt,
634 const bool pushInside
635)
637 // Takes a bb and a point on/close to the edge of the bb and pushes the
638 // point inside by a small fraction.
639
640 // Get local length scale.
641 const vector perturbVec = perturbTol_*bb.span();
642
643 point perturbedPt(pt);
644
645 // Modify all components which are close to any face of the bb to be
646 // well inside/outside them.
647
648 if (pushInside)
649 {
650 for (direction dir = 0; dir < vector::nComponents; dir++)
651 {
652 if (mag(pt[dir]-bb.min()[dir]) < mag(perturbVec[dir]))
653 {
654 // Close to 'left' side. Push well beyond left side.
655 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
656 perturbedPt[dir] = bb.min()[dir] + perturbDist;
657 }
658 else if (mag(pt[dir]-bb.max()[dir]) < mag(perturbVec[dir]))
659 {
660 // Close to 'right' side. Push well beyond right side.
661 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
662 perturbedPt[dir] = bb.max()[dir] - perturbDist;
663 }
664 }
665 }
666 else
667 {
668 for (direction dir = 0; dir < vector::nComponents; dir++)
669 {
670 if (mag(pt[dir]-bb.min()[dir]) < mag(perturbVec[dir]))
671 {
672 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
673 perturbedPt[dir] = bb.min()[dir] - perturbDist;
674 }
675 else if (mag(pt[dir]-bb.max()[dir]) < mag(perturbVec[dir]))
676 {
677 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
678 perturbedPt[dir] = bb.max()[dir] + perturbDist;
679 }
680 }
681 }
682
683 if (debug)
684 {
685 if (pushInside != bb.contains(perturbedPt))
686 {
687 auto fatal = FatalErrorInFunction;
688
689 fatal
690 << "pushed point:" << pt
691 << " to:" << perturbedPt
692 << " wanted side:" << pushInside
693 << " obtained side:" << bb.contains(perturbedPt)
694 << " of bb:" << bb << nl;
695
696 if (debug > 1)
697 {
698 fatal << abort(FatalError);
699 }
700 }
701 }
702
703 return perturbedPt;
704}
705
706
707template<class Type>
709(
710 const treeBoundBox& bb,
711 const direction faceID,
712 const point& pt,
713 const bool pushInside
714)
715{
716 // Takes a bb and a point on the edge of the bb and pushes the point
717 // outside by a small fraction.
718
719 // Get local length scale.
720 const vector perturbVec = perturbTol_*bb.span();
721
722 point perturbedPt(pt);
723
724 // Modify all components which are close to any face of the bb to be
725 // well outside them.
726
727 if (faceID == 0)
728 {
730 << abort(FatalError);
731 }
732
733 if (faceID & treeBoundBox::LEFTBIT)
734 {
735 if (pushInside)
736 {
737 perturbedPt[0] = bb.min()[0] + (perturbVec[0] + ROOTVSMALL);
738 }
739 else
740 {
741 perturbedPt[0] = bb.min()[0] - (perturbVec[0] + ROOTVSMALL);
742 }
743 }
744 else if (faceID & treeBoundBox::RIGHTBIT)
745 {
746 if (pushInside)
747 {
748 perturbedPt[0] = bb.max()[0] - (perturbVec[0] + ROOTVSMALL);
749 }
750 else
751 {
752 perturbedPt[0] = bb.max()[0] + (perturbVec[0] + ROOTVSMALL);
753 }
754 }
755
756 if (faceID & treeBoundBox::BOTTOMBIT)
757 {
758 if (pushInside)
759 {
760 perturbedPt[1] = bb.min()[1] + (perturbVec[1] + ROOTVSMALL);
761 }
762 else
763 {
764 perturbedPt[1] = bb.min()[1] - (perturbVec[1] + ROOTVSMALL);
765 }
766 }
767 else if (faceID & treeBoundBox::TOPBIT)
768 {
769 if (pushInside)
770 {
771 perturbedPt[1] = bb.max()[1] - (perturbVec[1] + ROOTVSMALL);
772 }
773 else
774 {
775 perturbedPt[1] = bb.max()[1] + (perturbVec[1] + ROOTVSMALL);
776 }
777 }
778
779 if (faceID & treeBoundBox::BACKBIT)
780 {
781 if (pushInside)
782 {
783 perturbedPt[2] = bb.min()[2] + (perturbVec[2] + ROOTVSMALL);
784 }
785 else
786 {
787 perturbedPt[2] = bb.min()[2] - (perturbVec[2] + ROOTVSMALL);
788 }
789 }
790 else if (faceID & treeBoundBox::FRONTBIT)
791 {
792 if (pushInside)
793 {
794 perturbedPt[2] = bb.max()[2] - (perturbVec[2] + ROOTVSMALL);
795 }
796 else
797 {
798 perturbedPt[2] = bb.max()[2] + (perturbVec[2] + ROOTVSMALL);
799 }
800 }
801
802 if (debug)
803 {
804 if (pushInside != bb.contains(perturbedPt))
805 {
806 auto fatal = FatalErrorInFunction;
807
808 fatal
809 << "pushed point:" << pt << " on face:" << faceString(faceID)
810 << " to:" << perturbedPt
811 << " wanted side:" << pushInside
812 << " obtained side:" << bb.contains(perturbedPt)
813 << " of bb:" << bb << nl;
814
815 if (debug > 1)
816 {
817 fatal << abort(FatalError);
818 }
819 }
820 }
821
822 return perturbedPt;
823}
824
825
826template<class Type>
828(
829 const treeBoundBox& bb,
830 const vector& dir, // end-start
831 const point& pt
832)
833{
834 if (debug)
835 {
836 if (bb.posBits(pt) != 0)
837 {
838 auto fatal = FatalErrorInFunction;
839
840 fatal
841 << " bb:" << bb << endl
842 << "does not contain point " << pt << nl;
843
844 if (debug > 1)
845 {
846 fatal << abort(FatalError);
847 }
848 }
849 }
850
851
852 // Handle two cases:
853 // - point exactly on multiple faces. Push away from all but one.
854 // - point on a single face. Push away from edges of face.
855
856 direction ptFaceID = bb.faceBits(pt);
857
858 direction nFaces = 0;
859 FixedList<direction, 3> faceIndices;
860
861 if (ptFaceID & treeBoundBox::LEFTBIT)
862 {
863 faceIndices[nFaces++] = treeBoundBox::LEFT;
864 }
865 else if (ptFaceID & treeBoundBox::RIGHTBIT)
866 {
867 faceIndices[nFaces++] = treeBoundBox::RIGHT;
868 }
869
870 if (ptFaceID & treeBoundBox::BOTTOMBIT)
871 {
872 faceIndices[nFaces++] = treeBoundBox::BOTTOM;
873 }
874 else if (ptFaceID & treeBoundBox::TOPBIT)
875 {
876 faceIndices[nFaces++] = treeBoundBox::TOP;
877 }
878
879 if (ptFaceID & treeBoundBox::BACKBIT)
880 {
881 faceIndices[nFaces++] = treeBoundBox::BACK;
882 }
883 else if (ptFaceID & treeBoundBox::FRONTBIT)
884 {
885 faceIndices[nFaces++] = treeBoundBox::FRONT;
886 }
887
888
889 // Determine the face we want to keep the point on
890
891 direction keepFaceID;
892
893 if (nFaces == 0)
894 {
895 // Return original point
896 return pt;
897 }
898 else if (nFaces == 1)
899 {
900 // Point is on a single face
901 keepFaceID = faceIndices[0];
902 }
903 else
904 {
905 // Determine best face out of faceIndices[0 .. nFaces-1].
906 // The best face is the one most perpendicular to the ray direction.
907
908 keepFaceID = faceIndices[0];
909 scalar maxInproduct = mag(treeBoundBox::faceNormals[keepFaceID] & dir);
910
911 for (direction i = 1; i < nFaces; i++)
912 {
913 direction face = faceIndices[i];
914 scalar s = mag(treeBoundBox::faceNormals[face] & dir);
915 if (s > maxInproduct)
916 {
917 maxInproduct = s;
918 keepFaceID = face;
919 }
920 }
921 }
922
923
924 // 1. Push point into bb, away from all corners
925
926 point facePoint(pushPoint(bb, pt, true));
927 direction faceID = 0;
928
929 // 2. Snap it back onto the preferred face
930
931 if (keepFaceID == treeBoundBox::LEFT)
932 {
933 facePoint.x() = bb.min().x();
934 faceID = treeBoundBox::LEFTBIT;
935 }
936 else if (keepFaceID == treeBoundBox::RIGHT)
937 {
938 facePoint.x() = bb.max().x();
939 faceID = treeBoundBox::RIGHTBIT;
940 }
941 else if (keepFaceID == treeBoundBox::BOTTOM)
942 {
943 facePoint.y() = bb.min().y();
944 faceID = treeBoundBox::BOTTOMBIT;
945 }
946 else if (keepFaceID == treeBoundBox::TOP)
947 {
948 facePoint.y() = bb.max().y();
949 faceID = treeBoundBox::TOPBIT;
950 }
951 else if (keepFaceID == treeBoundBox::BACK)
952 {
953 facePoint.z() = bb.min().z();
954 faceID = treeBoundBox::BACKBIT;
955 }
956 else if (keepFaceID == treeBoundBox::FRONT)
957 {
958 facePoint.z() = bb.max().z();
959 faceID = treeBoundBox::FRONTBIT;
960 }
961
962
963 if (debug)
964 {
965 if (faceID != bb.faceBits(facePoint))
966 {
967 auto fatal = FatalErrorInFunction;
968
969 fatal
970 << "Pushed point from " << pt
971 << " on face:" << ptFaceID << " of bb:" << bb << nl
972 << "onto " << facePoint
973 << " on face:" << faceID
974 << " which is not consistent with geometric face "
975 << bb.faceBits(facePoint) << nl;
976
977 if (debug > 1)
978 {
979 fatal << abort(FatalError);
980 }
981 }
982 if (bb.posBits(facePoint) != 0)
983 {
984 auto fatal = FatalErrorInFunction;
985
986 fatal
987 << " bb:" << bb << nl
988 << "does not contain perturbed point "
989 << facePoint << nl;
990
991 if (debug > 1)
992 {
993 fatal << abort(FatalError);
994 }
995 }
996 }
997
998
999 return facePoint;
1000}
1001
1002
1003template<class Type>
1005(
1006 const label nodeI,
1007 const direction octant,
1008
1009 label& parentNodeI,
1010 label& parentOctant
1011) const
1012{
1013 parentNodeI = nodes_[nodeI].parent_;
1014
1015 if (parentNodeI == -1)
1016 {
1017 // Reached edge of tree
1018 return false;
1019 }
1020
1021 const node& parentNode = nodes_[parentNodeI];
1022
1023 // Find octant nodeI is in.
1024 parentOctant = 255;
1025
1026 for (direction i = 0; i < parentNode.subNodes_.size(); i++)
1027 {
1028 labelBits index = parentNode.subNodes_[i];
1029
1030 if (isNode(index) && getNode(index) == nodeI)
1031 {
1032 parentOctant = i;
1033 break;
1034 }
1035 }
1036
1037 if (parentOctant == 255)
1038 {
1040 << "Problem: no parent found for octant:" << octant
1041 << " node:" << nodeI
1042 << abort(FatalError);
1043 }
1044
1045 return true;
1046}
1047
1048
1049
1050template<class Type>
1052(
1053 const point& facePoint,
1054 const direction faceID, // face(s) that facePoint is on
1055 label& nodeI,
1056 direction& octant
1057) const
1058{
1059 // Walk tree to neighbouring node. Gets current position as node and octant
1060 // in this node and walks in the direction given by the facePointBits
1061 // (combination of treeBoundBox::LEFTBIT, TOPBIT etc.) Returns false if
1062 // edge of tree hit.
1063
1064 label oldNodeI = nodeI;
1065 direction oldOctant = octant;
1066
1067 // Find out how to test for octant. Say if we want to go left we need
1068 // to traverse up the tree until we hit a node where our octant is
1069 // on the right.
1070
1071 // Coordinate direction to test
1072 const direction X = treeBoundBox::RIGHTHALF;
1073 const direction Y = treeBoundBox::TOPHALF;
1074 const direction Z = treeBoundBox::FRONTHALF;
1075
1076 direction octantMask = 0;
1077 direction wantedValue = 0;
1078
1079 if ((faceID & treeBoundBox::LEFTBIT) != 0)
1080 {
1081 // We want to go left so check if in right octant (i.e. x-bit is set)
1082 octantMask |= X;
1083 wantedValue |= X;
1084 }
1085 else if ((faceID & treeBoundBox::RIGHTBIT) != 0)
1086 {
1087 octantMask |= X; // wantedValue already 0
1088 }
1089
1090 if ((faceID & treeBoundBox::BOTTOMBIT) != 0)
1091 {
1092 // Want to go down so check for y-bit set.
1093 octantMask |= Y;
1094 wantedValue |= Y;
1095 }
1096 else if ((faceID & treeBoundBox::TOPBIT) != 0)
1097 {
1098 // Want to go up so check for y-bit not set.
1099 octantMask |= Y;
1100 }
1101
1102 if ((faceID & treeBoundBox::BACKBIT) != 0)
1103 {
1104 octantMask |= Z;
1105 wantedValue |= Z;
1106 }
1107 else if ((faceID & treeBoundBox::FRONTBIT) != 0)
1108 {
1109 octantMask |= Z;
1110 }
1111
1112 // So now we have the coordinate directions in the octant we need to check
1113 // and the resulting value.
1114
1115 /*
1116 // +---+---+
1117 // | | |
1118 // | | |
1119 // | | |
1120 // +---+-+-+
1121 // | | | |
1122 // | a+-+-+
1123 // | |\| |
1124 // +---+-+-+
1125 // \
1126 //
1127 // e.g. ray is at (a) in octant 0(or 4) with faceIDs : LEFTBIT+TOPBIT.
1128 // If we would be in octant 1(or 5) we could go to the correct octant
1129 // in the same node by just flipping the x and y bits (exoring).
1130 // But if we are not in octant 1/5 we have to go up until we are.
1131 // In general for leftbit+topbit:
1132 // - we have to check for x and y : octantMask = 011
1133 // - the checked bits have to be : wantedValue = ?01
1134 */
1135
1136 //Pout<< "For point " << facePoint << endl;
1137
1138 // Go up until we have chance to cross to the wanted direction
1139 while (wantedValue != (octant & octantMask))
1140 {
1141 // Go up to the parent.
1142
1143 // Remove the directions that are not on the boundary of the parent.
1144 // See diagram above
1145 if (wantedValue & X) // && octantMask&X
1146 {
1147 // Looking for right octant.
1148 if (octant & X)
1149 {
1150 // My octant is one of the left ones so punch out the x check
1151 octantMask &= ~X;
1152 wantedValue &= ~X;
1153 }
1154 }
1155 else
1156 {
1157 if (!(octant & X))
1158 {
1159 // My octant is right but I want to go left.
1160 octantMask &= ~X;
1161 wantedValue &= ~X;
1162 }
1163 }
1164
1165 if (wantedValue & Y)
1166 {
1167 if (octant & Y)
1168 {
1169 octantMask &= ~Y;
1170 wantedValue &= ~Y;
1171 }
1172 }
1173 else
1174 {
1175 if (!(octant & Y))
1176 {
1177 octantMask &= ~Y;
1178 wantedValue &= ~Y;
1179 }
1180 }
1181
1182 if (wantedValue & Z)
1183 {
1184 if (octant & Z)
1185 {
1186 octantMask &= ~Z;
1187 wantedValue &= ~Z;
1188 }
1189 }
1190 else
1191 {
1192 if (!(octant & Z))
1193 {
1194 octantMask &= ~Z;
1195 wantedValue &= ~Z;
1196 }
1197 }
1198
1199
1200 label parentNodeI;
1201 label parentOctant;
1202 walkToParent(nodeI, octant, parentNodeI, parentOctant);
1203
1204 if (parentNodeI == -1)
1205 {
1206 // Reached edge of tree
1207 return false;
1208 }
1209
1210 //Pout<< " walked from node:" << nodeI << " octant:" << octant
1211 // << " bb:" << nodes_[nodeI].bb_.subBbox(octant) << endl
1212 // << " to:" << parentNodeI << " octant:" << parentOctant
1213 // << " bb:" << nodes_[parentNodeI].bb_.subBbox(parentOctant)
1214 // << endl;
1215 //
1216 //Pout<< " octantMask:" << octantMask
1217 // << " wantedValue:" << wantedValue << endl;
1218
1219 nodeI = parentNodeI;
1220 octant = parentOctant;
1221 }
1222
1223 // So now we hit the correct parent node. Switch to the correct octant.
1224 // Is done by jumping to the 'other half' so if e.g. in x direction in
1225 // right half we now jump to the left half.
1226 octant ^= octantMask;
1227
1228 //Pout<< " to node:" << nodeI << " octant:" << octant
1229 // << " subBb:" <<subBbox(nodeI, octant) << endl;
1230
1231
1232 if (debug)
1233 {
1234 const treeBoundBox subBb(subBbox(nodeI, octant));
1235
1236 if (!subBb.contains(facePoint))
1237 {
1238 auto fatal = FatalErrorInFunction;
1239
1240 fatal
1241 << "When searching for " << facePoint
1242 << " ended up in node:" << nodeI
1243 << " octant:" << octant
1244 << " with bb:" << subBb << nl;
1245
1246 if (debug > 1)
1247 {
1248 fatal << abort(FatalError);
1249 }
1250 }
1251 }
1252
1253
1254 // See if we need to travel down. Note that we already go into the
1255 // the first level ourselves (instead of having findNode decide)
1256 labelBits index = nodes_[nodeI].subNodes_[octant];
1257
1258 if (isNode(index))
1259 {
1260 labelBits node = findNode(getNode(index), facePoint);
1261
1262 nodeI = getNode(node);
1263 octant = getOctant(node);
1264 }
1265
1266
1267 if (debug)
1268 {
1269 const treeBoundBox subBb(subBbox(nodeI, octant));
1270
1271 if (nodeI == oldNodeI && octant == oldOctant)
1272 {
1273 auto fatal = FatalErrorInFunction;
1274
1275 fatal
1276 << "Did not go to neighbour when searching for " << facePoint
1277 << nl
1278 << " starting from face:" << faceString(faceID)
1279 << " node:" << nodeI
1280 << " octant:" << octant
1281 << " bb:" << subBb << nl;
1282
1283 if (debug > 1)
1284 {
1285 fatal << abort(FatalError);
1286 }
1287 }
1288
1289 if (!subBb.contains(facePoint))
1290 {
1291 auto fatal = FatalErrorInFunction;
1292
1293 fatal
1294 << "When searching for " << facePoint
1295 << " ended up in node:" << nodeI
1296 << " octant:" << octant
1297 << " bb:" << subBb << nl;
1298
1299 if (debug > 1)
1300 {
1301 fatal << abort(FatalError);
1302 }
1303 }
1304 }
1305
1306
1307 return true;
1308}
1309
1310
1311template<class Type>
1313(
1314 const direction faceID
1315)
1316{
1317 word desc;
1318
1319 if (faceID == 0)
1320 {
1321 desc = "noFace";
1322 }
1323 if (faceID & treeBoundBox::LEFTBIT)
1324 {
1325 if (!desc.empty()) desc += "+";
1326 desc += "left";
1327 }
1328 if (faceID & treeBoundBox::RIGHTBIT)
1329 {
1330 if (!desc.empty()) desc += "+";
1331 desc += "right";
1332 }
1333 if (faceID & treeBoundBox::BOTTOMBIT)
1334 {
1335 if (!desc.empty()) desc += "+";
1336 desc += "bottom";
1337 }
1338 if (faceID & treeBoundBox::TOPBIT)
1339 {
1340 if (!desc.empty()) desc += "+";
1341 desc += "top";
1342 }
1343 if (faceID & treeBoundBox::BACKBIT)
1344 {
1345 if (!desc.empty()) desc += "+";
1346 desc += "back";
1347 }
1348 if (faceID & treeBoundBox::FRONTBIT)
1349 {
1350 if (!desc.empty()) desc += "+";
1351 desc += "front";
1352 }
1353 return desc;
1354}
1355
1356
1357template<class Type>
1359(
1360 const bool findAny,
1361 const point& treeStart,
1362 const vector& treeVec,
1363
1364 const point& start,
1365 const point& end,
1366 const label nodeI,
1367 const direction octant,
1368
1369 pointIndexHit& hitInfo,
1370 direction& hitBits
1371) const
1372{
1373 if (debug)
1374 {
1375 const treeBoundBox octantBb(subBbox(nodeI, octant));
1376
1377 if (octantBb.posBits(start) != 0)
1378 {
1379 auto fatal = FatalErrorInFunction;
1380
1381 fatal
1382 << "Node:" << nodeI << " octant:" << octant
1383 << " bb:" << octantBb << nl
1384 << "does not contain point " << start << nl;
1385
1386 if (debug > 1)
1387 {
1388 fatal << abort(FatalError);
1389 }
1390 }
1391 }
1392
1393
1394 const node& nod = nodes_[nodeI];
1395
1396 labelBits index = nod.subNodes_[octant];
1397
1398 if (isContent(index))
1399 {
1400 const labelList& indices = *(contents_[getContent(index)]);
1401
1402 if (indices.size())
1403 {
1404 if (findAny)
1405 {
1406 // Find any intersection
1407
1408 forAll(indices, elemI)
1409 {
1410 label shapeI = indices[elemI];
1411
1412 point pt;
1413 bool hit = shapes_.intersects(shapeI, start, end, pt);
1414
1415 // Note that intersection of shape might actually be
1416 // in a neighbouring box. For findAny this is not important.
1417 if (hit)
1418 {
1419 // Hit so pt is nearer than nearestPoint.
1420 // Update hit info
1421 hitInfo.setHit();
1422 hitInfo.setIndex(shapeI);
1423 hitInfo.setPoint(pt);
1424 return;
1425 }
1426 }
1427 }
1428 else
1429 {
1430 // Find nearest intersection
1431
1432 const treeBoundBox octantBb(subBbox(nodeI, octant));
1433
1434 point nearestPoint(end);
1435
1436 forAll(indices, elemI)
1437 {
1438 label shapeI = indices[elemI];
1439
1440 point pt;
1441 bool hit = shapes_.intersects
1442 (
1443 shapeI,
1444 start,
1445 nearestPoint,
1446 pt
1447 );
1448
1449 // Note that intersection of shape might actually be
1450 // in a neighbouring box. Since we need to maintain strict
1451 // (findAny=false) ordering skip such an intersection. It
1452 // will be found when we are doing the next box.
1453
1454 if (hit && octantBb.contains(pt))
1455 {
1456 // Hit so pt is nearer than nearestPoint.
1457 nearestPoint = pt;
1458 // Update hit info
1459 hitInfo.setHit();
1460 hitInfo.setIndex(shapeI);
1461 hitInfo.setPoint(pt);
1462 }
1463 }
1464
1465 if (hitInfo.hit())
1466 {
1467 // Found intersected shape.
1468 return;
1469 }
1470 }
1471 }
1472 }
1473
1474 // Nothing intersected in this node
1475 // Traverse to end of node. Do by ray tracing back from end and finding
1476 // intersection with bounding box of node.
1477 // start point is now considered to be inside bounding box.
1478
1479 const treeBoundBox octantBb(subBbox(nodeI, octant));
1480
1481 point pt;
1482 bool intersected = octantBb.intersects
1483 (
1484 end, //treeStart,
1485 (start-end), //treeVec,
1486
1487 end,
1488 start,
1489
1490 pt,
1491 hitBits
1492 );
1493
1494
1495 if (intersected)
1496 {
1497 // Return miss. Set misspoint to face.
1498 hitInfo.setPoint(pt);
1499 }
1500 else
1501 {
1502 // Rare case: did not find intersection of ray with octantBb. Can
1503 // happen if end is on face/edge of octantBb. Do like
1504 // lagrangian and re-track with perturbed vector (slightly pushed out
1505 // of bounding box)
1506
1507 point perturbedEnd(pushPoint(octantBb, end, false));
1508
1509 traverseNode
1510 (
1511 findAny,
1512 treeStart,
1513 treeVec,
1514 start,
1515 perturbedEnd,
1516 nodeI,
1517 octant,
1518
1519 hitInfo,
1520 hitBits
1521 );
1522 }
1523}
1524
1525
1526template<class Type>
1528(
1529 const bool findAny,
1530 const point& treeStart,
1531 const point& treeEnd,
1532 const label startNodeI,
1533 const direction startOctant,
1534 const bool verbose
1535) const
1536{
1537 const vector treeVec(treeEnd - treeStart);
1538
1539 // Current node as parent+octant
1540 label nodeI = startNodeI;
1541 direction octant = startOctant;
1542
1543 if (verbose)
1544 {
1545 Pout<< "findLine : treeStart:" << treeStart
1546 << " treeEnd:" << treeEnd << endl
1547 << "node:" << nodeI
1548 << " octant:" << octant
1549 << " bb:" << subBbox(nodeI, octant) << endl;
1550 }
1551
1552 // Current position. Initialize to miss
1553 pointIndexHit hitInfo(false, treeStart, -1);
1554
1555 //while (true)
1556 label i = 0;
1557 for (; i < 100000; i++)
1558 {
1559 // Ray-trace to end of current node. Updates point (either on triangle
1560 // in case of hit or on node bounding box in case of miss)
1561
1562 const treeBoundBox octantBb(subBbox(nodeI, octant));
1563
1564 // Make sure point is away from any edges/corners
1565 point startPoint
1566 (
1567 pushPointIntoFace
1568 (
1569 octantBb,
1570 treeVec,
1571 hitInfo.rawPoint()
1572 )
1573 );
1574
1575 if (verbose)
1576 {
1577 Pout<< "iter:" << i
1578 << " at current:" << hitInfo.rawPoint()
1579 << " (perturbed:" << startPoint << ")" << endl
1580 << " node:" << nodeI
1581 << " octant:" << octant
1582 << " bb:" << subBbox(nodeI, octant) << endl;
1583 }
1584
1585
1586
1587
1588 // Faces of current bounding box current point is on
1589 direction hitFaceID = 0;
1590
1591 traverseNode
1592 (
1593 findAny,
1594 treeStart,
1595 treeVec,
1596
1597 startPoint, // Note: pass in copy since hitInfo
1598 // also used in return value.
1599 treeEnd, // pass in overall end so is nicely outside bb
1600 nodeI,
1601 octant,
1602
1603 hitInfo,
1604 hitFaceID
1605 );
1606
1607 // Did we hit a triangle?
1608 if (hitInfo.hit())
1609 {
1610 break;
1611 }
1612
1613 if (hitFaceID == 0 || hitInfo.rawPoint() == treeEnd)
1614 {
1615 // endpoint inside the tree. Return miss.
1616 break;
1617 }
1618
1619 // Create a point on other side of face.
1620 point perturbedPoint
1621 (
1622 pushPoint
1623 (
1624 octantBb,
1625 hitFaceID,
1626 hitInfo.rawPoint(),
1627 false // push outside of octantBb
1628 )
1629 );
1630
1631 if (verbose)
1632 {
1633 Pout<< " iter:" << i
1634 << " hit face:" << faceString(hitFaceID)
1635 << " at:" << hitInfo.rawPoint() << nl
1636 << " node:" << nodeI
1637 << " octant:" << octant
1638 << " bb:" << subBbox(nodeI, octant) << nl
1639 << " walking to neighbour containing:" << perturbedPoint
1640 << endl;
1641 }
1642
1643
1644 // Nothing hit so we are on face of bounding box (given as node+octant+
1645 // position bits). Traverse to neighbouring node. Use slightly perturbed
1646 // point.
1647
1648 bool ok = walkToNeighbour
1649 (
1650 perturbedPoint,
1651 hitFaceID, // face(s) that hitInfo is on
1652
1653 nodeI,
1654 octant
1655 );
1656
1657 if (!ok)
1658 {
1659 // Hit the edge of the tree. Return miss.
1660 break;
1661 }
1662
1663 if (verbose)
1664 {
1665 const treeBoundBox octantBb(subBbox(nodeI, octant));
1666 Pout<< " walked for point:" << hitInfo.rawPoint() << endl
1667 << " to neighbour node:" << nodeI
1668 << " octant:" << octant
1669 << " face:" << faceString(octantBb.faceBits(hitInfo.rawPoint()))
1670 << " of octantBb:" << octantBb << endl
1671 << endl;
1672 }
1673 }
1674
1675 if (i == 100000)
1676 {
1677 // Probably in loop.
1678 if (!verbose)
1679 {
1680 // Redo intersection but now with verbose flag switched on.
1681 return findLine
1682 (
1683 findAny,
1684 treeStart,
1685 treeEnd,
1686 startNodeI,
1687 startOctant,
1688 true //verbose
1689 );
1690 }
1691 if (debug)
1692 {
1694 << "Got stuck in loop raytracing from:" << treeStart
1695 << " to:" << treeEnd << endl
1696 << "inside top box:" << subBbox(startNodeI, startOctant)
1697 << abort(FatalError);
1698 }
1699 else
1700 {
1702 << "Got stuck in loop raytracing from:" << treeStart
1703 << " to:" << treeEnd << endl
1704 << "inside top box:" << subBbox(startNodeI, startOctant)
1705 << endl;
1706 }
1707 }
1708
1709 return hitInfo;
1710}
1711
1712
1713template<class Type>
1715(
1716 const bool findAny,
1717 const point& start,
1718 const point& end
1719) const
1720{
1721 pointIndexHit hitInfo;
1722
1723 if (nodes_.size())
1724 {
1725 const treeBoundBox& treeBb = nodes_[0].bb_;
1726
1727 // No effort is made to deal with points which are on edge of tree
1728 // bounding box for now.
1729
1730 direction startBit = treeBb.posBits(start);
1731 direction endBit = treeBb.posBits(end);
1732
1733 if ((startBit & endBit) != 0)
1734 {
1735 // Both start and end outside domain and in same block.
1736 return pointIndexHit(false, Zero, -1);
1737 }
1738
1739
1740 // Trim segment to treeBb
1741
1742 point trackStart(start);
1743 point trackEnd(end);
1744
1745 if (startBit != 0)
1746 {
1747 // Track start to inside domain.
1748 if (!treeBb.intersects(start, end, trackStart))
1749 {
1750 return pointIndexHit(false, Zero, -1);
1751 }
1752 }
1753
1754 if (endBit != 0)
1755 {
1756 // Track end to inside domain.
1757 if (!treeBb.intersects(end, trackStart, trackEnd))
1758 {
1759 return pointIndexHit(false, Zero, -1);
1760 }
1761 }
1762
1763
1764 // Find lowest level tree node that start is in.
1765 labelBits index = findNode(0, trackStart);
1766
1767 label parentNodeI = getNode(index);
1768 direction octant = getOctant(index);
1769
1770 hitInfo = findLine
1771 (
1772 findAny,
1773 trackStart,
1774 trackEnd,
1775 parentNodeI,
1776 octant
1777 );
1778 }
1779
1780 return hitInfo;
1781}
1782
1783
1784template<class Type>
1786(
1787 const label nodeI,
1788 const treeBoundBox& searchBox,
1789 labelHashSet& elements
1790) const
1791{
1792 const node& nod = nodes_[nodeI];
1793 const treeBoundBox& nodeBb = nod.bb_;
1794
1795 for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
1796 {
1797 labelBits index = nod.subNodes_[octant];
1798
1799 if (isNode(index))
1800 {
1801 const treeBoundBox& subBb = nodes_[getNode(index)].bb_;
1802
1803 if (subBb.overlaps(searchBox))
1804 {
1805 findBox(getNode(index), searchBox, elements);
1806 }
1807 }
1808 else if (isContent(index))
1809 {
1810 const treeBoundBox subBb(nodeBb.subBbox(octant));
1811
1812 if (subBb.overlaps(searchBox))
1813 {
1814 const labelList& indices = *(contents_[getContent(index)]);
1815
1816 forAll(indices, i)
1817 {
1818 label shapeI = indices[i];
1819
1820 if (shapes_.overlaps(shapeI, searchBox))
1821 {
1822 elements.insert(shapeI);
1823 }
1824 }
1825 }
1826 }
1827 }
1828}
1829
1830
1831template<class Type>
1833(
1834 const label nodeI,
1835 const point& centre,
1836 const scalar radiusSqr,
1837 labelHashSet& elements
1838) const
1839{
1840 const node& nod = nodes_[nodeI];
1841 const treeBoundBox& nodeBb = nod.bb_;
1842
1843 for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
1844 {
1845 labelBits index = nod.subNodes_[octant];
1846
1847 if (isNode(index))
1848 {
1849 const treeBoundBox& subBb = nodes_[getNode(index)].bb_;
1850
1851 if (subBb.overlaps(centre, radiusSqr))
1852 {
1853 findSphere(getNode(index), centre, radiusSqr, elements);
1854 }
1855 }
1856 else if (isContent(index))
1857 {
1858 const treeBoundBox subBb(nodeBb.subBbox(octant));
1859
1860 if (subBb.overlaps(centre, radiusSqr))
1861 {
1862 const labelList& indices = *(contents_[getContent(index)]);
1863
1864 forAll(indices, i)
1865 {
1866 label shapeI = indices[i];
1867
1868 if (shapes_.overlaps(shapeI, centre, radiusSqr))
1869 {
1870 elements.insert(shapeI);
1871 }
1872 }
1873 }
1874 }
1875 }
1876}
1877
1878
1879template<class Type>
1880template<class CompareOp>
1882(
1883 const scalar nearDist,
1884 const bool okOrder,
1885 const dynamicIndexedOctree<Type>& tree1,
1886 const labelBits index1,
1887 const treeBoundBox& bb1,
1888 const dynamicIndexedOctree<Type>& tree2,
1889 const labelBits index2,
1890 const treeBoundBox& bb2,
1891 CompareOp& cop
1892)
1893{
1894 const vector nearSpan(nearDist, nearDist, nearDist);
1895
1896 if (tree1.isNode(index1))
1897 {
1898 const node& nod1 = tree1.nodes()[tree1.getNode(index1)];
1899 const treeBoundBox searchBox
1900 (
1901 bb1.min()-nearSpan,
1902 bb1.max()+nearSpan
1903 );
1904
1905 if (tree2.isNode(index2))
1906 {
1907 if (bb2.overlaps(searchBox))
1908 {
1909 const node& nod2 = tree2.nodes()[tree2.getNode(index2)];
1910
1911 for (direction i2 = 0; i2 < nod2.subNodes_.size(); i2++)
1912 {
1913 labelBits subIndex2 = nod2.subNodes_[i2];
1914 const treeBoundBox subBb2
1915 (
1916 tree2.isNode(subIndex2)
1917 ? tree2.nodes()[tree2.getNode(subIndex2)].bb_
1918 : bb2.subBbox(i2)
1919 );
1920
1921 findNear
1922 (
1923 nearDist,
1924 !okOrder,
1925 tree2,
1926 subIndex2,
1927 subBb2,
1928 tree1,
1929 index1,
1930 bb1,
1931 cop
1932 );
1933 }
1934 }
1935 }
1936 else if (tree2.isContent(index2))
1937 {
1938 // index2 is leaf, index1 not yet.
1939 for (direction i1 = 0; i1 < nod1.subNodes_.size(); i1++)
1940 {
1941 labelBits subIndex1 = nod1.subNodes_[i1];
1942 const treeBoundBox subBb1
1943 (
1944 tree1.isNode(subIndex1)
1945 ? tree1.nodes()[tree1.getNode(subIndex1)].bb_
1946 : bb1.subBbox(i1)
1947 );
1948
1949 findNear
1950 (
1951 nearDist,
1952 !okOrder,
1953 tree2,
1954 index2,
1955 bb2,
1956 tree1,
1957 subIndex1,
1958 subBb1,
1959 cop
1960 );
1961 }
1962 }
1963 }
1964 else if (tree1.isContent(index1))
1965 {
1966 const treeBoundBox searchBox
1967 (
1968 bb1.min()-nearSpan,
1969 bb1.max()+nearSpan
1970 );
1971
1972 if (tree2.isNode(index2))
1973 {
1974 const node& nod2 =
1975 tree2.nodes()[tree2.getNode(index2)];
1976
1977 if (bb2.overlaps(searchBox))
1978 {
1979 for (direction i2 = 0; i2 < nod2.subNodes_.size(); i2++)
1980 {
1981 labelBits subIndex2 = nod2.subNodes_[i2];
1982 const treeBoundBox subBb2
1983 (
1984 tree2.isNode(subIndex2)
1985 ? tree2.nodes()[tree2.getNode(subIndex2)].bb_
1986 : bb2.subBbox(i2)
1987 );
1988
1989 findNear
1990 (
1991 nearDist,
1992 !okOrder,
1993 tree2,
1994 subIndex2,
1995 subBb2,
1996 tree1,
1997 index1,
1998 bb1,
1999 cop
2000 );
2001 }
2002 }
2003 }
2004 else if (tree2.isContent(index2))
2005 {
2006 // Both are leaves. Check n^2.
2007
2008 const labelList& indices1 =
2009 tree1.contents()[tree1.getContent(index1)];
2010 const labelList& indices2 =
2011 tree2.contents()[tree2.getContent(index2)];
2012
2013 forAll(indices1, i)
2014 {
2015 label shape1 = indices1[i];
2016
2017 forAll(indices2, j)
2018 {
2019 label shape2 = indices2[j];
2020
2021 if ((&tree1 != &tree2) || (shape1 != shape2))
2022 {
2023 if (okOrder)
2024 {
2025 cop
2026 (
2027 nearDist,
2028 tree1.shapes(),
2029 shape1,
2030 tree2.shapes(),
2031 shape2
2032 );
2033 }
2034 else
2035 {
2036 cop
2037 (
2038 nearDist,
2039 tree2.shapes(),
2040 shape2,
2041 tree1.shapes(),
2042 shape1
2043 );
2044 }
2045 }
2046 }
2047 }
2048 }
2049 }
2050}
2051
2052
2053template<class Type>
2055(
2056 const labelBits index
2057) const
2058{
2059 label nElems = 0;
2060
2061 if (isNode(index))
2062 {
2063 // tree node.
2064 label nodeI = getNode(index);
2065
2066 const node& nod = nodes_[nodeI];
2067
2068 for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
2069 {
2070 nElems += countElements(nod.subNodes_[octant]);
2071 }
2072 }
2073 else if (isContent(index))
2074 {
2075 nElems += contents_[getContent(index)]->size();
2076 }
2077 else
2078 {
2079 // empty node
2080 }
2081
2082 return nElems;
2083}
2084
2085
2086template<class Type>
2088(
2089 const label nodeI,
2090 const direction octant
2091) const
2092{
2093 OFstream str
2094 (
2095 "node" + Foam::name(nodeI) + "_octant" + Foam::name(octant) + ".obj"
2096 );
2097
2098 labelBits index = nodes_[nodeI].subNodes_[octant];
2099
2100 treeBoundBox subBb;
2101
2102 if (isNode(index))
2103 {
2104 subBb = nodes_[getNode(index)].bb_;
2105 }
2106 else if (isContent(index) || isEmpty(index))
2107 {
2108 subBb = nodes_[nodeI].bb_.subBbox(octant);
2109 }
2110
2111 Pout<< "dumpContentNode : writing node:" << nodeI << " octant:" << octant
2112 << " to " << str.name() << endl;
2113
2114 // Dump bounding box
2115 pointField bbPoints(subBb.points());
2116
2117 forAll(bbPoints, i)
2118 {
2119 const point& pt = bbPoints[i];
2120
2121 str<< "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
2122 }
2123
2124 forAll(treeBoundBox::edges, i)
2125 {
2126 const edge& e = treeBoundBox::edges[i];
2127
2128 str<< "l " << e[0] + 1 << ' ' << e[1] + 1 << nl;
2129 }
2130}
2131
2132
2133// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
2134
2135template<class Type>
2137(
2138 const Type& shapes,
2139 const treeBoundBox& bb,
2140 const label maxLevels, // maximum number of levels
2141 const scalar maxLeafRatio,
2142 const scalar maxDuplicity
2143)
2144:
2145 shapes_(shapes),
2146 bb_(bb),
2147 maxLevels_(maxLevels),
2148 nLevelsMax_(0),
2149 maxLeafRatio_(maxLeafRatio),
2150 minSize_(label(maxLeafRatio)),
2151 maxDuplicity_(maxDuplicity),
2152 nodes_(label(shapes.size() / maxLeafRatio_)),
2153 contents_(label(shapes.size() / maxLeafRatio_)),
2154 nodeTypes_(0)
2155{
2156 if (shapes_.size() == 0)
2157 {
2158 return;
2159 }
2160
2161 insert(0, shapes_.size());
2162
2163 if (debug)
2164 {
2165 writeTreeInfo();
2166 }
2167}
2168
2169
2170// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
2171
2172template<class Type>
2174{
2175 return perturbTol_;
2176}
2177
2178
2179template<class Type>
2181(
2182 const point& sample,
2183 const scalar startDistSqr
2184) const
2185{
2186 scalar nearestDistSqr = startDistSqr;
2187 label nearestShapeI = -1;
2188 point nearestPoint = Zero;
2189
2190 if (nodes_.size())
2191 {
2192 findNearest
2193 (
2194 0,
2195 sample,
2196
2197 nearestDistSqr,
2198 nearestShapeI,
2199 nearestPoint
2200 );
2201 }
2202
2203 return pointIndexHit(nearestShapeI != -1, nearestPoint, nearestShapeI);
2204}
2205
2206
2207template<class Type>
2209(
2210 const linePointRef& ln,
2211 treeBoundBox& tightest,
2212 point& linePoint
2213) const
2214{
2215 label nearestShapeI = -1;
2216 point nearestPoint;
2217
2218 if (nodes_.size())
2219 {
2220 findNearest
2221 (
2222 0,
2223 ln,
2224
2225 tightest,
2226 nearestShapeI,
2227 linePoint,
2228 nearestPoint
2229 );
2230 }
2231 else
2232 {
2233 nearestPoint = Zero;
2234 }
2235
2236 return pointIndexHit(nearestShapeI != -1, nearestPoint, nearestShapeI);
2237}
2238
2239
2240template<class Type>
2242(
2243 const point& start,
2244 const point& end
2245) const
2246{
2247 return findLine(false, start, end);
2248}
2249
2250
2251template<class Type>
2253(
2254 const point& start,
2255 const point& end
2256) const
2257{
2258 return findLine(true, start, end);
2259}
2260
2261
2262template<class Type>
2264(
2265 const treeBoundBox& searchBox
2266) const
2267{
2268 if (nodes_.empty())
2269 {
2270 return labelList();
2271 }
2272
2273 // Storage for labels of shapes inside bb. Size estimate.
2274 labelHashSet elements(shapes_.size() / 100);
2275
2276 findBox(0, searchBox, elements);
2277
2278 return elements.toc();
2279}
2280
2281
2282template<class Type>
2284(
2285 const point& centre,
2286 const scalar radiusSqr
2287) const
2288{
2289 if (nodes_.empty())
2290 {
2291 return labelList();
2292 }
2293
2294 // Storage for labels of shapes inside bb. Size estimate.
2295 labelHashSet elements(shapes_.size() / 100);
2296
2297 findSphere(0, centre, radiusSqr, elements);
2298
2299 return elements.toc();
2300}
2301
2302
2303template<class Type>
2305(
2306 const label nodeI,
2307 const point& sample
2308) const
2309{
2310 if (nodes_.empty())
2311 {
2312 // Empty tree. Return what?
2313 return nodePlusOctant(nodeI, 0);
2314 }
2315
2316 const node& nod = nodes_[nodeI];
2317
2318 if (debug)
2319 {
2320 if (!nod.bb_.contains(sample))
2321 {
2323 << "Cannot find " << sample << " in node " << nodeI
2324 << abort(FatalError);
2325 }
2326 }
2327
2328 direction octant = nod.bb_.subOctant(sample);
2329
2330 labelBits index = nod.subNodes_[octant];
2331
2332 if (isNode(index))
2333 {
2334 // Recurse
2335 return findNode(getNode(index), sample);
2336 }
2337 else if (isContent(index))
2338 {
2339 // Content. Return treenode+octant
2340 return nodePlusOctant(nodeI, octant);
2341 }
2342 else
2343 {
2344 // Empty. Return treenode+octant
2345 return nodePlusOctant(nodeI, octant);
2346 }
2347}
2348
2349
2350template<class Type>
2352(
2353 const point& sample
2354) const
2355{
2356 labelBits index = findNode(0, sample);
2357
2358 const node& nod = nodes_[getNode(index)];
2359
2360 labelBits contentIndex = nod.subNodes_[getOctant(index)];
2361
2362 // Need to check for the presence of content, in-case the node is empty
2363 if (isContent(contentIndex))
2364 {
2365 const labelList& indices = *(contents_[getContent(contentIndex)]);
2366
2367 forAll(indices, elemI)
2368 {
2369 label shapeI = indices[elemI];
2370
2371 if (shapes_.contains(shapeI, sample))
2372 {
2373 return shapeI;
2374 }
2375 }
2376 }
2377
2378 return -1;
2379}
2380
2381
2382template<class Type>
2384(
2385 const point& sample
2386) const
2387{
2388 labelBits index = findNode(0, sample);
2389
2390 const node& nod = nodes_[getNode(index)];
2391
2392 labelBits contentIndex = nod.subNodes_[getOctant(index)];
2393
2394 // Need to check for the presence of content, in-case the node is empty
2395 if (isContent(contentIndex))
2396 {
2397 return *(contents_[getContent(contentIndex)]);
2398 }
2399
2400 return labelList::null();
2401}
2402
2403
2404template<class Type>
2406(
2407 const point& sample
2408) const
2409{
2410 if (nodes_.empty())
2411 {
2412 return volumeType::UNKNOWN;
2413 }
2414
2415 if (nodeTypes_.size() != 8*nodes_.size())
2416 {
2417 // Calculate type for every octant of node.
2418
2419 nodeTypes_.setSize(8*nodes_.size());
2420 nodeTypes_ = volumeType::UNKNOWN;
2421
2422 calcVolumeType(0);
2423
2424 if (debug)
2425 {
2426 label nUNKNOWN = 0;
2427 label nMIXED = 0;
2428 label nINSIDE = 0;
2429 label nOUTSIDE = 0;
2430
2431 forAll(nodeTypes_, i)
2432 {
2433 volumeType type = volumeType::type(nodeTypes_.get(i));
2434
2436 {
2437 nUNKNOWN++;
2438 }
2439 else if (type == volumeType::MIXED)
2440 {
2441 nMIXED++;
2442 }
2443 else if (type == volumeType::INSIDE)
2444 {
2445 nINSIDE++;
2446 }
2447 else if (type == volumeType::OUTSIDE)
2448 {
2449 nOUTSIDE++;
2450 }
2451 else
2452 {
2454 }
2455 }
2456
2457 Pout<< "dynamicIndexedOctree<Type>::getVolumeType : "
2458 << " bb:" << bb()
2459 << " nodes_:" << nodes_.size()
2460 << " nodeTypes_:" << nodeTypes_.size()
2461 << " nUNKNOWN:" << nUNKNOWN
2462 << " nMIXED:" << nMIXED
2463 << " nINSIDE:" << nINSIDE
2464 << " nOUTSIDE:" << nOUTSIDE
2465 << endl;
2466 }
2467 }
2468
2469 return getVolumeType(0, sample);
2470}
2471
2472
2473template<class Type>
2474template<class CompareOp>
2476(
2477 const scalar nearDist,
2478 const dynamicIndexedOctree<Type>& tree2,
2479 CompareOp& cop
2480) const
2481{
2482 findNear
2483 (
2484 nearDist,
2485 true,
2486 *this,
2487 nodePlusOctant(0, 0),
2488 bb(),
2489 tree2,
2490 nodePlusOctant(0, 0),
2491 tree2.bb(),
2492 cop
2493 );
2494}
2495
2496
2497template<class Type>
2498bool Foam::dynamicIndexedOctree<Type>::insert(label startIndex, label endIndex)
2499{
2500 if (startIndex == endIndex)
2501 {
2502 return false;
2503 }
2504
2505 if (nodes_.empty())
2506 {
2507 contents_.append
2508 (
2510 (
2511 new DynamicList<label>(1)
2512 )
2513 );
2514
2515 contents_[0]->append(0);
2516
2517 // Create topnode.
2518 node topNode = divide(bb_, 0, -1, 0);
2519
2520 nodes_.append(topNode);
2521
2522 startIndex++;
2523 }
2524
2525 bool success = true;
2526
2527 for (label pI = startIndex; pI < endIndex; ++pI)
2528 {
2529 label nLevels = 1;
2530
2531 if (!insertIndex(0, pI, nLevels))
2532 {
2533 success = false;
2534 }
2535
2536 nLevelsMax_ = max(nLevels, nLevelsMax_);
2537 }
2538
2539 return success;
2540}
2541
2542
2543template<class Type>
2545(
2546 const label nodIndex,
2547 const label index,
2548 label& nLevels
2549)
2550{
2551 bool shapeInserted = false;
2552
2553 for (direction octant = 0; octant < 8; octant++)
2554 {
2555 const labelBits& subNodeLabel = nodes_[nodIndex].subNodes_[octant];
2556
2557 if (isNode(subNodeLabel))
2558 {
2559 const treeBoundBox& subBb = nodes_[getNode(subNodeLabel)].bb_;
2560
2561 if (shapes().overlaps(index, subBb))
2562 {
2563 nLevels++;
2564
2565 if (insertIndex(getNode(subNodeLabel), index, nLevels))
2566 {
2567 shapeInserted = true;
2568 }
2569 }
2570 }
2571 else if (isContent(subNodeLabel))
2572 {
2573 const treeBoundBox subBb = nodes_[nodIndex].bb_.subBbox(octant);
2574
2575 if (shapes().overlaps(index, subBb))
2576 {
2577 const label contentI = getContent(subNodeLabel);
2578
2579 contents_[contentI]->append(index);
2580
2581 recursiveSubDivision
2582 (
2583 subBb,
2584 contentI,
2585 nodIndex,
2586 octant,
2587 nLevels
2588 );
2589
2590 shapeInserted = true;
2591 }
2592 }
2593 else
2594 {
2595 const treeBoundBox subBb = nodes_[nodIndex].bb_.subBbox(octant);
2596
2597 if (shapes().overlaps(index, subBb))
2598 {
2599 label sz = contents_.size();
2600
2601 contents_.append
2602 (
2604 );
2605
2606 contents_[sz]->append(index);
2607
2608 nodes_[nodIndex].subNodes_[octant]
2609 = contentPlusOctant(sz, octant);
2610 }
2611
2612 shapeInserted = true;
2613 }
2614 }
2615
2616 return shapeInserted;
2617}
2618
2619
2620template<class Type>
2622{
2623 if (nodes_.empty())
2624 {
2625 return false;
2626 }
2627
2628 removeIndex(0, index);
2629
2630 return true;
2631}
2632
2633
2634template<class Type>
2636(
2637 const label nodIndex,
2638 const label index
2639)
2640{
2641 label totalContents = 0;
2642
2643 for (direction octant = 0; octant < 8; octant++)
2644 {
2645 const labelBits& subNodeLabel = nodes_[nodIndex].subNodes_[octant];
2646
2647 if (isNode(subNodeLabel))
2648 {
2649 const treeBoundBox& subBb = nodes_[getNode(subNodeLabel)].bb_;
2650
2651 if (shapes().overlaps(index, subBb))
2652 {
2653 // Recursive function.
2654 label childContentsSize
2655 = removeIndex(getNode(subNodeLabel), index);
2656
2657 totalContents += childContentsSize;
2658
2659 if (childContentsSize == 0)
2660 {
2661 nodes_[nodIndex].subNodes_[octant]
2662 = emptyPlusOctant(octant);
2663 }
2664 }
2665 else
2666 {
2667 // Increment this so that the node will not be marked as empty
2668 totalContents++;
2669 }
2670 }
2671 else if (isContent(subNodeLabel))
2672 {
2673 const treeBoundBox subBb = nodes_[nodIndex].bb_.subBbox(octant);
2674
2675 const label contentI = getContent(subNodeLabel);
2676
2677 if (shapes().overlaps(index, subBb))
2678 {
2679 DynamicList<label>& contentList = *(contents_[contentI]);
2680
2681 DynamicList<label> newContent(contentList.size());
2682
2683 forAll(contentList, pI)
2684 {
2685 const label oldIndex = contentList[pI];
2686
2687 if (oldIndex != index)
2688 {
2689 newContent.append(oldIndex);
2690 }
2691 }
2692
2693 newContent.shrink();
2694
2695 if (newContent.size() == 0)
2696 {
2697 // Set to empty.
2698 nodes_[nodIndex].subNodes_[octant]
2699 = emptyPlusOctant(octant);
2700 }
2701
2702 contentList.transfer(newContent);
2703 }
2704
2705 totalContents += contents_[contentI]->size();
2706 }
2707 else
2708 {
2709 // Empty, do nothing.
2710 }
2711 }
2712
2713 return totalContents;
2714}
2715
2716
2717template<class Type>
2719(
2721 const bool printContents,
2722 const label nodeI
2723) const
2724{
2725 const node& nod = nodes_[nodeI];
2726 const treeBoundBox& bb = nod.bb_;
2727
2728 os << "nodeI:" << nodeI << " bb:" << bb << nl
2729 << "parent:" << nod.parent_ << nl
2730 << "n:" << countElements(nodePlusOctant(nodeI, 0)) << nl;
2731
2732 for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
2733 {
2734 const treeBoundBox subBb(bb.subBbox(octant));
2735
2736 labelBits index = nod.subNodes_[octant];
2737
2738 if (isNode(index))
2739 {
2740 // tree node.
2741 label subNodeI = getNode(index);
2742
2743 os << "octant:" << octant
2744 << " node: n:" << countElements(index)
2745 << " bb:" << subBb << endl;
2746
2747 string oldPrefix = os.prefix();
2748 os.prefix() = " " + oldPrefix;
2749
2750 print(os, printContents, subNodeI);
2751
2752 os.prefix() = oldPrefix;
2753 }
2754 else if (isContent(index))
2755 {
2756 const labelList& indices = *(contents_[getContent(index)]);
2757
2758 if (false) //debug)
2759 {
2760 writeOBJ(nodeI, octant);
2761 }
2762
2763 os << "octant:" << octant
2764 << " content: n:" << indices.size()
2765 << " bb:" << subBb;
2766
2767 if (printContents)
2768 {
2769 os << " contents:";
2770 forAll(indices, j)
2771 {
2772 os << ' ' << indices[j];
2773 }
2774 }
2775 os << endl;
2776 }
2777 else
2778 {
2779 if (false) //debug)
2780 {
2781 writeOBJ(nodeI, octant);
2782 }
2783
2784 os << "octant:" << octant << " empty:" << subBb << endl;
2785 }
2786 }
2787}
2788
2789
2790template<class Type>
2792{
2793 label nEntries = 0;
2794 forAll(contents_, i)
2795 {
2796 nEntries += contents_[i]->size();
2797 }
2798
2799 Pout<< "indexedOctree<Type>::indexedOctree"
2800 << " : finished construction of tree of:" << shapes().typeName
2801 << nl
2802 << " bounding box: " << this->bb() << nl
2803 << " shapes: " << shapes().size() << nl
2804 << " treeNodes: " << nodes_.size() << nl
2805 << " nEntries: " << nEntries << nl
2806 << " levels/maxLevels: " << nLevelsMax_ << "/" << maxLevels_ << nl
2807 << " minSize: " << minSize_ << nl
2808 << " per treeLeaf: "
2809 << scalar(nEntries)/contents_.size() << nl
2810 << " per shape (duplicity):"
2811 << scalar(nEntries)/shapes().size() << nl
2812 << endl;
2813}
2814
2815
2816template<class Type>
2818{
2819 os << *this;
2820
2821 return os.good();
2822}
2823
2824
2825template<class Type>
2828{
2829 os << t.bb() << token::SPACE << t.nodes() << endl;
2830
2831 forAll(t.contents(), cI)
2832 {
2833 os << t.contents()[cI]() << endl;
2834 }
2835
2836 return os;
2837}
2838
2839
2840// ************************************************************************* //
bool success
Various functions to operate on Lists.
Y[inertIndex] max(0.0)
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
DynamicList< T, SizeMin > & shrink()
Shrink the allocated space to the number of elements used.
Definition: DynamicListI.H:434
void transfer(List< T > &list)
Transfer contents of the argument List into this.
Definition: DynamicListI.H:467
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
A 1D vector of objects of type <T> with a fixed length <N>.
Definition: FixedList.H:81
Minimal example by using system/controlDict.functions:
List< Key > toc() const
The table of contents (the keys) in unsorted order.
Definition: HashTable.C:122
bool good() const noexcept
True if next operation might succeed.
Definition: IOstream.H:233
static const List< label > & null()
Return a null List.
Definition: ListI.H:109
virtual const fileName & name() const
Get the name of the stream.
Definition: OSstream.H:107
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:62
int overlaps
Flag to control which overlap calculations are performed.
Definition: PDRparams.H:97
unsigned int remove()
Remove and return the last element.
Definition: PackedListI.H:709
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:66
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
const Cmpt & z() const
Access to the vector z component.
Definition: VectorI.H:85
const Cmpt & y() const
Access to the vector y component.
Definition: VectorI.H:79
const Cmpt & x() const
Access to the vector x component.
Definition: VectorI.H:73
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:91
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:97
vector span() const
The bounding box span (from minimum to maximum)
Definition: boundBoxI.H:127
Tree node. Has up pointer and down pointers.
FixedList< labelBits, 8 > subNodes_
IDs of the 8 nodes on all sides of the mid point.
label parent_
Parent node (index into nodes_ of tree)
treeBoundBox bb_
Bounding box of this node.
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.
const labelList & findIndices(const point &) const
Find the shape indices that occupy the result of findNode.
static scalar & perturbTol()
Get the perturbation tolerance.
const treeBoundBox & bb() const
Top bounding box.
label removeIndex(const label nodIndex, const label index)
const contentListList & contents() const
List of all contents (referenced by those nodes that are.
static volumeType getSide(const vector &outsideNormal, const vector &vec)
Helper function to return the side. Returns outside if.
const List< node > & nodes() const
List of all nodes.
bool insert(label startIndex, label endIndex)
Insert a new object into the tree.
labelBits findNode(const label nodeI, const point &) const
Find deepest node (as parent+octant) containing point. Starts.
bool insertIndex(const label nodIndex, const label index, label &nLevels)
pointIndexHit findLineAny(const point &start, const point &end) const
Find any intersection of line between start and end.
virtual bool write()
Write the output fields.
A 29bits label and 3bits direction packed into single label.
Definition: labelBits.H:54
A line primitive.
Definition: line.H:68
scalar print()
Print to screen.
static constexpr direction nComponents
Number of components in bool is 1.
Definition: bool.H:98
Version of OSstream that prints a prefix on each line.
@ SPACE
Space [isspace].
Definition: token.H:125
Standard boundBox with extra functionality for use in octree.
Definition: treeBoundBox.H:89
bool overlaps(const boundBox &bb) const
Overlaps other bounding box?
Definition: boundBoxI.H:221
void searchOrder(const point &pt, FixedList< direction, 8 > &octantOrder) const
Calculates optimal order to look for nearest to point.
treeBoundBox subBbox(const direction octant) const
Sub-box of given octant. Midpoint calculated.
Definition: treeBoundBox.C:106
bool contains(const vector &dir, const point &) const
Contains point (inside or on edge) and moving in direction.
Definition: treeBoundBox.C:320
direction subOctant(const point &pt) const
Returns octant number given point and the calculated midpoint.
Definition: treeBoundBoxI.H:84
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
An enumeration wrapper for classification of a location as being inside/outside of a volume.
Definition: volumeType.H:61
type
Volume classification types.
Definition: volumeType.H:66
@ OUTSIDE
A location outside the volume.
Definition: volumeType.H:69
@ MIXED
A location that is partly inside and outside.
Definition: volumeType.H:70
@ UNKNOWN
Unknown state.
Definition: volumeType.H:67
@ INSIDE
A location inside the volume.
Definition: volumeType.H:68
A class for handling words, derived from Foam::string.
Definition: word.H:68
PtrList< volScalarField > & Y
const volScalarField & p0
Definition: EEqn.H:36
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
OBJstream os(runTime.globalPath()/outputName)
const labelList nFaces(UPstream::listGatherValues< label >(aMesh.nFaces()))
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))
#define WarningInFunction
Report a warning using Foam::Warning.
rAUs append(new volScalarField(IOobject::groupName("rAU", phase1.name()), 1.0/(U1Eqn.A()+byDt(max(phase1.residualAlpha() - alpha1, scalar(0)) *rho1))))
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
static unsigned getOctant(const point &p)
List< label > labelList
A List of labels.
Definition: List.H:66
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
Definition: pointIndexHit.H:46
label facePoint(const int facei, const block &block, const label i, const label j)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
vector point
Point is a vector.
Definition: point.H:43
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
Definition: MSwindows.C:598
Ostream & operator<<(Ostream &, const boundaryPatch &p)
Write boundaryPatch as dictionary entries (without surrounding braces)
Definition: boundaryPatch.C:83
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
DynamicList< autoPtr< DynamicList< label > > > contentListList
errorManip< error > abort(error &err)
Definition: errorManip.H:144
uint8_t direction
Definition: direction.H:56
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
error FatalError
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
void divide(FieldField< Field, Type > &f, const FieldField< Field, Type > &f1, const FieldField< Field, scalar > &f2)
bool ln(const fileName &src, const fileName &dst)
Create a softlink. dst should not exist. Returns true if successful.
Definition: MSwindows.C:933
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333