src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.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 "meshRefinement.H"
  28 #include "volMesh.H"
  29 #include "volFields.H"
  30 #include "surfaceMesh.H"
  31 #include "syncTools.H"
  32 #include "Time.H"
  33 #include "refinementHistory.H"
  34 #include "refinementSurfaces.H"
  35 #include "decompositionMethod.H"
  36 #include "regionSplit.H"
  37 #include "fvMeshDistribute.H"
  38 #include "indirectPrimitivePatch.H"
  39 #include "polyTopoChange.H"
  40 #include "removeCells.H"
  41 #include "mapDistributePolyMesh.H"
  42 #include "localPointRegion.H"
  43 #include "pointMesh.H"
  44 #include "pointFields.H"
  45 #include "slipPointPatchFields.H"
  46 #include "fixedValuePointPatchFields.H"
  47 #include "globalPointPatchFields.H"
  48 #include "calculatedPointPatchFields.H"
  49 #include "processorPointPatch.H"
  50 #include "globalIndex.H"
  51 #include "meshTools.H"
  52 #include "OFstream.H"
  53
  54 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  55
  56 namespace Foam
  57 {
  58     defineTypeNameAndDebug(meshRefinement, 0);
  59 }
  60
  61 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  62
  63 void Foam::meshRefinement::calcNeighbourData
  64 (
  65     labelList& neiLevel,
  66     pointField& neiCc
  67 )  const
  68 {
  69     const labelList& cellLevel = meshCutter_.cellLevel();
  70     const pointField& cellCentres = mesh_.cellCentres();
  71
  72     label nBoundaryFaces = mesh_.nFaces() - mesh_.nInternalFaces();
  73
  74     if (neiLevel.size() != nBoundaryFaces || neiCc.size() != nBoundaryFaces)
  75     {
  76         FatalErrorIn("meshRefinement::calcNeighbour(..)") << "nBoundaries:"
  77             << nBoundaryFaces << " neiLevel:" << neiLevel.size()
  78             << abort(FatalError);
  79     }
  80
  81     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
  82
  83     forAll(patches, patchI)
  84     {
  85         const polyPatch& pp = patches[patchI];
  86
  87         const unallocLabelList& faceCells = pp.faceCells();
  88         const vectorField::subField faceCentres = pp.faceCentres();
  89
  90         label bFaceI = pp.start()-mesh_.nInternalFaces();
  91
  92         if (pp.coupled())
  93         {
  94             forAll(faceCells, i)
  95             {
  96                 neiLevel[bFaceI] = cellLevel[faceCells[i]];
  97                 neiCc[bFaceI] = cellCentres[faceCells[i]];
  98                 bFaceI++;
  99             }
 100         }
 101         else
 102         {
 103             forAll(faceCells, i)
 104             {
 105                 neiLevel[bFaceI] = cellLevel[faceCells[i]];
 106                 neiCc[bFaceI] = faceCentres[i];
 107                 bFaceI++;
 108             }
 109         }
 110     }
 111
 112     // Swap coupled boundaries. Apply separation to cc since is coordinate.
 113     syncTools::swapBoundaryFaceList(mesh_, neiCc, true);
 114     syncTools::swapBoundaryFaceList(mesh_, neiLevel, false);
 115 }
 116
 117
 118 // Find intersections of edges (given by their two endpoints) with surfaces.
 119 // Returns first intersection if there are more than one.
 120 void Foam::meshRefinement::updateIntersections(const labelList& changedFaces)
 121 {
 122     const pointField& cellCentres = mesh_.cellCentres();
 123
 124     label nTotEdges = returnReduce(surfaceIndex_.size(), sumOp<label>());
 125     label nChangedFaces = returnReduce(changedFaces.size(), sumOp<label>());
 126
 127     Info<< "Edge intersection testing:" << nl
 128         << "    Number of edges             : " << nTotEdges << nl
 129         << "    Number of edges to retest   : " << nChangedFaces
 130         << endl;
 131
 132     // Get boundary face centre and level. Coupled aware.
 133     labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
 134     pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
 135     calcNeighbourData(neiLevel, neiCc);
 136
 137     // Collect segments we want to test for
 138     pointField start(changedFaces.size());
 139     pointField end(changedFaces.size());
 140
 141     forAll(changedFaces, i)
 142     {
 143         label faceI = changedFaces[i];
 144         label own = mesh_.faceOwner()[faceI];
 145
 146         start[i] = cellCentres[own];
 147         if (mesh_.isInternalFace(faceI))
 148         {
 149             end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
 150         }
 151         else
 152         {
 153             end[i] = neiCc[faceI-mesh_.nInternalFaces()];
 154         }
 155     }
 156
 157     // Do tests in one go
 158     labelList surfaceHit;
 159     {
 160         labelList surfaceLevel;
 161         surfaces_.findHigherIntersection
 162         (
 163             start,
 164             end,
 165             labelList(start.size(), -1),    // accept any intersection
 166             surfaceHit,
 167             surfaceLevel
 168         );
 169     }
 170
 171     // Keep just surface hit
 172     forAll(surfaceHit, i)
 173     {
 174         surfaceIndex_[changedFaces[i]] = surfaceHit[i];
 175     }
 176
 177     // Make sure both sides have same information. This should be
 178     // case in general since same vectors but just to make sure.
 179     syncTools::syncFaceList(mesh_, surfaceIndex_, maxEqOp<label>(), false);
 180
 181     label nHits = countHits();
 182     label nTotHits = returnReduce(nHits, sumOp<label>());
 183
 184     Info<< "    Number of intersected edges : " << nTotHits << endl;
 185
 186     // Set files to same time as mesh
 187     setInstance(mesh_.facesInstance());
 188 }
 189
 190
 191 void Foam::meshRefinement::checkData()
 192 {
 193     Pout<< "meshRefinement::checkData() : Checking refinement structure."
 194         << endl;
 195     meshCutter_.checkMesh();
 196
 197     Pout<< "meshRefinement::checkData() : Checking refinement levels."
 198         << endl;
 199     meshCutter_.checkRefinementLevels(1, labelList(0));
 200
 201
 202     Pout<< "meshRefinement::checkData() : Checking synchronization."
 203         << endl;
 204
 205     // Check face centres
 206     {
 207         // Boundary face centres
 208         pointField::subList boundaryFc
 209         (
 210             mesh_.faceCentres(),
 211             mesh_.nFaces()-mesh_.nInternalFaces(),
 212             mesh_.nInternalFaces()
 213         );
 214
 215         // Get neighbouring face centres
 216         pointField neiBoundaryFc(boundaryFc);
 217         syncTools::swapBoundaryFaceList
 218         (
 219             mesh_,
 220             neiBoundaryFc,
 221             true
 222         );
 223
 224         // Compare
 225         testSyncBoundaryFaceList
 226         (
 227             mergeDistance_,
 228             "testing faceCentres : ",
 229             boundaryFc,
 230             neiBoundaryFc
 231         );
 232     }
 233     // Check meshRefinement
 234     {
 235         // Get boundary face centre and level. Coupled aware.
 236         labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
 237         pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
 238         calcNeighbourData(neiLevel, neiCc);
 239
 240         // Collect segments we want to test for
 241         pointField start(mesh_.nFaces());
 242         pointField end(mesh_.nFaces());
 243
 244         forAll(start, faceI)
 245         {
 246             start[faceI] = mesh_.cellCentres()[mesh_.faceOwner()[faceI]];
 247
 248             if (mesh_.isInternalFace(faceI))
 249             {
 250                 end[faceI] = mesh_.cellCentres()[mesh_.faceNeighbour()[faceI]];
 251             }
 252             else
 253             {
 254                 end[faceI] = neiCc[faceI-mesh_.nInternalFaces()];
 255             }
 256         }
 257
 258         // Do tests in one go
 259         labelList surfaceHit;
 260         {
 261             labelList surfaceLevel;
 262             surfaces_.findHigherIntersection
 263             (
 264                 start,
 265                 end,
 266                 labelList(start.size(), -1),    // accept any intersection
 267                 surfaceHit,
 268                 surfaceLevel
 269             );
 270         }
 271
 272         // Check
 273         forAll(surfaceHit, faceI)
 274         {
 275             if (surfaceHit[faceI] != surfaceIndex_[faceI])
 276             {
 277                 if (mesh_.isInternalFace(faceI))
 278                 {
 279                     WarningIn("meshRefinement::checkData()")
 280                         << "Internal face:" << faceI
 281                         << " fc:" << mesh_.faceCentres()[faceI]
 282                         << " cached surfaceIndex_:" << surfaceIndex_[faceI]
 283                         << " current:" << surfaceHit[faceI]
 284                         << " ownCc:"
 285                         << mesh_.cellCentres()[mesh_.faceOwner()[faceI]]
 286                         << " neiCc:"
 287                         << mesh_.cellCentres()[mesh_.faceNeighbour()[faceI]]
 288                         << endl;
 289                 }
 290                 else
 291                 {
 292                     WarningIn("meshRefinement::checkData()")
 293                         << "Boundary face:" << faceI
 294                         << " fc:" << mesh_.faceCentres()[faceI]
 295                         << " cached surfaceIndex_:" << surfaceIndex_[faceI]
 296                         << " current:" << surfaceHit[faceI]
 297                         << " ownCc:"
 298                         << mesh_.cellCentres()[mesh_.faceOwner()[faceI]]
 299                         << endl;
 300                 }
 301             }
 302         }
 303     }
 304     {
 305         labelList::subList boundarySurface
 306         (
 307             surfaceIndex_,
 308             mesh_.nFaces()-mesh_.nInternalFaces(),
 309             mesh_.nInternalFaces()
 310         );
 311
 312         labelList neiBoundarySurface(boundarySurface);
 313         syncTools::swapBoundaryFaceList
 314         (
 315             mesh_,
 316             neiBoundarySurface,
 317             false
 318         );
 319
 320         // Compare
 321         testSyncBoundaryFaceList
 322         (
 323             0,                              // tolerance
 324             "testing surfaceIndex() : ",
 325             boundarySurface,
 326             neiBoundarySurface
 327         );
 328     }
 329
 330
 331     // Find duplicate faces
 332     Pout<< "meshRefinement::checkData() : Counting duplicate faces."
 333         << endl;
 334
 335     labelList duplicateFace
 336     (
 337         localPointRegion::findDuplicateFaces
 338         (
 339             mesh_,
 340             identity(mesh_.nFaces()-mesh_.nInternalFaces())
 341           + mesh_.nInternalFaces()
 342         )
 343     );
 344
 345     // Count
 346     {
 347         label nDup = 0;
 348
 349         forAll(duplicateFace, i)
 350         {
 351             if (duplicateFace[i] != -1)
 352             {
 353                 nDup++;
 354             }
 355         }
 356         nDup /= 2;  // will have counted both faces of duplicate
 357         Pout<< "meshRefinement::checkData() : Found " << nDup
 358             << " duplicate pairs of faces." << endl;
 359     }
 360 }
 361
 362
 363 void Foam::meshRefinement::setInstance(const fileName& inst)
 364 {
 365     meshCutter_.setInstance(inst);
 366     surfaceIndex_.instance() = inst;
 367 }
 368
 369
 370 // Remove cells. Put exposedFaces (output of getExposedFaces(cellsToRemove))
 371 // into exposedPatchIDs.
 372 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRemoveCells
 373 (
 374     const labelList& cellsToRemove,
 375     const labelList& exposedFaces,
 376     const labelList& exposedPatchIDs,
 377     removeCells& cellRemover
 378 )
 379 {
 380     polyTopoChange meshMod(mesh_);
 381
 382     // Arbitrary: put exposed faces into last patch.
 383     cellRemover.setRefinement
 384     (
 385         cellsToRemove,
 386         exposedFaces,
 387         exposedPatchIDs,
 388         meshMod
 389     );
 390
 391     // Change the mesh (no inflation)
 392     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
 393
 394     // Update fields
 395     mesh_.updateMesh(map);
 396
 397     // Move mesh (since morphing might not do this)
 398     if (map().hasMotionPoints())
 399     {
 400         mesh_.movePoints(map().preMotionPoints());
 401     }
 402     else
 403     {
 404         // Delete mesh volumes. No other way to do this?
 405         mesh_.clearOut();
 406     }
 407
 408     // Update local mesh data
 409     cellRemover.updateMesh(map);
 410
 411     // Update intersections. Recalculate intersections for exposed faces.
 412     labelList newExposedFaces = renumber
 413     (
 414         map().reverseFaceMap(),
 415         exposedFaces
 416     );
 417
 418     //Pout<< "removeCells : updating intersections for "
 419     //    << newExposedFaces.size() << " newly exposed faces." << endl;
 420
 421     updateMesh(map, newExposedFaces);
 422
 423     return map;
 424 }
 425
 426
 427 // Determine for multi-processor regions the lowest numbered cell on the lowest
 428 // numbered processor.
 429 void Foam::meshRefinement::getRegionMaster
 430 (
 431     const boolList& blockedFace,
 432     const regionSplit& globalRegion,
 433     Map<label>& regionToMaster
 434 ) const
 435 {
 436     globalIndex globalCells(mesh_.nCells());
 437
 438     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 439
 440     forAll(patches, patchI)
 441     {
 442         const polyPatch& pp = patches[patchI];
 443
 444         if (isA<processorPolyPatch>(pp))
 445         {
 446             forAll(pp, i)
 447             {
 448                 label faceI = pp.start()+i;
 449
 450                 if (!blockedFace[faceI])
 451                 {
 452                     // Only if there is a connection to the neighbour
 453                     // will there be a multi-domain region; if not through
 454                     // this face then through another.
 455
 456                     label cellI = mesh_.faceOwner()[faceI];
 457                     label globalCellI = globalCells.toGlobal(cellI);
 458
 459                     Map<label>::iterator iter =
 460                         regionToMaster.find(globalRegion[cellI]);
 461
 462                     if (iter != regionToMaster.end())
 463                     {
 464                         label master = iter();
 465                         iter() = min(master, globalCellI);
 466                     }
 467                     else
 468                     {
 469                         regionToMaster.insert
 470                         (
 471                             globalRegion[cellI],
 472                             globalCellI
 473                         );
 474                     }
 475                 }
 476             }
 477         }
 478     }
 479
 480     // Do reduction
 481     Pstream::mapCombineGather(regionToMaster, minEqOp<label>());
 482     Pstream::mapCombineScatter(regionToMaster);
 483 }
 484
 485
 486 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 487
 488 // Construct from components
 489 Foam::meshRefinement::meshRefinement
 490 (
 491     fvMesh& mesh,
 492     const scalar mergeDistance,
 493     const refinementSurfaces& surfaces,
 494     const shellSurfaces& shells
 495 )
 496 :
 497     mesh_(mesh),
 498     mergeDistance_(mergeDistance),
 499     surfaces_(surfaces),
 500     shells_(shells),
 501     meshCutter_
 502     (
 503         mesh,
 504         labelIOList
 505         (
 506             IOobject
 507             (
 508                 "cellLevel",
 509                 mesh_.facesInstance(),
 510                 fvMesh::meshSubDir,
 511                 mesh,
 512                 IOobject::READ_IF_PRESENT,
 513                 IOobject::NO_WRITE,
 514                 false
 515             ),
 516             labelList(mesh_.nCells(), 0)
 517         ),
 518         labelIOList
 519         (
 520             IOobject
 521             (
 522                 "pointLevel",
 523                 mesh_.facesInstance(),
 524                 fvMesh::meshSubDir,
 525                 mesh,
 526                 IOobject::READ_IF_PRESENT,
 527                 IOobject::NO_WRITE,
 528                 false
 529             ),
 530             labelList(mesh_.nPoints(), 0)
 531         ),
 532         refinementHistory
 533         (
 534             IOobject
 535             (
 536                 "refinementHistory",
 537                 mesh_.facesInstance(),
 538                 fvMesh::meshSubDir,
 539                 mesh_,
 540                 IOobject::NO_READ,
 541                 IOobject::NO_WRITE,
 542                 false
 543             ),
 544             List<refinementHistory::splitCell8>(0),
 545             labelList(0)
 546         )   // no unrefinement
 547     ),
 548     surfaceIndex_
 549     (
 550         IOobject
 551         (
 552             "surfaceIndex",
 553             mesh_.facesInstance(),
 554             fvMesh::meshSubDir,
 555             mesh_,
 556             IOobject::NO_READ,
 557             IOobject::NO_WRITE,
 558             false
 559         ),
 560         labelList(mesh_.nFaces(), -1)
 561     ),
 562     userFaceData_(0)
 563 {
 564     // recalculate intersections for all faces
 565     updateIntersections(identity(mesh_.nFaces()));
 566 }
 567
 568
 569 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 570
 571 Foam::label Foam::meshRefinement::countHits() const
 572 {
 573     label nHits = 0;
 574
 575     forAll(surfaceIndex_, faceI)
 576     {
 577         if (surfaceIndex_[faceI] >= 0)
 578         {
 579             nHits++;
 580         }
 581     }
 582     return nHits;
 583 }
 584
 585
 586 // Determine distribution to move connected regions onto one processor.
 587 Foam::labelList Foam::meshRefinement::decomposeCombineRegions
 588 (
 589     const boolList& blockedFace,
 590     const List<labelPair>& explicitConnections,
 591     decompositionMethod& decomposer
 592 ) const
 593 {
 594     // Determine global regions, separated by blockedFaces
 595     regionSplit globalRegion(mesh_, blockedFace, explicitConnections);
 596
 597     // Now globalRegion has global region per cell. Problem is that
 598     // the region might span multiple domains so we want to get
 599     // a "region master" per domain. Note that multi-processor
 600     // regions can only occur on cells on coupled patches.
 601
 602
 603     // Determine per coupled region the master cell (lowest numbered cell
 604     // on lowest numbered processor)
 605     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 606
 607     Map<label> regionToMaster(mesh_.nFaces()-mesh_.nInternalFaces());
 608     getRegionMaster(blockedFace, globalRegion, regionToMaster);
 609
 610
 611     // Global cell numbering engine
 612     globalIndex globalCells(mesh_.nCells());
 613
 614
 615     // Determine cell centre per region
 616     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 617     // Now we can divide regions into
 618     // - single-processor regions (almost all)
 619     //   - allocate local region + coordinate for it
 620     // - multi-processor for which I am master
 621     //   - allocate local region + coordinate for it
 622     // - multi-processor for which I am not master
 623     //   - do not allocate region.
 624     //   - but inherit new distribution from master.
 625
 626     Map<label> globalToLocalRegion(mesh_.nCells());
 627     DynamicList<point> localCc(mesh_.nCells()/10);
 628
 629     forAll(globalRegion, cellI)
 630     {
 631         Map<label>::const_iterator fndMaster =
 632             regionToMaster.find(globalRegion[cellI]);
 633
 634         if (fndMaster != regionToMaster.end())
 635         {
 636             // Multi-processor region.
 637             if (globalCells.toGlobal(cellI) == fndMaster())
 638             {
 639                 // I am master. Allocate region for me.
 640                 globalToLocalRegion.insert(globalRegion[cellI], localCc.size());
 641                 localCc.append(mesh_.cellCentres()[cellI]);
 642             }
 643         }
 644         else
 645         {
 646             // Single processor region.
 647             Map<label>::iterator iter =
 648                 globalToLocalRegion.find(globalRegion[cellI]);
 649
 650             if (iter == globalToLocalRegion.end())
 651             {
 652                 globalToLocalRegion.insert(globalRegion[cellI], localCc.size());
 653                 localCc.append(mesh_.cellCentres()[cellI]);
 654             }
 655         }
 656     }
 657     localCc.shrink();
 658     pointField localPoints;
 659     localPoints.transfer(localCc);
 660
 661
 662     // Call decomposer on localCc
 663     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 664
 665     labelList localDistribution = decomposer.decompose(localPoints);
 666
 667
 668     // Distribute the destination processor for coupled regions
 669     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 670
 671     Map<label> regionToDist(regionToMaster.size());
 672
 673     forAllConstIter(Map<label>, regionToMaster, iter)
 674     {
 675         if (globalCells.isLocal(iter()))
 676         {
 677             // Master cell is local.
 678             regionToDist.insert
 679             (
 680                 iter.key(),
 681                 localDistribution[globalToLocalRegion[iter.key()]]
 682             );
 683         }
 684         else
 685         {
 686             // Master cell is not on this processor. Make sure it is overridden.
 687             regionToDist.insert(iter.key(), labelMax);
 688         }
 689     }
 690     Pstream::mapCombineGather(regionToDist, minEqOp<label>());
 691     Pstream::mapCombineScatter(regionToDist);
 692
 693
 694     // Convert region-based decomposition back to cell-based one
 695     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 696
 697     labelList distribution(mesh_.nCells());
 698
 699     forAll(globalRegion, cellI)
 700     {
 701         Map<label>::const_iterator fndMaster =
 702             regionToDist.find(globalRegion[cellI]);
 703
 704         if (fndMaster != regionToDist.end())
 705         {
 706             // Special handling
 707             distribution[cellI] = fndMaster();
 708         }
 709         else
 710         {
 711             // region is local to the processor.
 712             label localRegionI = globalToLocalRegion[globalRegion[cellI]];
 713
 714             distribution[cellI] = localDistribution[localRegionI];
 715         }
 716     }
 717     return distribution;
 718 }
 719
 720
 721 Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
 722 (
 723     const bool keepZoneFaces,
 724     const bool keepBaffles,
 725     decompositionMethod& decomposer,
 726     fvMeshDistribute& distributor
 727 )
 728 {
 729     autoPtr<mapDistributePolyMesh> map;
 730
 731     if (Pstream::parRun())
 732     {
 733         //Info<< nl
 734         //    << "Doing final balancing" << nl
 735         //    << "---------------------" << nl
 736         //    << endl;
 737         //
 738         //if (debug_)
 739         //{
 740         //    const_cast<Time&>(mesh_.time())++;
 741         //}
 742
 743         // Wanted distribution
 744         labelList distribution;
 745
 746         if (keepZoneFaces || keepBaffles)
 747         {
 748             // Faces where owner and neighbour are not 'connected' so can
 749             // go to different processors.
 750             boolList blockedFace(mesh_.nFaces(), true);
 751             // Pairs of baffles
 752             List<labelPair> couples;
 753
 754             if (keepZoneFaces)
 755             {
 756                 label nNamed = surfaces().getNamedSurfaces().size();
 757
 758                 Info<< "Found " << nNamed << " surfaces with faceZones."
 759                     << " Applying special decomposition to keep those together."
 760                     << endl;
 761
 762                 // Determine decomposition to keep/move surface zones
 763                 // on one processor. The reason is that snapping will make these
 764                 // into baffles, move and convert them back so if they were
 765                 // proc boundaries after baffling&moving the points might be no
 766                 // longer snychronised so recoupling will fail. To prevent this
 767                 // keep owner&neighbour of such a surface zone on the same
 768                 // processor.
 769
 770                 const wordList& fzNames = surfaces().faceZoneNames();
 771                 const faceZoneMesh& fZones = mesh_.faceZones();
 772
 773                 // Get faces whose owner and neighbour should stay together,
 774                 // i.e. they are not 'blocked'.
 775
 776                 label nZoned = 0;
 777
 778                 forAll(fzNames, surfI)
 779                 {
 780                     if (fzNames[surfI].size() > 0)
 781                     {
 782                         // Get zone
 783                         label zoneI = fZones.findZoneID(fzNames[surfI]);
 784
 785                         const faceZone& fZone = fZones[zoneI];
 786
 787                         forAll(fZone, i)
 788                         {
 789                             if (blockedFace[fZone[i]])
 790                             {
 791                                 blockedFace[fZone[i]] = false;
 792                                 nZoned++;
 793                             }
 794                         }
 795                     }
 796                 }
 797                 Info<< "Found " << returnReduce(nZoned, sumOp<label>())
 798                     << " zoned faces to keep together."
 799                     << endl;
 800             }
 801
 802             if (keepBaffles)
 803             {
 804                 // Get boundary baffles that need to stay together.
 805                 couples = getDuplicateFaces   // all baffles
 806                 (
 807                     identity(mesh_.nFaces()-mesh_.nInternalFaces())
 808                    +mesh_.nInternalFaces()
 809                 );
 810
 811                 Info<< "Found " << returnReduce(couples.size(), sumOp<label>())
 812                     << " baffles to keep together."
 813                     << endl;
 814             }
 815
 816             distribution = decomposeCombineRegions
 817             (
 818                 blockedFace,
 819                 couples,
 820                 decomposer
 821             );
 822
 823             labelList nProcCells(distributor.countCells(distribution));
 824             Pstream::listCombineGather(nProcCells, plusEqOp<label>());
 825             Pstream::listCombineScatter(nProcCells);
 826
 827             Info<< "Calculated decomposition:" << endl;
 828             forAll(nProcCells, procI)
 829             {
 830                 Info<< "    " << procI << '\t' << nProcCells[procI] << endl;
 831             }
 832             Info<< endl;
 833         }
 834         else
 835         {
 836             // Normal decomposition
 837             distribution = decomposer.decompose(mesh_.cellCentres());
 838         }
 839
 840         if (debug)
 841         {
 842             Pout<< "Wanted distribution:"
 843                 << distributor.countCells(distribution)
 844                 << endl;
 845         }
 846         // Do actual sending/receiving of mesh
 847         map = distributor.distribute(distribution);
 848
 849         // Update numbering of meshRefiner
 850         distribute(map);
 851     }
 852     return map;
 853 }
 854
 855
 856 // Helper function to get intersected faces
 857 Foam::labelList Foam::meshRefinement::intersectedFaces() const
 858 {
 859     // Mark all faces that will become baffles
 860
 861     label nBoundaryFaces = 0;
 862
 863     forAll(surfaceIndex_, faceI)
 864     {
 865         if (surfaceIndex_[faceI] != -1)
 866         {
 867             nBoundaryFaces++;
 868         }
 869     }
 870
 871     labelList surfaceFaces(nBoundaryFaces);
 872     nBoundaryFaces = 0;
 873
 874     forAll(surfaceIndex_, faceI)
 875     {
 876         if (surfaceIndex_[faceI] != -1)
 877         {
 878             surfaceFaces[nBoundaryFaces++] = faceI;
 879         }
 880     }
 881     return surfaceFaces;
 882 }
 883
 884
 885 // Helper function to get points used by faces
 886 Foam::labelList Foam::meshRefinement::intersectedPoints
 887 (
 888 //    const labelList& globalToPatch
 889 ) const
 890 {
 891     const faceList& faces = mesh_.faces();
 892
 893     // Mark all points on faces that will become baffles
 894     PackedList<1> isBoundaryPoint(mesh_.nPoints(), 0u);
 895     label nBoundaryPoints = 0;
 896
 897     forAll(surfaceIndex_, faceI)
 898     {
 899         if (surfaceIndex_[faceI] != -1)
 900         {
 901             const face& f = faces[faceI];
 902
 903             forAll(f, fp)
 904             {
 905                 if (isBoundaryPoint.set(f[fp], 1u))
 906                 {
 907                     nBoundaryPoints++;
 908                 }
 909             }
 910         }
 911     }
 912
 913     //// Insert all meshed patches.
 914     //forAll(globalToPatch, i)
 915     //{
 916     //    label patchI = globalToPatch[i];
 917     //
 918     //    if (patchI != -1)
 919     //    {
 920     //        const polyPatch& pp = mesh_.boundaryMesh()[patchI];
 921     //
 922     //        label faceI = pp.start();
 923     //
 924     //        forAll(pp, i)
 925     //        {
 926     //            const face& f = faces[faceI];
 927     //
 928     //            forAll(f, fp)
 929     //            {
 930     //                if (isBoundaryPoint.set(f[fp], 1u))
 931     //                    nBoundaryPoints++;
 932     //                }
 933     //            }
 934     //            faceI++;
 935     //        }
 936     //    }
 937     //}
 938
 939
 940     // Pack
 941     labelList boundaryPoints(nBoundaryPoints);
 942     nBoundaryPoints = 0;
 943     forAll(isBoundaryPoint, pointI)
 944     {
 945         if (isBoundaryPoint.get(pointI) == 1u)
 946         {
 947             boundaryPoints[nBoundaryPoints++] = pointI;
 948         }
 949     }
 950
 951     return boundaryPoints;
 952 }
 953
 954
 955 Foam::labelList Foam::meshRefinement::addedPatches
 956 (
 957     const labelList& globalToPatch
 958 )
 959 {
 960     labelList patchIDs(globalToPatch.size());
 961     label addedI = 0;
 962
 963     forAll(globalToPatch, i)
 964     {
 965         if (globalToPatch[i] != -1)
 966         {
 967             patchIDs[addedI++] = globalToPatch[i];
 968         }
 969     }
 970     patchIDs.setSize(addedI);
 971
 972     return patchIDs;
 973 }
 974
 975
 976 //- Create patch from set of patches
 977 Foam::autoPtr<Foam::indirectPrimitivePatch> Foam::meshRefinement::makePatch
 978 (
 979     const polyMesh& mesh,
 980     const labelList& patchIDs
 981 )
 982 {
 983     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 984
 985     // Count faces.
 986     label nFaces = 0;
 987
 988     forAll(patchIDs, i)
 989     {
 990         const polyPatch& pp = patches[patchIDs[i]];
 991
 992         nFaces += pp.size();
 993     }
 994
 995     // Collect faces.
 996     labelList addressing(nFaces);
 997     nFaces = 0;
 998
 999     forAll(patchIDs, i)
1000     {
1001         const polyPatch& pp = patches[patchIDs[i]];
1002
1003         label meshFaceI = pp.start();
1004
1005         forAll(pp, i)
1006         {
1007             addressing[nFaces++] = meshFaceI++;
1008         }
1009     }
1010
1011     return autoPtr<indirectPrimitivePatch>
1012     (
1013         new indirectPrimitivePatch
1014         (
1015             IndirectList<face>(mesh.faces(), addressing),
1016             mesh.points()
1017         )
1018     );
1019 }
1020
1021
1022 // Construct pointVectorField with correct boundary conditions
1023 Foam::tmp<Foam::pointVectorField> Foam::meshRefinement::makeDisplacementField
1024 (
1025     const pointMesh& pMesh,
1026     const labelList& adaptPatchIDs
1027 )
1028 {
1029     const polyMesh& mesh = pMesh();
1030
1031     // Construct displacement field.
1032     const pointBoundaryMesh& pointPatches = pMesh.boundary();
1033
1034     wordList patchFieldTypes
1035     (
1036         pointPatches.size(),
1037         slipPointPatchVectorField::typeName
1038     );
1039
1040     forAll(adaptPatchIDs, i)
1041     {
1042         patchFieldTypes[adaptPatchIDs[i]] =
1043             fixedValuePointPatchVectorField::typeName;
1044     }
1045
1046     forAll(pointPatches, patchI)
1047     {
1048         if (isA<globalPointPatch>(pointPatches[patchI]))
1049         {
1050             patchFieldTypes[patchI] = globalPointPatchVectorField::typeName;
1051         }
1052         else if (isA<processorPointPatch>(pointPatches[patchI]))
1053         {
1054             patchFieldTypes[patchI] = calculatedPointPatchVectorField::typeName;
1055         }
1056     }
1057
1058     tmp<pointVectorField> tfld
1059     (
1060         new pointVectorField
1061         (
1062             IOobject
1063             (
1064                 "pointDisplacement",
1065                 mesh.time().timeName(),
1066                 mesh,
1067                 IOobject::NO_READ,
1068                 IOobject::AUTO_WRITE
1069             ),
1070             pMesh,
1071             dimensionedVector("displacement", dimLength, vector::zero),
1072             patchFieldTypes
1073         )
1074     );
1075     return tfld;
1076 }
1077
1078
1079 void Foam::meshRefinement::checkCoupledFaceZones(const polyMesh& mesh)
1080 {
1081     const faceZoneMesh& fZones = mesh.faceZones();
1082
1083     // Check any zones are present anywhere and in same order
1084
1085     {
1086         List<wordList> zoneNames(Pstream::nProcs());
1087         zoneNames[Pstream::myProcNo()] = fZones.names();
1088         Pstream::gatherList(zoneNames);
1089         Pstream::scatterList(zoneNames);
1090         // All have same data now. Check.
1091         forAll(zoneNames, procI)
1092         {
1093             if (procI != Pstream::myProcNo())
1094             {
1095                 if (zoneNames[procI] != zoneNames[Pstream::myProcNo()])
1096                 {
1097                     FatalErrorIn
1098                     (
1099                         "meshRefinement::checkCoupledFaceZones(const polyMesh&)"
1100                     )   << "faceZones are not synchronised on processors." << nl
1101                         << "Processor " << procI << " has faceZones "
1102                         << zoneNames[procI] << nl
1103                         << "Processor " << Pstream::myProcNo()
1104                         << " has faceZones "
1105                         << zoneNames[Pstream::myProcNo()] << nl
1106                         << exit(FatalError);
1107                 }
1108             }
1109         }
1110     }
1111
1112     // Check that coupled faces are present on both sides.
1113
1114     labelList faceToZone(mesh.nFaces()-mesh.nInternalFaces(), -1);
1115
1116     forAll(fZones, zoneI)
1117     {
1118         const faceZone& fZone = fZones[zoneI];
1119
1120         forAll(fZone, i)
1121         {
1122             label bFaceI = fZone[i]-mesh.nInternalFaces();
1123
1124             if (bFaceI >= 0)
1125             {
1126                 if (faceToZone[bFaceI] == -1)
1127                 {
1128                     faceToZone[bFaceI] = zoneI;
1129                 }
1130                 else if (faceToZone[bFaceI] == zoneI)
1131                 {
1132                     FatalErrorIn
1133                     (
1134                         "meshRefinement::checkCoupledFaceZones(const polyMesh&)"
1135                     )   << "Face " << fZone[i] << " in zone "
1136                         << fZone.name()
1137                         << " is twice in zone!"
1138                         << abort(FatalError);
1139                 }
1140                 else
1141                 {
1142                     FatalErrorIn
1143                     (
1144                         "meshRefinement::checkCoupledFaceZones(const polyMesh&)"
1145                     )   << "Face " << fZone[i] << " in zone "
1146                         << fZone.name()
1147                         << " is also in zone "
1148                         << fZones[faceToZone[bFaceI]].name()
1149                         << abort(FatalError);
1150                 }
1151             }
1152         }
1153     }
1154
1155     labelList neiFaceToZone(faceToZone);
1156     syncTools::swapBoundaryFaceList(mesh, neiFaceToZone, false);
1157
1158     forAll(faceToZone, i)
1159     {
1160         if (faceToZone[i] != neiFaceToZone[i])
1161         {
1162             FatalErrorIn
1163             (
1164                 "meshRefinement::checkCoupledFaceZones(const polyMesh&)"
1165             )   << "Face " << mesh.nInternalFaces()+i
1166                 << " is in zone " << faceToZone[i]
1167                 << ", its coupled face is in zone " << neiFaceToZone[i]
1168                 << abort(FatalError);
1169         }
1170     }
1171 }
1172
1173
1174 // Adds patch if not yet there. Returns patchID.
1175 Foam::label Foam::meshRefinement::addPatch
1176 (
1177     fvMesh& mesh,
1178     const word& patchName,
1179     const word& patchType
1180 )
1181 {
1182     polyBoundaryMesh& polyPatches =
1183         const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
1184
1185     label patchI = polyPatches.findPatchID(patchName);
1186     if (patchI != -1)
1187     {
1188         if (polyPatches[patchI].type() == patchType)
1189         {
1190             // Already there
1191             return patchI;
1192         }
1193         //else
1194         //{
1195         //    FatalErrorIn
1196         //    (
1197         //        "meshRefinement::addPatch(fvMesh&, const word&, const word&)"
1198         //    )   << "Patch " << patchName << " already exists but with type "
1199         //        << patchType << nl
1200         //        << "Current patch names:" << polyPatches.names()
1201         //        << exit(FatalError);
1202         //}
1203     }
1204
1205
1206     label insertPatchI = polyPatches.size();
1207     label startFaceI = mesh.nFaces();
1208
1209     forAll(polyPatches, patchI)
1210     {
1211         const polyPatch& pp = polyPatches[patchI];
1212
1213         if (isA<processorPolyPatch>(pp))
1214         {
1215             insertPatchI = patchI;
1216             startFaceI = pp.start();
1217             break;
1218         }
1219     }
1220
1221
1222     // Below is all quite a hack. Feel free to change once there is a better
1223     // mechanism to insert and reorder patches.
1224
1225     // Clear local fields and e.g. polyMesh parallelInfo.
1226     mesh.clearOut();
1227
1228     label sz = polyPatches.size();
1229
1230     fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
1231
1232     // Add polyPatch at the end
1233     polyPatches.setSize(sz+1);
1234     polyPatches.set
1235     (
1236         sz,
1237         polyPatch::New
1238         (
1239             patchType,
1240             patchName,
1241             0,              // size
1242             startFaceI,
1243             insertPatchI,
1244             polyPatches
1245         )
1246     );
1247     fvPatches.setSize(sz+1);
1248     fvPatches.set
1249     (
1250         sz,
1251         fvPatch::New
1252         (
1253             polyPatches[sz],  // point to newly added polyPatch
1254             mesh.boundary()
1255         )
1256     );
1257
1258     addPatchFields<volScalarField>
1259     (
1260         mesh,
1261         calculatedFvPatchField<scalar>::typeName
1262     );
1263     addPatchFields<volVectorField>
1264     (
1265         mesh,
1266         calculatedFvPatchField<vector>::typeName
1267     );
1268     addPatchFields<volSphericalTensorField>
1269     (
1270         mesh,
1271         calculatedFvPatchField<sphericalTensor>::typeName
1272     );
1273     addPatchFields<volSymmTensorField>
1274     (
1275         mesh,
1276         calculatedFvPatchField<symmTensor>::typeName
1277     );
1278     addPatchFields<volTensorField>
1279     (
1280         mesh,
1281         calculatedFvPatchField<tensor>::typeName
1282     );
1283
1284     // Surface fields
1285
1286     addPatchFields<surfaceScalarField>
1287     (
1288         mesh,
1289         calculatedFvPatchField<scalar>::typeName
1290     );
1291     addPatchFields<surfaceVectorField>
1292     (
1293         mesh,
1294         calculatedFvPatchField<vector>::typeName
1295     );
1296     addPatchFields<surfaceSphericalTensorField>
1297     (
1298         mesh,
1299         calculatedFvPatchField<sphericalTensor>::typeName
1300     );
1301     addPatchFields<surfaceSymmTensorField>
1302     (
1303         mesh,
1304         calculatedFvPatchField<symmTensor>::typeName
1305     );
1306     addPatchFields<surfaceTensorField>
1307     (
1308         mesh,
1309         calculatedFvPatchField<tensor>::typeName
1310     );
1311
1312     // Create reordering list
1313     // patches before insert position stay as is
1314     labelList oldToNew(sz+1);
1315     for (label i = 0; i < insertPatchI; i++)
1316     {
1317         oldToNew[i] = i;
1318     }
1319     // patches after insert position move one up
1320     for (label i = insertPatchI; i < sz; i++)
1321     {
1322         oldToNew[i] = i+1;
1323     }
1324     // appended patch gets moved to insert position
1325     oldToNew[sz] = insertPatchI;
1326
1327     // Shuffle into place
1328     polyPatches.reorder(oldToNew);
1329     fvPatches.reorder(oldToNew);
1330
1331     reorderPatchFields<volScalarField>(mesh, oldToNew);
1332     reorderPatchFields<volVectorField>(mesh, oldToNew);
1333     reorderPatchFields<volSphericalTensorField>(mesh, oldToNew);
1334     reorderPatchFields<volSymmTensorField>(mesh, oldToNew);
1335     reorderPatchFields<volTensorField>(mesh, oldToNew);
1336     reorderPatchFields<surfaceScalarField>(mesh, oldToNew);
1337     reorderPatchFields<surfaceVectorField>(mesh, oldToNew);
1338     reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
1339     reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
1340     reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
1341
1342     return insertPatchI;
1343 }
1344
1345
1346 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMeshRegions
1347 (
1348     const point& keepPoint
1349 )
1350 {
1351     // Determine connected regions. regionSplit is the labelList with the
1352     // region per cell.
1353     regionSplit cellRegion(mesh_);
1354
1355     label regionI = -1;
1356
1357     label cellI = mesh_.findCell(keepPoint);
1358
1359     if (cellI != -1)
1360     {
1361         regionI = cellRegion[cellI];
1362     }
1363
1364     reduce(regionI, maxOp<label>());
1365
1366     if (regionI == -1)
1367     {
1368         FatalErrorIn
1369         (
1370             "meshRefinement::splitMeshRegions(const point&)"
1371         )   << "Point " << keepPoint
1372             << " is not inside the mesh." << nl
1373             << "Bounding box of the mesh:" << mesh_.globalData().bb()
1374             << exit(FatalError);
1375     }
1376
1377     // Subset
1378     // ~~~~~~
1379
1380     // Get cells to remove
1381     DynamicList<label> cellsToRemove(mesh_.nCells());
1382     forAll(cellRegion, cellI)
1383     {
1384         if (cellRegion[cellI] != regionI)
1385         {
1386             cellsToRemove.append(cellI);
1387         }
1388     }
1389     cellsToRemove.shrink();
1390
1391     label nCellsToKeep = mesh_.nCells() - cellsToRemove.size();
1392     reduce(nCellsToKeep, sumOp<label>());
1393
1394     Info<< "Keeping all cells in region " << regionI
1395         << " containing point " << keepPoint << endl
1396         << "Selected for keeping : "
1397         << nCellsToKeep
1398         << " cells." << endl;
1399
1400
1401     // Remove cells
1402     removeCells cellRemover(mesh_);
1403
1404     labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
1405
1406     if (exposedFaces.size() > 0)
1407     {
1408         FatalErrorIn
1409         (
1410             "meshRefinement::splitMeshRegions(const point&)"
1411         )   << "Removing non-reachable cells should only expose boundary faces"
1412             << nl
1413             << "ExposedFaces:" << exposedFaces << abort(FatalError);
1414     }
1415
1416     return doRemoveCells
1417     (
1418         cellsToRemove,
1419         exposedFaces,
1420         labelList(exposedFaces.size(),-1),  // irrelevant since 0 size.
1421         cellRemover
1422     );
1423 }
1424
1425
1426 void Foam::meshRefinement::distribute(const mapDistributePolyMesh& map)
1427 {
1428     // mesh_ already distributed; distribute my member data
1429
1430     // surfaceQueries_ ok.
1431
1432     // refinement
1433     meshCutter_.distribute(map);
1434
1435     // surfaceIndex is face data.
1436     map.distributeFaceData(surfaceIndex_);
1437
1438     // maintainedFaces are indices of faces.
1439     forAll(userFaceData_, i)
1440     {
1441         map.distributeFaceData(userFaceData_[i].second());
1442     }
1443 }
1444
1445
1446 void Foam::meshRefinement::updateMesh
1447 (
1448     const mapPolyMesh& map,
1449     const labelList& changedFaces
1450 )
1451 {
1452     Map<label> dummyMap(0);
1453
1454     updateMesh(map, changedFaces, dummyMap, dummyMap, dummyMap);
1455 }
1456
1457
1458 void Foam::meshRefinement::storeData
1459 (
1460     const labelList& pointsToStore,
1461     const labelList& facesToStore,
1462     const labelList& cellsToStore
1463 )
1464 {
1465     // For now only meshCutter has storable/retrievable data.
1466     meshCutter_.storeData
1467     (
1468         pointsToStore,
1469         facesToStore,
1470         cellsToStore
1471     );
1472 }
1473
1474
1475 void Foam::meshRefinement::updateMesh
1476 (
1477     const mapPolyMesh& map,
1478     const labelList& changedFaces,
1479     const Map<label>& pointsToRestore,
1480     const Map<label>& facesToRestore,
1481     const Map<label>& cellsToRestore
1482 )
1483 {
1484     // For now only meshCutter has storable/retrievable data.
1485
1486     // Update numbering of cells/vertices.
1487     meshCutter_.updateMesh
1488     (
1489         map,
1490         pointsToRestore,
1491         facesToRestore,
1492         cellsToRestore
1493     );
1494
1495     // Update surfaceIndex
1496     updateList(map.faceMap(), -1, surfaceIndex_);
1497
1498     // Update cached intersection information
1499     updateIntersections(changedFaces);
1500
1501     // Update maintained faces
1502     forAll(userFaceData_, i)
1503     {
1504         labelList& data = userFaceData_[i].second();
1505
1506         if (userFaceData_[i].first() == KEEPALL)
1507         {
1508             // extend list with face-from-face data
1509             updateList(map.faceMap(), -1, data);
1510         }
1511         else if (userFaceData_[i].first() == MASTERONLY)
1512         {
1513             // keep master only
1514             labelList newFaceData(map.faceMap().size(), -1);
1515
1516             forAll(newFaceData, faceI)
1517             {
1518                 label oldFaceI = map.faceMap()[faceI];
1519
1520                 if (oldFaceI >= 0 && map.reverseFaceMap()[oldFaceI] == faceI)
1521                 {
1522                     newFaceData[faceI] = data[oldFaceI];
1523                 }
1524             }
1525             data.transfer(newFaceData);
1526         }
1527         else
1528         {
1529             // remove any face that has been refined i.e. referenced more than
1530             // once.
1531
1532             // 1. Determine all old faces that get referenced more than once.
1533             // These get marked with -1 in reverseFaceMap
1534             labelList reverseFaceMap(map.reverseFaceMap());
1535
1536             forAll(map.faceMap(), faceI)
1537             {
1538                 label oldFaceI = map.faceMap()[faceI];
1539
1540                 if (oldFaceI >= 0)
1541                 {
1542                     if (reverseFaceMap[oldFaceI] != faceI)
1543                     {
1544                         // faceI is slave face. Mark old face.
1545                         reverseFaceMap[oldFaceI] = -1;
1546                     }
1547                 }
1548             }
1549
1550             // 2. Map only faces with intact reverseFaceMap
1551             labelList newFaceData(map.faceMap().size(), -1);
1552             forAll(newFaceData, faceI)
1553             {
1554                 label oldFaceI = map.faceMap()[faceI];
1555
1556                 if (oldFaceI >= 0)
1557                 {
1558                     if (reverseFaceMap[oldFaceI] == faceI)
1559                     {
1560                         newFaceData[faceI] = data[oldFaceI];
1561                     }
1562                 }
1563             }
1564             data.transfer(newFaceData);
1565         }
1566     }
1567 }
1568
1569
1570 bool Foam::meshRefinement::write() const
1571 {
1572     bool writeOk =
1573         mesh_.write()
1574      && meshCutter_.write()
1575      && surfaceIndex_.write();
1576
1577     return writeOk;
1578 }
1579
1580
1581 void Foam::meshRefinement::printMeshInfo(const bool debug, const string& msg)
1582  const
1583 {
1584     const globalMeshData& pData = mesh_.globalData();
1585
1586     if (debug)
1587     {
1588         Pout<< msg.c_str()
1589             << " : cells(local):" << mesh_.nCells()
1590             << "  faces(local):" << mesh_.nFaces()
1591             << "  points(local):" << mesh_.nPoints()
1592             << endl;
1593     }
1594     else
1595     {
1596         Info<< msg.c_str()
1597             << " : cells:" << pData.nTotalCells()
1598             << "  faces:" << pData.nTotalFaces()
1599             << "  points:" << pData.nTotalPoints()
1600             << endl;
1601     }
1602
1603
1604     //if (debug)
1605     {
1606         const labelList& cellLevel = meshCutter_.cellLevel();
1607
1608         labelList nCells(gMax(cellLevel)+1, 0);
1609
1610         forAll(cellLevel, cellI)
1611         {
1612             nCells[cellLevel[cellI]]++;
1613         }
1614
1615         Pstream::listCombineGather(nCells, plusEqOp<label>());
1616         Pstream::listCombineScatter(nCells);
1617
1618         Info<< "Cells per refinement level:" << endl;
1619         forAll(nCells, levelI)
1620         {
1621             Info<< "    " << levelI << '\t' << nCells[levelI]
1622                 << endl;
1623         }
1624     }
1625 }
1626
1627
1628 void Foam::meshRefinement::dumpRefinementLevel() const
1629 {
1630     volScalarField volRefLevel
1631     (
1632         IOobject
1633         (
1634             "cellLevel",
1635             mesh_.time().timeName(),
1636             mesh_,
1637             IOobject::NO_READ,
1638             IOobject::AUTO_WRITE,
1639             false
1640         ),
1641         mesh_,
1642         dimensionedScalar("zero", dimless, 0)
1643     );
1644
1645     const labelList& cellLevel = meshCutter_.cellLevel();
1646
1647     forAll(volRefLevel, cellI)
1648     {
1649         volRefLevel[cellI] = cellLevel[cellI];
1650     }
1651
1652     volRefLevel.write();
1653
1654
1655     pointMesh pMesh(mesh_);
1656
1657     pointScalarField pointRefLevel
1658     (
1659         IOobject
1660         (
1661             "pointLevel",
1662             mesh_.time().timeName(),
1663             mesh_,
1664             IOobject::NO_READ,
1665             IOobject::NO_WRITE,
1666             false
1667         ),
1668         pMesh,
1669         dimensionedScalar("zero", dimless, 0)
1670     );
1671
1672     const labelList& pointLevel = meshCutter_.pointLevel();
1673
1674     forAll(pointRefLevel, pointI)
1675     {
1676         pointRefLevel[pointI] = pointLevel[pointI];
1677     }
1678
1679     pointRefLevel.write();
1680 }
1681
1682
1683 void Foam::meshRefinement::dumpIntersections(const fileName& prefix) const
1684 {
1685     {
1686         const pointField& cellCentres = mesh_.cellCentres();
1687
1688         OFstream str(prefix + "_edges.obj");
1689         label vertI = 0;
1690         Pout<< "meshRefinement::dumpIntersections :"
1691             << " Writing cellcentre-cellcentre intersections to file "
1692             << str.name() << endl;
1693
1694
1695         // Redo all intersections
1696         // ~~~~~~~~~~~~~~~~~~~~~~
1697
1698         // Get boundary face centre and level. Coupled aware.
1699         labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
1700         pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
1701         calcNeighbourData(neiLevel, neiCc);
1702
1703         labelList intersectionFaces(intersectedFaces());
1704
1705         // Collect segments we want to test for
1706         pointField start(intersectionFaces.size());
1707         pointField end(intersectionFaces.size());
1708
1709         forAll(intersectionFaces, i)
1710         {
1711             label faceI = intersectionFaces[i];
1712             start[i] = cellCentres[mesh_.faceOwner()[faceI]];
1713
1714             if (mesh_.isInternalFace(faceI))
1715             {
1716                 end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
1717             }
1718             else
1719             {
1720                 end[i] = neiCc[faceI-mesh_.nInternalFaces()];
1721             }
1722         }
1723
1724         // Do tests in one go
1725         labelList surfaceHit;
1726         List<pointIndexHit> surfaceHitInfo;
1727         surfaces_.findAnyIntersection
1728         (
1729             start,
1730             end,
1731             surfaceHit,
1732             surfaceHitInfo
1733         );
1734
1735         forAll(intersectionFaces, i)
1736         {
1737             if (surfaceHit[i] != -1)
1738             {
1739                 meshTools::writeOBJ(str, start[i]);
1740                 vertI++;
1741                 meshTools::writeOBJ(str, surfaceHitInfo[i].hitPoint());
1742                 vertI++;
1743                 meshTools::writeOBJ(str, end[i]);
1744                 vertI++;
1745                 str << "l " << vertI-2 << ' ' << vertI-1 << nl
1746                     << "l " << vertI-1 << ' ' << vertI << nl;
1747             }
1748         }
1749     }
1750
1751     // Convert to vtk format
1752     string cmd
1753     (
1754         "objToVTK " + prefix + "_edges.obj " + prefix + "_edges.vtk > /dev/null"
1755     );
1756     system(cmd.c_str());
1757
1758     Pout<< endl;
1759 }
1760
1761
1762 void Foam::meshRefinement::write
1763 (
1764     const label flag,
1765     const fileName& prefix
1766 ) const
1767 {
1768     if (flag & MESH)
1769     {
1770         write();
1771     }
1772     if (flag & SCALARLEVELS)
1773     {
1774         dumpRefinementLevel();
1775     }
1776     if (flag&OBJINTERSECTIONS && prefix.size()>0)
1777     {
1778         dumpIntersections(prefix);
1779     }
1780 }
1781
1782
1783 // ************************************************************************* //