src/meshTools/octree/treeNode.C

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