FreeFOAM The Cross-Platform CFD Toolkit
octree.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2011 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 Description
25 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "octree.H"
29 #include "treeLeaf.H"
30 #include "treeNode.H"
31 #include <OpenFOAM/long.H>
32 #include <OSspecific/cpuTime.H>
33 #include <OpenFOAM/linePointRef.H>
34 #include "pointIndexHit.H"
35 
36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37 
38 template <class Type>
40 {
41  if (type == UNKNOWN)
42  {
43  return "unknown";
44  }
45  else if (type == MIXED)
46  {
47  return "mixed";
48  }
49  else if (type == INSIDE)
50  {
51  return "inside";
52  }
53  else if (type == OUTSIDE)
54  {
55  return "outside";
56  }
57  else
58  {
59  FatalErrorIn("volType(const label)") << "type:" << type
60  << " unknown." << abort(FatalError);
61 
62  return "dummy";
63  }
64 }
65 
66 
67 // Determine inside/outside status of vector compared to geometry-based normal
68 template <class Type>
70 (
71  const vector& geomNormal,
72  const vector& vec
73 )
74 {
75  scalar sign = geomNormal & vec;
76 
77  if (sign >= 0)
78  {
79  return octree<Type>::OUTSIDE;
80  }
81  else
82  {
83  return octree<Type>::INSIDE;
84  }
85 }
86 
87 
88 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
89 
90 template <class Type>
92 (
93  const treeBoundBox& octreeBb,
94  const Type& shapes,
95  const label minNLevels,
96  const scalar maxLeafRatio,
97  const scalar maxShapeRatio
98 )
99 :
100  topNode_(new treeNode<Type>(octreeBb)),
101  shapes_(shapes),
102  octreeBb_(octreeBb),
103  maxLeafRatio_(maxLeafRatio),
104  maxShapeRatio_(maxShapeRatio),
105  minNLevels_(minNLevels),
106  deepestLevel_(0),
107  nEntries_(0),
108  nNodes_(0),
109  nLeaves_(0),
110  endIter_(*this, -1),
111  endConstIter_(*this, -1)
112 {
113  cpuTime timer;
114 
115  setNodes(nNodes() + 1);
116 
117  const label nShapes = shapes_.size();
118 
119  labelList indices(nShapes);
120  for (label i = 0; i < nShapes; i++)
121  {
122  indices[i] = i;
123  }
124 
125  // Create initial level (0) of subLeaves
126  if (debug & 1)
127  {
128  Pout<< "octree : --- Start of Level " << deepestLevel_
129  << " ----" << endl;
130  }
131  topNode_->distribute(0, *this, shapes_, indices);
132 
133  if (debug & 1)
134  {
135  printStats(Pout);
136  Pout<< "octree : --- End of Level " << deepestLevel_
137  << " ----" << endl;
138  }
139 
140  // Breadth first creation of tree
141  // Stop if: - level above minlevel and
142  // - less than so many cells per endpoint
143  // (so bottom level is fine enough)
144  // - every shape mentioned in only so many
145  // leaves. (if shape bb quite big it will end up
146  // in lots of leaves).
147  // - no change in number of leaves
148  // (happens if leafs fine enough)
149  // This guarantees that tree
150  // - is fine enough (if minLevel > 0)
151  // - has some guaranteed maximum size (maxShapeRatio)
152 
153  label oldNLeaves = -1; // make test below pass first time.
154  label oldNNodes = -1;
155  deepestLevel_ = 1;
156  while
157  (
158  (deepestLevel_ <= minNLevels_)
159  || (
160  (nEntries() > maxLeafRatio * nLeaves()) // shapes per leaf
161  && (nEntries() < maxShapeRatio * nShapes) // entries per shape
162  )
163  )
164  {
165  if (deepestLevel_ >= maxNLevels)
166  {
167  if (debug & 1)
168  {
169  Pout<< "octree : exiting since maxNLevels "
170  << maxNLevels << " reached" << endl;
171  }
172  break;
173  }
174 
175  if ((oldNLeaves == nLeaves()) && (oldNNodes == nNodes()))
176  {
177  if (debug & 1)
178  {
179  Pout<< "octree : exiting since nLeaves and nNodes do not change"
180  << endl;
181  }
182  break;
183  }
184  if (debug & 1)
185  {
186  Pout<< "octree : --- Start of Level " << deepestLevel_
187  << " ----" << endl;
188  }
189 
190  oldNLeaves = nLeaves();
191  oldNNodes = nNodes();
192 
193  topNode_->redistribute
194  (
195  1,
196  *this,
197  shapes_,
198  deepestLevel_
199  );
200 
201  if (debug & 1)
202  {
203  printStats(Pout);
204 
205  Pout<< "octree : --- End of Level " << deepestLevel_
206  << " ----" << endl;
207  }
208 
209  deepestLevel_++;
210  }
211 
212 
213  if (debug & 1)
214  {
215  Pout<< "octree : Constructed octree in = "
216  << timer.cpuTimeIncrement()
217  << " s\n" << endl << endl;
218  }
219 
220  // Set volume type of non-treeleaf nodes.
221  topNode_->setSubNodeType(0, *this, shapes_);
222 
223  if (debug & 1)
224  {
225  Pout<< "octree : Added node information to octree in = "
226  << timer.cpuTimeIncrement()
227  << " s\n" << endl << endl;
228  }
229 }
230 
231 
232 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
233 
234 template <class Type>
236 {
237  delete topNode_;
238 }
239 
240 
241 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
242 
243 template <class Type>
244 Foam::label Foam::octree<Type>::getSampleType(const point& sample) const
245 {
246  return topNode_->getSampleType(0, *this, shapes_, sample);
247 }
248 
249 
250 template <class Type>
251 Foam::label Foam::octree<Type>::find(const point& sample) const
252 {
253  return topNode_->find(shapes_, sample);
254 }
255 
256 
257 template <class Type>
259 (
260  const point& sample,
261  treeBoundBox& tightest
262 ) const
263 {
264  return topNode_->findTightest
265  (
266  shapes_,
267  sample,
268  tightest
269  );
270 }
271 
272 
273 template <class Type>
275 (
276  const point& sample,
277  treeBoundBox& tightest,
278  scalar& tightestDist
279 ) const
280 {
281  label tightestI = -1;
282 
283  if (debug & 4)
284  {
285  Pout<< "octree::findNearest : searching for nearest for "
286  << "sample:" << sample
287  << " tightest:" << tightest << endl;
288  }
289 
290  topNode_->findNearest
291  (
292  shapes_,
293  sample,
294  tightest,
295  tightestI,
296  tightestDist
297  );
298 
299  if (debug & 4)
300  {
301  Pout<< "octree::findNearest : found nearest for "
302  << "sample:" << sample << " with "
303  << " tightestI:" << tightestI
304  << " tightest:" << tightest
305  << " tightestDist:" << tightestDist
306  << endl;
307  }
308 
309  return tightestI;
310 }
311 
312 
313 template <class Type>
315 (
316  const linePointRef& ln,
317  treeBoundBox& tightest,
318  point& linePoint,
319  point& shapePoint
320 ) const
321 {
322  // Start off from miss with points at large distance apart.
323  label tightestI = -1;
324  linePoint = point(-GREAT, -GREAT, -GREAT);
325  shapePoint = point(GREAT, GREAT, GREAT);
326 
327  topNode_->findNearest
328  (
329  shapes_,
330  ln,
331  tightest,
332  tightestI,
333  linePoint,
334  shapePoint
335  );
336 
337  return tightestI;
338 }
339 
340 
341 template <class Type>
343 {
344  // Storage for labels of shapes inside bb. Size estimate.
345  labelHashSet elements(100);
346 
347  topNode_->findBox(shapes_, bb, elements);
348 
349  return elements.toc();
350 }
351 
352 
353 template <class Type>
355 (
356  const point& treeStart,
357  const point& treeEnd
358 ) const
359 {
360  // Initialize to a miss
361  pointIndexHit hitInfo(false, treeStart, -1);
362 
363  const vector dir(treeEnd - treeStart);
364 
365  // Current line segment to search
366  point start(treeStart);
367  point end(treeEnd);
368 
369  while (true)
370  {
371  // Find nearest treeLeaf intersected by line
372  point leafIntPoint;
373 
374  const treeLeaf<Type>* leafPtr = findLeafLine
375  (
376  start,
377  end,
378  leafIntPoint
379  );
380 
381  if (!leafPtr)
382  {
383  // Reached end of string of treeLeaves to be searched. Return
384  // whatever we have found so far.
385  break;
386  }
387 
388  // Inside treeLeaf find nearest intersection
389  scalar minS = GREAT;
390 
391  const labelList& indices = leafPtr->indices();
392 
393  forAll(indices, elemI)
394  {
395  label index = indices[elemI];
396 
397  point pt;
398  bool hit = shapes().intersects(index, start, end, pt);
399 
400  if (hit)
401  {
402  // Check whether intersection nearer p
403  scalar s = (pt - treeStart) & dir;
404 
405  if (s < minS)
406  {
407  minS = s;
408 
409  // Update hit info
410  hitInfo.setHit();
411  hitInfo.setPoint(pt);
412  hitInfo.setIndex(index);
413 
414  // Update segment to search
415  end = pt;
416  }
417  }
418  }
419 
420  if (hitInfo.hit())
421  {
422  // Found intersected shape.
423  break;
424  }
425 
426  // Start from end of current leaf
427  start = leafIntPoint;
428  }
429 
430  return hitInfo;
431 }
432 
433 
434 template <class Type>
436 (
437  const point& start,
438  const point& end
439 ) const
440 {
441  // Initialize to a miss
442  pointIndexHit hitInfo(false, start, -1);
443 
444  // Start of segment in current treeNode.
445  point p(start);
446  while (true)
447  {
448  // Find treeLeaf intersected by line
449  point leafIntPoint;
450 
451  const treeLeaf<Type>* leafPtr = findLeafLine(p, end, leafIntPoint);
452 
453  if (!leafPtr)
454  {
455  // Reached end of string of treeLeaves to be searched. Return
456  // whatever we have found so far.
457  break;
458  }
459 
460  // Inside treeLeaf find any intersection
461 
462  // Vector from endPoint to entry point of leaf.
463  const vector edgeVec(end - p);
464 
465  const labelList& indices = leafPtr->indices();
466 
467  forAll(indices, elemI)
468  {
469  label index = indices[elemI];
470 
471  point pt;
472  bool hit = shapes().intersects
473  (
474  index,
475  p,
476  end,
477  pt
478  );
479 
480  if (hit)
481  {
482  hitInfo.setHit();
483  hitInfo.setPoint(pt);
484  hitInfo.setIndex(index);
485 
486  break;
487  }
488  }
489 
490  if (hitInfo.hit())
491  {
492  // Found intersected shape.
493  break;
494  }
495 
496  // Start from end of current leaf
497  p = leafIntPoint;
498  }
499 
500  return hitInfo;
501 }
502 
503 
504 template <class Type>
506 (
507  const point& start,
508  const point& end,
509  point& leafIntPoint
510 ) const
511 {
512  // returns first found cube along line
513 
514  if (debug & 2)
515  {
516  Pout<< "octree::findLeafLine : searching for shapes on line "
517  << "start:" << start
518  << " end:" << end << endl;
519  }
520 
521  // If start is outside project onto top cube
522  if (octreeBb_.contains(start))
523  {
524  leafIntPoint = start;
525  }
526  else
527  {
528  if (!octreeBb_.intersects(start, end, leafIntPoint))
529  {
530  if (debug & 2)
531  {
532  Pout<< "octree::findLeafLine : start outside domain but does"
533  << " not intersect : start:"
534  << start << endl;
535  }
536  return NULL;
537  }
538 
539  if (debug & 2)
540  {
541  Pout<< "octree::findLeafLine : start propagated to inside"
542  " domain : "
543  << leafIntPoint << endl;
544  }
545  }
546 
547  // Normal action: find next intersection along line
548  const treeLeaf<Type>* leafPtr = topNode_->findLeafLine
549  (
550  0,
551  shapes_,
552  leafIntPoint,
553  end
554  );
555 
556  if (debug & 2)
557  {
558  Pout<< "returning from octree::findLeafLine with "
559  << "leafIntersection:" << leafIntPoint
560  << " leafPtr:" << long(leafPtr) << endl;
561  }
562 
563  return leafPtr;
564 }
565 
566 
567 template <class Type>
569 (
570  Ostream& os,
571  label& vertNo
572 ) const
573 {
574  scalar minx = octreeBb_.min().x();
575  scalar miny = octreeBb_.min().y();
576  scalar minz = octreeBb_.min().z();
577 
578  scalar maxx = octreeBb_.max().x();
579  scalar maxy = octreeBb_.max().y();
580  scalar maxz = octreeBb_.max().z();
581 
582  os << "v " << minx << " " << miny << " " << minz << endl;
583  os << "v " << maxx << " " << miny << " " << minz << endl;
584  os << "v " << maxx << " " << maxy << " " << minz << endl;
585  os << "v " << minx << " " << maxy << " " << minz << endl;
586 
587  os << "v " << minx << " " << miny << " " << maxz << endl;
588  os << "v " << maxx << " " << miny << " " << maxz << endl;
589  os << "v " << maxx << " " << maxy << " " << maxz << endl;
590  os << "v " << minx << " " << maxy << " " << maxz << endl;
591 
592  // Bottom face
593  os << "l " << vertNo + 1 << " " << vertNo + 2 << endl;
594  os << "l " << vertNo + 2 << " " << vertNo + 3 << endl;
595  os << "l " << vertNo + 3 << " " << vertNo + 4 << endl;
596  os << "l " << vertNo + 4 << " " << vertNo + 1 << endl;
597 
598  // Top face
599  os << "l " << vertNo + 5 << " " << vertNo + 6 << endl;
600  os << "l " << vertNo + 6 << " " << vertNo + 7 << endl;
601  os << "l " << vertNo + 7 << " " << vertNo + 8 << endl;
602  os << "l " << vertNo + 8 << " " << vertNo + 5 << endl;
603 
604  // Edges from bottom to top face
605  os << "l " << vertNo + 1 << " " << vertNo + 5 << endl;
606  os << "l " << vertNo + 2 << " " << vertNo + 6 << endl;
607  os << "l " << vertNo + 3 << " " << vertNo + 7 << endl;
608  os << "l " << vertNo + 4 << " " << vertNo + 8 << endl;
609 
610  vertNo += 8;
611 
612  topNode_->writeOBJ(os, 1, vertNo);
613 }
614 
615 
616 template <class Type>
618 {
619  os << "Statistics after iteration " << deepestLevel() << ':' << endl
620  << " nShapes :" << shapes().size() << endl
621  << " nNodes :" << nNodes() << endl
622  << " nLeaves :" << nLeaves() << endl
623  << " nEntries :" << nEntries() << endl;
624 
625  if (nLeaves() && shapes().size())
626  {
627  os
628  << " Cells per leaf :"
629  << scalar(nEntries())/nLeaves()
630  << nl
631  << " Every cell in :"
632  << scalar(nEntries())/shapes().size() << " cubes"
633  << endl;
634  }
635 }
636 
637 
638 // * * * * * * * * * * * * * * * STL iterator * * * * * * * * * * * * * * * //
639 
640 // Construct from a octree. Set index at end.
641 template <class Type>
643 :
644  octree_(oc),
645  curLeaf_(oc.nLeaves())
646 {
647  leaves_.setSize(0);
648 }
649 
650 
651 // Construct from octree. Set index.
652 template <class Type>
653 Foam::octree<Type>::iterator::iterator(octree<Type>& oc, label index)
654 :
655  octree_(oc),
656  curLeaf_(index)
657 {
658  if (index >= 0)
659  {
660  leaves_.setSize(oc.nLeaves());
661 
662  label leafIndex = 0;
663  octree_.topNode()->findLeaves(leaves_, leafIndex);
664 
665  if (leafIndex != oc.nLeaves())
666  {
668  (
669  "octree::iterator::iterator"
670  "(octree&, label)"
671  )
672  << "Traversal of tree returns : " << leafIndex << endl
673  << "Statistics of octree say : " << oc.nLeaves() << endl
674  << abort(FatalError);
675  }
676  }
677 }
678 
679 
680 template <class Type>
682 {
683  if ((curLeaf_ < 0) && (iter.curLeaf_ >= 0))
684  {
686  (
687  "octree::iterator::operator="
688  "(const iterator&)"
689  )
690  << "lhs : " << curLeaf_ << endl
691  << "rhs : " << iter.curLeaf_ << endl
692  << abort(FatalError);
693  }
694  curLeaf_ = iter.curLeaf_;
695 }
696 
697 
698 template <class Type>
700 {
701  label index1 =
702  (curLeaf_ >= 0 ? curLeaf_ : octree_.nLeaves());
703  label index2 =
704  (iter.curLeaf_ >= 0 ? iter.curLeaf_ : iter.octree_.nLeaves());
705 
706  return index1 == index2;
707 }
708 
709 
710 template <class Type>
712 {
713  return !(iterator::operator==(iter));
714 }
715 
716 
717 template <class Type>
719 {
720  return *leaves_[curLeaf_];
721 }
722 
723 
724 template <class Type>
727 {
728  curLeaf_++;
729  return *this;
730 }
731 
732 
733 template <class Type>
736 {
737  iterator tmp = *this;
738  ++*this;
739  return tmp;
740 }
741 
742 
743 template <class Type>
746 {
747  return iterator(*this, 0);
748 }
749 
750 
751 template <class Type>
752 const typename Foam::octree<Type>::iterator&
754 {
755  return octree<Type>::endIter_;
756 }
757 
758 
759 // * * * * * * * * * * * * * * STL const_iterator * * * * * * * * * * * * * //
760 
761 // Construct for a given octree
762 template <class Type>
764 :
765  octree_(oc),
766  curLeaf_(oc.nLeaves())
767 {
768  leaves_.setSize(oc.nLeaves());
769 }
770 
771 
772 // Construct for a given octree
773 template <class Type>
775 (
776  const octree<Type>& oc,
777  label index
778 )
779 :
780  octree_(oc),
781  curLeaf_(index)
782 {
783  if (index >= 0)
784  {
785  leaves_.setSize(oc.nLeaves());
786 
787  label leafIndex = 0;
788  octree_.topNode()->findLeaves(leaves_, leafIndex);
789 
790  if (leafIndex != oc.nLeaves())
791  {
793  (
794  "octree::const_iterator::const_iterator"
795  "(octree&, label)"
796  )
797  << "Traversal of tree returns : " << leafIndex << endl
798  << "Statistics of octree say : " << oc.nLeaves() << endl
799  << abort(FatalError);
800  }
801  }
802 }
803 
804 
805 template <class Type>
807 {
808  if ((curLeaf_ < 0) && (iter.curLeaf_ >= 0))
809  {
811  (
812  "octree::const_iterator::operator="
813  "(const const_iterator&)"
814  )
815  << "lhs : " << curLeaf_ << endl
816  << "rhs : " << iter.curLeaf_ << endl
817  << abort(FatalError);
818  }
819  curLeaf_ = iter.curLeaf_;
820  curLeaf_ = iter.curLeaf_;
821 }
822 
823 
824 template <class Type>
826 (
827  const const_iterator& iter
828 ) const
829 {
830  label index1 =
831  (curLeaf_ >= 0 ? curLeaf_ : octree_.nLeaves());
832  label index2 =
833  (iter.curLeaf_ >= 0 ? iter.curLeaf_ : iter.octree_.nLeaves());
834 
835  return index1 == index2;
836 }
837 
838 
839 template <class Type>
840 bool Foam::octree<Type>::const_iterator::operator!=
841 (
842  const const_iterator& iter
843 ) const
844 {
845  return !(const_iterator::operator==(iter));
846 }
847 
848 
849 template <class Type>
851 {
852  return *leaves_[curLeaf_];
853 }
854 
855 
856 template <class Type>
859 {
860  curLeaf_++;
861  return *this;
862 }
863 
864 
865 template <class Type>
868 {
869  const_iterator tmp = *this;
870  ++*this;
871  return tmp;
872 }
873 
874 
875 template <class Type>
878 {
879  return const_iterator(*this, 0);
880 }
881 
882 
883 template <class Type>
886 {
887  return const_iterator(*this, 0);
888 }
889 
890 
891 template <class Type>
894 {
896 }
897 
898 
899 template <class Type>
902 {
904 }
905 
906 
907 // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
908 
909 template <class Type>
910 Foam::Ostream& Foam::operator<<(Ostream& os, const octree<Type>& oc)
911 {
912  return os << token::BEGIN_LIST
913  //<< token::SPACE << oc.shapes_
914  << token::SPACE << oc.octreeBb_
915  << token::SPACE << oc.maxLeafRatio_
916  << token::SPACE << oc.maxShapeRatio_
917  << token::SPACE << oc.minNLevels_
918  << token::SPACE << oc.deepestLevel_
919  << token::SPACE << oc.nEntries_
920  << token::SPACE << oc.nNodes_
921  << token::SPACE << oc.nLeaves_
922  << token::SPACE << *oc.topNode_
924 }
925 
926 
927 // ************************ vim: set sw=4 sts=4 et: ************************ //