src/mesh/cfMesh/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         const 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++ < 100) );
 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& sPoints = 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 //    Info << "Finding nearest point for " << p << " size vec " << sizeVec << endl;
 240
 241     do
 242     {
 243         const boundBox bb(p - sizeVec, p + sizeVec);
 244
 245         neighbours.clear();
 246         findLeavesContainedInBox(bb, neighbours);
 247
 248 //        Info << "Iteration " << iterationI << " nu found boxes "
 249 //             << neighbours.size() << endl;
 250
 251         forAll(neighbours, neiI)
 252         {
 253             if( !neighbours[neiI]->hasContainedEdges() )
 254                 continue;
 255
 256             const VRWGraph& containedEdges =
 257                 neighbours[neiI]->slotPtr()->containedEdges_;
 258             const constRow ce =
 259                 containedEdges[neighbours[neiI]->containedEdges()];
 260
 261 //            Info << "Number of contained edges in box "
 262 //                 << neighbours[neiI]->cubeLabel()
 263 //                 << " are " << ce.size() << endl;
 264
 265             forAll(ce, eI)
 266             {
 267                 const label edgeI = ce[eI];
 268
 269                 //- find if the edge is in correct patches
 270                 bool correctPatches(true);
 271
 272                 forAllRow(edgeFaces, edgeI, efI)
 273                 {
 274                     const label facetI = edgeFaces(edgeI, efI);
 275
 276                     if( !regions.contains(surface_[facetI].region()) )
 277                     {
 278                         correctPatches = false;
 279                         break;
 280                     }
 281                 }
 282
 283                 if( !correctPatches )
 284                     continue;
 285
 286                 const edge& e = edges[edgeI];
 287                 const point& sp = sPoints[e.start()];
 288                 const point& ep = sPoints[e.end()];
 289                 const point np = help::nearestPointOnTheEdgeExact(sp, ep, p);
 290                 const scalar dSq = Foam::magSqr(np - p);
 291
 292                 if( dSq < distSq )
 293                 {
 294                     distSq = dSq;
 295                     edgePoint = np;
 296                     nearestEdge = ce[eI];
 297                     foundAnEdge = true;
 298                 }
 299             }
 300         }
 301
 302         if( !foundAnEdge )
 303             sizeVec *= 2.0;
 304
 305     } while( !foundAnEdge && (++iterationI < 3) );
 306
 307     return foundAnEdge;
 308 }
 309
 310 bool meshOctree::findNearestPointToEdge
 311 (
 312     point& nearest,
 313     scalar& distSq,
 314     label& nearestEdge,
 315     const FixedList<point, 2>& edgePoints,
 316     const FixedList<label, 2>& edgePointRegions
 317 ) const
 318 {
 319     const point c = 0.5 * (edgePoints[0] + edgePoints[1]);
 320     const scalar dst = mag(edgePoints[0] - edgePoints[1]);
 321     vector sizeVec(dst, dst, dst);
 322
 323     boundBox bb(c - 0.75 * sizeVec, c + 0.75 * sizeVec);
 324
 325     DynList<const meshOctreeCube*, 256> leavesInBox;
 326     findLeavesContainedInBox(bb, leavesInBox);
 327
 328     const VRWGraph& edgeFaces = surface_.edgeFacets();
 329     const pointField& points = surface_.points();
 330     const edgeLongList& surfaceEdges = surface_.edges();
 331
 332     distSq = VGREAT;
 333     nearestEdge = -1;
 334
 335     bool found(false);
 336
 337     forAll(leavesInBox, leafI)
 338     {
 339         if( !leavesInBox[leafI]->hasContainedEdges() )
 340             continue;
 341
 342         const VRWGraph& containedEdges =
 343             leavesInBox[leafI]->slotPtr()->containedEdges_;
 344         const constRow edges =
 345             containedEdges[leavesInBox[leafI]->containedEdges()];
 346
 347         forAll(edges, eI)
 348         {
 349             const constRow ef = edgeFaces[edges[eI]];
 350             if( ef.size() != 2 )
 351                 continue;
 352
 353             if(
 354                 (
 355                     (surface_[ef[0]].region() == edgePointRegions[0]) &&
 356                     (surface_[ef[1]].region() == edgePointRegions[1])
 357                 ) ||
 358                 (
 359                     (surface_[ef[1]].region() == edgePointRegions[0]) &&
 360                     (surface_[ef[0]].region() == edgePointRegions[1])
 361                 )
 362             )
 363             {
 364                 const edge& edg = surfaceEdges[edges[eI]];
 365
 366                 point nearestOnEdge, nearestOnLine;
 367                 if(
 368                     help::nearestEdgePointToTheLine
 369                     (
 370                         points[edg[0]],
 371                         points[edg[1]],
 372                         edgePoints[0],
 373                         edgePoints[1],
 374                         nearestOnEdge,
 375                         nearestOnLine
 376                     )
 377                 )
 378                 {
 379                     if( magSqr(nearestOnEdge - nearestOnLine) < distSq )
 380                     {
 381                         nearest = nearestOnEdge;
 382                         nearestEdge = edges[eI];
 383                         distSq = magSqr(nearestOnEdge - nearestOnLine);
 384                         found = true;
 385                     }
 386                 }
 387             }
 388         }
 389     }
 390
 391     return found;
 392 }
 393
 394 bool meshOctree::findNearestCorner
 395 (
 396     point& nearest,
 397     scalar& distSq,
 398     label& nearestPoint,
 399     const point& p,
 400     const DynList<label>& patches
 401 ) const
 402 {
 403
 404
 405     const label cLabel = findLeafContainingVertex(p);
 406     vector sizeVec;
 407     if( cLabel < 0 )
 408     {
 409         sizeVec.x() = sizeVec.y() = sizeVec.z() = searchRange_;
 410     }
 411     else
 412     {
 413         const scalar s = 0.75 * leaves_[cLabel]->size(rootBox_);
 414         sizeVec.x() = sizeVec.y() = sizeVec.z() = s;
 415     }
 416
 417     //- find nearest surface vertex to the point p
 418     bool found(false);
 419     label iterationI(0);
 420     DynList<const meshOctreeCube*, 256> neighbours;
 421
 422     const pointField& points = surface_.points();
 423     const VRWGraph& pEdges = surface_.pointEdges();
 424     const VRWGraph& eFacets = surface_.edgeFacets();
 425
 426     distSq = VGREAT;
 427     nearestPoint = -1;
 428
 429     do
 430     {
 431         const boundBox bb(p - sizeVec, p + sizeVec);
 432
 433         neighbours.clear();
 434         findLeavesContainedInBox(bb, neighbours);
 435         labelHashSet checkedPoint;
 436
 437         //- find nearest projection
 438         forAll(neighbours, neiI)
 439         {
 440             if( !neighbours[neiI]->hasContainedElements() )
 441                 continue;
 442
 443             const VRWGraph& ct =
 444                 neighbours[neiI]->slotPtr()->containedTriangles_;
 445             const constRow el = ct[neighbours[neiI]->containedElements()];
 446             forAll(el, tI)
 447             {
 448                 const labelledTri& tri = surface_[el[tI]];
 449
 450                 forAll(tri, pI)
 451                 {
 452                     const label spI = tri[pI];
 453
 454                     if( checkedPoint.found(spI) )
 455                         continue;
 456
 457                     checkedPoint.insert(spI);
 458
 459                     DynList<label> nodePatches;
 460                     label nEdges(0);
 461
 462                     forAllRow(pEdges, spI, i)
 463                     {
 464                         const label eI = pEdges(spI, i);
 465
 466                         if( eFacets.sizeOfRow(eI) != 2 )
 467                             break;
 468
 469                         if(
 470                             surface_[eFacets(eI, 0)].region() !=
 471                             surface_[eFacets(eI, 1)].region()
 472                         )
 473                         {
 474                             //- found an edge attached to this vertex
 475                             ++nEdges;
 476                             nodePatches.appendIfNotIn
 477                             (
 478                                 surface_[eFacets(eI, 0)].region()
 479                             );
 480                             nodePatches.appendIfNotIn
 481                             (
 482                                 surface_[eFacets(eI, 1)].region()
 483                             );
 484                         }
 485                     }
 486
 487                     if( nEdges > 2 )
 488                     {
 489                         //- check if all required patches
 490                         //- are present at this corner
 491                         nEdges = 0;
 492                         forAll(patches, i)
 493                         {
 494                             if( nodePatches.contains(patches[i]) )
 495                                 ++nEdges;
 496                         }
 497
 498                         if( nEdges >= patches.size() )
 499                         {
 500                             //- all patches are present, check the distance
 501                             const scalar dSq = Foam::magSqr(points[spI] - p);
 502
 503                             if( dSq < distSq )
 504                             {
 505                                 distSq = dSq;
 506                                 found = true;
 507                                 nearest = points[spI];
 508                                 nearestPoint = spI;
 509                             }
 510                         }
 511                     }
 512                 }
 513             }
 514         }
 515
 516         if( !found )
 517             sizeVec *= 2.0;
 518
 519     } while( !found && (iterationI++ < 3) );
 520
 521     return found;
 522 }
 523
 524 bool meshOctree::findNearestPointToPatches
 525 (
 526     point& nearest,
 527     scalar& distSq,
 528     const point& p,
 529     const DynList<label>& patches,
 530     const scalar tol
 531 ) const
 532 {
 533     if( patches.size() == 0 )
 534         return false;
 535
 536     nearest = p;
 537     scalar distSqApprox;
 538     label iter(0);
 539     while( iter++ < 40 )
 540     {
 541         point newP(vector::zero);
 542         forAll(patches, patchI)
 543         {
 544             point np;
 545             label nearestTri;
 546             this->findNearestSurfacePointInRegion
 547             (
 548                 np,
 549                 distSqApprox,
 550                 nearestTri,
 551                 patches[patchI],
 552                 nearest
 553             );
 554
 555             newP += np;
 556         }
 557
 558         newP /= patches.size();
 559         distSq = magSqr(newP - p);
 560
 561         if( Foam::magSqr(newP - nearest) < tol * distSq )
 562             break;
 563
 564         nearest = newP;
 565     }
 566
 567     return true;
 568 }
 569
 570 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 571
 572 } // End namespace Foam
 573
 574 // ************************************************************************* //