src/dynamicMesh/dynamicTopoFvMesh/lengthScaleEstimator/lengthScaleEstimator.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 Class
  25     lengthScaleEstimator
  26
  27 Description
  28     Implementation of lengthScaleEstimator members
  29
  30 Author
  31     Sandeep Menon
  32     University of Massachusetts Amherst
  33     All rights reserved
  34
  35 \*----------------------------------------------------------------------------*/
  36
  37 #include "fvc.H"
  38 #include "volFields.H"
  39 #include "lengthScaleEstimator.H"
  40 #include "processorPolyPatch.H"
  41
  42 namespace Foam
  43 {
  44
  45 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  46
  47 defineTypeNameAndDebug(lengthScaleEstimator, 0);
  48
  49 lengthScaleEstimator::ScaleFnPair lengthScaleEstimator::methods_[] =
  50 {
  51     ScaleFnPair("constant", &lengthScaleEstimator::constantScale),
  52     ScaleFnPair("direct", &lengthScaleEstimator::directScale),
  53     ScaleFnPair("inverse", &lengthScaleEstimator::inverseScale)
  54 };
  55
  56 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  57
  58 // Construct from polyMesh and dictionary
  59 lengthScaleEstimator::lengthScaleEstimator
  60 (
  61     const polyMesh& mesh
  62 )
  63 :
  64     mesh_(mesh),
  65     ratioMin_(0.0),
  66     ratioMax_(0.0),
  67     growthFactor_(1.0),
  68     minLengthScale_(VSMALL),
  69     maxLengthScale_(GREAT),
  70     curvatureDeviation_(0.0),
  71     spatialRes_(-1),
  72     proximityBins_(0),
  73     sliceThreshold_(VSMALL),
  74     sliceHoldOff_(0),
  75     sliceBoxes_(0),
  76     field_("none"),
  77     scale_(NULL),
  78     gradient_(false),
  79     fieldLength_(0.0),
  80     lowerRefineLevel_(0.001),
  81     upperRefineLevel_(0.999),
  82     meanScale_(-1.0),
  83     maxRefineLevel_(labelMax)
  84 {}
  85
  86 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  87
  88 lengthScaleEstimator::~lengthScaleEstimator()
  89 {}
  90
  91 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  92
  93 // Check for legitimacy of patches
  94 void lengthScaleEstimator::checkPatches(const wordList& patchList) const
  95 {
  96     const polyBoundaryMesh& boundary = mesh_.boundaryMesh();
  97
  98     forAll(patchList, wordI)
  99     {
 100         if (boundary.findPatchID(patchList[wordI]) < 0)
 101         {
 102             FatalErrorIn
 103             (
 104                 "void lengthScaleEstimator::checkPatches"
 105                 "(const wordList& patchList) const"
 106             )
 107                 << " Could not find patch: "
 108                 << patchList[wordI] << nl
 109                 << abort(FatalError);
 110         }
 111     }
 112 }
 113
 114
 115 // Prepare for proximity-based refinement, if necessary
 116 void lengthScaleEstimator::prepareProximityPatches()
 117 {
 118     if (!proximityPatches_.size())
 119     {
 120         return;
 121     }
 122
 123     if (debug)
 124     {
 125         Info << "Preparing patches for proximity-based refinement...";
 126     }
 127
 128     // Clear out existing lists.
 129     proximityBins_.clear();
 130     proximityBins_.setSize(997, labelList(0));
 131
 132     const polyBoundaryMesh& boundary = mesh_.boundaryMesh();
 133
 134     bool setSpatialRes = false;
 135
 136     // Loop through all proximity patches and spatially hash patch faces.
 137     forAll(boundary, patchI)
 138     {
 139         if
 140         (
 141             (proximityPatches_.found(boundary[patchI].name())) ||
 142             (boundary[patchI].type() == "symmetryPlane")
 143         )
 144         {
 145             const polyPatch& proxPatch = boundary[patchI];
 146
 147             // Construct a bounding-box of face centres.
 148             // Do not synchronize in parallel, since the patch
 149             // may not be present on all sub-domains.
 150             boundBox box(proxPatch.faceCentres(), false);
 151
 152             const point& bMin = box.min();
 153             const point& bMax = box.max();
 154
 155             // Further hashing requires this information.
 156             proxBoundBox_.min() = bMin;
 157             proxBoundBox_.max() = bMax;
 158
 159             // Build a list of face indices
 160             labelList faceIndices
 161             (
 162                 identity(proxPatch.size()) + proxPatch.start()
 163             );
 164
 165             // For spatial resolution, pick a face on this patch.
 166             if (!setSpatialRes)
 167             {
 168                 scalar eLength =
 169                 (
 170                     Foam::sqrt(mag(proxPatch.faceAreas()[0]))
 171                 );
 172
 173                 spatialRes_ =
 174                 (
 175                     label(::floor(mag(bMax - bMin) / (3.0 * eLength)))
 176                 );
 177
 178                 setSpatialRes = true;
 179             }
 180
 181             spatialHash
 182             (
 183                 proxPatch.faceCentres(),
 184                 faceIndices,
 185                 box,
 186                 spatialRes_,
 187                 proximityBins_
 188             );
 189         }
 190     }
 191
 192     if (debug)
 193     {
 194         Info << "Done." << endl;
 195     }
 196 }
 197
 198
 199 // Perform spatial hashing on a set of points
 200 void lengthScaleEstimator::spatialHash
 201 (
 202     const pointField& pointLocations,
 203     const labelList& pointIndices,
 204     const boundBox& box,
 205     const label resolution,
 206     labelListList& bins
 207 )
 208 {
 209     label binSize = bins.size(), nD = resolution;
 210
 211     const point& bMin = box.min();
 212     const point& bMax = box.max();
 213
 214     // Extend bounding-box dimensions a bit to avoid edge-effects.
 215     scalar ext = 0.02*(mag(bMax - bMin));
 216
 217     // Define an inverse grid-cell size.
 218     scalar xL = nD/(bMax.x() - bMin.x() + ext);
 219     scalar yL = nD/(bMax.y() - bMin.y() + ext);
 220     scalar zL = nD/(bMax.z() - bMin.z() + ext);
 221
 222     // Loop through all points and bin them.
 223     forAll(pointLocations, pointI)
 224     {
 225         // Translate to boundBox minimum.
 226         point p = pointLocations[pointI] - bMin;
 227
 228         // Hash the position.
 229         label i = label(mag(::floor(p.x()*xL)));
 230         label j = label(mag(::floor(p.y()*yL)));
 231         label k = label(mag(::floor(p.z()*zL)));
 232
 233         label pos = mag(((k*nD*nD)+(j*nD)+i) % binSize);
 234
 235         // Store the index.
 236         label newSize = bins[pos].size() + 1;
 237         bins[pos].setSize(newSize, pointIndices[pointI]);
 238     }
 239 }
 240
 241
 242 // Send length-scale info across processors
 243 void lengthScaleEstimator::writeLengthScaleInfo
 244 (
 245     const labelList& cellLevels,
 246     const scalarList& lengthScale
 247 )
 248 {
 249     const polyBoundaryMesh& boundary = mesh_.boundaryMesh();
 250
 251     // Clear existing buffers
 252     sendLblBuffer_.clear();
 253     recvLblBuffer_.clear();
 254     sendLvlBuffer_.clear();
 255     recvLvlBuffer_.clear();
 256     sendSclBuffer_.clear();
 257     recvSclBuffer_.clear();
 258
 259     // Number of sent faces across processors
 260     labelList nSendFaces(boundary.size(), 0);
 261     labelList nRecvFaces(boundary.size(), 0);
 262
 263     // Corresponding face labels
 264     sendLblBuffer_.setSize(boundary.size());
 265     sendLvlBuffer_.setSize(boundary.size());
 266     recvLblBuffer_.setSize(boundary.size());
 267     recvLvlBuffer_.setSize(boundary.size());
 268
 269     // Length-scales corresponding to face-labels
 270     sendSclBuffer_.setSize(boundary.size());
 271     recvSclBuffer_.setSize(boundary.size());
 272
 273     // Fill send buffers with cell-level and length-scale info.
 274     forAll(boundary, pI)
 275     {
 276         if (isA<processorPolyPatch>(boundary[pI]))
 277         {
 278             const labelList& fCells = boundary[pI].faceCells();
 279
 280             // Set the initial buffer size
 281             sendLblBuffer_[pI].setSize(boundary[pI].size(), 0);
 282             sendLvlBuffer_[pI].setSize(boundary[pI].size(), 0);
 283             sendSclBuffer_[pI].setSize(boundary[pI].size(), 0.0);
 284
 285             forAll(fCells, faceI)
 286             {
 287                 label cI = fCells[faceI];
 288
 289                 // Does the adjacent cell have a non-zero level?
 290                 if (cellLevels[cI] > 0)
 291                 {
 292                     // Fill the send buffer.
 293                     sendLblBuffer_[pI][nSendFaces[pI]] = faceI;
 294                     sendLvlBuffer_[pI][nSendFaces[pI]] = cellLevels[cI];
 295                     sendSclBuffer_[pI][nSendFaces[pI]] = lengthScale[cI];
 296
 297                     nSendFaces[pI]++;
 298                 }
 299             }
 300
 301             // Resize to actual value
 302             sendLblBuffer_[pI].setSize(nSendFaces[pI]);
 303             sendLvlBuffer_[pI].setSize(nSendFaces[pI]);
 304             sendSclBuffer_[pI].setSize(nSendFaces[pI]);
 305
 306             const processorPolyPatch& pp =
 307             (
 308                 refCast<const processorPolyPatch>(boundary[pI])
 309             );
 310
 311             label neiProcNo = pp.neighbProcNo();
 312
 313             // First perform a blocking send/receive of the number of faces.
 314             OPstream::write
 315             (
 316                 Pstream::nonBlocking,
 317                 neiProcNo,
 318                 reinterpret_cast<const char*>(&(nSendFaces[pI])),
 319                 sizeof(label)
 320             );
 321
 322             IPstream::read
 323             (
 324                 Pstream::nonBlocking,
 325                 neiProcNo,
 326                 reinterpret_cast<char*>(&(nRecvFaces[pI])),
 327                 sizeof(label)
 328             );
 329         }
 330     }
 331
 332     // Wait for all transfers to complete.
 333     OPstream::waitRequests();
 334     IPstream::waitRequests();
 335
 336     // Send info to neighbouring processors.
 337     forAll(boundary, patchI)
 338     {
 339         if (isA<processorPolyPatch>(boundary[patchI]))
 340         {
 341             const processorPolyPatch& pp =
 342             (
 343                 refCast<const processorPolyPatch>(boundary[patchI])
 344             );
 345
 346             label neiProcNo = pp.neighbProcNo();
 347
 348             if (debug > 4)
 349             {
 350                 Pout << " Processor patch " << patchI << ' ' << pp.name()
 351                      << " communicating with " << neiProcNo
 352                      << "  Sending: " << nSendFaces[patchI]
 353                      << "  Recving: " << nRecvFaces[patchI]
 354                      << endl;
 355             }
 356
 357             // Next, perform a non-blocking send/receive of buffers.
 358             // But only if the buffer size is non-zero.
 359             if (nSendFaces[patchI] != 0)
 360             {
 361                 OPstream::write
 362                 (
 363                     Pstream::nonBlocking,
 364                     neiProcNo,
 365                     reinterpret_cast<const char*>(&(sendLblBuffer_[patchI][0])),
 366                     sendLblBuffer_[patchI].size()*sizeof(label)
 367                 );
 368
 369                 OPstream::write
 370                 (
 371                     Pstream::nonBlocking,
 372                     neiProcNo,
 373                     reinterpret_cast<const char*>(&(sendLvlBuffer_[patchI][0])),
 374                     sendLblBuffer_[patchI].size()*sizeof(label)
 375                 );
 376
 377                 OPstream::write
 378                 (
 379                     Pstream::nonBlocking,
 380                     neiProcNo,
 381                     reinterpret_cast<const char*>(&(sendSclBuffer_[patchI][0])),
 382                     sendSclBuffer_[patchI].size()*sizeof(scalar)
 383                 );
 384             }
 385
 386             if (nRecvFaces[patchI] != 0)
 387             {
 388                 // Size the receive buffers
 389                 recvLblBuffer_[patchI].setSize(nRecvFaces[patchI], 0);
 390                 recvLvlBuffer_[patchI].setSize(nRecvFaces[patchI], 0);
 391                 recvSclBuffer_[patchI].setSize(nRecvFaces[patchI], 0.0);
 392
 393                 IPstream::read
 394                 (
 395                     Pstream::nonBlocking,
 396                     neiProcNo,
 397                     reinterpret_cast<char*>(&(recvLblBuffer_[patchI][0])),
 398                     nRecvFaces[patchI]*sizeof(label)
 399                 );
 400
 401                 IPstream::read
 402                 (
 403                     Pstream::nonBlocking,
 404                     neiProcNo,
 405                     reinterpret_cast<char*>(&(recvLvlBuffer_[patchI][0])),
 406                     nRecvFaces[patchI]*sizeof(label)
 407                 );
 408
 409                 IPstream::read
 410                 (
 411                     Pstream::nonBlocking,
 412                     neiProcNo,
 413                     reinterpret_cast<char*>(&(recvSclBuffer_[patchI][0])),
 414                     nRecvFaces[patchI]*sizeof(scalar)
 415                 );
 416             }
 417         }
 418     }
 419 }
 420
 421
 422 // Receive length-scale info across processors
 423 void lengthScaleEstimator::readLengthScaleInfo
 424 (
 425     const label level,
 426     label& visitedCells,
 427     labelList& cellLevels,
 428     UList<scalar>& lengthScale,
 429     labelHashSet& levelCells
 430 ) const
 431 {
 432     const polyBoundaryMesh& boundary = mesh_.boundaryMesh();
 433
 434     const labelList& own = mesh_.faceOwner();
 435     const labelList& nei = mesh_.faceNeighbour();
 436
 437     // Wait for all transfers to complete.
 438     OPstream::waitRequests();
 439     IPstream::waitRequests();
 440
 441     // Now re-visit cells and update length-scales.
 442     forAll(boundary, patchI)
 443     {
 444         if (recvLblBuffer_[patchI].size())
 445         {
 446             const labelList& fCells = boundary[patchI].faceCells();
 447
 448             forAll(recvLblBuffer_[patchI], i)
 449             {
 450                 label nLabel = recvLblBuffer_[patchI][i];
 451                 label pLevel = recvLvlBuffer_[patchI][i];
 452
 453                 label cI = fCells[nLabel];
 454
 455                 label& ngbLevel = cellLevels[cI];
 456
 457                 // For an unvisited cell, update the level
 458                 bool unvisited = false;
 459
 460                 if (ngbLevel == 0)
 461                 {
 462                     ngbLevel = pLevel + 1;
 463                     unvisited = true;
 464                 }
 465
 466                 // Visit all neighbours and re-calculate length-scale
 467                 const cell& cellCheck = mesh_.cells()[cI];
 468
 469                 scalar sumLength = 0.0;
 470                 label nTouchedNgb = 0;
 471
 472                 forAll(cellCheck, faceI)
 473                 {
 474                     label sLevel = -1, sLCell = -1;
 475                     label pF = boundary.whichPatch(cellCheck[faceI]);
 476
 477                     if (pF == -1)
 478                     {
 479                         // Internal face. Determine neighbour.
 480                         if (own[cellCheck[faceI]] == cI)
 481                         {
 482                             sLCell = nei[cellCheck[faceI]];
 483                         }
 484                         else
 485                         {
 486                             sLCell = own[cellCheck[faceI]];
 487                         }
 488
 489                         sLevel = cellLevels[sLCell];
 490
 491                         if ((sLevel < ngbLevel) && (sLevel > 0))
 492                         {
 493                             sumLength += lengthScale[sLCell];
 494
 495                             nTouchedNgb++;
 496                         }
 497                     }
 498                     else
 499                     if (isA<processorPolyPatch>(boundary[pF]))
 500                     {
 501                         // Determine the local index.
 502                         label local =
 503                         (
 504                             boundary[pF].whichFace(cellCheck[faceI])
 505                         );
 506
 507                         // Is this label present in the list?
 508                         label j = -1;
 509
 510                         if ((j = findIndex(recvLblBuffer_[pF], local)) > -1)
 511                         {
 512                             if
 513                             (
 514                                 (recvLvlBuffer_[pF][j] < ngbLevel) &&
 515                                 (recvLvlBuffer_[pF][j] > 0)
 516                             )
 517                             {
 518                                 sumLength += recvSclBuffer_[pF][j];
 519
 520                                 nTouchedNgb++;
 521                             }
 522                         }
 523                     }
 524                     else
 525                     if (!isFreePatch(pF))
 526                     {
 527                         sumLength +=
 528                         (
 529                             fixedLengthScale(cellCheck[faceI], pF, true)
 530                         );
 531
 532                         nTouchedNgb++;
 533                     }
 534                 }
 535
 536                 sumLength /= nTouchedNgb;
 537
 538                 // Scale the length and assign to this cell
 539                 if (level < maxRefineLevel_)
 540                 {
 541                     sumLength *= growthFactor_;
 542                 }
 543                 else
 544                 if (meanScale_ > 0.0)
 545                 {
 546                     // If a mean scale has been specified,
 547                     // override the value
 548                     sumLength = meanScale_;
 549                 }
 550
 551                 lengthScale[cI] = sumLength;
 552
 553                 if (unvisited)
 554                 {
 555                     levelCells.insert(cI);
 556
 557                     visitedCells++;
 558                 }
 559             }
 560         }
 561     }
 562 }
 563
 564
 565 // Use a constant length scale for field-based refinement
 566 scalar lengthScaleEstimator::constantScale(const scalar fieldValue) const
 567 {
 568     return fieldLength_;
 569 }
 570
 571
 572 // Use a direct-proportion length scale for field-based refinement
 573 scalar lengthScaleEstimator::directScale(const scalar fieldValue) const
 574 {
 575     if (fieldValue > upperRefineLevel_)
 576     {
 577         return meanScale_;
 578     }
 579
 580     scalar value = fieldValue;
 581
 582     // Bracket the value within limits
 583     value = min(value, upperRefineLevel_);
 584     value = max(value, lowerRefineLevel_);
 585
 586     scalar ratio
 587     (
 588         (value - lowerRefineLevel_) / (upperRefineLevel_ - lowerRefineLevel_)
 589     );
 590
 591     scalar scale
 592     (
 593         minLengthScale_ + ((maxLengthScale_ - minLengthScale_) * ratio)
 594     );
 595
 596     return scale;
 597 }
 598
 599
 600 // Use an inverse-proportion length scale for field-based refinement
 601 scalar lengthScaleEstimator::inverseScale(const scalar fieldValue) const
 602 {
 603     if (fieldValue < lowerRefineLevel_)
 604     {
 605         return meanScale_;
 606     }
 607
 608     scalar value = fieldValue;
 609
 610     // Bracket the value within limits
 611     value = min(value, upperRefineLevel_);
 612     value = max(value, lowerRefineLevel_);
 613
 614     scalar ratio
 615     (
 616         (value - lowerRefineLevel_) / (upperRefineLevel_ - lowerRefineLevel_)
 617     );
 618
 619     scalar scale
 620     (
 621         minLengthScale_ + ((maxLengthScale_ - minLengthScale_) * (1.0 - ratio))
 622     );
 623
 624     return scale;
 625 }
 626
 627
 628 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 629
 630 // Read edge refinement options from the dictionary
 631 void lengthScaleEstimator::readRefinementOptions
 632 (
 633     const dictionary& refineDict,
 634     bool reRead,
 635     bool mandatory
 636 )
 637 {
 638     // Save old values
 639     scalar oldRatioMax = ratioMax_;
 640     scalar oldRatioMin = ratioMin_;
 641     scalar oldGrowthFactor = growthFactor_;
 642     scalar oldMinScale = minLengthScale_;
 643     scalar oldMaxScale = maxLengthScale_;
 644
 645     ratioMax_ = readScalar(refineDict.lookup("bisectionRatio"));
 646     ratioMin_ = readScalar(refineDict.lookup("collapseRatio"));
 647     growthFactor_ = readScalar(refineDict.lookup("growthFactor"));
 648
 649     // Sanity check: Are ratios and growth-factor correctly specified?
 650     if
 651     (
 652         ratioMin_ > ratioMax_ ||
 653         ratioMin_ < VSMALL ||
 654         ratioMax_ < VSMALL ||
 655         growthFactor_ < 1.0
 656     )
 657     {
 658         FatalErrorIn
 659         (
 660             "void lengthScaleEstimator::readRefinementOptions"
 661             "(const dictionary&, bool, bool)"
 662         )
 663             << " Options are incorrectly specified." << nl
 664             << " ratioMin: " << ratioMin_ << nl
 665             << " ratioMax: " << ratioMax_ << nl
 666             << " growthFactor: " << growthFactor_ << nl
 667             << abort(FatalError);
 668     }
 669
 670     if (refineDict.found("minLengthScale") || mandatory)
 671     {
 672         minLengthScale_ = readScalar(refineDict.lookup("minLengthScale"));
 673     }
 674
 675     if (refineDict.found("maxLengthScale") || mandatory)
 676     {
 677         maxLengthScale_ = readScalar(refineDict.lookup("maxLengthScale"));
 678     }
 679
 680     // Sanity check: Are min/max length scales correctly specified?
 681     if (minLengthScale_ > maxLengthScale_)
 682     {
 683         FatalErrorIn
 684         (
 685             "void lengthScaleEstimator::readRefinementOptions"
 686             "(const dictionary&, bool, bool)"
 687         )
 688             << " Length-scales are incorrectly specified." << nl
 689             << " minLengthScale: " << minLengthScale_ << nl
 690             << " maxLengthScale: " << maxLengthScale_ << nl
 691             << abort(FatalError);
 692     }
 693
 694     if (reRead)
 695     {
 696         // Check if values have changed, and report it.
 697         if (mag(oldRatioMax - ratioMax_) > SMALL)
 698         {
 699             Info << "\tOld ratioMax: " << oldRatioMax << nl
 700                  << "\tNew ratioMax: " << ratioMax_ << endl;
 701         }
 702
 703         if (mag(oldRatioMin - ratioMin_) > SMALL)
 704         {
 705             Info << "\tOld ratioMin: " << oldRatioMin << nl
 706                  << "\tNew ratioMin: " << ratioMin_ << endl;
 707         }
 708
 709         if (mag(oldGrowthFactor - growthFactor_) > SMALL)
 710         {
 711             Info << "\tOld growthFactor: " << oldGrowthFactor << nl
 712                  << "\tNew growthFactor: " << growthFactor_ << endl;
 713         }
 714
 715         if (mag(oldMinScale - minLengthScale_) > SMALL)
 716         {
 717             Info << "\tOld minLengthScale: " << oldMinScale << nl
 718                  << "\tNew minLengthScale: " << minLengthScale_ << endl;
 719         }
 720
 721         if (mag(oldMaxScale - maxLengthScale_) > SMALL)
 722         {
 723             Info << "\tOld maxLengthScale: " << oldMaxScale << nl
 724                  << "\tNew maxLengthScale: " << maxLengthScale_ << endl;
 725         }
 726     }
 727
 728     if (refineDict.found("fixedLengthScalePatches") || mandatory)
 729     {
 730         fixedPatches_ = refineDict.subDict("fixedLengthScalePatches");
 731
 732         // Ensure that patches are legitimate.
 733         checkPatches(fixedPatches_.toc());
 734     }
 735
 736     if (refineDict.found("freeLengthScalePatches") || mandatory)
 737     {
 738         freePatches_ = refineDict.subDict("freeLengthScalePatches");
 739
 740         // Ensure that patches are legitimate.
 741         checkPatches(freePatches_.toc());
 742
 743         // Check if fixed and free patches are conflicting
 744         if (fixedPatches_.size() && freePatches_.size())
 745         {
 746             wordList fixedPatchList = fixedPatches_.toc();
 747             wordList freePatchList = freePatches_.toc();
 748
 749             forAll(fixedPatchList, wordI)
 750             {
 751                 if (findIndex(freePatchList, fixedPatchList[wordI]) > -1)
 752                 {
 753                     FatalErrorIn
 754                     (
 755                         "void lengthScaleEstimator::readRefinementOptions"
 756                         "(const dictionary&, bool, bool)"
 757                     )
 758                         << " Conflicting fixed/free patches." << nl
 759                         << " Patch: " << fixedPatchList[wordI] << nl
 760                         << abort(FatalError);
 761                 }
 762             }
 763         }
 764     }
 765
 766     if (refineDict.found("curvaturePatches") || mandatory)
 767     {
 768         curvaturePatches_ = refineDict.subDict("curvaturePatches");
 769
 770         // Ensure that patches are legitimate.
 771         checkPatches(curvaturePatches_.toc());
 772
 773         curvatureDeviation_ =
 774         (
 775             readScalar(refineDict.lookup("curvatureDeviation"))
 776         );
 777
 778         if (curvatureDeviation_ > 1.0 || curvatureDeviation_ < 0.0)
 779         {
 780             FatalErrorIn
 781             (
 782                 "void lengthScaleEstimator::readRefinementOptions"
 783                 "(const dictionary&, bool, bool)"
 784             )
 785                 << " Curvature deviation out of range [0..1]"
 786                 << abort(FatalError);
 787         }
 788     }
 789
 790     if (refineDict.found("proximityPatches") || mandatory)
 791     {
 792         proximityPatches_ =
 793         (
 794             refineDict.subDict("proximityPatches")
 795         );
 796
 797         // Ensure that patches are legitimate.
 798         checkPatches(proximityPatches_.toc());
 799
 800         // Check if a threshold for slicing has been specified.
 801         if (refineDict.found("sliceThreshold") || mandatory)
 802         {
 803             sliceThreshold_ =
 804             (
 805                 readScalar(refineDict.lookup("sliceThreshold"))
 806             );
 807
 808             // Cap the threshold value
 809             // sliceThreshold_ = Foam::max(sliceThreshold_, minLengthScale_);
 810         }
 811     }
 812
 813     // Check if edge-refinement is to be avoided on any patches
 814     if (refineDict.found("noModificationPatches") || mandatory)
 815     {
 816         wordList noModPatches =
 817         (
 818             refineDict.subDict("noModificationPatches").toc()
 819         );
 820
 821         // Ensure that patches are legitimate.
 822         checkPatches(noModPatches);
 823
 824         noModPatchIDs_.setSize(noModPatches.size());
 825
 826         label indexI = 0;
 827
 828         forAll(noModPatches, wordI)
 829         {
 830             const word& patchName = noModPatches[wordI];
 831
 832             noModPatchIDs_[indexI++] =
 833             (
 834                 mesh_.boundaryMesh().findPatchID(patchName)
 835             );
 836         }
 837     }
 838
 839     // Check if refinement is to be performed based on a field-value
 840     if (refineDict.found("fieldRefinement") || mandatory)
 841     {
 842         field_ = word(refineDict.lookup("fieldRefinement"));
 843
 844         word scaleMethod(refineDict.lookup("fieldScaleMethod"));
 845
 846         const label nMethods(sizeof(methods_) / sizeof(methods_[0]));
 847
 848         bool invalidMethod = true;
 849
 850         for (label i = 0; i < nMethods; i++)
 851         {
 852             if (scaleMethod == methods_[i].first())
 853             {
 854                 scale_ = methods_[i].second();
 855                 invalidMethod = false;
 856                 break;
 857             }
 858         }
 859
 860         if (invalidMethod)
 861         {
 862             Info << nl << " Available methods: " << endl;
 863
 864             for (label i = 0; i < nMethods; i++)
 865             {
 866                 Info << "    " << methods_[i].first() << nl;
 867             }
 868
 869             FatalErrorIn
 870             (
 871                 "void lengthScaleEstimator::readRefinementOptions"
 872                 "(const dictionary&, bool, bool)"
 873             )
 874                 << " Invalid method: " << scaleMethod << " specified." << nl
 875                 << abort(FatalError);
 876         }
 877
 878         // Define whether the gradient must be evaluated
 879         gradient_ = readBool(refineDict.lookup("fieldGradient"));
 880
 881         // Lookup a specified length-scale
 882         fieldLength_ = readScalar(refineDict.lookup("fieldLengthScale"));
 883
 884         lowerRefineLevel_ = readScalar(refineDict.lookup("lowerRefineLevel"));
 885         upperRefineLevel_ = readScalar(refineDict.lookup("upperRefineLevel"));
 886
 887         // Sanity check: Are refinement levels correctly specified?
 888         if (lowerRefineLevel_ > upperRefineLevel_)
 889         {
 890             FatalErrorIn
 891             (
 892                 "void lengthScaleEstimator::readRefinementOptions"
 893                 "(const dictionary&, bool, bool)"
 894             )
 895                 << " Refinement levels are incorrectly specified." << nl
 896                 << " lowerRefineLevel: " << lowerRefineLevel_ << nl
 897                 << " upperRefineLevel: " << upperRefineLevel_ << nl
 898                 << abort(FatalError);
 899         }
 900     }
 901
 902     // Check if a max refinement level has been specified
 903     if (refineDict.found("maxRefineLevel") || mandatory)
 904     {
 905         maxRefineLevel_ = readLabel(refineDict.lookup("maxRefineLevel"));
 906     }
 907
 908     // Update switch for mean-scale
 909     if (refineDict.found("meanScale") || mandatory)
 910     {
 911         meanScale_ = readScalar(refineDict.lookup("meanScale"));
 912     }
 913 }
 914
 915
 916 // Set explicitly coupled patch information
 917 void lengthScaleEstimator::setCoupledPatches
 918 (
 919     const dictionary& coupledPatches
 920 )
 921 {
 922     label indexI = 0;
 923
 924     // Determine master and slave patches
 925     wordList masterPatches(coupledPatches.size());
 926     wordList slavePatches(coupledPatches.size());
 927
 928     forAllConstIter(dictionary, coupledPatches, dIter)
 929     {
 930         const dictionary& dictI = dIter().dict();
 931
 932         masterPatches[indexI] = word(dictI.lookup("master"));
 933         slavePatches[indexI] = word(dictI.lookup("slave"));
 934
 935         indexI++;
 936     }
 937
 938     // Ensure that patches are legitimate.
 939     checkPatches(masterPatches);
 940
 941     // Check whether coupled patches are fixedPatches as well.
 942     forAll(masterPatches, wordI)
 943     {
 944         word pName(masterPatches[wordI]);
 945
 946         if (fixedPatches_.found(pName))
 947         {
 948             // Add the slave patch to the list as well.
 949             // If it already exists, over-ride the value.
 950             fixedPatches_.add
 951             (
 952                 slavePatches[wordI],
 953                 fixedPatches_[pName][0].scalarToken(),
 954                 true
 955             );
 956         }
 957     }
 958 }
 959
 960
 961 // Calculate the length scale field
 962 void lengthScaleEstimator::calculateLengthScale
 963 (
 964     UList<scalar>& lengthScale
 965 )
 966 {
 967     label level = 1, visitedCells = 0;
 968     labelList cellLevels(mesh_.nCells(), 0);
 969
 970     // Check for allocation
 971     if (lengthScale.size() != mesh_.nCells())
 972     {
 973         FatalErrorIn
 974         (
 975             "void lengthScaleEstimator::calculateLengthScale"
 976             "(UList<scalar>& lengthScale)"
 977         )
 978             << " Field is incorrectly sized." << nl
 979             << " Field size: " << lengthScale.size()
 980             << " nCells: " << mesh_.nCells()
 981             << abort(FatalError);
 982     }
 983
 984     // HashSet to keep track of cells in each level
 985     labelHashSet levelCells;
 986
 987     // Prepare for proximity-based refinement, if necessary
 988     prepareProximityPatches();
 989
 990     // Obtain the cellCells addressing list
 991     const labelList& own = mesh_.faceOwner();
 992     const labelListList& cc = mesh_.cellCells();
 993     const polyBoundaryMesh& boundary = mesh_.boundaryMesh();
 994
 995     forAll(boundary,patchI)
 996     {
 997         // Skip floating length-scale patches
 998         if (isFreePatch(patchI))
 999         {
1000             continue;
1001         }
1002
1003         const polyPatch& bdyPatch = boundary[patchI];
1004
1005         label pStart = bdyPatch.start();
1006
1007         forAll(bdyPatch,faceI)
1008         {
1009             label ownCell = own[pStart+faceI];
1010
1011             if (cellLevels[ownCell] != 0)
1012             {
1013                 continue;
1014             }
1015
1016             cellLevels[ownCell] = level;
1017
1018             lengthScale[ownCell] =
1019             (
1020                 fixedLengthScale(pStart+faceI, patchI, true) * growthFactor_
1021             );
1022
1023             levelCells.insert(ownCell);
1024
1025             visitedCells++;
1026         }
1027     }
1028
1029     // If a field has been specified, use that.
1030     if (field_ != "none")
1031     {
1032         const volScalarField* vFldPtr = NULL;
1033
1034         // Temporary pointer for gradient
1035         volScalarField* magGradPtr = NULL;
1036
1037         if (gradient_)
1038         {
1039             // Lookup various field types, and evaluate the gradient
1040             bool invalidObject = true;
1041
1042             // Evaluate using gradient scheme
1043             word gradName("grad(" + field_ + ')');
1044
1045             // Register field under a name that's unique
1046             word registerName("lengthScaleGradient(" + field_ + ')');
1047
1048             if (mesh_.objectRegistry::foundObject<volScalarField>(field_))
1049             {
1050                 const volScalarField& field =
1051                 (
1052                     mesh_.objectRegistry::lookupObject<volScalarField>(field_)
1053                 );
1054
1055                 magGradPtr =
1056                 (
1057                     new volScalarField
1058                     (
1059                         IOobject
1060                         (
1061                             registerName,
1062                             mesh_.time().timeName(),
1063                             mesh_,
1064                             IOobject::NO_READ,
1065                             IOobject::NO_WRITE,
1066                             false
1067                         ),
1068                         mag(fvc::grad(field, gradName))
1069                     )
1070                 );
1071
1072                 invalidObject = false;
1073             }
1074
1075             if (mesh_.objectRegistry::foundObject<volVectorField>(field_))
1076             {
1077                 const volVectorField& field =
1078                 (
1079                     mesh_.objectRegistry::lookupObject<volVectorField>(field_)
1080                 );
1081
1082                 magGradPtr =
1083                 (
1084                     new volScalarField
1085                     (
1086                         IOobject
1087                         (
1088                             registerName,
1089                             mesh_.time().timeName(),
1090                             mesh_,
1091                             IOobject::NO_READ,
1092                             IOobject::NO_WRITE,
1093                             false
1094                         ),
1095                         mag(fvc::grad(field, gradName))
1096                     )
1097                 );
1098
1099                 invalidObject = false;
1100             }
1101
1102             if (invalidObject)
1103             {
1104                 FatalErrorIn
1105                 (
1106                     "void lengthScaleEstimator::calculateLengthScale"
1107                     "(UList<scalar>& lengthScale)"
1108                 )
1109                     << " Invalid field: " << field_
1110                     << " specified for gradient." << nl
1111                     << abort(FatalError);
1112             }
1113
1114             vFldPtr = magGradPtr;
1115         }
1116         else
1117         {
1118             // Lookup a scalar field by default
1119             const volScalarField& field =
1120             (
1121                 mesh_.objectRegistry::lookupObject<volScalarField>(field_)
1122             );
1123
1124             vFldPtr = &field;
1125         }
1126
1127         const volScalarField& vFld = *vFldPtr;
1128
1129         const labelList& own = mesh_.faceOwner();
1130         const labelList& nei = mesh_.faceNeighbour();
1131
1132         forAll(nei, faceI)
1133         {
1134             label ownCell = own[faceI], neiCell = nei[faceI];
1135
1136             scalar fAvg = 0.5 * (vFld[ownCell] + vFld[neiCell]);
1137
1138             if ((fAvg > lowerRefineLevel_) && (fAvg < upperRefineLevel_))
1139             {
1140                 // Set the scale for cells on either side
1141                 if (!cellLevels[ownCell])
1142                 {
1143                     cellLevels[ownCell] = level;
1144                     lengthScale[ownCell] = (this->*scale_)(fAvg);
1145
1146                     levelCells.insert(ownCell);
1147
1148                     visitedCells++;
1149                 }
1150
1151                 if (!cellLevels[neiCell])
1152                 {
1153                     cellLevels[neiCell] = level;
1154                     lengthScale[neiCell] = (this->*scale_)(fAvg);
1155
1156                     levelCells.insert(neiCell);
1157
1158                     visitedCells++;
1159                 }
1160             }
1161         }
1162
1163         // Clean up gradient, if necessary
1164         if (gradient_)
1165         {
1166             deleteDemandDrivenData(magGradPtr);
1167         }
1168     }
1169
1170     bool doneWithSweeps = false;
1171
1172     // Perform multiple sweeps through the mesh...
1173     while (!doneWithSweeps)
1174     {
1175         if (Pstream::parRun())
1176         {
1177             writeLengthScaleInfo(cellLevels, lengthScale);
1178         }
1179
1180         // Loop through cells of the current level
1181         labelList currLvlCells = levelCells.toc();
1182         levelCells.clear();
1183
1184         // Loop through cells, and increment neighbour
1185         // cells of the current level
1186         forAll(currLvlCells,cellI)
1187         {
1188             // Obtain the cells neighbouring this one
1189             const labelList& cList = cc[currLvlCells[cellI]];
1190
1191             forAll(cList, indexI)
1192             {
1193                 label& ngbLevel = cellLevels[cList[indexI]];
1194
1195                 if (ngbLevel == 0)
1196                 {
1197                     ngbLevel = level + 1;
1198
1199                     // Compute the mean of the existing
1200                     // neighbour length-scales
1201                     const labelList& ncList = cc[cList[indexI]];
1202                     scalar sumLength = 0.0;
1203                     label nTouchedNgb = 0;
1204
1205                     forAll(ncList, indexJ)
1206                     {
1207                         label sLevel = cellLevels[ncList[indexJ]];
1208
1209                         if ((sLevel < ngbLevel) && (sLevel > 0))
1210                         {
1211                             sumLength += lengthScale[ncList[indexJ]];
1212
1213                             nTouchedNgb++;
1214                         }
1215                     }
1216
1217                     sumLength /= nTouchedNgb;
1218
1219                     // Scale the length and assign to this cell
1220                     if (level < maxRefineLevel_)
1221                     {
1222                         sumLength *= growthFactor_;
1223                     }
1224                     else
1225                     if (meanScale_ > 0.0)
1226                     {
1227                         // If a mean scale has been specified,
1228                         // override the value
1229                         sumLength = meanScale_;
1230                     }
1231
1232                     lengthScale[cList[indexI]] = sumLength;
1233
1234                     levelCells.insert(cList[indexI]);
1235
1236                     visitedCells++;
1237                 }
1238             }
1239         }
1240
1241         if (Pstream::parRun())
1242         {
1243             readLengthScaleInfo
1244             (
1245                 level,
1246                 visitedCells,
1247                 cellLevels,
1248                 lengthScale,
1249                 levelCells
1250             );
1251         }
1252
1253         if (debug > 2)
1254         {
1255             Pout << "Processed level: " << level << nl
1256                  << " Visited: " << visitedCells
1257                  << " out of " << mesh_.nCells() << endl;
1258         }
1259
1260         // Move on to the next level
1261         level++;
1262
1263         if (visitedCells >= mesh_.nCells())
1264         {
1265             doneWithSweeps = true;
1266         }
1267
1268         // Wait for everyone to complete.
1269         reduce(doneWithSweeps, andOp<bool>());
1270     }
1271
1272     if (debug)
1273     {
1274         Info << "Max Length Scale: " << maxLengthScale_ << endl;
1275         Info << "Length Scale sweeps: " << level << endl;
1276     }
1277
1278     // Check if everything went okay
1279     if (visitedCells != mesh_.nCells())
1280     {
1281         FatalErrorIn
1282         (
1283             "void lengthScaleEstimator::calculateLengthScale"
1284             "(UList<scalar>& lengthScale)"
1285         )
1286             << " Algorithm did not visit every cell in the mesh."
1287             << " Something's messed up." << nl
1288             << " Visited cells: " << visitedCells
1289             << " nCells: " << mesh_.nCells()
1290             << abort(FatalError);
1291     }
1292
1293     // Wait for transfers before continuing.
1294     if (Pstream::parRun())
1295     {
1296         OPstream::waitRequests();
1297         IPstream::waitRequests();
1298     }
1299 }
1300
1301
1302 } // End namespace Foam
1303
1304 // ************************************************************************* //