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