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