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