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