src/engine/engineTopoChangerMesh/pistonRefine/pistonRefine.C

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