src/meshTools/meshSearch/meshSearch.C

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