src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C

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