src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "dynamicRefineFvMesh.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "surfaceInterpolate.H"
  29 #include "volFields.H"
  30 #include "polyTopoChange.H"
  31 #include "surfaceFields.H"
  32 #include "syncTools.H"
  33 #include "pointFields.H"
  34
  35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  36
  37 namespace Foam
  38 {
  39     defineTypeNameAndDebug(dynamicRefineFvMesh, 0);
  40     addToRunTimeSelectionTable(dynamicFvMesh, dynamicRefineFvMesh, IOobject);
  41 }
  42
  43 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  44
  45 // the PackedBoolList::count method would probably be faster
  46 // since we are only checking for 'true' anyhow
  47 Foam::label Foam::dynamicRefineFvMesh::count
  48 (
  49     const PackedBoolList& l,
  50     const unsigned int val
  51 )
  52 {
  53     label n = 0;
  54     forAll(l, i)
  55     {
  56         if (l.get(i) == val)
  57         {
  58             n++;
  59         }
  60     }
  61     return n;
  62 }
  63
  64
  65 void Foam::dynamicRefineFvMesh::calculateProtectedCells
  66 (
  67     PackedBoolList& unrefineableCell
  68 ) const
  69 {
  70     if (protectedCell_.empty())
  71     {
  72         unrefineableCell.clear();
  73         return;
  74     }
  75
  76     const labelList& cellLevel = meshCutter_.cellLevel();
  77
  78     unrefineableCell = protectedCell_;
  79
  80     // Get neighbouring cell level
  81     labelList neiLevel(nFaces()-nInternalFaces());
  82
  83     for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
  84     {
  85         neiLevel[faceI-nInternalFaces()] = cellLevel[faceOwner()[faceI]];
  86     }
  87     syncTools::swapBoundaryFaceList(*this, neiLevel);
  88
  89
  90     while (true)
  91     {
  92         // Pick up faces on border of protected cells
  93         boolList seedFace(nFaces(), false);
  94
  95         forAll(faceNeighbour(), faceI)
  96         {
  97             label own = faceOwner()[faceI];
  98             bool ownProtected = unrefineableCell.get(own);
  99             label nei = faceNeighbour()[faceI];
 100             bool neiProtected = unrefineableCell.get(nei);
 101
 102             if (ownProtected && (cellLevel[nei] > cellLevel[own]))
 103             {
 104                 seedFace[faceI] = true;
 105             }
 106             else if (neiProtected && (cellLevel[own] > cellLevel[nei]))
 107             {
 108                 seedFace[faceI] = true;
 109             }
 110         }
 111         for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
 112         {
 113             label own = faceOwner()[faceI];
 114             bool ownProtected = unrefineableCell.get(own);
 115             if
 116             (
 117                 ownProtected
 118              && (neiLevel[faceI-nInternalFaces()] > cellLevel[own])
 119             )
 120             {
 121                 seedFace[faceI] = true;
 122             }
 123         }
 124
 125         syncTools::syncFaceList(*this, seedFace, orEqOp<bool>());
 126
 127
 128         // Extend unrefineableCell
 129         bool hasExtended = false;
 130
 131         for (label faceI = 0; faceI < nInternalFaces(); faceI++)
 132         {
 133             if (seedFace[faceI])
 134             {
 135                 label own = faceOwner()[faceI];
 136                 if (unrefineableCell.get(own) == 0)
 137                 {
 138                     unrefineableCell.set(own, 1);
 139                     hasExtended = true;
 140                 }
 141
 142                 label nei = faceNeighbour()[faceI];
 143                 if (unrefineableCell.get(nei) == 0)
 144                 {
 145                     unrefineableCell.set(nei, 1);
 146                     hasExtended = true;
 147                 }
 148             }
 149         }
 150         for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
 151         {
 152             if (seedFace[faceI])
 153             {
 154                 label own = faceOwner()[faceI];
 155                 if (unrefineableCell.get(own) == 0)
 156                 {
 157                     unrefineableCell.set(own, 1);
 158                     hasExtended = true;
 159                 }
 160             }
 161         }
 162
 163         if (!returnReduce(hasExtended, orOp<bool>()))
 164         {
 165             break;
 166         }
 167     }
 168 }
 169
 170
 171 void Foam::dynamicRefineFvMesh::readDict()
 172 {
 173     dictionary refineDict
 174     (
 175         IOdictionary
 176         (
 177             IOobject
 178             (
 179                 "dynamicMeshDict",
 180                 time().constant(),
 181                 *this,
 182                 IOobject::MUST_READ_IF_MODIFIED,
 183                 IOobject::NO_WRITE,
 184                 false
 185             )
 186         ).subDict(typeName + "Coeffs")
 187     );
 188
 189     List<Pair<word> > fluxVelocities = List<Pair<word> >
 190     (
 191         refineDict.lookup("correctFluxes")
 192     );
 193     // Rework into hashtable.
 194     correctFluxes_.resize(fluxVelocities.size());
 195     forAll(fluxVelocities, i)
 196     {
 197         correctFluxes_.insert(fluxVelocities[i][0], fluxVelocities[i][1]);
 198     }
 199
 200     dumpLevel_ = Switch(refineDict.lookup("dumpLevel"));
 201 }
 202
 203
 204 // Refines cells, maps fields and recalculates (an approximate) flux
 205 Foam::autoPtr<Foam::mapPolyMesh>
 206 Foam::dynamicRefineFvMesh::refine
 207 (
 208     const labelList& cellsToRefine
 209 )
 210 {
 211     // Mesh changing engine.
 212     polyTopoChange meshMod(*this);
 213
 214     // Play refinement commands into mesh changer.
 215     meshCutter_.setRefinement(cellsToRefine, meshMod);
 216
 217     // Create mesh (with inflation), return map from old to new mesh.
 218     //autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, true);
 219     autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, false);
 220
 221     Info<< "Refined from "
 222         << returnReduce(map().nOldCells(), sumOp<label>())
 223         << " to " << globalData().nTotalCells() << " cells." << endl;
 224
 225     if (debug)
 226     {
 227         // Check map.
 228         for (label faceI = 0; faceI < nInternalFaces(); faceI++)
 229         {
 230             label oldFaceI = map().faceMap()[faceI];
 231
 232             if (oldFaceI >= nInternalFaces())
 233             {
 234                 FatalErrorIn("dynamicRefineFvMesh::refine(const labelList&)")
 235                     << "New internal face:" << faceI
 236                     << " fc:" << faceCentres()[faceI]
 237                     << " originates from boundary oldFace:" << oldFaceI
 238                     << abort(FatalError);
 239             }
 240         }
 241     }
 242
 243     // Update fields
 244     updateMesh(map);
 245     // Move mesh
 246     /*
 247     pointField newPoints;
 248     if (map().hasMotionPoints())
 249     {
 250         newPoints = map().preMotionPoints();
 251     }
 252     else
 253     {
 254         newPoints = points();
 255     }
 256     movePoints(newPoints);
 257     */
 258
 259     // Correct the flux for modified/added faces. All the faces which only
 260     // have been renumbered will already have been handled by the mapping.
 261     {
 262         const labelList& faceMap = map().faceMap();
 263         const labelList& reverseFaceMap = map().reverseFaceMap();
 264
 265         // Storage for any master faces. These will be the original faces
 266         // on the coarse cell that get split into four (or rather the
 267         // master face gets modified and three faces get added from the master)
 268         labelHashSet masterFaces(4*cellsToRefine.size());
 269
 270         forAll(faceMap, faceI)
 271         {
 272             label oldFaceI = faceMap[faceI];
 273
 274             if (oldFaceI >= 0)
 275             {
 276                 label masterFaceI = reverseFaceMap[oldFaceI];
 277
 278                 if (masterFaceI < 0)
 279                 {
 280                     FatalErrorIn
 281                     (
 282                         "dynamicRefineFvMesh::refine(const labelList&)"
 283                     )   << "Problem: should not have removed faces"
 284                         << " when refining."
 285                         << nl << "face:" << faceI << abort(FatalError);
 286                 }
 287                 else if (masterFaceI != faceI)
 288                 {
 289                     masterFaces.insert(masterFaceI);
 290                 }
 291             }
 292         }
 293         if (debug)
 294         {
 295             Info<< "Found " << returnReduce(masterFaces.size(), sumOp<label>())
 296                 << " split faces " << endl;
 297         }
 298
 299         HashTable<const surfaceScalarField*> fluxes
 300         (
 301             lookupClass<surfaceScalarField>()
 302         );
 303         forAllConstIter(HashTable<const surfaceScalarField*>, fluxes, iter)
 304         {
 305             if (!correctFluxes_.found(iter.key()))
 306             {
 307                 WarningIn("dynamicRefineFvMesh::refine(const labelList&)")
 308                     << "Cannot find surfaceScalarField " << iter.key()
 309                     << " in user-provided flux mapping table "
 310                     << correctFluxes_ << endl
 311                     << "    The flux mapping table is used to recreate the"
 312                     << " flux on newly created faces." << endl
 313                     << "    Either add the entry if it is a flux or use ("
 314                     << iter.key() << " none) to suppress this warning."
 315                     << endl;
 316                 continue;
 317             }
 318
 319             const word& UName = correctFluxes_[iter.key()];
 320
 321             if (UName == "none")
 322             {
 323                 continue;
 324             }
 325
 326             if (debug)
 327             {
 328                 Info<< "Mapping flux " << iter.key()
 329                     << " using interpolated flux " << UName
 330                     << endl;
 331             }
 332
 333             surfaceScalarField& phi = const_cast<surfaceScalarField&>(*iter());
 334
 335             const surfaceScalarField phiU
 336             (
 337                 fvc::interpolate
 338                 (
 339                     lookupObject<volVectorField>(UName)
 340                 )
 341               & Sf()
 342             );
 343
 344             // Recalculate new internal faces.
 345             for (label faceI = 0; faceI < nInternalFaces(); faceI++)
 346             {
 347                 label oldFaceI = faceMap[faceI];
 348
 349                 if (oldFaceI == -1)
 350                 {
 351                     // Inflated/appended
 352                     phi[faceI] = phiU[faceI];
 353                 }
 354                 else if (reverseFaceMap[oldFaceI] != faceI)
 355                 {
 356                     // face-from-masterface
 357                     phi[faceI] = phiU[faceI];
 358                 }
 359             }
 360
 361             // Recalculate new boundary faces.
 362             forAll(phi.boundaryField(), patchI)
 363             {
 364                 fvsPatchScalarField& patchPhi = phi.boundaryField()[patchI];
 365                 const fvsPatchScalarField& patchPhiU =
 366                     phiU.boundaryField()[patchI];
 367
 368                 label faceI = patchPhi.patch().start();
 369
 370                 forAll(patchPhi, i)
 371                 {
 372                     label oldFaceI = faceMap[faceI];
 373
 374                     if (oldFaceI == -1)
 375                     {
 376                         // Inflated/appended
 377                         patchPhi[i] = patchPhiU[i];
 378                     }
 379                     else if (reverseFaceMap[oldFaceI] != faceI)
 380                     {
 381                         // face-from-masterface
 382                         patchPhi[i] = patchPhiU[i];
 383                     }
 384
 385                     faceI++;
 386                 }
 387             }
 388
 389             // Update master faces
 390             forAllConstIter(labelHashSet, masterFaces, iter)
 391             {
 392                 label faceI = iter.key();
 393
 394                 if (isInternalFace(faceI))
 395                 {
 396                     phi[faceI] = phiU[faceI];
 397                 }
 398                 else
 399                 {
 400                     label patchI = boundaryMesh().whichPatch(faceI);
 401                     label i = faceI - boundaryMesh()[patchI].start();
 402
 403                     const fvsPatchScalarField& patchPhiU =
 404                         phiU.boundaryField()[patchI];
 405
 406                     fvsPatchScalarField& patchPhi =
 407                         phi.boundaryField()[patchI];
 408
 409                     if (patchPhi.size() > 0)
 410                     {
 411                         patchPhi[i] = patchPhiU[i];
 412                     }
 413                 }
 414             }
 415         }
 416     }
 417
 418     // Update numbering of cells/vertices.
 419     meshCutter_.updateMesh(map);
 420
 421     // Update numbering of protectedCell_
 422     if (protectedCell_.size())
 423     {
 424         PackedBoolList newProtectedCell(nCells());
 425
 426         forAll(newProtectedCell, cellI)
 427         {
 428             label oldCellI = map().cellMap()[cellI];
 429             newProtectedCell.set(cellI, protectedCell_.get(oldCellI));
 430         }
 431         protectedCell_.transfer(newProtectedCell);
 432     }
 433
 434     // Debug: Check refinement levels (across faces only)
 435     meshCutter_.checkRefinementLevels(-1, labelList(0));
 436
 437     return map;
 438 }
 439
 440
 441 // Combines previously split cells, maps fields and recalculates
 442 // (an approximate) flux
 443 Foam::autoPtr<Foam::mapPolyMesh>
 444 Foam::dynamicRefineFvMesh::unrefine
 445 (
 446     const labelList& splitPoints
 447 )
 448 {
 449     polyTopoChange meshMod(*this);
 450
 451     // Play refinement commands into mesh changer.
 452     meshCutter_.setUnrefinement(splitPoints, meshMod);
 453
 454
 455     // Save information on faces that will be combined
 456     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 457
 458     // Find the faceMidPoints on cells to be combined.
 459     // for each face resulting of split of face into four store the
 460     // midpoint
 461     Map<label> faceToSplitPoint(3*splitPoints.size());
 462
 463     {
 464         forAll(splitPoints, i)
 465         {
 466             label pointI = splitPoints[i];
 467
 468             const labelList& pEdges = pointEdges()[pointI];
 469
 470             forAll(pEdges, j)
 471             {
 472                 label otherPointI = edges()[pEdges[j]].otherVertex(pointI);
 473
 474                 const labelList& pFaces = pointFaces()[otherPointI];
 475
 476                 forAll(pFaces, pFaceI)
 477                 {
 478                     faceToSplitPoint.insert(pFaces[pFaceI], otherPointI);
 479                 }
 480             }
 481         }
 482     }
 483
 484
 485     // Change mesh and generate map.
 486     //autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, true);
 487     autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, false);
 488
 489     Info<< "Unrefined from "
 490         << returnReduce(map().nOldCells(), sumOp<label>())
 491         << " to " << globalData().nTotalCells() << " cells."
 492         << endl;
 493
 494     // Update fields
 495     updateMesh(map);
 496
 497
 498     // Move mesh
 499     /*
 500     pointField newPoints;
 501     if (map().hasMotionPoints())
 502     {
 503         newPoints = map().preMotionPoints();
 504     }
 505     else
 506     {
 507         newPoints = points();
 508     }
 509     movePoints(newPoints);
 510     */
 511
 512     // Correct the flux for modified faces.
 513     {
 514         const labelList& reversePointMap = map().reversePointMap();
 515         const labelList& reverseFaceMap = map().reverseFaceMap();
 516
 517         HashTable<const surfaceScalarField*> fluxes
 518         (
 519             lookupClass<surfaceScalarField>()
 520         );
 521         forAllConstIter(HashTable<const surfaceScalarField*>, fluxes, iter)
 522         {
 523             if (!correctFluxes_.found(iter.key()))
 524             {
 525                 WarningIn("dynamicRefineFvMesh::refine(const labelList&)")
 526                     << "Cannot find surfaceScalarField " << iter.key()
 527                     << " in user-provided flux mapping table "
 528                     << correctFluxes_ << endl
 529                     << "    The flux mapping table is used to recreate the"
 530                     << " flux on newly created faces." << endl
 531                     << "    Either add the entry if it is a flux or use ("
 532                     << iter.key() << " none) to suppress this warning."
 533                     << endl;
 534                 continue;
 535             }
 536
 537             const word& UName = correctFluxes_[iter.key()];
 538
 539             if (UName == "none")
 540             {
 541                 continue;
 542             }
 543
 544             if (debug)
 545             {
 546                 Info<< "Mapping flux " << iter.key()
 547                     << " using interpolated flux " << UName
 548                     << endl;
 549             }
 550
 551             surfaceScalarField& phi = const_cast<surfaceScalarField&>(*iter());
 552             const surfaceScalarField phiU
 553             (
 554                 fvc::interpolate
 555                 (
 556                     lookupObject<volVectorField>(UName)
 557                 )
 558               & Sf()
 559             );
 560
 561
 562             forAllConstIter(Map<label>, faceToSplitPoint, iter)
 563             {
 564                 label oldFaceI = iter.key();
 565                 label oldPointI = iter();
 566
 567                 if (reversePointMap[oldPointI] < 0)
 568                 {
 569                     // midpoint was removed. See if face still exists.
 570                     label faceI = reverseFaceMap[oldFaceI];
 571
 572                     if (faceI >= 0)
 573                     {
 574                         if (isInternalFace(faceI))
 575                         {
 576                             phi[faceI] = phiU[faceI];
 577                         }
 578                         else
 579                         {
 580                             label patchI = boundaryMesh().whichPatch(faceI);
 581                             label i = faceI - boundaryMesh()[patchI].start();
 582
 583                             const fvsPatchScalarField& patchPhiU =
 584                                 phiU.boundaryField()[patchI];
 585
 586                             fvsPatchScalarField& patchPhi =
 587                                 phi.boundaryField()[patchI];
 588
 589                             patchPhi[i] = patchPhiU[i];
 590                         }
 591                     }
 592                 }
 593             }
 594         }
 595     }
 596
 597
 598     // Update numbering of cells/vertices.
 599     meshCutter_.updateMesh(map);
 600
 601     // Update numbering of protectedCell_
 602     if (protectedCell_.size())
 603     {
 604         PackedBoolList newProtectedCell(nCells());
 605
 606         forAll(newProtectedCell, cellI)
 607         {
 608             label oldCellI = map().cellMap()[cellI];
 609             if (oldCellI >= 0)
 610             {
 611                 newProtectedCell.set(cellI, protectedCell_.get(oldCellI));
 612             }
 613         }
 614         protectedCell_.transfer(newProtectedCell);
 615     }
 616
 617     // Debug: Check refinement levels (across faces only)
 618     meshCutter_.checkRefinementLevels(-1, labelList(0));
 619
 620     return map;
 621 }
 622
 623
 624 // Get max of connected point
 625 Foam::scalarField
 626 Foam::dynamicRefineFvMesh::maxPointField(const scalarField& pFld) const
 627 {
 628     scalarField vFld(nCells(), -GREAT);
 629
 630     forAll(pointCells(), pointI)
 631     {
 632         const labelList& pCells = pointCells()[pointI];
 633
 634         forAll(pCells, i)
 635         {
 636             vFld[pCells[i]] = max(vFld[pCells[i]], pFld[pointI]);
 637         }
 638     }
 639     return vFld;
 640 }
 641
 642
 643 // Get min of connected cell
 644 Foam::scalarField
 645 Foam::dynamicRefineFvMesh::minCellField(const volScalarField& vFld) const
 646 {
 647     scalarField pFld(nPoints(), GREAT);
 648
 649     forAll(pointCells(), pointI)
 650     {
 651         const labelList& pCells = pointCells()[pointI];
 652
 653         forAll(pCells, i)
 654         {
 655             pFld[pointI] = min(pFld[pointI], vFld[pCells[i]]);
 656         }
 657     }
 658     return pFld;
 659 }
 660
 661
 662 // Simple (non-parallel) interpolation by averaging.
 663 Foam::scalarField
 664 Foam::dynamicRefineFvMesh::cellToPoint(const scalarField& vFld) const
 665 {
 666     scalarField pFld(nPoints());
 667
 668     forAll(pointCells(), pointI)
 669     {
 670         const labelList& pCells = pointCells()[pointI];
 671
 672         scalar sum = 0.0;
 673         forAll(pCells, i)
 674         {
 675             sum += vFld[pCells[i]];
 676         }
 677         pFld[pointI] = sum/pCells.size();
 678     }
 679     return pFld;
 680 }
 681
 682
 683 // Calculate error. Is < 0 or distance to minLevel, maxLevel
 684 Foam::scalarField Foam::dynamicRefineFvMesh::error
 685 (
 686     const scalarField& fld,
 687     const scalar minLevel,
 688     const scalar maxLevel
 689 ) const
 690 {
 691     scalarField c(fld.size(), -1);
 692
 693     forAll(fld, i)
 694     {
 695         scalar err = min(fld[i]-minLevel, maxLevel-fld[i]);
 696
 697         if (err >= 0)
 698         {
 699             c[i] = err;
 700         }
 701     }
 702     return c;
 703 }
 704
 705
 706 void Foam::dynamicRefineFvMesh::selectRefineCandidates
 707 (
 708     const scalar lowerRefineLevel,
 709     const scalar upperRefineLevel,
 710     const scalarField& vFld,
 711     PackedBoolList& candidateCell
 712 ) const
 713 {
 714     // Get error per cell. Is -1 (not to be refined) to >0 (to be refined,
 715     // higher more desirable to be refined).
 716     scalarField cellError
 717     (
 718         maxPointField
 719         (
 720             error
 721             (
 722                 cellToPoint(vFld),
 723                 lowerRefineLevel,
 724                 upperRefineLevel
 725             )
 726         )
 727     );
 728
 729     // Mark cells that are candidates for refinement.
 730     forAll(cellError, cellI)
 731     {
 732         if (cellError[cellI] > 0)
 733         {
 734             candidateCell.set(cellI, 1);
 735         }
 736     }
 737 }
 738
 739
 740 Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
 741 (
 742     const label maxCells,
 743     const label maxRefinement,
 744     const PackedBoolList& candidateCell
 745 ) const
 746 {
 747     // Every refined cell causes 7 extra cells
 748     label nTotToRefine = (maxCells - globalData().nTotalCells()) / 7;
 749
 750     const labelList& cellLevel = meshCutter_.cellLevel();
 751
 752     // Mark cells that cannot be refined since they would trigger refinement
 753     // of protected cells (since 2:1 cascade)
 754     PackedBoolList unrefineableCell;
 755     calculateProtectedCells(unrefineableCell);
 756
 757     // Count current selection
 758     label nCandidates = returnReduce(count(candidateCell, 1), sumOp<label>());
 759
 760     // Collect all cells
 761     DynamicList<label> candidates(nCells());
 762
 763     if (nCandidates < nTotToRefine)
 764     {
 765         forAll(candidateCell, cellI)
 766         {
 767             if
 768             (
 769                 cellLevel[cellI] < maxRefinement
 770              && candidateCell.get(cellI)
 771              && (
 772                     unrefineableCell.empty()
 773                  || !unrefineableCell.get(cellI)
 774                 )
 775             )
 776             {
 777                 candidates.append(cellI);
 778             }
 779         }
 780     }
 781     else
 782     {
 783         // Sort by error? For now just truncate.
 784         for (label level = 0; level < maxRefinement; level++)
 785         {
 786             forAll(candidateCell, cellI)
 787             {
 788                 if
 789                 (
 790                     cellLevel[cellI] == level
 791                  && candidateCell.get(cellI)
 792                  && (
 793                         unrefineableCell.empty()
 794                      || !unrefineableCell.get(cellI)
 795                     )
 796                 )
 797                 {
 798                     candidates.append(cellI);
 799                 }
 800             }
 801
 802             if (returnReduce(candidates.size(), sumOp<label>()) > nTotToRefine)
 803             {
 804                 break;
 805             }
 806         }
 807     }
 808
 809     // Guarantee 2:1 refinement after refinement
 810     labelList consistentSet
 811     (
 812         meshCutter_.consistentRefinement
 813         (
 814             candidates.shrink(),
 815             true               // Add to set to guarantee 2:1
 816         )
 817     );
 818
 819     Info<< "Selected " << returnReduce(consistentSet.size(), sumOp<label>())
 820         << " cells for refinement out of " << globalData().nTotalCells()
 821         << "." << endl;
 822
 823     return consistentSet;
 824 }
 825
 826
 827 Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
 828 (
 829     const scalar unrefineLevel,
 830     const PackedBoolList& markedCell,
 831     const scalarField& pFld
 832 ) const
 833 {
 834     // All points that can be unrefined
 835     const labelList splitPoints(meshCutter_.getSplitPoints());
 836
 837     DynamicList<label> newSplitPoints(splitPoints.size());
 838
 839     forAll(splitPoints, i)
 840     {
 841         label pointI = splitPoints[i];
 842
 843         if (pFld[pointI] < unrefineLevel)
 844         {
 845             // Check that all cells are not marked
 846             const labelList& pCells = pointCells()[pointI];
 847
 848             bool hasMarked = false;
 849
 850             forAll(pCells, pCellI)
 851             {
 852                 if (markedCell.get(pCells[pCellI]))
 853                 {
 854                     hasMarked = true;
 855                     break;
 856                 }
 857             }
 858
 859             if (!hasMarked)
 860             {
 861                 newSplitPoints.append(pointI);
 862             }
 863         }
 864     }
 865
 866
 867     newSplitPoints.shrink();
 868
 869     // Guarantee 2:1 refinement after unrefinement
 870     labelList consistentSet
 871     (
 872         meshCutter_.consistentUnrefinement
 873         (
 874             newSplitPoints,
 875             false
 876         )
 877     );
 878     Info<< "Selected " << returnReduce(consistentSet.size(), sumOp<label>())
 879         << " split points out of a possible "
 880         << returnReduce(splitPoints.size(), sumOp<label>())
 881         << "." << endl;
 882
 883     return consistentSet;
 884 }
 885
 886
 887 void Foam::dynamicRefineFvMesh::extendMarkedCells
 888 (
 889     PackedBoolList& markedCell
 890 ) const
 891 {
 892     // Mark faces using any marked cell
 893     boolList markedFace(nFaces(), false);
 894
 895     forAll(markedCell, cellI)
 896     {
 897         if (markedCell.get(cellI))
 898         {
 899             const cell& cFaces = cells()[cellI];
 900
 901             forAll(cFaces, i)
 902             {
 903                 markedFace[cFaces[i]] = true;
 904             }
 905         }
 906     }
 907
 908     syncTools::syncFaceList(*this, markedFace, orEqOp<bool>());
 909
 910     // Update cells using any markedFace
 911     for (label faceI = 0; faceI < nInternalFaces(); faceI++)
 912     {
 913         if (markedFace[faceI])
 914         {
 915             markedCell.set(faceOwner()[faceI], 1);
 916             markedCell.set(faceNeighbour()[faceI], 1);
 917         }
 918     }
 919     for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
 920     {
 921         if (markedFace[faceI])
 922         {
 923             markedCell.set(faceOwner()[faceI], 1);
 924         }
 925     }
 926 }
 927
 928
 929 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 930
 931 Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
 932 :
 933     dynamicFvMesh(io),
 934     meshCutter_(*this),
 935     dumpLevel_(false),
 936     nRefinementIterations_(0),
 937     protectedCell_(nCells(), 0)
 938 {
 939     // Read static part of dictionary
 940     readDict();
 941
 942
 943     const labelList& cellLevel = meshCutter_.cellLevel();
 944     const labelList& pointLevel = meshCutter_.pointLevel();
 945
 946     // Set cells that should not be refined.
 947     // This is currently any cell which does not have 8 anchor points or
 948     // uses any face which does not have 4 anchor points.
 949     // Note: do not use cellPoint addressing
 950
 951     // Count number of points <= cellLevel
 952     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 953
 954     labelList nAnchors(nCells(), 0);
 955
 956     label nProtected = 0;
 957
 958     forAll(pointCells(), pointI)
 959     {
 960         const labelList& pCells = pointCells()[pointI];
 961
 962         forAll(pCells, i)
 963         {
 964             label cellI = pCells[i];
 965
 966             if (!protectedCell_.get(cellI))
 967             {
 968                 if (pointLevel[pointI] <= cellLevel[cellI])
 969                 {
 970                     nAnchors[cellI]++;
 971
 972                     if (nAnchors[cellI] > 8)
 973                     {
 974                         protectedCell_.set(cellI, 1);
 975                         nProtected++;
 976                     }
 977                 }
 978             }
 979         }
 980     }
 981
 982
 983     // Count number of points <= faceLevel
 984     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 985     // Bit tricky since proc face might be one more refined than the owner since
 986     // the coupled one is refined.
 987
 988     {
 989         labelList neiLevel(nFaces());
 990
 991         for (label faceI = 0; faceI < nInternalFaces(); faceI++)
 992         {
 993             neiLevel[faceI] = cellLevel[faceNeighbour()[faceI]];
 994         }
 995         for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
 996         {
 997             neiLevel[faceI] = cellLevel[faceOwner()[faceI]];
 998         }
 999         syncTools::swapFaceList(*this, neiLevel);
1000
1001
1002         boolList protectedFace(nFaces(), false);
1003
1004         forAll(faceOwner(), faceI)
1005         {
1006             label faceLevel = max
1007             (
1008                 cellLevel[faceOwner()[faceI]],
1009                 neiLevel[faceI]
1010             );
1011
1012             const face& f = faces()[faceI];
1013
1014             label nAnchors = 0;
1015
1016             forAll(f, fp)
1017             {
1018                 if (pointLevel[f[fp]] <= faceLevel)
1019                 {
1020                     nAnchors++;
1021
1022                     if (nAnchors > 4)
1023                     {
1024                         protectedFace[faceI] = true;
1025                         break;
1026                     }
1027                 }
1028             }
1029         }
1030
1031         syncTools::syncFaceList(*this, protectedFace, orEqOp<bool>());
1032
1033         for (label faceI = 0; faceI < nInternalFaces(); faceI++)
1034         {
1035             if (protectedFace[faceI])
1036             {
1037                 protectedCell_.set(faceOwner()[faceI], 1);
1038                 nProtected++;
1039                 protectedCell_.set(faceNeighbour()[faceI], 1);
1040                 nProtected++;
1041             }
1042         }
1043         for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
1044         {
1045             if (protectedFace[faceI])
1046             {
1047                 protectedCell_.set(faceOwner()[faceI], 1);
1048                 nProtected++;
1049             }
1050         }
1051     }
1052
1053     if (returnReduce(nProtected, sumOp<label>()) == 0)
1054     {
1055         protectedCell_.clear();
1056     }
1057 }
1058
1059
1060 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
1061
1062 Foam::dynamicRefineFvMesh::~dynamicRefineFvMesh()
1063 {}
1064
1065
1066 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
1067
1068 bool Foam::dynamicRefineFvMesh::update()
1069 {
1070     // Re-read dictionary. Choosen since usually -small so trivial amount
1071     // of time compared to actual refinement. Also very useful to be able
1072     // to modify on-the-fly.
1073     dictionary refineDict
1074     (
1075         IOdictionary
1076         (
1077             IOobject
1078             (
1079                 "dynamicMeshDict",
1080                 time().constant(),
1081                 *this,
1082                 IOobject::MUST_READ_IF_MODIFIED,
1083                 IOobject::NO_WRITE,
1084                 false
1085             )
1086         ).subDict(typeName + "Coeffs")
1087     );
1088
1089     label refineInterval = readLabel(refineDict.lookup("refineInterval"));
1090
1091     bool hasChanged = false;
1092
1093     if (refineInterval == 0)
1094     {
1095         changing(hasChanged);
1096
1097         return false;
1098     }
1099     else if (refineInterval < 0)
1100     {
1101         FatalErrorIn("dynamicRefineFvMesh::update()")
1102             << "Illegal refineInterval " << refineInterval << nl
1103             << "The refineInterval setting in the dynamicMeshDict should"
1104             << " be >= 1." << nl
1105             << exit(FatalError);
1106     }
1107
1108
1109
1110
1111     // Note: cannot refine at time 0 since no V0 present since mesh not
1112     //       moved yet.
1113
1114     if (time().timeIndex() > 0 && time().timeIndex() % refineInterval == 0)
1115     {
1116         label maxCells = readLabel(refineDict.lookup("maxCells"));
1117
1118         if (maxCells <= 0)
1119         {
1120             FatalErrorIn("dynamicRefineFvMesh::update()")
1121                 << "Illegal maximum number of cells " << maxCells << nl
1122                 << "The maxCells setting in the dynamicMeshDict should"
1123                 << " be > 0." << nl
1124                 << exit(FatalError);
1125         }
1126
1127         label maxRefinement = readLabel(refineDict.lookup("maxRefinement"));
1128
1129         if (maxRefinement <= 0)
1130         {
1131             FatalErrorIn("dynamicRefineFvMesh::update()")
1132                 << "Illegal maximum refinement level " << maxRefinement << nl
1133                 << "The maxCells setting in the dynamicMeshDict should"
1134                 << " be > 0." << nl
1135                 << exit(FatalError);
1136         }
1137
1138         const word fieldName(refineDict.lookup("field"));
1139
1140         const volScalarField& vFld = lookupObject<volScalarField>(fieldName);
1141
1142         const scalar lowerRefineLevel =
1143             readScalar(refineDict.lookup("lowerRefineLevel"));
1144         const scalar upperRefineLevel =
1145             readScalar(refineDict.lookup("upperRefineLevel"));
1146         const scalar unrefineLevel =
1147             readScalar(refineDict.lookup("unrefineLevel"));
1148         const label nBufferLayers =
1149             readLabel(refineDict.lookup("nBufferLayers"));
1150
1151         // Cells marked for refinement or otherwise protected from unrefinement.
1152         PackedBoolList refineCell(nCells());
1153
1154         if (globalData().nTotalCells() < maxCells)
1155         {
1156             // Determine candidates for refinement (looking at field only)
1157             selectRefineCandidates
1158             (
1159                 lowerRefineLevel,
1160                 upperRefineLevel,
1161                 vFld,
1162                 refineCell
1163             );
1164
1165             // Select subset of candidates. Take into account max allowable
1166             // cells, refinement level, protected cells.
1167             labelList cellsToRefine
1168             (
1169                 selectRefineCells
1170                 (
1171                     maxCells,
1172                     maxRefinement,
1173                     refineCell
1174                 )
1175             );
1176
1177             label nCellsToRefine = returnReduce
1178             (
1179                 cellsToRefine.size(), sumOp<label>()
1180             );
1181
1182             if (nCellsToRefine > 0)
1183             {
1184                 // Refine/update mesh and map fields
1185                 autoPtr<mapPolyMesh> map = refine(cellsToRefine);
1186
1187                 // Update refineCell. Note that some of the marked ones have
1188                 // not been refined due to constraints.
1189                 {
1190                     const labelList& cellMap = map().cellMap();
1191                     const labelList& reverseCellMap = map().reverseCellMap();
1192
1193                     PackedBoolList newRefineCell(cellMap.size());
1194
1195                     forAll(cellMap, cellI)
1196                     {
1197                         label oldCellI = cellMap[cellI];
1198
1199                         if (oldCellI < 0)
1200                         {
1201                             newRefineCell.set(cellI, 1);
1202                         }
1203                         else if (reverseCellMap[oldCellI] != cellI)
1204                         {
1205                             newRefineCell.set(cellI, 1);
1206                         }
1207                         else
1208                         {
1209                             newRefineCell.set(cellI, refineCell.get(oldCellI));
1210                         }
1211                     }
1212                     refineCell.transfer(newRefineCell);
1213                 }
1214
1215                 // Extend with a buffer layer to prevent neighbouring points
1216                 // being unrefined.
1217                 for (label i = 0; i < nBufferLayers; i++)
1218                 {
1219                     extendMarkedCells(refineCell);
1220                 }
1221
1222                 hasChanged = true;
1223             }
1224         }
1225
1226
1227         {
1228             // Select unrefineable points that are not marked in refineCell
1229             labelList pointsToUnrefine
1230             (
1231                 selectUnrefinePoints
1232                 (
1233                     unrefineLevel,
1234                     refineCell,
1235                     minCellField(vFld)
1236                 )
1237             );
1238
1239             label nSplitPoints = returnReduce
1240             (
1241                 pointsToUnrefine.size(),
1242                 sumOp<label>()
1243             );
1244
1245             if (nSplitPoints > 0)
1246             {
1247                 // Refine/update mesh
1248                 unrefine(pointsToUnrefine);
1249
1250                 hasChanged = true;
1251             }
1252         }
1253
1254
1255         if ((nRefinementIterations_ % 10) == 0)
1256         {
1257             // Compact refinement history occassionally (how often?).
1258             // Unrefinement causes holes in the refinementHistory.
1259             const_cast<refinementHistory&>(meshCutter().history()).compact();
1260         }
1261         nRefinementIterations_++;
1262     }
1263
1264     changing(hasChanged);
1265
1266     return hasChanged;
1267 }
1268
1269
1270 bool Foam::dynamicRefineFvMesh::writeObject
1271 (
1272     IOstream::streamFormat fmt,
1273     IOstream::versionNumber ver,
1274     IOstream::compressionType cmp
1275 ) const
1276 {
1277     // Force refinement data to go to the current time directory.
1278     const_cast<hexRef8&>(meshCutter_).setInstance(time().timeName());
1279
1280     bool writeOk =
1281     (
1282         dynamicFvMesh::writeObjects(fmt, ver, cmp)
1283      && meshCutter_.write()
1284     );
1285
1286     if (dumpLevel_)
1287     {
1288         volScalarField scalarCellLevel
1289         (
1290             IOobject
1291             (
1292                 "cellLevel",
1293                 time().timeName(),
1294                 *this,
1295                 IOobject::NO_READ,
1296                 IOobject::AUTO_WRITE,
1297                 false
1298             ),
1299             *this,
1300             dimensionedScalar("level", dimless, 0)
1301         );
1302
1303         const labelList& cellLevel = meshCutter_.cellLevel();
1304
1305         forAll(cellLevel, cellI)
1306         {
1307             scalarCellLevel[cellI] = cellLevel[cellI];
1308         }
1309
1310         writeOk = writeOk && scalarCellLevel.write();
1311     }
1312
1313     return writeOk;
1314 }
1315
1316
1317 // ************************************************************************* //