FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

treeNode.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 Description
00025 
00026 \*---------------------------------------------------------------------------*/
00027 
00028 #include "treeNode.H"
00029 #include "octree.H"
00030 #include "treeLeaf.H"
00031 #include "treeBoundBox.H"
00032 #include <OpenFOAM/long.H>
00033 #include <OpenFOAM/linePointRef.H>
00034 
00035 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00036 
00037 template <class Type>
00038 const Foam::label Foam::treeNode<Type>::leafOffset(100);
00039 
00040 
00041 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00042 
00043 template <class Type>
00044 void Foam::treeNode<Type>::setAsNode(const label octant)
00045 {
00046     subNodeTypes_ |= (0x1 << octant);
00047 }
00048 
00049 
00050 template <class Type>
00051 void Foam::treeNode<Type>::setAsLeaf(const label octant)
00052 {
00053     subNodeTypes_ &= ~(0x1 << octant);
00054 }
00055 
00056 
00057 // Set pointer to sub node
00058 template <class Type>
00059 void Foam::treeNode<Type>::setNodePtr
00060 (
00061     const label octant,
00062     treeElem<Type>* treeNodePtr
00063 )
00064 {
00065     setAsNode(octant);
00066     subNodes_[octant] = treeNodePtr;
00067 }
00068 
00069 
00070 // Set pointer to sub leaf
00071 template <class Type>
00072 void Foam::treeNode<Type>::setLeafPtr
00073 (
00074     const label octant,
00075     treeElem<Type>* treeLeafPtr
00076 )
00077 {
00078     setAsLeaf(octant);
00079     subNodes_[octant] = treeLeafPtr;
00080 }
00081 
00082 
00083 template <class Type>
00084 void Foam::treeNode<Type>::setVolType
00085 (
00086     const label octant,
00087     const label type
00088 )
00089 {
00090     if ((type < 0) || (type > 3))
00091     {
00092         FatalErrorIn("treeNode<Type>::setVolType(const label, const label)")
00093             << "Type " << type << " not within range 0..3" << endl;
00094     }
00095 
00096     // Clear out two bits at position 2*octant
00097     volType_ &= ~(0x3 << 2*octant);
00098 
00099     // Add the two bits of type
00100     volType_ |= (type << 2*octant);
00101 }
00102 
00103 
00104 template <class Type>
00105 void Foam::treeNode<Type>::space(Ostream& os, const label n)
00106 {
00107     for (label i=0; i<n; i++)
00108     {
00109         os<< ' ';
00110     }
00111 }
00112 
00113 
00114 // look in single octant starting from <start>
00115 template <class Type>
00116 const Foam::treeLeaf<Type>* Foam::treeNode<Type>::findLeafLineOctant
00117 (
00118     const int level,
00119     const Type& shapes,
00120     const label octant,
00121     const vector& direction,
00122     point& start,
00123     const point& end
00124 ) const
00125 {
00126     static const char* functionName =
00127         "treeNode<Type>::findLeafLineOctant"
00128         "(const int, const Type&, const label, const vector&,"
00129         " point&, const point&)";
00130 
00131     if (debug & 2)
00132     {
00133         space(Pout, 2*level);
00134         Pout<< "findLeafLineOctant : bb:" << this->bb()
00135             << "  start:" << start
00136             << "  end:" << end
00137             << "  mid:" << midpoint()
00138             << " Searching octant:" << octant
00139             << endl;
00140     }
00141 
00142     if (subNodes()[octant])
00143     {
00144         if (isNode(octant))
00145         {
00146             // Node: recurse into subnodes
00147             const treeNode<Type>* subNodePtr = getNodePtr(octant);
00148 
00149             if (subNodePtr->bb().contains(direction, start))
00150             {
00151                 // Search on lower level
00152                 const treeLeaf<Type>* subLeafPtr = subNodePtr->findLeafLine
00153                 (
00154                     level + 1,
00155                     shapes,
00156                     start,
00157                     end
00158                 );
00159 
00160                 if (debug & 2)
00161                 {
00162                     space(Pout, 2*level);
00163                     Pout<< "findLeafLineOctant : bb:" << this->bb()
00164                         << " returning from sub treeNode"
00165                         << " with start:" << start << "  subLeaf:"
00166                         << long(subLeafPtr) << endl;
00167                 }
00168 
00169                 return subLeafPtr;
00170             }
00171             else
00172             {
00173                 FatalErrorIn(functionName)
00174                     << "Sub node " << subNodePtr->bb()
00175                     << " at octant " << octant
00176                     << " does not contain start " << start
00177                     << abort(FatalError);
00178             }
00179         }
00180         else
00181         {
00182             // Leaf
00183             const treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
00184 
00185             if (subLeafPtr->bb().contains(direction, start))
00186             {
00187                 // Step to end of subleaf bb
00188                 point tmp;
00189                 if
00190                 (
00191                     !subLeafPtr->bb().intersects
00192                     (
00193                         end,
00194                         start,
00195                         tmp
00196                     )
00197                 )
00198                 {
00199                     FatalErrorIn(functionName)
00200                         << "Sub leaf contains start " << start
00201                         << " but line does not intersect its bb "
00202                         << subLeafPtr->bb()
00203                         << abort(FatalError);
00204                 }
00205                 start = tmp;
00206 
00207                 if (debug & 2)
00208                 {
00209                     space(Pout, 2*level);
00210                     Pout<< "findLeafLineOctant : returning from intersecting"
00211                         << " treeLeaf " << subLeafPtr->bb()
00212                         << " with start:" << start << "  subLeaf:"
00213                         << long(subLeafPtr) << endl;
00214                 }
00215 
00216                 return subLeafPtr;
00217             }
00218             else
00219             {
00220                 FatalErrorIn(functionName)
00221                     << "Sub leaf " << subLeafPtr->bb()
00222                     << " at octant " << octant
00223                     << " does not contain start " << start
00224                     << abort(FatalError);
00225             }
00226         }
00227     }
00228     else
00229     {
00230         // Empty subNode. Transfer across.
00231         const treeBoundBox emptyBb = this->bb().subBbox(midpoint(), octant);
00232 
00233         if (emptyBb.contains(direction, start))
00234         {
00235             if (debug & 2)
00236             {
00237                 space(Pout, 2*level);
00238                 Pout<< "findLeafLineOctant : Empty node. Octant:" << octant
00239                     << "  start:" << start
00240                     << "  bb:" << this->bb()
00241                     << "  emptyBb:" << emptyBb << endl;
00242             }
00243 
00244             // Update start by clipping to emptyBb
00245             point tmp;
00246             if
00247             (
00248                 !emptyBb.intersects
00249                 (
00250                     end,
00251                     start,
00252                     tmp
00253                 )
00254             )
00255             {
00256                 FatalErrorIn(functionName)
00257                     << "Empty node contains start " << start
00258                     << " but line does not intersect its (calculated)"
00259                     << " bb " << emptyBb
00260                     << endl << "This might be due to truncation error"
00261                     << abort(FatalError);
00262             }
00263             start = tmp;
00264 
00265             if (debug & 2)
00266             {
00267                 space(Pout, 2*level);
00268                 Pout<< "findLeafLineOctant : returning from intersecting with"
00269                     << " empty " << emptyBb
00270                     << " with start:" << start << "  subLeaf:" << 0 << endl;
00271             }
00272 
00273             return NULL;
00274         }
00275         else
00276         {
00277             FatalErrorIn(functionName)
00278                 << "Empty node " << emptyBb
00279                 << " at octant " << octant
00280                 << " does not contain start " << start
00281                 << abort(FatalError);
00282         }
00283     }
00284 
00285     FatalErrorIn(functionName)
00286         << "Octant " << octant << " of cube " << this->bb()
00287         << " does not contain start " << start
00288         << abort(FatalError);
00289 
00290     return NULL;
00291 }
00292 
00293 
00294 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00295 
00296 // Construct from components
00297 template <class Type>
00298 Foam::treeNode<Type>::treeNode(const treeBoundBox& bb)
00299 :
00300     treeElem<Type>(bb),
00301     treeNodeName(),
00302     mid_(bb.midpoint()),
00303     subNodeTypes_(0),
00304     volType_(0)
00305 {
00306     for (label octantI=0; octantI<8; octantI++)
00307     {
00308         subNodes_[octantI] = NULL;
00309         setVolType(octantI, octree<Type>::UNKNOWN);
00310     }
00311 }
00312 
00313 
00314 // Construct from Istream
00315 template <class Type>
00316 Foam::treeNode<Type>::treeNode(Istream& is)
00317 {
00318     is >> *this;
00319 }
00320 
00321 
00322 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
00323 
00324 template <class Type>
00325 Foam::treeNode<Type>::~treeNode()
00326 {
00327     for (int octant=0; octant<8; octant++)
00328     {
00329         if (subNodes()[octant])
00330         {
00331             if (isNode(octant))
00332             {
00333                 delete getNodePtr(octant);
00334             }
00335             else
00336             {
00337                 delete getLeafPtr(octant);
00338             }
00339         }
00340     }
00341 }
00342 
00343 
00344 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
00345 
00346 // Distributes cells to subLeaves
00347 template <class Type>
00348 void Foam::treeNode<Type>::distribute
00349 (
00350     const label level,
00351     octree<Type>& top,
00352     const Type& shapes,
00353     const labelList& indices
00354 )
00355 {
00356     if (debug & 1)
00357     {
00358         space(Pout, level);
00359         Pout<< "treeNode::distributing " << indices.size() << endl;
00360     }
00361 
00362     // Create subLeaves if necessary
00363     for (label octant=0; octant<8; octant++)
00364     {
00365         if (subNodes()[octant])
00366         {
00367             printNode(Pout, level);
00368             FatalErrorIn
00369             (
00370                 "treeNode<Type>::distribute(const label, octree<Type>&, "
00371                 "const Type&, const labelList&)"
00372             )   << "subNode already available at octant:" << octant
00373                 << abort(FatalError);
00374         }
00375         else
00376         {
00377             treeLeaf<Type>* subLeafPtr = new treeLeaf<Type>
00378             (
00379                 this->bb().subBbox(midpoint(), octant),
00380                 indices.size()
00381             );
00382 
00383             top.setLeaves(top.nLeaves() + 1);
00384             setLeafPtr(octant, subLeafPtr);
00385         }
00386     }
00387 
00388 
00389     // add cells to correct sub leaf
00390     forAll(indices, i)
00391     {
00392         const label shapei = indices[i];
00393 
00394         for (label octant=0; octant<8; octant++)
00395         {
00396             treeLeaf<Type>* leafPtr = getLeafPtr(octant);
00397 
00398             if (shapes.overlaps(shapei, leafPtr->bb()))
00399             {
00400                 if (debug == 1)
00401                 {
00402                     space(Pout, level);
00403                     Pout<< "inserting " << shapei;
00404                     shapes.write(Pout, shapei);
00405                     Pout<< " into " << leafPtr->bb() << endl;
00406                 }
00407                 leafPtr->insert(shapei);
00408                 top.setEntries(top.nEntries() + 1);
00409             }
00410         }
00411     }
00412 
00413     // Trim size of subLeaves
00414     for (label octant=0; octant<8; octant++)
00415     {
00416         treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
00417 
00418         if (subLeafPtr->size() == 0)
00419         {
00420             // Contains no data. Delete.
00421             setLeafPtr(octant, NULL);
00422             delete subLeafPtr;
00423             top.setLeaves(top.nLeaves() - 1);
00424         }
00425         else
00426         {
00427             // Trim to actual size.
00428             subLeafPtr->trim();
00429         }
00430     }
00431 
00432     if (debug & 1)
00433     {
00434         space(Pout, level);
00435         Pout<< "end of treeNode::distribute" << endl;
00436     }
00437 }
00438 
00439 
00440 // Descends to refineLevel and checks the subLeaves for redistribution
00441 template <class Type>
00442 void Foam::treeNode<Type>::redistribute
00443 (
00444     const label level,
00445     octree<Type>& top,
00446     const Type& shapes,
00447     const label refineLevel
00448 )
00449 {
00450     if (debug & 1)
00451     {
00452         space(Pout, level);
00453         Pout<< "treeNode::redistribute with level:" << level
00454             << "  refineLevel:" << refineLevel << endl;
00455     }
00456 
00457     // Descend to correct level
00458     if (level < refineLevel)
00459     {
00460         for (label octant=0; octant<8; octant++)
00461         {
00462             if (subNodes()[octant])
00463             {
00464                 if (isNode(octant))
00465                 {
00466                     getNodePtr(octant)->redistribute
00467                     (
00468                         level+1,
00469                         top,
00470                         shapes,
00471                         refineLevel
00472                     );
00473                 }
00474             }
00475         }
00476     }
00477     else
00478     {
00479         // Reached correct (should also be deepest) level of treeNode
00480         if (debug & 1)
00481         {
00482             space(Pout, level);
00483             Pout<< "treeNode::redistribute : now at correct level" << endl;
00484         }
00485 
00486         // handle redistribution of sub leaves
00487         for (label octant=0; octant<8; octant++)
00488         {
00489             if (subNodes()[octant])
00490             {
00491                 if (isNode(octant))
00492                 {
00493                     FatalErrorIn
00494                     (
00495                         "treeNode<Type>::redistribute(const int, octree& top,"
00496                         "const int, const treeBoundBox&)"
00497                     )   << "found treeNode instead of treeLeaf" << endl
00498                         << abort(FatalError);
00499                 }
00500                 else
00501                 {
00502                     treeLeaf<Type>* leafPtr = getLeafPtr(octant);
00503 
00504                     treeLeaf<Type>* newSubPtr = leafPtr->redistribute
00505                     (
00506                         level,
00507                         top,
00508                         shapes
00509                     );
00510 
00511                     if (newSubPtr && (newSubPtr != leafPtr))
00512                     {
00513                         // redistribute has created nodePtr
00514                         // so delete no longer used subPtr and update info
00515                         if (debug & 1)
00516                         {
00517                             Pout<< "deleting "
00518                                 << top.nEntries() - leafPtr->size()
00519                                 << " entries" << endl;
00520                         }
00521                         top.setEntries(top.nEntries() - leafPtr->size());
00522 
00523                         delete leafPtr;
00524 
00525                         top.setLeaves(top.nLeaves() - 1);
00526 
00527                         setNodePtr(octant, newSubPtr);
00528                     }
00529                 }
00530             }
00531         }
00532         if (debug & 1)
00533         {
00534             space(Pout, level);
00535             Pout<< "end of treeNode::redistribute for correct level" << endl;
00536         }
00537     }
00538 
00539     if (debug & 1)
00540     {
00541         space(Pout, level);
00542         Pout<< "return from treeNode::redistribute with bb:" << this->bb()
00543             << endl;
00544     }
00545 }
00546 
00547 
00548 // Set type of node.
00549 template <class Type>
00550 Foam::label Foam::treeNode<Type>::setSubNodeType
00551 (
00552     const label level,
00553     octree<Type>& top,
00554     const Type& shapes
00555 )
00556 {
00557     if (debug & 4)
00558     {
00559         space(Pout, level);
00560         Pout<< "treeNode::setSubNodeType with level:" << level
00561             << "   bb:" << this->bb() << endl;
00562     }
00563 
00564     label myType = -1;
00565 
00566     for (label octant=0; octant<8; octant++)
00567     {
00568         label subType = -1;
00569 
00570         if (subNodes()[octant])
00571         {
00572             if (isNode(octant))
00573             {
00574                 subType = getNodePtr(octant)->setSubNodeType
00575                 (
00576                     level+1,
00577                     top,
00578                     shapes
00579                 );
00580             }
00581             else
00582             {
00583                 subType = getLeafPtr(octant)->setSubNodeType
00584                 (
00585                     level+1,
00586                     top,
00587                     shapes
00588                 );
00589             }
00590         }
00591         else
00592         {
00593             // No data in this one. Set type for octant acc. to its bounding
00594             // box.
00595             const treeBoundBox subBb = this->bb().subBbox(midpoint(), octant);
00596 
00597             subType = shapes.getSampleType(top, subBb.midpoint());
00598         }
00599 
00600         if (debug & 4)
00601         {
00602             space(Pout, level);
00603             Pout<< "treeNode::setSubNodeType : setting octant with bb:"
00604                 << this->bb().subBbox(midpoint(), octant)
00605                 << "  to type:" << octree<Type>::volType(subType) << endl;
00606         }
00607         setVolType(octant, subType);
00608 
00609         // Combine sub node types into type for treeNode. Result is 'mixed' if
00610         // types differ among subnodes.
00611         if (myType == -1)
00612         {
00613             myType = subType;
00614         }
00615         else if (subType != myType)
00616         {
00617             myType = octree<Type>::MIXED;
00618         }
00619     }
00620 
00621     if (debug & 4)
00622     {
00623         space(Pout, level);
00624         Pout<< "return from treeNode::setSubNodeType with type:"
00625             << octree<Type>::volType(myType)
00626             << "  bb:" << this->bb() << endl;
00627     }
00628 
00629     return myType;
00630 }
00631 
00632 
00633 // Get type of node.
00634 template <class Type>
00635 Foam::label Foam::treeNode<Type>::getSampleType
00636 (
00637     const label level,
00638     const octree<Type>& top,
00639     const Type& shapes,
00640     const point& sample
00641 ) const
00642 {
00643     if (debug & 4)
00644     {
00645         space(Pout, level);
00646         Pout<< "treeNode::getSampleType with level:" << level
00647             << " bb:" << this->bb() << "  sample:" << sample << endl;
00648     }
00649 
00650     // Determine octant of bb. If on edge just use whichever octant.
00651     bool onEdge = false;
00652 
00653     label octant = this->bb().subOctant(midpoint(), sample, onEdge);
00654 
00655     label type = getVolType(octant);
00656 
00657     if (type == octree<Type>::MIXED)
00658     {
00659         // At this level multiple sub types. Recurse to resolve.
00660         if (subNodes()[octant])
00661         {
00662             if (isNode(octant))
00663             {
00664                 // Node: recurse into subnodes
00665                 type = getNodePtr(octant)->getSampleType
00666                 (
00667                     level + 1,
00668                     top,
00669                     shapes,
00670                     sample
00671                 );
00672             }
00673             else
00674             {
00675                 // Leaf
00676                 type = getLeafPtr(octant)->getSampleType
00677                 (
00678                     level + 1,
00679                     top,
00680                     shapes,
00681                     sample
00682                 );
00683             }
00684         }
00685         else
00686         {
00687             // Problem: empty subnode should have a type
00688             FatalErrorIn
00689             (
00690                 "treeNode<Type>::getSampleType"
00691                 "(const label, octree<Type>&, const Type&, const point&)"
00692             )   << "Empty node bb:" << this->bb().subBbox(midpoint(), octant)
00693                 << " has non-mixed type:"
00694                 << octree<Type>::volType(type)
00695                 << abort(FatalError);
00696         }
00697     }
00698 
00699     if (type == octree<Type>::MIXED)
00700     {
00701         FatalErrorIn
00702         (
00703             "treeNode<Type>::getSampleType"
00704             "(const label, octree<Type>&, const Type&, const point&)"
00705         )   << "Type is MIXED when searching for " << sample
00706             << " at level " << this->bb() << endl
00707             << "This probably is because the octree has not been constructed"
00708             << " with search facility." << exit(FatalError);
00709     }
00710 
00711     if (debug & 4)
00712     {
00713         space(Pout, level);
00714         Pout<< "return from treeNode::getSampleType with type:"
00715             << octree<Type>::volType(type)
00716             << "  bb:" << this->bb()
00717             << "  sample:" << sample << endl;
00718     }
00719     return type;
00720 }
00721 
00722 
00723 template <class Type>
00724 Foam::label Foam::treeNode<Type>::find
00725 (
00726     const Type& shapes,
00727     const point& sample
00728 ) const
00729 {
00730     // Find octant of sample. Don't care if on edge (since any item on edge
00731     // will have been inserted in both subcubes)
00732     bool onEdge = false;
00733 
00734     label octant = this->bb().subOctant(midpoint(), sample, onEdge);
00735 
00736     if (subNodes()[octant])
00737     {
00738         if (isNode(octant))
00739         {
00740             // Node: recurse into subnodes
00741             return getNodePtr(octant)->find(shapes, sample);
00742         }
00743         else
00744         {
00745             // Leaf: let leaf::find handle this
00746             return getLeafPtr(octant)->find(shapes, sample);
00747         }
00748     }
00749     return -1;
00750 }
00751 
00752 
00753 template <class Type>
00754 bool Foam::treeNode<Type>::findTightest
00755 (
00756     const Type& shapes,
00757     const point& sample,
00758     treeBoundBox& tightest
00759 ) const
00760 {
00761     bool changed = false;
00762     bool onEdge = false;
00763     // Estimate for best place to start searching
00764     label sampleOctant = this->bb().subOctant(midpoint(), sample, onEdge);
00765 
00766     // Go into all suboctants (one containing sample first) and update tightest.
00767     // Order of visiting is if e.g. sampleOctant = 5:
00768     //  5 1 2 3 4 0 6 7
00769     for (label octantI=0; octantI<8; octantI++)
00770     {
00771         label octant;
00772         if (octantI == 0)
00773         {
00774             // Use sampleOctant first
00775             octant = sampleOctant;
00776         }
00777         else if (octantI == sampleOctant)
00778         {
00779             octant = 0;
00780         }
00781         else
00782         {
00783             octant = octantI;
00784         }
00785 
00786         if (subNodes()[octant])
00787         {
00788             if (isNode(octant))
00789             {
00790                 // Node: recurse into subnodes
00791                 const treeNode<Type>* subNodePtr = getNodePtr(octant);
00792 
00793                 if (subNodePtr->bb().overlaps(tightest))
00794                 {
00795                     // there might be a better fit inside this subNode
00796                     changed |= subNodePtr->findTightest
00797                     (
00798                         shapes,
00799                         sample,
00800                         tightest
00801                     );
00802                 }
00803             }
00804             else
00805             {
00806                 // Leaf: let leaf::find handle this
00807                 const treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
00808 
00809                 if (subLeafPtr->bb().overlaps(tightest))
00810                 {
00811                     // there might be a better fit inside this subLeaf
00812                     changed |= subLeafPtr->findTightest
00813                     (
00814                         shapes,
00815                         sample,
00816                         tightest
00817                     );
00818                 }
00819             }
00820         }
00821     }
00822 
00823     return changed;
00824 }
00825 
00826 
00827 template <class Type>
00828 bool Foam::treeNode<Type>::findNearest
00829 (
00830     const Type& shapes,
00831     const point& sample,
00832     treeBoundBox& tightest,
00833     label& tightestI,
00834     scalar& tightestDist
00835 ) const
00836 {
00837     if (debug & 8)
00838     {
00839         Pout<< "In findNearest with sample:" << sample << " cube:"
00840             << this->bb() << " tightest:" << tightest << endl;
00841     }
00842 
00843     bool changed = false;
00844     bool onEdge = false;
00845     // Estimate for best place to start searching
00846     label sampleOctant = this->bb().subOctant(midpoint(), sample, onEdge);
00847 
00848     // Go into all suboctants (one containing sample first) and update tightest.
00849     // Order of visiting is if e.g. sampleOctant = 5:
00850     //  5 1 2 3 4 0 6 7
00851     for (label octantI=0; octantI<8; octantI++)
00852     {
00853         label octant;
00854         if (octantI == 0)
00855         {
00856             // Use sampleOctant first
00857             octant = sampleOctant;
00858         }
00859         else if (octantI == sampleOctant)
00860         {
00861             octant = 0;
00862         }
00863         else
00864         {
00865             octant = octantI;
00866         }
00867 
00868         if (subNodes()[octant])
00869         {
00870             if (isNode(octant))
00871             {
00872                 // Node
00873                 const treeNode<Type>* subNodePtr = getNodePtr(octant);
00874 
00875                 if (subNodePtr->bb().overlaps(tightest))
00876                 {
00877                     // there might be a better fit inside this subNode
00878                     changed |= subNodePtr->findNearest
00879                     (
00880                         shapes,
00881                         sample,
00882                         tightest,
00883                         tightestI,
00884                         tightestDist
00885                     );
00886                 }
00887             }
00888             else
00889             {
00890                 // Leaf: let leaf::find handle this
00891                 const treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
00892 
00893                 if (subLeafPtr->bb().overlaps(tightest))
00894                 {
00895                     // there might be a better fit inside this subNode
00896                     changed |= subLeafPtr->findNearest
00897                     (
00898                         shapes,
00899                         sample,
00900                         tightest,
00901                         tightestI,
00902                         tightestDist
00903                     );
00904                 }
00905             }
00906         }
00907     }
00908 
00909     if (debug & 8)
00910     {
00911         Pout<< "Exiting findNearest for sample:" << sample << " cube:"
00912             << this->bb() << " tightestI:" << tightestI << endl;
00913     }
00914 
00915     return changed;
00916 }
00917 
00918 
00919 template <class Type>
00920 bool Foam::treeNode<Type>::findNearest
00921 (
00922     const Type& shapes,
00923     const linePointRef& ln,
00924     treeBoundBox& tightest,
00925     label& tightestI,
00926     point& linePoint,   // nearest point on line
00927     point& shapePoint   // nearest point on shape
00928 ) const
00929 {
00930     bool changed = false;
00931     bool onEdge = false;
00932     // Estimate for best place to start searching
00933     label sampleOctant = this->bb().subOctant(midpoint(), ln.centre(), onEdge);
00934 
00935     // Go into all suboctants (one containing sample first) and update tightest.
00936     // Order of visiting is if e.g. sampleOctant = 5:
00937     //  5 1 2 3 4 0 6 7
00938     for (label octantI=0; octantI<8; octantI++)
00939     {
00940         label octant;
00941         if (octantI == 0)
00942         {
00943             // Use sampleOctant first
00944             octant = sampleOctant;
00945         }
00946         else if (octantI == sampleOctant)
00947         {
00948             octant = 0;
00949         }
00950         else
00951         {
00952             octant = octantI;
00953         }
00954 
00955         if (subNodes()[octant])
00956         {
00957             if (isNode(octant))
00958             {
00959                 // Node
00960                 const treeNode<Type>* subNodePtr = getNodePtr(octant);
00961 
00962                 if (subNodePtr->bb().overlaps(tightest))
00963                 {
00964                     // there might be a better fit inside this subNode
00965                     changed |= subNodePtr->findNearest
00966                     (
00967                         shapes,
00968                         ln,
00969                         tightest,
00970                         tightestI,
00971                         linePoint,
00972                         shapePoint
00973                     );
00974                 }
00975             }
00976             else
00977             {
00978                 // Leaf: let leaf::find handle this
00979                 const treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
00980 
00981                 if (subLeafPtr->bb().overlaps(tightest))
00982                 {
00983                     // there might be a better fit inside this subNode
00984                     changed |= subLeafPtr->findNearest
00985                     (
00986                         shapes,
00987                         ln,
00988                         tightest,
00989                         tightestI,
00990                         linePoint,
00991                         shapePoint
00992                     );
00993                 }
00994             }
00995         }
00996     }
00997 
00998     return changed;
00999 }
01000 
01001 
01002 template <class Type>
01003 bool Foam::treeNode<Type>::findBox
01004 (
01005     const Type& shapes,
01006     const boundBox& box,
01007     labelHashSet& elements
01008 ) const
01009 {
01010     bool changed = false;
01011     bool onEdge = false;
01012     // Estimate for best place to start searching
01013     label sampleOctant = this->bb().subOctant
01014     (
01015         midpoint(),
01016         box.midpoint(),
01017         onEdge
01018     );
01019 
01020     // Go into all suboctants (one containing sample first) and update tightest.
01021     // Order of visiting is if e.g. sampleOctant = 5:
01022     //  5 1 2 3 4 0 6 7
01023     for (label octantI=0; octantI<8; octantI++)
01024     {
01025         label octant;
01026         if (octantI == 0)
01027         {
01028             // Use sampleOctant first
01029             octant = sampleOctant;
01030         }
01031         else if (octantI == sampleOctant)
01032         {
01033             octant = 0;
01034         }
01035         else
01036         {
01037             octant = octantI;
01038         }
01039 
01040         if (subNodes()[octant])
01041         {
01042             if (isNode(octant))
01043             {
01044                 // Node
01045                 const treeNode<Type>* subNodePtr = getNodePtr(octant);
01046 
01047                 if (subNodePtr->bb().overlaps(box))
01048                 {
01049                     // Visit sub node.
01050                     changed |= subNodePtr->findBox(shapes, box, elements);
01051                 }
01052             }
01053             else
01054             {
01055                 // Leaf: let leaf::find handle this
01056                 const treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
01057 
01058                 if (subLeafPtr->bb().overlaps(box))
01059                 {
01060                     // Visit sub leaf.
01061                     changed |= subLeafPtr->findBox(shapes, box, elements);
01062                 }
01063             }
01064         }
01065     }
01066 
01067     return changed;
01068 }
01069 
01070 
01071 // look from <start> in current cube (given by this->bb()).
01072 template <class Type>
01073 const Foam::treeLeaf<Type>* Foam::treeNode<Type>::findLeafLine
01074 (
01075     const int level,
01076     const Type& shapes,
01077     point& start,
01078     const point& end
01079 ) const
01080 {
01081     if (debug & 2)
01082     {
01083         space(Pout, 2*level);
01084         Pout<< "findLeafLine : bb:" << this->bb() << "  mid:" << midpoint()
01085             << "  start:" << start << endl;
01086     }
01087 
01088     scalar typDim = this->bb().avgDim();
01089 
01090     const vector direction = end - start;
01091 
01092     // Loop on current level until start has been updated to be outside
01093     // of this->bb(). Note that max only four subcubes can be crossed so this is
01094     // check on whether there are any truncation error problems.
01095 
01096     label iter = 0;
01097 
01098     while (true)
01099     {
01100         if (!this->bb().contains(direction, start))
01101         {
01102             if (debug & 2)
01103             {
01104                 space(Pout, 2*level);
01105                 Pout<< "findLeafLine : Start not inside bb " << this->bb()
01106                     << ". Returning with start:" << start << "  subLeaf:"
01107                     << 0 << endl;
01108             }
01109             return NULL;
01110         }
01111 
01112         // Check if start and <end> equal
01113         if ((mag(start - end)/typDim) < SMALL)
01114         {
01115             if (debug & 2)
01116             {
01117                 space(Pout, 2*level);
01118                 Pout<< "findLeafLine : start equals end"
01119                     << ". Returning with start:" << start << "  subLeaf:"
01120                     << 0 << endl;
01121             }
01122             return NULL;
01123         }
01124 
01125         if (iter >= 4)
01126         {
01127             // Too many iterations. Is hanging. Handle outside of loop.
01128             break;
01129         }
01130 
01131         bool onEdge = false;
01132         label octant = this->bb().subOctant
01133         (
01134             midpoint(), direction, start, onEdge
01135         );
01136 
01137         // Try finding non-empty treeleaf in octant
01138         const treeLeaf<Type>* leafPtr = findLeafLineOctant
01139         (
01140             level,
01141             shapes,
01142             octant,
01143             direction,
01144             start,
01145             end
01146         );
01147 
01148         if (leafPtr)
01149         {
01150             // Found treeLeaf -> return
01151             if (debug & 2)
01152             {
01153                 space(Pout, 2*level);
01154                 Pout<< "findLeafLine : Found treeLeaf"
01155                         << ". Returning with start:" << start << "  subLeaf:"
01156                         << long(leafPtr) << endl;
01157             }
01158 
01159             return leafPtr;
01160         }
01161 
01162         iter++;
01163     }
01164 
01165     // Check if is hanging. Max 4 octants can be crossed by a straight line
01166     FatalErrorIn
01167     (
01168         "treeNode<Type>::findLeafLine"
01169         "(const label, octree<Type>&, point&,"
01170         " const point&)"
01171     )   << "Did not leave bb " << this->bb()
01172         << " after " << iter
01173         << " iterations of updating starting point."
01174         << "start:" << start << "  end:" << end
01175         << abort(FatalError);
01176 
01177     return NULL;
01178 }
01179 
01180 
01181 template <class Type>
01182 void Foam::treeNode<Type>::findLeaves
01183 (
01184     List<treeLeaf<Type>*>& leafArray,
01185     label& leafIndex
01186 ) const
01187 {
01188     // Go into all sub boxes
01189     for (label octant=0; octant<8; octant++)
01190     {
01191         if (subNodes()[octant])
01192         {
01193             if (isNode(octant))
01194             {
01195                 // Node: recurse into subnodes
01196                 const treeNode<Type>* subNodePtr = getNodePtr(octant);
01197                 subNodePtr->findLeaves(leafArray, leafIndex);
01198             }
01199             else
01200             {
01201                 // Leaf: store
01202                 treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
01203                 leafArray[leafIndex++] = subLeafPtr;
01204             }
01205         }
01206     }
01207 }
01208 
01209 
01210 template <class Type>
01211 void Foam::treeNode<Type>::findLeaves
01212 (
01213     List<const treeLeaf<Type>*>& leafArray,
01214     label& leafIndex
01215 ) const
01216 {
01217     // Go into all sub boxes
01218     for (label octant=0; octant<8; octant++)
01219     {
01220         if (subNodes()[octant])
01221         {
01222             if (isNode(octant))
01223             {
01224                 // Node: recurse into subnodes
01225                 const treeNode<Type>* subNodePtr = getNodePtr(octant);
01226                 subNodePtr->findLeaves(leafArray, leafIndex);
01227             }
01228             else
01229             {
01230                 // Leaf: store
01231                 treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
01232                 leafArray[leafIndex++] = subLeafPtr;
01233             }
01234         }
01235     }
01236 }
01237 
01238 
01239 template <class Type>
01240 void Foam::treeNode<Type>::printNode
01241 (
01242     Ostream& os,
01243     const label level
01244 ) const
01245 {
01246     space(os, 2*level);
01247 
01248     os << "node:" << this->bb() << endl;
01249 
01250     for (label octant=0; octant<8; octant++)
01251     {
01252         label type = getVolType(octant);
01253 
01254         string typeString = octree<Type>::volType(type);
01255 
01256         if (!subNodes_[octant])
01257         {
01258             space(os, level);
01259             os << octant << ":" << typeString << " : null" << endl;
01260         }
01261         else if (isNode(octant))
01262         {
01263             space(os, level);
01264             os << octant << ":" << typeString << " : node" << endl;
01265             getNodePtr(octant)->printNode(os, level+1);
01266         }
01267         else
01268         {
01269             space(os, level);
01270             os << octant << ":" << typeString << " : leaf" << endl;
01271 
01272             treeLeaf<Type>* leafPtr = getLeafPtr(octant);
01273             leafPtr->printLeaf(os, level+1);
01274         }
01275     }
01276 }
01277 
01278 
01279 template <class Type>
01280 void Foam::treeNode<Type>::writeOBJ
01281 (
01282     Ostream& os,
01283     const label level,
01284     label& vertNo
01285 ) const
01286 {
01287     point midPoint(this->bb().midpoint());
01288 
01289     label midVertNo = vertNo;
01290     os << "v " << midPoint.x() << " " << midPoint.y() << " "
01291        << midPoint.z() << endl;
01292     vertNo++;
01293 
01294     for (label octant=0; octant<8; octant++)
01295     {
01296         if (subNodes_[octant])
01297         {
01298             if (isNode(octant))
01299             {
01300                 treeNode<Type>* nodePtr = getNodePtr(octant);
01301 
01302                 point subMidPoint(nodePtr->bb().midpoint());
01303                 os << "v " << subMidPoint.x() << " " << subMidPoint.y() << " "
01304                    << subMidPoint.z() << endl;
01305                 os << "l " << midVertNo + 1<< " " << vertNo + 1 << endl;
01306                 vertNo++;
01307 
01308                 nodePtr->writeOBJ(os, level+1, vertNo);
01309             }
01310             else
01311             {
01312                 treeLeaf<Type>* leafPtr = getLeafPtr(octant);
01313 
01314                 point subMidPoint(leafPtr->bb().midpoint());
01315                 os << "v " << subMidPoint.x() << " " << subMidPoint.y() << " "
01316                    << subMidPoint.z() << endl;
01317                 os << "l " << midVertNo + 1<< " " << vertNo + 1 << endl;
01318                 vertNo++;
01319 
01320                 //leafPtr->writeOBJ(os, level+1, vertNo);
01321             }
01322         }
01323     }
01324 }
01325 
01326 
01327 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
01328 
01329 template <class Type>
01330 Foam::Istream& Foam::operator>>(Istream& is, treeNode<Type>& oc)
01331 {
01332     for (label octant = 0; octant < 8; octant++)
01333     {
01334         oc.subNodes_[octant] = NULL;
01335     }
01336 
01337     is >> oc.bb();
01338 
01339     label nPtrs;
01340 
01341     // Read number of entries folllowing
01342     is >> nPtrs;
01343 
01344     is.readBegin("treeNode");
01345     for (label octant = 0; octant < nPtrs; octant++)
01346     {
01347         label index;
01348         is >> index;
01349 
01350         if (index >= treeNode<Type>::leafOffset)
01351         {
01352             // Leaf recognized by out of range index
01353             treeLeaf<Type>* leafPtr = new treeLeaf<Type>(is);
01354             oc.setLeafPtr(index - treeNode<Type>::leafOffset, leafPtr);
01355         }
01356         else
01357         {
01358             oc.setNodePtr(index, new treeNode<Type>(is));
01359         }
01360     }
01361 
01362     // Read end of treeNode list
01363     is.readEnd("treeNode");
01364 
01365     // Check state of Istream
01366     is.check("Istream& operator>>(Istream&, treeNode&)");
01367 
01368     return is;
01369 }
01370 
01371 
01372 template <class Type>
01373 Foam::Ostream& Foam::operator<<(Ostream& os, const treeNode<Type>& tn)
01374 {
01375     // Count valid subnodes:
01376     //   - treeNode
01377     //   - treeLeafs with non-zero cell list.
01378     label nPtrs = 0;
01379     for (label octant = 0; octant < 8; octant++)
01380     {
01381         if (tn.subNodes_[octant])
01382         {
01383             if (tn.isNode(octant) || tn.getLeafPtr(octant)->indices().size())
01384             {
01385                 nPtrs++;
01386             }
01387         }
01388     }
01389 
01390 
01391     // output subnodes as list of length nPtrs
01392     os << token::SPACE << tn.bb() << token::SPACE << nPtrs
01393        << token::SPACE << token::BEGIN_LIST;
01394 
01395     for (label octant = 0; octant < 8; octant++)
01396     {
01397         if (tn.subNodes_[octant])
01398         {
01399             if (tn.isNode(octant))
01400             {
01401                 const treeNode<Type>* subNodePtr = tn.getNodePtr(octant);
01402 
01403                 // Node: output index, value
01404                 os  << token::SPACE << octant << token::SPACE << *subNodePtr
01405                     << token::NL;
01406             }
01407             else if (tn.getLeafPtr(octant)->indices().size())
01408             {
01409                 // treeLeaf: mark by putting index invalid
01410                 const treeLeaf<Type>* subLeafPtr = tn.getLeafPtr(octant);
01411 
01412                 os  << token::SPACE << octant + treeNode<Type>::leafOffset
01413                     << token::SPACE << *subLeafPtr
01414                     << token::NL;
01415             }
01416         }
01417     }
01418 
01419     os  << token::SPACE << token::END_LIST;
01420 
01421     return os;
01422 }
01423 
01424 
01425 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines