src/dynamicMesh/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C

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