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