src/meshLibrary/utilities/octrees/meshOctree/meshOctreeFindNearestSurfacePoint.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | cfMesh: A library for mesh generation
   4    \\    /   O peration     |
   5     \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
   6      \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of cfMesh.
  10
  11     cfMesh 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     cfMesh 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 cfMesh.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include "meshOctree.H"
  29 #include "triSurf.H"
  30 #include "demandDrivenData.H"
  31 #include "helperFunctions.H"
  32 #include "HashSet.H"
  33
  34 // #define DEBUGSearch
  35
  36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40
  41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  42
  43 void meshOctree::findNearestSurfacePoint
  44 (
  45     point& nearest,
  46     scalar& distSq,
  47     label& nearestTriangle,
  48     label& region,
  49     const point& p
  50 ) const
  51 {
  52     region = -1;
  53     nearestTriangle = 1;
  54
  55     const label cLabel = findLeafContainingVertex(p);
  56     vector sizeVec;
  57     if( cLabel < 0 )
  58     {
  59         sizeVec.x() = sizeVec.y() = sizeVec.z() = searchRange_;
  60     }
  61     else
  62     {
  63         const scalar s = 0.75 * leaves_[cLabel]->size(rootBox_);
  64         sizeVec.x() = sizeVec.y() = sizeVec.z() = s;
  65     }
  66
  67     //- find nearest surface vertex to the point p
  68     bool found(false);
  69     label iterationI(0);
  70     DynList<const meshOctreeCube*, 256> neighbours;
  71
  72     distSq = VGREAT;
  73
  74     do
  75     {
  76         boundBox bb(p - sizeVec, p + sizeVec);
  77
  78         neighbours.clear();
  79         findLeavesContainedInBox(bb, neighbours);
  80
  81         //- find nearest projection
  82         forAll(neighbours, neiI)
  83         {
  84             if( !neighbours[neiI]->hasContainedElements() )
  85                 continue;
  86
  87             const VRWGraph& ct =
  88                 neighbours[neiI]->slotPtr()->containedTriangles_;
  89             const constRow el = ct[neighbours[neiI]->containedElements()];
  90             forAll(el, tI)
  91             {
  92                 const point p0 =
  93                     help::nearestPointOnTheTriangle(el[tI], surface_, p);
  94
  95                 const scalar dSq = Foam::magSqr(p0 - p);
  96                 if( dSq < distSq )
  97                 {
  98                     distSq = dSq;
  99                     nearest = p0;
 100                     nearestTriangle = el[tI];
 101                     region = surface_[el[tI]].region();
 102                     found = true;
 103                 }
 104             }
 105         }
 106
 107         if( !found )
 108             sizeVec *= 2.0;
 109
 110     } while( !found && (iterationI++ < 5) );
 111
 112     # ifdef DEBUGSearch
 113     forAll(surface_, triI)
 114     {
 115         const point pp = help::nearestPointOnTheTriangle(triI, surface_, p);
 116
 117         if( distSq - magSqr(pp - p) > SMALL )
 118             Pout << "Point " << p << " current nearest " << nearest
 119                  << " closer point " << pp << endl;
 120     }
 121     # endif
 122
 123     if( (!found || (region < 0)) && !Pstream::parRun() )
 124     {
 125         Warning << "Could not find a boundary region for vertex " << p << endl;
 126         Warning << "Found " << found << " and region " << region << endl;
 127     }
 128 }
 129
 130 void meshOctree::findNearestSurfacePointInRegion
 131 (
 132     point& nearest,
 133     scalar& distSq,
 134     label& nearestTriangle,
 135     const label region,
 136     const point& p
 137 ) const
 138 {
 139     const label cLabel = findLeafContainingVertex(p);
 140     vector sizeVec;
 141     if( cLabel < 0 )
 142     {
 143         sizeVec.x() = sizeVec.y() = sizeVec.z() = searchRange_;
 144     }
 145     else
 146     {
 147         const scalar s = 0.75 * leaves_[cLabel]->size(rootBox_);
 148         sizeVec.x() = sizeVec.y() = sizeVec.z() = s;
 149     }
 150
 151     //- find nearest surface vertex to the point p
 152     bool found(false);
 153     label iterationI(0);
 154     DynList<const meshOctreeCube*, 256> neighbours;
 155     nearestTriangle = -1;
 156
 157     distSq = VGREAT;
 158
 159     do
 160     {
 161         const boundBox bb(p - sizeVec, p + sizeVec);
 162
 163         neighbours.clear();
 164         findLeavesContainedInBox(bb, neighbours);
 165
 166         //- find nearest projection
 167         forAll(neighbours, neiI)
 168         {
 169             if( !neighbours[neiI]->hasContainedElements() )
 170                 continue;
 171
 172             const VRWGraph& ct =
 173                 neighbours[neiI]->slotPtr()->containedTriangles_;
 174             const constRow el =
 175                 ct[neighbours[neiI]->containedElements()];
 176             forAll(el, tI)
 177             {
 178                 if( surface_[el[tI]].region() != region )
 179                     continue;
 180
 181                 const point p0 =
 182                     help::nearestPointOnTheTriangle(el[tI], surface_, p);
 183
 184                 const scalar dSq = Foam::magSqr(p0 - p);
 185                 if( dSq < distSq )
 186                 {
 187                     distSq = dSq;
 188                     nearest = p0;
 189                     nearestTriangle = el[tI];
 190                     found = true;
 191                 }
 192             }
 193         }
 194
 195         if( !found )
 196             sizeVec *= 2.0;
 197
 198     } while( !found && (iterationI++ < 5) );
 199
 200     if( (!found || (region < 0)) && !Pstream::parRun() )
 201         Warning << "Could not find a boundary region for vertex " << p << endl;
 202 }
 203
 204 bool meshOctree::findNearestEdgePoint
 205 (
 206     point& edgePoint,
 207     scalar& distSq,
 208     label& nearestEdge,
 209     const point& p,
 210     const DynList<label>& regions
 211 ) const
 212 {
 213     //- find the estimate for the searching range
 214     const label cLabel = findLeafContainingVertex(p);
 215     vector sizeVec;
 216     if( cLabel < 0 )
 217     {
 218         sizeVec.x() = sizeVec.y() = sizeVec.z() = searchRange_;
 219     }
 220     else
 221     {
 222         const scalar s = 0.75 * leaves_[cLabel]->size(rootBox_);
 223         sizeVec.x() = sizeVec.y() = sizeVec.z() = s;
 224     }
 225
 226     DynList<const meshOctreeCube*, 256> neighbours;
 227
 228     const pointField& sp = surface_.points();
 229     const edgeLongList& edges = surface_.edges();
 230     const VRWGraph& edgeFaces = surface_.edgeFacets();
 231
 232     edgePoint = p;
 233     bool foundAnEdge(false);
 234     label iterationI(0);
 235
 236     distSq = VGREAT;
 237     nearestEdge = -1;
 238
 239     do
 240     {
 241         boundBox bb(p - sizeVec, p + sizeVec);
 242
 243         neighbours.clear();
 244         findLeavesContainedInBox(bb, neighbours);
 245
 246         forAll(neighbours, neiI)
 247         {
 248             if( !neighbours[neiI]->hasContainedEdges() )
 249                 continue;
 250
 251             const VRWGraph& containedEdges =
 252                 neighbours[neiI]->slotPtr()->containedEdges_;
 253             const constRow ce =
 254                 containedEdges[neighbours[neiI]->containedEdges()];
 255
 256             forAll(ce, eI)
 257             {
 258                 //- find if the edge is in correct patches
 259                 bool correctPatches(true);
 260
 261                 forAllRow(edgeFaces, ce[eI], efI)
 262                 {
 263                     const label facetI = edgeFaces(ce[eI], efI);
 264
 265                     if( !regions.contains(surface_[facetI].region()) )
 266                     {
 267                         correctPatches = false;
 268                         break;
 269                     }
 270                 }
 271
 272                 if( !correctPatches )
 273                     continue;
 274
 275                 const point s = sp[edges[ce[eI]].start()];
 276                 const point e = sp[edges[ce[eI]].end()];
 277                 const point np = help::nearestPointOnTheEdgeExact(s, e, p);
 278                 const scalar dSq = Foam::magSqr(np - p);
 279
 280                 if( dSq < distSq )
 281                 {
 282                     distSq = dSq;
 283                     edgePoint = np;
 284                     nearestEdge = ce[eI];
 285                     foundAnEdge = true;
 286                 }
 287             }
 288         }
 289
 290         if( !foundAnEdge )
 291             sizeVec *= 2.0;
 292
 293     } while( !foundAnEdge && (++iterationI < 5) );
 294
 295     return foundAnEdge;
 296 }
 297
 298 bool meshOctree::findNearestPointToEdge
 299 (
 300     point& nearest,
 301     scalar& distSq,
 302     label& nearestEdge,
 303     const FixedList<point, 2>& edgePoints,
 304     const FixedList<label, 2>& edgePointRegions
 305 ) const
 306 {
 307     const point c = 0.5 * (edgePoints[0] + edgePoints[1]);
 308     const scalar dst = mag(edgePoints[0] - edgePoints[1]);
 309     vector sizeVec(dst, dst, dst);
 310
 311     boundBox bb(c - 0.75 * sizeVec, c + 0.75 * sizeVec);
 312
 313     DynList<const meshOctreeCube*, 256> leavesInBox;
 314     findLeavesContainedInBox(bb, leavesInBox);
 315
 316     const VRWGraph& edgeFaces = surface_.edgeFacets();
 317     const pointField& points = surface_.points();
 318     const edgeLongList& surfaceEdges = surface_.edges();
 319
 320     distSq = VGREAT;
 321     nearestEdge = -1;
 322
 323     bool found(false);
 324
 325     forAll(leavesInBox, leafI)
 326     {
 327         if( !leavesInBox[leafI]->hasContainedEdges() )
 328             continue;
 329
 330         const VRWGraph& containedEdges =
 331             leavesInBox[leafI]->slotPtr()->containedEdges_;
 332         const constRow edges =
 333             containedEdges[leavesInBox[leafI]->containedEdges()];
 334
 335         forAll(edges, eI)
 336         {
 337             const constRow ef = edgeFaces[edges[eI]];
 338             if( ef.size() != 2 )
 339                 continue;
 340
 341             if(
 342                 (
 343                     (surface_[ef[0]].region() == edgePointRegions[0]) &&
 344                     (surface_[ef[1]].region() == edgePointRegions[1])
 345                 ) ||
 346                 (
 347                     (surface_[ef[1]].region() == edgePointRegions[0]) &&
 348                     (surface_[ef[0]].region() == edgePointRegions[1])
 349                 )
 350             )
 351             {
 352                 const edge& edg = surfaceEdges[edges[eI]];
 353
 354                 point nearestOnEdge, nearestOnLine;
 355                 if(
 356                     help::nearestEdgePointToTheLine
 357                     (
 358                         points[edg[0]],
 359                         points[edg[1]],
 360                         edgePoints[0],
 361                         edgePoints[1],
 362                         nearestOnEdge,
 363                         nearestOnLine
 364                     )
 365                 )
 366                 {
 367                     if( magSqr(nearestOnEdge - nearestOnLine) < distSq )
 368                     {
 369                         nearest = nearestOnEdge;
 370                         nearestEdge = edges[eI];
 371                         distSq = magSqr(nearestOnEdge - nearestOnLine);
 372                         found = true;
 373                     }
 374                 }
 375             }
 376         }
 377     }
 378
 379     return found;
 380 }
 381
 382 bool meshOctree::findNearestCorner
 383 (
 384     point& nearest,
 385     scalar& distSq,
 386     label& nearestPoint,
 387     const point& p,
 388     const DynList<label>& patches
 389 ) const
 390 {
 391
 392
 393     const label cLabel = findLeafContainingVertex(p);
 394     vector sizeVec;
 395     if( cLabel < 0 )
 396     {
 397         sizeVec.x() = sizeVec.y() = sizeVec.z() = searchRange_;
 398     }
 399     else
 400     {
 401         const scalar s = 0.75 * leaves_[cLabel]->size(rootBox_);
 402         sizeVec.x() = sizeVec.y() = sizeVec.z() = s;
 403     }
 404
 405     //- find nearest surface vertex to the point p
 406     bool found(false);
 407     label iterationI(0);
 408     DynList<const meshOctreeCube*, 256> neighbours;
 409
 410     const pointField& points = surface_.points();
 411     const VRWGraph& pEdges = surface_.pointEdges();
 412     const VRWGraph& eFacets = surface_.edgeFacets();
 413
 414     distSq = VGREAT;
 415     nearestPoint = -1;
 416
 417     do
 418     {
 419         boundBox bb(p - sizeVec, p + sizeVec);
 420
 421         neighbours.clear();
 422         findLeavesContainedInBox(bb, neighbours);
 423         labelHashSet checkedPoint;
 424
 425         //- find nearest projection
 426         forAll(neighbours, neiI)
 427         {
 428             if( !neighbours[neiI]->hasContainedElements() )
 429                 continue;
 430
 431             const VRWGraph& ct =
 432                 neighbours[neiI]->slotPtr()->containedTriangles_;
 433             const constRow el = ct[neighbours[neiI]->containedElements()];
 434             forAll(el, tI)
 435             {
 436                 const labelledTri& tri = surface_[el[tI]];
 437
 438                 forAll(tri, pI)
 439                 {
 440                     const label spI = tri[pI];
 441
 442                     if( checkedPoint.found(spI) )
 443                         continue;
 444
 445                     checkedPoint.insert(spI);
 446
 447                     DynList<label> nodePatches;
 448                     label nEdges(0);
 449
 450                     forAllRow(pEdges, spI, i)
 451                     {
 452                         const label eI = pEdges(spI, i);
 453
 454                         if( eFacets.sizeOfRow(eI) != 2 )
 455                             break;
 456
 457                         if(
 458                             surface_[eFacets(eI, 0)].region() !=
 459                             surface_[eFacets(eI, 1)].region()
 460                         )
 461                         {
 462                             //- found an edge attached to this vertex
 463                             ++nEdges;
 464                             nodePatches.appendIfNotIn
 465                             (
 466                                 surface_[eFacets(eI, 0)].region()
 467                             );
 468                             nodePatches.appendIfNotIn
 469                             (
 470                                 surface_[eFacets(eI, 1)].region()
 471                             );
 472                         }
 473                     }
 474
 475                     if( nEdges > 2 )
 476                     {
 477                         //- check if all required patches
 478                         //- are present at this corner
 479                         nEdges = 0;
 480                         forAll(patches, i)
 481                         {
 482                             if( nodePatches.contains(patches[i]) )
 483                                 ++nEdges;
 484                         }
 485
 486                         if( nEdges >= patches.size() )
 487                         {
 488                             //- all patches are present, check the distance
 489                             const scalar dSq = Foam::magSqr(points[spI] - p);
 490
 491                             if( dSq < distSq )
 492                             {
 493                                 distSq = dSq;
 494                                 found = true;
 495                                 nearest = points[spI];
 496                                 nearestPoint = spI;
 497                             }
 498                         }
 499                     }
 500                 }
 501             }
 502         }
 503
 504         if( !found )
 505             sizeVec *= 2.0;
 506
 507     } while( !found && (iterationI++ < 5) );
 508
 509     return found;
 510 }
 511
 512 bool meshOctree::findNearestPointToPatches
 513 (
 514     point& nearest,
 515     scalar& distSq,
 516     const point& p,
 517     const DynList<label>& patches,
 518     const scalar tol
 519 ) const
 520 {
 521     if( patches.size() == 0 )
 522         return false;
 523
 524     point mapPoint(p);
 525     scalar dSq(VGREAT);
 526
 527     bool found(false);
 528     if( patches.size() == 2 )
 529     {
 530         label nse;
 531         found = findNearestEdgePoint(mapPoint, dSq, nse, p, patches);
 532     }
 533     else if( patches.size() > 2 )
 534     {
 535         label nsp;
 536         found = findNearestCorner(mapPoint, dSq, nsp, p, patches);
 537     }
 538
 539     point mapPointApprox(p);
 540     scalar distSqApprox;
 541     label iter(0);
 542     while( iter++ < 20 )
 543     {
 544         point newP(vector::zero);
 545         forAll(patches, patchI)
 546         {
 547             point np;
 548             label nearestTri;
 549             this->findNearestSurfacePointInRegion
 550             (
 551                 np,
 552                 distSqApprox,
 553                 nearestTri,
 554                 patches[patchI],
 555                 mapPointApprox
 556             );
 557
 558             newP += np;
 559         }
 560
 561         newP /= patches.size();
 562         if( Foam::magSqr(newP - mapPointApprox) < tol * dSq )
 563             break;
 564
 565         mapPointApprox = newP;
 566     }
 567
 568     distSq = Foam::magSqr(mapPointApprox - p);
 569     nearest = mapPointApprox;
 570
 571     if( found && (dSq < 1.5 * distSq) )
 572     {
 573         nearest = mapPoint;
 574         distSq = dSq;
 575     }
 576
 577     return true;
 578 }
 579
 580 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 581
 582 } // End namespace Foam
 583
 584 // ************************************************************************* //