src/foam/interpolations/GGIInterpolation/GGIInterpolationWeights.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     Mass-conservative face interpolation of face data between two
  26     primitivePatches
  27
  28 Author
  29     Hrvoje Jasak, Wikki Ltd.  All rights reserved
  30
  31 Modification by:
  32     Martin Beaudoin, Hydro-Quebec, (2008)
  33
  34 \*---------------------------------------------------------------------------*/
  35
  36 #include "GGIInterpolationTemplate.H"
  37 #include "objectHit.H"
  38 #include "boolList.H"
  39 #include "DynamicList.H"
  40 #include "dimensionedConstants.H"
  41
  42 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  43
  44 namespace Foam
  45 {
  46
  47 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  48
  49 template<class MasterPatch, class SlavePatch>
  50 const Foam::debug::tolerancesSwitch
  51 GGIInterpolation<MasterPatch, SlavePatch>::areaErrorTol_
  52 (
  53     "GGIAreaErrorTol",
  54     1.0e-8,
  55     "Minimum GGI face to face intersection area. The smallest accepted GGI weighting factor."
  56 );
  57
  58 template<class MasterPatch, class SlavePatch>
  59 const Foam::debug::tolerancesSwitch
  60 GGIInterpolation<MasterPatch, SlavePatch>::featureCosTol_
  61 (
  62     "GGIFeatureCosTol",
  63     0.8,
  64     "Minimum cosine value between 2 GGI patch neighbouring facet normals."
  65 );
  66
  67 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  68
  69 template<class MasterPatch, class SlavePatch>
  70 void GGIInterpolation<MasterPatch, SlavePatch>::calcAddressing() const
  71 {
  72     if
  73     (
  74         masterAddrPtr_
  75      || masterWeightsPtr_
  76      || slaveAddrPtr_
  77      || slaveWeightsPtr_
  78      || uncoveredMasterAddrPtr_
  79      || uncoveredSlaveAddrPtr_
  80     )
  81     {
  82         FatalErrorIn
  83         (
  84             "void GGIInterpolation<MasterPatch, SlavePatch>::"
  85             "calcAddressing() const"
  86         )   << "Addressing already calculated"
  87             << abort(FatalError);
  88     }
  89
  90     if (debug)
  91     {
  92         InfoIn
  93         (
  94             "void GGIInterpolation<MasterPatch, SlavePatch>::"
  95             "calcAddressing() const"
  96         )   << "Evaluation of GGI weighting factors:" << endl;
  97     }
  98
  99     // Create the dynamic lists to hold the addressing
 100
 101     // The final master/slave list, after filtering out the "false" neighbours
 102     List<DynamicList<label, 8> > masterNeighbors(masterPatch_.size());
 103     List<DynamicList<label, 8> > slaveNeighbors(slavePatch_.size());
 104
 105     // The weights
 106     List<DynamicList<scalar, 8> > masterNeighborsWeights(masterPatch_.size());
 107     List<DynamicList<scalar, 8> > slaveNeighborsWeights(slavePatch_.size());
 108
 109     // First, find a rough estimate of each slave and master facet
 110     // neighborhood by filtering out all the faces located outside of
 111     // an Axis-Aligned Bounding Box (AABB).  Warning: This algorithm
 112     // is based on the evaluation of AABB boxes, which is pretty fast;
 113     // but still the complexity of the algorithm is n^2, wich is
 114     // pretty bad for GGI patches composed of 100,000 of facets...  So
 115     // here is the place where we could certainly gain major speedup
 116     // for larger meshes.
 117
 118     // The candidates master neighbours
 119     // Choice of algorithm:
 120     // 1) Axis-aligned bounding box
 121     // 2) Octree search with bounding box
 122     // 3) 3-D vector distance
 123     // 4) n-Squared search
 124     labelListList candidateMasterNeighbors;
 125
 126     if (reject_ == AABB)
 127     {
 128          findNeighboursAABB(candidateMasterNeighbors);
 129     }
 130     else if (reject_ == BB_OCTREE)
 131     {
 132          findNeighboursBBOctree(candidateMasterNeighbors);
 133     }
 134     else if (reject_ == THREE_D_DISTANCE)
 135     {
 136          findNeighbours3D(candidateMasterNeighbors);
 137     }
 138     else if (reject_ == N_SQUARED)
 139     {
 140         candidateMasterNeighbors.setSize(masterPatch_.size());
 141
 142         // Mark N-squared search
 143         labelList nSquaredList(slavePatch_.size());
 144         forAll (nSquaredList, i)
 145         {
 146             nSquaredList[i] = i;
 147         }
 148
 149         forAll (candidateMasterNeighbors, j)
 150         {
 151             candidateMasterNeighbors[j] = nSquaredList;
 152         }
 153     }
 154
 155     // Next, we move to the 2D world.  We project each slave and
 156     // master face onto a local plane defined by the master face
 157     // normal.  We filter out a few false neighbors using the
 158     // Separating Axes Theorem
 159
 160     // It is in this local plane that we will refine our list of
 161     // neighbors.  So for a given a neighbor face, we need as many
 162     // projections as there are neighbors closeby.
 163
 164     const pointField& masterPatchPoints = masterPatch_.points();
 165     const vectorField masterPatchNormals = masterPatch_.faceNormals();
 166
 167     // Store the polygon made by projecting the face points onto the
 168     // face normal
 169     // The master faces polygons
 170     List<pointField>  masterFace2DPolygon(masterPatch_.size());
 171
 172     // Tolerance factor for the Separation of Axes Theorem == distErrorTol_
 173
 174     forAll (masterPatch_, faceMi)
 175     {
 176         // First, we make sure that all the master faces points are
 177         // recomputed onto the 2D plane defined by the master faces
 178         // normals.
 179         // For triangles, this is useless, but for N-gons
 180         // with more than 3 points, this is essential.
 181         // The intersection between the master and slave faces will be
 182         // done in these 2D reference frames
 183
 184         // A few basic information to keep close-by
 185         vector currentMasterFaceNormal = masterPatchNormals[faceMi];
 186         vector currentMasterFaceCentre =
 187             masterPatch_[faceMi].centre(masterPatchPoints);
 188
 189         scalarField facePolygonErrorProjection;
 190
 191         // Project the master faces points onto the normal face plane to
 192         // form a flattened polygon
 193         masterFace2DPolygon[faceMi] =
 194             projectPointsOnPlane
 195             (
 196                 masterPatch_[faceMi].points(masterPatchPoints),
 197                 currentMasterFaceCentre,
 198                 currentMasterFaceNormal,
 199                 facePolygonErrorProjection
 200             );
 201
 202         // Next we compute an orthonormal basis (u, v, w) aligned with
 203         // the face normal for doing the 3D to 2D projection.
 204         //
 205         // "w" is aligned on the face normal.  We need to select a "u"
 206         // direction, it can be anything as long as it lays on the
 207         // projection plane.  We chose to use the direction from the
 208         // master face center to the most distant projected master face
 209         // point on the plane.  Finally, we get "v" by evaluating the
 210         // cross-product w^u = v.  And we make sure that u, v, and w are
 211         // normalized.
 212         //
 213         //
 214         //                                                       u  =  vector from face center to most distant projected master face point.
 215         //                                               /       .
 216         //           ^y                                / |       .      .w = normal to master face
 217         //           |                               /   |       .    .
 218         //           |                             /     |       .  .
 219         //           |                            |      |       .
 220         //           |                            |      /        .
 221         //           |                            |    /           .
 222         //           |                            |  /              .
 223         //           ---------> x                 |/                 .
 224         //          /                                                 v = w^u
 225         //         /
 226         //        /
 227         //       z
 228         //
 229         //
 230
 231         orthoNormalBasis uvw =
 232             computeOrthonormalBasis
 233             (
 234                 currentMasterFaceCentre,
 235                 currentMasterFaceNormal,
 236                 masterFace2DPolygon[faceMi]
 237             );
 238
 239         // Recompute the master polygon into this orthoNormalBasis
 240         // We should only see a rotation along the normal of the face here
 241         List<point2D> masterPointsInUV;
 242         scalarField masterErrorProjectionAlongW;
 243
 244         masterPointsInUV =
 245             projectPoints3Dto2D
 246             (
 247                 uvw,
 248                 currentMasterFaceCentre,
 249                 masterFace2DPolygon[faceMi],
 250                 masterErrorProjectionAlongW   // Should be at zero all the way
 251             );
 252
 253         // Compute the surface area of the polygon;
 254         // We need this for computing the weighting factors
 255         scalar surfaceAreaMasterPointsInUV = area2D(masterPointsInUV);
 256
 257         // Check if polygon is CW.. Should not, it should be CCW; but
 258         // better and cheaper to check here
 259         if (surfaceAreaMasterPointsInUV < 0.0)
 260         {
 261             reverse(masterPointsInUV);
 262             surfaceAreaMasterPointsInUV = -surfaceAreaMasterPointsInUV;
 263
 264             // Just generate a warning until we can verify this is a non issue
 265             InfoIn
 266             (
 267                 "void GGIInterpolation<MasterPatch, SlavePatch>::"
 268                 "calcAddressing()"
 269             )   << "The master projected polygon was CW instead of CCW.  "
 270                 << "This is strange..."  << endl;
 271         }
 272
 273         // Next, project the candidate master neighbours faces points
 274         // onto the same plane using the new orthonormal basis
 275         const labelList& curCMN = candidateMasterNeighbors[faceMi];
 276
 277         forAll (curCMN, neighbI)
 278         {
 279             // For each points, compute the dot product with u,v,w.  The
 280             // [u,v] component will gives us the 2D cordinates we are
 281             // looking for for doing the 2D intersection The w component
 282             // is basically the projection error normal to the projection
 283             // plane
 284
 285             // NB: this polygon is most certainly CW w/r to the uvw
 286             // axis because of the way the normals are oriented on
 287             // each side of the GGI interface... We will switch the
 288             // polygon to CCW in due time...
 289             List<point2D> neighbPointsInUV;
 290             scalarField neighbErrorProjectionAlongW;
 291
 292             // We use the xyz points directly, with a possible transformation
 293             pointField curSlaveFacePoints =
 294                 slavePatch_[curCMN[neighbI]].points(slavePatch_.points());
 295
 296             if (doTransform())
 297             {
 298                 // Transform points to master plane
 299                 if (forwardT_.size() == 1)
 300                 {
 301                     transform
 302                     (
 303                         curSlaveFacePoints,
 304                         forwardT_[0],
 305                         curSlaveFacePoints
 306                     );
 307                 }
 308                 else
 309                 {
 310                     transform
 311                     (
 312                         curSlaveFacePoints,
 313                         forwardT_[curCMN[neighbI]],
 314                         curSlaveFacePoints
 315                     );
 316                 }
 317             }
 318
 319             // Apply the translation offset in order to keep the
 320             // neighbErrorProjectionAlongW values to a minimum
 321             if (doSeparation())
 322             {
 323                 if (forwardSep_.size() == 1)
 324                 {
 325                     curSlaveFacePoints += forwardSep_[0];
 326                 }
 327                 else
 328                 {
 329                     curSlaveFacePoints += forwardSep_[curCMN[neighbI]];
 330                 }
 331             }
 332
 333             neighbPointsInUV =
 334                 projectPoints3Dto2D
 335                 (
 336                     uvw,
 337                     currentMasterFaceCentre,
 338                     curSlaveFacePoints,
 339                     neighbErrorProjectionAlongW
 340                 );
 341
 342             // We are now ready to filter out the "bad" neighbours.
 343             // For this, we will apply the Separating Axes Theorem
 344             // http://en.wikipedia.org/wiki/Separating_axis_theorem.
 345
 346             // This will be the second and last quick reject test.
 347             // We will use the 2D projected points for both the master
 348             // patch and its neighbour candidates
 349             if
 350             (
 351                 detect2dPolygonsOverlap
 352                 (
 353                     masterPointsInUV,
 354                     neighbPointsInUV,
 355                     sqrt(areaErrorTol_()) // distErrorTol
 356                 )
 357             )
 358             {
 359                 // We have an overlap between the master face and this
 360                 // neighbor face.
 361                 label faceMaster = faceMi;
 362                 label faceSlave  = curCMN[neighbI];
 363
 364                 // Compute the surface area of the neighbour polygon;
 365                 // We need this for computing the weighting factors
 366                 scalar surfaceAreaNeighbPointsInUV = area2D(neighbPointsInUV);
 367
 368                 // Check for CW polygons. It most certainly is, and
 369                 // the polygon intersection algorithms are expecting
 370                 // to work with CCW point ordering for the polygons
 371                 if (surfaceAreaNeighbPointsInUV < 0.0)
 372                 {
 373                     reverse(neighbPointsInUV);
 374                     surfaceAreaNeighbPointsInUV = -surfaceAreaNeighbPointsInUV;
 375                 }
 376
 377
 378                 // We compute the intersection area using the
 379                 // Sutherland-Hodgman algorithm.  Of course, if the
 380                 // intersection area is 0, that would constitute the last and
 381                 // final reject test, but it would also be an indication that
 382                 // our 2 previous rejection tests are a bit laxed...  or that
 383                 // maybe we are in presence of concave polygons....
 384                 scalar intersectionArea =
 385                     polygonIntersection
 386                     (
 387                         masterPointsInUV,
 388                         neighbPointsInUV
 389                     );
 390
 391                 if (intersectionArea > VSMALL) // Or > areaErrorTol_ ???
 392                 {
 393                     // We compute the GGI weights based on this
 394                     // intersection area, and on the individual face
 395                     // area on each side of the GGI.
 396
 397                     // Since all the intersection have been computed
 398                     // in the projected UV space we need to compute
 399                     // the weights using the surface area from the
 400                     // faces projection as well. That way, we make
 401                     // sure all our factors will sum up to 1.0.
 402
 403                     masterNeighbors[faceMaster].append(faceSlave);
 404                     slaveNeighbors[faceSlave].append(faceMaster);
 405
 406                     masterNeighborsWeights[faceMaster].append
 407                     (
 408                         intersectionArea/surfaceAreaMasterPointsInUV
 409                     );
 410
 411                     slaveNeighborsWeights[faceSlave].append
 412                     (
 413                         intersectionArea/surfaceAreaNeighbPointsInUV
 414                     );
 415                 }
 416                 else
 417                 {
 418                     WarningIn
 419                     (
 420                         "GGIInterpolation<MasterPatch, SlavePatch>::"
 421                         "calcAddressing()"
 422                     )   << "polygonIntersection is returning a "
 423                         << "zero surface area between " << nl
 424                         << "     Master face: " << faceMi
 425                         << " and Neighbour face: " << curCMN[neighbI]
 426                         << " intersection area = " << intersectionArea << nl
 427                         << "Please check the two quick-check algorithms for "
 428                         << "GGIInterpolation.  Something is  missing." << endl;
 429                 }
 430             }
 431         }
 432
 433         // We went through all the possible neighbors for this face.
 434     }
 435
 436
 437     // Allocate the member attributes and pack addressing
 438     masterAddrPtr_ = new labelListList(masterPatch_.size());
 439     labelListList& ma  = *masterAddrPtr_;
 440
 441     masterWeightsPtr_ = new scalarListList(masterPatch_.size());
 442     scalarListList& maW = *masterWeightsPtr_;
 443
 444     forAll (ma, mfI)
 445     {
 446         ma[mfI].transfer(masterNeighbors[mfI].shrink());
 447         maW[mfI].transfer(masterNeighborsWeights[mfI].shrink());
 448     }
 449
 450     slaveAddrPtr_ = new labelListList(slavePatch_.size());
 451     labelListList& sa = *slaveAddrPtr_;
 452
 453     slaveWeightsPtr_ = new scalarListList(slavePatch_.size());
 454     scalarListList& saW = *slaveWeightsPtr_;
 455
 456     forAll (sa, sfI)
 457     {
 458         sa[sfI].transfer(slaveNeighbors[sfI].shrink());
 459         saW[sfI].transfer(slaveNeighborsWeights[sfI].shrink());
 460     }
 461
 462     // Now that the neighbourhood is known, let's go hunting for
 463     // non-overlapping faces
 464     uncoveredMasterAddrPtr_ =
 465         new labelList
 466         (
 467             findNonOverlappingFaces(maW, masterNonOverlapFaceTol_)
 468         );
 469
 470     uncoveredSlaveAddrPtr_ =
 471         new labelList
 472         (
 473             findNonOverlappingFaces(saW, slaveNonOverlapFaceTol_)
 474         );
 475
 476     // Rescaling the weighting factors so they will sum up to 1.0
 477     // See the comment for the method ::rescaleWeightingFactors() for
 478     // more information.  By default, we always rescale.  But for some
 479     // special kind of GGI interpolation, like the mixingPlaneGGI,
 480     // then we need the brute values, so no rescaling in that
 481     // case. Hence the little flag rescaleGGIWeightingFactors_
 482
 483     if (rescaleGGIWeightingFactors_)
 484     {
 485         rescaleWeightingFactors();
 486     }
 487 }
 488
 489
 490 // Rescaling the weighting factors so they will sum up to 1.0 This is
 491 // necessary for the slave weighting factors because intersection with
 492 // master neighbours are usually computed from different projection
 493 // plane, so the weighting factor don't quite sum up to 1.0 For the
 494 // slave weighting factors, we are usually talking of delta of the
 495 // order of 10e-6 here.
 496 //
 497 // For the master weighting factor, this is another story. We truly
 498 // expect that the master weighting factors will exactly sum up to 1.0
 499 // if all the neighbours are properly identified.
 500 //
 501 // However, for concentric circular geometry, if the circumferantial
 502 // resolution is too coarse, we will end up with some part of the face
 503 // surface that are not taken into account because they do not
 504 // physically overlap any neighbours.  For example, think of 2
 505 // concentric circular patches, slightly rotated one relatively to the
 506 // other.  A good case: the ercoftac conical diffuser, Case0...  GGI
 507 // located right between the cylindrical and conical parts, rotate the
 508 // cylindrical by 15 degrees.  For this case, we will need to devise a
 509 // decent strategy in order to intelligently take care of these
 510 // "missing weights"
 511 //
 512 // The purpose of the ::rescaleWeightingFactors() method is mainly for
 513 // this.
 514 template<class MasterPatch, class SlavePatch>
 515 void GGIInterpolation<MasterPatch, SlavePatch>::rescaleWeightingFactors() const
 516 {
 517     scalarListList& maW = *masterWeightsPtr_;
 518     scalarListList& saW = *slaveWeightsPtr_;
 519
 520     // Memorize the largest correction needed in order to provide some
 521     // basic info to the user
 522     scalar largestSWC = 0;
 523     scalar sumSWC = 0;
 524     scalar curSWC = 0;
 525
 526     scalar largestMWC = 0;
 527     scalar sumMWC = 0;
 528     scalar curMWC = 0;
 529
 530     // Rescaling the slave weights
 531     if (debug)
 532     {
 533         if
 534         (
 535             uncoveredMasterFaces().size() > 0
 536          || uncoveredSlaveFaces().size() > 0
 537         )
 538         {
 539             InfoIn
 540             (
 541                 "void GGIInterpolation<MasterPatch, SlavePatch>::"
 542                 "rescaleWeightingFactors() const"
 543             )   << "Uncovered faces found.  On master: "
 544                 << uncoveredMasterFaces().size()
 545                 << " on slave: " << uncoveredSlaveFaces().size() << endl;
 546         }
 547       }
 548
 549     forAll (saW, saWi)
 550     {
 551         scalar slaveWeightSum = Foam::sum(saW[saWi]);
 552
 553         if (saW[saWi].size() > 0)
 554         {
 555             saW[saWi] = saW[saWi]/slaveWeightSum;
 556
 557             // Some book-keeping
 558             curSWC = mag(1.0 - slaveWeightSum);
 559             largestSWC = Foam::max(largestSWC, curSWC);
 560
 561             sumSWC += curSWC;
 562         }
 563     }
 564
 565     // Rescaling the master weights
 566     forAll (maW, maWi)
 567     {
 568         scalar masterWeightSum = Foam::sum(maW[maWi]);
 569
 570         if (maW[maWi].size() > 0)
 571         {
 572             maW[maWi] = maW[maWi]/masterWeightSum;
 573
 574             // Some book-keeping
 575             curMWC = mag(1.0 - masterWeightSum);
 576             largestMWC = Foam::max(largestMWC, curMWC);
 577
 578             sumMWC += curMWC;
 579         }
 580     }
 581
 582     if (debug)
 583     {
 584         if (saW.size() > 0 && maW.size() > 0)
 585         {
 586             Info<< "  Largest slave weighting factor correction : "
 587                 << largestSWC
 588                 << " average: " << sumSWC/saW.size() << nl
 589                 << "  Largest master weighting factor correction: "
 590                 << largestMWC
 591                 << " average: " << sumMWC/maW.size() << endl;
 592         }
 593     }
 594 }
 595
 596
 597 // Find non-overlapping faces from both master and slave patches
 598 // The default non-overlapping criteria is total absence of neighbours.
 599 // Later on, ths criteria will be based on minimum surface intersection, or
 600 // minimum weight factor
 601 template<class MasterPatch, class SlavePatch>
 602 tmp<labelField>
 603 GGIInterpolation<MasterPatch, SlavePatch>::findNonOverlappingFaces
 604 (
 605     const scalarListList& patchWeights,
 606     const scalar& nonOverlapFaceTol   //  = min sum of the neighbour weights
 607 ) const
 608 {
 609     tmp<labelField> tpatchFaceNonOverlapAddr(new labelField());
 610     labelField& patchFaceNonOverlapAddr = tpatchFaceNonOverlapAddr();
 611
 612     DynamicList<label, 64> patchFaceNonOverlap(patchWeights.size());
 613
 614     // Scan the list of patch weights, looking for empty lists
 615     forAll (patchWeights, paWi)
 616     {
 617         scalar sumWeightsFace = sum(patchWeights[paWi]);
 618
 619         if (sumWeightsFace <= nonOverlapFaceTol)
 620         {
 621             // Store local index.
 622             patchFaceNonOverlap.append(paWi);
 623         }
 624     }
 625
 626     if (patchFaceNonOverlap.size() > 0)
 627     {
 628         patchFaceNonOverlapAddr.transfer(patchFaceNonOverlap.shrink());
 629     }
 630
 631     if (debug)
 632     {
 633         InfoIn("GGIInterpolation::findNonOverlappingFaces")
 634             << "   : Found " << patchFaceNonOverlapAddr.size()
 635             << " non-overlapping faces for this GGI patch" << endl;
 636     }
 637
 638     return tpatchFaceNonOverlapAddr;
 639 }
 640
 641 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 642
 643 } // End namespace Foam
 644
 645 // ************************************************************************* //