00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #include "octree.H"
00029 #include "treeLeaf.H"
00030 #include "treeNode.H"
00031 #include <OpenFOAM/long.H>
00032 #include <OSspecific/cpuTime.H>
00033 #include <OpenFOAM/linePointRef.H>
00034 #include "pointIndexHit.H"
00035
00036
00037
00038 template <class Type>
00039 Foam::string Foam::octree<Type>::volType(const label type)
00040 {
00041 if (type == UNKNOWN)
00042 {
00043 return "unknown";
00044 }
00045 else if (type == MIXED)
00046 {
00047 return "mixed";
00048 }
00049 else if (type == INSIDE)
00050 {
00051 return "inside";
00052 }
00053 else if (type == OUTSIDE)
00054 {
00055 return "outside";
00056 }
00057 else
00058 {
00059 FatalErrorIn("volType(const label)") << "type:" << type
00060 << " unknown." << abort(FatalError);
00061
00062 return "dummy";
00063 }
00064 }
00065
00066
00067
00068 template <class Type>
00069 Foam::label Foam::octree<Type>::getVolType
00070 (
00071 const vector& geomNormal,
00072 const vector& vec
00073 )
00074 {
00075 scalar sign = geomNormal & vec;
00076
00077 if (sign >= 0)
00078 {
00079 return octree<Type>::OUTSIDE;
00080 }
00081 else
00082 {
00083 return octree<Type>::INSIDE;
00084 }
00085 }
00086
00087
00088
00089
00090 template <class Type>
00091 Foam::octree<Type>::octree
00092 (
00093 const treeBoundBox& octreeBb,
00094 const Type& shapes,
00095 const label minNLevels,
00096 const scalar maxLeafRatio,
00097 const scalar maxShapeRatio
00098 )
00099 :
00100 topNode_(new treeNode<Type>(octreeBb)),
00101 shapes_(shapes),
00102 octreeBb_(octreeBb),
00103 maxLeafRatio_(maxLeafRatio),
00104 maxShapeRatio_(maxShapeRatio),
00105 minNLevels_(minNLevels),
00106 deepestLevel_(0),
00107 nEntries_(0),
00108 nNodes_(0),
00109 nLeaves_(0),
00110 endIter_(*this, -1),
00111 endConstIter_(*this, -1)
00112 {
00113 cpuTime timer;
00114
00115 setNodes(nNodes() + 1);
00116
00117 const label nShapes = shapes_.size();
00118
00119 labelList indices(nShapes);
00120 for (label i = 0; i < nShapes; i++)
00121 {
00122 indices[i] = i;
00123 }
00124
00125
00126 if (debug & 1)
00127 {
00128 Pout<< "octree : --- Start of Level " << deepestLevel_
00129 << " ----" << endl;
00130 }
00131 topNode_->distribute(0, *this, shapes_, indices);
00132
00133 if (debug & 1)
00134 {
00135 printStats(Pout);
00136 Pout<< "octree : --- End of Level " << deepestLevel_
00137 << " ----" << endl;
00138 }
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153 label oldNLeaves = -1;
00154 label oldNNodes = -1;
00155 deepestLevel_ = 1;
00156 while
00157 (
00158 (deepestLevel_ <= minNLevels_)
00159 || (
00160 (nEntries() > maxLeafRatio * nLeaves())
00161 && (nEntries() < maxShapeRatio * nShapes)
00162 )
00163 )
00164 {
00165 if (deepestLevel_ >= maxNLevels)
00166 {
00167 if (debug & 1)
00168 {
00169 Pout<< "octree : exiting since maxNLevels "
00170 << maxNLevels << " reached" << endl;
00171 }
00172 break;
00173 }
00174
00175 if ((oldNLeaves == nLeaves()) && (oldNNodes == nNodes()))
00176 {
00177 if (debug & 1)
00178 {
00179 Pout<< "octree : exiting since nLeaves and nNodes do not change"
00180 << endl;
00181 }
00182 break;
00183 }
00184 if (debug & 1)
00185 {
00186 Pout<< "octree : --- Start of Level " << deepestLevel_
00187 << " ----" << endl;
00188 }
00189
00190 oldNLeaves = nLeaves();
00191 oldNNodes = nNodes();
00192
00193 topNode_->redistribute
00194 (
00195 1,
00196 *this,
00197 shapes_,
00198 deepestLevel_
00199 );
00200
00201 if (debug & 1)
00202 {
00203 printStats(Pout);
00204
00205 Pout<< "octree : --- End of Level " << deepestLevel_
00206 << " ----" << endl;
00207 }
00208
00209 deepestLevel_++;
00210 }
00211
00212
00213 if (debug & 1)
00214 {
00215 Pout<< "octree : Constructed octree in = "
00216 << timer.cpuTimeIncrement()
00217 << " s\n" << endl << endl;
00218 }
00219
00220
00221 topNode_->setSubNodeType(0, *this, shapes_);
00222
00223 if (debug & 1)
00224 {
00225 Pout<< "octree : Added node information to octree in = "
00226 << timer.cpuTimeIncrement()
00227 << " s\n" << endl << endl;
00228 }
00229 }
00230
00231
00232
00233
00234 template <class Type>
00235 Foam::octree<Type>::~octree()
00236 {
00237 delete topNode_;
00238 }
00239
00240
00241
00242
00243 template <class Type>
00244 Foam::label Foam::octree<Type>::getSampleType(const point& sample) const
00245 {
00246 return topNode_->getSampleType(0, *this, shapes_, sample);
00247 }
00248
00249
00250 template <class Type>
00251 Foam::label Foam::octree<Type>::find(const point& sample) const
00252 {
00253 return topNode_->find(shapes_, sample);
00254 }
00255
00256
00257 template <class Type>
00258 bool Foam::octree<Type>::findTightest
00259 (
00260 const point& sample,
00261 treeBoundBox& tightest
00262 ) const
00263 {
00264 return topNode_->findTightest
00265 (
00266 shapes_,
00267 sample,
00268 tightest
00269 );
00270 }
00271
00272
00273 template <class Type>
00274 Foam::label Foam::octree<Type>::findNearest
00275 (
00276 const point& sample,
00277 treeBoundBox& tightest,
00278 scalar& tightestDist
00279 ) const
00280 {
00281 label tightestI = -1;
00282
00283 if (debug & 4)
00284 {
00285 Pout<< "octree::findNearest : searching for nearest for "
00286 << "sample:" << sample
00287 << " tightest:" << tightest << endl;
00288 }
00289
00290 topNode_->findNearest
00291 (
00292 shapes_,
00293 sample,
00294 tightest,
00295 tightestI,
00296 tightestDist
00297 );
00298
00299 if (debug & 4)
00300 {
00301 Pout<< "octree::findNearest : found nearest for "
00302 << "sample:" << sample << " with "
00303 << " tightestI:" << tightestI
00304 << " tightest:" << tightest
00305 << " tightestDist:" << tightestDist
00306 << endl;
00307 }
00308
00309 return tightestI;
00310 }
00311
00312
00313 template <class Type>
00314 Foam::label Foam::octree<Type>::findNearest
00315 (
00316 const linePointRef& ln,
00317 treeBoundBox& tightest,
00318 point& linePoint,
00319 point& shapePoint
00320 ) const
00321 {
00322
00323 label tightestI = -1;
00324 linePoint = point(-GREAT, -GREAT, -GREAT);
00325 shapePoint = point(GREAT, GREAT, GREAT);
00326
00327 topNode_->findNearest
00328 (
00329 shapes_,
00330 ln,
00331 tightest,
00332 tightestI,
00333 linePoint,
00334 shapePoint
00335 );
00336
00337 return tightestI;
00338 }
00339
00340
00341 template <class Type>
00342 Foam::labelList Foam::octree<Type>::findBox(const boundBox& bb) const
00343 {
00344
00345 labelHashSet elements(100);
00346
00347 topNode_->findBox(shapes_, bb, elements);
00348
00349 return elements.toc();
00350 }
00351
00352
00353 template <class Type>
00354 Foam::pointIndexHit Foam::octree<Type>::findLine
00355 (
00356 const point& treeStart,
00357 const point& treeEnd
00358 ) const
00359 {
00360
00361 pointIndexHit hitInfo(false, treeStart, -1);
00362
00363 const vector dir(treeEnd - treeStart);
00364
00365
00366 point start(treeStart);
00367 point end(treeEnd);
00368
00369 while (true)
00370 {
00371
00372 point leafIntPoint;
00373
00374 const treeLeaf<Type>* leafPtr = findLeafLine
00375 (
00376 start,
00377 end,
00378 leafIntPoint
00379 );
00380
00381 if (!leafPtr)
00382 {
00383
00384
00385 break;
00386 }
00387
00388
00389 scalar minS = GREAT;
00390
00391 const labelList& indices = leafPtr->indices();
00392
00393 forAll(indices, elemI)
00394 {
00395 label index = indices[elemI];
00396
00397 point pt;
00398 bool hit = shapes().intersects(index, start, end, pt);
00399
00400 if (hit)
00401 {
00402
00403 scalar s = (pt - treeStart) & dir;
00404
00405 if (s < minS)
00406 {
00407 minS = s;
00408
00409
00410 hitInfo.setHit();
00411 hitInfo.setPoint(pt);
00412 hitInfo.setIndex(index);
00413
00414
00415 end = pt;
00416 }
00417 }
00418 }
00419
00420 if (hitInfo.hit())
00421 {
00422
00423 break;
00424 }
00425
00426
00427 start = leafIntPoint;
00428 }
00429
00430 return hitInfo;
00431 }
00432
00433
00434 template <class Type>
00435 Foam::pointIndexHit Foam::octree<Type>::findLineAny
00436 (
00437 const point& start,
00438 const point& end
00439 ) const
00440 {
00441
00442 pointIndexHit hitInfo(false, start, -1);
00443
00444
00445 point p(start);
00446 while (true)
00447 {
00448
00449 point leafIntPoint;
00450
00451 const treeLeaf<Type>* leafPtr = findLeafLine(p, end, leafIntPoint);
00452
00453 if (!leafPtr)
00454 {
00455
00456
00457 break;
00458 }
00459
00460
00461
00462
00463 const vector edgeVec(end - p);
00464
00465 const labelList& indices = leafPtr->indices();
00466
00467 forAll(indices, elemI)
00468 {
00469 label index = indices[elemI];
00470
00471 point pt;
00472 bool hit = shapes().intersects
00473 (
00474 index,
00475 p,
00476 end,
00477 pt
00478 );
00479
00480 if (hit)
00481 {
00482 hitInfo.setHit();
00483 hitInfo.setPoint(pt);
00484 hitInfo.setIndex(index);
00485
00486 break;
00487 }
00488 }
00489
00490 if (hitInfo.hit())
00491 {
00492
00493 break;
00494 }
00495
00496
00497 p = leafIntPoint;
00498 }
00499
00500 return hitInfo;
00501 }
00502
00503
00504 template <class Type>
00505 const Foam::treeLeaf<Type>* Foam::octree<Type>::findLeafLine
00506 (
00507 const point& start,
00508 const point& end,
00509 point& leafIntPoint
00510 ) const
00511 {
00512
00513
00514 if (debug & 2)
00515 {
00516 Pout<< "octree::findLeafLine : searching for shapes on line "
00517 << "start:" << start
00518 << " end:" << end << endl;
00519 }
00520
00521
00522 if (octreeBb_.contains(start))
00523 {
00524 leafIntPoint = start;
00525 }
00526 else
00527 {
00528 if (!octreeBb_.intersects(start, end, leafIntPoint))
00529 {
00530 if (debug & 2)
00531 {
00532 Pout<< "octree::findLeafLine : start outside domain but does"
00533 << " not intersect : start:"
00534 << start << endl;
00535 }
00536 return NULL;
00537 }
00538
00539 if (debug & 2)
00540 {
00541 Pout<< "octree::findLeafLine : start propagated to inside"
00542 " domain : "
00543 << leafIntPoint << endl;
00544 }
00545 }
00546
00547
00548 const treeLeaf<Type>* leafPtr = topNode_->findLeafLine
00549 (
00550 0,
00551 shapes_,
00552 leafIntPoint,
00553 end
00554 );
00555
00556 if (debug & 2)
00557 {
00558 Pout<< "returning from octree::findLeafLine with "
00559 << "leafIntersection:" << leafIntPoint
00560 << " leafPtr:" << long(leafPtr) << endl;
00561 }
00562
00563 return leafPtr;
00564 }
00565
00566
00567 template <class Type>
00568 void Foam::octree<Type>::writeOBJ
00569 (
00570 Ostream& os,
00571 label& vertNo
00572 ) const
00573 {
00574 scalar minx = octreeBb_.min().x();
00575 scalar miny = octreeBb_.min().y();
00576 scalar minz = octreeBb_.min().z();
00577
00578 scalar maxx = octreeBb_.max().x();
00579 scalar maxy = octreeBb_.max().y();
00580 scalar maxz = octreeBb_.max().z();
00581
00582 os << "v " << minx << " " << miny << " " << minz << endl;
00583 os << "v " << maxx << " " << miny << " " << minz << endl;
00584 os << "v " << maxx << " " << maxy << " " << minz << endl;
00585 os << "v " << minx << " " << maxy << " " << minz << endl;
00586
00587 os << "v " << minx << " " << miny << " " << maxz << endl;
00588 os << "v " << maxx << " " << miny << " " << maxz << endl;
00589 os << "v " << maxx << " " << maxy << " " << maxz << endl;
00590 os << "v " << minx << " " << maxy << " " << maxz << endl;
00591
00592
00593 os << "l " << vertNo + 1 << " " << vertNo + 2 << endl;
00594 os << "l " << vertNo + 2 << " " << vertNo + 3 << endl;
00595 os << "l " << vertNo + 3 << " " << vertNo + 4 << endl;
00596 os << "l " << vertNo + 4 << " " << vertNo + 1 << endl;
00597
00598
00599 os << "l " << vertNo + 5 << " " << vertNo + 6 << endl;
00600 os << "l " << vertNo + 6 << " " << vertNo + 7 << endl;
00601 os << "l " << vertNo + 7 << " " << vertNo + 8 << endl;
00602 os << "l " << vertNo + 8 << " " << vertNo + 5 << endl;
00603
00604
00605 os << "l " << vertNo + 1 << " " << vertNo + 5 << endl;
00606 os << "l " << vertNo + 2 << " " << vertNo + 6 << endl;
00607 os << "l " << vertNo + 3 << " " << vertNo + 7 << endl;
00608 os << "l " << vertNo + 4 << " " << vertNo + 8 << endl;
00609
00610 vertNo += 8;
00611
00612 topNode_->writeOBJ(os, 1, vertNo);
00613 }
00614
00615
00616 template <class Type>
00617 void Foam::octree<Type>::printStats(Ostream& os) const
00618 {
00619 os << "Statistics after iteration " << deepestLevel() << ':' << endl
00620 << " nShapes :" << shapes().size() << endl
00621 << " nNodes :" << nNodes() << endl
00622 << " nLeaves :" << nLeaves() << endl
00623 << " nEntries :" << nEntries() << endl;
00624
00625 if (nLeaves() && shapes().size())
00626 {
00627 os
00628 << " Cells per leaf :"
00629 << scalar(nEntries())/nLeaves()
00630 << nl
00631 << " Every cell in :"
00632 << scalar(nEntries())/shapes().size() << " cubes"
00633 << endl;
00634 }
00635 }
00636
00637
00638
00639
00640
00641 template <class Type>
00642 Foam::octree<Type>::iterator::iterator(octree<Type>& oc)
00643 :
00644 octree_(oc),
00645 curLeaf_(oc.nLeaves())
00646 {
00647 leaves_.setSize(0);
00648 }
00649
00650
00651
00652 template <class Type>
00653 Foam::octree<Type>::iterator::iterator(octree<Type>& oc, label index)
00654 :
00655 octree_(oc),
00656 curLeaf_(index)
00657 {
00658 if (index >= 0)
00659 {
00660 leaves_.setSize(oc.nLeaves());
00661
00662 label leafIndex = 0;
00663 octree_.topNode()->findLeaves(leaves_, leafIndex);
00664
00665 if (leafIndex != oc.nLeaves())
00666 {
00667 FatalErrorIn
00668 (
00669 "octree::iterator::iterator"
00670 "(octree&, label)"
00671 )
00672 << "Traversal of tree returns : " << leafIndex << endl
00673 << "Statistics of octree say : " << oc.nLeaves() << endl
00674 << abort(FatalError);
00675 }
00676 }
00677 }
00678
00679
00680 template <class Type>
00681 void Foam::octree<Type>::iterator::operator=(const iterator& iter)
00682 {
00683 if ((curLeaf_ < 0) && (iter.curLeaf_ >= 0))
00684 {
00685 FatalErrorIn
00686 (
00687 "octree::iterator::operator="
00688 "(const iterator&)"
00689 )
00690 << "lhs : " << curLeaf_ << endl
00691 << "rhs : " << iter.curLeaf_ << endl
00692 << abort(FatalError);
00693 }
00694 curLeaf_ = iter.curLeaf_;
00695 }
00696
00697
00698 template <class Type>
00699 bool Foam::octree<Type>::iterator::operator==(const iterator& iter) const
00700 {
00701 label index1 =
00702 (curLeaf_ >= 0 ? curLeaf_ : octree_.nLeaves());
00703 label index2 =
00704 (iter.curLeaf_ >= 0 ? iter.curLeaf_ : iter.octree_.nLeaves());
00705
00706 return index1 == index2;
00707 }
00708
00709
00710 template <class Type>
00711 bool Foam::octree<Type>::iterator::operator!=(const iterator& iter) const
00712 {
00713 return !(iterator::operator==(iter));
00714 }
00715
00716
00717 template <class Type>
00718 Foam::treeLeaf<Type>& Foam::octree<Type>::iterator::operator*()
00719 {
00720 return *leaves_[curLeaf_];
00721 }
00722
00723
00724 template <class Type>
00725 typename Foam::octree<Type>::iterator&
00726 Foam::octree<Type>::iterator::operator++()
00727 {
00728 curLeaf_++;
00729 return *this;
00730 }
00731
00732
00733 template <class Type>
00734 typename Foam::octree<Type>::iterator
00735 Foam::octree<Type>::iterator::operator++(int)
00736 {
00737 iterator tmp = *this;
00738 ++*this;
00739 return tmp;
00740 }
00741
00742
00743 template <class Type>
00744 typename Foam::octree<Type>::iterator
00745 Foam::octree<Type>::begin()
00746 {
00747 return iterator(*this, 0);
00748 }
00749
00750
00751 template <class Type>
00752 const typename Foam::octree<Type>::iterator&
00753 Foam::octree<Type>::end()
00754 {
00755 return octree<Type>::endIter_;
00756 }
00757
00758
00759
00760
00761
00762 template <class Type>
00763 Foam::octree<Type>::const_iterator::const_iterator(const octree<Type>& oc)
00764 :
00765 octree_(oc),
00766 curLeaf_(oc.nLeaves())
00767 {
00768 leaves_.setSize(oc.nLeaves());
00769 }
00770
00771
00772
00773 template <class Type>
00774 Foam::octree<Type>::const_iterator::const_iterator
00775 (
00776 const octree<Type>& oc,
00777 label index
00778 )
00779 :
00780 octree_(oc),
00781 curLeaf_(index)
00782 {
00783 if (index >= 0)
00784 {
00785 leaves_.setSize(oc.nLeaves());
00786
00787 label leafIndex = 0;
00788 octree_.topNode()->findLeaves(leaves_, leafIndex);
00789
00790 if (leafIndex != oc.nLeaves())
00791 {
00792 FatalErrorIn
00793 (
00794 "octree::const_iterator::const_iterator"
00795 "(octree&, label)"
00796 )
00797 << "Traversal of tree returns : " << leafIndex << endl
00798 << "Statistics of octree say : " << oc.nLeaves() << endl
00799 << abort(FatalError);
00800 }
00801 }
00802 }
00803
00804
00805 template <class Type>
00806 void Foam::octree<Type>::const_iterator::operator=(const const_iterator& iter)
00807 {
00808 if ((curLeaf_ < 0) && (iter.curLeaf_ >= 0))
00809 {
00810 FatalErrorIn
00811 (
00812 "octree::const_iterator::operator="
00813 "(const const_iterator&)"
00814 )
00815 << "lhs : " << curLeaf_ << endl
00816 << "rhs : " << iter.curLeaf_ << endl
00817 << abort(FatalError);
00818 }
00819 curLeaf_ = iter.curLeaf_;
00820 curLeaf_ = iter.curLeaf_;
00821 }
00822
00823
00824 template <class Type>
00825 bool Foam::octree<Type>::const_iterator::operator==
00826 (
00827 const const_iterator& iter
00828 ) const
00829 {
00830 label index1 =
00831 (curLeaf_ >= 0 ? curLeaf_ : octree_.nLeaves());
00832 label index2 =
00833 (iter.curLeaf_ >= 0 ? iter.curLeaf_ : iter.octree_.nLeaves());
00834
00835 return index1 == index2;
00836 }
00837
00838
00839 template <class Type>
00840 bool Foam::octree<Type>::const_iterator::operator!=
00841 (
00842 const const_iterator& iter
00843 ) const
00844 {
00845 return !(const_iterator::operator==(iter));
00846 }
00847
00848
00849 template <class Type>
00850 const Foam::treeLeaf<Type>& Foam::octree<Type>::const_iterator::operator*()
00851 {
00852 return *leaves_[curLeaf_];
00853 }
00854
00855
00856 template <class Type>
00857 typename Foam::octree<Type>::const_iterator&
00858 Foam::octree<Type>::const_iterator::operator++()
00859 {
00860 curLeaf_++;
00861 return *this;
00862 }
00863
00864
00865 template <class Type>
00866 typename Foam::octree<Type>::const_iterator
00867 Foam::octree<Type>::const_iterator::operator++(int)
00868 {
00869 const_iterator tmp = *this;
00870 ++*this;
00871 return tmp;
00872 }
00873
00874
00875 template <class Type>
00876 typename Foam::octree<Type>::const_iterator
00877 Foam::octree<Type>::begin() const
00878 {
00879 return const_iterator(*this, 0);
00880 }
00881
00882
00883 template <class Type>
00884 typename Foam::octree<Type>::const_iterator
00885 Foam::octree<Type>::cbegin() const
00886 {
00887 return const_iterator(*this, 0);
00888 }
00889
00890
00891 template <class Type>
00892 const typename Foam::octree<Type>::const_iterator&
00893 Foam::octree<Type>::end() const
00894 {
00895 return octree<Type>::endConstIter_;
00896 }
00897
00898
00899 template <class Type>
00900 const typename Foam::octree<Type>::const_iterator&
00901 Foam::octree<Type>::cend() const
00902 {
00903 return octree<Type>::endConstIter_;
00904 }
00905
00906
00907
00908
00909 template <class Type>
00910 Foam::Ostream& Foam::operator<<(Ostream& os, const octree<Type>& oc)
00911 {
00912 return os << token::BEGIN_LIST
00913
00914 << token::SPACE << oc.octreeBb_
00915 << token::SPACE << oc.maxLeafRatio_
00916 << token::SPACE << oc.maxShapeRatio_
00917 << token::SPACE << oc.minNLevels_
00918 << token::SPACE << oc.deepestLevel_
00919 << token::SPACE << oc.nEntries_
00920 << token::SPACE << oc.nNodes_
00921 << token::SPACE << oc.nLeaves_
00922 << token::SPACE << *oc.topNode_
00923 << token::SPACE << token::END_LIST;
00924 }
00925
00926
00927