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