src/foam/interpolations/GGIInterpolation/GGIInterpolationQuickRejectTests.C

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