src/meshLibrary/utilities/helperFunctions/helperFunctionsGeometryQueriesI.H

   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 "error.H"
  29 #include "helperFunctionsGeometryQueries.H"
  30 #include "helperFunctionsTopologyManipulation.H"
  31 #include "edgeList.H"
  32 #include "pointField.H"
  33 #include "boolList.H"
  34 #include "triSurf.H"
  35 #include "matrix3D.H"
  36 #include "HashSet.H"
  37 #include "tetrahedron.H"
  38 #include "boundBox.H"
  39 #include "Map.H"
  40
  41 //#define DEBUG_pMesh
  42
  43 namespace Foam
  44 {
  45
  46 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
  47
  48 namespace help
  49 {
  50
  51 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
  52
  53 template<class ListType>
  54 inline bool isnan(const ListType& l)
  55 {
  56     forAll(l, i)
  57         if( l[i] != l[i] )
  58             return true;
  59
  60     return false;
  61 }
  62
  63 template<class ListType>
  64 bool isinf(const ListType& l)
  65 {
  66     forAll(l, i)
  67         if( (l[i] < -VGREAT) || (l[i] > VGREAT) )
  68             return true;
  69
  70     return false;
  71 }
  72
  73 inline bool isSharedEdgeConvex
  74 (
  75     const pointField& points,
  76     const face& f1,
  77     const face& f2
  78 )
  79 {
  80     face triOwn(3);
  81     face triNei(3);
  82
  83     forAll(f1, pI)
  84     {
  85         label pos = f2.which(f1[pI]);
  86         if( pos == -1 )
  87             continue;
  88
  89         triNei[0] = f2[pos];
  90         triNei[1] = f2.nextLabel(pos);
  91         triNei[2] = f2.prevLabel(pos);
  92
  93         triOwn[0] = f1[pI];
  94         triOwn[1] = f1.nextLabel(pI);
  95         triOwn[2] = f1.prevLabel(pI);
  96
  97         scalar vol(0.0);
  98
  99         forAll(triOwn, pJ)
 100         {
 101             if( triNei.which(triOwn[pJ]) == -1 )
 102             {
 103                 tetrahedron<point, point> tet
 104                 (
 105                     points[triNei[0]],
 106                     points[triNei[1]],
 107                     points[triNei[2]],
 108                     points[triOwn[pJ]]
 109                 );
 110
 111                 vol = tet.mag();
 112                 break;
 113             }
 114         }
 115
 116         if( vol > -VSMALL )
 117             return false;
 118     }
 119
 120     return true;
 121 }
 122
 123 inline faceList mergePatchFaces
 124 (
 125     const List< DynList<label> >& pfcs,
 126     const pointField& polyPoints
 127 )
 128 {
 129     //- merge faces which share a common edge
 130     faceList patchFaces(pfcs.size());
 131     label counter(0);
 132     forAll(pfcs, faceI)
 133         if( pfcs[faceI].size() > 2 )
 134         {
 135             const DynList<label>& f = pfcs[faceI];
 136             face f_(f.size());
 137             forAll(f_, fJ)
 138                 f_[fJ] = f[fJ];
 139
 140             patchFaces[counter++] = f_;
 141         }
 142
 143     patchFaces.setSize(counter);
 144
 145     bool merged;
 146     do
 147     {
 148         faceList mergedFaces(patchFaces.size());
 149         boolList currentlyMerged(patchFaces.size(), false);
 150
 151         counter = 0;
 152         merged = false;
 153
 154         for(label nI=0;nI<(patchFaces.size()-1);nI++)
 155         {
 156             vector n0 = patchFaces[nI].normal(polyPoints);
 157             n0 /= mag(n0);
 158
 159             for(label nJ=nI+1;nJ<patchFaces.size();nJ++)
 160             {
 161                 vector n1 = patchFaces[nI].normal(polyPoints);
 162                 n1 /= mag(n1);
 163                 if(
 164                     help::shareAnEdge(patchFaces[nI], patchFaces[nJ]) &&
 165                     mag(n0 & n1) > 0.95
 166                 )
 167                 {
 168                     merged = true;
 169                     currentlyMerged[nI] = currentlyMerged[nJ] = true;
 170                     mergedFaces[counter++] =
 171                         help::mergeTwoFaces
 172                         (
 173                             patchFaces[nI],
 174                             patchFaces[nJ]
 175                         );
 176
 177                     break;
 178                 }
 179             }
 180
 181             if( merged ) break;
 182         }
 183
 184         forAll(patchFaces, pfI)
 185             if( !currentlyMerged[pfI] )
 186                 mergedFaces.newElmt(counter++) = patchFaces[pfI];
 187
 188         if( merged )
 189         {
 190             patchFaces.setSize(counter);
 191             for(label k=0;k<counter;k++)
 192                 patchFaces[k] = mergedFaces[k];
 193         }
 194
 195     } while( merged );
 196
 197     return patchFaces;
 198 }
 199
 200 inline bool vertexOnLine(const point& p, const edge& e, const pointField& ep)
 201 {
 202     vector v = e.vec(ep);
 203     v /= mag(v);
 204
 205     vector pv = p - ep[e.start()];
 206     pv /= mag(pv);
 207
 208     if( mag(pv & v) > (1.0-SMALL) )
 209         return true;
 210
 211     return false;
 212 }
 213
 214 inline bool vertexInPlane(const point& p, const plane& pl)
 215 {
 216     const vector& n = pl.normal();
 217     const point& fp = pl.refPoint();
 218
 219     vector d = p - fp;
 220     if( mag(d) > VSMALL )
 221         d /= mag(d);
 222
 223     if( mag(d & n) < SMALL )
 224         return true;
 225
 226     return false;
 227 }
 228
 229 inline bool planeIntersectsEdge
 230 (
 231     const point& start,
 232     const point& end,
 233     const plane& pl,
 234     point& intersection
 235 )
 236 {
 237     const vector v = end - start;
 238
 239     const vector& n = pl.normal();
 240     const point& fp = pl.refPoint();
 241
 242     if( (n & (v/mag(v))) < SMALL )
 243         return false;
 244
 245     const scalar t((n & (fp - start)) / (n & v));
 246
 247     if( (t > -SMALL) && (t < (1.0+SMALL)) )
 248     {
 249         intersection = start + v * t;
 250         return true;
 251     }
 252
 253     return false;
 254 }
 255
 256 inline bool pointInTetrahedron
 257 (
 258     const point& p,
 259     const tetrahedron<point, point>& tet
 260 )
 261 {
 262     const vector v0 = tet.a() - tet.d();
 263     const vector v1 = tet.b() - tet.d();
 264     const vector v2 = tet.c() - tet.d();
 265     const vector sp = p - tet.d();
 266
 267     matrix3D mat;
 268     FixedList<scalar, 3> source;
 269     for(label i=0;i<3;++i)
 270     {
 271         mat[i][0] = v0[i];
 272         mat[i][1] = v1[i];
 273         mat[i][2] = v2[i];
 274         source[i] = sp[i];
 275     }
 276
 277     //- check the determinant of the transformation
 278     const scalar det = mat.determinant();
 279
 280     if( mag(det) < VSMALL )
 281         return false;
 282
 283     //- get the coordinates of the point in the barycentric corrdinate system
 284     const scalar u0 = mat.solveFirst(source);
 285
 286     if( (u0 < -SMALL) || (u0 > (1.0+SMALL)) )
 287         return false;
 288
 289     const scalar u1 = mat.solveSecond(source);
 290
 291     if( (u1 < -SMALL) || ((u0+u1) > (1.0+SMALL)) )
 292         return false;
 293
 294     const scalar u2 = mat.solveThird(source);
 295
 296     if( (u2 < -SMALL) || (u2 > (1.0+SMALL)) )
 297         return false;
 298
 299
 300     const scalar u3 = 1.0 - u0 - u1 - u2;
 301
 302     if( (u3 < -SMALL) || (u3 > (1.0+SMALL)) )
 303         return false;
 304
 305     return true;
 306 }
 307
 308 inline bool nearestEdgePointToTheLine
 309 (
 310     const point& edgePoint0,
 311     const point& edgePoint1,
 312     const point& lp0,
 313     const point& lp1,
 314     point& nearestOnEdge,
 315     point& nearestOnLine
 316 )
 317 {
 318     const vector v = lp1 - lp0;
 319     const vector d = lp0 - edgePoint0;
 320     const vector e = edgePoint1 - edgePoint0;
 321
 322     const scalar vMag = mag(v);
 323     if( vMag < VSMALL )
 324         return false;
 325
 326     const scalar eMag = mag(e);
 327     if( eMag < VSMALL )
 328     {
 329         nearestOnEdge = edgePoint0;
 330         nearestOnLine = nearestPointOnTheEdge(lp0, lp1, nearestOnEdge);
 331         return true;
 332     }
 333
 334     if( mag((v/vMag) & (e/eMag)) > (1.0 - SMALL) )
 335         return false;
 336
 337     tensor mat(tensor::zero);
 338     mat.xx() = (v&v);
 339     mat.xy() = mat.yx() = -1.0 * (v&e);
 340     mat.yy() = (e&e);
 341     mat.zz() = SMALL;
 342
 343     vector source(vector::zero);
 344     source[0] = -1.0 * (d&v);
 345     source[1] = (d&e);
 346
 347     const vector sol = (inv(mat) & source);
 348
 349     nearestOnLine = lp0 + v * sol[0];
 350     if( sol[1] > 1.0 )
 351     {
 352         nearestOnEdge = edgePoint1;
 353     }
 354     else if( sol[1] < 0.0 )
 355     {
 356         nearestOnEdge = edgePoint0;
 357     }
 358     else
 359     {
 360         nearestOnEdge = edgePoint0 + e * sol[1];
 361     }
 362
 363     return true;
 364 }
 365
 366 inline point nearestPointOnTheTriangle
 367 (
 368     const label tI,
 369     const triSurf& surface,
 370     const point& p
 371 )
 372 {
 373     const labelledTri& ltri = surface[tI];
 374     const pointField& points = surface.points();
 375
 376     const vector v0 = points[ltri[1]] - points[ltri[0]];
 377     const vector v1 = points[ltri[2]] - points[ltri[0]];
 378     const vector v2 = p - points[ltri[0]];
 379
 380     const scalar dot00 = (v0 & v0);
 381     const scalar dot01 = (v0 & v1);
 382     const scalar dot02 = (v0 & v2);
 383     const scalar dot11 = (v1 & v1);
 384     const scalar dot12 = (v1 & v2);
 385
 386     // Compute barycentric coordinates
 387     const scalar det = dot00 * dot11 - dot01 * dot01;
 388
 389     if( mag(det) < VSMALL )
 390     {
 391         point nearest(p);
 392         scalar dist(VGREAT);
 393         const edgeList edges = ltri.edges();
 394         forAll(edges, eI)
 395         {
 396             const edge& e = edges[eI];
 397             const point np =
 398                 nearestPointOnTheEdge
 399                 (
 400                     points[e.start()],
 401                     points[e.end()],
 402                     p
 403                 );
 404
 405             if( magSqr(p - np) < dist )
 406             {
 407                 nearest = np;
 408                 dist = magSqr(p - np);
 409             }
 410         }
 411
 412         return nearest;
 413     }
 414
 415     const scalar u = (dot11 * dot02 - dot01 * dot12) / det;
 416     const scalar v = (dot00 * dot12 - dot01 * dot02) / det;
 417
 418     const point pProj = points[ltri[0]] + u * v0 + v * v1;
 419
 420     //- Check if point is in triangle
 421     if( (u >= -SMALL) && (v >= -SMALL) && ((u + v) <= (1.0+SMALL)) )
 422     {
 423         return pProj;
 424     }
 425     else
 426     {
 427         if( u < -SMALL )
 428         {
 429             const edge e(ltri[0], ltri[2]);
 430             const scalar ed = ((pProj - points[e.start()]) & v1) / magSqr(v1);
 431
 432             if( ed > 1.0 )
 433             {
 434                 return points[e.end()];
 435             }
 436             else if( ed < 0.0 )
 437             {
 438                 return points[e.start()];
 439             }
 440             else
 441             {
 442                 return (points[e.start()] + v1 * ed);
 443             }
 444         }
 445         else if( v < -SMALL )
 446         {
 447             const edge e(ltri[0], ltri[1]);
 448             const scalar ed = ((pProj - points[e.start()]) & v0) / magSqr(v0);
 449
 450             if( ed > 1.0 )
 451             {
 452                 return points[e.end()];
 453             }
 454             else if( ed < 0.0 )
 455             {
 456                 return points[e.start()];
 457             }
 458             else
 459             {
 460                 return (points[e.start()] + v0 * ed);
 461             }
 462         }
 463         else
 464         {
 465             const edge e(ltri[2], ltri[1]);
 466             const vector ev = e.vec(points);
 467             const scalar ed = ((pProj - points[e.start()]) & ev) / magSqr(ev);
 468
 469             if( ed > 1.0 )
 470             {
 471                 return points[e.end()];
 472             }
 473             else if( ed < 0.0 )
 474             {
 475                 return points[e.start()];
 476             }
 477             else
 478             {
 479                 return (points[e.start()] + ev * ed);
 480             }
 481         }
 482     }
 483
 484     return pProj;
 485 }
 486
 487 inline bool triLineIntersection
 488 (
 489     const triangle<point, point>& tria,
 490     const point& lineStart,
 491     const point& lineEnd,
 492     point& intersection
 493 )
 494 {
 495     const point& p0 = tria.a();
 496     const vector v(lineStart - lineEnd);
 497     const vector v0 = tria.b() - p0;
 498     const vector v1 = tria.c() - p0;
 499     const vector sp = lineStart - p0;
 500
 501     matrix3D mat;
 502     FixedList<scalar, 3> source;
 503     for(label i=0;i<3;++i)
 504     {
 505         mat[i][0] = v0[i];
 506         mat[i][1] = v1[i];
 507         mat[i][2] = v[i];
 508         source[i] = sp[i];
 509     }
 510
 511     const scalar det = mat.determinant();
 512
 513     if( mag(det) < SMALL )
 514         return false;
 515
 516     const scalar t = mat.solveThird(source);
 517
 518     if( (t < -SMALL) || (t > (1.0+SMALL)) )
 519         return false;
 520
 521     const scalar u0 = mat.solveFirst(source);
 522
 523     if( u0 < -SMALL )
 524         return false;
 525
 526     const scalar u1 = mat.solveSecond(source);
 527
 528     if( (u1 < -SMALL) || ((u0+u1) > (1.0+SMALL)) )
 529         return false;
 530
 531     intersection = lineStart - t * v;
 532     return true;
 533 }
 534
 535 inline bool triLineIntersection
 536 (
 537     const triSurf& surface,
 538     const label tI,
 539     const point& s,
 540     const point& e,
 541     point& intersection
 542 )
 543 {
 544     const pointField& pts = surface.points();
 545     const labelledTri& tri = surface[tI];
 546
 547     const triangle<point, point> tria
 548     (
 549         pts[tri[0]],
 550         pts[tri[1]],
 551         pts[tri[2]]
 552     );
 553
 554     return triLineIntersection(tria, s, e, intersection);
 555 }
 556
 557 inline bool boundBoxLineIntersection
 558 (
 559     const point& s,
 560     const point& e,
 561     const boundBox& bb
 562 )
 563 {
 564     scalar tMax(1.0+SMALL), tMin(-SMALL);
 565
 566     const vector v = e - s;
 567     const scalar d = mag(v);
 568
 569     //- check if the vector has length
 570     if( d < VSMALL )
 571     {
 572         if( bb.contains(s) )
 573             return true;
 574
 575         return false;
 576     }
 577
 578     const point& pMin = bb.min();
 579     const point& pMax = bb.max();
 580
 581     //- check coordinates
 582     for(label dir=0;dir<3;++dir)
 583     {
 584         const scalar vd = v[dir];
 585         const scalar sd = s[dir];
 586
 587         if( mag(vd) > (SMALL * d) )
 588         {
 589             if( vd >= 0.0 )
 590             {
 591                 tMin = Foam::max(tMin, (pMin[dir] - sd) / vd);
 592                 tMax = Foam::min(tMax, (pMax[dir] - sd) / vd);
 593             }
 594             else
 595             {
 596                 tMin = Foam::max(tMin, (pMax[dir] - sd) / vd);
 597                 tMax = Foam::min(tMax, (pMin[dir] - sd) / vd);
 598             }
 599         }
 600         else if( (sd < pMin[dir]) || (sd > pMax[dir]) )
 601         {
 602             return false;
 603         }
 604     }
 605
 606     if( (tMax - tMin) > -SMALL )
 607         return true;
 608
 609     return false;
 610 }
 611
 612 inline bool doFaceAndTriangleIntersect
 613 (
 614     const triSurf& surface,
 615     const label triI,
 616     const face& f,
 617     const pointField& facePoints
 618 )
 619 {
 620     const pointField& triPoints = surface.points();
 621
 622     const point centre = f.centre(facePoints);
 623     point intersection;
 624
 625     //- check if any triangle edge intersects the face
 626     const labelledTri& tri = surface[triI];
 627
 628     forAll(tri, eI)
 629     {
 630         const point& s = triPoints[tri[eI]];
 631         const point& e = triPoints[tri[(eI+1)%3]];
 632
 633         forAll(f, pI)
 634         {
 635             const triangle<point, point> tria
 636             (
 637                 facePoints[f[pI]],
 638                 facePoints[f.nextLabel(pI)],
 639                 centre
 640             );
 641
 642             const bool currIntersection =
 643                 help::triLineIntersection
 644                 (
 645                     tria,
 646                     s,
 647                     e,
 648                     intersection
 649                 );
 650
 651             if( currIntersection )
 652                 return true;
 653         }
 654     }
 655
 656     //- check if any face edges intersect the triangle
 657     forAll(f, pI)
 658     {
 659         const point& s = facePoints[f[pI]];
 660         const point& e = facePoints[f.nextLabel(pI)];
 661
 662         const bool intersected =
 663             help::triLineIntersection
 664             (
 665                 surface,
 666                 triI,
 667                 s,
 668                 e,
 669                 intersection
 670             );
 671
 672         if( intersected )
 673             return true;
 674     }
 675
 676     return false;
 677 }
 678
 679 inline bool pointInsideFace
 680 (
 681     const point& p,
 682     const face& f,
 683     const vector& n,
 684     const pointField& fp
 685 )
 686 {
 687     const edgeList fe = f.edges();
 688     forAll(f, pI)
 689     {
 690         if( mag(p - fp[f[pI]]) < SMALL )
 691             return true;
 692
 693         vector pv = p - fp[f[pI]];
 694         pv /= mag(pv);
 695
 696         vector lv = n ^ fe[pI].vec(fp);
 697         lv /= mag(lv);
 698
 699         const scalar d = pv & lv;
 700         if( d < -SMALL )
 701         {
 702             return false;
 703         }
 704     }
 705
 706     return true;
 707 }
 708
 709 inline bool pointInsideFace
 710 (
 711     const point& p,
 712     const face& f,
 713     const pointField& fp
 714 )
 715 {
 716     vector n = f.normal(fp);
 717     n /= mag(n);
 718
 719     const edgeList fe = f.edges();
 720     forAll(f, pI)
 721     {
 722         if( mag(p - fp[f[pI]]) < SMALL )
 723         return true;
 724
 725         vector pv = p - fp[f[pI]];
 726         pv /= mag(pv);
 727
 728         vector lv = n ^ fe[pI].vec(fp);
 729         lv /= mag(lv);
 730
 731         const scalar d = pv & lv;
 732         if( d < -SMALL )
 733         {
 734             return false;
 735         }
 736     }
 737
 738     return true;
 739 }
 740
 741 inline point nearestPointOnTheEdge
 742 (
 743     const point& edgePoint0,
 744     const point& edgePoint1,
 745     const point& p
 746 )
 747 {
 748     const vector e = edgePoint1 - edgePoint0;
 749     const scalar d = mag(e);
 750     const vector k = p - edgePoint0;
 751
 752     if( d < VSMALL )
 753         return edgePoint0;
 754
 755     return edgePoint0 + ((e / (d*d)) * (e & k));
 756 }
 757
 758 inline point nearestPointOnTheEdgeExact
 759 (
 760     const point& edgePoint0,
 761     const point& edgePoint1,
 762     const point& p
 763 )
 764 {
 765     const vector e = edgePoint1 - edgePoint0;
 766     const scalar d = mag(e);
 767     const vector k = p - edgePoint0;
 768
 769     if( d < VSMALL )
 770         return edgePoint0;
 771
 772     const scalar t = (e & k) / (d * d);
 773     if( t > 1.0 )
 774     {
 775         return edgePoint1;
 776     }
 777     else if( t < 0.0 )
 778     {
 779         return edgePoint0;
 780     }
 781
 782     return edgePoint0 + (e * t);
 783 }
 784
 785 inline scalar distanceOfPointFromTheEdge
 786 (
 787     const point& edgePoint0,
 788     const point& edgePoint1,
 789     const point& p
 790 )
 791 {
 792     return mag(nearestPointOnTheEdge(edgePoint0, edgePoint1, p) - p);
 793 }
 794
 795 inline label numberOfFaceGroups
 796 (
 797     const labelHashSet& containedElements,
 798     const point& centre,
 799     const scalar range,
 800     const triSurf& surface
 801 )
 802 {
 803     const pointField& points = surface.points();
 804     const edgeLongList& edges = surface.edges();
 805     const VRWGraph& faceEdges = surface.facetEdges();
 806     const VRWGraph& edgeFaces = surface.edgeFacets();
 807
 808     labelHashSet triaInRange(containedElements.size());
 809     const scalar rangeSq = range * range;
 810     forAllConstIter(labelHashSet, containedElements, it)
 811     {
 812         const point p = nearestPointOnTheTriangle(it.key(), surface, centre);
 813         if( magSqr(p - centre) < rangeSq )
 814             triaInRange.insert(it.key());
 815     }
 816
 817     Map<label> elGroup(triaInRange.size());
 818     Map<label> testEdge(triaInRange.size());
 819
 820     label nGroups(0);
 821
 822     DynList<label> front;
 823
 824     forAllConstIter(labelHashSet, triaInRange, it)
 825         if( !elGroup.found(it.key()) )
 826         {
 827             front.clear();
 828             front.append(it.key());
 829             elGroup.insert(it.key(), nGroups);
 830
 831             while( front.size() )
 832             {
 833                 const label fLabel = front.removeLastElement();
 834
 835                 forAllRow(faceEdges, fLabel, feI)
 836                 {
 837                     const label edgeI = faceEdges(fLabel, feI);
 838
 839                     //- check if the edge intersects the bounding box
 840                     if( testEdge.found(edgeI) )
 841                     {
 842                         if( !testEdge[edgeI] )
 843                             continue;
 844                     }
 845                     else
 846                     {
 847                         const point& s = points[edges[edgeI][0]];
 848                         const point& e = points[edges[edgeI][1]];
 849                         const point np = nearestPointOnTheEdge(s, e, centre);
 850                         if( magSqr(np - centre) < rangeSq )
 851                         {
 852                             testEdge.insert(edgeI, 1);
 853                         }
 854                         else
 855                         {
 856                             testEdge.insert(edgeI, 0);
 857                             continue;
 858                         }
 859                     }
 860
 861                     forAllRow(edgeFaces, edgeI, efI)
 862                     {
 863                         const label nei = edgeFaces(edgeI, efI);
 864                         if
 865                         (
 866                             triaInRange.found(nei) &&
 867                             !elGroup.found(nei)
 868                         )
 869                         {
 870                             elGroup.insert(nei, nGroups);
 871                             front.append(nei);
 872                         }
 873                     }
 874                 }
 875             }
 876
 877             ++nGroups;
 878         }
 879
 880     return nGroups;
 881 }
 882
 883 inline label numberOfEdgeGroups
 884 (
 885     const labelHashSet& containedEdges,
 886     const point& centre,
 887     const scalar range,
 888     const triSurf& surface
 889 )
 890 {
 891     const pointField& points = surface.points();
 892     const edgeLongList& edges = surface.edges();
 893     const VRWGraph& pointEdges = surface.pointEdges();
 894
 895     const scalar rangeSq = range * range;
 896     labelHashSet edgesInRange(containedEdges.size());
 897     forAllConstIter(labelHashSet, containedEdges, it)
 898     {
 899         const edge& e = edges[it.key()];
 900         const point& sp = points[e[0]];
 901         const point& ep = points[e[1]];
 902
 903         const point p = nearestPointOnTheEdgeExact(sp, ep, centre);
 904         if( magSqr(p - centre) < rangeSq )
 905             edgesInRange.insert(it.key());
 906     }
 907
 908     Map<label> elGroup(edgesInRange.size());
 909     Map<label> pointTest(edgesInRange.size());
 910
 911     label nGroups(0);
 912
 913     DynList<label> front;
 914
 915     forAllConstIter(labelHashSet, edgesInRange, it)
 916         if( !elGroup.found(it.key()) )
 917         {
 918             front.clear();
 919             front.append(it.key());
 920             elGroup.insert(it.key(), nGroups);
 921
 922             while( front.size() )
 923             {
 924                 const label eLabel = front.removeLastElement();
 925                 const edge& e = edges[eLabel];
 926
 927                 for(label i=0;i<2;++i)
 928                 {
 929                     if( pointTest.found(e[i]) )
 930                     {
 931                         if( !pointTest[e[i]] )
 932                             continue;
 933                     }
 934                     else
 935                     {
 936                         if( magSqr(points[e[i]] - centre) < rangeSq )
 937                         {
 938                             pointTest.insert(e[i], 1);
 939                         }
 940                         else
 941                         {
 942                             pointTest.insert(e[i], 0);
 943                             continue;
 944                         }
 945                     }
 946
 947                     forAllRow(pointEdges, e[i], peI)
 948                     {
 949                         const label nei = pointEdges(e[i], peI);
 950
 951                         if( edgesInRange.found(nei) && !elGroup.found(nei) )
 952                         {
 953                             elGroup.insert(nei, nGroups);
 954                             front.append(nei);
 955                         }
 956                     }
 957                 }
 958             }
 959
 960             ++nGroups;
 961         }
 962
 963     return nGroups;
 964 }
 965
 966 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
 967
 968 } // End namespace help
 969
 970 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
 971
 972 } // End namespace Foam
 973
 974 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //