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

indexedOctree.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 \*---------------------------------------------------------------------------*/
00025 
00026 #include "indexedOctree.H"
00027 #include <OpenFOAM/linePointRef.H>
00028 #include <meshTools/meshTools.H>
00029 #include <OpenFOAM/OFstream.H>
00030 
00031 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00032 
00033 template <class Type>
00034 Foam::scalar Foam::indexedOctree<Type>::perturbTol_ = 10*SMALL;
00035 
00036 
00037 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00038 
00039 // Does bb intersect a sphere around sample? Or is any corner point of bb
00040 // closer than nearestDistSqr to sample.
00041 template <class Type>
00042 bool Foam::indexedOctree<Type>::overlaps
00043 (
00044     const point& p0,
00045     const point& p1,
00046     const scalar nearestDistSqr,
00047     const point& sample
00048 )
00049 {
00050     // Find out where sample is in relation to bb.
00051     // Find nearest point on bb.
00052     scalar distSqr = 0;
00053 
00054     for (direction dir = 0; dir < vector::nComponents; dir++)
00055     {
00056         scalar d0 = p0[dir] - sample[dir];
00057         scalar d1 = p1[dir] - sample[dir];
00058 
00059         if ((d0 > 0) != (d1 > 0))
00060         {
00061             // sample inside both extrema. This component does not add any
00062             // distance.
00063         }
00064         else if (mag(d0) < mag(d1))
00065         {
00066             distSqr += d0*d0;
00067         }
00068         else
00069         {
00070             distSqr += d1*d1;
00071         }
00072 
00073         if (distSqr > nearestDistSqr)
00074         {
00075             return false;
00076         }
00077     }
00078 
00079     return true;
00080 }
00081 
00082 
00083 // Does bb intersect a sphere around sample? Or is any corner point of bb
00084 // closer than nearestDistSqr to sample.
00085 template <class Type>
00086 bool Foam::indexedOctree<Type>::overlaps
00087 (
00088     const treeBoundBox& parentBb,
00089     const direction octant,
00090     const scalar nearestDistSqr,
00091     const point& sample
00092 )
00093 {
00094     //- Accelerated version of
00095     //     treeBoundBox subBb(parentBb.subBbox(mid, octant))
00096     //     overlaps
00097     //     (
00098     //          subBb.min(),
00099     //          subBb.max(),
00100     //          nearestDistSqr,
00101     //          sample
00102     //     )
00103 
00104     const point& min = parentBb.min();
00105     const point& max = parentBb.max();
00106 
00107     point other;
00108 
00109     if (octant & treeBoundBox::RIGHTHALF)
00110     {
00111         other.x() = max.x();
00112     }
00113     else
00114     {
00115         other.x() = min.x();
00116     }
00117 
00118     if (octant & treeBoundBox::TOPHALF)
00119     {
00120         other.y() = max.y();
00121     }
00122     else
00123     {
00124         other.y() = min.y();
00125     }
00126 
00127     if (octant & treeBoundBox::FRONTHALF)
00128     {
00129         other.z() = max.z();
00130     }
00131     else
00132     {
00133         other.z() = min.z();
00134     }
00135 
00136     const point mid(0.5*(min+max));
00137 
00138     return overlaps(mid, other, nearestDistSqr, sample);
00139 }
00140 
00141 
00142 //
00143 // Construction helper routines
00144 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00145 //
00146 
00147 // Split list of indices into 8 bins
00148 template <class Type>
00149 void Foam::indexedOctree<Type>::divide
00150 (
00151     const labelList& indices,
00152     const treeBoundBox& bb,
00153     labelListList& result
00154 ) const
00155 {
00156     List<DynamicList<label> > subIndices(8);
00157     for (direction octant = 0; octant < subIndices.size(); octant++)
00158     {
00159         subIndices[octant].setCapacity(indices.size()/8);
00160     }
00161 
00162     // Precalculate bounding boxes.
00163     FixedList<treeBoundBox, 8> subBbs;
00164     for (direction octant = 0; octant < subBbs.size(); octant++)
00165     {
00166         subBbs[octant] = bb.subBbox(octant);
00167     }
00168 
00169     forAll(indices, i)
00170     {
00171         label shapeI = indices[i];
00172 
00173         for (direction octant = 0; octant < 8; octant++)
00174         {
00175             if (shapes_.overlaps(shapeI, subBbs[octant]))
00176             {
00177                 subIndices[octant].append(shapeI);
00178             }
00179         }
00180     }
00181 
00182     result.setSize(8);
00183     for (direction octant = 0; octant < subIndices.size(); octant++)
00184     {
00185         result[octant].transfer(subIndices[octant]);
00186     }
00187 }
00188 
00189 
00190 // Subdivide the (content) node.
00191 template <class Type>
00192 typename Foam::indexedOctree<Type>::node
00193 Foam::indexedOctree<Type>::divide
00194 (
00195     const treeBoundBox& bb,
00196     DynamicList<labelList>& contents,
00197     const label contentI
00198 ) const
00199 {
00200     const labelList& indices = contents[contentI];
00201 
00202     node nod;
00203 
00204     if
00205     (
00206         bb.min()[0] >= bb.max()[0]
00207      || bb.min()[1] >= bb.max()[1]
00208      || bb.min()[2] >= bb.max()[2]
00209     )
00210     {
00211         FatalErrorIn("indexedOctree<Type>::divide(..)")
00212             << "Badly formed bounding box:" << bb
00213             << abort(FatalError);
00214     }
00215 
00216     nod.bb_ = bb;
00217     nod.parent_ = -1;
00218 
00219     labelListList dividedIndices(8);
00220     divide(indices, bb, dividedIndices);
00221 
00222     // Have now divided the indices into 8 (possibly empty) subsets.
00223     // Replace current contentI with the first (non-empty) subset.
00224     // Append the rest.
00225     bool replaced = false;
00226 
00227     for (direction octant = 0; octant < dividedIndices.size(); octant++)
00228     {
00229         labelList& subIndices = dividedIndices[octant];
00230 
00231         if (subIndices.size())
00232         {
00233             if (!replaced)
00234             {
00235                 contents[contentI].transfer(subIndices);
00236                 nod.subNodes_[octant] = contentPlusOctant(contentI, octant);
00237                 replaced = true;
00238             }
00239             else
00240             {
00241                 // Store at end of contents.
00242                 // note dummy append + transfer trick
00243                 label sz = contents.size();
00244                 contents.append(labelList(0));
00245                 contents[sz].transfer(subIndices);
00246                 nod.subNodes_[octant] = contentPlusOctant(sz, octant);
00247             }
00248         }
00249         else
00250         {
00251             // Mark node as empty
00252             nod.subNodes_[octant] = emptyPlusOctant(octant);
00253         }
00254     }
00255 
00256     return nod;
00257 }
00258 
00259 
00260 // Split any contents node with more than minSize elements.
00261 template <class Type>
00262 void Foam::indexedOctree<Type>::splitNodes
00263 (
00264     const label minSize,
00265     DynamicList<indexedOctree<Type>::node>& nodes,
00266     DynamicList<labelList>& contents
00267 ) const
00268 {
00269     label currentSize = nodes.size();
00270 
00271     // Take care to loop only over old nodes.
00272     // Also we loop over the same DynamicList which gets modified and
00273     // moved so make sure not to keep any references!
00274     for (label nodeI = 0; nodeI < currentSize; nodeI++)
00275     {
00276         for
00277         (
00278             direction octant = 0;
00279             octant < nodes[nodeI].subNodes_.size();
00280             octant++
00281         )
00282         {
00283             labelBits index = nodes[nodeI].subNodes_[octant];
00284 
00285             if (isNode(index))
00286             {
00287                 // tree node. Leave intact.
00288             }
00289             else if (isContent(index))
00290             {
00291                 label contentI = getContent(index);
00292 
00293                 if (contents[contentI].size() > minSize)
00294                 {
00295                     // Create node for content.
00296 
00297                     // Find the bounding box for the subnode
00298                     const node& nod = nodes[nodeI];
00299                     const treeBoundBox bb(nod.bb_.subBbox(octant));
00300 
00301                     node subNode(divide(bb, contents, contentI));
00302                     subNode.parent_ = nodeI;
00303                     label sz = nodes.size();
00304                     nodes.append(subNode);
00305                     nodes[nodeI].subNodes_[octant] = nodePlusOctant(sz, octant);
00306                 }
00307             }
00308         }
00309     }
00310 }
00311 
00312 
00313 // Reorder contents to be in same order as nodes. Returns number of nodes on
00314 // the compactLevel.
00315 template <class Type>
00316 Foam::label Foam::indexedOctree<Type>::compactContents
00317 (
00318     DynamicList<node>& nodes,
00319     DynamicList<labelList>& contents,
00320     const label compactLevel,
00321     const label nodeI,
00322     const label level,
00323 
00324     List<labelList>& compactedContents,
00325     label& compactI
00326 )
00327 {
00328     const node& nod = nodes[nodeI];
00329 
00330     label nNodes = 0;
00331 
00332     if (level < compactLevel)
00333     {
00334         for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
00335         {
00336             labelBits index = nod.subNodes_[octant];
00337 
00338             if (isNode(index))
00339             {
00340                 nNodes += compactContents
00341                 (
00342                     nodes,
00343                     contents,
00344                     compactLevel,
00345                     getNode(index),
00346                     level+1,
00347                     compactedContents,
00348                     compactI
00349                 );
00350             }
00351         }
00352     }
00353     else if (level == compactLevel)
00354     {
00355         // Compact all content on this level
00356         for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
00357         {
00358             labelBits index = nod.subNodes_[octant];
00359 
00360             if (isContent(index))
00361             {
00362                 label contentI = getContent(index);
00363 
00364                 compactedContents[compactI].transfer(contents[contentI]);
00365 
00366                 // Subnode is at compactI. Adapt nodeI to point to it
00367                 nodes[nodeI].subNodes_[octant] =
00368                     contentPlusOctant(compactI, octant);
00369 
00370                 compactI++;
00371             }
00372             else if (isNode(index))
00373             {
00374                 nNodes++;
00375             }
00376         }
00377     }
00378     return nNodes;
00379 }
00380 
00381 
00382 // Pre-calculates wherever possible the volume status per node/subnode.
00383 // Recurses to determine status of lowest level boxes. Level above is
00384 // combination of octants below.
00385 template <class Type>
00386 typename Foam::indexedOctree<Type>::volumeType
00387 Foam::indexedOctree<Type>::calcVolumeType
00388 (
00389     const label nodeI
00390 ) const
00391 {
00392     const node& nod = nodes_[nodeI];
00393 
00394     volumeType myType = UNKNOWN;
00395 
00396     for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
00397     {
00398         volumeType subType;
00399 
00400         labelBits index = nod.subNodes_[octant];
00401 
00402         if (isNode(index))
00403         {
00404             // tree node. Recurse.
00405             subType = calcVolumeType(getNode(index));
00406         }
00407         else if (isContent(index))
00408         {
00409             // Contents. Depending on position in box might be on either
00410             // side.
00411             subType = MIXED;
00412         }
00413         else
00414         {
00415             // No data in this octant. Set type for octant acc. to the mid
00416             // of its bounding box.
00417             const treeBoundBox subBb = nod.bb_.subBbox(octant);
00418 
00419             subType = volumeType
00420             (
00421                 shapes_.getVolumeType(*this, subBb.midpoint())
00422             );
00423         }
00424 
00425         // Store octant type
00426         nodeTypes_.set((nodeI<<3)+octant, subType);
00427 
00428         // Combine sub node types into type for treeNode. Result is 'mixed' if
00429         // types differ among subnodes.
00430         if (myType == UNKNOWN)
00431         {
00432             myType = subType;
00433         }
00434         else if (subType != myType)
00435         {
00436             myType = MIXED;
00437         }
00438     }
00439     return myType;
00440 }
00441 
00442 
00443 template <class Type>
00444 typename Foam::indexedOctree<Type>::volumeType
00445 Foam::indexedOctree<Type>::getVolumeType
00446 (
00447     const label nodeI,
00448     const point& sample
00449 ) const
00450 {
00451     const node& nod = nodes_[nodeI];
00452 
00453     direction octant = nod.bb_.subOctant(sample);
00454 
00455     volumeType octantType = volumeType(nodeTypes_.get((nodeI<<3)+octant));
00456 
00457     if (octantType == INSIDE)
00458     {
00459         return octantType;
00460     }
00461     else if (octantType == OUTSIDE)
00462     {
00463         return octantType;
00464     }
00465     else if (octantType == UNKNOWN)
00466     {
00467         // Can happen for e.g. non-manifold surfaces.
00468         return octantType;
00469     }
00470     else if (octantType == MIXED)
00471     {
00472         labelBits index = nod.subNodes_[octant];
00473 
00474         if (isNode(index))
00475         {
00476             // Recurse
00477             volumeType subType = getVolumeType(getNode(index), sample);
00478 
00479             return subType;
00480         }
00481         else if (isContent(index))
00482         {
00483             // Content. Defer to shapes.
00484             return volumeType(shapes_.getVolumeType(*this, sample));
00485         }
00486         else
00487         {
00488             // Empty node. Cannot have 'mixed' as its type since not divided
00489             // up and has no items inside it.
00490             FatalErrorIn
00491             (
00492                 "indexedOctree<Type>::getVolumeType"
00493                 "(const label, const point&)"
00494             )   << "Sample:" << sample << " node:" << nodeI
00495                 << " with bb:" << nodes_[nodeI].bb_ << nl
00496                 << "Empty subnode has invalid volume type MIXED."
00497                 << abort(FatalError);
00498 
00499             return UNKNOWN;
00500         }
00501     }
00502     else
00503     {
00504         FatalErrorIn
00505         (
00506             "indexedOctree<Type>::getVolumeType"
00507             "(const label, const point&)"
00508         )   << "Sample:" << sample << " at node:" << nodeI
00509             << " octant:" << octant
00510             << " with bb:" << nod.bb_.subBbox(octant) << nl
00511             << "Node has invalid volume type " << octantType
00512             << abort(FatalError);
00513 
00514         return UNKNOWN;
00515     }
00516 }
00517 
00518 
00519 template <class Type>
00520 typename Foam::indexedOctree<Type>::volumeType
00521 Foam::indexedOctree<Type>::getSide
00522 (
00523     const vector& outsideNormal,
00524     const vector& vec
00525 )
00526 {
00527     if ((outsideNormal&vec) >= 0)
00528     {
00529         return OUTSIDE;
00530     }
00531     else
00532     {
00533         return INSIDE;
00534     }
00535 }
00536 
00537 
00538 //
00539 // Query routines
00540 // ~~~~~~~~~~~~~~
00541 //
00542 
00543 // Find nearest point starting from nodeI
00544 template <class Type>
00545 void Foam::indexedOctree<Type>::findNearest
00546 (
00547     const label nodeI,
00548     const point& sample,
00549 
00550     scalar& nearestDistSqr,
00551     label& nearestShapeI,
00552     point& nearestPoint
00553 ) const
00554 {
00555     const node& nod = nodes_[nodeI];
00556 
00557     // Determine order to walk through octants
00558     FixedList<direction, 8> octantOrder;
00559     nod.bb_.searchOrder(sample, octantOrder);
00560 
00561     // Go into all suboctants (one containing sample first) and update nearest.
00562     for (direction i = 0; i < 8; i++)
00563     {
00564         direction octant = octantOrder[i];
00565 
00566         labelBits index = nod.subNodes_[octant];
00567 
00568         if (isNode(index))
00569         {
00570             label subNodeI = getNode(index);
00571 
00572             const treeBoundBox& subBb = nodes_[subNodeI].bb_;
00573 
00574             if (overlaps(subBb.min(), subBb.max(), nearestDistSqr, sample))
00575             {
00576                 findNearest
00577                 (
00578                     subNodeI,
00579                     sample,
00580 
00581                     nearestDistSqr,
00582                     nearestShapeI,
00583                     nearestPoint
00584                 );
00585             }
00586         }
00587         else if (isContent(index))
00588         {
00589             if
00590             (
00591                 overlaps
00592                 (
00593                     nod.bb_,
00594                     octant,
00595                     nearestDistSqr,
00596                     sample
00597                 )
00598             )
00599             {
00600                 shapes_.findNearest
00601                 (
00602                     contents_[getContent(index)],
00603                     sample,
00604 
00605                     nearestDistSqr,
00606                     nearestShapeI,
00607                     nearestPoint
00608                 );
00609             }
00610         }
00611     }
00612 }
00613 
00614 
00615 // Find nearest point to line.
00616 template <class Type>
00617 void Foam::indexedOctree<Type>::findNearest
00618 (
00619     const label nodeI,
00620     const linePointRef& ln,
00621 
00622     treeBoundBox& tightest,
00623     label& nearestShapeI,
00624     point& linePoint,
00625     point& nearestPoint
00626 ) const
00627 {
00628     const node& nod = nodes_[nodeI];
00629     const treeBoundBox& nodeBb = nod.bb_;
00630 
00631     // Determine order to walk through octants
00632     FixedList<direction, 8> octantOrder;
00633     nod.bb_.searchOrder(ln.centre(), octantOrder);
00634 
00635     // Go into all suboctants (one containing sample first) and update nearest.
00636     for (direction i = 0; i < 8; i++)
00637     {
00638         direction octant = octantOrder[i];
00639 
00640         labelBits index = nod.subNodes_[octant];
00641 
00642         if (isNode(index))
00643         {
00644             const treeBoundBox& subBb = nodes_[getNode(index)].bb_;
00645 
00646             if (subBb.overlaps(tightest))
00647             {
00648                 findNearest
00649                 (
00650                     getNode(index),
00651                     ln,
00652 
00653                     tightest,
00654                     nearestShapeI,
00655                     linePoint,
00656                     nearestPoint
00657                 );
00658             }
00659         }
00660         else if (isContent(index))
00661         {
00662             const treeBoundBox subBb(nodeBb.subBbox(octant));
00663 
00664             if (subBb.overlaps(tightest))
00665             {
00666                 shapes_.findNearest
00667                 (
00668                     contents_[getContent(index)],
00669                     ln,
00670 
00671                     tightest,
00672                     nearestShapeI,
00673                     linePoint,
00674                     nearestPoint
00675                 );
00676             }
00677         }
00678     }
00679 }
00680 
00681 
00682 template <class Type>
00683 Foam::treeBoundBox Foam::indexedOctree<Type>::subBbox
00684 (
00685     const label parentNodeI,
00686     const direction octant
00687 ) const
00688 {
00689     // Get type of node at octant
00690     const node& nod = nodes_[parentNodeI];
00691     labelBits index = nod.subNodes_[octant];
00692 
00693     if (isNode(index))
00694     {
00695         // Use stored bb
00696         return nodes_[getNode(index)].bb_;
00697     }
00698     else
00699     {
00700         // Calculate subBb
00701         return nod.bb_.subBbox(octant);
00702     }
00703 }
00704 
00705 
00706 // Takes a bb and a point on/close to the edge of the bb and pushes the point
00707 // inside by a small fraction.
00708 template <class Type>
00709 Foam::point Foam::indexedOctree<Type>::pushPoint
00710 (
00711     const treeBoundBox& bb,
00712     const point& pt,
00713     const bool pushInside
00714 )
00715 {
00716     // Get local length scale.
00717     const vector perturbVec = perturbTol_*(bb.span());
00718 
00719     point perturbedPt(pt);
00720 
00721     // Modify all components which are close to any face of the bb to be
00722     // well inside/outside them.
00723 
00724     if (pushInside)
00725     {
00726         for (direction dir = 0; dir < vector::nComponents; dir++)
00727         {
00728             if (mag(pt[dir]-bb.min()[dir]) < mag(perturbVec[dir]))
00729             {
00730                 // Close to 'left' side. Push well beyond left side.
00731                 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
00732                 perturbedPt[dir] = bb.min()[dir] + perturbDist;
00733             }
00734             else if (mag(pt[dir]-bb.max()[dir]) < mag(perturbVec[dir]))
00735             {
00736                 // Close to 'right' side. Push well beyond right side.
00737                 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
00738                 perturbedPt[dir] = bb.max()[dir] - perturbDist;
00739             }
00740         }
00741     }
00742     else
00743     {
00744         for (direction dir = 0; dir < vector::nComponents; dir++)
00745         {
00746             if (mag(pt[dir]-bb.min()[dir]) < mag(perturbVec[dir]))
00747             {
00748                 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
00749                 perturbedPt[dir] = bb.min()[dir] - perturbDist;
00750             }
00751             else if (mag(pt[dir]-bb.max()[dir]) < mag(perturbVec[dir]))
00752             {
00753                 scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
00754                 perturbedPt[dir] = bb.max()[dir] + perturbDist;
00755             }
00756         }
00757     }
00758 
00759     if (debug)
00760     {
00761         if (pushInside != bb.contains(perturbedPt))
00762         {
00763             FatalErrorIn("indexedOctree<Type>::pushPoint(..)")
00764                 << "pushed point:" << pt
00765                 << " to:" << perturbedPt
00766                 << " wanted side:" << pushInside
00767                 << " obtained side:" << bb.contains(perturbedPt)
00768                 << " of bb:" << bb
00769                 << abort(FatalError);
00770         }
00771     }
00772 
00773     return perturbedPt;
00774 }
00775 
00776 
00777 // Takes a bb and a point on the edge of the bb and pushes the point
00778 // outside by a small fraction.
00779 template <class Type>
00780 Foam::point Foam::indexedOctree<Type>::pushPoint
00781 (
00782     const treeBoundBox& bb,
00783     const direction faceID,
00784     const point& pt,
00785     const bool pushInside
00786 )
00787 {
00788     // Get local length scale.
00789     const vector perturbVec = perturbTol_*bb.span();
00790 
00791     point perturbedPt(pt);
00792 
00793     // Modify all components which are close to any face of the bb to be
00794     // well outside them.
00795 
00796     if (faceID == 0)
00797     {
00798         FatalErrorIn("indexedOctree<Type>::pushPoint(..)")
00799             << abort(FatalError);
00800     }
00801 
00802     if (faceID & treeBoundBox::LEFTBIT)
00803     {
00804         if (pushInside)
00805         {
00806             perturbedPt[0] = bb.min()[0] + (perturbVec[0] + ROOTVSMALL);
00807         }
00808         else
00809         {
00810             perturbedPt[0] = bb.min()[0] - (perturbVec[0] + ROOTVSMALL);
00811         }
00812     }
00813     else if (faceID & treeBoundBox::RIGHTBIT)
00814     {
00815         if (pushInside)
00816         {
00817             perturbedPt[0] = bb.max()[0] - (perturbVec[0] + ROOTVSMALL);
00818         }
00819         else
00820         {
00821             perturbedPt[0] = bb.max()[0] + (perturbVec[0] + ROOTVSMALL);
00822         }
00823     }
00824 
00825     if (faceID & treeBoundBox::BOTTOMBIT)
00826     {
00827         if (pushInside)
00828         {
00829             perturbedPt[1] = bb.min()[1] + (perturbVec[1] + ROOTVSMALL);
00830         }
00831         else
00832         {
00833             perturbedPt[1] = bb.min()[1] - (perturbVec[1] + ROOTVSMALL);
00834         }
00835     }
00836     else if (faceID & treeBoundBox::TOPBIT)
00837     {
00838         if (pushInside)
00839         {
00840             perturbedPt[1] = bb.max()[1] - (perturbVec[1] + ROOTVSMALL);
00841         }
00842         else
00843         {
00844             perturbedPt[1] = bb.max()[1] + (perturbVec[1] + ROOTVSMALL);
00845         }
00846     }
00847 
00848     if (faceID & treeBoundBox::BACKBIT)
00849     {
00850         if (pushInside)
00851         {
00852             perturbedPt[2] = bb.min()[2] + (perturbVec[2] + ROOTVSMALL);
00853         }
00854         else
00855         {
00856             perturbedPt[2] = bb.min()[2] - (perturbVec[2] + ROOTVSMALL);
00857         }
00858     }
00859     else if (faceID & treeBoundBox::FRONTBIT)
00860     {
00861         if (pushInside)
00862         {
00863             perturbedPt[2] = bb.max()[2] - (perturbVec[2] + ROOTVSMALL);
00864         }
00865         else
00866         {
00867             perturbedPt[2] = bb.max()[2] + (perturbVec[2] + ROOTVSMALL);
00868         }
00869     }
00870 
00871     if (debug)
00872     {
00873         if (pushInside != bb.contains(perturbedPt))
00874         {
00875             FatalErrorIn("indexedOctree<Type>::pushPoint(..)")
00876                 << "pushed point:" << pt << " on face:" << faceString(faceID)
00877                 << " to:" << perturbedPt
00878                 << " wanted side:" << pushInside
00879                 << " obtained side:" << bb.contains(perturbedPt)
00880                 << " of bb:" << bb
00881                 << abort(FatalError);
00882         }
00883     }
00884 
00885     return perturbedPt;
00886 }
00887 
00888 
00889 // Guarantees that if pt is on a face it gets perturbed so it is away
00890 // from the face edges.
00891 // If pt is not on a face does nothing.
00892 template <class Type>
00893 Foam::point Foam::indexedOctree<Type>::pushPointIntoFace
00894 (
00895     const treeBoundBox& bb,
00896     const vector& dir,          // end-start
00897     const point& pt
00898 )
00899 {
00900     if (debug)
00901     {
00902         if (bb.posBits(pt) != 0)
00903         {
00904             FatalErrorIn("indexedOctree<Type>::pushPointIntoFace(..)")
00905                 << " bb:" << bb << endl
00906                 << "does not contain point " << pt << abort(FatalError);
00907         }
00908     }
00909 
00910 
00911     // Handle two cases:
00912     // - point exactly on multiple faces. Push away from all but one.
00913     // - point on a single face. Push away from edges of face.
00914 
00915     direction ptFaceID = bb.faceBits(pt);
00916 
00917     direction nFaces = 0;
00918     FixedList<direction, 3> faceIndices;
00919 
00920     if (ptFaceID & treeBoundBox::LEFTBIT)
00921     {
00922         faceIndices[nFaces++] = treeBoundBox::LEFT;
00923     }
00924     else if (ptFaceID & treeBoundBox::RIGHTBIT)
00925     {
00926         faceIndices[nFaces++] = treeBoundBox::RIGHT;
00927     }
00928 
00929     if (ptFaceID & treeBoundBox::BOTTOMBIT)
00930     {
00931         faceIndices[nFaces++] = treeBoundBox::BOTTOM;
00932     }
00933     else if (ptFaceID & treeBoundBox::TOPBIT)
00934     {
00935         faceIndices[nFaces++] = treeBoundBox::TOP;
00936     }
00937 
00938     if (ptFaceID & treeBoundBox::BACKBIT)
00939     {
00940         faceIndices[nFaces++] = treeBoundBox::BACK;
00941     }
00942     else if (ptFaceID & treeBoundBox::FRONTBIT)
00943     {
00944         faceIndices[nFaces++] = treeBoundBox::FRONT;
00945     }
00946 
00947 
00948     // Determine the face we want to keep the point on
00949 
00950     direction keepFaceID;
00951 
00952     if (nFaces == 0)
00953     {
00954         // Return original point
00955         return pt;
00956     }
00957     else if (nFaces == 1)
00958     {
00959         // Point is on a single face
00960         keepFaceID = faceIndices[0];
00961     }
00962     else
00963     {
00964         // Determine best face out of faceIndices[0 .. nFaces-1].
00965         // The best face is the one most perpendicular to the ray direction.
00966 
00967         keepFaceID = faceIndices[0];
00968         scalar maxInproduct = mag(treeBoundBox::faceNormals[keepFaceID] & dir);
00969 
00970         for (direction i = 1; i < nFaces; i++)
00971         {
00972             direction face = faceIndices[i];
00973             scalar s = mag(treeBoundBox::faceNormals[face] & dir);
00974             if (s > maxInproduct)
00975             {
00976                 maxInproduct = s;
00977                 keepFaceID = face;
00978             }
00979         }
00980     }
00981 
00982 
00983     // 1. Push point into bb, away from all corners
00984 
00985     point facePoint(pushPoint(bb, pt, true));
00986     direction faceID = 0;
00987 
00988     // 2. Snap it back onto the preferred face
00989 
00990     if (keepFaceID == treeBoundBox::LEFT)
00991     {
00992         facePoint.x() = bb.min().x();
00993         faceID = treeBoundBox::LEFTBIT;
00994     }
00995     else if (keepFaceID == treeBoundBox::RIGHT)
00996     {
00997         facePoint.x() = bb.max().x();
00998         faceID = treeBoundBox::RIGHTBIT;
00999     }
01000     else if (keepFaceID == treeBoundBox::BOTTOM)
01001     {
01002         facePoint.y() = bb.min().y();
01003         faceID = treeBoundBox::BOTTOMBIT;
01004     }
01005     else if (keepFaceID == treeBoundBox::TOP)
01006     {
01007         facePoint.y() = bb.max().y();
01008         faceID = treeBoundBox::TOPBIT;
01009     }
01010     else if (keepFaceID == treeBoundBox::BACK)
01011     {
01012         facePoint.z() = bb.min().z();
01013         faceID = treeBoundBox::BACKBIT;
01014     }
01015     else if (keepFaceID == treeBoundBox::FRONT)
01016     {
01017         facePoint.z() = bb.max().z();
01018         faceID = treeBoundBox::FRONTBIT;
01019     }
01020 
01021 
01022     if (debug)
01023     {
01024         if (faceID != bb.faceBits(facePoint))
01025         {
01026             FatalErrorIn("indexedOctree<Type>::pushPointIntoFace(..)")
01027                 << "Pushed point from " << pt
01028                 << " on face:" << ptFaceID << " of bb:" << bb << endl
01029                 << "onto " << facePoint
01030                 << " on face:" << faceID
01031                 << " which is not consistent with geometric face "
01032                 << bb.faceBits(facePoint)
01033                 << abort(FatalError);
01034         }
01035         if (bb.posBits(facePoint) != 0)
01036         {
01037             FatalErrorIn("indexedOctree<Type>::pushPointIntoFace(..)")
01038                 << " bb:" << bb << endl
01039                 << "does not contain perturbed point "
01040                 << facePoint << abort(FatalError);
01041         }
01042     }
01043 
01044 
01045     return facePoint;
01046 }
01047 
01048 
01051 //template <class Type>
01052 //void Foam::indexedOctree<Type>::checkMultipleFaces
01053 //(
01054 //    const treeBoundBox& bb,
01055 //    const vector& dir,          // end-start
01056 //    pointIndexHit& faceHitInfo,
01057 //    direction& faceID
01058 //)
01059 //{
01060 //    // Do the quick elimination of no or one face.
01061 //    if
01062 //    (
01063 //        (faceID == 0)
01064 //     || (faceID == treeBoundBox::LEFTBIT)
01065 //     || (faceID == treeBoundBox::RIGHTBIT)
01066 //     || (faceID == treeBoundBox::BOTTOMBIT)
01067 //     || (faceID == treeBoundBox::TOPBIT)
01068 //     || (faceID == treeBoundBox::BACKBIT)
01069 //     || (faceID == treeBoundBox::FRONTBIT)
01070 //    )
01071 //    {
01072 //        return;
01073 //    }
01074 //
01075 //
01076 //    // Check the direction of vector w.r.t. faces being intersected.
01077 //    FixedList<scalar, 6> inproducts(-GREAT);
01078 //
01079 //    direction nFaces = 0;
01080 //
01081 //    if (faceID & treeBoundBox::LEFTBIT)
01082 //    {
01083 //        inproducts[treeBoundBox::LEFT] = mag
01084 //        (
01085 //            treeBoundBox::faceNormals[treeBoundBox::LEFT]
01086 //          & dir
01087 //        );
01088 //        nFaces++;
01089 //    }
01090 //    if (faceID & treeBoundBox::RIGHTBIT)
01091 //    {
01092 //        inproducts[treeBoundBox::RIGHT] = mag
01093 //        (
01094 //            treeBoundBox::faceNormals[treeBoundBox::RIGHT]
01095 //          & dir
01096 //        );
01097 //        nFaces++;
01098 //    }
01099 //
01100 //    if (faceID & treeBoundBox::BOTTOMBIT)
01101 //    {
01102 //        inproducts[treeBoundBox::BOTTOM] = mag
01103 //        (
01104 //            treeBoundBox::faceNormals[treeBoundBox::BOTTOM]
01105 //          & dir
01106 //        );
01107 //        nFaces++;
01108 //    }
01109 //    if (faceID & treeBoundBox::TOPBIT)
01110 //    {
01111 //        inproducts[treeBoundBox::TOP] = mag
01112 //        (
01113 //            treeBoundBox::faceNormals[treeBoundBox::TOP]
01114 //          & dir
01115 //        );
01116 //        nFaces++;
01117 //    }
01118 //
01119 //    if (faceID & treeBoundBox::BACKBIT)
01120 //    {
01121 //        inproducts[treeBoundBox::BACK] = mag
01122 //        (
01123 //            treeBoundBox::faceNormals[treeBoundBox::BACK]
01124 //          & dir
01125 //        );
01126 //        nFaces++;
01127 //    }
01128 //    if (faceID & treeBoundBox::FRONTBIT)
01129 //    {
01130 //        inproducts[treeBoundBox::FRONT] = mag
01131 //        (
01132 //            treeBoundBox::faceNormals[treeBoundBox::FRONT]
01133 //          & dir
01134 //        );
01135 //        nFaces++;
01136 //    }
01137 //
01138 //    if (nFaces == 0 || nFaces == 1 || nFaces > 3)
01139 //    {
01140 //        FatalErrorIn("indexedOctree<Type>::checkMultipleFaces(..)")
01141 //            << "Problem : nFaces:" << nFaces << abort(FatalError);
01142 //    }
01143 //
01144 //    // Keep point on most perpendicular face; shift it away from the aligned
01145 //    // ones.
01146 //    // E.g. line hits top and left face:
01147 //    //     a
01148 //    // ----+----+
01149 //    //     |    |
01150 //    //     |    |
01151 //    //     +----+
01152 //    // Shift point down (away from top):
01153 //    //
01154 //    //    a+----+
01155 //    // ----|    |
01156 //    //     |    |
01157 //    //     +----+
01158 //
01159 //    label maxIndex = -1;
01160 //    scalar maxInproduct = -GREAT;
01161 //
01162 //    for (direction i = 0; i < 6; i++)
01163 //    {
01164 //        if (inproducts[i] > maxInproduct)
01165 //        {
01166 //            maxInproduct = inproducts[i];
01167 //            maxIndex = i;
01168 //        }
01169 //    }
01170 //
01171 //    if (maxIndex == -1)
01172 //    {
01173 //        FatalErrorIn("indexedOctree<Type>::checkMultipleFaces(..)")
01174 //            << "Problem maxIndex:" << maxIndex << " inproducts:" << inproducts
01175 //            << abort(FatalError);
01176 //    }
01177 //
01178 //    const point oldPoint(faceHitInfo.rawPoint());
01179 //    const direction oldFaceID = faceID;
01180 //
01181 //    // 1. Push point into bb, away from all corners
01182 //
01183 //    faceHitInfo.rawPoint() = pushPoint(bb, oldFaceID, oldPoint, true);
01184 //
01185 //    // 2. Snap it back onto the preferred face
01186 //
01187 //    if (maxIndex == treeBoundBox::LEFT)
01188 //    {
01189 //        faceHitInfo.rawPoint().x() = bb.min().x();
01190 //        faceID = treeBoundBox::LEFTBIT;
01191 //    }
01192 //    else if (maxIndex == treeBoundBox::RIGHT)
01193 //    {
01194 //        faceHitInfo.rawPoint().x() = bb.max().x();
01195 //        faceID = treeBoundBox::RIGHTBIT;
01196 //    }
01197 //    else if (maxIndex == treeBoundBox::BOTTOM)
01198 //    {
01199 //        faceHitInfo.rawPoint().y() = bb.min().y();
01200 //        faceID = treeBoundBox::BOTTOMBIT;
01201 //    }
01202 //    else if (maxIndex == treeBoundBox::TOP)
01203 //    {
01204 //        faceHitInfo.rawPoint().y() = bb.max().y();
01205 //        faceID = treeBoundBox::TOPBIT;
01206 //    }
01207 //    else if (maxIndex == treeBoundBox::BACK)
01208 //    {
01209 //        faceHitInfo.rawPoint().z() = bb.min().z();
01210 //        faceID = treeBoundBox::BACKBIT;
01211 //    }
01212 //    else if (maxIndex == treeBoundBox::FRONT)
01213 //    {
01214 //        faceHitInfo.rawPoint().z() = bb.max().z();
01215 //        faceID = treeBoundBox::FRONTBIT;
01216 //    }
01217 //
01218 //    Pout<< "From ray:" << dir
01219 //        << " from point:" << oldPoint
01220 //        << " on faces:" << faceString(oldFaceID)
01221 //        << " of bb:" << bb
01222 //        << " with inprods:" << inproducts
01223 //        << " maxIndex:" << maxIndex << endl
01224 //        << "perturbed to point:" << faceHitInfo.rawPoint()
01225 //        << " on face:" << faceString(faceID)
01226 //        << endl;
01227 //
01228 //
01229 //    if (debug)
01230 //    {
01231 //        if (faceID != bb.faceBits(faceHitInfo.rawPoint()))
01232 //        {
01233 //            FatalErrorIn("indexedOctree<Type>::checkMultipleFaces(..)")
01234 //                << "Pushed point from " << oldPoint
01235 //                << " on face:" << oldFaceID << " of bb:" << bb << endl
01236 //                << "onto " << faceHitInfo.rawPoint()
01237 //                << " on face:" << faceID
01238 //                << " which is not consistent with geometric face "
01239 //                <<  bb.faceBits(faceHitInfo.rawPoint())
01240 //                << abort(FatalError);
01241 //        }
01242 //    }
01243 //}
01244 
01245 
01246 // Get parent node and octant. Return false if top of tree reached.
01247 template <class Type>
01248 bool Foam::indexedOctree<Type>::walkToParent
01249 (
01250     const label nodeI,
01251     const direction octant,
01252 
01253     label& parentNodeI,
01254     label& parentOctant
01255 ) const
01256 {
01257     parentNodeI = nodes_[nodeI].parent_;
01258 
01259     if (parentNodeI == -1)
01260     {
01261         // Reached edge of tree
01262         return false;
01263     }
01264 
01265     const node& parentNode = nodes_[parentNodeI];
01266 
01267     // Find octant nodeI is in.
01268     parentOctant = 255;
01269 
01270     for (direction i = 0; i < parentNode.subNodes_.size(); i++)
01271     {
01272         labelBits index = parentNode.subNodes_[i];
01273 
01274         if (isNode(index) && getNode(index) == nodeI)
01275         {
01276             parentOctant = i;
01277             break;
01278         }
01279     }
01280 
01281     if (parentOctant == 255)
01282     {
01283         FatalErrorIn("walkToParent(..)")
01284             << "Problem: no parent found for octant:" << octant
01285             << " node:" << nodeI
01286             << abort(FatalError);
01287     }
01288 
01289     return true;
01290 }
01291 
01292 
01293 // Walk tree to neighbouring node. Gets current position as
01294 // node and octant in this node and walks in the direction given by
01295 // the facePointBits (combination of treeBoundBox::LEFTBIT, TOPBIT etc.)
01296 // Returns false if edge of tree hit.
01297 template <class Type>
01298 bool Foam::indexedOctree<Type>::walkToNeighbour
01299 (
01300     const point& facePoint,
01301     const direction faceID,  // face(s) that facePoint is on
01302     label& nodeI,
01303     direction& octant
01304 ) const
01305 {
01306     label oldNodeI = nodeI;
01307     direction oldOctant = octant;
01308 
01309     // Find out how to test for octant. Say if we want to go left we need
01310     // to traverse up the tree until we hit a node where our octant is
01311     // on the right.
01312 
01313     // Coordinate direction to test
01314     const direction X = treeBoundBox::RIGHTHALF;
01315     const direction Y = treeBoundBox::TOPHALF;
01316     const direction Z = treeBoundBox::FRONTHALF;
01317 
01318     direction octantMask = 0;
01319     direction wantedValue = 0;
01320 
01321     if ((faceID & treeBoundBox::LEFTBIT) != 0)
01322     {
01323         // We want to go left so check if in right octant (i.e. x-bit is set)
01324         octantMask |= X;
01325         wantedValue |= X;
01326     }
01327     else if ((faceID & treeBoundBox::RIGHTBIT) != 0)
01328     {
01329         octantMask |= X;  // wantedValue already 0
01330     }
01331 
01332     if ((faceID & treeBoundBox::BOTTOMBIT) != 0)
01333     {
01334         // Want to go down so check for y-bit set.
01335         octantMask |= Y;
01336         wantedValue |= Y;
01337     }
01338     else if ((faceID & treeBoundBox::TOPBIT) != 0)
01339     {
01340         // Want to go up so check for y-bit not set.
01341         octantMask |= Y;
01342     }
01343 
01344     if ((faceID & treeBoundBox::BACKBIT) != 0)
01345     {
01346         octantMask |= Z;
01347         wantedValue |= Z;
01348     }
01349     else if ((faceID & treeBoundBox::FRONTBIT) != 0)
01350     {
01351         octantMask |= Z;
01352     }
01353 
01354     // So now we have the coordinate directions in the octant we need to check
01355     // and the resulting value.
01356 
01357     /*
01358     // +---+---+
01359     // |   |   |
01360     // |   |   |
01361     // |   |   |
01362     // +---+-+-+
01363     // |   | | |
01364     // |  a+-+-+
01365     // |   |\| |
01366     // +---+-+-+
01367     //        \
01368     //
01369     // e.g. ray is at (a) in octant 0(or 4) with faceIDs : LEFTBIT+TOPBIT.
01370     // If we would be in octant 1(or 5) we could go to the correct octant
01371     // in the same node by just flipping the x and y bits (exoring).
01372     // But if we are not in octant 1/5 we have to go up until we are.
01373     // In general for leftbit+topbit:
01374     // - we have to check for x and y : octantMask  = 011
01375     // - the checked bits have to be  : wantedValue = ?01
01376     */
01377 
01378     //Pout<< "For point " << facePoint << endl;
01379 
01380     // Go up until we have chance to cross to the wanted direction
01381     while (wantedValue != (octant & octantMask))
01382     {
01383         // Go up to the parent.
01384 
01385         // Remove the directions that are not on the boundary of the parent.
01386         // See diagram above
01387         if (wantedValue & X)    // && octantMask&X
01388         {
01389             // Looking for right octant.
01390             if (octant & X)
01391             {
01392                 // My octant is one of the left ones so punch out the x check
01393                 octantMask &= ~X;
01394                 wantedValue &= ~X;
01395             }
01396         }
01397         else
01398         {
01399             if (!(octant & X))
01400             {
01401                 // My octant is right but I want to go left.
01402                 octantMask &= ~X;
01403                 wantedValue &= ~X;
01404             }
01405         }
01406 
01407         if (wantedValue & Y)
01408         {
01409             if (octant & Y)
01410             {
01411                 octantMask &= ~Y;
01412                 wantedValue &= ~Y;
01413             }
01414         }
01415         else
01416         {
01417             if (!(octant & Y))
01418             {
01419                 octantMask &= ~Y;
01420                 wantedValue &= ~Y;
01421             }
01422         }
01423 
01424         if (wantedValue & Z)
01425         {
01426             if (octant & Z)
01427             {
01428                 octantMask &= ~Z;
01429                 wantedValue &= ~Z;
01430             }
01431         }
01432         else
01433         {
01434             if (!(octant & Z))
01435             {
01436                 octantMask &= ~Z;
01437                 wantedValue &= ~Z;
01438             }
01439         }
01440 
01441 
01442         label parentNodeI;
01443         label parentOctant;
01444         walkToParent(nodeI, octant, parentNodeI, parentOctant);
01445 
01446         if (parentNodeI == -1)
01447         {
01448             // Reached edge of tree
01449             return false;
01450         }
01451 
01452         //Pout<< "    walked from node:" << nodeI << " octant:" << octant
01453         //    << " bb:" << nodes_[nodeI].bb_.subBbox(octant) << endl
01454         //    << "    to:" << parentNodeI << " octant:" << parentOctant
01455         //    << " bb:" << nodes_[parentNodeI].bb_.subBbox(parentOctant)
01456         //    << endl;
01457         //
01458         //Pout<< "    octantMask:" << octantMask
01459         //    << " wantedValue:" << wantedValue << endl;
01460 
01461         nodeI = parentNodeI;
01462         octant = parentOctant;
01463     }
01464 
01465     // So now we hit the correct parent node. Switch to the correct octant.
01466     // Is done by jumping to the 'other half' so if e.g. in x direction in
01467     // right half we now jump to the left half.
01468     octant ^= octantMask;
01469 
01470     //Pout<< "    to node:" << nodeI << " octant:" << octant
01471     //    << " subBb:" <<subBbox(nodeI, octant) << endl;
01472 
01473 
01474     if (debug)
01475     {
01476         const treeBoundBox subBb(subBbox(nodeI, octant));
01477 
01478         if (!subBb.contains(facePoint))
01479         {
01480             FatalErrorIn("indexedOctree<Type>::walkToNeighbour(..)")
01481                 << "When searching for " << facePoint
01482                 << " ended up in node:" << nodeI
01483                 << " octant:" << octant
01484                 << " with bb:" << subBb
01485                 << abort(FatalError);
01486         }
01487     }
01488 
01489 
01490     // See if we need to travel down. Note that we already go into the
01491     // the first level ourselves (instead of having findNode decide)
01492     labelBits index = nodes_[nodeI].subNodes_[octant];
01493 
01494     if (isNode(index))
01495     {
01496         labelBits node = findNode(getNode(index), facePoint);
01497 
01498         nodeI = getNode(node);
01499         octant = getOctant(node);
01500     }
01501 
01502 
01503     if (debug)
01504     {
01505         const treeBoundBox subBb(subBbox(nodeI, octant));
01506 
01507         if (nodeI == oldNodeI && octant == oldOctant)
01508         {
01509             FatalErrorIn("indexedOctree<Type>::walkToNeighbour(..)")
01510                 << "Did not go to neighbour when searching for " << facePoint
01511                 << endl
01512                 << "    starting from face:" << faceString(faceID)
01513                 << " node:" << nodeI
01514                 << " octant:" << octant
01515                 << " bb:" << subBb
01516                 << abort(FatalError);
01517         }
01518 
01519         if (!subBb.contains(facePoint))
01520         {
01521             FatalErrorIn("indexedOctree<Type>::walkToNeighbour(..)")
01522                 << "When searching for " << facePoint
01523                 << " ended up in node:" << nodeI
01524                 << " octant:" << octant
01525                 << " bb:" << subBb
01526                 << abort(FatalError);
01527         }
01528     }
01529 
01530 
01531     return true;
01532 }
01533 
01534 
01535 template <class Type>
01536 Foam::word Foam::indexedOctree<Type>::faceString
01537 (
01538     const direction faceID
01539 )
01540 {
01541     word desc;
01542 
01543     if (faceID == 0)
01544     {
01545         desc = "noFace";
01546     }
01547     if (faceID & treeBoundBox::LEFTBIT)
01548     {
01549         if (!desc.empty()) desc += "+";
01550         desc += "left";
01551     }
01552     if (faceID & treeBoundBox::RIGHTBIT)
01553     {
01554         if (!desc.empty()) desc += "+";
01555         desc += "right";
01556     }
01557     if (faceID & treeBoundBox::BOTTOMBIT)
01558     {
01559         if (!desc.empty()) desc += "+";
01560         desc += "bottom";
01561     }
01562     if (faceID & treeBoundBox::TOPBIT)
01563     {
01564         if (!desc.empty()) desc += "+";
01565         desc += "top";
01566     }
01567     if (faceID & treeBoundBox::BACKBIT)
01568     {
01569         if (!desc.empty()) desc += "+";
01570         desc += "back";
01571     }
01572     if (faceID & treeBoundBox::FRONTBIT)
01573     {
01574         if (!desc.empty()) desc += "+";
01575         desc += "front";
01576     }
01577     return desc;
01578 }
01579 
01580 
01581 // Traverse a node. If intersects a triangle return first intersection point:
01582 //  hitInfo.index = index of shape
01583 //  hitInfo.point = point on shape
01584 // Else return a miss and the bounding box face hit:
01585 //  hitInfo.point = coordinate of intersection of ray with bounding box
01586 //  hitBits  = posbits of point on bounding box
01587 template <class Type>
01588 void Foam::indexedOctree<Type>::traverseNode
01589 (
01590     const bool findAny,
01591     const point& treeStart,
01592     const vector& treeVec,
01593 
01594     const point& start,
01595     const point& end,
01596     const label nodeI,
01597     const direction octant,
01598 
01599     pointIndexHit& hitInfo,
01600     direction& hitBits
01601 ) const
01602 {
01603     if (debug)
01604     {
01605         const treeBoundBox octantBb(subBbox(nodeI, octant));
01606 
01607         if (octantBb.posBits(start) != 0)
01608         {
01609             FatalErrorIn("indexedOctree<Type>::traverseNode(..)")
01610                 << "Node:" << nodeI << " octant:" << octant
01611                 << " bb:" << octantBb << endl
01612                 << "does not contain point " << start << abort(FatalError);
01613         }
01614     }
01615 
01616 
01617     const node& nod = nodes_[nodeI];
01618 
01619     labelBits index = nod.subNodes_[octant];
01620 
01621     if (isContent(index))
01622     {
01623         const labelList& indices = contents_[getContent(index)];
01624 
01625         if (findAny)
01626         {
01627             // Find any intersection
01628 
01629             forAll(indices, elemI)
01630             {
01631                 label shapeI = indices[elemI];
01632 
01633                 point pt;
01634                 bool hit = shapes_.intersects(shapeI, start, end, pt);
01635 
01636                 if (hit)
01637                 {
01638                     // Hit so pt is nearer than nearestPoint.
01639                     // Update hit info
01640                     hitInfo.setHit();
01641                     hitInfo.setIndex(shapeI);
01642                     hitInfo.setPoint(pt);
01643                     return;
01644                 }
01645             }
01646         }
01647         else
01648         {
01649             // Find nearest intersection.
01650 
01651             point nearestPoint(end);
01652 
01653             forAll(indices, elemI)
01654             {
01655                 label shapeI = indices[elemI];
01656 
01657                 point pt;
01658                 bool hit = shapes_.intersects(shapeI, start, nearestPoint, pt);
01659 
01660                 if (hit)
01661                 {
01662                     // Hit so pt is nearer than nearestPoint.
01663                     nearestPoint = pt;
01664                     // Update hit info
01665                     hitInfo.setHit();
01666                     hitInfo.setIndex(shapeI);
01667                     hitInfo.setPoint(pt);
01668                 }
01669             }
01670 
01671             if (hitInfo.hit())
01672             {
01673                 // Found intersected shape.
01674                 return;
01675             }
01676         }
01677     }
01678 
01679     // Nothing intersected in this node
01680     // Traverse to end of node. Do by ray tracing back from end and finding
01681     // intersection with bounding box of node.
01682     // start point is now considered to be inside bounding box.
01683 
01684     const treeBoundBox octantBb(subBbox(nodeI, octant));
01685 
01686     point pt;
01687     bool intersected = octantBb.intersects
01688     (
01689         end,            //treeStart,
01690         (start-end),    //treeVec,
01691 
01692         end,
01693         start,
01694 
01695         pt,
01696         hitBits
01697     );
01698 
01699 
01700     if (intersected)
01701     {
01702         // Return miss. Set misspoint to face.
01703         hitInfo.setPoint(pt);
01704     }
01705     else
01706     {
01707         // Rare case: did not find intersection of ray with octantBb. Can
01708         // happen if end is on face/edge of octantBb. Do like
01709         // lagrangian and re-track with perturbed vector (slightly pushed out
01710         // of bounding box)
01711 
01712         point perturbedEnd(pushPoint(octantBb, end, false));
01713 
01714 
01715         //if (debug)
01716         {
01717             // Dump octantBb to obj
01718             writeOBJ(nodeI, octant);
01719             // Dump ray to obj as well
01720             {
01721                 OFstream str("ray.obj");
01722                 meshTools::writeOBJ(str, start);
01723                 meshTools::writeOBJ(str, end);
01724                 str << "l 1 2" << nl;
01725             }
01726             WarningIn("indexedOctree<Type>::traverseNode(..)")
01727                 << "Did not intersect ray from endpoint:" << end
01728                 << " to startpoint:" << start
01729                 << " with bounding box:" << octantBb << nl
01730                 << "Re-intersecting with perturbed endpoint:" << perturbedEnd
01731                 << endl;
01732         }
01733 
01734         traverseNode
01735         (
01736             findAny,
01737             treeStart,
01738             treeVec,
01739             start,
01740             perturbedEnd,
01741             nodeI,
01742             octant,
01743 
01744             hitInfo,
01745             hitBits
01746         );
01747     }
01748 }
01749 
01750 
01751 // Find first intersection
01752 template <class Type>
01753 Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
01754 (
01755     const bool findAny,
01756     const point& treeStart,
01757     const point& treeEnd,
01758     const label startNodeI,
01759     const direction startOctant,
01760     const bool verbose
01761 ) const
01762 {
01763     const vector treeVec(treeEnd - treeStart);
01764 
01765     // Current node as parent+octant
01766     label nodeI = startNodeI;
01767     direction octant = startOctant;
01768 
01769     if (verbose)
01770     {
01771         Pout<< "findLine : treeStart:" << treeStart
01772             << " treeEnd:" << treeEnd << endl
01773             << "node:" << nodeI
01774             << " octant:" << octant
01775             << " bb:" << subBbox(nodeI, octant) << endl;
01776     }
01777 
01778     // Current position. Initialize to miss
01779     pointIndexHit hitInfo(false, treeStart, -1);
01780 
01781     //while (true)
01782     label i = 0;
01783     for (; i < 100000; i++)
01784     {
01785         // Ray-trace to end of current node. Updates point (either on triangle
01786         // in case of hit or on node bounding box in case of miss)
01787 
01788         const treeBoundBox octantBb(subBbox(nodeI, octant));
01789 
01790         // Make sure point is away from any edges/corners
01791         point startPoint
01792         (
01793             pushPointIntoFace
01794             (
01795                 octantBb,
01796                 treeVec,
01797                 hitInfo.rawPoint()
01798             )
01799         );
01800 
01801         if (verbose)
01802         {
01803             Pout<< "iter:" << i
01804                 << " at current:" << hitInfo.rawPoint()
01805                 << " (perturbed:" << startPoint << ")" << endl
01806                 << "    node:" << nodeI
01807                 << " octant:" << octant
01808                 << " bb:" << subBbox(nodeI, octant) << endl;
01809         }
01810 
01811 
01812 
01813 
01814         // Faces of current bounding box current point is on
01815         direction hitFaceID = 0;
01816 
01817         traverseNode
01818         (
01819             findAny,
01820             treeStart,
01821             treeVec,
01822 
01823             startPoint,     // Note: pass in copy since hitInfo
01824                             // also used in return value.
01825             treeEnd,        // pass in overall end so is nicely outside bb
01826             nodeI,
01827             octant,
01828 
01829             hitInfo,
01830             hitFaceID
01831         );
01832 
01833         // Did we hit a triangle?
01834         if (hitInfo.hit())
01835         {
01836             break;
01837         }
01838 
01839         if (hitFaceID == 0 || hitInfo.rawPoint() == treeEnd)
01840         {
01841             // endpoint inside the tree. Return miss.
01842             break;
01843         }
01844 
01845         // Create a point on other side of face.
01846         point perturbedPoint
01847         (
01848             pushPoint
01849             (
01850                 octantBb,
01851                 hitFaceID,
01852                 hitInfo.rawPoint(),
01853                 false                   // push outside of octantBb
01854             )
01855         );
01856 
01857         if (verbose)
01858         {
01859             Pout<< "    iter:" << i
01860                 << " hit face:" << faceString(hitFaceID)
01861                 << " at:" << hitInfo.rawPoint() << nl
01862                 << "    node:" << nodeI
01863                 << " octant:" << octant
01864                 << " bb:" << subBbox(nodeI, octant) << nl
01865                 << "    walking to neighbour containing:" << perturbedPoint
01866                 << endl;
01867         }
01868 
01869 
01870         // Nothing hit so we are on face of bounding box (given as node+octant+
01871         // position bits). Traverse to neighbouring node. Use slightly perturbed
01872         // point.
01873 
01874         bool ok = walkToNeighbour
01875         (
01876             perturbedPoint,
01877             hitFaceID,  // face(s) that hitInfo is on
01878 
01879             nodeI,
01880             octant
01881         );
01882 
01883         if (!ok)
01884         {
01885             // Hit the edge of the tree. Return miss.
01886             break;
01887         }
01888 
01889         if (verbose)
01890         {
01891             const treeBoundBox octantBb(subBbox(nodeI, octant));
01892             Pout<< "    walked for point:" << hitInfo.rawPoint() << endl
01893                 << "    to neighbour node:" << nodeI
01894                 << " octant:" << octant
01895                 << " face:" << faceString(octantBb.faceBits(hitInfo.rawPoint()))
01896                 << " of octantBb:" << octantBb << endl
01897                 << endl;
01898         }
01899     }
01900 
01901     if (i == 100000)
01902     {
01903         // Probably in loop.
01904         if (!verbose)
01905         {
01906             // Redo intersection but now with verbose flag switched on.
01907             return findLine
01908             (
01909                 findAny,
01910                 treeStart,
01911                 treeEnd,
01912                 startNodeI,
01913                 startOctant,
01914                 true            //verbose
01915             );
01916         }
01917         if (debug)
01918         {
01919             FatalErrorIn("indexedOctree<Type>::findLine(..)")
01920                 << "Got stuck in loop raytracing from:" << treeStart
01921                 << " to:" << treeEnd << endl
01922                 << "inside top box:" << subBbox(startNodeI, startOctant)
01923                 << abort(FatalError);
01924         }
01925         else
01926         {
01927             WarningIn("indexedOctree<Type>::findLine(..)")
01928                 << "Got stuck in loop raytracing from:" << treeStart
01929                 << " to:" << treeEnd << endl
01930                 << "inside top box:" << subBbox(startNodeI, startOctant)
01931                 << endl;
01932         }
01933     }
01934 
01935     return hitInfo;
01936 }
01937 
01938 
01939 // Find first intersection
01940 template <class Type>
01941 Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
01942 (
01943     const bool findAny,
01944     const point& start,
01945     const point& end
01946 ) const
01947 {
01948     pointIndexHit hitInfo;
01949 
01950     if (nodes_.size())
01951     {
01952         const treeBoundBox& treeBb = nodes_[0].bb_;
01953 
01954         // No effort is made to deal with points which are on edge of tree
01955         // bounding box for now.
01956 
01957         direction startBit = treeBb.posBits(start);
01958         direction endBit = treeBb.posBits(end);
01959 
01960         if ((startBit & endBit) != 0)
01961         {
01962             // Both start and end outside domain and in same block.
01963             return pointIndexHit(false, vector::zero, -1);
01964         }
01965 
01966 
01967         // Trim segment to treeBb
01968 
01969         point trackStart(start);
01970         point trackEnd(end);
01971 
01972         if (startBit != 0)
01973         {
01974             // Track start to inside domain.
01975             if (!treeBb.intersects(start, end, trackStart))
01976             {
01977                 return pointIndexHit(false, vector::zero, -1);
01978             }
01979         }
01980 
01981         if (endBit != 0)
01982         {
01983             // Track end to inside domain.
01984             if (!treeBb.intersects(end, trackStart, trackEnd))
01985             {
01986                 return pointIndexHit(false, vector::zero, -1);
01987             }
01988         }
01989 
01990 
01991         // Find lowest level tree node that start is in.
01992         labelBits index = findNode(0, trackStart);
01993 
01994         label parentNodeI = getNode(index);
01995         direction octant = getOctant(index);
01996 
01997         hitInfo = findLine
01998         (
01999             findAny,
02000             trackStart,
02001             trackEnd,
02002             parentNodeI,
02003             octant
02004         );
02005     }
02006 
02007     return hitInfo;
02008 }
02009 
02010 
02011 template <class Type>
02012 void Foam::indexedOctree<Type>::findBox
02013 (
02014     const label nodeI,
02015     const treeBoundBox& searchBox,
02016     labelHashSet& elements
02017 ) const
02018 {
02019     const node& nod = nodes_[nodeI];
02020     const treeBoundBox& nodeBb = nod.bb_;
02021 
02022     for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
02023     {
02024         labelBits index = nod.subNodes_[octant];
02025 
02026         if (isNode(index))
02027         {
02028             const treeBoundBox& subBb = nodes_[getNode(index)].bb_;
02029 
02030             if (subBb.overlaps(searchBox))
02031             {
02032                 findBox(getNode(index), searchBox, elements);
02033             }
02034         }
02035         else if (isContent(index))
02036         {
02037             const treeBoundBox subBb(nodeBb.subBbox(octant));
02038 
02039             if (subBb.overlaps(searchBox))
02040             {
02041                 const labelList& indices = contents_[getContent(index)];
02042 
02043                 forAll(indices, i)
02044                 {
02045                     label shapeI = indices[i];
02046 
02047                     if (shapes_.overlaps(shapeI, searchBox))
02048                     {
02049                         elements.insert(shapeI);
02050                     }
02051                 }
02052             }
02053         }
02054     }
02055 }
02056 
02057 
02058 // Number of elements in node.
02059 template <class Type>
02060 Foam::label Foam::indexedOctree<Type>::countElements
02061 (
02062     const labelBits index
02063 ) const
02064 {
02065     label nElems = 0;
02066 
02067     if (isNode(index))
02068     {
02069         // tree node.
02070         label nodeI = getNode(index);
02071 
02072         const node& nod = nodes_[nodeI];
02073 
02074         for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
02075         {
02076             nElems += countElements(nod.subNodes_[octant]);
02077         }
02078     }
02079     else if (isContent(index))
02080     {
02081         nElems += contents_[getContent(index)].size();
02082     }
02083     else
02084     {
02085         // empty node
02086     }
02087 
02088     return nElems;
02089 }
02090 
02091 
02092 template <class Type>
02093 void Foam::indexedOctree<Type>::writeOBJ
02094 (
02095     const label nodeI,
02096     const direction octant
02097 ) const
02098 {
02099     OFstream str
02100     (
02101         "node" + Foam::name(nodeI) + "_octant" + Foam::name(octant) + ".obj"
02102     );
02103 
02104     labelBits index = nodes_[nodeI].subNodes_[octant];
02105 
02106     treeBoundBox subBb;
02107 
02108     if (isNode(index))
02109     {
02110         subBb = nodes_[getNode(index)].bb_;
02111     }
02112     else if (isContent(index))
02113     {
02114         subBb = nodes_[nodeI].bb_.subBbox(octant);
02115     }
02116 
02117     Pout<< "dumpContentNode : writing node:" << nodeI << " octant:" << octant
02118         << " to " << str.name() << endl;
02119 
02120     label vertI = 0;
02121 
02122     // Dump bounding box
02123     pointField bbPoints(subBb.points());
02124 
02125     label pointVertI = vertI;
02126     forAll(bbPoints, i)
02127     {
02128         meshTools::writeOBJ(str, bbPoints[i]);
02129         vertI++;
02130     }
02131 
02132     forAll(treeBoundBox::edges, i)
02133     {
02134         const edge& e = treeBoundBox::edges[i];
02135 
02136         str<< "l " << e[0]+pointVertI+1 << ' ' << e[1]+pointVertI+1 << nl;
02137     }
02138 }
02139 
02140 
02141 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
02142 
02143 template <class Type>
02144 Foam::indexedOctree<Type>::indexedOctree(const Type& shapes)
02145 :
02146     shapes_(shapes),
02147     nodes_(0),
02148     contents_(0),
02149     nodeTypes_(0)
02150 {}
02151 
02152 
02153 template <class Type>
02154 Foam::indexedOctree<Type>::indexedOctree
02155 (
02156     const Type& shapes,
02157     const List<node>& nodes,
02158     const labelListList& contents
02159 )
02160 :
02161     shapes_(shapes),
02162     nodes_(nodes),
02163     contents_(contents),
02164     nodeTypes_(0)
02165 {}
02166 
02167 
02168 template <class Type>
02169 Foam::indexedOctree<Type>::indexedOctree
02170 (
02171     const Type& shapes,
02172     const treeBoundBox& bb,
02173     const label maxLevels,          // maximum number of levels
02174     const scalar maxLeafRatio,
02175     const scalar maxDuplicity
02176 )
02177 :
02178     shapes_(shapes),
02179     nodes_(0),
02180     contents_(0),
02181     nodeTypes_(0)
02182 {
02183     if (debug)
02184     {
02185         Pout<< "indexedOctree<Type>::indexedOctree:" << nl
02186             << "    shapes:" << shapes.size() << nl
02187             << "    bb:" << bb << nl
02188             << endl;
02189     }
02190 
02191     if (shapes.size() == 0)
02192     {
02193         return;
02194     }
02195 
02196     // Start off with one node with all shapes in it.
02197     DynamicList<node> nodes(label(shapes.size() / maxLeafRatio));
02198     DynamicList<labelList> contents(label(shapes.size() / maxLeafRatio));
02199     contents.append(identity(shapes.size()));
02200 
02201     // Create topnode.
02202     node topNode(divide(bb, contents, 0));
02203     nodes.append(topNode);
02204 
02205 
02206     // Now have all contents at level 1. Create levels by splitting levels
02207     // above.
02208 
02209     label nLevels = 1;
02210 
02211     for (; nLevels < maxLevels; nLevels++)
02212     {
02213         // Count number of references into shapes (i.e. contents)
02214         label nEntries = 0;
02215         forAll(contents, i)
02216         {
02217             nEntries += contents[i].size();
02218         }
02219 
02220         if (debug)
02221         {
02222             Pout<< "indexedOctree<Type>::indexedOctree:" << nl
02223                 << "    nLevels:" << nLevels << nl
02224                 << "    nEntries per treeLeaf:" << nEntries/contents.size()
02225                 << nl
02226                 << "    nEntries per shape (duplicity):"
02227                 << nEntries/shapes.size()
02228                 << nl
02229                 << endl;
02230         }
02231 
02232         if
02233         (
02234             //nEntries < maxLeafRatio*contents.size()
02235          // ||
02236             nEntries > maxDuplicity*shapes.size()
02237         )
02238         {
02239             break;
02240         }
02241 
02242 
02243         // Split nodes with more than maxLeafRatio elements
02244         label nOldNodes = nodes.size();
02245         splitNodes
02246         (
02247             label(maxLeafRatio),
02248             nodes,
02249             contents
02250         );
02251 
02252         if (nOldNodes == nodes.size())
02253         {
02254             break;
02255         }
02256     }
02257 
02258     // Shrink
02259     nodes.shrink();
02260     contents.shrink();
02261 
02262 
02263     // Compact such that deeper level contents are always after the
02264     // ones for a shallower level. This way we can slice a coarser level
02265     // off the tree.
02266     contents_.setSize(contents.size());
02267     label compactI = 0;
02268 
02269     label level = 0;
02270 
02271     while (true)
02272     {
02273         compactContents
02274         (
02275             nodes,
02276             contents,
02277             level,
02278             0,
02279             0,
02280             contents_,
02281             compactI
02282         );
02283 
02284         if (compactI == contents_.size())
02285         {
02286             // Transferred all contents to contents_ (in order breadth first)
02287             break;
02288         }
02289 
02290         level++;
02291     }
02292     nodes_.transfer(nodes);
02293     nodes.clear();
02294 
02295     if (debug)
02296     {
02297         label nEntries = 0;
02298         forAll(contents_, i)
02299         {
02300             nEntries += contents_[i].size();
02301         }
02302 
02303         Pout<< "indexedOctree<Type>::indexedOctree"
02304             << " : finished construction of tree of:" << shapes.typeName
02305             << nl
02306             << "    bb:" << this->bb() << nl
02307             << "    shapes:" << shapes.size() << nl
02308             << "    nLevels:" << nLevels << nl
02309             << "    treeNodes:" << nodes_.size() << nl
02310             << "    nEntries:" << nEntries << nl
02311             << "        per treeLeaf:"
02312             << scalar(nEntries)/contents.size() << nl
02313             << "        per shape (duplicity):"
02314             << scalar(nEntries)/shapes.size() << nl
02315             << endl;
02316     }
02317 }
02318 
02319 
02320 template <class Type>
02321 Foam::indexedOctree<Type>::indexedOctree
02322 (
02323     const Type& shapes,
02324     Istream& is
02325 )
02326 :
02327     shapes_(shapes),
02328     nodes_(is),
02329     contents_(is),
02330     nodeTypes_(0)
02331 {}
02332 
02333 
02334 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
02335 
02336 template <class Type>
02337 Foam::scalar& Foam::indexedOctree<Type>::perturbTol()
02338 {
02339     return perturbTol_;
02340 }
02341 
02342 
02343 template <class Type>
02344 Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
02345 (
02346     const point& sample,
02347     const scalar startDistSqr
02348 ) const
02349 {
02350     scalar nearestDistSqr = startDistSqr;
02351     label nearestShapeI = -1;
02352     point nearestPoint;
02353 
02354     if (nodes_.size())
02355     {
02356         findNearest
02357         (
02358             0,
02359             sample,
02360 
02361             nearestDistSqr,
02362             nearestShapeI,
02363             nearestPoint
02364         );
02365     }
02366     else
02367     {
02368         nearestPoint = vector::zero;
02369     }
02370 
02371     return pointIndexHit(nearestShapeI != -1, nearestPoint, nearestShapeI);
02372 }
02373 
02374 
02375 template <class Type>
02376 Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
02377 (
02378     const linePointRef& ln,
02379     treeBoundBox& tightest,
02380     point& linePoint
02381 ) const
02382 {
02383     label nearestShapeI = -1;
02384     point nearestPoint;
02385 
02386     if (nodes_.size())
02387     {
02388         findNearest
02389         (
02390             0,
02391             ln,
02392 
02393             tightest,
02394             nearestShapeI,
02395             linePoint,
02396             nearestPoint
02397         );
02398     }
02399     else
02400     {
02401         nearestPoint = vector::zero;
02402     }
02403 
02404     return pointIndexHit(nearestShapeI != -1, nearestPoint, nearestShapeI);
02405 }
02406 
02407 
02408 // Find nearest intersection
02409 template <class Type>
02410 Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
02411 (
02412     const point& start,
02413     const point& end
02414 ) const
02415 {
02416     return findLine(false, start, end);
02417 }
02418 
02419 
02420 // Find nearest intersection
02421 template <class Type>
02422 Foam::pointIndexHit Foam::indexedOctree<Type>::findLineAny
02423 (
02424     const point& start,
02425     const point& end
02426 ) const
02427 {
02428     return findLine(true, start, end);
02429 }
02430 
02431 
02432 template <class Type>
02433 Foam::labelList Foam::indexedOctree<Type>::findBox
02434 (
02435     const treeBoundBox& searchBox
02436 ) const
02437 {
02438     // Storage for labels of shapes inside bb. Size estimate.
02439     labelHashSet elements(shapes_.size() / 100);
02440 
02441     if (nodes_.size())
02442     {
02443         findBox(0, searchBox, elements);
02444     }
02445 
02446     return elements.toc();
02447 }
02448 
02449 
02450 // Find node (as parent+octant) containing point
02451 template <class Type>
02452 Foam::labelBits Foam::indexedOctree<Type>::findNode
02453 (
02454     const label nodeI,
02455     const point& sample
02456 ) const
02457 {
02458     if (nodes_.empty())
02459     {
02460         // Empty tree. Return what?
02461         return nodePlusOctant(nodeI, 0);
02462     }
02463 
02464     const node& nod = nodes_[nodeI];
02465 
02466     if (debug)
02467     {
02468         if (!nod.bb_.contains(sample))
02469         {
02470             FatalErrorIn("findNode(..)")
02471                 << "Cannot find " << sample << " in node " << nodeI
02472                 << abort(FatalError);
02473         }
02474     }
02475 
02476     direction octant = nod.bb_.subOctant(sample);
02477 
02478     labelBits index = nod.subNodes_[octant];
02479 
02480     if (isNode(index))
02481     {
02482         // Recurse
02483         return findNode(getNode(index), sample);
02484     }
02485     else if (isContent(index))
02486     {
02487         // Content. Return treenode+octant
02488         return nodePlusOctant(nodeI, octant);
02489     }
02490     else
02491     {
02492         // Empty. Return treenode+octant
02493         return nodePlusOctant(nodeI, octant);
02494     }
02495 }
02496 
02497 
02498 // Determine type (inside/outside/mixed) per node.
02499 template <class Type>
02500 typename Foam::indexedOctree<Type>::volumeType
02501 Foam::indexedOctree<Type>::getVolumeType
02502 (
02503     const point& sample
02504 ) const
02505 {
02506     if (nodes_.empty())
02507     {
02508         return UNKNOWN;
02509     }
02510 
02511     if (nodeTypes_.size() != 8*nodes_.size())
02512     {
02513         // Calculate type for every octant of node.
02514 
02515         nodeTypes_.setSize(8*nodes_.size());
02516         nodeTypes_ = UNKNOWN;
02517 
02518         calcVolumeType(0);
02519 
02520         if (debug)
02521         {
02522             label nUNKNOWN = 0;
02523             label nMIXED = 0;
02524             label nINSIDE = 0;
02525             label nOUTSIDE = 0;
02526 
02527             forAll(nodeTypes_, i)
02528             {
02529                 volumeType type = volumeType(nodeTypes_.get(i));
02530 
02531                 if (type == UNKNOWN)
02532                 {
02533                     nUNKNOWN++;
02534                 }
02535                 else if (type == MIXED)
02536                 {
02537                     nMIXED++;
02538                 }
02539                 else if (type == INSIDE)
02540                 {
02541                     nINSIDE++;
02542                 }
02543                 else if (type == OUTSIDE)
02544                 {
02545                     nOUTSIDE++;
02546                 }
02547                 else
02548                 {
02549                     FatalErrorIn("getVolumeType") << abort(FatalError);
02550                 }
02551             }
02552 
02553             Pout<< "indexedOctree<Type>::getVolumeType : "
02554                 << " bb:" << bb()
02555                 << " nodes_:" << nodes_.size()
02556                 << " nodeTypes_:" << nodeTypes_.size()
02557                 << " nUNKNOWN:" << nUNKNOWN
02558                 << " nMIXED:" << nMIXED
02559                 << " nINSIDE:" << nINSIDE
02560                 << " nOUTSIDE:" << nOUTSIDE
02561                 << endl;
02562         }
02563     }
02564 
02565     return getVolumeType(0, sample);
02566 }
02567 
02568 
02569 // Print contents of nodeI
02570 template <class Type>
02571 void Foam::indexedOctree<Type>::print
02572 (
02573     prefixOSstream& os,
02574     const bool printContents,
02575     const label nodeI
02576 ) const
02577 {
02578     const node& nod = nodes_[nodeI];
02579     const treeBoundBox& bb = nod.bb_;
02580 
02581     os  << "nodeI:" << nodeI << " bb:" << bb << nl
02582         << "parent:" << nod.parent_ << nl
02583         << "n:" << countElements(nodePlusOctant(nodeI, 0)) << nl;
02584 
02585     for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
02586     {
02587         const treeBoundBox subBb(bb.subBbox(octant));
02588 
02589         labelBits index = nod.subNodes_[octant];
02590 
02591         if (isNode(index))
02592         {
02593             // tree node.
02594             label subNodeI = getNode(index);
02595 
02596             os  << "octant:" << octant
02597                 << " node: n:" << countElements(index)
02598                 << " bb:" << subBb << endl;
02599 
02600             string oldPrefix = os.prefix();
02601             os.prefix() = "  " + oldPrefix;
02602 
02603             print(os, printContents, subNodeI);
02604 
02605             os.prefix() = oldPrefix;
02606         }
02607         else if (isContent(index))
02608         {
02609             const labelList& indices = contents_[getContent(index)];
02610 
02611             os  << "octant:" << octant
02612                 << " content: n:" << indices.size()
02613                 << " bb:" << subBb;
02614 
02615             if (printContents)
02616             {
02617                 os << " contents:";
02618                 forAll(indices, j)
02619                 {
02620                     os  << ' ' << indices[j];
02621                 }
02622             }
02623             os  << endl;
02624         }
02625         else
02626         {
02627             os  << "octant:" << octant << " empty:" << subBb << endl;
02628         }
02629     }
02630 }
02631 
02632 
02633 // Print contents of nodeI
02634 template <class Type>
02635 bool Foam::indexedOctree<Type>::write(Ostream& os) const
02636 {
02637     os << *this;
02638 
02639     return os.good();
02640 }
02641 
02642 
02643 template <class Type>
02644 Foam::Ostream& Foam::operator<<(Ostream& os, const indexedOctree<Type>& t)
02645 {
02646     return
02647         os  << t.bb() << token::SPACE << t.nodes()
02648             << token::SPACE << t.contents();
02649 }
02650 
02651 
02652 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines