src/mesh/cfMesh/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 "Pair.H"
  40 #include "Map.H"
  41
  42 namespace Foam
  43 {
  44
  45 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
  46
  47 namespace help
  48 {
  49
  50 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
  51
  52 template<class ListType>
  53 inline bool isnan(const ListType& l)
  54 {
  55     forAll(l, i)
  56         if( l[i] != l[i] )
  57             return true;
  58
  59     return false;
  60 }
  61
  62 template<class ListType>
  63 bool isinf(const ListType& l)
  64 {
  65     forAll(l, i)
  66         if( (l[i] < -VGREAT) || (l[i] > VGREAT) )
  67             return true;
  68
  69     return false;
  70 }
  71
  72 template<class Face1, class Face2>
  73 inline bool isSharedEdgeConvex
  74 (
  75     const pointField& points,
  76     const Face1& f1,
  77     const Face2& f2
  78 )
  79 {
  80     DynList<label, 3> triOwn(3);
  81     DynList<label, 3> triNei(3);
  82
  83     forAll(f1, pI)
  84     {
  85         label pos(-1);
  86         forAll(f2, pJ)
  87             if( f2[pJ] == f1[pI] )
  88             {
  89                 pos = pJ;
  90                 break;
  91             }
  92
  93         if( pos < 0 )
  94             continue;
  95
  96         triNei[0] = f2[pos];
  97         triNei[1] = f2[f2.fcIndex(pos)];
  98         triNei[2] = f2[f2.rcIndex(pos)];
  99
 100         triOwn[0] = f1[pI];
 101         triOwn[1] = f1[f1.fcIndex(pI)];
 102         triOwn[2] = f1[f1.rcIndex(pI)];
 103
 104         scalar vol(0.0);
 105
 106         forAll(triOwn, pJ)
 107         {
 108             if( !triNei.contains(triOwn[pJ]) )
 109             {
 110                 tetrahedron<point, point> tet
 111                 (
 112                     points[triNei[0]],
 113                     points[triNei[1]],
 114                     points[triNei[2]],
 115                     points[triOwn[pJ]]
 116                 );
 117
 118                 vol = tet.mag();
 119                 break;
 120             }
 121         }
 122
 123         if( vol > -VSMALL )
 124             return false;
 125     }
 126
 127     return true;
 128 }
 129
 130 template<class Face1, class Face2>
 131 inline scalar angleBetweenFaces
 132 (
 133     const pointField& points,
 134     const Face1& f1,
 135     const Face2& f2
 136 )
 137 {
 138     DynList<label, 3> triOwn(3);
 139     DynList<label, 3> triNei(3);
 140
 141     scalar angle(0.0);
 142     label counter(0);
 143
 144     forAll(f1, pI)
 145     {
 146         label pos(-1);
 147         forAll(f2, pJ)
 148             if( f2[pJ] == f1[pI] )
 149             {
 150                 pos = pJ;
 151                 break;
 152             }
 153
 154         if( pos < 0 )
 155             continue;
 156
 157         triNei[0] = f2[pos];
 158         triNei[1] = f2[f2.fcIndex(pos)];
 159         triNei[2] = f2[f2.rcIndex(pos)];
 160
 161         triOwn[0] = f1[pI];
 162         triOwn[1] = f1[f1.fcIndex(pI)];
 163         triOwn[2] = f1[f1.rcIndex(pI)];
 164
 165         scalar vol(0.0);
 166
 167         forAll(triOwn, pJ)
 168         {
 169             if( !triNei.contains(triOwn[pJ]) )
 170             {
 171                 tetrahedron<point, point> tet
 172                 (
 173                     points[triNei[0]],
 174                     points[triNei[1]],
 175                     points[triNei[2]],
 176                     points[triOwn[pJ]]
 177                 );
 178
 179                 vol = tet.mag();
 180                 break;
 181             }
 182         }
 183
 184         vector nOwn
 185         (
 186             (points[triOwn[1]] - points[triOwn[0]]) ^
 187             (points[triOwn[2]] - points[triOwn[0]])
 188         );
 189         nOwn /= (mag(nOwn) + VSMALL);
 190
 191         vector nNei
 192         (
 193             (points[triNei[1]] - points[triNei[0]]) ^
 194             (points[triNei[2]] - points[triNei[0]])
 195         );
 196         nNei /= (mag(nNei) + VSMALL);
 197
 198         const scalar dot = Foam::max(-1.0, Foam::min(nOwn & nNei, 1.0));
 199
 200         if( vol > -VSMALL )
 201         {
 202             //- the angle is in the interval [Pi, 2Pi>
 203             const scalar ang = Foam::acos(dot);
 204
 205             angle += ang + M_PI;
 206             ++counter;
 207         }
 208         else
 209         {
 210             //- the angle is in the interval [0, Pi>
 211             const scalar ang = Foam::acos(-dot);
 212
 213             angle += ang;
 214             ++counter;
 215         }
 216     }
 217
 218     if( counter == 0 )
 219     {
 220         FatalErrorIn
 221         (
 222             "scalar angleBetweenFaces"
 223             "(const pointField&, const face&, const face&)"
 224         ) << "Faces " << f1 << " and " << f2
 225           << " do no share an edge" << abort(FatalError);
 226     }
 227
 228     return angle / counter;
 229 }
 230
 231 inline faceList mergePatchFaces
 232 (
 233     const List< DynList<label> >& pfcs,
 234     const pointField& polyPoints
 235 )
 236 {
 237     //- merge faces which share a common edge
 238     faceList patchFaces(pfcs.size());
 239     label counter(0);
 240     forAll(pfcs, faceI)
 241         if( pfcs[faceI].size() > 2 )
 242         {
 243             const DynList<label>& f = pfcs[faceI];
 244             face f_(f.size());
 245             forAll(f_, fJ)
 246                 f_[fJ] = f[fJ];
 247
 248             patchFaces[counter++] = f_;
 249         }
 250
 251     patchFaces.setSize(counter);
 252
 253     bool merged;
 254     do
 255     {
 256         faceList mergedFaces(patchFaces.size());
 257         boolList currentlyMerged(patchFaces.size(), false);
 258
 259         counter = 0;
 260         merged = false;
 261
 262         for(label nI=0;nI<(patchFaces.size()-1);nI++)
 263         {
 264             vector n0 = patchFaces[nI].normal(polyPoints);
 265             n0 /= mag(n0);
 266
 267             for(label nJ=nI+1;nJ<patchFaces.size();nJ++)
 268             {
 269                 vector n1 = patchFaces[nI].normal(polyPoints);
 270                 n1 /= mag(n1);
 271                 if(
 272                     help::shareAnEdge(patchFaces[nI], patchFaces[nJ]) &&
 273                     mag(n0 & n1) > 0.95
 274                 )
 275                 {
 276                     merged = true;
 277                     currentlyMerged[nI] = currentlyMerged[nJ] = true;
 278                     mergedFaces[counter++] =
 279                         help::mergeTwoFaces
 280                         (
 281                             patchFaces[nI],
 282                             patchFaces[nJ]
 283                         );
 284
 285                     break;
 286                 }
 287             }
 288
 289             if( merged ) break;
 290         }
 291
 292         forAll(patchFaces, pfI)
 293             if( !currentlyMerged[pfI] )
 294                 mergedFaces.newElmt(counter++) = patchFaces[pfI];
 295
 296         if( merged )
 297         {
 298             patchFaces.setSize(counter);
 299             for(label k=0;k<counter;k++)
 300                 patchFaces[k] = mergedFaces[k];
 301         }
 302
 303     } while( merged );
 304
 305     return patchFaces;
 306 }
 307
 308 inline bool vertexOnLine(const point& p, const edge& e, const pointField& ep)
 309 {
 310     vector v = e.vec(ep);
 311     v /= mag(v);
 312
 313     vector pv = p - ep[e.start()];
 314     pv /= mag(pv);
 315
 316     if( mag(pv & v) > (1.0-SMALL) )
 317         return true;
 318
 319     return false;
 320 }
 321
 322 inline bool vertexInPlane(const point& p, const plane& pl)
 323 {
 324     const vector& n = pl.normal();
 325     const point& fp = pl.refPoint();
 326
 327     vector d = p - fp;
 328     if( mag(d) > VSMALL )
 329         d /= mag(d);
 330
 331     if( mag(d & n) < SMALL )
 332         return true;
 333
 334     return false;
 335 }
 336
 337 inline bool planeIntersectsEdge
 338 (
 339     const point& start,
 340     const point& end,
 341     const plane& pl,
 342     point& intersection
 343 )
 344 {
 345     const vector v = end - start;
 346
 347     const vector& n = pl.normal();
 348     const point& fp = pl.refPoint();
 349
 350     if( (n & (v/mag(v))) < SMALL )
 351         return false;
 352
 353     const scalar t((n & (fp - start)) / (n & v));
 354
 355     if( (t > -SMALL) && (t < (1.0+SMALL)) )
 356     {
 357         intersection = start + v * t;
 358         return true;
 359     }
 360
 361     return false;
 362 }
 363
 364 inline bool pointInTetrahedron
 365 (
 366     const point& p,
 367     const tetrahedron<point, point>& tet
 368 )
 369 {
 370     const vector v0 = tet.a() - tet.d();
 371     const vector v1 = tet.b() - tet.d();
 372     const vector v2 = tet.c() - tet.d();
 373     const vector sp = p - tet.d();
 374
 375     matrix3D mat;
 376     FixedList<scalar, 3> source;
 377     for(label i=0;i<3;++i)
 378     {
 379         mat[i][0] = v0[i];
 380         mat[i][1] = v1[i];
 381         mat[i][2] = v2[i];
 382         source[i] = sp[i];
 383     }
 384
 385     //- check the determinant of the transformation
 386     const scalar det = mat.determinant();
 387
 388     if( mag(det) < VSMALL )
 389         return false;
 390
 391     //- get the coordinates of the point in the barycentric corrdinate system
 392     const scalar u0 = mat.solveFirst(source);
 393
 394     if( (u0 < -SMALL) || (u0 > (1.0+SMALL)) )
 395         return false;
 396
 397     const scalar u1 = mat.solveSecond(source);
 398
 399     if( (u1 < -SMALL) || ((u0+u1) > (1.0+SMALL)) )
 400         return false;
 401
 402     const scalar u2 = mat.solveThird(source);
 403
 404     if( (u2 < -SMALL) || (u2 > (1.0+SMALL)) )
 405         return false;
 406
 407
 408     const scalar u3 = 1.0 - u0 - u1 - u2;
 409
 410     if( (u3 < -SMALL) || (u3 > (1.0+SMALL)) )
 411         return false;
 412
 413     return true;
 414 }
 415
 416 inline bool nearestEdgePointToTheLine
 417 (
 418     const point& edgePoint0,
 419     const point& edgePoint1,
 420     const point& lp0,
 421     const point& lp1,
 422     point& nearestOnEdge,
 423     point& nearestOnLine
 424 )
 425 {
 426     const vector v = lp1 - lp0;
 427     const vector d = lp0 - edgePoint0;
 428     const vector e = edgePoint1 - edgePoint0;
 429
 430     const scalar vMag = mag(v);
 431     if( vMag < VSMALL )
 432         return false;
 433
 434     const scalar eMag = mag(e);
 435     if( eMag < VSMALL )
 436     {
 437         nearestOnEdge = edgePoint0;
 438         nearestOnLine = nearestPointOnTheEdge(lp0, lp1, nearestOnEdge);
 439         return true;
 440     }
 441
 442     if( mag((v/vMag) & (e/eMag)) > (1.0 - SMALL) )
 443         return false;
 444
 445     tensor mat(tensor::zero);
 446     mat.xx() = (v&v);
 447     mat.xy() = mat.yx() = -1.0 * (v&e);
 448     mat.yy() = (e&e);
 449     mat.zz() = SMALL;
 450
 451     vector source(vector::zero);
 452     source[0] = -1.0 * (d&v);
 453     source[1] = (d&e);
 454
 455     const vector sol = (inv(mat) & source);
 456
 457     nearestOnLine = lp0 + v * sol[0];
 458     if( sol[1] > 1.0 )
 459     {
 460         nearestOnEdge = edgePoint1;
 461     }
 462     else if( sol[1] < 0.0 )
 463     {
 464         nearestOnEdge = edgePoint0;
 465     }
 466     else
 467     {
 468         nearestOnEdge = edgePoint0 + e * sol[1];
 469     }
 470
 471     return true;
 472 }
 473
 474 inline point nearestPointOnTheTriangle
 475 (
 476     const triangle<point, point>& tri,
 477     const point& p
 478 )
 479 {
 480     const vector v0 = tri.b() - tri.a();
 481     const vector v1 = tri.c() - tri.a();
 482     const vector v2 = p - tri.a();
 483
 484     const scalar dot00 = (v0 & v0);
 485     const scalar dot01 = (v0 & v1);
 486     const scalar dot02 = (v0 & v2);
 487     const scalar dot11 = (v1 & v1);
 488     const scalar dot12 = (v1 & v2);
 489
 490     // Compute barycentric coordinates
 491     const scalar det = dot00 * dot11 - dot01 * dot01;
 492
 493     if( mag(det) < VSMALL )
 494     {
 495         point nearest(p);
 496         scalar dist(VGREAT);
 497
 498         FixedList<Pair<point>, 3> edges;
 499         edges[0] = Pair<point>(tri.a(), tri.b());
 500         edges[1] = Pair<point>(tri.b(), tri.c());
 501         edges[2] = Pair<point>(tri.c(), tri.a());
 502
 503         forAll(edges, eI)
 504         {
 505             const Pair<point>& e = edges[eI];
 506             const point np =
 507                 nearestPointOnTheEdge
 508                 (
 509                     e.first(),
 510                     e.second(),
 511                     p
 512                 );
 513
 514             if( magSqr(p - np) < dist )
 515             {
 516                 nearest = np;
 517                 dist = magSqr(p - np);
 518             }
 519         }
 520
 521         return nearest;
 522     }
 523
 524     const scalar u = (dot11 * dot02 - dot01 * dot12) / det;
 525     const scalar v = (dot00 * dot12 - dot01 * dot02) / det;
 526
 527     const point pProj = tri.a() + u * v0 + v * v1;
 528
 529     //- Check if point is in triangle
 530     if( (u >= -SMALL) && (v >= -SMALL) && ((u + v) <= (1.0+SMALL)) )
 531     {
 532         return pProj;
 533     }
 534     else
 535     {
 536         if( u < -SMALL )
 537         {
 538             const scalar ed = ((pProj - tri.a()) & v1) / (magSqr(v1) + VSMALL);
 539
 540             if( ed > 1.0 )
 541             {
 542                 return tri.c();
 543             }
 544             else if( ed < 0.0 )
 545             {
 546                 return tri.a();
 547             }
 548             else
 549             {
 550                 return (tri.a() + v1 * ed);
 551             }
 552         }
 553         else if( v < -SMALL )
 554         {
 555             const scalar ed = ((pProj - tri.a()) & v0) / (magSqr(v0) + VSMALL);
 556
 557             if( ed > 1.0 )
 558             {
 559                 return tri.b();
 560             }
 561             else if( ed < 0.0 )
 562             {
 563                 return tri.a();
 564             }
 565             else
 566             {
 567                 return (tri.a() + v0 * ed);
 568             }
 569         }
 570         else
 571         {
 572             const vector ev = tri.b() - tri.c();
 573             const scalar ed = ((pProj - tri.c()) & ev) / (magSqr(ev) + VSMALL);
 574
 575             if( ed > 1.0 )
 576             {
 577                 return tri.b();
 578             }
 579             else if( ed < 0.0 )
 580             {
 581                 return tri.c();
 582             }
 583             else
 584             {
 585                 return (tri.c() + ev * ed);
 586             }
 587         }
 588     }
 589
 590     return pProj;
 591 }
 592
 593 inline point nearestPointOnTheTriangle
 594 (
 595     const label tI,
 596     const triSurf& surface,
 597     const point& p
 598 )
 599 {
 600     const labelledTri& ltri = surface[tI];
 601     const pointField& points = surface.points();
 602
 603     const triangle<point, point> tri
 604     (
 605         points[ltri[0]],
 606         points[ltri[1]],
 607         points[ltri[2]]
 608     );
 609
 610     return nearestPointOnTheTriangle(tri, p);
 611 }
 612
 613 inline bool findMinimizerPoint
 614 (
 615     const DynList<point>& origins,
 616     const DynList<vector>& normals,
 617     point& pMin
 618 )
 619 {
 620     if( origins.size() != normals.size() )
 621         FatalErrorIn
 622         (
 623             "inline bool findMinimizerPoint"
 624             "(const DynList<point>&, const DynList<vector>&, point&)"
 625         ) << "Size of normals and origins do not match" << abort(FatalError);
 626
 627     tensor mat(tensor::zero);
 628     vector source(vector::zero);
 629
 630     forAll(origins, i)
 631     {
 632         //- use the normalized vector
 633         vector n = normals[i];
 634         n /= (mag(n) + VSMALL);
 635
 636         const tensor t = n * n;
 637
 638         mat += t;
 639
 640         source += (origins[i] & n) * n;
 641     }
 642
 643     const scalar determinant = Foam::mag(Foam::det(mat));
 644     if( determinant < SMALL )
 645         return false;
 646
 647     pMin = (inv(mat, determinant) & source);
 648
 649     return true;
 650 }
 651
 652 inline bool triLineIntersection
 653 (
 654     const triangle<point, point>& tria,
 655     const point& lineStart,
 656     const point& lineEnd,
 657     point& intersection
 658 )
 659 {
 660     const point& p0 = tria.a();
 661     const vector v(lineStart - lineEnd);
 662     const vector v0 = tria.b() - p0;
 663     const vector v1 = tria.c() - p0;
 664     const vector sp = lineStart - p0;
 665
 666     matrix3D mat;
 667     FixedList<scalar, 3> source;
 668     for(label i=0;i<3;++i)
 669     {
 670         mat[i][0] = v0[i];
 671         mat[i][1] = v1[i];
 672         mat[i][2] = v[i];
 673         source[i] = sp[i];
 674     }
 675
 676     const scalar det = mat.determinant();
 677
 678     if( mag(det) < SMALL )
 679         return false;
 680
 681     const scalar t = mat.solveThird(source);
 682
 683     if( (t < -SMALL) || (t > (1.0+SMALL)) )
 684         return false;
 685
 686     const scalar u0 = mat.solveFirst(source);
 687
 688     if( u0 < -SMALL )
 689         return false;
 690
 691     const scalar u1 = mat.solveSecond(source);
 692
 693     if( (u1 < -SMALL) || ((u0+u1) > (1.0+SMALL)) )
 694         return false;
 695
 696     intersection = lineStart - t * v;
 697     return true;
 698 }
 699
 700 inline bool triLineIntersection
 701 (
 702     const triSurf& surface,
 703     const label tI,
 704     const point& s,
 705     const point& e,
 706     point& intersection
 707 )
 708 {
 709     const pointField& pts = surface.points();
 710     const labelledTri& tri = surface[tI];
 711
 712     const triangle<point, point> tria
 713     (
 714         pts[tri[0]],
 715         pts[tri[1]],
 716         pts[tri[2]]
 717     );
 718
 719     return triLineIntersection(tria, s, e, intersection);
 720 }
 721
 722 inline bool boundBoxLineIntersection
 723 (
 724     const point& s,
 725     const point& e,
 726     const boundBox& bb
 727 )
 728 {
 729     scalar tMax(1.0+SMALL), tMin(-SMALL);
 730
 731     const vector v = e - s;
 732     const scalar d = mag(v);
 733
 734     //- check if the vector has length
 735     if( d < VSMALL )
 736     {
 737         if( bb.contains(s) )
 738             return true;
 739
 740         return false;
 741     }
 742
 743     const point& pMin = bb.min();
 744     const point& pMax = bb.max();
 745
 746     //- check coordinates
 747     for(label dir=0;dir<3;++dir)
 748     {
 749         const scalar vd = v[dir];
 750         const scalar sd = s[dir];
 751
 752         if( mag(vd) > (SMALL * d) )
 753         {
 754             if( vd >= 0.0 )
 755             {
 756                 tMin = Foam::max(tMin, (pMin[dir] - sd) / vd);
 757                 tMax = Foam::min(tMax, (pMax[dir] - sd) / vd);
 758             }
 759             else
 760             {
 761                 tMin = Foam::max(tMin, (pMax[dir] - sd) / vd);
 762                 tMax = Foam::min(tMax, (pMin[dir] - sd) / vd);
 763             }
 764         }
 765         else if( (sd < pMin[dir]) || (sd > pMax[dir]) )
 766         {
 767             return false;
 768         }
 769     }
 770
 771     if( (tMax - tMin) > -SMALL )
 772         return true;
 773
 774     return false;
 775 }
 776
 777 inline bool lineFaceIntersection
 778 (
 779     const point& sp,
 780     const point& ep,
 781     const face& f,
 782     const pointField& fp,
 783     point& intersection
 784 )
 785 {
 786     const point c = f.centre(fp);
 787
 788     forAll(f, pI)
 789     {
 790         const triangle<point, point> tri
 791         (
 792             fp[f[pI]],
 793             fp[f.nextLabel(pI)],
 794             c
 795         );
 796
 797         if( triLineIntersection(tri, sp, ep, intersection) )
 798             return true;
 799     }
 800
 801     return false;
 802 }
 803
 804 inline bool doFaceAndTriangleIntersect
 805 (
 806     const triSurf& surface,
 807     const label triI,
 808     const face& f,
 809     const pointField& facePoints
 810 )
 811 {
 812     const pointField& triPoints = surface.points();
 813
 814     const point centre = f.centre(facePoints);
 815     point intersection;
 816
 817     //- check if any triangle edge intersects the face
 818     const labelledTri& tri = surface[triI];
 819
 820     forAll(tri, eI)
 821     {
 822         const point& s = triPoints[tri[eI]];
 823         const point& e = triPoints[tri[(eI+1)%3]];
 824
 825         forAll(f, pI)
 826         {
 827             const triangle<point, point> tria
 828             (
 829                 facePoints[f[pI]],
 830                 facePoints[f.nextLabel(pI)],
 831                 centre
 832             );
 833
 834             const bool currIntersection =
 835                 help::triLineIntersection
 836                 (
 837                     tria,
 838                     s,
 839                     e,
 840                     intersection
 841                 );
 842
 843             if( currIntersection )
 844                 return true;
 845         }
 846     }
 847
 848     //- check if any face edges intersect the triangle
 849     forAll(f, pI)
 850     {
 851         const point& s = facePoints[f[pI]];
 852         const point& e = facePoints[f.nextLabel(pI)];
 853
 854         const bool intersected =
 855             help::triLineIntersection
 856             (
 857                 surface,
 858                 triI,
 859                 s,
 860                 e,
 861                 intersection
 862             );
 863
 864         if( intersected )
 865             return true;
 866     }
 867
 868     return false;
 869 }
 870
 871 inline bool doEdgesOverlap
 872 (
 873     const point& e0p0,
 874     const point& e0p1,
 875     const point& e1p0,
 876     const point& e1p1,
 877     FixedList<point, 2>& overlappingPart,
 878     const scalar distTol,
 879     const scalar cosTol
 880 )
 881 {
 882     if( distTol < 0.0 )
 883     {
 884         WarningIn
 885         (
 886             "inline bool doEdgesOverlap(const point, const point&,"
 887             "const point&, const point&, const scalar, const scalar,"
 888             "FixedList<point, 2>& overlappingPart)"
 889         ) << "Distance is not specified" << endl;
 890
 891         return false;
 892     }
 893
 894     vector e0 = e0p1 - e0p0;
 895     const scalar de0 = mag(e0);
 896     e0 /= (de0 + VSMALL);
 897
 898     vector e1 = e1p1 - e1p0;
 899     const scalar de1 = mag(e1);
 900     e1 /= (de1 + VSMALL);
 901
 902     //- check the angle deviation between the two vectors
 903     if( mag(e0 & e1) < cosTol )
 904         return false;
 905
 906     scalar t00 = (e1p0 - e0p0) & e0;
 907     scalar t01 = (e1p1 - e0p0) & e0;
 908
 909     scalar t10 = (e0p0 - e1p0) & e1;
 910     scalar t11 = (e0p1 - e1p0) & e1;
 911
 912     //- check if points are colinear within the tolerance
 913     if
 914     (
 915         (magSqr(e0p0 + t00*e0 - e1p0) <= distTol*distTol) ||
 916         (magSqr(e0p0 + t01*e0 - e1p1) <= distTol*distTol) ||
 917         (magSqr(e1p0 + t10*e1 - e0p0) <= distTol*distTol) ||
 918         (magSqr(e1p0 + t11*e1 - e0p1) <= distTol*distTol)
 919     )
 920     {
 921         vector vec = e0 + (((e0 & e1) > 0.0) ? e1 : -e1);
 922         vec /= (mag(vec) + VSMALL);
 923         const point origin = 0.25 * (e0p0 + e0p1 + e1p0 + e1p1);
 924
 925         //- calculate parameters on the regression line
 926         t00 = (e0p0 - origin) & vec;
 927         t01 = (e0p1 - origin) & vec;
 928         t10 = (e1p0 - origin) & vec;
 929         t11 = (e1p1 - origin) & vec;
 930
 931         //- check interval overlapping over the line
 932         const scalar t0Min = Foam::min(t00, t01);
 933         const scalar t0Max = Foam::max(t00, t01);
 934         const scalar t1Min = Foam::min(t10, t11);
 935         const scalar t1Max = Foam::max(t10, t11);
 936
 937         if( t1Min < t0Max )
 938         {
 939             overlappingPart[0] = origin + t1Min * vec;
 940             overlappingPart[1] = origin + t0Max * vec;
 941
 942             return true;
 943         }
 944         else if( t0Min < t1Max )
 945         {
 946             overlappingPart[0] = origin + t0Min * vec;
 947             overlappingPart[1] = origin + t1Max * vec;
 948
 949             return true;
 950         }
 951     }
 952
 953     return false;
 954 }
 955
 956 inline bool doTrianglesOverlap
 957 (
 958     const triangle<point, point>& tri0,
 959     const triangle<point, point>& tri1,
 960     DynList<point>& overlappingPolygon,
 961     const scalar distTol,
 962     const scalar cosTol
 963 )
 964 {
 965     if( distTol < 0.0 )
 966     {
 967         WarningIn
 968         (
 969             "inline bool doTrianglesOverlap(const triangle<point, point>&,"
 970             " const triangle<point, point>&, DynList<point>&,"
 971             " const scalar, const scalar)"
 972         ) << "Distance is not specified" << endl;
 973
 974         return false;
 975     }
 976
 977     vector n0 = tri0.normal();
 978     const scalar dn0 = mag(n0);
 979     n0 /= (dn0 + VSMALL);
 980
 981     if( dn0 < VSMALL )
 982         return false;
 983
 984     vector n1 = tri1.normal();
 985     const scalar dn1 = mag(n1);
 986     n1 /= (dn1 + VSMALL);
 987
 988     if( dn1 < VSMALL )
 989         return false;
 990
 991     //- check the angle deviation between the two vectors
 992     if( (mag(n0 & n1) < cosTol) && (dn0 >= VSMALL) && (dn1 >= VSMALL) )
 993         return false;
 994
 995     //- check if the two nearest points are within tolerance
 996     if
 997     (
 998         (mag((tri1.a() - tri0.a()) & n0) < distTol) ||
 999         (mag((tri1.b() - tri0.a()) & n0) < distTol) ||
1000         (mag((tri1.c() - tri0.a()) & n0) < distTol) ||
1001         (mag((tri0.a() - tri1.a()) & n1) < distTol) ||
1002         (mag((tri0.b() - tri1.a()) & n1) < distTol) ||
1003         (mag((tri0.c() - tri1.a()) & n1) < distTol)
1004     )
1005     {
1006         vector vec = n0 + ((n0 & n1) >= 0.0 ? n1 : -n1);
1007         vec /= (mag(vec) + VSMALL);
1008         const point origin = 0.5 * (tri0.centre() + tri1.centre());
1009
1010         overlappingPolygon.clear();
1011
1012         //- calculate max distance point from the origin
1013         point bestPoint = tri0.a();
1014         scalar distSq = magSqr(tri0.a() - origin);
1015
1016         scalar dSq = magSqr(tri0.b() - origin);
1017         if( dSq > distSq )
1018         {
1019             distSq = dSq;
1020             bestPoint = tri0.b();
1021         }
1022
1023         dSq = magSqr(tri0.c() - origin);
1024         if( dSq > distSq )
1025         {
1026             distSq = dSq;
1027             bestPoint = tri0.c();
1028         }
1029
1030         dSq = magSqr(tri1.a() - origin);
1031         if( dSq > distSq )
1032         {
1033             distSq = dSq;
1034             bestPoint = tri1.a();
1035         }
1036
1037         dSq = magSqr(tri1.b() - origin);
1038         if( dSq > distSq )
1039         {
1040             distSq = dSq;
1041             bestPoint = tri1.b();
1042         }
1043
1044         dSq = magSqr(tri1.c() - origin);
1045         if( dSq > distSq )
1046         {
1047             distSq = dSq;
1048             bestPoint = tri1.c();
1049         }
1050
1051         if( distSq < VSMALL )
1052             return false;
1053
1054         //- transform into planar coordinates
1055         vector x = (bestPoint - origin) - ((bestPoint - origin) & vec) * vec;
1056         x /= (mag(x) + VSMALL);
1057         vector y = vec ^ x;
1058
1059         DynList<point2D, 6> poly2D(3);
1060         poly2D[0] = point2D((tri0.a() - origin) & x, (tri0.a() - origin) & y);
1061         poly2D[1] = point2D((tri0.b() - origin) & x, (tri0.b() - origin) & y);
1062         poly2D[2] = point2D((tri0.c() - origin) & x, (tri0.c() - origin) & y);
1063
1064         FixedList<point2D, 3> t1Proj;
1065         t1Proj[0] = point2D((tri1.a() - origin) & x, (tri1.a() - origin) & y);
1066         t1Proj[1] = point2D((tri1.b() - origin) & x, (tri1.b() - origin) & y);
1067         t1Proj[2] = point2D((tri1.c() - origin) & x, (tri1.c() - origin) & y);
1068
1069         forAll(t1Proj, eI)
1070         {
1071             const vector2D vec = t1Proj[(eI+1)%3] - t1Proj[eI];
1072
1073             DynList<scalar, 6> distance(poly2D.size());
1074             forAll(poly2D, pI)
1075             {
1076                 const vector2D pVec = poly2D[pI] - t1Proj[eI];
1077                 distance[pI] = vec.y() * pVec.x() - vec.x() * pVec.y();
1078             }
1079
1080             DynList<point2D, 6> newPoly2D;
1081
1082             forAll(distance, pI)
1083             {
1084                 if( distance[pI] >= 0.0 )
1085                 {
1086                     newPoly2D.append(poly2D[pI]);
1087
1088                     if( distance.fcElement(pI) < 0.0 )
1089                     {
1090                         //- this is very sensitive to floaing point tolerances
1091                         const point2D newP =
1092                             (
1093                                 mag(distance[pI]) * poly2D.fcElement(pI) +
1094                                 mag(distance.fcElement(pI)) * poly2D[pI]
1095                             ) /
1096                             (
1097                                 mag(distance[pI]) +
1098                                 mag(distance.fcElement(pI)) +
1099                                 VSMALL
1100                             );
1101
1102                         newPoly2D.append(newP);
1103                     }
1104                 }
1105                 else if( distance.fcElement(pI) >= 0.0 )
1106                 {
1107                     //- this is very sensitive to floaing point tolerances
1108                     const point2D newP =
1109                         (
1110                             mag(distance[pI]) * poly2D.fcElement(pI) +
1111                             mag(distance.fcElement(pI)) * poly2D[pI]
1112                         ) /
1113                         (
1114                             mag(distance[pI]) +
1115                             mag(distance.fcElement(pI)) +
1116                             VSMALL
1117                         );
1118
1119                     newPoly2D.append(newP);
1120                 }
1121             }
1122
1123             poly2D = newPoly2D;
1124         }
1125
1126         //- check if the overlapping polygon exists
1127         if( poly2D.size() == 0 )
1128         {
1129             overlappingPolygon.clear();
1130             return false;
1131         }
1132
1133         //- fill the overlapping polygon
1134         overlappingPolygon.setSize(poly2D.size());
1135         forAll(poly2D, pI)
1136         {
1137             const point2D& pp = poly2D[pI];
1138             overlappingPolygon[pI] = origin + x * pp.x() + y * pp.y();
1139         }
1140
1141         return true;
1142     }
1143
1144     overlappingPolygon.clear();
1145     return false;
1146 }
1147
1148 inline bool doTrianglesIntersect
1149 (
1150     const triangle<point, point> &tri0,
1151     const triangle<point, point> &tri1,
1152     const scalar distTol
1153 )
1154 {
1155     if( distTol < 0.0 )
1156     {
1157         WarningIn
1158         (
1159             "inline bool doTrianglesIntersect(const triangle<point, point>&,"
1160             " const triangle<point, point>&,"
1161             " const scalar, const scalar)"
1162         ) << "Distance is not specified" << endl;
1163
1164         return false;
1165     }
1166
1167     //- find distances of points from the second triangle
1168     point np = nearestPointOnTheTriangle(tri1, tri0.a());
1169     if( magSqr(np - tri0.a()) < distTol*distTol )
1170         return true;
1171     np = nearestPointOnTheTriangle(tri1, tri0.b());
1172     if( magSqr(np - tri0.b()) < distTol*distTol )
1173         return true;
1174     np = nearestPointOnTheTriangle(tri1, tri0.c());
1175     if( magSqr(np - tri0.c()) < distTol*distTol )
1176         return true;
1177
1178     //- find distances of points from the first triangle
1179     np = nearestPointOnTheTriangle(tri0, tri1.a());
1180     if( magSqr(np - tri1.a()) < distTol*distTol )
1181         return true;
1182     np = nearestPointOnTheTriangle(tri0, tri1.b());
1183     if( magSqr(np - tri1.b()) < distTol*distTol )
1184         return true;
1185     np = nearestPointOnTheTriangle(tri0, tri1.c());
1186     if( magSqr(np - tri1.c()) < distTol*distTol )
1187         return true;
1188
1189     //- find edge intersections
1190     point intersection;
1191     if( triLineIntersection(tri0, tri1.a(), tri1.b(), intersection) )
1192         return true;
1193     if( triLineIntersection(tri0, tri1.b(), tri1.c(), intersection) )
1194         return true;
1195     if( triLineIntersection(tri0, tri1.c(), tri1.a(), intersection) )
1196         return true;
1197
1198     if( triLineIntersection(tri1, tri0.a(), tri0.b(), intersection) )
1199         return true;
1200     if( triLineIntersection(tri1, tri0.b(), tri0.c(), intersection) )
1201         return true;
1202     if( triLineIntersection(tri1, tri0.c(), tri0.a(), intersection) )
1203         return true;
1204
1205     return false;
1206 }
1207
1208 inline bool doTrianglesIntersect
1209 (
1210     const triangle<point, point>& tri0,
1211     const triangle<point, point>& tri1,
1212     DynList<point> &intersectionPoints,
1213     const scalar distTol
1214 )
1215 {
1216     if( distTol < 0.0 )
1217     {
1218         WarningIn
1219         (
1220             "inline bool doTrianglesIntersect(const triangle<point, point>&,"
1221             " const triangle<point, point>&, DynList<point>&,"
1222             " const scalar, const scalar)"
1223         ) << "Distance is not specified" << endl;
1224
1225         return false;
1226     }
1227
1228     vector n0 = tri0.normal();
1229     const scalar dn0 = mag(n0);
1230     n0 /= (dn0 + VSMALL);
1231
1232     vector n1 = tri1.normal();
1233     const scalar dn1 = mag(n1);
1234     n1 /= (dn1 + VSMALL);
1235
1236     //- distance of the points of the first triangle from the plane
1237     //- of the second triangle
1238     FixedList<scalar, 3> distancesTri0;
1239     distancesTri0[0] = (tri0.a() - tri1.a()) & n1;
1240     distancesTri0[1] = (tri0.b() - tri1.a()) & n1;
1241     distancesTri0[2] = (tri0.c() - tri1.a()) & n1;
1242
1243     forAll(distancesTri0, i)
1244         if( mag(distancesTri0[i]) < distTol )
1245             distancesTri0[i] = 0.0;
1246
1247     bool hasPositive(false), hasNegative(false), hasZero(false);
1248     DynList<point, 2> intersectionLine0;
1249     forAll(distancesTri0, pI)
1250     {
1251         if( distancesTri0[pI] >= distTol )
1252         {
1253             hasPositive = true;
1254         }
1255         else if( distancesTri0[pI] <= -distTol )
1256         {
1257             hasNegative = true;
1258         }
1259         else
1260         {
1261             hasZero = true;
1262         }
1263     }
1264
1265     //- find points on the intersection line
1266     if
1267     (
1268         (hasPositive && !hasNegative && !hasZero) ||
1269         (hasNegative && !hasPositive && !hasZero)
1270     )
1271     {
1272         //- there can be no intersection
1273         intersectionPoints.clear();
1274         return false;
1275     }
1276     else if
1277     (
1278         (hasPositive && (hasNegative || hasZero)) ||
1279         (hasNegative && (hasPositive || hasZero))
1280     )
1281     {
1282         //- find points on the intersection line
1283         if( mag(distancesTri0[0]) < distTol)
1284         {
1285             intersectionLine0.append(tri0.a());
1286         }
1287         else if( distancesTri0[0] * distancesTri0[1] < 0.0 )
1288         {
1289             intersectionLine0.append
1290             (
1291                 (tri0.a() * distancesTri0[1] - tri0.b() * distancesTri0[0]) /
1292                 (distancesTri0[0] - distancesTri0[1])
1293             );
1294         }
1295
1296         if( mag(distancesTri0[1]) < distTol )
1297         {
1298             intersectionLine0.append(tri0.b());
1299         }
1300         else if( distancesTri0[1] * distancesTri0[2] < 0.0 )
1301         {
1302             intersectionLine0.append
1303             (
1304                 (tri0.b() * distancesTri0[2] - tri0.c() * distancesTri0[1]) /
1305                 (distancesTri0[1] - distancesTri0[2])
1306             );
1307         }
1308
1309         if( mag(distancesTri0[2]) < distTol )
1310         {
1311             intersectionLine0.append(tri0.c());
1312         }
1313         else if( distancesTri0[2] * distancesTri0[0] < 0.0 )
1314         {
1315             intersectionLine0.append
1316             (
1317                 (tri0.c() * distancesTri0[0] - tri0.a() * distancesTri0[2]) /
1318                 (distancesTri0[2] - distancesTri0[0])
1319             );
1320         }
1321     }
1322     else
1323     {
1324         //- these triangles are in the same plane
1325         //- check if they overlap
1326         return doTrianglesOverlap(tri0, tri1, intersectionPoints, distTol);
1327     }
1328
1329     //- distance of the points of the second triangle from the plane
1330     //- of the first triangle
1331     FixedList<scalar, 3> distancesTri1;
1332     distancesTri1[0] = (tri1.a() - tri0.a()) & n0;
1333     distancesTri1[1] = (tri1.b() - tri0.a()) & n0;
1334     distancesTri1[2] = (tri1.c() - tri0.a()) & n0;
1335
1336     forAll(distancesTri1, i)
1337         if( mag(distancesTri1[i]) < distTol )
1338             distancesTri1[i] = 0.0;
1339
1340     hasPositive = false;
1341     hasNegative = false;
1342     hasZero = false;
1343
1344     DynList<point, 2> intersectionLine1;
1345     forAll(distancesTri1, pI)
1346     {
1347         if( distancesTri1[pI] >= distTol )
1348         {
1349             hasPositive = true;
1350         }
1351         else if( distancesTri1[pI] <= -distTol )
1352         {
1353             hasNegative = true;
1354         }
1355         else
1356         {
1357             hasZero = true;
1358         }
1359     }
1360
1361     if
1362     (
1363         (hasPositive && !hasNegative && !hasZero) ||
1364         (hasNegative && !hasPositive && !hasZero)
1365     )
1366     {
1367         //- there can be no intersection
1368         intersectionPoints.clear();
1369         return false;
1370     }
1371     else if
1372     (
1373         (hasPositive && (hasNegative || hasZero)) ||
1374         (hasNegative && (hasPositive || hasZero))
1375     )
1376     {
1377         //- find points on the intersection line
1378         if( mag(distancesTri1[0]) < distTol)
1379         {
1380             intersectionLine1.append(tri1.a());
1381         }
1382         else if( distancesTri1[0] * distancesTri1[1] < 0.0 )
1383         {
1384             intersectionLine1.append
1385             (
1386                 (tri1.a() * distancesTri1[1] - tri1.b() * distancesTri1[0]) /
1387                 (distancesTri1[0] - distancesTri1[1])
1388             );
1389         }
1390
1391         if( mag(distancesTri1[1]) < distTol)
1392         {
1393             intersectionLine1.append(tri1.b());
1394         }
1395         else if( distancesTri1[1] * distancesTri1[2] < 0.0 )
1396         {
1397             intersectionLine1.append
1398             (
1399                 (tri1.b() * distancesTri1[2] - tri1.c() * distancesTri1[1]) /
1400                 (distancesTri1[1] - distancesTri1[2])
1401             );
1402         }
1403
1404         if( mag(distancesTri1[2]) < distTol)
1405         {
1406             intersectionLine1.append(tri1.c());
1407         }
1408         else if( distancesTri1[2] * distancesTri1[0] < 0.0 )
1409             intersectionLine1.append
1410             (
1411                 (tri1.c() * distancesTri1[0] - tri1.a() * distancesTri1[2]) /
1412                 (distancesTri1[2] - distancesTri1[0])
1413             );
1414     }
1415     else
1416     {
1417         //- these triangles are in the same plane
1418         //- check if they overlap
1419         return doTrianglesOverlap(tri0, tri1, intersectionPoints, distTol);
1420     }
1421
1422     //- calculate the direction of the intersection line
1423     vector vec = n0 ^ n1;
1424
1425     //- n0 and n1 are nearly coplanar, try overlap intersection
1426     if( magSqr(vec) < SMALL )
1427     {
1428         //- these triangles are in the same plane
1429         return doTrianglesOverlap(tri0, tri1, intersectionPoints, distTol);
1430     }
1431
1432     //- calculate intervals on the intersection line
1433     scalar t0Min(VGREAT), t0Max(-VGREAT);
1434     point p0Min(vector::zero), p0Max(vector::zero);
1435     forAll(intersectionLine0, i)
1436     {
1437         const scalar t = intersectionLine0[i] & vec;
1438
1439         if( t < t0Min )
1440         {
1441             t0Min = t;
1442             p0Min = intersectionLine0[i];
1443         }
1444         if( t > t0Max )
1445         {
1446             t0Max = t;
1447             p0Max = intersectionLine0[i];
1448         }
1449     }
1450
1451     scalar t1Min(VGREAT), t1Max(-VGREAT);
1452     point p1Min(vector::zero), p1Max(vector::zero);
1453     forAll(intersectionLine1, i)
1454     {
1455         const scalar t = intersectionLine1[i] & vec;
1456
1457         if( t < t1Min )
1458         {
1459             t1Min = t;
1460             p1Min = intersectionLine1[i];
1461         }
1462         if( t > t1Max )
1463         {
1464             t1Max = t;
1465             p1Max = intersectionLine1[i];
1466         }
1467     }
1468
1469     if( (t1Min <= t0Max) && (t1Max >= t0Min) )
1470     {
1471         intersectionPoints.setSize(2);
1472         intersectionPoints[0] = p1Min;
1473         intersectionPoints[1] = p0Max;
1474
1475         return true;
1476     }
1477     else if( (t0Min <= t1Max) && (t0Max >= t1Min) )
1478     {
1479         intersectionPoints.setSize(2);
1480         intersectionPoints[0] = p0Min;
1481         intersectionPoints[1] = p1Max;
1482
1483         return true;
1484     }
1485
1486     intersectionPoints.setSize(0);
1487     return false;
1488 }
1489
1490 inline bool pointInsideFace
1491 (
1492     const point& p,
1493     const face& f,
1494     const vector& n,
1495     const pointField& fp,
1496     const scalar distTol
1497 )
1498 {
1499     const edgeList fe = f.edges();
1500     forAll(f, pI)
1501     {
1502         if( mag(p - fp[f[pI]]) < distTol )
1503             return true;
1504
1505         vector pv = p - fp[f[pI]];
1506         pv /= mag(pv);
1507
1508         vector lv = n ^ fe[pI].vec(fp);
1509         lv /= mag(lv);
1510
1511         const scalar d = pv & lv;
1512         if( d < -distTol )
1513         {
1514             return false;
1515         }
1516     }
1517
1518     return true;
1519 }
1520
1521 inline bool pointInsideFace
1522 (
1523     const point& p,
1524     const face& f,
1525     const pointField& fp,
1526     const scalar distTol
1527 )
1528 {
1529     const point c = f.centre(fp);
1530     const scalar tolSqr = sqr(distTol);
1531
1532     forAll(f, eI)
1533     {
1534         const edge fe = f.faceEdge(eI);
1535
1536         triangle<point, point> tri
1537         (
1538             fp[fe[0]],
1539             fp[fe[1]],
1540             c
1541         );
1542
1543         const point np = nearestPointOnTheTriangle(tri, p);
1544
1545         if( magSqr(np - p) <= tolSqr )
1546             return true;
1547     }
1548
1549     return false;
1550 }
1551
1552 inline bool isFaceConvexAndOk
1553 (
1554     const face& f,
1555     const pointField& fp,
1556     DynList<bool>& OkPoints
1557 )
1558 {
1559     bool valid(true);
1560
1561     vector normal = f.normal(fp);
1562     const scalar magN = mag(normal);
1563
1564     //- face has zero area. All points are inverted
1565     if( magN < VSMALL )
1566     {
1567         OkPoints.setSize(f.size());
1568         OkPoints = false;
1569
1570         return false;
1571     }
1572
1573     normal /= (magN + VSMALL);
1574
1575     //- project points into the plane formed by the face centre
1576     //- and the face normal
1577     plane pl(f.centre(fp), normal);
1578
1579     DynList<point> projFace(f.size());
1580     forAll(f, pI)
1581         projFace[pI] = pl.nearestPoint(fp[f[pI]]);
1582
1583     //- calculate normalised edge vectors
1584     DynList<vector> edgeVecs(f.size());
1585     forAll(f, eI)
1586     {
1587         vector& v = edgeVecs[eI];
1588
1589         const edge e = f.faceEdge(eI);
1590         v = e.vec(fp);
1591         v /= (mag(v) + VSMALL);
1592     }
1593
1594     //- check if the face is convex
1595     OkPoints.setSize(f.size());
1596     forAll(f, pI)
1597     {
1598         const vector& pv = edgeVecs[f.rcIndex(pI)];
1599         const vector& nv = edgeVecs[pI];
1600
1601         if( ((pv ^ nv) & normal) >= -0.05 )
1602         {
1603             OkPoints[pI] = true;
1604         }
1605         else
1606         {
1607             OkPoints[pI] = false;
1608             valid = false;
1609         }
1610     }
1611
1612     return valid;
1613 }
1614
1615 inline point nearestPointOnTheEdge
1616 (
1617     const point& edgePoint0,
1618     const point& edgePoint1,
1619     const point& p
1620 )
1621 {
1622     const vector e = edgePoint1 - edgePoint0;
1623     const scalar d = mag(e);
1624     const vector k = p - edgePoint0;
1625
1626     if( d < ROOTVSMALL )
1627         return edgePoint0;
1628
1629     return edgePoint0 + ((e / (d*d)) * (e & k));
1630 }
1631
1632 inline point nearestPointOnTheEdgeExact
1633 (
1634     const point& edgePoint0,
1635     const point& edgePoint1,
1636     const point& p
1637 )
1638 {
1639     const vector e = edgePoint1 - edgePoint0;
1640     const scalar dSq = magSqr(e);
1641     const vector k = p - edgePoint0;
1642
1643     if( dSq < VSMALL )
1644         return edgePoint0;
1645
1646     const scalar t = (e & k) / (dSq + VSMALL);
1647     if( t > 1.0 )
1648     {
1649         return edgePoint1;
1650     }
1651     else if( t < 0.0 )
1652     {
1653         return edgePoint0;
1654     }
1655
1656     return edgePoint0 + (e * t);
1657 }
1658
1659 inline scalar distanceOfPointFromTheEdge
1660 (
1661     const point& edgePoint0,
1662     const point& edgePoint1,
1663     const point& p
1664 )
1665 {
1666     return mag(nearestPointOnTheEdge(edgePoint0, edgePoint1, p) - p);
1667 }
1668
1669 inline label numberOfFaceGroups
1670 (
1671     const labelHashSet& containedElements,
1672     const point& centre,
1673     const scalar range,
1674     const triSurf& surface
1675 )
1676 {
1677     const pointField& points = surface.points();
1678     const edgeLongList& edges = surface.edges();
1679     const VRWGraph& faceEdges = surface.facetEdges();
1680     const VRWGraph& edgeFaces = surface.edgeFacets();
1681
1682     labelHashSet triaInRange(containedElements.size());
1683     const scalar rangeSq = range * range;
1684     forAllConstIter(labelHashSet, containedElements, it)
1685     {
1686         const point p = nearestPointOnTheTriangle(it.key(), surface, centre);
1687         if( magSqr(p - centre) < rangeSq )
1688             triaInRange.insert(it.key());
1689     }
1690
1691     Map<label> elGroup(triaInRange.size());
1692     Map<label> testEdge(triaInRange.size());
1693
1694     label nGroups(0);
1695
1696     DynList<label> front;
1697
1698     forAllConstIter(labelHashSet, triaInRange, it)
1699         if( !elGroup.found(it.key()) )
1700         {
1701             front.clear();
1702             front.append(it.key());
1703             elGroup.insert(it.key(), nGroups);
1704
1705             while( front.size() )
1706             {
1707                 const label fLabel = front.removeLastElement();
1708
1709                 forAllRow(faceEdges, fLabel, feI)
1710                 {
1711                     const label edgeI = faceEdges(fLabel, feI);
1712
1713                     //- check if the edge intersects the bounding box
1714                     if( testEdge.found(edgeI) )
1715                     {
1716                         if( !testEdge[edgeI] )
1717                             continue;
1718                     }
1719                     else
1720                     {
1721                         const point& s = points[edges[edgeI][0]];
1722                         const point& e = points[edges[edgeI][1]];
1723                         const point np = nearestPointOnTheEdge(s, e, centre);
1724                         if( magSqr(np - centre) < rangeSq )
1725                         {
1726                             testEdge.insert(edgeI, 1);
1727                         }
1728                         else
1729                         {
1730                             testEdge.insert(edgeI, 0);
1731                             continue;
1732                         }
1733                     }
1734
1735                     forAllRow(edgeFaces, edgeI, efI)
1736                     {
1737                         const label nei = edgeFaces(edgeI, efI);
1738                         if
1739                         (
1740                             triaInRange.found(nei) &&
1741                             !elGroup.found(nei)
1742                         )
1743                         {
1744                             elGroup.insert(nei, nGroups);
1745                             front.append(nei);
1746                         }
1747                     }
1748                 }
1749             }
1750
1751             ++nGroups;
1752         }
1753
1754     return nGroups;
1755 }
1756
1757 inline label numberOfEdgeGroups
1758 (
1759     const labelHashSet& containedEdges,
1760     const point& centre,
1761     const scalar range,
1762     const triSurf& surface
1763 )
1764 {
1765     const pointField& points = surface.points();
1766     const edgeLongList& edges = surface.edges();
1767     const VRWGraph& pointEdges = surface.pointEdges();
1768
1769     const scalar rangeSq = range * range;
1770     labelHashSet edgesInRange(containedEdges.size());
1771     forAllConstIter(labelHashSet, containedEdges, it)
1772     {
1773         const edge& e = edges[it.key()];
1774         const point& sp = points[e[0]];
1775         const point& ep = points[e[1]];
1776
1777         const point p = nearestPointOnTheEdgeExact(sp, ep, centre);
1778         if( magSqr(p - centre) < rangeSq )
1779             edgesInRange.insert(it.key());
1780     }
1781
1782     Map<label> elGroup(edgesInRange.size());
1783     Map<label> pointTest(edgesInRange.size());
1784
1785     label nGroups(0);
1786
1787     DynList<label> front;
1788
1789     forAllConstIter(labelHashSet, edgesInRange, it)
1790         if( !elGroup.found(it.key()) )
1791         {
1792             front.clear();
1793             front.append(it.key());
1794             elGroup.insert(it.key(), nGroups);
1795
1796             while( front.size() )
1797             {
1798                 const label eLabel = front.removeLastElement();
1799                 const edge& e = edges[eLabel];
1800
1801                 for(label i=0;i<2;++i)
1802                 {
1803                     if( pointTest.found(e[i]) )
1804                     {
1805                         if( !pointTest[e[i]] )
1806                             continue;
1807                     }
1808                     else
1809                     {
1810                         if( magSqr(points[e[i]] - centre) < rangeSq )
1811                         {
1812                             pointTest.insert(e[i], 1);
1813                         }
1814                         else
1815                         {
1816                             pointTest.insert(e[i], 0);
1817                             continue;
1818                         }
1819                     }
1820
1821                     forAllRow(pointEdges, e[i], peI)
1822                     {
1823                         const label nei = pointEdges(e[i], peI);
1824
1825                         if( edgesInRange.found(nei) && !elGroup.found(nei) )
1826                         {
1827                             elGroup.insert(nei, nGroups);
1828                             front.append(nei);
1829                         }
1830                     }
1831                 }
1832             }
1833
1834             ++nGroups;
1835         }
1836
1837     return nGroups;
1838 }
1839
1840 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
1841
1842 } // End namespace help
1843
1844 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
1845
1846 } // End namespace Foam
1847
1848 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //