src/meshTools/meshSearch/meshSearch.C

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