src/meshTools/meshSearch/meshSearch.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "meshSearch.H"
  27 #include "polyMesh.H"
  28 #include "indexedOctree.H"
  29 #include "DynamicList.H"
  30 #include "demandDrivenData.H"
  31 #include "treeDataPolyMeshCell.H"
  32 #include "treeDataFace.H"
  33
  34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  35
  36 defineTypeNameAndDebug(Foam::meshSearch, 0);
  37
  38 Foam::scalar Foam::meshSearch::tol_ = 1E-3;
  39
  40
  41 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  42
  43 bool Foam::meshSearch::findNearer
  44 (
  45     const point& sample,
  46     const pointField& points,
  47     label& nearestI,
  48     scalar& nearestDistSqr
  49 )
  50 {
  51     bool nearer = false;
  52
  53     forAll(points, pointI)
  54     {
  55         scalar distSqr = magSqr(points[pointI] - sample);
  56
  57         if (distSqr < nearestDistSqr)
  58         {
  59             nearestDistSqr = distSqr;
  60             nearestI = pointI;
  61             nearer = true;
  62         }
  63     }
  64
  65     return nearer;
  66 }
  67
  68
  69 bool Foam::meshSearch::findNearer
  70 (
  71     const point& sample,
  72     const pointField& points,
  73     const labelList& indices,
  74     label& nearestI,
  75     scalar& nearestDistSqr
  76 )
  77 {
  78     bool nearer = false;
  79
  80     forAll(indices, i)
  81     {
  82         label pointI = indices[i];
  83
  84         scalar distSqr = magSqr(points[pointI] - sample);
  85
  86         if (distSqr < nearestDistSqr)
  87         {
  88             nearestDistSqr = distSqr;
  89             nearestI = pointI;
  90             nearer = true;
  91         }
  92     }
  93
  94     return nearer;
  95 }
  96
  97
  98 // tree based searching
  99 Foam::label Foam::meshSearch::findNearestCellTree(const point& location) const
 100 {
 101     const indexedOctree<treeDataPolyMeshCell>& tree = cellTree();
 102
 103     pointIndexHit info = tree.findNearest
 104     (
 105         location,
 106         magSqr(tree.bb().max()-tree.bb().min())
 107     );
 108
 109     if (!info.hit())
 110     {
 111         info = tree.findNearest(location, Foam::sqr(GREAT));
 112     }
 113     return info.index();
 114 }
 115
 116
 117 // linear searching
 118 Foam::label Foam::meshSearch::findNearestCellLinear(const point& location) const
 119 {
 120     const vectorField& centres = mesh_.cellCentres();
 121
 122     label nearestIndex = 0;
 123     scalar minProximity = magSqr(centres[nearestIndex] - location);
 124
 125     findNearer
 126     (
 127         location,
 128         centres,
 129         nearestIndex,
 130         minProximity
 131     );
 132
 133     return nearestIndex;
 134 }
 135
 136
 137 // walking from seed
 138 Foam::label Foam::meshSearch::findNearestCellWalk
 139 (
 140     const point& location,
 141     const label seedCellI
 142 ) const
 143 {
 144     if (seedCellI < 0)
 145     {
 146         FatalErrorIn
 147         (
 148             "meshSearch::findNearestCellWalk(const point&, const label)"
 149         )   << "illegal seedCell:" << seedCellI << exit(FatalError);
 150     }
 151
 152     // Walk in direction of face that decreases distance
 153
 154     label curCellI = seedCellI;
 155     scalar distanceSqr = magSqr(mesh_.cellCentres()[curCellI] - location);
 156
 157     bool closer;
 158
 159     do
 160     {
 161         // Try neighbours of curCellI
 162         closer = findNearer
 163         (
 164             location,
 165             mesh_.cellCentres(),
 166             mesh_.cellCells()[curCellI],
 167             curCellI,
 168             distanceSqr
 169         );
 170     } while (closer);
 171
 172     return curCellI;
 173 }
 174
 175
 176 // tree based searching
 177 Foam::label Foam::meshSearch::findNearestFaceTree(const point& location) const
 178 {
 179     // Search nearest cell centre.
 180     const indexedOctree<treeDataPolyMeshCell>& tree = cellTree();
 181
 182     // Search with decent span
 183     pointIndexHit info = tree.findNearest
 184     (
 185         location,
 186         magSqr(tree.bb().max()-tree.bb().min())
 187     );
 188
 189     if (!info.hit())
 190     {
 191         // Search with desparate span
 192         info = tree.findNearest(location, Foam::sqr(GREAT));
 193     }
 194
 195
 196     // Now check any of the faces of the nearest cell
 197     const vectorField& centres = mesh_.faceCentres();
 198     const cell& ownFaces = mesh_.cells()[info.index()];
 199
 200     label nearestFaceI = ownFaces[0];
 201     scalar minProximity = magSqr(centres[nearestFaceI] - location);
 202
 203     findNearer
 204     (
 205         location,
 206         centres,
 207         ownFaces,
 208         nearestFaceI,
 209         minProximity
 210     );
 211
 212     return nearestFaceI;
 213 }
 214
 215
 216 // linear searching
 217 Foam::label Foam::meshSearch::findNearestFaceLinear(const point& location) const
 218 {
 219     const vectorField& centres = mesh_.faceCentres();
 220
 221     label nearestFaceI = 0;
 222     scalar minProximity = magSqr(centres[nearestFaceI] - location);
 223
 224     findNearer
 225     (
 226         location,
 227         centres,
 228         nearestFaceI,
 229         minProximity
 230     );
 231
 232     return nearestFaceI;
 233 }
 234
 235
 236 // walking from seed
 237 Foam::label Foam::meshSearch::findNearestFaceWalk
 238 (
 239     const point& location,
 240     const label seedFaceI
 241 ) const
 242 {
 243     if (seedFaceI < 0)
 244     {
 245         FatalErrorIn
 246         (
 247             "meshSearch::findNearestFaceWalk(const point&, const label)"
 248         )   << "illegal seedFace:" << seedFaceI << exit(FatalError);
 249     }
 250
 251     const vectorField& centres = mesh_.faceCentres();
 252
 253
 254     // Walk in direction of face that decreases distance
 255
 256     label curFaceI = seedFaceI;
 257     scalar distanceSqr = magSqr(centres[curFaceI] - location);
 258
 259     while (true)
 260     {
 261         label betterFaceI = curFaceI;
 262
 263         findNearer
 264         (
 265             location,
 266             centres,
 267             mesh_.cells()[mesh_.faceOwner()[curFaceI]],
 268             betterFaceI,
 269             distanceSqr
 270         );
 271
 272         if (mesh_.isInternalFace(curFaceI))
 273         {
 274             findNearer
 275             (
 276                 location,
 277                 centres,
 278                 mesh_.cells()[mesh_.faceNeighbour()[curFaceI]],
 279                 betterFaceI,
 280                 distanceSqr
 281             );
 282         }
 283
 284         if (betterFaceI == curFaceI)
 285         {
 286             break;
 287         }
 288
 289         curFaceI = betterFaceI;
 290     }
 291
 292     return curFaceI;
 293 }
 294
 295
 296 Foam::label Foam::meshSearch::findCellLinear(const point& location) const
 297 {
 298     bool cellFound = false;
 299     label n = 0;
 300
 301     label cellI = -1;
 302
 303     while ((!cellFound) && (n < mesh_.nCells()))
 304     {
 305         if (pointInCell(location, n))
 306         {
 307             cellFound = true;
 308             cellI = n;
 309         }
 310         else
 311         {
 312             n++;
 313         }
 314     }
 315     if (cellFound)
 316     {
 317         return cellI;
 318     }
 319     else
 320     {
 321         return -1;
 322     }
 323 }
 324
 325
 326 // walking from seed
 327 Foam::label Foam::meshSearch::findCellWalk
 328 (
 329     const point& location,
 330     const label seedCellI
 331 ) const
 332 {
 333     if (seedCellI < 0)
 334     {
 335         FatalErrorIn
 336         (
 337             "meshSearch::findCellWalk(const point&, const label)"
 338         )   << "illegal seedCell:" << seedCellI << exit(FatalError);
 339     }
 340
 341     if (pointInCell(location, seedCellI))
 342     {
 343         return seedCellI;
 344     }
 345
 346     // Walk in direction of face that decreases distance
 347     label curCellI = seedCellI;
 348     scalar nearestDistSqr = magSqr(mesh_.cellCentres()[curCellI] - location);
 349
 350     while(true)
 351     {
 352         // Try neighbours of curCellI
 353
 354         const cell& cFaces = mesh_.cells()[curCellI];
 355
 356         label nearestCellI = -1;
 357
 358         forAll(cFaces, i)
 359         {
 360             label faceI = cFaces[i];
 361
 362             if (mesh_.isInternalFace(faceI))
 363             {
 364                 label cellI = mesh_.faceOwner()[faceI];
 365                 if (cellI == curCellI)
 366                 {
 367                     cellI = mesh_.faceNeighbour()[faceI];
 368                 }
 369
 370                 // Check if this is the correct cell
 371                 if (pointInCell(location, cellI))
 372                 {
 373                     return cellI;
 374                 }
 375
 376                 // Also calculate the nearest cell
 377                 scalar distSqr = magSqr(mesh_.cellCentres()[cellI] - location);
 378
 379                 if (distSqr < nearestDistSqr)
 380                 {
 381                     nearestDistSqr = distSqr;
 382                     nearestCellI = cellI;
 383                 }
 384             }
 385         }
 386
 387         if (nearestCellI == -1)
 388         {
 389             return -1;
 390         }
 391
 392         // Continue with the nearest cell
 393         curCellI = nearestCellI;
 394     }
 395
 396     return -1;
 397 }
 398
 399
 400 Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
 401 (
 402     const point& location,
 403     const label seedFaceI
 404 ) const
 405 {
 406     if (seedFaceI < 0)
 407     {
 408         FatalErrorIn
 409         (
 410             "meshSearch::findNearestBoundaryFaceWalk"
 411             "(const point&, const label)"
 412         )   << "illegal seedFace:" << seedFaceI << exit(FatalError);
 413     }
 414
 415     // Start off from seedFaceI
 416
 417     label curFaceI = seedFaceI;
 418
 419     const face& f =  mesh_.faces()[curFaceI];
 420
 421     scalar minDist = f.nearestPoint
 422     (
 423         location,
 424         mesh_.points()
 425     ).distance();
 426
 427     bool closer;
 428
 429     do
 430     {
 431         closer = false;
 432
 433         // Search through all neighbouring boundary faces by going
 434         // across edges
 435
 436         label lastFaceI = curFaceI;
 437
 438         const labelList& myEdges = mesh_.faceEdges()[curFaceI];
 439
 440         forAll(myEdges, myEdgeI)
 441         {
 442             const labelList& neighbours = mesh_.edgeFaces()[myEdges[myEdgeI]];
 443
 444             // Check any face which uses edge, is boundary face and
 445             // is not curFaceI itself.
 446
 447             forAll(neighbours, nI)
 448             {
 449                 label faceI = neighbours[nI];
 450
 451                 if
 452                 (
 453                     (faceI >= mesh_.nInternalFaces())
 454                  && (faceI != lastFaceI)
 455                 )
 456                 {
 457                     const face& f =  mesh_.faces()[faceI];
 458
 459                     pointHit curHit = f.nearestPoint
 460                     (
 461                         location,
 462                         mesh_.points()
 463                     );
 464
 465                     // If the face is closer, reset current face and distance
 466                     if (curHit.distance() < minDist)
 467                     {
 468                         minDist = curHit.distance();
 469                         curFaceI = faceI;
 470                         closer = true;  // a closer neighbour has been found
 471                     }
 472                 }
 473             }
 474         }
 475     } while (closer);
 476
 477     return curFaceI;
 478 }
 479
 480
 481 Foam::vector Foam::meshSearch::offset
 482 (
 483     const point& bPoint,
 484     const label bFaceI,
 485     const vector& dir
 486 ) const
 487 {
 488     // Get the neighbouring cell
 489     label ownerCellI = mesh_.faceOwner()[bFaceI];
 490
 491     const point& c = mesh_.cellCentres()[ownerCellI];
 492
 493     // Typical dimension: distance from point on face to cell centre
 494     scalar typDim = mag(c - bPoint);
 495
 496     return tol_*typDim*dir;
 497 }
 498
 499
 500 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 501
 502 // Construct from components
 503 Foam::meshSearch::meshSearch(const polyMesh& mesh, const bool faceDecomp)
 504 :
 505     mesh_(mesh),
 506     faceDecomp_(faceDecomp)
 507 {}
 508
 509
 510 // Construct with a custom bounding box
 511 Foam::meshSearch::meshSearch
 512 (
 513     const polyMesh& mesh,
 514     const treeBoundBox& bb,
 515     const bool faceDecomp
 516 )
 517 :
 518     mesh_(mesh),
 519     faceDecomp_(faceDecomp)
 520 {
 521     overallBbPtr_.reset(new treeBoundBox(bb));
 522 }
 523
 524
 525 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 526
 527 Foam::meshSearch::~meshSearch()
 528 {
 529     clearOut();
 530 }
 531
 532
 533 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 534
 535 const Foam::indexedOctree<Foam::treeDataFace>& Foam::meshSearch::boundaryTree()
 536  const
 537 {
 538     if (!boundaryTreePtr_.valid())
 539     {
 540         //
 541         // Construct tree
 542         //
 543
 544         if (!overallBbPtr_.valid())
 545         {
 546             Random rndGen(261782);
 547             overallBbPtr_.reset
 548             (
 549                 new treeBoundBox(mesh_.points())
 550             );
 551
 552             treeBoundBox& overallBb = overallBbPtr_();
 553             // Extend slightly and make 3D
 554             overallBb = overallBb.extend(rndGen, 1E-4);
 555             overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 556             overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 557         }
 558
 559         // all boundary faces (not just walls)
 560         labelList bndFaces(mesh_.nFaces()-mesh_.nInternalFaces());
 561         forAll(bndFaces, i)
 562         {
 563             bndFaces[i] = mesh_.nInternalFaces() + i;
 564         }
 565
 566         boundaryTreePtr_.reset
 567         (
 568             new indexedOctree<treeDataFace>
 569             (
 570                 treeDataFace    // all information needed to search faces
 571                 (
 572                     false,                      // do not cache bb
 573                     mesh_,
 574                     bndFaces                    // boundary faces only
 575                 ),
 576                 overallBbPtr_(),                // overall search domain
 577                 8,                              // maxLevel
 578                 10,                             // leafsize
 579                 3.0                             // duplicity
 580             )
 581         );
 582     }
 583
 584     return boundaryTreePtr_();
 585 }
 586
 587
 588 const Foam::indexedOctree<Foam::treeDataPolyMeshCell>&
 589 Foam::meshSearch::cellTree() const
 590 {
 591     if (!cellTreePtr_.valid())
 592     {
 593         //
 594         // Construct tree
 595         //
 596
 597         if (!overallBbPtr_.valid())
 598         {
 599             Random rndGen(261782);
 600             overallBbPtr_.reset
 601             (
 602                 new treeBoundBox(mesh_.points())
 603             );
 604
 605             treeBoundBox& overallBb = overallBbPtr_();
 606             // Extend slightly and make 3D
 607             overallBb = overallBb.extend(rndGen, 1E-4);
 608             overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 609             overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 610         }
 611
 612         cellTreePtr_.reset
 613         (
 614             new indexedOctree<treeDataPolyMeshCell>
 615             (
 616                 treeDataPolyMeshCell
 617                 (
 618                     false,  // not cache bb
 619                     mesh_
 620                 ),
 621                 overallBbPtr_(),
 622                 8,              // maxLevel
 623                 10,             // leafsize
 624                 6.0             // duplicity
 625             )
 626         );
 627     }
 628
 629     return cellTreePtr_();
 630 }
 631
 632
 633 // Is the point in the cell
 634 // Works by checking if there is a face inbetween the point and the cell
 635 // centre.
 636 // Check for internal uses proper face decomposition or just average normal.
 637 bool Foam::meshSearch::pointInCell(const point& p, label cellI) const
 638 {
 639     if (faceDecomp_)
 640     {
 641         const point& ctr = mesh_.cellCentres()[cellI];
 642
 643         vector dir(p - ctr);
 644         scalar magDir = mag(dir);
 645
 646         // Check if any faces are hit by ray from cell centre to p.
 647         // If none -> p is in cell.
 648         const labelList& cFaces = mesh_.cells()[cellI];
 649
 650         // Make sure half_ray does not pick up any faces on the wrong
 651         // side of the ray.
 652         scalar oldTol = intersection::setPlanarTol(0.0);
 653
 654         forAll(cFaces, i)
 655         {
 656             label faceI = cFaces[i];
 657
 658             pointHit inter = mesh_.faces()[faceI].ray
 659             (
 660                 ctr,
 661                 dir,
 662                 mesh_.points(),
 663                 intersection::HALF_RAY,
 664                 intersection::VECTOR
 665             );
 666
 667             if (inter.hit())
 668             {
 669                 scalar dist = inter.distance();
 670
 671                 if (dist < magDir)
 672                 {
 673                     // Valid hit. Hit face so point is not in cell.
 674                     intersection::setPlanarTol(oldTol);
 675
 676                     return false;
 677                 }
 678             }
 679         }
 680
 681         intersection::setPlanarTol(oldTol);
 682
 683         // No face inbetween point and cell centre so point is inside.
 684         return true;
 685     }
 686     else
 687     {
 688         const labelList& f = mesh_.cells()[cellI];
 689         const labelList& owner = mesh_.faceOwner();
 690         const vectorField& cf = mesh_.faceCentres();
 691         const vectorField& Sf = mesh_.faceAreas();
 692
 693         forAll(f, facei)
 694         {
 695             label nFace = f[facei];
 696             vector proj = p - cf[nFace];
 697             vector normal = Sf[nFace];
 698             if (owner[nFace] == cellI)
 699             {
 700                 if ((normal & proj) > 0)
 701                 {
 702                     return false;
 703                 }
 704             }
 705             else
 706             {
 707                 if ((normal & proj) < 0)
 708                 {
 709                     return false;
 710                 }
 711             }
 712         }
 713
 714         return true;
 715     }
 716 }
 717
 718
 719 Foam::label Foam::meshSearch::findNearestCell
 720 (
 721     const point& location,
 722     const label seedCellI,
 723     const bool useTreeSearch
 724 ) const
 725 {
 726     if (seedCellI == -1)
 727     {
 728         if (useTreeSearch)
 729         {
 730             return findNearestCellTree(location);
 731         }
 732         else
 733         {
 734             return findNearestCellLinear(location);
 735         }
 736     }
 737     else
 738     {
 739         return findNearestCellWalk(location, seedCellI);
 740     }
 741 }
 742
 743
 744 Foam::label Foam::meshSearch::findNearestFace
 745 (
 746     const point& location,
 747     const label seedFaceI,
 748     const bool useTreeSearch
 749 ) const
 750 {
 751     if (seedFaceI == -1)
 752     {
 753         if (useTreeSearch)
 754         {
 755             return findNearestFaceTree(location);
 756         }
 757         else
 758         {
 759             return findNearestFaceLinear(location);
 760         }
 761     }
 762     else
 763     {
 764         return findNearestFaceWalk(location, seedFaceI);
 765     }
 766 }
 767
 768
 769 Foam::label Foam::meshSearch::findCell
 770 (
 771     const point& location,
 772     const label seedCellI,
 773     const bool useTreeSearch
 774 ) const
 775 {
 776     // Find the nearest cell centre to this location
 777     if (seedCellI == -1)
 778     {
 779         if (useTreeSearch)
 780         {
 781             return cellTree().findInside(location);
 782         }
 783         else
 784         {
 785              return findCellLinear(location);
 786         }
 787     }
 788     else
 789     {
 790         return findCellWalk(location, seedCellI);
 791     }
 792 }
 793
 794
 795 Foam::label Foam::meshSearch::findNearestBoundaryFace
 796 (
 797     const point& location,
 798     const label seedFaceI,
 799     const bool useTreeSearch
 800 ) const
 801 {
 802     if (seedFaceI == -1)
 803     {
 804         if (useTreeSearch)
 805         {
 806             const indexedOctree<treeDataFace>& tree = boundaryTree();
 807
 808             pointIndexHit info = boundaryTree().findNearest
 809             (
 810                 location,
 811                 magSqr(tree.bb().max()-tree.bb().min())
 812             );
 813
 814             if (!info.hit())
 815             {
 816                 info = boundaryTree().findNearest
 817                 (
 818                     location,
 819                     Foam::sqr(GREAT)
 820                 );
 821             }
 822
 823             return tree.shapes().faceLabels()[info.index()];
 824         }
 825         else
 826         {
 827             scalar minDist = GREAT;
 828
 829             label minFaceI = -1;
 830
 831             for
 832             (
 833                 label faceI = mesh_.nInternalFaces();
 834                 faceI < mesh_.nFaces();
 835                 faceI++
 836             )
 837             {
 838                 const face& f =  mesh_.faces()[faceI];
 839
 840                 pointHit curHit =
 841                     f.nearestPoint
 842                     (
 843                         location,
 844                         mesh_.points()
 845                     );
 846
 847                 if (curHit.distance() < minDist)
 848                 {
 849                     minDist = curHit.distance();
 850                     minFaceI = faceI;
 851                 }
 852             }
 853             return minFaceI;
 854         }
 855     }
 856     else
 857     {
 858         return findNearestBoundaryFaceWalk(location, seedFaceI);
 859     }
 860 }
 861
 862
 863 Foam::pointIndexHit Foam::meshSearch::intersection
 864 (
 865     const point& pStart,
 866     const point& pEnd
 867 ) const
 868 {
 869     pointIndexHit curHit = boundaryTree().findLine(pStart, pEnd);
 870
 871     if (curHit.hit())
 872     {
 873         // Change index into octreeData into face label
 874         curHit.setIndex(boundaryTree().shapes().faceLabels()[curHit.index()]);
 875     }
 876     return curHit;
 877 }
 878
 879
 880 Foam::List<Foam::pointIndexHit> Foam::meshSearch::intersections
 881 (
 882     const point& pStart,
 883     const point& pEnd
 884 ) const
 885 {
 886     DynamicList<pointIndexHit> hits;
 887
 888     vector edgeVec = pEnd - pStart;
 889     edgeVec /= mag(edgeVec);
 890
 891     point pt = pStart;
 892
 893     pointIndexHit bHit;
 894     do
 895     {
 896         bHit = intersection(pt, pEnd);
 897
 898         if (bHit.hit())
 899         {
 900             hits.append(bHit);
 901
 902             const vector& area = mesh_.faceAreas()[bHit.index()];
 903
 904             scalar typDim = Foam::sqrt(mag(area));
 905
 906             if ((mag(bHit.hitPoint() - pEnd)/typDim) < SMALL)
 907             {
 908                 break;
 909             }
 910
 911             // Restart from hitPoint shifted a little bit in the direction
 912             // of the destination
 913
 914             pt =
 915                 bHit.hitPoint()
 916               + offset(bHit.hitPoint(), bHit.index(), edgeVec);
 917         }
 918
 919     } while (bHit.hit());
 920
 921
 922     hits.shrink();
 923
 924     return hits;
 925 }
 926
 927
 928 bool Foam::meshSearch::isInside(const point& p) const
 929 {
 930     return
 931         boundaryTree().getVolumeType(p)
 932      == indexedOctree<treeDataFace>::INSIDE;
 933 }
 934
 935
 936 // Delete all storage
 937 void Foam::meshSearch::clearOut()
 938 {
 939     boundaryTreePtr_.clear();
 940     cellTreePtr_.clear();
 941     overallBbPtr_.clear();
 942 }
 943
 944
 945 void Foam::meshSearch::correct()
 946 {
 947     clearOut();
 948 }
 949
 950
 951 // ************************************************************************* //