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