src/meshTools/meshSearch/meshSearch.C

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