src/foam/interpolations/MixingPlaneInterpolation/MixingPlaneProfile.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     MixingPlane interpolation functions
  26
  27 Author
  28     Martin Beaudoin, Hydro-Quebec, 2009.  All rights reserved
  29
  30 Contributor
  31     Hrvoje Jasak, Wikki Ltd.
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 namespace Foam
  36 {
  37
  38 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  39
  40 template<class MasterPatch, class SlavePatch>
  41 tmp<pointField>
  42 MixingPlaneInterpolation<MasterPatch, SlavePatch>::computeProfileFromHistograms
  43 (
  44     const profileHistogram& masterHisto,
  45     const profileHistogram& slaveHisto,
  46     const scalar halfSizeBin
  47 ) const
  48 {
  49     // Find min, max bounds
  50     scalar histoMinValue =
  51         Foam::min
  52         (
  53             masterHisto.begin()->first,
  54             slaveHisto.begin()->first
  55         );
  56
  57     scalar histoMaxValue =
  58         Foam::max
  59         (
  60             (--masterHisto.end())->first,
  61             (--slaveHisto.end())->first
  62         );
  63
  64     // Next, we compare both histograms, leap-frogging from one histogram
  65     // to the other, everytime jumping to the next largest value from
  66     // a given position
  67     std::list<point> leapFrogProfile(0);
  68
  69     scalar curRvalue = histoMinValue;
  70
  71     if (masterHisto.find(curRvalue) != masterHisto.end())
  72     {
  73         leapFrogProfile.push_back
  74         (
  75             *(masterHisto.find(curRvalue)->second.begin())
  76         );
  77     }
  78     else
  79     {
  80         leapFrogProfile.push_back
  81          (
  82              *(slaveHisto.find(curRvalue)->second.begin())
  83          );
  84     }
  85
  86     do
  87     {
  88         profileHistogram::const_iterator nextMaster =
  89             masterHisto.lower_bound(curRvalue + halfSizeBin - SMALL);
  90
  91         profileHistogram::const_iterator nextSlave  =
  92             slaveHisto.lower_bound(curRvalue + halfSizeBin - SMALL);
  93
  94         if
  95         (
  96             nextMaster == masterHisto.end()
  97          || nextSlave  == slaveHisto.end()
  98         )
  99         {
 100             // We are done
 101             if (curRvalue != histoMaxValue)
 102             {
 103                 if (masterHisto.find(histoMaxValue) != masterHisto.end())
 104                 {
 105                     leapFrogProfile.push_back
 106                     (
 107                         *(masterHisto.find(histoMaxValue)->second.begin())
 108                     );
 109                 }
 110                 else
 111                 {
 112                     leapFrogProfile.push_back
 113                     (
 114                         *(slaveHisto.find(histoMaxValue)->second.begin())
 115                     );
 116                 }
 117             }
 118             break;
 119         }
 120
 121         // Leap frog to the next largest delta
 122         // Please note that we are not taking into account the number of
 123         // values at the specific map index (the attribute called 'second'
 124         // of the map container)
 125         // So by leap-frogging this way, we might ran into the situation
 126         // where the next jump will bring us to a radius value shared by only
 127         // one patch face point, and we might oversee a slightly smaller
 128         // radius shared by 100s of points.
 129         // There is a definite advantage of sticking to existing radius values
 130         // for the profile because this will tend to keep the number of
 131         // GGI face neighbours low.
 132         // I am not sure if there is a nice solution to this problem.... MB
 133         //
 134         curRvalue = max(nextMaster->first, nextSlave->first);
 135
 136         if (masterHisto.find(curRvalue) != masterHisto.end())
 137         {
 138             leapFrogProfile.push_back
 139             (
 140                 *(masterHisto.find(curRvalue)->second.begin())
 141             );
 142         }
 143         else
 144         {
 145             leapFrogProfile.push_back
 146             (
 147                 *(slaveHisto.find(curRvalue)->second.begin())
 148             );
 149         }
 150     } while (curRvalue != histoMaxValue);
 151
 152     // Re-package the data into pointField
 153     tmp<pointField> tprofile(new pointField(leapFrogProfile.size()));
 154
 155     pointField& profile = tprofile();
 156
 157     label pI = 0;
 158
 159     forAllIter (std::list<point>, leapFrogProfile, lI)
 160     {
 161         profile[pI++] = *lI;
 162     }
 163
 164     return tprofile;
 165 }
 166
 167
 168 template<class MasterPatch, class SlavePatch>
 169 void
 170 MixingPlaneInterpolation<MasterPatch, SlavePatch>::updateProfileHistogram
 171 (
 172     profileHistogram& histo,
 173     const point& profileCoord,  // 3D point reference
 174     const direction dir,        // Sorting dimension 0: x, 1: y, 2: z
 175     const scalar halfSizeBin    // half size of min width for histogram bins
 176 ) const
 177 {
 178     bool foundNewBin = true;
 179
 180     scalar keyValue = profileCoord.component(dir);
 181
 182     forAllIter (profileHistogram, histo, histoI)
 183     {
 184         if
 185         (
 186             keyValue >= histoI->first - halfSizeBin
 187          && keyValue < histoI->first + halfSizeBin
 188         )
 189         {
 190             foundNewBin = false;
 191             histoI->second.push_back(profileCoord);
 192             break;
 193         }
 194     }
 195
 196     if (foundNewBin)
 197     {
 198         std::list<point> initValue;
 199
 200         initValue.push_back(profileCoord);
 201         histo.insert
 202         (
 203             std::pair<scalar, std::list<point> >(keyValue, initValue)
 204         );
 205     }
 206 }
 207
 208
 209 template<class MasterPatch, class SlavePatch>
 210 tmp<pointField>
 211 MixingPlaneInterpolation<MasterPatch, SlavePatch>::calcProfile() const
 212 {
 213     // First, transform patch points over one single global profile
 214     // into local coordinates
 215
 216     pointField masterGlobalProfile =
 217         cs_.localPosition(masterPatch_.localPoints());
 218
 219     pointField slaveGlobalProfile =
 220         cs_.localPosition(slavePatch_.localPoints());
 221
 222     // Collapse all points in the sweep direction
 223
 224     const direction sweepDir = sweepAxisSwitch();
 225
 226     masterGlobalProfile.replace(sweepDir, 0);
 227     slaveGlobalProfile.replace(sweepDir, 0);
 228
 229     if (debug)
 230     {
 231         InfoIn
 232         (
 233             "tmp<pointField> MixingPlaneInterpolation::calcProfile()"
 234         )   << "masterGlobalProfile: " << masterGlobalProfile << nl
 235             << "slaveGlobalProfile: " << slaveGlobalProfile << endl;
 236     }
 237
 238     // Find the smallest edge length from both patches.
 239     // This length will control the size of the bins for the histograms
 240     // we are about to build.
 241     scalar masterMinEdgeLength = GREAT;
 242     const edgeList& masterEdgeList = masterPatch_.edges();
 243
 244     forAll (masterEdgeList, mEi)
 245     {
 246         masterMinEdgeLength =
 247             Foam::min
 248             (
 249                 masterMinEdgeLength,
 250                 masterEdgeList[mEi].mag(masterPatch_.localPoints())
 251             );
 252     }
 253
 254     scalar slaveMinEdgeLength = GREAT;
 255     const edgeList& slaveEdgeList = slavePatch_.edges();
 256
 257     forAll (slaveEdgeList, sEi)
 258     {
 259         slaveMinEdgeLength =
 260             Foam::min
 261             (
 262                 slaveMinEdgeLength,
 263                 slaveEdgeList[sEi].mag(slavePatch_.localPoints())
 264             );
 265     }
 266
 267     // There is no point classifying data to a resolution smaller than the
 268     // largest of the two minimum edges found.
 269     // This will drive the size of the smallest bin needed to construct the
 270     // z and r histograms
 271     scalar halfMinSizeBin =
 272         Foam::max(masterMinEdgeLength, slaveMinEdgeLength)/2.0;
 273
 274     if (debug)
 275     {
 276         InfoIn
 277         (
 278             "tmp<pointField> MixingPlaneInterpolation::calcProfile()"
 279         )   << "halfMinSizeBin: " << halfMinSizeBin << endl;
 280     }
 281
 282     // Build master and slave histogram
 283     profileHistogram masterHistogram;
 284     profileHistogram slaveHistogram;
 285
 286     // stackingDir switch
 287     const direction stackingDir = stackAxisSwitch();
 288
 289     // Master side
 290     forAll (masterGlobalProfile, mI)
 291     {
 292         updateProfileHistogram
 293         (
 294             masterHistogram,
 295             masterGlobalProfile[mI],
 296             stackingDir,
 297             halfMinSizeBin
 298         );
 299     }
 300
 301     // Shadow side
 302     forAll (slaveGlobalProfile, sI)
 303     {
 304         updateProfileHistogram
 305         (
 306             slaveHistogram,
 307             slaveGlobalProfile[sI],
 308             stackingDir,
 309             halfMinSizeBin
 310         );
 311     }
 312
 313     if (debug > 1)
 314     {
 315         // Write histograms
 316         forAllIter (profileHistogram, masterHistogram, zHi)
 317         {
 318             Info<< "master histo (z, n): (" << zHi->first << " "
 319                 << static_cast<int>(zHi->second.size()) << ")" << endl;
 320         }
 321
 322         forAllIter (profileHistogram, slaveHistogram, zHi)
 323         {
 324             Info<< "slave histo (z, n): (" << zHi->first << " "
 325                 << static_cast<int>(zHi->second.size()) << ")" << endl;
 326         }
 327     }
 328
 329     pointField profileBeforeClip;
 330
 331     // Select which type of profile we need
 332     switch (discretisationType_)
 333     {
 334         // To Do: add uniform mixing plane
 335
 336         case MASTER_PATCH:
 337         {
 338             profileBeforeClip = computeProfileFromHistograms
 339             (
 340                 masterHistogram,
 341                 masterHistogram,
 342                 halfMinSizeBin
 343             );
 344         }
 345         break;
 346
 347         case SLAVE_PATCH:
 348         {
 349             profileBeforeClip = computeProfileFromHistograms
 350             (
 351                 slaveHistogram,
 352                 slaveHistogram,
 353                 halfMinSizeBin
 354             );
 355         }
 356         break;
 357
 358         case BOTH_PATCHES:
 359         {
 360             profileBeforeClip = computeProfileFromHistograms
 361             (
 362                 masterHistogram,
 363                 slaveHistogram,
 364                 halfMinSizeBin
 365             );
 366         }
 367         break;
 368
 369         case UNIFORM:
 370         case USER_DEFINED:
 371         default:
 372         {
 373             FatalErrorIn
 374             (
 375                 "tmp<pointField> MixingPlaneInterpolation<MasterPatch, "
 376                 "SlavePatch>::calcProfile() const"
 377             )   << "Bad type of mixing plane discretisation: "
 378                 << discretisationNames_[discretisationType_]
 379                 << nl << "Available types are: " << nl
 380                 << MixingPlaneInterpolationName::discretisationNames_
 381                 << abort(FatalError);
 382         }
 383     }
 384
 385     // Remove interpolationProfile_ points located outside of either
 386     // master/slave patch boundingBox,
 387     // with the exception of the first and last profile points
 388
 389     // Martin to work here. HJ, 27/Jan/2011
 390
 391     boundBox masterBB
 392     (
 393         masterGlobalProfile,
 394         false
 395     );
 396
 397     boundBox slaveBB
 398     (
 399         slaveGlobalProfile,
 400         false
 401     );
 402
 403     // Expand the bounding box in the sweepAxis-wise direction
 404     // Note: All points are collapsed to zero in sweepDir
 405     // It is sufficient to expand the box by 1 in this direction
 406
 407     masterBB.min().replace(sweepDir, -1);
 408     masterBB.max().replace(sweepDir, 1);
 409
 410     slaveBB.min().replace(sweepDir, -1);
 411     slaveBB.max().replace(sweepDir, 1);
 412
 413     boundBox globSpanBB
 414     (
 415         point(Foam::min(masterBB.min(), slaveBB.min()))
 416       - point(SMALL, SMALL, SMALL),
 417         point(Foam::max(masterBB.max(), slaveBB.max()))
 418       + point(SMALL, SMALL, SMALL)
 419     );
 420
 421     if (debug)
 422     {
 423         InfoIn
 424         (
 425             "tmp<pointField> MixingPlaneInterpolation<MasterPatch, "
 426             "SlavePatch>::calcProfile() const"
 427         )   << setprecision(12) << nl
 428             << "masterBB: " << masterBB << nl
 429             << "slaveBB: " <<  slaveBB << nl
 430             << "globSpanBB: " <<  globSpanBB << nl
 431             << "initial profile values: " << profileBeforeClip << endl;
 432     }
 433
 434     // Iterate through profile, removing points located
 435     // outside of either master/slave BB
 436
 437     tmp<pointField> tprofile(new pointField(profileBeforeClip.size()));
 438     pointField& profile = tprofile();
 439     label curIndex = 0;
 440
 441     forAll (profileBeforeClip, pI)
 442     {
 443         if (globSpanBB.contains(profileBeforeClip[pI]))
 444         {
 445             // We keep that profile point
 446             profile[curIndex] = profileBeforeClip[pI];
 447             curIndex++;  // Next slot
 448         }
 449         else
 450         {
 451             if (debug)
 452             {
 453                 InfoIn
 454                 (
 455                     "MixingPlaneInterpolation"
 456                     "<MasterPatch, SlavePatch>::"
 457                     "removeNonOverlappedProfilePoints"
 458                 )   << setprecision(12) << "   Removing point: "
 459                     << profileBeforeClip[pI] << endl;
 460             }
 461         }
 462     }
 463
 464     profile.setSize(curIndex);
 465
 466     if (profile.size() < 2)
 467     {
 468         FatalErrorIn
 469         (
 470             "tmp<pointField> MixingPlaneInterpolation<MasterPatch, "
 471             "SlavePatch>::calcProfile() const"
 472         )   << "Lost all points in profile: " << profile
 473             << abort(FatalError);
 474     }
 475
 476     if (debug)
 477     {
 478         InfoIn
 479         (
 480             "MixingPlaneInterpolation<MasterPatch, SlavePatch>::"
 481             "removeNonOverlappedProfilePoints"
 482         )   << "cleaned-up profile values: " << profile << endl;
 483     }
 484
 485     return tprofile;
 486 }
 487
 488
 489 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 490
 491 template<class MasterPatch, class SlavePatch>
 492 const Foam::pointField&
 493 MixingPlaneInterpolation<MasterPatch, SlavePatch>::interpolationProfile() const
 494 {
 495     if (interpolationProfile_.size() == 0)
 496     {
 497         // Not a user-defined profile: calculate as per algorithm
 498         interpolationProfile_ = calcProfile();
 499     }
 500
 501     return interpolationProfile_;
 502 }
 503
 504
 505 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 506
 507 } // End namespace Foam
 508
 509 // ************************************************************************* //