src/meshTools/octree/treeNode.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 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 "treeNode.H"
  27 #include "octree.H"
  28 #include "treeLeaf.H"
  29 #include "treeBoundBox.H"
  30 #include "long.H"
  31 #include "linePointRef.H"
  32
  33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  34
  35 template <class Type>
  36 const Foam::label Foam::treeNode<Type>::leafOffset(100);
  37
  38
  39 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  40
  41 template <class Type>
  42 void Foam::treeNode<Type>::setAsNode(const label octant)
  43 {
  44     subNodeTypes_ |= (0x1 << octant);
  45 }
  46
  47
  48 template <class Type>
  49 void Foam::treeNode<Type>::setAsLeaf(const label octant)
  50 {
  51     subNodeTypes_ &= ~(0x1 << octant);
  52 }
  53
  54
  55 // Set pointer to sub node
  56 template <class Type>
  57 void Foam::treeNode<Type>::setNodePtr
  58 (
  59     const label octant,
  60     treeElem<Type>* treeNodePtr
  61 )
  62 {
  63     setAsNode(octant);
  64     subNodes_[octant] = treeNodePtr;
  65 }
  66
  67
  68 // Set pointer to sub leaf
  69 template <class Type>
  70 void Foam::treeNode<Type>::setLeafPtr
  71 (
  72     const label octant,
  73     treeElem<Type>* treeLeafPtr
  74 )
  75 {
  76     setAsLeaf(octant);
  77     subNodes_[octant] = treeLeafPtr;
  78 }
  79
  80
  81 template <class Type>
  82 void Foam::treeNode<Type>::setVolType
  83 (
  84     const label octant,
  85     const label type
  86 )
  87 {
  88     if ((type < 0) || (type > 3))
  89     {
  90         FatalErrorIn("treeNode<Type>::setVolType(const label, const label)")
  91             << "Type " << type << " not within range 0..3" << endl;
  92     }
  93
  94     // Clear out two bits at position 2*octant
  95     volType_ &= ~(0x3 << 2*octant);
  96
  97     // Add the two bits of type
  98     volType_ |= (type << 2*octant);
  99 }
 100
 101
 102 template <class Type>
 103 void Foam::treeNode<Type>::space(Ostream& os, const label n)
 104 {
 105     for (label i=0; i<n; i++)
 106     {
 107         os  << ' ';
 108     }
 109 }
 110
 111
 112 // look in single octant starting from <start>
 113 template <class Type>
 114 const Foam::treeLeaf<Type>* Foam::treeNode<Type>::findLeafLineOctant
 115 (
 116     const int level,
 117     const Type& shapes,
 118     const label octant,
 119     const vector& direction,
 120     point& start,
 121     const point& end
 122 ) const
 123 {
 124     static const char* functionName =
 125         "treeNode<Type>::findLeafLineOctant"
 126         "(const int, const Type&, const label, const vector&,"
 127         " point&, const point&)";
 128
 129     if (debug & 2)
 130     {
 131         space(Pout, 2*level);
 132         Pout<< "findLeafLineOctant : bb:" << this->bb()
 133             << "  start:" << start
 134             << "  end:" << end
 135             << "  mid:" << midpoint()
 136             << " Searching octant:" << octant
 137             << endl;
 138     }
 139
 140     if (subNodes()[octant])
 141     {
 142         if (isNode(octant))
 143         {
 144             // Node: recurse into subnodes
 145             const treeNode<Type>* subNodePtr = getNodePtr(octant);
 146
 147             if (subNodePtr->bb().contains(direction, start))
 148             {
 149                 // Search on lower level
 150                 const treeLeaf<Type>* subLeafPtr = subNodePtr->findLeafLine
 151                 (
 152                     level + 1,
 153                     shapes,
 154                     start,
 155                     end
 156                 );
 157
 158                 if (debug & 2)
 159                 {
 160                     space(Pout, 2*level);
 161                     Pout<< "findLeafLineOctant : bb:" << this->bb()
 162                         << " returning from sub treeNode"
 163                         << " with start:" << start << "  subLeaf:"
 164                         << long(subLeafPtr) << endl;
 165                 }
 166
 167                 return subLeafPtr;
 168             }
 169             else
 170             {
 171                 FatalErrorIn(functionName)
 172                     << "Sub node " << subNodePtr->bb()
 173                     << " at octant " << octant
 174                     << " does not contain start " << start
 175                     << abort(FatalError);
 176             }
 177         }
 178         else
 179         {
 180             // Leaf
 181             const treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
 182
 183             if (subLeafPtr->bb().contains(direction, start))
 184             {
 185                 // Step to end of subleaf bb
 186                 point tmp;
 187                 if
 188                 (
 189                     !subLeafPtr->bb().intersects
 190                     (
 191                         end,
 192                         start,
 193                         tmp
 194                     )
 195                 )
 196                 {
 197                     FatalErrorIn(functionName)
 198                         << "Sub leaf contains start " << start
 199                         << " but line does not intersect its bb "
 200                         << subLeafPtr->bb()
 201                         << abort(FatalError);
 202                 }
 203                 start = tmp;
 204
 205                 if (debug & 2)
 206                 {
 207                     space(Pout, 2*level);
 208                     Pout<< "findLeafLineOctant : returning from intersecting"
 209                         << " treeLeaf " << subLeafPtr->bb()
 210                         << " with start:" << start << "  subLeaf:"
 211                         << long(subLeafPtr) << endl;
 212                 }
 213
 214                 return subLeafPtr;
 215             }
 216             else
 217             {
 218                 FatalErrorIn(functionName)
 219                     << "Sub leaf " << subLeafPtr->bb()
 220                     << " at octant " << octant
 221                     << " does not contain start " << start
 222                     << abort(FatalError);
 223             }
 224         }
 225     }
 226     else
 227     {
 228         // Empty subNode. Transfer across.
 229         const treeBoundBox emptyBb = this->bb().subBbox(midpoint(), octant);
 230
 231         if (emptyBb.contains(direction, start))
 232         {
 233             if (debug & 2)
 234             {
 235                 space(Pout, 2*level);
 236                 Pout<< "findLeafLineOctant : Empty node. Octant:" << octant
 237                     << "  start:" << start
 238                     << "  bb:" << this->bb()
 239                     << "  emptyBb:" << emptyBb << endl;
 240             }
 241
 242             // Update start by clipping to emptyBb
 243             point tmp;
 244             if
 245             (
 246                 !emptyBb.intersects
 247                 (
 248                     end,
 249                     start,
 250                     tmp
 251                 )
 252             )
 253             {
 254                 FatalErrorIn(functionName)
 255                     << "Empty node contains start " << start
 256                     << " but line does not intersect its (calculated)"
 257                     << " bb " << emptyBb
 258                     << endl << "This might be due to truncation error"
 259                     << abort(FatalError);
 260             }
 261             start = tmp;
 262
 263             if (debug & 2)
 264             {
 265                 space(Pout, 2*level);
 266                 Pout<< "findLeafLineOctant : returning from intersecting with"
 267                     << " empty " << emptyBb
 268                     << " with start:" << start << "  subLeaf:" << 0 << endl;
 269             }
 270
 271             return NULL;
 272         }
 273         else
 274         {
 275             FatalErrorIn(functionName)
 276                 << "Empty node " << emptyBb
 277                 << " at octant " << octant
 278                 << " does not contain start " << start
 279                 << abort(FatalError);
 280         }
 281     }
 282
 283     FatalErrorIn(functionName)
 284         << "Octant " << octant << " of cube " << this->bb()
 285         << " does not contain start " << start
 286         << abort(FatalError);
 287
 288     return NULL;
 289 }
 290
 291
 292 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 293
 294 // Construct from components
 295 template <class Type>
 296 Foam::treeNode<Type>::treeNode(const treeBoundBox& bb)
 297 :
 298     treeElem<Type>(bb),
 299     treeNodeName(),
 300     mid_(bb.midpoint()),
 301     subNodeTypes_(0),
 302     volType_(0)
 303 {
 304     for (label octantI=0; octantI<8; octantI++)
 305     {
 306         subNodes_[octantI] = NULL;
 307         setVolType(octantI, octree<Type>::UNKNOWN);
 308     }
 309 }
 310
 311
 312 // Construct from Istream
 313 template <class Type>
 314 Foam::treeNode<Type>::treeNode(Istream& is)
 315 {
 316     is >> *this;
 317 }
 318
 319
 320 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 321
 322 template <class Type>
 323 Foam::treeNode<Type>::~treeNode()
 324 {
 325     for (int octant=0; octant<8; octant++)
 326     {
 327         if (subNodes()[octant])
 328         {
 329             if (isNode(octant))
 330             {
 331                 delete getNodePtr(octant);
 332             }
 333             else
 334             {
 335                 delete getLeafPtr(octant);
 336             }
 337         }
 338     }
 339 }
 340
 341
 342 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 343
 344 // Distributes cells to subLeaves
 345 template <class Type>
 346 void Foam::treeNode<Type>::distribute
 347 (
 348     const label level,
 349     octree<Type>& top,
 350     const Type& shapes,
 351     const labelList& indices
 352 )
 353 {
 354     if (debug & 1)
 355     {
 356         space(Pout, level);
 357         Pout<< "treeNode::distributing " << indices.size() << endl;
 358     }
 359
 360     // Create subLeaves if necessary
 361     for (label octant=0; octant<8; octant++)
 362     {
 363         if (subNodes()[octant])
 364         {
 365             printNode(Pout, level);
 366             FatalErrorIn
 367             (
 368                 "treeNode<Type>::distribute(const label, octree<Type>&, "
 369                 "const Type&, const labelList&)"
 370             )   << "subNode already available at octant:" << octant
 371                 << abort(FatalError);
 372         }
 373         else
 374         {
 375             treeLeaf<Type>* subLeafPtr = new treeLeaf<Type>
 376             (
 377                 this->bb().subBbox(midpoint(), octant),
 378                 indices.size()
 379             );
 380
 381             top.setLeaves(top.nLeaves() + 1);
 382             setLeafPtr(octant, subLeafPtr);
 383         }
 384     }
 385
 386
 387     // add cells to correct sub leaf
 388     forAll(indices, i)
 389     {
 390         const label shapei = indices[i];
 391
 392         for (label octant=0; octant<8; octant++)
 393         {
 394             treeLeaf<Type>* leafPtr = getLeafPtr(octant);
 395
 396             if (shapes.overlaps(shapei, leafPtr->bb()))
 397             {
 398                 if (debug == 1)
 399                 {
 400                     space(Pout, level);
 401                     Pout<< "inserting " << shapei;
 402                     shapes.write(Pout, shapei);
 403                     Pout<< " into " << leafPtr->bb() << endl;
 404                 }
 405                 leafPtr->insert(shapei);
 406                 top.setEntries(top.nEntries() + 1);
 407             }
 408         }
 409     }
 410
 411     // Trim size of subLeaves
 412     for (label octant=0; octant<8; octant++)
 413     {
 414         treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
 415
 416         if (subLeafPtr->size() == 0)
 417         {
 418             // Contains no data. Delete.
 419             setLeafPtr(octant, NULL);
 420             delete subLeafPtr;
 421             top.setLeaves(top.nLeaves() - 1);
 422         }
 423         else
 424         {
 425             // Trim to actual size.
 426             subLeafPtr->trim();
 427         }
 428     }
 429
 430     if (debug & 1)
 431     {
 432         space(Pout, level);
 433         Pout<< "end of treeNode::distribute" << endl;
 434     }
 435 }
 436
 437
 438 // Descends to refineLevel and checks the subLeaves for redistribution
 439 template <class Type>
 440 void Foam::treeNode<Type>::redistribute
 441 (
 442     const label level,
 443     octree<Type>& top,
 444     const Type& shapes,
 445     const label refineLevel
 446 )
 447 {
 448     if (debug & 1)
 449     {
 450         space(Pout, level);
 451         Pout<< "treeNode::redistribute with level:" << level
 452             << "  refineLevel:" << refineLevel << endl;
 453     }
 454
 455     // Descend to correct level
 456     if (level < refineLevel)
 457     {
 458         for (label octant=0; octant<8; octant++)
 459         {
 460             if (subNodes()[octant])
 461             {
 462                 if (isNode(octant))
 463                 {
 464                     getNodePtr(octant)->redistribute
 465                     (
 466                         level+1,
 467                         top,
 468                         shapes,
 469                         refineLevel
 470                     );
 471                 }
 472             }
 473         }
 474     }
 475     else
 476     {
 477         // Reached correct (should also be deepest) level of treeNode
 478         if (debug & 1)
 479         {
 480             space(Pout, level);
 481             Pout<< "treeNode::redistribute : now at correct level" << endl;
 482         }
 483
 484         // handle redistribution of sub leaves
 485         for (label octant=0; octant<8; octant++)
 486         {
 487             if (subNodes()[octant])
 488             {
 489                 if (isNode(octant))
 490                 {
 491                     FatalErrorIn
 492                     (
 493                         "treeNode<Type>::redistribute(const int, octree& top,"
 494                         "const int, const treeBoundBox&)"
 495                     )   << "found treeNode instead of treeLeaf" << endl
 496                         << abort(FatalError);
 497                 }
 498                 else
 499                 {
 500                     treeLeaf<Type>* leafPtr = getLeafPtr(octant);
 501
 502                     treeLeaf<Type>* newSubPtr = leafPtr->redistribute
 503                     (
 504                         level,
 505                         top,
 506                         shapes
 507                     );
 508
 509                     if (newSubPtr && (newSubPtr != leafPtr))
 510                     {
 511                         // redistribute has created nodePtr
 512                         // so delete no longer used subPtr and update info
 513                         if (debug & 1)
 514                         {
 515                             Pout<< "deleting "
 516                                 << top.nEntries() - leafPtr->size()
 517                                 << " entries" << endl;
 518                         }
 519                         top.setEntries(top.nEntries() - leafPtr->size());
 520
 521                         delete leafPtr;
 522
 523                         top.setLeaves(top.nLeaves() - 1);
 524
 525                         setNodePtr(octant, newSubPtr);
 526                     }
 527                 }
 528             }
 529         }
 530         if (debug & 1)
 531         {
 532             space(Pout, level);
 533             Pout<< "end of treeNode::redistribute for correct level" << endl;
 534         }
 535     }
 536
 537     if (debug & 1)
 538     {
 539         space(Pout, level);
 540         Pout<< "return from treeNode::redistribute with bb:" << this->bb()
 541             << endl;
 542     }
 543 }
 544
 545
 546 // Set type of node.
 547 template <class Type>
 548 Foam::label Foam::treeNode<Type>::setSubNodeType
 549 (
 550     const label level,
 551     octree<Type>& top,
 552     const Type& shapes
 553 )
 554 {
 555     if (debug & 4)
 556     {
 557         space(Pout, level);
 558         Pout<< "treeNode::setSubNodeType with level:" << level
 559             << "   bb:" << this->bb() << endl;
 560     }
 561
 562     label myType = -1;
 563
 564     for (label octant=0; octant<8; octant++)
 565     {
 566         label subType = -1;
 567
 568         if (subNodes()[octant])
 569         {
 570             if (isNode(octant))
 571             {
 572                 subType = getNodePtr(octant)->setSubNodeType
 573                 (
 574                     level+1,
 575                     top,
 576                     shapes
 577                 );
 578             }
 579             else
 580             {
 581                 subType = getLeafPtr(octant)->setSubNodeType
 582                 (
 583                     level+1,
 584                     top,
 585                     shapes
 586                 );
 587             }
 588         }
 589         else
 590         {
 591             // No data in this one. Set type for octant acc. to its bounding
 592             // box.
 593             const treeBoundBox subBb = this->bb().subBbox(midpoint(), octant);
 594
 595             subType = shapes.getSampleType(top, subBb.midpoint());
 596         }
 597
 598         if (debug & 4)
 599         {
 600             space(Pout, level);
 601             Pout<< "treeNode::setSubNodeType : setting octant with bb:"
 602                 << this->bb().subBbox(midpoint(), octant)
 603                 << "  to type:" << octree<Type>::volType(subType) << endl;
 604         }
 605         setVolType(octant, subType);
 606
 607         // Combine sub node types into type for treeNode. Result is 'mixed' if
 608         // types differ among subnodes.
 609         if (myType == -1)
 610         {
 611             myType = subType;
 612         }
 613         else if (subType != myType)
 614         {
 615             myType = octree<Type>::MIXED;
 616         }
 617     }
 618
 619     if (debug & 4)
 620     {
 621         space(Pout, level);
 622         Pout<< "return from treeNode::setSubNodeType with type:"
 623             << octree<Type>::volType(myType)
 624             << "  bb:" << this->bb() << endl;
 625     }
 626
 627     return myType;
 628 }
 629
 630
 631 // Get type of node.
 632 template <class Type>
 633 Foam::label Foam::treeNode<Type>::getSampleType
 634 (
 635     const label level,
 636     const octree<Type>& top,
 637     const Type& shapes,
 638     const point& sample
 639 ) const
 640 {
 641     if (debug & 4)
 642     {
 643         space(Pout, level);
 644         Pout<< "treeNode::getSampleType with level:" << level
 645             << " bb:" << this->bb() << "  sample:" << sample << endl;
 646     }
 647
 648     // Determine octant of bb. If on edge just use whichever octant.
 649     bool onEdge = false;
 650
 651     label octant = this->bb().subOctant(midpoint(), sample, onEdge);
 652
 653     label type = getVolType(octant);
 654
 655     if (type == octree<Type>::MIXED)
 656     {
 657         // At this level multiple sub types. Recurse to resolve.
 658         if (subNodes()[octant])
 659         {
 660             if (isNode(octant))
 661             {
 662                 // Node: recurse into subnodes
 663                 type = getNodePtr(octant)->getSampleType
 664                 (
 665                     level + 1,
 666                     top,
 667                     shapes,
 668                     sample
 669                 );
 670             }
 671             else
 672             {
 673                 // Leaf
 674                 type = getLeafPtr(octant)->getSampleType
 675                 (
 676                     level + 1,
 677                     top,
 678                     shapes,
 679                     sample
 680                 );
 681             }
 682         }
 683         else
 684         {
 685             // Problem: empty subnode should have a type
 686             FatalErrorIn
 687             (
 688                 "treeNode<Type>::getSampleType"
 689                 "(const label, octree<Type>&, const Type&, const point&)"
 690             )   << "Empty node bb:" << this->bb().subBbox(midpoint(), octant)
 691                 << " has non-mixed type:"
 692                 << octree<Type>::volType(type)
 693                 << abort(FatalError);
 694         }
 695     }
 696
 697     if (type == octree<Type>::MIXED)
 698     {
 699         FatalErrorIn
 700         (
 701             "treeNode<Type>::getSampleType"
 702             "(const label, octree<Type>&, const Type&, const point&)"
 703         )   << "Type is MIXED when searching for " << sample
 704             << " at level " << this->bb() << endl
 705             << "This probably is because the octree has not been constructed"
 706             << " with search facility." << exit(FatalError);
 707     }
 708
 709     if (debug & 4)
 710     {
 711         space(Pout, level);
 712         Pout<< "return from treeNode::getSampleType with type:"
 713             << octree<Type>::volType(type)
 714             << "  bb:" << this->bb()
 715             << "  sample:" << sample << endl;
 716     }
 717     return type;
 718 }
 719
 720
 721 template <class Type>
 722 Foam::label Foam::treeNode<Type>::find
 723 (
 724     const Type& shapes,
 725     const point& sample
 726 ) const
 727 {
 728     // Find octant of sample. Don't care if on edge (since any item on edge
 729     // will have been inserted in both subcubes)
 730     bool onEdge = false;
 731
 732     label octant = this->bb().subOctant(midpoint(), sample, onEdge);
 733
 734     if (subNodes()[octant])
 735     {
 736         if (isNode(octant))
 737         {
 738             // Node: recurse into subnodes
 739             return getNodePtr(octant)->find(shapes, sample);
 740         }
 741         else
 742         {
 743             // Leaf: let leaf::find handle this
 744             return getLeafPtr(octant)->find(shapes, sample);
 745         }
 746     }
 747     return -1;
 748 }
 749
 750
 751 template <class Type>
 752 bool Foam::treeNode<Type>::findTightest
 753 (
 754     const Type& shapes,
 755     const point& sample,
 756     treeBoundBox& tightest
 757 ) const
 758 {
 759     bool changed = false;
 760     bool onEdge = false;
 761     // Estimate for best place to start searching
 762     label sampleOctant = this->bb().subOctant(midpoint(), sample, onEdge);
 763
 764     // Go into all suboctants (one containing sample first) and update tightest.
 765     // Order of visiting is if e.g. sampleOctant = 5:
 766     //  5 1 2 3 4 0 6 7
 767     for (label octantI=0; octantI<8; octantI++)
 768     {
 769         label octant;
 770         if (octantI == 0)
 771         {
 772             // Use sampleOctant first
 773             octant = sampleOctant;
 774         }
 775         else if (octantI == sampleOctant)
 776         {
 777             octant = 0;
 778         }
 779         else
 780         {
 781             octant = octantI;
 782         }
 783
 784         if (subNodes()[octant])
 785         {
 786             if (isNode(octant))
 787             {
 788                 // Node: recurse into subnodes
 789                 const treeNode<Type>* subNodePtr = getNodePtr(octant);
 790
 791                 if (subNodePtr->bb().overlaps(tightest))
 792                 {
 793                     // there might be a better fit inside this subNode
 794                     changed |= subNodePtr->findTightest
 795                     (
 796                         shapes,
 797                         sample,
 798                         tightest
 799                     );
 800                 }
 801             }
 802             else
 803             {
 804                 // Leaf: let leaf::find handle this
 805                 const treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
 806
 807                 if (subLeafPtr->bb().overlaps(tightest))
 808                 {
 809                     // there might be a better fit inside this subLeaf
 810                     changed |= subLeafPtr->findTightest
 811                     (
 812                         shapes,
 813                         sample,
 814                         tightest
 815                     );
 816                 }
 817             }
 818         }
 819     }
 820
 821     return changed;
 822 }
 823
 824
 825 template <class Type>
 826 bool Foam::treeNode<Type>::findNearest
 827 (
 828     const Type& shapes,
 829     const point& sample,
 830     treeBoundBox& tightest,
 831     label& tightestI,
 832     scalar& tightestDist
 833 ) const
 834 {
 835     if (debug & 8)
 836     {
 837         Pout<< "In findNearest with sample:" << sample << " cube:"
 838             << this->bb() << " tightest:" << tightest << endl;
 839     }
 840
 841     bool changed = false;
 842     bool onEdge = false;
 843     // Estimate for best place to start searching
 844     label sampleOctant = this->bb().subOctant(midpoint(), sample, onEdge);
 845
 846     // Go into all suboctants (one containing sample first) and update tightest.
 847     // Order of visiting is if e.g. sampleOctant = 5:
 848     //  5 1 2 3 4 0 6 7
 849     for (label octantI=0; octantI<8; octantI++)
 850     {
 851         label octant;
 852         if (octantI == 0)
 853         {
 854             // Use sampleOctant first
 855             octant = sampleOctant;
 856         }
 857         else if (octantI == sampleOctant)
 858         {
 859             octant = 0;
 860         }
 861         else
 862         {
 863             octant = octantI;
 864         }
 865
 866         if (subNodes()[octant])
 867         {
 868             if (isNode(octant))
 869             {
 870                 // Node
 871                 const treeNode<Type>* subNodePtr = getNodePtr(octant);
 872
 873                 if (subNodePtr->bb().overlaps(tightest))
 874                 {
 875                     // there might be a better fit inside this subNode
 876                     changed |= subNodePtr->findNearest
 877                     (
 878                         shapes,
 879                         sample,
 880                         tightest,
 881                         tightestI,
 882                         tightestDist
 883                     );
 884                 }
 885             }
 886             else
 887             {
 888                 // Leaf: let leaf::find handle this
 889                 const treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
 890
 891                 if (subLeafPtr->bb().overlaps(tightest))
 892                 {
 893                     // there might be a better fit inside this subNode
 894                     changed |= subLeafPtr->findNearest
 895                     (
 896                         shapes,
 897                         sample,
 898                         tightest,
 899                         tightestI,
 900                         tightestDist
 901                     );
 902                 }
 903             }
 904         }
 905     }
 906
 907     if (debug & 8)
 908     {
 909         Pout<< "Exiting findNearest for sample:" << sample << " cube:"
 910             << this->bb() << " tightestI:" << tightestI << endl;
 911     }
 912
 913     return changed;
 914 }
 915
 916
 917 template <class Type>
 918 bool Foam::treeNode<Type>::findNearest
 919 (
 920     const Type& shapes,
 921     const linePointRef& ln,
 922     treeBoundBox& tightest,
 923     label& tightestI,
 924     point& linePoint,   // nearest point on line
 925     point& shapePoint   // nearest point on shape
 926 ) const
 927 {
 928     bool changed = false;
 929     bool onEdge = false;
 930     // Estimate for best place to start searching
 931     label sampleOctant = this->bb().subOctant(midpoint(), ln.centre(), onEdge);
 932
 933     // Go into all suboctants (one containing sample first) and update tightest.
 934     // Order of visiting is if e.g. sampleOctant = 5:
 935     //  5 1 2 3 4 0 6 7
 936     for (label octantI=0; octantI<8; octantI++)
 937     {
 938         label octant;
 939         if (octantI == 0)
 940         {
 941             // Use sampleOctant first
 942             octant = sampleOctant;
 943         }
 944         else if (octantI == sampleOctant)
 945         {
 946             octant = 0;
 947         }
 948         else
 949         {
 950             octant = octantI;
 951         }
 952
 953         if (subNodes()[octant])
 954         {
 955             if (isNode(octant))
 956             {
 957                 // Node
 958                 const treeNode<Type>* subNodePtr = getNodePtr(octant);
 959
 960                 if (subNodePtr->bb().overlaps(tightest))
 961                 {
 962                     // there might be a better fit inside this subNode
 963                     changed |= subNodePtr->findNearest
 964                     (
 965                         shapes,
 966                         ln,
 967                         tightest,
 968                         tightestI,
 969                         linePoint,
 970                         shapePoint
 971                     );
 972                 }
 973             }
 974             else
 975             {
 976                 // Leaf: let leaf::find handle this
 977                 const treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
 978
 979                 if (subLeafPtr->bb().overlaps(tightest))
 980                 {
 981                     // there might be a better fit inside this subNode
 982                     changed |= subLeafPtr->findNearest
 983                     (
 984                         shapes,
 985                         ln,
 986                         tightest,
 987                         tightestI,
 988                         linePoint,
 989                         shapePoint
 990                     );
 991                 }
 992             }
 993         }
 994     }
 995
 996     return changed;
 997 }
 998
 999
1000 template <class Type>
1001 bool Foam::treeNode<Type>::findBox
1002 (
1003     const Type& shapes,
1004     const boundBox& box,
1005     labelHashSet& elements
1006 ) const
1007 {
1008     bool changed = false;
1009     bool onEdge = false;
1010     // Estimate for best place to start searching
1011     label sampleOctant = this->bb().subOctant
1012     (
1013         midpoint(),
1014         box.midpoint(),
1015         onEdge
1016     );
1017
1018     // Go into all suboctants (one containing sample first) and update tightest.
1019     // Order of visiting is if e.g. sampleOctant = 5:
1020     //  5 1 2 3 4 0 6 7
1021     for (label octantI=0; octantI<8; octantI++)
1022     {
1023         label octant;
1024         if (octantI == 0)
1025         {
1026             // Use sampleOctant first
1027             octant = sampleOctant;
1028         }
1029         else if (octantI == sampleOctant)
1030         {
1031             octant = 0;
1032         }
1033         else
1034         {
1035             octant = octantI;
1036         }
1037
1038         if (subNodes()[octant])
1039         {
1040             if (isNode(octant))
1041             {
1042                 // Node
1043                 const treeNode<Type>* subNodePtr = getNodePtr(octant);
1044
1045                 if (subNodePtr->bb().overlaps(box))
1046                 {
1047                     // Visit sub node.
1048                     changed |= subNodePtr->findBox(shapes, box, elements);
1049                 }
1050             }
1051             else
1052             {
1053                 // Leaf: let leaf::find handle this
1054                 const treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
1055
1056                 if (subLeafPtr->bb().overlaps(box))
1057                 {
1058                     // Visit sub leaf.
1059                     changed |= subLeafPtr->findBox(shapes, box, elements);
1060                 }
1061             }
1062         }
1063     }
1064
1065     return changed;
1066 }
1067
1068
1069 // look from <start> in current cube (given by this->bb()).
1070 template <class Type>
1071 const Foam::treeLeaf<Type>* Foam::treeNode<Type>::findLeafLine
1072 (
1073     const int level,
1074     const Type& shapes,
1075     point& start,
1076     const point& end
1077 ) const
1078 {
1079     if (debug & 2)
1080     {
1081         space(Pout, 2*level);
1082         Pout<< "findLeafLine : bb:" << this->bb() << "  mid:" << midpoint()
1083             << "  start:" << start << endl;
1084     }
1085
1086     scalar typDim = this->bb().avgDim();
1087
1088     const vector direction = end - start;
1089
1090     // Loop on current level until start has been updated to be outside
1091     // of this->bb(). Note that max only four subcubes can be crossed so this is
1092     // check on whether there are any truncation error problems.
1093
1094     label iter = 0;
1095
1096     while (true)
1097     {
1098         if (!this->bb().contains(direction, start))
1099         {
1100             if (debug & 2)
1101             {
1102                 space(Pout, 2*level);
1103                 Pout<< "findLeafLine : Start not inside bb " << this->bb()
1104                     << ". Returning with start:" << start << "  subLeaf:"
1105                     << 0 << endl;
1106             }
1107             return NULL;
1108         }
1109
1110         // Check if start and <end> equal
1111         if ((mag(start - end)/typDim) < SMALL)
1112         {
1113             if (debug & 2)
1114             {
1115                 space(Pout, 2*level);
1116                 Pout<< "findLeafLine : start equals end"
1117                     << ". Returning with start:" << start << "  subLeaf:"
1118                     << 0 << endl;
1119             }
1120             return NULL;
1121         }
1122
1123         if (iter >= 4)
1124         {
1125             // Too many iterations. Is hanging. Handle outside of loop.
1126             break;
1127         }
1128
1129         bool onEdge = false;
1130         label octant = this->bb().subOctant
1131         (
1132             midpoint(), direction, start, onEdge
1133         );
1134
1135         // Try finding non-empty treeleaf in octant
1136         const treeLeaf<Type>* leafPtr = findLeafLineOctant
1137         (
1138             level,
1139             shapes,
1140             octant,
1141             direction,
1142             start,
1143             end
1144         );
1145
1146         if (leafPtr)
1147         {
1148             // Found treeLeaf -> return
1149             if (debug & 2)
1150             {
1151                 space(Pout, 2*level);
1152                 Pout<< "findLeafLine : Found treeLeaf"
1153                         << ". Returning with start:" << start << "  subLeaf:"
1154                         << long(leafPtr) << endl;
1155             }
1156
1157             return leafPtr;
1158         }
1159
1160         iter++;
1161     }
1162
1163     // Check if is hanging. Max 4 octants can be crossed by a straight line
1164     FatalErrorIn
1165     (
1166         "treeNode<Type>::findLeafLine"
1167         "(const label, octree<Type>&, point&,"
1168         " const point&)"
1169     )   << "Did not leave bb " << this->bb()
1170         << " after " << iter
1171         << " iterations of updating starting point."
1172         << "start:" << start << "  end:" << end
1173         << abort(FatalError);
1174
1175     return NULL;
1176 }
1177
1178
1179 template <class Type>
1180 void Foam::treeNode<Type>::findLeaves
1181 (
1182     List<treeLeaf<Type>*>& leafArray,
1183     label& leafIndex
1184 ) const
1185 {
1186     // Go into all sub boxes
1187     for (label octant=0; octant<8; octant++)
1188     {
1189         if (subNodes()[octant])
1190         {
1191             if (isNode(octant))
1192             {
1193                 // Node: recurse into subnodes
1194                 const treeNode<Type>* subNodePtr = getNodePtr(octant);
1195                 subNodePtr->findLeaves(leafArray, leafIndex);
1196             }
1197             else
1198             {
1199                 // Leaf: store
1200                 treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
1201                 leafArray[leafIndex++] = subLeafPtr;
1202             }
1203         }
1204     }
1205 }
1206
1207
1208 template <class Type>
1209 void Foam::treeNode<Type>::findLeaves
1210 (
1211     List<const treeLeaf<Type>*>& leafArray,
1212     label& leafIndex
1213 ) const
1214 {
1215     // Go into all sub boxes
1216     for (label octant=0; octant<8; octant++)
1217     {
1218         if (subNodes()[octant])
1219         {
1220             if (isNode(octant))
1221             {
1222                 // Node: recurse into subnodes
1223                 const treeNode<Type>* subNodePtr = getNodePtr(octant);
1224                 subNodePtr->findLeaves(leafArray, leafIndex);
1225             }
1226             else
1227             {
1228                 // Leaf: store
1229                 treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
1230                 leafArray[leafIndex++] = subLeafPtr;
1231             }
1232         }
1233     }
1234 }
1235
1236
1237 template <class Type>
1238 void Foam::treeNode<Type>::printNode
1239 (
1240     Ostream& os,
1241     const label level
1242 ) const
1243 {
1244     space(os, 2*level);
1245
1246     os << "node:" << this->bb() << endl;
1247
1248     for (label octant=0; octant<8; octant++)
1249     {
1250         label type = getVolType(octant);
1251
1252         string typeString = octree<Type>::volType(type);
1253
1254         if (!subNodes_[octant])
1255         {
1256             space(os, level);
1257             os << octant << ":" << typeString << " : null" << endl;
1258         }
1259         else if (isNode(octant))
1260         {
1261             space(os, level);
1262             os << octant << ":" << typeString << " : node" << endl;
1263             getNodePtr(octant)->printNode(os, level+1);
1264         }
1265         else
1266         {
1267             space(os, level);
1268             os << octant << ":" << typeString << " : leaf" << endl;
1269
1270             treeLeaf<Type>* leafPtr = getLeafPtr(octant);
1271             leafPtr->printLeaf(os, level+1);
1272         }
1273     }
1274 }
1275
1276
1277 template <class Type>
1278 void Foam::treeNode<Type>::writeOBJ
1279 (
1280     Ostream& os,
1281     const label level,
1282     label& vertNo
1283 ) const
1284 {
1285     point midPoint(this->bb().midpoint());
1286
1287     label midVertNo = vertNo;
1288     os << "v " << midPoint.x() << " " << midPoint.y() << " "
1289        << midPoint.z() << endl;
1290     vertNo++;
1291
1292     for (label octant=0; octant<8; octant++)
1293     {
1294         if (subNodes_[octant])
1295         {
1296             if (isNode(octant))
1297             {
1298                 treeNode<Type>* nodePtr = getNodePtr(octant);
1299
1300                 point subMidPoint(nodePtr->bb().midpoint());
1301                 os << "v " << subMidPoint.x() << " " << subMidPoint.y() << " "
1302                    << subMidPoint.z() << endl;
1303                 os << "l " << midVertNo + 1<< " " << vertNo + 1 << endl;
1304                 vertNo++;
1305
1306                 nodePtr->writeOBJ(os, level+1, vertNo);
1307             }
1308             else
1309             {
1310                 treeLeaf<Type>* leafPtr = getLeafPtr(octant);
1311
1312                 point subMidPoint(leafPtr->bb().midpoint());
1313                 os << "v " << subMidPoint.x() << " " << subMidPoint.y() << " "
1314                    << subMidPoint.z() << endl;
1315                 os << "l " << midVertNo + 1<< " " << vertNo + 1 << endl;
1316                 vertNo++;
1317
1318                 //leafPtr->writeOBJ(os, level+1, vertNo);
1319             }
1320         }
1321     }
1322 }
1323
1324
1325 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
1326
1327 template <class Type>
1328 Foam::Istream& Foam::operator>>(Istream& is, treeNode<Type>& oc)
1329 {
1330     for (label octant = 0; octant < 8; octant++)
1331     {
1332         oc.subNodes_[octant] = NULL;
1333     }
1334
1335     is >> oc.bb();
1336
1337     label nPtrs;
1338
1339     // Read number of entries folllowing
1340     is >> nPtrs;
1341
1342     is.readBegin("treeNode");
1343     for (label octant = 0; octant < nPtrs; octant++)
1344     {
1345         label index;
1346         is >> index;
1347
1348         if (index >= treeNode<Type>::leafOffset)
1349         {
1350             // Leaf recognized by out of range index
1351             treeLeaf<Type>* leafPtr = new treeLeaf<Type>(is);
1352             oc.setLeafPtr(index - treeNode<Type>::leafOffset, leafPtr);
1353         }
1354         else
1355         {
1356             oc.setNodePtr(index, new treeNode<Type>(is));
1357         }
1358     }
1359
1360     // Read end of treeNode list
1361     is.readEnd("treeNode");
1362
1363     // Check state of Istream
1364     is.check("Istream& operator>>(Istream&, treeNode&)");
1365
1366     return is;
1367 }
1368
1369
1370 template <class Type>
1371 Foam::Ostream& Foam::operator<<(Ostream& os, const treeNode<Type>& tn)
1372 {
1373     // Count valid subnodes:
1374     //   - treeNode
1375     //   - treeLeafs with non-zero cell list.
1376     label nPtrs = 0;
1377     for (label octant = 0; octant < 8; octant++)
1378     {
1379         if (tn.subNodes_[octant])
1380         {
1381             if (tn.isNode(octant) || tn.getLeafPtr(octant)->indices().size())
1382             {
1383                 nPtrs++;
1384             }
1385         }
1386     }
1387
1388
1389     // output subnodes as list of length nPtrs
1390     os << token::SPACE << tn.bb() << token::SPACE << nPtrs
1391        << token::SPACE << token::BEGIN_LIST;
1392
1393     for (label octant = 0; octant < 8; octant++)
1394     {
1395         if (tn.subNodes_[octant])
1396         {
1397             if (tn.isNode(octant))
1398             {
1399                 const treeNode<Type>* subNodePtr = tn.getNodePtr(octant);
1400
1401                 // Node: output index, value
1402                 os  << token::SPACE << octant << token::SPACE << *subNodePtr
1403                     << token::NL;
1404             }
1405             else if (tn.getLeafPtr(octant)->indices().size())
1406             {
1407                 // treeLeaf: mark by putting index invalid
1408                 const treeLeaf<Type>* subLeafPtr = tn.getLeafPtr(octant);
1409
1410                 os  << token::SPACE << octant + treeNode<Type>::leafOffset
1411                     << token::SPACE << *subLeafPtr
1412                     << token::NL;
1413             }
1414         }
1415     }
1416
1417     os  << token::SPACE << token::END_LIST;
1418
1419     return os;
1420 }
1421
1422
1423 // ************************************************************************* //