src/foam/interpolations/GGIInterpolation/GGIInterpolationQuickRejectTests.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     |
   5     \\  /    A nd           | For copyright notice see file Copyright
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     3D and 2D quick reject tests for filtering out the neighborhood
  26     of the GGI master patch faces
  27
  28 Author
  29     Martin Beaudoin, Hydro-Quebec, (2008)
  30
  31 \*---------------------------------------------------------------------------*/
  32
  33 #include "boundBox.H"
  34 #include "plane.H"
  35 #include "transformField.H"
  36 #include "octree.H"
  37 #include "octreeDataBoundBox.H"
  38
  39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  40
  41 namespace Foam
  42 {
  43
  44 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  45
  46 template<class MasterPatch, class SlavePatch>
  47 const Foam::debug::tolerancesSwitch
  48 GGIInterpolation<MasterPatch, SlavePatch>::faceBoundBoxExtendSpanFraction_
  49 (
  50     "GGIFaceBoundBoxExtendSpanFraction",
  51     1.0e-2,
  52     "GGI faces bounding box expansion factor. Add robustness for quick-search algo. Keep it to a few percent."
  53 );
  54
  55 template<class MasterPatch, class SlavePatch>
  56 const Foam::debug::optimisationSwitch
  57 GGIInterpolation<MasterPatch, SlavePatch>::octreeSearchMinNLevel_
  58 (
  59     "GGIOctreeSearchMinNLevel",
  60     3,
  61     "GGI neighbouring facets octree-based search: minNlevel parameter for octree"
  62 );
  63
  64 template<class MasterPatch, class SlavePatch>
  65 const Foam::debug::optimisationSwitch
  66 GGIInterpolation<MasterPatch, SlavePatch>::octreeSearchMaxLeafRatio_
  67 (
  68     "GGIOctreeSearchMaxLeafRatio",
  69     3,
  70     "GGI neighbouring facets octree-based search: maxLeafRatio parameter for octree"
  71 );
  72
  73 template<class MasterPatch, class SlavePatch>
  74 const Foam::debug::optimisationSwitch
  75 GGIInterpolation<MasterPatch, SlavePatch>::octreeSearchMaxShapeRatio_
  76 (
  77     "GGIOctreeSearchMaxShapeRatio",
  78     1,
  79     "GGI neighbouring facets octree-based search: maxShapeRatio parameter for octree"
  80 );
  81
  82
  83 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  84
  85 // From: http://www.gamasutra.com/features/20000330/bobic_02.htm
  86 // One of the most primitive ways of doing collision detection is to
  87 // approximate each object or a part of the object with a sphere, and
  88 // then check whether spheres intersect each other. This method is
  89 // widely used even today because it’s computationally inexpensive.
  90 // We merely check whether the distance between the centers of two
  91 // spheres is less than the sum of the two radii (which indicates that
  92 // a collision has occurred). Even better, if we calculate whether the
  93 // distance squared is less than the sum of the radii squared, then we
  94 // eliminate that nasty square root in our distance calculation, but
  95 // only if we take care of this:
  96 //
  97 // In the interval [0,1], sqr(x) underestimate x, so we will have
  98 // false negative for the neighbors, this is bad...  In the interval
  99 // ]1, infinite], sqr(x) will overestimate x, so we will get false
 100 // positive, which is not that bad for a quick reject test
 101 //
 102 // So we will check the values of the sqr(x) before doing comparisons,
 103 // if sqr(x) is < 1, then we will pay for the sqrt. If sqr(x) > 1,
 104 // then we only compare squared values.
 105 //
 106 // Of course, we will end up comparing squared and and no-squared
 107 // values, which is a bit weird; but we must remember that this is a
 108 // quick reject test, and nothing more.
 109 //
 110 // Overall, we will overestimate the number of neighbours, which is
 111 // much more preferable than to underestimate. We will have to use
 112 // another test to remove the false positives (Separating Axis Theorem
 113 // test)
 114 //
 115 // This constitutes the first good quick reject test before going into
 116 // more computing intensive tasks with the calculation of the GGI
 117 // weights
 118
 119 template<class MasterPatch, class SlavePatch>
 120 void GGIInterpolation<MasterPatch, SlavePatch>::findNeighbours3D
 121 (
 122     labelListList& result
 123 ) const
 124 {
 125     List<DynamicList<label> > candidateMasterNeighbors(masterPatch_.size());
 126
 127     // First, compute the face center and the sphere radius (squared)
 128     // of the slave patch faces so we will not have to recompute this
 129     // n times
 130     vectorField slaveFaceCentre(slavePatch_.size());
 131     scalarField slaveRadius2(slavePatch_.size());
 132
 133     forAll (slavePatch_, faceSi)
 134     {
 135         // Grab points from slave face
 136         pointField curFacePoints =
 137             slavePatch_[faceSi].points(slavePatch_.points());
 138
 139         slaveFaceCentre[faceSi] =
 140             slavePatch_[faceSi].centre(slavePatch_.points());
 141
 142         if (doTransform())
 143         {
 144             // Transform points and normal to master plane
 145
 146             if (forwardT_.size() == 1)
 147             {
 148                 // Constant transform
 149                 transform
 150                 (
 151                     curFacePoints,
 152                     forwardT_[0],
 153                     curFacePoints
 154                 );
 155
 156                 slaveFaceCentre[faceSi] = transform
 157                 (
 158                     forwardT_[0],
 159                     slaveFaceCentre[faceSi]
 160                 );
 161             }
 162             else
 163             {
 164                 transform
 165                 (
 166                     curFacePoints,
 167                     forwardT_[faceSi],
 168                     curFacePoints
 169                 );
 170
 171                 slaveFaceCentre[faceSi] = transform
 172                 (
 173                     forwardT_[faceSi],
 174                     slaveFaceCentre[faceSi]
 175                 );
 176             }
 177         }
 178
 179         boundBox bbSlave(curFacePoints, false);
 180
 181         scalar tmpValue = Foam::magSqr(bbSlave.max() - bbSlave.min())/4.0;
 182
 183         // We will compare squared distances, so save the sqrt() if value > 1.0.
 184         if (tmpValue < 1.0)
 185         {
 186             // Take the sqrt, otherwise, we underestimate the radius
 187             slaveRadius2[faceSi] = sqrt(tmpValue);
 188         }
 189         else
 190         {
 191             slaveRadius2[faceSi] = tmpValue;
 192         }
 193     }
 194
 195     // Next, we search for each master face a list of potential neighbors
 196     forAll (masterPatch_, faceMi)
 197     {
 198         // For each masterPatch faces, compute the bounding box. With
 199         // this, we compute the radius of the bounding sphere for this
 200         // face
 201         boundBox bbMaster
 202         (
 203             masterPatch_[faceMi].points(masterPatch_.points()),
 204             false
 205         );
 206
 207         // We will compare squared distances, so save the sqrt() if > 1.0.
 208         scalar masterRadius2 =
 209             Foam::magSqr(bbMaster.max() - bbMaster.min())/4.0;
 210
 211         // Take the sqrt, otherwise, we underestimate the radius
 212         if (masterRadius2 < 1.0)
 213         {
 214             masterRadius2 = sqrt(masterRadius2);
 215         }
 216
 217         // This is the inner loop, so the face centres and face radii
 218         // for the slave faces are already pre-computed
 219         forAll (slavePatch_, faceSi)
 220         {
 221             // We test if the squared distance between the two face
 222             // centers is less than the sum of the 2 squared radii from
 223             // each face bounding sphere.
 224
 225             // If so, we have found possibly 2 neighbours
 226             scalar distFaceCenters =
 227                 Foam::magSqr
 228                 (
 229                     masterPatch_[faceMi].centre(masterPatch_.points())
 230                   - slaveFaceCentre[faceSi]
 231                 );
 232
 233             if (distFaceCenters < 1.0)
 234             {
 235                 distFaceCenters = sqrt(distFaceCenters);
 236             }
 237
 238             if (distFaceCenters < (masterRadius2 + slaveRadius2[faceSi]))
 239             {
 240                 candidateMasterNeighbors[faceMi].append(faceSi);
 241             }
 242         }
 243     }
 244
 245     // Repack the list
 246     result.setSize(masterPatch_.size());
 247
 248     forAll (result, i)
 249     {
 250         result[i].transfer(candidateMasterNeighbors[i].shrink());
 251     }
 252 }
 253
 254
 255 // This algorithm find the faces in proximity of another face based
 256 // on the AABB (Axis Aligned Boundary Box. These tests are done in 3D
 257 // space.
 258 //
 259 // For a given face and a potential neighbours, we construct first an
 260 // "augmented BB" by substracting/adding the slave delta BB to the
 261 // BBmin/BBmax of the master face. Then, we compare if the "augmented
 262 // BB" intersects with the potential neighbour BB.
 263 //
 264 // This is a fairly quick reject test, based only on additions and
 265 // comparisons... But still a n^2 test... very costly
 266 //
 267 // Furthermore, before selecting a given face as a potential
 268 // neighbour, we evaluate the featureCos of both master and slave
 269 // face normals. We reject immediatly all the faces that might fall
 270 // into the bounding box, but where the normals are too dissimilar.
 271
 272 template<class MasterPatch, class SlavePatch>
 273 void GGIInterpolation<MasterPatch, SlavePatch>::findNeighboursAABB
 274 (
 275     labelListList& result
 276 ) const
 277 {
 278     List<DynamicList<label> > candidateMasterNeighbors(masterPatch_.size());
 279
 280     // Grab the master patch faces bounding boxes
 281     List<boundBox> masterPatchBB(masterPatch_.size());
 282
 283     forAll (masterPatch_, faceMi)
 284     {
 285         boundBox bbMaster
 286         (
 287             masterPatch_[faceMi].points(masterPatch_.points()),
 288             false
 289         );
 290
 291         masterPatchBB[faceMi] = bbMaster;
 292     }
 293
 294     // Grab the slave patch faces bounding boxes, plus compute its
 295     // extend or delta
 296     List<boundBox> slavePatchBB(slavePatch_.size());
 297     pointField deltaBBSlave(slavePatch_.size());
 298
 299     // We expect that any neighbour face to face intersection will fall
 300     // within augmented BB.
 301     vectorField slaveFaceBBminThickness(slavePatch_.size(), vector::zero);
 302
 303     const faceList& slaveLocalFaces = slavePatch_.localFaces();
 304     vectorField slaveNormals = slavePatch_.faceNormals();
 305     const pointField& slaveLocalPoints = slavePatch_.localPoints();
 306
 307     // Transform slave normals to master plane if needed
 308     if (doTransform())
 309     {
 310         if (forwardT_.size() == 1)
 311         {
 312             transform(slaveNormals, forwardT_[0], slaveNormals);
 313         }
 314         else
 315         {
 316             transform(slaveNormals, forwardT_, slaveNormals);
 317         }
 318     }
 319
 320     forAll (slaveFaceBBminThickness, sI)
 321     {
 322         scalar maxEdgeLength = 0.0;
 323
 324         // Let's use the length of the longest edge from each faces
 325         edgeList el = slaveLocalFaces[sI].edges();
 326
 327         forAll (el, elI)
 328         {
 329             scalar edgeLength = el[elI].mag(slaveLocalPoints);
 330             maxEdgeLength = Foam::max(edgeLength, maxEdgeLength);
 331         }
 332
 333         // Make sure our offset is positive. Ugly, but cheap
 334         vector posNormal = cmptMag(slaveNormals[sI]);
 335
 336         slaveFaceBBminThickness[sI] = posNormal*maxEdgeLength;
 337     }
 338
 339     // Iterate over slave patch faces, compute its bounding box,
 340     // using a possible transformation and separation for cyclic patches
 341     forAll (slavePatch_, faceSi)
 342     {
 343         pointField curFacePoints =
 344             slavePatch_[faceSi].points(slavePatch_.points());
 345
 346         if (doTransform())
 347         {
 348             if (forwardT_.size() == 1)
 349             {
 350                 transform(curFacePoints, forwardT_[0], curFacePoints);
 351             }
 352             else
 353             {
 354                 transform(curFacePoints, forwardT_[faceSi], curFacePoints);
 355             }
 356         }
 357
 358         if (doSeparation())
 359         {
 360             if (forwardSep_.size() == 1)
 361             {
 362                 curFacePoints += forwardSep_[0];
 363             }
 364             else
 365             {
 366                 curFacePoints += forwardSep_[faceSi];
 367             }
 368         }
 369
 370         slavePatchBB[faceSi] = boundBox(curFacePoints, false);
 371
 372         // We compute the extent of the slave face BB.
 373         // Plus, we boost it a little bit, just to stay clear
 374         // of floating point numerical issues when doing intersections
 375         // Let's boost by 10%.
 376         deltaBBSlave[faceSi] =
 377             1.1*
 378             (
 379                 slavePatchBB[faceSi].max()
 380               - slavePatchBB[faceSi].min()
 381               + slaveFaceBBminThickness[faceSi]
 382             );
 383     }
 384
 385     // Visit each master patch face BB,
 386     // augment it with the info from the slave patch face BB
 387     // then, compute the intersection
 388     const vectorField& masterFaceNormals = masterPatch_.faceNormals();
 389
 390     forAll (masterPatchBB, faceMi)
 391     {
 392         forAll (slavePatchBB, faceSi)
 393         {
 394             // Compute the augmented AABB
 395             boundBox augmentedBBMaster
 396             (
 397                 masterPatchBB[faceMi].min() - deltaBBSlave[faceSi],
 398                 masterPatchBB[faceMi].max() + deltaBBSlave[faceSi]
 399             );
 400
 401             if (augmentedBBMaster.overlaps(slavePatchBB[faceSi]))
 402             {
 403                 // Compute featureCos between the two face normals
 404                 // before adding to the list of candidates
 405                 scalar featureCos =
 406                     masterFaceNormals[faceMi] & slaveNormals[faceSi];
 407
 408                 if (mag(featureCos) > featureCosTol_)
 409                 {
 410                     candidateMasterNeighbors[faceMi].append(faceSi);
 411                 }
 412             }
 413         }
 414     }
 415
 416     // Repack the list
 417     result.setSize(masterPatch_.size());
 418
 419     forAll (result, i)
 420     {
 421         result[i].transfer(candidateMasterNeighbors[i].shrink());
 422     }
 423 }
 424
 425 // This algorithm find the faces in proximity of another face based
 426 // on the face BB (Bounding Box) and an octree of bounding boxes.
 427 template<class MasterPatch, class SlavePatch>
 428 void GGIInterpolation<MasterPatch, SlavePatch>::findNeighboursBBOctree
 429 (
 430     labelListList& result
 431 ) const
 432 {
 433     List<DynamicList<label> > candidateMasterNeighbors(masterPatch_.size());
 434
 435     // Initialize the list of master patch faces bounding box
 436     treeBoundBoxList lmasterFaceBB(masterPatch_.size());
 437
 438     forAll (masterPatch_, faceMi)
 439     {
 440         pointField facePoints
 441         (
 442             masterPatch_[faceMi].points(masterPatch_.points())
 443         );
 444
 445         // Construct face BB with an extension of face span defined by the
 446         //  global tolerance factor faceBoundBoxExtendSpanFraction_
 447         // (1% by default)
 448         treeBoundBox bbFaceMaster(facePoints);
 449
 450         lmasterFaceBB[faceMi] =
 451             bbFaceMaster.extend(faceBoundBoxExtendSpanFraction_);
 452     }
 453
 454     // Initialize the list of slave patch faces bounding box
 455     treeBoundBoxList lslaveFaceBB(slavePatch_.size());
 456
 457     forAll (slavePatch_, faceSi)
 458     {
 459         pointField facePoints
 460         (
 461             slavePatch_[faceSi].points(slavePatch_.points())
 462         );
 463
 464         // possible transformation and separation for cyclic patches
 465         if (doTransform())
 466         {
 467             if (forwardT_.size() == 1)
 468             {
 469                 transform(facePoints, forwardT_[0], facePoints);
 470             }
 471             else
 472             {
 473                 transform(facePoints, forwardT_[faceSi], facePoints);
 474             }
 475         }
 476
 477         if (doSeparation())
 478         {
 479             if (forwardSep_.size() == 1)
 480             {
 481                 facePoints += forwardSep_[0];
 482             }
 483             else
 484             {
 485                 facePoints += forwardSep_[faceSi];
 486             }
 487         }
 488
 489         // Construct face BB with an extension of face span defined by the
 490         //  global tolerance factor faceBoundBoxExtendSpanFraction_
 491         // (1% by default)
 492         treeBoundBox bbFaceSlave(facePoints);
 493
 494         lslaveFaceBB[faceSi] =
 495             bbFaceSlave.extend(faceBoundBoxExtendSpanFraction_);
 496     }
 497
 498     // Create the slave octreeData, using the boundBox flavor
 499     octreeDataBoundBox slaveDataBB(lslaveFaceBB);
 500
 501     // Overall slave patch BB
 502     treeBoundBox slaveOverallBB(slavePatch_.points());
 503
 504     // Create the slave patch octree
 505
 506     octree<octreeDataBoundBox> slavePatchOctree
 507     (
 508         slaveOverallBB,              // overall search domain
 509         slaveDataBB,
 510         octreeSearchMinNLevel_,      // min number of levels
 511         octreeSearchMaxLeafRatio_,   // max avg. size of leaves
 512         octreeSearchMaxShapeRatio_   // max avg. duplicity.
 513     );
 514
 515     const vectorField& masterFaceNormals = masterPatch_.faceNormals();
 516     vectorField slaveNormals = slavePatch_.faceNormals();
 517
 518     // Transform slave normals to master plane if needed
 519     if (doTransform())
 520     {
 521         if (forwardT_.size() == 1)
 522         {
 523             transform(slaveNormals, forwardT_[0], slaveNormals);
 524         }
 525         else
 526         {
 527             transform(slaveNormals, forwardT_, slaveNormals);
 528         }
 529     }
 530
 531     // Visit each master patch face BB and find the potential neighbours
 532     // using the slave patch octree
 533
 534     forAll (lmasterFaceBB, faceMi)
 535     {
 536         // List of candidate neighbours
 537         labelList overlappedFaces  =
 538             slavePatchOctree.findBox(lmasterFaceBB[faceMi]);
 539
 540         forAll (overlappedFaces, ovFi)
 541         {
 542             label faceSi = overlappedFaces[ovFi];
 543
 544             // Compute and verify featureCos between the two face normals
 545             // before adding to the list of candidates
 546             scalar featureCos =
 547                 masterFaceNormals[faceMi] & slaveNormals[faceSi];
 548
 549             if (mag(featureCos) > featureCosTol_)
 550             {
 551                 candidateMasterNeighbors[faceMi].append(faceSi);
 552             }
 553         }
 554     }
 555
 556     // Repack the list
 557     result.setSize(masterPatch_.size());
 558
 559     forAll (result, i)
 560     {
 561         result[i].transfer(candidateMasterNeighbors[i].shrink());
 562     }
 563
 564     return;
 565 }
 566
 567
 568
 569 // Projects a list of points onto a plane located at planeOrig,
 570 // oriented along planeNormal.  Return the projected points in a
 571 // pointField, and the normal distance of each points from the
 572 // projection plane
 573 template<class MasterPatch, class SlavePatch>
 574 tmp<pointField> GGIInterpolation<MasterPatch, SlavePatch>::projectPointsOnPlane
 575 (
 576     const pointField& lpoints,
 577     const vector& planeOrig,
 578     const vector& planeDirection,
 579     scalarField& distanceProjection
 580 ) const
 581 {
 582     tmp<pointField> tprojectedPoints(new pointField(lpoints.size()));
 583     pointField& projectedPoints = tprojectedPoints();
 584
 585     vector normalVector = planeDirection/(mag(planeDirection) + VSMALL);
 586
 587     scalarField dist(lpoints.size(), 0.0);
 588
 589     if (lpoints.size() > 3 && mag(normalVector) > SMALL)
 590     {
 591         // Construct the plane
 592         plane projectionPlane(planeOrig, normalVector);
 593
 594         forAll (lpoints, pointI)
 595         {
 596             projectedPoints[pointI] =
 597                 projectionPlane.nearestPoint(lpoints[pointI]);
 598
 599             dist[pointI] = projectionPlane.distance(lpoints[pointI]);
 600         }
 601     }
 602     else
 603     {
 604           // Triangle... nothing to project, the points are already
 605           // located on the right plane
 606         projectedPoints = lpoints;
 607     }
 608
 609     // Transfer the projection distance
 610     distanceProjection = dist;
 611
 612     return tprojectedPoints;
 613 }
 614
 615
 616 // Compute an orthonormal basis (u, v, w) where: w is aligned on the
 617 // normalVector u is the direction from normalVectorCentre to the most
 618 // distant point in the list pointsOnPlane v = w^u
 619 //
 620 // Everything is normalized
 621 template<class MasterPatch, class SlavePatch>
 622 typename GGIInterpolation<MasterPatch, SlavePatch>::orthoNormalBasis
 623 GGIInterpolation<MasterPatch, SlavePatch>::computeOrthonormalBasis
 624 (
 625     const vector& normalVectorCentre,
 626     const vector& normalVector,
 627     const pointField& pointsOnPlane
 628 ) const
 629 {
 630     // The orthonormal basis uvw
 631     orthoNormalBasis uvw;
 632
 633     vector u = pTraits<vector>::zero;
 634     vector v = pTraits<vector>::zero;
 635     vector w = normalVector;
 636
 637     // Normalized w
 638     w /= mag(w) + VSMALL;
 639
 640     // Find the projected point from the master face that is the most distant
 641     scalar longestDistanceFromCenter = 0;
 642
 643     forAll (pointsOnPlane, pointI)
 644     {
 645         vector delta = pointsOnPlane[pointI] - normalVectorCentre;
 646
 647         if (mag(delta) > longestDistanceFromCenter)
 648         {
 649             longestDistanceFromCenter = mag(delta);
 650             u = delta; // This is our next candidate for the u direction
 651         }
 652     }
 653
 654     // Normalized u
 655     u /= mag(u) + VSMALL;
 656
 657     // Compute v from u and w
 658     v = w ^ u;    // v = w^u;
 659
 660     // We got ourselves a new orthonormal basis to play with
 661     uvw[0] = u;
 662     uvw[1] = v;
 663     uvw[2] = w;
 664
 665     return uvw;
 666 }
 667
 668
 669 // Projection of 3D points onto a 2D UV plane defined by an
 670 // orthonormal basis We project onto the uv plane. Could be
 671 // parametrized if we ever need to project onto other uvw planes
 672 template<class MasterPatch, class SlavePatch>
 673 List<point2D> GGIInterpolation<MasterPatch, SlavePatch>::projectPoints3Dto2D
 674 (
 675     const orthoNormalBasis& orthoBase,
 676     const vector& orthoBaseOffset,
 677     const pointField& pointsIn3D,
 678     scalarField& distanceProjection
 679 ) const
 680 {
 681     List<point2D> pointsIn2D(pointsIn3D.size());
 682     scalarField  dist(pointsIn3D.size(), 0.0);
 683
 684     pointField pointsIn3DTranslated = pointsIn3D - orthoBaseOffset;
 685
 686     // Project onto the uv plane. The distance from the plane is
 687     // computed with w
 688     forAll (pointsIn3D, pointsI)
 689     {
 690         // u component
 691         pointsIn2D[pointsI][0] =
 692             pointsIn3DTranslated[pointsI] & orthoBase[0];
 693
 694         // v component
 695         pointsIn2D[pointsI][1] =
 696             pointsIn3DTranslated[pointsI] & orthoBase[1];
 697
 698         // w component = error above projection plane
 699         dist[pointsI] = pointsIn3DTranslated[pointsI] & orthoBase[2];
 700     }
 701
 702     distanceProjection = dist;
 703
 704     return pointsIn2D;
 705 }
 706
 707
 708 // Projection of 2D points onto a 1D normalized direction vector
 709 template<class MasterPatch, class SlavePatch>
 710 scalarField GGIInterpolation<MasterPatch, SlavePatch>::projectPoints2Dto1D
 711 (
 712     const vector2D&      normalizedProjectionDir,
 713     const point2D&       normalizedProjectionDirOffset,
 714     const List<point2D>& lPoints2D
 715     ) const
 716 {
 717     scalarField pointsIn1D(lPoints2D.size()); // Return values.
 718
 719     forAll (lPoints2D, ptsI)
 720     {
 721         pointsIn1D[ptsI] =
 722             normalizedProjectionDir & (lPoints2D[ptsI]
 723           - normalizedProjectionDirOffset);
 724     }
 725
 726     return pointsIn1D;
 727 }
 728
 729
 730 // Polygon overlap or polygon collision detection using the Separating
 731 // Axes Theorem.  We expect this algorithm to possibly call himself a
 732 // second time: code needs to be re-entrant!!!
 733 //
 734 // First time it is called, we test polygon2 agains each edges of
 735 // polygon1 If we still detect an overlap, then the algorith will call
 736 // itself to test polygon1 against polygon2 in order to find a possibe
 737 // separating Axes.  The flag firstCall will be used for detecting if
 738 // we are on the first or second call, so we can exit properly after
 739 // the second call.
 740 template<class MasterPatch, class SlavePatch>
 741 bool GGIInterpolation<MasterPatch, SlavePatch>::detect2dPolygonsOverlap
 742 (
 743     const List<point2D>& poly1,
 744     const List<point2D>& poly2,
 745     const scalar& tolFactor,
 746     const bool firstCall
 747 ) const
 748 {
 749     bool isOverlapping = true;
 750
 751     for
 752     (
 753         label istart1 = 0, iend1 = poly1.size() - 1;
 754         istart1 < poly1.size();
 755         iend1 = istart1, istart1++
 756     )
 757     {
 758         // Reference edge from polygon1
 759         point2D curEdge1 = poly1[istart1] - poly1[iend1];
 760
 761         point2D origEdge1 = poly1[iend1];
 762
 763         // normalPerpDirection1 & curEdge1 == 0
 764         vector2D normalPerpDirection1(curEdge1[1], - curEdge1[0]);
 765
 766         // Normalize
 767         normalPerpDirection1 /= mag(normalPerpDirection1) + VSMALL;
 768
 769         // Project poly1 onto normalPerpDirection1
 770         scalarField poly1PointsProjection = projectPoints2Dto1D
 771         (
 772             normalPerpDirection1,
 773             origEdge1,
 774             poly1
 775         );
 776
 777         // Project poly2 onto normalPerpDirection1
 778         scalarField poly2PointsProjection = projectPoints2Dto1D
 779         (
 780             normalPerpDirection1,
 781             origEdge1,
 782             poly2
 783         );
 784
 785         // Grab range of the projected polygons
 786         scalar p1_min = Foam::min(poly1PointsProjection);
 787         scalar p1_max = Foam::max(poly1PointsProjection);
 788         scalar p2_min = Foam::min(poly2PointsProjection);
 789         scalar p2_max = Foam::max(poly2PointsProjection);
 790
 791         // Here are all the possible situations for the range overlap
 792         // detection, and a simple test that discriminates them all
 793         //
 794         //
 795         //    P1 -------------                       p1_min + epsilon < p2_max  &&  p2_min + epsilon < p1_max  == true
 796         //    P2 -------------
 797         //
 798         //    P2 -------------                       p1_min + epsilon < p2_max  &&  p2_min + epsilon < p1_max  == true
 799         //    P1 -------------
 800         //
 801         //
 802         //    P1  -------------                      p1_min + epsilon < p2_max  &&  p2_min + epsilon < p1_max  == true
 803         //         P2  ----
 804         //
 805         //         P1  ----
 806         //    P2  -------------                      p1_min + epsilon < p2_max  &&  p2_min + epsilon < p1_max  == true
 807         //
 808         //
 809         //    P1  -------------                      p1_min + epsilon < p2_max  &&  p2_min + epsilon < p1_max  == true
 810         //              P2  ---------
 811         //
 812         //    P2  -------------                      p1_min + epsilon < p2_max  &&  p2_min + epsilon < p1_max  == true
 813         //              P1  ---------
 814         //
 815         //
 816         //                        P1 -------------   p1_min + epsilon < p2_max  &&  p2_min + epsilon < p1_max  == false
 817         //    P2  -------------
 818         //
 819         //
 820         //                        P2 -------------   p1_min + epsilon < p2_max  &&  p2_min + epsilon < p1_max  == false
 821         //    P1  -------------
 822         //
 823         //
 824         //
 825         //    P1  -------------                      p1_min + epsilon < p2_max  &&  p2_min + epsilon < p1_max  == false
 826         //                     ------------- P2
 827         //
 828         //
 829         //    P2  -------------                      p1_min + epsilon < p2_max  &&  p2_min + epsilon < p1_max  == false
 830         //                     ------------- P1
 831         //
 832         //
 833         //
 834         // What value should we give to epsilon???
 835         //
 836         // The intervals [p1_min, p1_max] and [p2_min, p2_max] are
 837         // basically the dimension of one side of the BB for a given
 838         // polygon, for a given "orientation" of the polygon.
 839         //
 840         // This means that epsilon^2 is roughly the size of the minimal
 841         // surface area intersecting the 2 polygons that one accept to
 842         // discard.
 843         //
 844         // So, if for instance we fix epsilon = 10e-3 * min
 845         // (range_of_polygon1, range_of_polygon2), this means that we
 846         // accept to discard an intersecting area roughly 10e-6 times the
 847         // surface of the smallest polygon, so 1 PPM.
 848         //
 849         // Which means also that our GGI weighting factors will never be
 850         //  smaller than roughly 10e-6, because this is the fraction of
 851         // the intersection surface area we choose to discard.
 852
 853         scalar _epsilon = tolFactor*
 854             Foam::min
 855             (
 856                 mag(p1_max - p1_min),
 857                 mag(p2_max - p2_min)
 858             );
 859
 860         // Evaluate the presence or not of an overlapping
 861         if
 862         (
 863             !((p1_min + _epsilon < p2_max)
 864           && (p2_min + _epsilon < p1_max))
 865         )
 866         {
 867             isOverlapping = false; // We have found a separating axis!
 868             break;
 869         }
 870     }
 871
 872     if (isOverlapping && firstCall)
 873     {
 874         // We have not found any separating axes by exploring from
 875         // poly1, let's switch by exploring from poly2 instead
 876         isOverlapping =
 877             detect2dPolygonsOverlap(poly2, poly1, tolFactor, false);
 878     }
 879
 880     return isOverlapping;
 881 }
 882
 883
 884 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 885
 886 } // End namespace Foam
 887
 888 // ************************************************************************* //