FreeFOAM The Cross-Platform CFD Toolkit
indexedOctree.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) 1991-2010 OpenCFD Ltd.
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 "indexedOctree.H"
27 #include <OpenFOAM/linePointRef.H>
28 #include <meshTools/meshTools.H>
29 #include <OpenFOAM/OFstream.H>
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
33 template <class Type>
34 Foam::scalar Foam::indexedOctree<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 labelList& indices,
152  const treeBoundBox& bb,
153  labelListList& result
154 ) const
155 {
156  List<DynamicList<label> > subIndices(8);
157  for (direction octant = 0; octant < subIndices.size(); octant++)
158  {
159  subIndices[octant].setCapacity(indices.size()/8);
160  }
161 
162  // Precalculate bounding boxes.
163  FixedList<treeBoundBox, 8> subBbs;
164  for (direction octant = 0; octant < subBbs.size(); octant++)
165  {
166  subBbs[octant] = bb.subBbox(octant);
167  }
168 
169  forAll(indices, i)
170  {
171  label shapeI = indices[i];
172 
173  for (direction octant = 0; octant < 8; octant++)
174  {
175  if (shapes_.overlaps(shapeI, subBbs[octant]))
176  {
177  subIndices[octant].append(shapeI);
178  }
179  }
180  }
181 
182  result.setSize(8);
183  for (direction octant = 0; octant < subIndices.size(); octant++)
184  {
185  result[octant].transfer(subIndices[octant]);
186  }
187 }
188 
189 
190 // Subdivide the (content) node.
191 template <class Type>
194 (
195  const treeBoundBox& bb,
196  DynamicList<labelList>& contents,
197  const label contentI
198 ) const
199 {
200  const labelList& 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("indexedOctree<Type>::divide(..)")
212  << "Badly formed bounding box:" << bb
213  << abort(FatalError);
214  }
215 
216  nod.bb_ = bb;
217  nod.parent_ = -1;
218 
219  labelListList 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  labelList& 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  replaced = true;
238  }
239  else
240  {
241  // Store at end of contents.
242  // note dummy append + transfer trick
243  label sz = contents.size();
244  contents.append(labelList(0));
245  contents[sz].transfer(subIndices);
246  nod.subNodes_[octant] = contentPlusOctant(sz, octant);
247  }
248  }
249  else
250  {
251  // Mark node as empty
252  nod.subNodes_[octant] = emptyPlusOctant(octant);
253  }
254  }
255 
256  return nod;
257 }
258 
259 
260 // Split any contents node with more than minSize elements.
261 template <class Type>
263 (
264  const label minSize,
265  DynamicList<indexedOctree<Type>::node>& nodes,
266  DynamicList<labelList>& contents
267 ) const
268 {
269  label currentSize = nodes.size();
270 
271  // Take care to loop only over old nodes.
272  // Also we loop over the same DynamicList which gets modified and
273  // moved so make sure not to keep any references!
274  for (label nodeI = 0; nodeI < currentSize; nodeI++)
275  {
276  for
277  (
278  direction octant = 0;
279  octant < nodes[nodeI].subNodes_.size();
280  octant++
281  )
282  {
283  labelBits index = nodes[nodeI].subNodes_[octant];
284 
285  if (isNode(index))
286  {
287  // tree node. Leave intact.
288  }
289  else if (isContent(index))
290  {
291  label contentI = getContent(index);
292 
293  if (contents[contentI].size() > minSize)
294  {
295  // Create node for content.
296 
297  // Find the bounding box for the subnode
298  const node& nod = nodes[nodeI];
299  const treeBoundBox bb(nod.bb_.subBbox(octant));
300 
301  node subNode(divide(bb, contents, contentI));
302  subNode.parent_ = nodeI;
303  label sz = nodes.size();
304  nodes.append(subNode);
305  nodes[nodeI].subNodes_[octant] = nodePlusOctant(sz, octant);
306  }
307  }
308  }
309  }
310 }
311 
312 
313 // Reorder contents to be in same order as nodes. Returns number of nodes on
314 // the compactLevel.
315 template <class Type>
317 (
318  DynamicList<node>& nodes,
319  DynamicList<labelList>& contents,
320  const label compactLevel,
321  const label nodeI,
322  const label level,
323 
324  List<labelList>& compactedContents,
325  label& compactI
326 )
327 {
328  const node& nod = nodes[nodeI];
329 
330  label nNodes = 0;
331 
332  if (level < compactLevel)
333  {
334  for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
335  {
336  labelBits index = nod.subNodes_[octant];
337 
338  if (isNode(index))
339  {
340  nNodes += compactContents
341  (
342  nodes,
343  contents,
344  compactLevel,
345  getNode(index),
346  level+1,
347  compactedContents,
348  compactI
349  );
350  }
351  }
352  }
353  else if (level == compactLevel)
354  {
355  // Compact all content on this level
356  for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
357  {
358  labelBits index = nod.subNodes_[octant];
359 
360  if (isContent(index))
361  {
362  label contentI = getContent(index);
363 
364  compactedContents[compactI].transfer(contents[contentI]);
365 
366  // Subnode is at compactI. Adapt nodeI to point to it
367  nodes[nodeI].subNodes_[octant] =
368  contentPlusOctant(compactI, octant);
369 
370  compactI++;
371  }
372  else if (isNode(index))
373  {
374  nNodes++;
375  }
376  }
377  }
378  return nNodes;
379 }
380 
381 
382 // Pre-calculates wherever possible the volume status per node/subnode.
383 // Recurses to determine status of lowest level boxes. Level above is
384 // combination of octants below.
385 template <class Type>
388 (
389  const label nodeI
390 ) const
391 {
392  const node& nod = nodes_[nodeI];
393 
394  volumeType myType = UNKNOWN;
395 
396  for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
397  {
398  volumeType subType;
399 
400  labelBits index = nod.subNodes_[octant];
401 
402  if (isNode(index))
403  {
404  // tree node. Recurse.
405  subType = calcVolumeType(getNode(index));
406  }
407  else if (isContent(index))
408  {
409  // Contents. Depending on position in box might be on either
410  // side.
411  subType = MIXED;
412  }
413  else
414  {
415  // No data in this octant. Set type for octant acc. to the mid
416  // of its bounding box.
417  const treeBoundBox subBb = nod.bb_.subBbox(octant);
418 
419  subType = volumeType
420  (
421  shapes_.getVolumeType(*this, subBb.midpoint())
422  );
423  }
424 
425  // Store octant type
426  nodeTypes_.set((nodeI<<3)+octant, subType);
427 
428  // Combine sub node types into type for treeNode. Result is 'mixed' if
429  // types differ among subnodes.
430  if (myType == UNKNOWN)
431  {
432  myType = subType;
433  }
434  else if (subType != myType)
435  {
436  myType = MIXED;
437  }
438  }
439  return myType;
440 }
441 
442 
443 template <class Type>
446 (
447  const label nodeI,
448  const point& sample
449 ) const
450 {
451  const node& nod = nodes_[nodeI];
452 
453  direction octant = nod.bb_.subOctant(sample);
454 
455  volumeType octantType = volumeType(nodeTypes_.get((nodeI<<3)+octant));
456 
457  if (octantType == INSIDE)
458  {
459  return octantType;
460  }
461  else if (octantType == OUTSIDE)
462  {
463  return octantType;
464  }
465  else if (octantType == UNKNOWN)
466  {
467  // Can happen for e.g. non-manifold surfaces.
468  return octantType;
469  }
470  else if (octantType == MIXED)
471  {
472  labelBits index = nod.subNodes_[octant];
473 
474  if (isNode(index))
475  {
476  // Recurse
477  volumeType subType = getVolumeType(getNode(index), sample);
478 
479  return subType;
480  }
481  else if (isContent(index))
482  {
483  // Content. Defer to shapes.
484  return volumeType(shapes_.getVolumeType(*this, sample));
485  }
486  else
487  {
488  // Empty node. Cannot have 'mixed' as its type since not divided
489  // up and has no items inside it.
491  (
492  "indexedOctree<Type>::getVolumeType"
493  "(const label, const point&)"
494  ) << "Sample:" << sample << " node:" << nodeI
495  << " with bb:" << nodes_[nodeI].bb_ << nl
496  << "Empty subnode has invalid volume type MIXED."
497  << abort(FatalError);
498 
499  return UNKNOWN;
500  }
501  }
502  else
503  {
505  (
506  "indexedOctree<Type>::getVolumeType"
507  "(const label, const point&)"
508  ) << "Sample:" << sample << " at node:" << nodeI
509  << " octant:" << octant
510  << " with bb:" << nod.bb_.subBbox(octant) << nl
511  << "Node has invalid volume type " << octantType
512  << abort(FatalError);
513 
514  return UNKNOWN;
515  }
516 }
517 
518 
519 template <class Type>
522 (
523  const vector& outsideNormal,
524  const vector& vec
525 )
526 {
527  if ((outsideNormal&vec) >= 0)
528  {
529  return OUTSIDE;
530  }
531  else
532  {
533  return INSIDE;
534  }
535 }
536 
537 
538 //
539 // Query routines
540 // ~~~~~~~~~~~~~~
541 //
542 
543 // Find nearest point starting from nodeI
544 template <class Type>
546 (
547  const label nodeI,
548  const point& sample,
549 
550  scalar& nearestDistSqr,
551  label& nearestShapeI,
552  point& nearestPoint
553 ) const
554 {
555  const node& nod = nodes_[nodeI];
556 
557  // Determine order to walk through octants
558  FixedList<direction, 8> octantOrder;
559  nod.bb_.searchOrder(sample, octantOrder);
560 
561  // Go into all suboctants (one containing sample first) and update nearest.
562  for (direction i = 0; i < 8; i++)
563  {
564  direction octant = octantOrder[i];
565 
566  labelBits index = nod.subNodes_[octant];
567 
568  if (isNode(index))
569  {
570  label subNodeI = getNode(index);
571 
572  const treeBoundBox& subBb = nodes_[subNodeI].bb_;
573 
574  if (overlaps(subBb.min(), subBb.max(), nearestDistSqr, sample))
575  {
576  findNearest
577  (
578  subNodeI,
579  sample,
580 
581  nearestDistSqr,
582  nearestShapeI,
583  nearestPoint
584  );
585  }
586  }
587  else if (isContent(index))
588  {
589  if
590  (
591  overlaps
592  (
593  nod.bb_,
594  octant,
595  nearestDistSqr,
596  sample
597  )
598  )
599  {
600  shapes_.findNearest
601  (
602  contents_[getContent(index)],
603  sample,
604 
605  nearestDistSqr,
606  nearestShapeI,
607  nearestPoint
608  );
609  }
610  }
611  }
612 }
613 
614 
615 // Find nearest point to line.
616 template <class Type>
618 (
619  const label nodeI,
620  const linePointRef& ln,
621 
622  treeBoundBox& tightest,
623  label& nearestShapeI,
624  point& linePoint,
625  point& nearestPoint
626 ) const
627 {
628  const node& nod = nodes_[nodeI];
629  const treeBoundBox& nodeBb = nod.bb_;
630 
631  // Determine order to walk through octants
632  FixedList<direction, 8> octantOrder;
633  nod.bb_.searchOrder(ln.centre(), octantOrder);
634 
635  // Go into all suboctants (one containing sample first) and update nearest.
636  for (direction i = 0; i < 8; i++)
637  {
638  direction octant = octantOrder[i];
639 
640  labelBits index = nod.subNodes_[octant];
641 
642  if (isNode(index))
643  {
644  const treeBoundBox& subBb = nodes_[getNode(index)].bb_;
645 
646  if (subBb.overlaps(tightest))
647  {
648  findNearest
649  (
650  getNode(index),
651  ln,
652 
653  tightest,
654  nearestShapeI,
655  linePoint,
656  nearestPoint
657  );
658  }
659  }
660  else if (isContent(index))
661  {
662  const treeBoundBox subBb(nodeBb.subBbox(octant));
663 
664  if (subBb.overlaps(tightest))
665  {
666  shapes_.findNearest
667  (
668  contents_[getContent(index)],
669  ln,
670 
671  tightest,
672  nearestShapeI,
673  linePoint,
674  nearestPoint
675  );
676  }
677  }
678  }
679 }
680 
681 
682 template <class Type>
684 (
685  const label parentNodeI,
686  const direction octant
687 ) const
688 {
689  // Get type of node at octant
690  const node& nod = nodes_[parentNodeI];
691  labelBits index = nod.subNodes_[octant];
692 
693  if (isNode(index))
694  {
695  // Use stored bb
696  return nodes_[getNode(index)].bb_;
697  }
698  else
699  {
700  // Calculate subBb
701  return nod.bb_.subBbox(octant);
702  }
703 }
704 
705 
706 // Takes a bb and a point on/close to the edge of the bb and pushes the point
707 // inside by a small fraction.
708 template <class Type>
710 (
711  const treeBoundBox& bb,
712  const point& pt,
713  const bool pushInside
714 )
715 {
716  // Get local length scale.
717  const vector perturbVec = perturbTol_*(bb.span());
718 
719  point perturbedPt(pt);
720 
721  // Modify all components which are close to any face of the bb to be
722  // well inside/outside them.
723 
724  if (pushInside)
725  {
726  for (direction dir = 0; dir < vector::nComponents; dir++)
727  {
728  if (mag(pt[dir]-bb.min()[dir]) < mag(perturbVec[dir]))
729  {
730  // Close to 'left' side. Push well beyond left side.
731  scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
732  perturbedPt[dir] = bb.min()[dir] + perturbDist;
733  }
734  else if (mag(pt[dir]-bb.max()[dir]) < mag(perturbVec[dir]))
735  {
736  // Close to 'right' side. Push well beyond right side.
737  scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
738  perturbedPt[dir] = bb.max()[dir] - perturbDist;
739  }
740  }
741  }
742  else
743  {
744  for (direction dir = 0; dir < vector::nComponents; dir++)
745  {
746  if (mag(pt[dir]-bb.min()[dir]) < mag(perturbVec[dir]))
747  {
748  scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
749  perturbedPt[dir] = bb.min()[dir] - perturbDist;
750  }
751  else if (mag(pt[dir]-bb.max()[dir]) < mag(perturbVec[dir]))
752  {
753  scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
754  perturbedPt[dir] = bb.max()[dir] + perturbDist;
755  }
756  }
757  }
758 
759  if (debug)
760  {
761  if (pushInside != bb.contains(perturbedPt))
762  {
763  FatalErrorIn("indexedOctree<Type>::pushPoint(..)")
764  << "pushed point:" << pt
765  << " to:" << perturbedPt
766  << " wanted side:" << pushInside
767  << " obtained side:" << bb.contains(perturbedPt)
768  << " of bb:" << bb
769  << abort(FatalError);
770  }
771  }
772 
773  return perturbedPt;
774 }
775 
776 
777 // Takes a bb and a point on the edge of the bb and pushes the point
778 // outside by a small fraction.
779 template <class Type>
781 (
782  const treeBoundBox& bb,
783  const direction faceID,
784  const point& pt,
785  const bool pushInside
786 )
787 {
788  // Get local length scale.
789  const vector perturbVec = perturbTol_*bb.span();
790 
791  point perturbedPt(pt);
792 
793  // Modify all components which are close to any face of the bb to be
794  // well outside them.
795 
796  if (faceID == 0)
797  {
798  FatalErrorIn("indexedOctree<Type>::pushPoint(..)")
799  << abort(FatalError);
800  }
801 
802  if (faceID & treeBoundBox::LEFTBIT)
803  {
804  if (pushInside)
805  {
806  perturbedPt[0] = bb.min()[0] + (perturbVec[0] + ROOTVSMALL);
807  }
808  else
809  {
810  perturbedPt[0] = bb.min()[0] - (perturbVec[0] + ROOTVSMALL);
811  }
812  }
813  else if (faceID & treeBoundBox::RIGHTBIT)
814  {
815  if (pushInside)
816  {
817  perturbedPt[0] = bb.max()[0] - (perturbVec[0] + ROOTVSMALL);
818  }
819  else
820  {
821  perturbedPt[0] = bb.max()[0] + (perturbVec[0] + ROOTVSMALL);
822  }
823  }
824 
825  if (faceID & treeBoundBox::BOTTOMBIT)
826  {
827  if (pushInside)
828  {
829  perturbedPt[1] = bb.min()[1] + (perturbVec[1] + ROOTVSMALL);
830  }
831  else
832  {
833  perturbedPt[1] = bb.min()[1] - (perturbVec[1] + ROOTVSMALL);
834  }
835  }
836  else if (faceID & treeBoundBox::TOPBIT)
837  {
838  if (pushInside)
839  {
840  perturbedPt[1] = bb.max()[1] - (perturbVec[1] + ROOTVSMALL);
841  }
842  else
843  {
844  perturbedPt[1] = bb.max()[1] + (perturbVec[1] + ROOTVSMALL);
845  }
846  }
847 
848  if (faceID & treeBoundBox::BACKBIT)
849  {
850  if (pushInside)
851  {
852  perturbedPt[2] = bb.min()[2] + (perturbVec[2] + ROOTVSMALL);
853  }
854  else
855  {
856  perturbedPt[2] = bb.min()[2] - (perturbVec[2] + ROOTVSMALL);
857  }
858  }
859  else if (faceID & treeBoundBox::FRONTBIT)
860  {
861  if (pushInside)
862  {
863  perturbedPt[2] = bb.max()[2] - (perturbVec[2] + ROOTVSMALL);
864  }
865  else
866  {
867  perturbedPt[2] = bb.max()[2] + (perturbVec[2] + ROOTVSMALL);
868  }
869  }
870 
871  if (debug)
872  {
873  if (pushInside != bb.contains(perturbedPt))
874  {
875  FatalErrorIn("indexedOctree<Type>::pushPoint(..)")
876  << "pushed point:" << pt << " on face:" << faceString(faceID)
877  << " to:" << perturbedPt
878  << " wanted side:" << pushInside
879  << " obtained side:" << bb.contains(perturbedPt)
880  << " of bb:" << bb
881  << abort(FatalError);
882  }
883  }
884 
885  return perturbedPt;
886 }
887 
888 
889 // Guarantees that if pt is on a face it gets perturbed so it is away
890 // from the face edges.
891 // If pt is not on a face does nothing.
892 template <class Type>
894 (
895  const treeBoundBox& bb,
896  const vector& dir, // end-start
897  const point& pt
898 )
899 {
900  if (debug)
901  {
902  if (bb.posBits(pt) != 0)
903  {
904  FatalErrorIn("indexedOctree<Type>::pushPointIntoFace(..)")
905  << " bb:" << bb << endl
906  << "does not contain point " << pt << abort(FatalError);
907  }
908  }
909 
910 
911  // Handle two cases:
912  // - point exactly on multiple faces. Push away from all but one.
913  // - point on a single face. Push away from edges of face.
914 
915  direction ptFaceID = bb.faceBits(pt);
916 
917  direction nFaces = 0;
918  FixedList<direction, 3> faceIndices;
919 
920  if (ptFaceID & treeBoundBox::LEFTBIT)
921  {
922  faceIndices[nFaces++] = treeBoundBox::LEFT;
923  }
924  else if (ptFaceID & treeBoundBox::RIGHTBIT)
925  {
926  faceIndices[nFaces++] = treeBoundBox::RIGHT;
927  }
928 
929  if (ptFaceID & treeBoundBox::BOTTOMBIT)
930  {
931  faceIndices[nFaces++] = treeBoundBox::BOTTOM;
932  }
933  else if (ptFaceID & treeBoundBox::TOPBIT)
934  {
935  faceIndices[nFaces++] = treeBoundBox::TOP;
936  }
937 
938  if (ptFaceID & treeBoundBox::BACKBIT)
939  {
940  faceIndices[nFaces++] = treeBoundBox::BACK;
941  }
942  else if (ptFaceID & treeBoundBox::FRONTBIT)
943  {
944  faceIndices[nFaces++] = treeBoundBox::FRONT;
945  }
946 
947 
948  // Determine the face we want to keep the point on
949 
950  direction keepFaceID;
951 
952  if (nFaces == 0)
953  {
954  // Return original point
955  return pt;
956  }
957  else if (nFaces == 1)
958  {
959  // Point is on a single face
960  keepFaceID = faceIndices[0];
961  }
962  else
963  {
964  // Determine best face out of faceIndices[0 .. nFaces-1].
965  // The best face is the one most perpendicular to the ray direction.
966 
967  keepFaceID = faceIndices[0];
968  scalar maxInproduct = mag(treeBoundBox::faceNormals[keepFaceID] & dir);
969 
970  for (direction i = 1; i < nFaces; i++)
971  {
972  direction face = faceIndices[i];
973  scalar s = mag(treeBoundBox::faceNormals[face] & dir);
974  if (s > maxInproduct)
975  {
976  maxInproduct = s;
977  keepFaceID = face;
978  }
979  }
980  }
981 
982 
983  // 1. Push point into bb, away from all corners
984 
985  point facePoint(pushPoint(bb, pt, true));
986  direction faceID = 0;
987 
988  // 2. Snap it back onto the preferred face
989 
990  if (keepFaceID == treeBoundBox::LEFT)
991  {
992  facePoint.x() = bb.min().x();
993  faceID = treeBoundBox::LEFTBIT;
994  }
995  else if (keepFaceID == treeBoundBox::RIGHT)
996  {
997  facePoint.x() = bb.max().x();
998  faceID = treeBoundBox::RIGHTBIT;
999  }
1000  else if (keepFaceID == treeBoundBox::BOTTOM)
1001  {
1002  facePoint.y() = bb.min().y();
1003  faceID = treeBoundBox::BOTTOMBIT;
1004  }
1005  else if (keepFaceID == treeBoundBox::TOP)
1006  {
1007  facePoint.y() = bb.max().y();
1008  faceID = treeBoundBox::TOPBIT;
1009  }
1010  else if (keepFaceID == treeBoundBox::BACK)
1011  {
1012  facePoint.z() = bb.min().z();
1013  faceID = treeBoundBox::BACKBIT;
1014  }
1015  else if (keepFaceID == treeBoundBox::FRONT)
1016  {
1017  facePoint.z() = bb.max().z();
1018  faceID = treeBoundBox::FRONTBIT;
1019  }
1020 
1021 
1022  if (debug)
1023  {
1024  if (faceID != bb.faceBits(facePoint))
1025  {
1026  FatalErrorIn("indexedOctree<Type>::pushPointIntoFace(..)")
1027  << "Pushed point from " << pt
1028  << " on face:" << ptFaceID << " of bb:" << bb << endl
1029  << "onto " << facePoint
1030  << " on face:" << faceID
1031  << " which is not consistent with geometric face "
1032  << bb.faceBits(facePoint)
1033  << abort(FatalError);
1034  }
1035  if (bb.posBits(facePoint) != 0)
1036  {
1037  FatalErrorIn("indexedOctree<Type>::pushPointIntoFace(..)")
1038  << " bb:" << bb << endl
1039  << "does not contain perturbed point "
1040  << facePoint << abort(FatalError);
1041  }
1042  }
1043 
1044 
1045  return facePoint;
1046 }
1047 
1048 
1051 //template <class Type>
1052 //void Foam::indexedOctree<Type>::checkMultipleFaces
1053 //(
1054 // const treeBoundBox& bb,
1055 // const vector& dir, // end-start
1056 // pointIndexHit& faceHitInfo,
1057 // direction& faceID
1058 //)
1059 //{
1060 // // Do the quick elimination of no or one face.
1061 // if
1062 // (
1063 // (faceID == 0)
1064 // || (faceID == treeBoundBox::LEFTBIT)
1065 // || (faceID == treeBoundBox::RIGHTBIT)
1066 // || (faceID == treeBoundBox::BOTTOMBIT)
1067 // || (faceID == treeBoundBox::TOPBIT)
1068 // || (faceID == treeBoundBox::BACKBIT)
1069 // || (faceID == treeBoundBox::FRONTBIT)
1070 // )
1071 // {
1072 // return;
1073 // }
1074 //
1075 //
1076 // // Check the direction of vector w.r.t. faces being intersected.
1077 // FixedList<scalar, 6> inproducts(-GREAT);
1078 //
1079 // direction nFaces = 0;
1080 //
1081 // if (faceID & treeBoundBox::LEFTBIT)
1082 // {
1083 // inproducts[treeBoundBox::LEFT] = mag
1084 // (
1085 // treeBoundBox::faceNormals[treeBoundBox::LEFT]
1086 // & dir
1087 // );
1088 // nFaces++;
1089 // }
1090 // if (faceID & treeBoundBox::RIGHTBIT)
1091 // {
1092 // inproducts[treeBoundBox::RIGHT] = mag
1093 // (
1094 // treeBoundBox::faceNormals[treeBoundBox::RIGHT]
1095 // & dir
1096 // );
1097 // nFaces++;
1098 // }
1099 //
1100 // if (faceID & treeBoundBox::BOTTOMBIT)
1101 // {
1102 // inproducts[treeBoundBox::BOTTOM] = mag
1103 // (
1104 // treeBoundBox::faceNormals[treeBoundBox::BOTTOM]
1105 // & dir
1106 // );
1107 // nFaces++;
1108 // }
1109 // if (faceID & treeBoundBox::TOPBIT)
1110 // {
1111 // inproducts[treeBoundBox::TOP] = mag
1112 // (
1113 // treeBoundBox::faceNormals[treeBoundBox::TOP]
1114 // & dir
1115 // );
1116 // nFaces++;
1117 // }
1118 //
1119 // if (faceID & treeBoundBox::BACKBIT)
1120 // {
1121 // inproducts[treeBoundBox::BACK] = mag
1122 // (
1123 // treeBoundBox::faceNormals[treeBoundBox::BACK]
1124 // & dir
1125 // );
1126 // nFaces++;
1127 // }
1128 // if (faceID & treeBoundBox::FRONTBIT)
1129 // {
1130 // inproducts[treeBoundBox::FRONT] = mag
1131 // (
1132 // treeBoundBox::faceNormals[treeBoundBox::FRONT]
1133 // & dir
1134 // );
1135 // nFaces++;
1136 // }
1137 //
1138 // if (nFaces == 0 || nFaces == 1 || nFaces > 3)
1139 // {
1140 // FatalErrorIn("indexedOctree<Type>::checkMultipleFaces(..)")
1141 // << "Problem : nFaces:" << nFaces << abort(FatalError);
1142 // }
1143 //
1144 // // Keep point on most perpendicular face; shift it away from the aligned
1145 // // ones.
1146 // // E.g. line hits top and left face:
1147 // // a
1148 // // ----+----+
1149 // // | |
1150 // // | |
1151 // // +----+
1152 // // Shift point down (away from top):
1153 // //
1154 // // a+----+
1155 // // ----| |
1156 // // | |
1157 // // +----+
1158 //
1159 // label maxIndex = -1;
1160 // scalar maxInproduct = -GREAT;
1161 //
1162 // for (direction i = 0; i < 6; i++)
1163 // {
1164 // if (inproducts[i] > maxInproduct)
1165 // {
1166 // maxInproduct = inproducts[i];
1167 // maxIndex = i;
1168 // }
1169 // }
1170 //
1171 // if (maxIndex == -1)
1172 // {
1173 // FatalErrorIn("indexedOctree<Type>::checkMultipleFaces(..)")
1174 // << "Problem maxIndex:" << maxIndex << " inproducts:" << inproducts
1175 // << abort(FatalError);
1176 // }
1177 //
1178 // const point oldPoint(faceHitInfo.rawPoint());
1179 // const direction oldFaceID = faceID;
1180 //
1181 // // 1. Push point into bb, away from all corners
1182 //
1183 // faceHitInfo.rawPoint() = pushPoint(bb, oldFaceID, oldPoint, true);
1184 //
1185 // // 2. Snap it back onto the preferred face
1186 //
1187 // if (maxIndex == treeBoundBox::LEFT)
1188 // {
1189 // faceHitInfo.rawPoint().x() = bb.min().x();
1190 // faceID = treeBoundBox::LEFTBIT;
1191 // }
1192 // else if (maxIndex == treeBoundBox::RIGHT)
1193 // {
1194 // faceHitInfo.rawPoint().x() = bb.max().x();
1195 // faceID = treeBoundBox::RIGHTBIT;
1196 // }
1197 // else if (maxIndex == treeBoundBox::BOTTOM)
1198 // {
1199 // faceHitInfo.rawPoint().y() = bb.min().y();
1200 // faceID = treeBoundBox::BOTTOMBIT;
1201 // }
1202 // else if (maxIndex == treeBoundBox::TOP)
1203 // {
1204 // faceHitInfo.rawPoint().y() = bb.max().y();
1205 // faceID = treeBoundBox::TOPBIT;
1206 // }
1207 // else if (maxIndex == treeBoundBox::BACK)
1208 // {
1209 // faceHitInfo.rawPoint().z() = bb.min().z();
1210 // faceID = treeBoundBox::BACKBIT;
1211 // }
1212 // else if (maxIndex == treeBoundBox::FRONT)
1213 // {
1214 // faceHitInfo.rawPoint().z() = bb.max().z();
1215 // faceID = treeBoundBox::FRONTBIT;
1216 // }
1217 //
1218 // Pout<< "From ray:" << dir
1219 // << " from point:" << oldPoint
1220 // << " on faces:" << faceString(oldFaceID)
1221 // << " of bb:" << bb
1222 // << " with inprods:" << inproducts
1223 // << " maxIndex:" << maxIndex << endl
1224 // << "perturbed to point:" << faceHitInfo.rawPoint()
1225 // << " on face:" << faceString(faceID)
1226 // << endl;
1227 //
1228 //
1229 // if (debug)
1230 // {
1231 // if (faceID != bb.faceBits(faceHitInfo.rawPoint()))
1232 // {
1233 // FatalErrorIn("indexedOctree<Type>::checkMultipleFaces(..)")
1234 // << "Pushed point from " << oldPoint
1235 // << " on face:" << oldFaceID << " of bb:" << bb << endl
1236 // << "onto " << faceHitInfo.rawPoint()
1237 // << " on face:" << faceID
1238 // << " which is not consistent with geometric face "
1239 // << bb.faceBits(faceHitInfo.rawPoint())
1240 // << abort(FatalError);
1241 // }
1242 // }
1243 //}
1244 
1245 
1246 // Get parent node and octant. Return false if top of tree reached.
1247 template <class Type>
1249 (
1250  const label nodeI,
1251  const direction octant,
1252 
1253  label& parentNodeI,
1254  label& parentOctant
1255 ) const
1256 {
1257  parentNodeI = nodes_[nodeI].parent_;
1258 
1259  if (parentNodeI == -1)
1260  {
1261  // Reached edge of tree
1262  return false;
1263  }
1264 
1265  const node& parentNode = nodes_[parentNodeI];
1266 
1267  // Find octant nodeI is in.
1268  parentOctant = 255;
1269 
1270  for (direction i = 0; i < parentNode.subNodes_.size(); i++)
1271  {
1272  labelBits index = parentNode.subNodes_[i];
1273 
1274  if (isNode(index) && getNode(index) == nodeI)
1275  {
1276  parentOctant = i;
1277  break;
1278  }
1279  }
1280 
1281  if (parentOctant == 255)
1282  {
1283  FatalErrorIn("walkToParent(..)")
1284  << "Problem: no parent found for octant:" << octant
1285  << " node:" << nodeI
1286  << abort(FatalError);
1287  }
1288 
1289  return true;
1290 }
1291 
1292 
1293 // Walk tree to neighbouring node. Gets current position as
1294 // node and octant in this node and walks in the direction given by
1295 // the facePointBits (combination of treeBoundBox::LEFTBIT, TOPBIT etc.)
1296 // Returns false if edge of tree hit.
1297 template <class Type>
1299 (
1300  const point& facePoint,
1301  const direction faceID, // face(s) that facePoint is on
1302  label& nodeI,
1303  direction& octant
1304 ) const
1305 {
1306  label oldNodeI = nodeI;
1307  direction oldOctant = octant;
1308 
1309  // Find out how to test for octant. Say if we want to go left we need
1310  // to traverse up the tree until we hit a node where our octant is
1311  // on the right.
1312 
1313  // Coordinate direction to test
1314  const direction X = treeBoundBox::RIGHTHALF;
1315  const direction Y = treeBoundBox::TOPHALF;
1316  const direction Z = treeBoundBox::FRONTHALF;
1317 
1318  direction octantMask = 0;
1319  direction wantedValue = 0;
1320 
1321  if ((faceID & treeBoundBox::LEFTBIT) != 0)
1322  {
1323  // We want to go left so check if in right octant (i.e. x-bit is set)
1324  octantMask |= X;
1325  wantedValue |= X;
1326  }
1327  else if ((faceID & treeBoundBox::RIGHTBIT) != 0)
1328  {
1329  octantMask |= X; // wantedValue already 0
1330  }
1331 
1332  if ((faceID & treeBoundBox::BOTTOMBIT) != 0)
1333  {
1334  // Want to go down so check for y-bit set.
1335  octantMask |= Y;
1336  wantedValue |= Y;
1337  }
1338  else if ((faceID & treeBoundBox::TOPBIT) != 0)
1339  {
1340  // Want to go up so check for y-bit not set.
1341  octantMask |= Y;
1342  }
1343 
1344  if ((faceID & treeBoundBox::BACKBIT) != 0)
1345  {
1346  octantMask |= Z;
1347  wantedValue |= Z;
1348  }
1349  else if ((faceID & treeBoundBox::FRONTBIT) != 0)
1350  {
1351  octantMask |= Z;
1352  }
1353 
1354  // So now we have the coordinate directions in the octant we need to check
1355  // and the resulting value.
1356 
1357  /*
1358  // +---+---+
1359  // | | |
1360  // | | |
1361  // | | |
1362  // +---+-+-+
1363  // | | | |
1364  // | a+-+-+
1365  // | |\| |
1366  // +---+-+-+
1367  // \
1368  //
1369  // e.g. ray is at (a) in octant 0(or 4) with faceIDs : LEFTBIT+TOPBIT.
1370  // If we would be in octant 1(or 5) we could go to the correct octant
1371  // in the same node by just flipping the x and y bits (exoring).
1372  // But if we are not in octant 1/5 we have to go up until we are.
1373  // In general for leftbit+topbit:
1374  // - we have to check for x and y : octantMask = 011
1375  // - the checked bits have to be : wantedValue = ?01
1376  */
1377 
1378  //Pout<< "For point " << facePoint << endl;
1379 
1380  // Go up until we have chance to cross to the wanted direction
1381  while (wantedValue != (octant & octantMask))
1382  {
1383  // Go up to the parent.
1384 
1385  // Remove the directions that are not on the boundary of the parent.
1386  // See diagram above
1387  if (wantedValue & X) // && octantMask&X
1388  {
1389  // Looking for right octant.
1390  if (octant & X)
1391  {
1392  // My octant is one of the left ones so punch out the x check
1393  octantMask &= ~X;
1394  wantedValue &= ~X;
1395  }
1396  }
1397  else
1398  {
1399  if (!(octant & X))
1400  {
1401  // My octant is right but I want to go left.
1402  octantMask &= ~X;
1403  wantedValue &= ~X;
1404  }
1405  }
1406 
1407  if (wantedValue & Y)
1408  {
1409  if (octant & Y)
1410  {
1411  octantMask &= ~Y;
1412  wantedValue &= ~Y;
1413  }
1414  }
1415  else
1416  {
1417  if (!(octant & Y))
1418  {
1419  octantMask &= ~Y;
1420  wantedValue &= ~Y;
1421  }
1422  }
1423 
1424  if (wantedValue & Z)
1425  {
1426  if (octant & Z)
1427  {
1428  octantMask &= ~Z;
1429  wantedValue &= ~Z;
1430  }
1431  }
1432  else
1433  {
1434  if (!(octant & Z))
1435  {
1436  octantMask &= ~Z;
1437  wantedValue &= ~Z;
1438  }
1439  }
1440 
1441 
1442  label parentNodeI;
1443  label parentOctant;
1444  walkToParent(nodeI, octant, parentNodeI, parentOctant);
1445 
1446  if (parentNodeI == -1)
1447  {
1448  // Reached edge of tree
1449  return false;
1450  }
1451 
1452  //Pout<< " walked from node:" << nodeI << " octant:" << octant
1453  // << " bb:" << nodes_[nodeI].bb_.subBbox(octant) << endl
1454  // << " to:" << parentNodeI << " octant:" << parentOctant
1455  // << " bb:" << nodes_[parentNodeI].bb_.subBbox(parentOctant)
1456  // << endl;
1457  //
1458  //Pout<< " octantMask:" << octantMask
1459  // << " wantedValue:" << wantedValue << endl;
1460 
1461  nodeI = parentNodeI;
1462  octant = parentOctant;
1463  }
1464 
1465  // So now we hit the correct parent node. Switch to the correct octant.
1466  // Is done by jumping to the 'other half' so if e.g. in x direction in
1467  // right half we now jump to the left half.
1468  octant ^= octantMask;
1469 
1470  //Pout<< " to node:" << nodeI << " octant:" << octant
1471  // << " subBb:" <<subBbox(nodeI, octant) << endl;
1472 
1473 
1474  if (debug)
1475  {
1476  const treeBoundBox subBb(subBbox(nodeI, octant));
1477 
1478  if (!subBb.contains(facePoint))
1479  {
1480  FatalErrorIn("indexedOctree<Type>::walkToNeighbour(..)")
1481  << "When searching for " << facePoint
1482  << " ended up in node:" << nodeI
1483  << " octant:" << octant
1484  << " with bb:" << subBb
1485  << abort(FatalError);
1486  }
1487  }
1488 
1489 
1490  // See if we need to travel down. Note that we already go into the
1491  // the first level ourselves (instead of having findNode decide)
1492  labelBits index = nodes_[nodeI].subNodes_[octant];
1493 
1494  if (isNode(index))
1495  {
1496  labelBits node = findNode(getNode(index), facePoint);
1497 
1498  nodeI = getNode(node);
1499  octant = getOctant(node);
1500  }
1501 
1502 
1503  if (debug)
1504  {
1505  const treeBoundBox subBb(subBbox(nodeI, octant));
1506 
1507  if (nodeI == oldNodeI && octant == oldOctant)
1508  {
1509  FatalErrorIn("indexedOctree<Type>::walkToNeighbour(..)")
1510  << "Did not go to neighbour when searching for " << facePoint
1511  << endl
1512  << " starting from face:" << faceString(faceID)
1513  << " node:" << nodeI
1514  << " octant:" << octant
1515  << " bb:" << subBb
1516  << abort(FatalError);
1517  }
1518 
1519  if (!subBb.contains(facePoint))
1520  {
1521  FatalErrorIn("indexedOctree<Type>::walkToNeighbour(..)")
1522  << "When searching for " << facePoint
1523  << " ended up in node:" << nodeI
1524  << " octant:" << octant
1525  << " bb:" << subBb
1526  << abort(FatalError);
1527  }
1528  }
1529 
1530 
1531  return true;
1532 }
1533 
1534 
1535 template <class Type>
1537 (
1538  const direction faceID
1539 )
1540 {
1541  word desc;
1542 
1543  if (faceID == 0)
1544  {
1545  desc = "noFace";
1546  }
1547  if (faceID & treeBoundBox::LEFTBIT)
1548  {
1549  if (!desc.empty()) desc += "+";
1550  desc += "left";
1551  }
1552  if (faceID & treeBoundBox::RIGHTBIT)
1553  {
1554  if (!desc.empty()) desc += "+";
1555  desc += "right";
1556  }
1557  if (faceID & treeBoundBox::BOTTOMBIT)
1558  {
1559  if (!desc.empty()) desc += "+";
1560  desc += "bottom";
1561  }
1562  if (faceID & treeBoundBox::TOPBIT)
1563  {
1564  if (!desc.empty()) desc += "+";
1565  desc += "top";
1566  }
1567  if (faceID & treeBoundBox::BACKBIT)
1568  {
1569  if (!desc.empty()) desc += "+";
1570  desc += "back";
1571  }
1572  if (faceID & treeBoundBox::FRONTBIT)
1573  {
1574  if (!desc.empty()) desc += "+";
1575  desc += "front";
1576  }
1577  return desc;
1578 }
1579 
1580 
1581 // Traverse a node. If intersects a triangle return first intersection point:
1582 // hitInfo.index = index of shape
1583 // hitInfo.point = point on shape
1584 // Else return a miss and the bounding box face hit:
1585 // hitInfo.point = coordinate of intersection of ray with bounding box
1586 // hitBits = posbits of point on bounding box
1587 template <class Type>
1589 (
1590  const bool findAny,
1591  const point& treeStart,
1592  const vector& treeVec,
1593 
1594  const point& start,
1595  const point& end,
1596  const label nodeI,
1597  const direction octant,
1598 
1599  pointIndexHit& hitInfo,
1600  direction& hitBits
1601 ) const
1602 {
1603  if (debug)
1604  {
1605  const treeBoundBox octantBb(subBbox(nodeI, octant));
1606 
1607  if (octantBb.posBits(start) != 0)
1608  {
1609  FatalErrorIn("indexedOctree<Type>::traverseNode(..)")
1610  << "Node:" << nodeI << " octant:" << octant
1611  << " bb:" << octantBb << endl
1612  << "does not contain point " << start << abort(FatalError);
1613  }
1614  }
1615 
1616 
1617  const node& nod = nodes_[nodeI];
1618 
1619  labelBits index = nod.subNodes_[octant];
1620 
1621  if (isContent(index))
1622  {
1623  const labelList& indices = contents_[getContent(index)];
1624 
1625  if (findAny)
1626  {
1627  // Find any intersection
1628 
1629  forAll(indices, elemI)
1630  {
1631  label shapeI = indices[elemI];
1632 
1633  point pt;
1634  bool hit = shapes_.intersects(shapeI, start, end, pt);
1635 
1636  if (hit)
1637  {
1638  // Hit so pt is nearer than nearestPoint.
1639  // Update hit info
1640  hitInfo.setHit();
1641  hitInfo.setIndex(shapeI);
1642  hitInfo.setPoint(pt);
1643  return;
1644  }
1645  }
1646  }
1647  else
1648  {
1649  // Find nearest intersection.
1650 
1651  point nearestPoint(end);
1652 
1653  forAll(indices, elemI)
1654  {
1655  label shapeI = indices[elemI];
1656 
1657  point pt;
1658  bool hit = shapes_.intersects(shapeI, start, nearestPoint, pt);
1659 
1660  if (hit)
1661  {
1662  // Hit so pt is nearer than nearestPoint.
1663  nearestPoint = pt;
1664  // Update hit info
1665  hitInfo.setHit();
1666  hitInfo.setIndex(shapeI);
1667  hitInfo.setPoint(pt);
1668  }
1669  }
1670 
1671  if (hitInfo.hit())
1672  {
1673  // Found intersected shape.
1674  return;
1675  }
1676  }
1677  }
1678 
1679  // Nothing intersected in this node
1680  // Traverse to end of node. Do by ray tracing back from end and finding
1681  // intersection with bounding box of node.
1682  // start point is now considered to be inside bounding box.
1683 
1684  const treeBoundBox octantBb(subBbox(nodeI, octant));
1685 
1686  point pt;
1687  bool intersected = octantBb.intersects
1688  (
1689  end, //treeStart,
1690  (start-end), //treeVec,
1691 
1692  end,
1693  start,
1694 
1695  pt,
1696  hitBits
1697  );
1698 
1699 
1700  if (intersected)
1701  {
1702  // Return miss. Set misspoint to face.
1703  hitInfo.setPoint(pt);
1704  }
1705  else
1706  {
1707  // Rare case: did not find intersection of ray with octantBb. Can
1708  // happen if end is on face/edge of octantBb. Do like
1709  // lagrangian and re-track with perturbed vector (slightly pushed out
1710  // of bounding box)
1711 
1712  point perturbedEnd(pushPoint(octantBb, end, false));
1713 
1714 
1715  //if (debug)
1716  {
1717  // Dump octantBb to obj
1718  writeOBJ(nodeI, octant);
1719  // Dump ray to obj as well
1720  {
1721  OFstream str("ray.obj");
1722  meshTools::writeOBJ(str, start);
1723  meshTools::writeOBJ(str, end);
1724  str << "l 1 2" << nl;
1725  }
1726  WarningIn("indexedOctree<Type>::traverseNode(..)")
1727  << "Did not intersect ray from endpoint:" << end
1728  << " to startpoint:" << start
1729  << " with bounding box:" << octantBb << nl
1730  << "Re-intersecting with perturbed endpoint:" << perturbedEnd
1731  << endl;
1732  }
1733 
1734  traverseNode
1735  (
1736  findAny,
1737  treeStart,
1738  treeVec,
1739  start,
1740  perturbedEnd,
1741  nodeI,
1742  octant,
1743 
1744  hitInfo,
1745  hitBits
1746  );
1747  }
1748 }
1749 
1750 
1751 // Find first intersection
1752 template <class Type>
1754 (
1755  const bool findAny,
1756  const point& treeStart,
1757  const point& treeEnd,
1758  const label startNodeI,
1759  const direction startOctant,
1760  const bool verbose
1761 ) const
1762 {
1763  const vector treeVec(treeEnd - treeStart);
1764 
1765  // Current node as parent+octant
1766  label nodeI = startNodeI;
1767  direction octant = startOctant;
1768 
1769  if (verbose)
1770  {
1771  Pout<< "findLine : treeStart:" << treeStart
1772  << " treeEnd:" << treeEnd << endl
1773  << "node:" << nodeI
1774  << " octant:" << octant
1775  << " bb:" << subBbox(nodeI, octant) << endl;
1776  }
1777 
1778  // Current position. Initialize to miss
1779  pointIndexHit hitInfo(false, treeStart, -1);
1780 
1781  //while (true)
1782  label i = 0;
1783  for (; i < 100000; i++)
1784  {
1785  // Ray-trace to end of current node. Updates point (either on triangle
1786  // in case of hit or on node bounding box in case of miss)
1787 
1788  const treeBoundBox octantBb(subBbox(nodeI, octant));
1789 
1790  // Make sure point is away from any edges/corners
1791  point startPoint
1792  (
1793  pushPointIntoFace
1794  (
1795  octantBb,
1796  treeVec,
1797  hitInfo.rawPoint()
1798  )
1799  );
1800 
1801  if (verbose)
1802  {
1803  Pout<< "iter:" << i
1804  << " at current:" << hitInfo.rawPoint()
1805  << " (perturbed:" << startPoint << ")" << endl
1806  << " node:" << nodeI
1807  << " octant:" << octant
1808  << " bb:" << subBbox(nodeI, octant) << endl;
1809  }
1810 
1811 
1812 
1813 
1814  // Faces of current bounding box current point is on
1815  direction hitFaceID = 0;
1816 
1817  traverseNode
1818  (
1819  findAny,
1820  treeStart,
1821  treeVec,
1822 
1823  startPoint, // Note: pass in copy since hitInfo
1824  // also used in return value.
1825  treeEnd, // pass in overall end so is nicely outside bb
1826  nodeI,
1827  octant,
1828 
1829  hitInfo,
1830  hitFaceID
1831  );
1832 
1833  // Did we hit a triangle?
1834  if (hitInfo.hit())
1835  {
1836  break;
1837  }
1838 
1839  if (hitFaceID == 0 || hitInfo.rawPoint() == treeEnd)
1840  {
1841  // endpoint inside the tree. Return miss.
1842  break;
1843  }
1844 
1845  // Create a point on other side of face.
1846  point perturbedPoint
1847  (
1848  pushPoint
1849  (
1850  octantBb,
1851  hitFaceID,
1852  hitInfo.rawPoint(),
1853  false // push outside of octantBb
1854  )
1855  );
1856 
1857  if (verbose)
1858  {
1859  Pout<< " iter:" << i
1860  << " hit face:" << faceString(hitFaceID)
1861  << " at:" << hitInfo.rawPoint() << nl
1862  << " node:" << nodeI
1863  << " octant:" << octant
1864  << " bb:" << subBbox(nodeI, octant) << nl
1865  << " walking to neighbour containing:" << perturbedPoint
1866  << endl;
1867  }
1868 
1869 
1870  // Nothing hit so we are on face of bounding box (given as node+octant+
1871  // position bits). Traverse to neighbouring node. Use slightly perturbed
1872  // point.
1873 
1874  bool ok = walkToNeighbour
1875  (
1876  perturbedPoint,
1877  hitFaceID, // face(s) that hitInfo is on
1878 
1879  nodeI,
1880  octant
1881  );
1882 
1883  if (!ok)
1884  {
1885  // Hit the edge of the tree. Return miss.
1886  break;
1887  }
1888 
1889  if (verbose)
1890  {
1891  const treeBoundBox octantBb(subBbox(nodeI, octant));
1892  Pout<< " walked for point:" << hitInfo.rawPoint() << endl
1893  << " to neighbour node:" << nodeI
1894  << " octant:" << octant
1895  << " face:" << faceString(octantBb.faceBits(hitInfo.rawPoint()))
1896  << " of octantBb:" << octantBb << endl
1897  << endl;
1898  }
1899  }
1900 
1901  if (i == 100000)
1902  {
1903  // Probably in loop.
1904  if (!verbose)
1905  {
1906  // Redo intersection but now with verbose flag switched on.
1907  return findLine
1908  (
1909  findAny,
1910  treeStart,
1911  treeEnd,
1912  startNodeI,
1913  startOctant,
1914  true //verbose
1915  );
1916  }
1917  if (debug)
1918  {
1919  FatalErrorIn("indexedOctree<Type>::findLine(..)")
1920  << "Got stuck in loop raytracing from:" << treeStart
1921  << " to:" << treeEnd << endl
1922  << "inside top box:" << subBbox(startNodeI, startOctant)
1923  << abort(FatalError);
1924  }
1925  else
1926  {
1927  WarningIn("indexedOctree<Type>::findLine(..)")
1928  << "Got stuck in loop raytracing from:" << treeStart
1929  << " to:" << treeEnd << endl
1930  << "inside top box:" << subBbox(startNodeI, startOctant)
1931  << endl;
1932  }
1933  }
1934 
1935  return hitInfo;
1936 }
1937 
1938 
1939 // Find first intersection
1940 template <class Type>
1942 (
1943  const bool findAny,
1944  const point& start,
1945  const point& end
1946 ) const
1947 {
1948  pointIndexHit hitInfo;
1949 
1950  if (nodes_.size())
1951  {
1952  const treeBoundBox& treeBb = nodes_[0].bb_;
1953 
1954  // No effort is made to deal with points which are on edge of tree
1955  // bounding box for now.
1956 
1957  direction startBit = treeBb.posBits(start);
1958  direction endBit = treeBb.posBits(end);
1959 
1960  if ((startBit & endBit) != 0)
1961  {
1962  // Both start and end outside domain and in same block.
1963  return pointIndexHit(false, vector::zero, -1);
1964  }
1965 
1966 
1967  // Trim segment to treeBb
1968 
1969  point trackStart(start);
1970  point trackEnd(end);
1971 
1972  if (startBit != 0)
1973  {
1974  // Track start to inside domain.
1975  if (!treeBb.intersects(start, end, trackStart))
1976  {
1977  return pointIndexHit(false, vector::zero, -1);
1978  }
1979  }
1980 
1981  if (endBit != 0)
1982  {
1983  // Track end to inside domain.
1984  if (!treeBb.intersects(end, trackStart, trackEnd))
1985  {
1986  return pointIndexHit(false, vector::zero, -1);
1987  }
1988  }
1989 
1990 
1991  // Find lowest level tree node that start is in.
1992  labelBits index = findNode(0, trackStart);
1993 
1994  label parentNodeI = getNode(index);
1995  direction octant = getOctant(index);
1996 
1997  hitInfo = findLine
1998  (
1999  findAny,
2000  trackStart,
2001  trackEnd,
2002  parentNodeI,
2003  octant
2004  );
2005  }
2006 
2007  return hitInfo;
2008 }
2009 
2010 
2011 template <class Type>
2013 (
2014  const label nodeI,
2015  const treeBoundBox& searchBox,
2016  labelHashSet& elements
2017 ) const
2018 {
2019  const node& nod = nodes_[nodeI];
2020  const treeBoundBox& nodeBb = nod.bb_;
2021 
2022  for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
2023  {
2024  labelBits index = nod.subNodes_[octant];
2025 
2026  if (isNode(index))
2027  {
2028  const treeBoundBox& subBb = nodes_[getNode(index)].bb_;
2029 
2030  if (subBb.overlaps(searchBox))
2031  {
2032  findBox(getNode(index), searchBox, elements);
2033  }
2034  }
2035  else if (isContent(index))
2036  {
2037  const treeBoundBox subBb(nodeBb.subBbox(octant));
2038 
2039  if (subBb.overlaps(searchBox))
2040  {
2041  const labelList& indices = contents_[getContent(index)];
2042 
2043  forAll(indices, i)
2044  {
2045  label shapeI = indices[i];
2046 
2047  if (shapes_.overlaps(shapeI, searchBox))
2048  {
2049  elements.insert(shapeI);
2050  }
2051  }
2052  }
2053  }
2054  }
2055 }
2056 
2057 
2058 // Number of elements in node.
2059 template <class Type>
2061 (
2062  const labelBits index
2063 ) const
2064 {
2065  label nElems = 0;
2066 
2067  if (isNode(index))
2068  {
2069  // tree node.
2070  label nodeI = getNode(index);
2071 
2072  const node& nod = nodes_[nodeI];
2073 
2074  for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
2075  {
2076  nElems += countElements(nod.subNodes_[octant]);
2077  }
2078  }
2079  else if (isContent(index))
2080  {
2081  nElems += contents_[getContent(index)].size();
2082  }
2083  else
2084  {
2085  // empty node
2086  }
2087 
2088  return nElems;
2089 }
2090 
2091 
2092 template <class Type>
2094 (
2095  const label nodeI,
2096  const direction octant
2097 ) const
2098 {
2099  OFstream str
2100  (
2101  "node" + Foam::name(nodeI) + "_octant" + Foam::name(octant) + ".obj"
2102  );
2103 
2104  labelBits index = nodes_[nodeI].subNodes_[octant];
2105 
2106  treeBoundBox subBb;
2107 
2108  if (isNode(index))
2109  {
2110  subBb = nodes_[getNode(index)].bb_;
2111  }
2112  else if (isContent(index))
2113  {
2114  subBb = nodes_[nodeI].bb_.subBbox(octant);
2115  }
2116 
2117  Pout<< "dumpContentNode : writing node:" << nodeI << " octant:" << octant
2118  << " to " << str.name() << endl;
2119 
2120  label vertI = 0;
2121 
2122  // Dump bounding box
2123  pointField bbPoints(subBb.points());
2124 
2125  label pointVertI = vertI;
2126  forAll(bbPoints, i)
2127  {
2128  meshTools::writeOBJ(str, bbPoints[i]);
2129  vertI++;
2130  }
2131 
2132  forAll(treeBoundBox::edges, i)
2133  {
2134  const edge& e = treeBoundBox::edges[i];
2135 
2136  str<< "l " << e[0]+pointVertI+1 << ' ' << e[1]+pointVertI+1 << nl;
2137  }
2138 }
2139 
2140 
2141 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
2142 
2143 template <class Type>
2145 :
2146  shapes_(shapes),
2147  nodes_(0),
2148  contents_(0),
2149  nodeTypes_(0)
2150 {}
2151 
2152 
2153 template <class Type>
2156  const Type& shapes,
2157  const List<node>& nodes,
2158  const labelListList& contents
2159 )
2160 :
2161  shapes_(shapes),
2162  nodes_(nodes),
2163  contents_(contents),
2164  nodeTypes_(0)
2165 {}
2166 
2167 
2168 template <class Type>
2171  const Type& shapes,
2172  const treeBoundBox& bb,
2173  const label maxLevels, // maximum number of levels
2174  const scalar maxLeafRatio,
2175  const scalar maxDuplicity
2176 )
2177 :
2178  shapes_(shapes),
2179  nodes_(0),
2180  contents_(0),
2181  nodeTypes_(0)
2182 {
2183  if (debug)
2184  {
2185  Pout<< "indexedOctree<Type>::indexedOctree:" << nl
2186  << " shapes:" << shapes.size() << nl
2187  << " bb:" << bb << nl
2188  << endl;
2189  }
2190 
2191  if (shapes.size() == 0)
2192  {
2193  return;
2194  }
2195 
2196  // Start off with one node with all shapes in it.
2197  DynamicList<node> nodes(label(shapes.size() / maxLeafRatio));
2198  DynamicList<labelList> contents(label(shapes.size() / maxLeafRatio));
2199  contents.append(identity(shapes.size()));
2200 
2201  // Create topnode.
2202  node topNode(divide(bb, contents, 0));
2203  nodes.append(topNode);
2204 
2205 
2206  // Now have all contents at level 1. Create levels by splitting levels
2207  // above.
2208 
2209  label nLevels = 1;
2210 
2211  for (; nLevels < maxLevels; nLevels++)
2212  {
2213  // Count number of references into shapes (i.e. contents)
2214  label nEntries = 0;
2215  forAll(contents, i)
2216  {
2217  nEntries += contents[i].size();
2218  }
2219 
2220  if (debug)
2221  {
2222  Pout<< "indexedOctree<Type>::indexedOctree:" << nl
2223  << " nLevels:" << nLevels << nl
2224  << " nEntries per treeLeaf:" << nEntries/contents.size()
2225  << nl
2226  << " nEntries per shape (duplicity):"
2227  << nEntries/shapes.size()
2228  << nl
2229  << endl;
2230  }
2231 
2232  if
2233  (
2234  //nEntries < maxLeafRatio*contents.size()
2235  // ||
2236  nEntries > maxDuplicity*shapes.size()
2237  )
2238  {
2239  break;
2240  }
2241 
2242 
2243  // Split nodes with more than maxLeafRatio elements
2244  label nOldNodes = nodes.size();
2245  splitNodes
2246  (
2247  label(maxLeafRatio),
2248  nodes,
2249  contents
2250  );
2251 
2252  if (nOldNodes == nodes.size())
2253  {
2254  break;
2255  }
2256  }
2257 
2258  // Shrink
2259  nodes.shrink();
2260  contents.shrink();
2261 
2262 
2263  // Compact such that deeper level contents are always after the
2264  // ones for a shallower level. This way we can slice a coarser level
2265  // off the tree.
2266  contents_.setSize(contents.size());
2267  label compactI = 0;
2268 
2269  label level = 0;
2270 
2271  while (true)
2272  {
2273  compactContents
2274  (
2275  nodes,
2276  contents,
2277  level,
2278  0,
2279  0,
2280  contents_,
2281  compactI
2282  );
2283 
2284  if (compactI == contents_.size())
2285  {
2286  // Transferred all contents to contents_ (in order breadth first)
2287  break;
2288  }
2289 
2290  level++;
2291  }
2292  nodes_.transfer(nodes);
2293  nodes.clear();
2294 
2295  if (debug)
2296  {
2297  label nEntries = 0;
2298  forAll(contents_, i)
2299  {
2300  nEntries += contents_[i].size();
2301  }
2302 
2303  Pout<< "indexedOctree<Type>::indexedOctree"
2304  << " : finished construction of tree of:" << shapes.typeName
2305  << nl
2306  << " bb:" << this->bb() << nl
2307  << " shapes:" << shapes.size() << nl
2308  << " nLevels:" << nLevels << nl
2309  << " treeNodes:" << nodes_.size() << nl
2310  << " nEntries:" << nEntries << nl
2311  << " per treeLeaf:"
2312  << scalar(nEntries)/contents.size() << nl
2313  << " per shape (duplicity):"
2314  << scalar(nEntries)/shapes.size() << nl
2315  << endl;
2316  }
2317 }
2318 
2319 
2320 template <class Type>
2323  const Type& shapes,
2324  Istream& is
2325 )
2326 :
2327  shapes_(shapes),
2328  nodes_(is),
2329  contents_(is),
2330  nodeTypes_(0)
2331 {}
2332 
2333 
2334 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
2335 
2336 template <class Type>
2338 {
2339  return perturbTol_;
2340 }
2341 
2342 
2343 template <class Type>
2346  const point& sample,
2347  const scalar startDistSqr
2348 ) const
2349 {
2350  scalar nearestDistSqr = startDistSqr;
2351  label nearestShapeI = -1;
2352  point nearestPoint;
2353 
2354  if (nodes_.size())
2355  {
2356  findNearest
2357  (
2358  0,
2359  sample,
2360 
2361  nearestDistSqr,
2362  nearestShapeI,
2363  nearestPoint
2364  );
2365  }
2366  else
2367  {
2368  nearestPoint = vector::zero;
2369  }
2370 
2371  return pointIndexHit(nearestShapeI != -1, nearestPoint, nearestShapeI);
2372 }
2373 
2374 
2375 template <class Type>
2378  const linePointRef& ln,
2379  treeBoundBox& tightest,
2380  point& linePoint
2381 ) const
2382 {
2383  label nearestShapeI = -1;
2384  point nearestPoint;
2385 
2386  if (nodes_.size())
2387  {
2388  findNearest
2389  (
2390  0,
2391  ln,
2392 
2393  tightest,
2394  nearestShapeI,
2395  linePoint,
2396  nearestPoint
2397  );
2398  }
2399  else
2400  {
2401  nearestPoint = vector::zero;
2402  }
2403 
2404  return pointIndexHit(nearestShapeI != -1, nearestPoint, nearestShapeI);
2405 }
2406 
2407 
2408 // Find nearest intersection
2409 template <class Type>
2412  const point& start,
2413  const point& end
2414 ) const
2415 {
2416  return findLine(false, start, end);
2417 }
2418 
2419 
2420 // Find nearest intersection
2421 template <class Type>
2424  const point& start,
2425  const point& end
2426 ) const
2427 {
2428  return findLine(true, start, end);
2429 }
2430 
2431 
2432 template <class Type>
2435  const treeBoundBox& searchBox
2436 ) const
2437 {
2438  // Storage for labels of shapes inside bb. Size estimate.
2439  labelHashSet elements(shapes_.size() / 100);
2440 
2441  if (nodes_.size())
2442  {
2443  findBox(0, searchBox, elements);
2444  }
2445 
2446  return elements.toc();
2447 }
2448 
2449 
2450 // Find node (as parent+octant) containing point
2451 template <class Type>
2454  const label nodeI,
2455  const point& sample
2456 ) const
2457 {
2458  if (nodes_.empty())
2459  {
2460  // Empty tree. Return what?
2461  return nodePlusOctant(nodeI, 0);
2462  }
2463 
2464  const node& nod = nodes_[nodeI];
2465 
2466  if (debug)
2467  {
2468  if (!nod.bb_.contains(sample))
2469  {
2470  FatalErrorIn("findNode(..)")
2471  << "Cannot find " << sample << " in node " << nodeI
2472  << abort(FatalError);
2473  }
2474  }
2475 
2476  direction octant = nod.bb_.subOctant(sample);
2477 
2478  labelBits index = nod.subNodes_[octant];
2479 
2480  if (isNode(index))
2481  {
2482  // Recurse
2483  return findNode(getNode(index), sample);
2484  }
2485  else if (isContent(index))
2486  {
2487  // Content. Return treenode+octant
2488  return nodePlusOctant(nodeI, octant);
2489  }
2490  else
2491  {
2492  // Empty. Return treenode+octant
2493  return nodePlusOctant(nodeI, octant);
2494  }
2495 }
2496 
2497 
2498 // Determine type (inside/outside/mixed) per node.
2499 template <class Type>
2503  const point& sample
2504 ) const
2505 {
2506  if (nodes_.empty())
2507  {
2508  return UNKNOWN;
2509  }
2510 
2511  if (nodeTypes_.size() != 8*nodes_.size())
2512  {
2513  // Calculate type for every octant of node.
2514 
2515  nodeTypes_.setSize(8*nodes_.size());
2516  nodeTypes_ = UNKNOWN;
2517 
2518  calcVolumeType(0);
2519 
2520  if (debug)
2521  {
2522  label nUNKNOWN = 0;
2523  label nMIXED = 0;
2524  label nINSIDE = 0;
2525  label nOUTSIDE = 0;
2526 
2527  forAll(nodeTypes_, i)
2528  {
2529  volumeType type = volumeType(nodeTypes_.get(i));
2530 
2531  if (type == UNKNOWN)
2532  {
2533  nUNKNOWN++;
2534  }
2535  else if (type == MIXED)
2536  {
2537  nMIXED++;
2538  }
2539  else if (type == INSIDE)
2540  {
2541  nINSIDE++;
2542  }
2543  else if (type == OUTSIDE)
2544  {
2545  nOUTSIDE++;
2546  }
2547  else
2548  {
2549  FatalErrorIn("getVolumeType") << abort(FatalError);
2550  }
2551  }
2552 
2553  Pout<< "indexedOctree<Type>::getVolumeType : "
2554  << " bb:" << bb()
2555  << " nodes_:" << nodes_.size()
2556  << " nodeTypes_:" << nodeTypes_.size()
2557  << " nUNKNOWN:" << nUNKNOWN
2558  << " nMIXED:" << nMIXED
2559  << " nINSIDE:" << nINSIDE
2560  << " nOUTSIDE:" << nOUTSIDE
2561  << endl;
2562  }
2563  }
2564 
2565  return getVolumeType(0, sample);
2566 }
2567 
2568 
2569 // Print contents of nodeI
2570 template <class Type>
2573  prefixOSstream& os,
2574  const bool printContents,
2575  const label nodeI
2576 ) const
2577 {
2578  const node& nod = nodes_[nodeI];
2579  const treeBoundBox& bb = nod.bb_;
2580 
2581  os << "nodeI:" << nodeI << " bb:" << bb << nl
2582  << "parent:" << nod.parent_ << nl
2583  << "n:" << countElements(nodePlusOctant(nodeI, 0)) << nl;
2584 
2585  for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
2586  {
2587  const treeBoundBox subBb(bb.subBbox(octant));
2588 
2589  labelBits index = nod.subNodes_[octant];
2590 
2591  if (isNode(index))
2592  {
2593  // tree node.
2594  label subNodeI = getNode(index);
2595 
2596  os << "octant:" << octant
2597  << " node: n:" << countElements(index)
2598  << " bb:" << subBb << endl;
2599 
2600  string oldPrefix = os.prefix();
2601  os.prefix() = " " + oldPrefix;
2602 
2603  print(os, printContents, subNodeI);
2604 
2605  os.prefix() = oldPrefix;
2606  }
2607  else if (isContent(index))
2608  {
2609  const labelList& indices = contents_[getContent(index)];
2610 
2611  os << "octant:" << octant
2612  << " content: n:" << indices.size()
2613  << " bb:" << subBb;
2614 
2615  if (printContents)
2616  {
2617  os << " contents:";
2618  forAll(indices, j)
2619  {
2620  os << ' ' << indices[j];
2621  }
2622  }
2623  os << endl;
2624  }
2625  else
2626  {
2627  os << "octant:" << octant << " empty:" << subBb << endl;
2628  }
2629  }
2630 }
2631 
2632 
2633 // Print contents of nodeI
2634 template <class Type>
2636 {
2637  os << *this;
2638 
2639  return os.good();
2640 }
2641 
2642 
2643 template <class Type>
2644 Foam::Ostream& Foam::operator<<(Ostream& os, const indexedOctree<Type>& t)
2645 {
2646  return
2647  os << t.bb() << token::SPACE << t.nodes()
2648  << token::SPACE << t.contents();
2649 }
2650 
2651 
2652 // ************************ vim: set sw=4 sts=4 et: ************************ //