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-2009 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 #include "geomDecomp.H"
  54 #include "Random.H"
  55 #include "searchableSurfaces.H"
  56 #include "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     labelHashSet addedPatchIDSet(meshedPatches());
  88
  89     forAll(patches, patchI)
  90     {
  91         const polyPatch& pp = patches[patchI];
  92
  93         const unallocLabelList& faceCells = pp.faceCells();
  94         const vectorField::subField faceCentres = pp.faceCentres();
  95         const vectorField::subField faceAreas = pp.faceAreas();
  96
  97         label bFaceI = pp.start()-mesh_.nInternalFaces();
  98
  99         if (pp.coupled())
 100         {
 101             forAll(faceCells, i)
 102             {
 103                 neiLevel[bFaceI] = cellLevel[faceCells[i]];
 104                 neiCc[bFaceI] = cellCentres[faceCells[i]];
 105                 bFaceI++;
 106             }
 107         }
 108         else if (addedPatchIDSet.found(patchI))
 109         {
 110             // Face was introduced from cell-cell intersection. Try to
 111             // reconstruct other side cell(centre). Three possibilities:
 112             // - cells same size.
 113             // - preserved cell smaller. Not handled.
 114             // - preserved cell larger.
 115             forAll(faceCells, i)
 116             {
 117                 // Extrapolate the face centre.
 118                 vector fn = faceAreas[i];
 119                 fn /= mag(fn)+VSMALL;
 120
 121                 label own = faceCells[i];
 122                 label ownLevel = cellLevel[own];
 123                 label faceLevel = meshCutter_.getAnchorLevel(pp.start()+i);
 124
 125                 // Normal distance from face centre to cell centre
 126                 scalar d = ((faceCentres[i] - cellCentres[own]) & fn);
 127                 if (faceLevel > ownLevel)
 128                 {
 129                     // Other cell more refined. Adjust normal distance
 130                     d *= 0.5;
 131                 }
 132                 neiLevel[bFaceI] = cellLevel[ownLevel];
 133                 // Calculate other cell centre by extrapolation
 134                 neiCc[bFaceI] = faceCentres[i] + d*fn;
 135                 bFaceI++;
 136             }
 137         }
 138         else
 139         {
 140             forAll(faceCells, i)
 141             {
 142                 neiLevel[bFaceI] = cellLevel[faceCells[i]];
 143                 neiCc[bFaceI] = faceCentres[i];
 144                 bFaceI++;
 145             }
 146         }
 147     }
 148
 149     // Swap coupled boundaries. Apply separation to cc since is coordinate.
 150     syncTools::swapBoundaryFaceList(mesh_, neiCc, true);
 151     syncTools::swapBoundaryFaceList(mesh_, neiLevel, false);
 152 }
 153
 154
 155 // Find intersections of edges (given by their two endpoints) with surfaces.
 156 // Returns first intersection if there are more than one.
 157 void Foam::meshRefinement::updateIntersections(const labelList& changedFaces)
 158 {
 159     const pointField& cellCentres = mesh_.cellCentres();
 160
 161     // Stats on edges to test. Count proc faces only once.
 162     PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
 163
 164     {
 165         label nMasterFaces = 0;
 166         forAll(isMasterFace, faceI)
 167         {
 168             if (isMasterFace.get(faceI) == 1)
 169             {
 170                 nMasterFaces++;
 171             }
 172         }
 173         reduce(nMasterFaces, sumOp<label>());
 174
 175         label nChangedFaces = 0;
 176         forAll(changedFaces, i)
 177         {
 178             if (isMasterFace.get(changedFaces[i]) == 1)
 179             {
 180                 nChangedFaces++;
 181             }
 182         }
 183         reduce(nChangedFaces, sumOp<label>());
 184
 185         Info<< "Edge intersection testing:" << nl
 186             << "    Number of edges             : " << nMasterFaces << nl
 187             << "    Number of edges to retest   : " << nChangedFaces
 188             << endl;
 189     }
 190
 191
 192     // Get boundary face centre and level. Coupled aware.
 193     labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
 194     pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
 195     calcNeighbourData(neiLevel, neiCc);
 196
 197     // Collect segments we want to test for
 198     pointField start(changedFaces.size());
 199     pointField end(changedFaces.size());
 200
 201     forAll(changedFaces, i)
 202     {
 203         label faceI = changedFaces[i];
 204         label own = mesh_.faceOwner()[faceI];
 205
 206         start[i] = cellCentres[own];
 207         if (mesh_.isInternalFace(faceI))
 208         {
 209             end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
 210         }
 211         else
 212         {
 213             end[i] = neiCc[faceI-mesh_.nInternalFaces()];
 214         }
 215     }
 216
 217     // Do tests in one go
 218     labelList surfaceHit;
 219     {
 220         labelList surfaceLevel;
 221         surfaces_.findHigherIntersection
 222         (
 223             start,
 224             end,
 225             labelList(start.size(), -1),    // accept any intersection
 226             surfaceHit,
 227             surfaceLevel
 228         );
 229     }
 230
 231     // Keep just surface hit
 232     forAll(surfaceHit, i)
 233     {
 234         surfaceIndex_[changedFaces[i]] = surfaceHit[i];
 235     }
 236
 237     // Make sure both sides have same information. This should be
 238     // case in general since same vectors but just to make sure.
 239     syncTools::syncFaceList(mesh_, surfaceIndex_, maxEqOp<label>(), false);
 240
 241     label nHits = countHits();
 242     label nTotHits = returnReduce(nHits, sumOp<label>());
 243
 244     Info<< "    Number of intersected edges : " << nTotHits << endl;
 245
 246     // Set files to same time as mesh
 247     setInstance(mesh_.facesInstance());
 248 }
 249
 250
 251 void Foam::meshRefinement::checkData()
 252 {
 253     Pout<< "meshRefinement::checkData() : Checking refinement structure."
 254         << endl;
 255     meshCutter_.checkMesh();
 256
 257     Pout<< "meshRefinement::checkData() : Checking refinement levels."
 258         << endl;
 259     meshCutter_.checkRefinementLevels(1, labelList(0));
 260
 261
 262     Pout<< "meshRefinement::checkData() : Checking synchronization."
 263         << endl;
 264
 265     // Check face centres
 266     {
 267         // Boundary face centres
 268         pointField::subList boundaryFc
 269         (
 270             mesh_.faceCentres(),
 271             mesh_.nFaces()-mesh_.nInternalFaces(),
 272             mesh_.nInternalFaces()
 273         );
 274
 275         // Get neighbouring face centres
 276         pointField neiBoundaryFc(boundaryFc);
 277         syncTools::swapBoundaryFaceList
 278         (
 279             mesh_,
 280             neiBoundaryFc,
 281             true
 282         );
 283
 284         // Compare
 285         testSyncBoundaryFaceList
 286         (
 287             mergeDistance_,
 288             "testing faceCentres : ",
 289             boundaryFc,
 290             neiBoundaryFc
 291         );
 292     }
 293     // Check meshRefinement
 294     {
 295         // Get boundary face centre and level. Coupled aware.
 296         labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
 297         pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
 298         calcNeighbourData(neiLevel, neiCc);
 299
 300         // Collect segments we want to test for
 301         pointField start(mesh_.nFaces());
 302         pointField end(mesh_.nFaces());
 303
 304         forAll(start, faceI)
 305         {
 306             start[faceI] = mesh_.cellCentres()[mesh_.faceOwner()[faceI]];
 307
 308             if (mesh_.isInternalFace(faceI))
 309             {
 310                 end[faceI] = mesh_.cellCentres()[mesh_.faceNeighbour()[faceI]];
 311             }
 312             else
 313             {
 314                 end[faceI] = neiCc[faceI-mesh_.nInternalFaces()];
 315             }
 316         }
 317
 318         // Do tests in one go
 319         labelList surfaceHit;
 320         {
 321             labelList surfaceLevel;
 322             surfaces_.findHigherIntersection
 323             (
 324                 start,
 325                 end,
 326                 labelList(start.size(), -1),    // accept any intersection
 327                 surfaceHit,
 328                 surfaceLevel
 329             );
 330         }
 331
 332         // Check
 333         forAll(surfaceHit, faceI)
 334         {
 335             if (surfaceHit[faceI] != surfaceIndex_[faceI])
 336             {
 337                 if (mesh_.isInternalFace(faceI))
 338                 {
 339                     WarningIn("meshRefinement::checkData()")
 340                         << "Internal face:" << faceI
 341                         << " fc:" << mesh_.faceCentres()[faceI]
 342                         << " cached surfaceIndex_:" << surfaceIndex_[faceI]
 343                         << " current:" << surfaceHit[faceI]
 344                         << " ownCc:"
 345                         << mesh_.cellCentres()[mesh_.faceOwner()[faceI]]
 346                         << " neiCc:"
 347                         << mesh_.cellCentres()[mesh_.faceNeighbour()[faceI]]
 348                         << endl;
 349                 }
 350                 else
 351                 {
 352                     WarningIn("meshRefinement::checkData()")
 353                         << "Boundary face:" << faceI
 354                         << " fc:" << mesh_.faceCentres()[faceI]
 355                         << " cached surfaceIndex_:" << surfaceIndex_[faceI]
 356                         << " current:" << surfaceHit[faceI]
 357                         << " ownCc:"
 358                         << mesh_.cellCentres()[mesh_.faceOwner()[faceI]]
 359                         << endl;
 360                 }
 361             }
 362         }
 363     }
 364     {
 365         labelList::subList boundarySurface
 366         (
 367             surfaceIndex_,
 368             mesh_.nFaces()-mesh_.nInternalFaces(),
 369             mesh_.nInternalFaces()
 370         );
 371
 372         labelList neiBoundarySurface(boundarySurface);
 373         syncTools::swapBoundaryFaceList
 374         (
 375             mesh_,
 376             neiBoundarySurface,
 377             false
 378         );
 379
 380         // Compare
 381         testSyncBoundaryFaceList
 382         (
 383             0,                              // tolerance
 384             "testing surfaceIndex() : ",
 385             boundarySurface,
 386             neiBoundarySurface
 387         );
 388     }
 389
 390
 391     // Find duplicate faces
 392     Pout<< "meshRefinement::checkData() : Counting duplicate faces."
 393         << endl;
 394
 395     labelList duplicateFace
 396     (
 397         localPointRegion::findDuplicateFaces
 398         (
 399             mesh_,
 400             identity(mesh_.nFaces()-mesh_.nInternalFaces())
 401           + mesh_.nInternalFaces()
 402         )
 403     );
 404
 405     // Count
 406     {
 407         label nDup = 0;
 408
 409         forAll(duplicateFace, i)
 410         {
 411             if (duplicateFace[i] != -1)
 412             {
 413                 nDup++;
 414             }
 415         }
 416         nDup /= 2;  // will have counted both faces of duplicate
 417         Pout<< "meshRefinement::checkData() : Found " << nDup
 418             << " duplicate pairs of faces." << endl;
 419     }
 420 }
 421
 422
 423 void Foam::meshRefinement::setInstance(const fileName& inst)
 424 {
 425     meshCutter_.setInstance(inst);
 426     surfaceIndex_.instance() = inst;
 427 }
 428
 429
 430 // Remove cells. Put exposedFaces (output of getExposedFaces(cellsToRemove))
 431 // into exposedPatchIDs.
 432 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRemoveCells
 433 (
 434     const labelList& cellsToRemove,
 435     const labelList& exposedFaces,
 436     const labelList& exposedPatchIDs,
 437     removeCells& cellRemover
 438 )
 439 {
 440     polyTopoChange meshMod(mesh_);
 441
 442     // Arbitrary: put exposed faces into last patch.
 443     cellRemover.setRefinement
 444     (
 445         cellsToRemove,
 446         exposedFaces,
 447         exposedPatchIDs,
 448         meshMod
 449     );
 450
 451     // Change the mesh (no inflation)
 452     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
 453
 454     // Update fields
 455     mesh_.updateMesh(map);
 456
 457     // Move mesh (since morphing might not do this)
 458     if (map().hasMotionPoints())
 459     {
 460         mesh_.movePoints(map().preMotionPoints());
 461     }
 462     else
 463     {
 464         // Delete mesh volumes. No other way to do this?
 465         mesh_.clearOut();
 466     }
 467
 468     if (overwrite_)
 469     {
 470         mesh_.setInstance(oldInstance_);
 471     }
 472
 473     // Update local mesh data
 474     cellRemover.updateMesh(map);
 475
 476     // Update intersections. Recalculate intersections for exposed faces.
 477     labelList newExposedFaces = renumber
 478     (
 479         map().reverseFaceMap(),
 480         exposedFaces
 481     );
 482
 483     //Pout<< "removeCells : updating intersections for "
 484     //    << newExposedFaces.size() << " newly exposed faces." << endl;
 485
 486     updateMesh(map, newExposedFaces);
 487
 488     return map;
 489 }
 490
 491
 492 // Determine for multi-processor regions the lowest numbered cell on the lowest
 493 // numbered processor.
 494 void Foam::meshRefinement::getCoupledRegionMaster
 495 (
 496     const globalIndex& globalCells,
 497     const boolList& blockedFace,
 498     const regionSplit& globalRegion,
 499     Map<label>& regionToMaster
 500 ) const
 501 {
 502     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 503
 504     forAll(patches, patchI)
 505     {
 506         const polyPatch& pp = patches[patchI];
 507
 508         if (isA<processorPolyPatch>(pp))
 509         {
 510             forAll(pp, i)
 511             {
 512                 label faceI = pp.start()+i;
 513
 514                 if (!blockedFace[faceI])
 515                 {
 516                     // Only if there is a connection to the neighbour
 517                     // will there be a multi-domain region; if not through
 518                     // this face then through another.
 519
 520                     label cellI = mesh_.faceOwner()[faceI];
 521                     label globalCellI = globalCells.toGlobal(cellI);
 522
 523                     Map<label>::iterator iter =
 524                         regionToMaster.find(globalRegion[cellI]);
 525
 526                     if (iter != regionToMaster.end())
 527                     {
 528                         label master = iter();
 529                         iter() = min(master, globalCellI);
 530                     }
 531                     else
 532                     {
 533                         regionToMaster.insert
 534                         (
 535                             globalRegion[cellI],
 536                             globalCellI
 537                         );
 538                     }
 539                 }
 540             }
 541         }
 542     }
 543
 544     // Do reduction
 545     Pstream::mapCombineGather(regionToMaster, minEqOp<label>());
 546     Pstream::mapCombineScatter(regionToMaster);
 547 }
 548
 549
 550 void Foam::meshRefinement::calcLocalRegions
 551 (
 552     const globalIndex& globalCells,
 553     const labelList& globalRegion,
 554     const Map<label>& coupledRegionToMaster,
 555
 556     Map<label>& globalToLocalRegion,
 557     pointField& localPoints
 558 ) const
 559 {
 560     globalToLocalRegion.resize(globalRegion.size());
 561     DynamicList<point> localCc(globalRegion.size()/2);
 562
 563     forAll(globalRegion, cellI)
 564     {
 565         Map<label>::const_iterator fndMaster =
 566             coupledRegionToMaster.find(globalRegion[cellI]);
 567
 568         if (fndMaster != coupledRegionToMaster.end())
 569         {
 570             // Multi-processor region.
 571             if (globalCells.toGlobal(cellI) == fndMaster())
 572             {
 573                 // I am master. Allocate region for me.
 574                 globalToLocalRegion.insert(globalRegion[cellI], localCc.size());
 575                 localCc.append(mesh_.cellCentres()[cellI]);
 576             }
 577         }
 578         else
 579         {
 580             // Single processor region.
 581             if (globalToLocalRegion.insert(globalRegion[cellI], localCc.size()))
 582             {
 583                 localCc.append(mesh_.cellCentres()[cellI]);
 584             }
 585         }
 586     }
 587
 588     localPoints.transfer(localCc);
 589
 590     if (localPoints.size() != globalToLocalRegion.size())
 591     {
 592         FatalErrorIn("calcLocalRegions(..)")
 593             << "localPoints:" << localPoints.size()
 594             << " globalToLocalRegion:" << globalToLocalRegion.size()
 595             << exit(FatalError);
 596     }
 597 }
 598
 599
 600 Foam::label Foam::meshRefinement::getShiftedRegion
 601 (
 602     const globalIndex& indexer,
 603     const Map<label>& globalToLocalRegion,
 604     const Map<label>& coupledRegionToShifted,
 605     const label globalRegion
 606 )
 607 {
 608     Map<label>::const_iterator iter =
 609         globalToLocalRegion.find(globalRegion);
 610
 611     if (iter != globalToLocalRegion.end())
 612     {
 613         // Region is 'owned' locally. Convert local region index into global.
 614         return indexer.toGlobal(iter());
 615     }
 616     else
 617     {
 618         return coupledRegionToShifted[globalRegion];
 619     }
 620 }
 621
 622
 623 // Add if not yet present
 624 void Foam::meshRefinement::addUnique(const label elem, labelList& lst)
 625 {
 626     if (findIndex(lst, elem) == -1)
 627     {
 628         label sz = lst.size();
 629         lst.setSize(sz+1);
 630         lst[sz] = elem;
 631     }
 632 }
 633
 634
 635 void Foam::meshRefinement::calcRegionRegions
 636 (
 637     const labelList& globalRegion,
 638     const Map<label>& globalToLocalRegion,
 639     const Map<label>& coupledRegionToMaster,
 640     labelListList& regionRegions
 641 ) const
 642 {
 643     // Global region indexing since we now know the shifted regions.
 644     globalIndex shiftIndexer(globalToLocalRegion.size());
 645
 646     // Redo the coupledRegionToMaster to be in shifted region indexing.
 647     Map<label> coupledRegionToShifted(coupledRegionToMaster.size());
 648     forAllConstIter(Map<label>, coupledRegionToMaster, iter)
 649     {
 650         label region = iter.key();
 651
 652         Map<label>::const_iterator fndRegion = globalToLocalRegion.find(region);
 653
 654         if (fndRegion != globalToLocalRegion.end())
 655         {
 656             // A local cell is master of this region. Get its shifted region.
 657             coupledRegionToShifted.insert
 658             (
 659                 region,
 660                 shiftIndexer.toGlobal(fndRegion())
 661             );
 662         }
 663         Pstream::mapCombineGather(coupledRegionToShifted, minEqOp<label>());
 664         Pstream::mapCombineScatter(coupledRegionToShifted);
 665     }
 666
 667
 668     // Determine region-region connectivity.
 669     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 670     // This is for all my regions (so my local ones or the ones I am master
 671     // of) the neighbouring region indices.
 672
 673
 674     // Transfer lists.
 675     PtrList<HashSet<edge, Hash<edge> > > regionConnectivity(Pstream::nProcs());
 676     forAll(regionConnectivity, procI)
 677     {
 678         if (procI != Pstream::myProcNo())
 679         {
 680             regionConnectivity.set
 681             (
 682                 procI,
 683                 new HashSet<edge, Hash<edge> >
 684                 (
 685                     coupledRegionToShifted.size()
 686                   / Pstream::nProcs()
 687                 )
 688             );
 689         }
 690     }
 691
 692
 693     // Connectivity. For all my local regions the connected regions.
 694     regionRegions.setSize(globalToLocalRegion.size());
 695
 696     // Add all local connectivity to regionRegions; add all non-local
 697     // to the transferlists.
 698     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
 699     {
 700         label ownRegion = globalRegion[mesh_.faceOwner()[faceI]];
 701         label neiRegion = globalRegion[mesh_.faceNeighbour()[faceI]];
 702
 703         if (ownRegion != neiRegion)
 704         {
 705             label shiftOwnRegion = getShiftedRegion
 706             (
 707                 shiftIndexer,
 708                 globalToLocalRegion,
 709                 coupledRegionToShifted,
 710                 ownRegion
 711             );
 712             label shiftNeiRegion = getShiftedRegion
 713             (
 714                 shiftIndexer,
 715                 globalToLocalRegion,
 716                 coupledRegionToShifted,
 717                 neiRegion
 718             );
 719
 720
 721             // Connection between two regions. Send to owner of region.
 722             // - is ownRegion 'owned' by me
 723             // - is neiRegion 'owned' by me
 724
 725             if (shiftIndexer.isLocal(shiftOwnRegion))
 726             {
 727                 label localI = shiftIndexer.toLocal(shiftOwnRegion);
 728                 addUnique(shiftNeiRegion, regionRegions[localI]);
 729             }
 730             else
 731             {
 732                 label masterProc = shiftIndexer.whichProcID(shiftOwnRegion);
 733                 regionConnectivity[masterProc].insert
 734                 (
 735                     edge(shiftOwnRegion, shiftNeiRegion)
 736                 );
 737             }
 738
 739             if (shiftIndexer.isLocal(shiftNeiRegion))
 740             {
 741                 label localI = shiftIndexer.toLocal(shiftNeiRegion);
 742                 addUnique(shiftOwnRegion, regionRegions[localI]);
 743             }
 744             else
 745             {
 746                 label masterProc = shiftIndexer.whichProcID(shiftNeiRegion);
 747                 regionConnectivity[masterProc].insert
 748                 (
 749                     edge(shiftOwnRegion, shiftNeiRegion)
 750                 );
 751             }
 752         }
 753     }
 754
 755
 756     // Send
 757     forAll(regionConnectivity, procI)
 758     {
 759         if (procI != Pstream::myProcNo())
 760         {
 761             OPstream str(Pstream::blocking, procI);
 762             str << regionConnectivity[procI];
 763         }
 764     }
 765     // Receive
 766     forAll(regionConnectivity, procI)
 767     {
 768         if (procI != Pstream::myProcNo())
 769         {
 770             IPstream str(Pstream::blocking, procI);
 771             str >> regionConnectivity[procI];
 772         }
 773     }
 774
 775     // Add to addressing.
 776     forAll(regionConnectivity, procI)
 777     {
 778         if (procI != Pstream::myProcNo())
 779         {
 780             for
 781             (
 782                 HashSet<edge, Hash<edge> >::const_iterator iter =
 783                     regionConnectivity[procI].begin();
 784                 iter != regionConnectivity[procI].end();
 785                 ++iter
 786             )
 787             {
 788                 const edge& e = iter.key();
 789
 790                 bool someLocal = false;
 791                 if (shiftIndexer.isLocal(e[0]))
 792                 {
 793                     label localI = shiftIndexer.toLocal(e[0]);
 794                     addUnique(e[1], regionRegions[localI]);
 795                     someLocal = true;
 796                 }
 797                 if (shiftIndexer.isLocal(e[1]))
 798                 {
 799                     label localI = shiftIndexer.toLocal(e[1]);
 800                     addUnique(e[0], regionRegions[localI]);
 801                     someLocal = true;
 802                 }
 803
 804                 if (!someLocal)
 805                 {
 806                     FatalErrorIn("calcRegionRegions(..)")
 807                         << "Received from processor " << procI
 808                         << " connection " << e
 809                         << " where none of the elements is local to me."
 810                         << abort(FatalError);
 811                 }
 812             }
 813         }
 814     }
 815 }
 816
 817
 818 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 819
 820 // Construct from components
 821 Foam::meshRefinement::meshRefinement
 822 (
 823     fvMesh& mesh,
 824     const scalar mergeDistance,
 825     const bool overwrite,
 826     const refinementSurfaces& surfaces,
 827     const shellSurfaces& shells
 828 )
 829 :
 830     mesh_(mesh),
 831     mergeDistance_(mergeDistance),
 832     overwrite_(overwrite),
 833     oldInstance_(mesh.pointsInstance()),
 834     surfaces_(surfaces),
 835     shells_(shells),
 836     meshCutter_
 837     (
 838         mesh,
 839         labelIOList
 840         (
 841             IOobject
 842             (
 843                 "cellLevel",
 844                 mesh_.facesInstance(),
 845                 fvMesh::meshSubDir,
 846                 mesh,
 847                 IOobject::READ_IF_PRESENT,
 848                 IOobject::NO_WRITE,
 849                 false
 850             ),
 851             labelList(mesh_.nCells(), 0)
 852         ),
 853         labelIOList
 854         (
 855             IOobject
 856             (
 857                 "pointLevel",
 858                 mesh_.facesInstance(),
 859                 fvMesh::meshSubDir,
 860                 mesh,
 861                 IOobject::READ_IF_PRESENT,
 862                 IOobject::NO_WRITE,
 863                 false
 864             ),
 865             labelList(mesh_.nPoints(), 0)
 866         ),
 867         refinementHistory
 868         (
 869             IOobject
 870             (
 871                 "refinementHistory",
 872                 mesh_.facesInstance(),
 873                 fvMesh::meshSubDir,
 874                 mesh_,
 875                 IOobject::NO_READ,
 876                 IOobject::NO_WRITE,
 877                 false
 878             ),
 879             List<refinementHistory::splitCell8>(0),
 880             labelList(0)
 881         )   // no unrefinement
 882     ),
 883     surfaceIndex_
 884     (
 885         IOobject
 886         (
 887             "surfaceIndex",
 888             mesh_.facesInstance(),
 889             fvMesh::meshSubDir,
 890             mesh_,
 891             IOobject::NO_READ,
 892             IOobject::NO_WRITE,
 893             false
 894         ),
 895         labelList(mesh_.nFaces(), -1)
 896     ),
 897     userFaceData_(0)
 898 {
 899     // recalculate intersections for all faces
 900     updateIntersections(identity(mesh_.nFaces()));
 901 }
 902
 903
 904 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 905
 906 Foam::label Foam::meshRefinement::countHits() const
 907 {
 908     // Stats on edges to test. Count proc faces only once.
 909     PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
 910
 911     label nHits = 0;
 912
 913     forAll(surfaceIndex_, faceI)
 914     {
 915         if (surfaceIndex_[faceI] >= 0 && isMasterFace.get(faceI) == 1)
 916         {
 917             nHits++;
 918         }
 919     }
 920     return nHits;
 921 }
 922
 923
 924 // Determine distribution to move connected regions onto one processor.
 925 Foam::labelList Foam::meshRefinement::decomposeCombineRegions
 926 (
 927     const boolList& blockedFace,
 928     const List<labelPair>& explicitConnections,
 929     decompositionMethod& decomposer
 930 ) const
 931 {
 932     // Determine global regions, separated by blockedFaces
 933     regionSplit globalRegion(mesh_, blockedFace, explicitConnections);
 934
 935     // Now globalRegion has global region per cell. Problem is that
 936     // the region might span multiple domains so we want to get
 937     // a "region master" per domain. Note that multi-processor
 938     // regions can only occur on cells on coupled patches.
 939     // Note: since the number of regions does not change by this the
 940     // process can be seen as just a shift of a region onto a single
 941     // processor.
 942
 943
 944     // Global cell numbering engine
 945     globalIndex globalCells(mesh_.nCells());
 946
 947
 948     // Determine per coupled region the master cell (lowest numbered cell
 949     // on lowest numbered processor)
 950     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 951     // (does not determine master for non-coupled=fully-local regions)
 952
 953     Map<label> coupledRegionToMaster(mesh_.nFaces()-mesh_.nInternalFaces());
 954     getCoupledRegionMaster
 955     (
 956         globalCells,
 957         blockedFace,
 958         globalRegion,
 959         coupledRegionToMaster
 960     );
 961
 962     // Determine my regions
 963     // ~~~~~~~~~~~~~~~~~~~~
 964     // These are the fully local ones or the coupled ones of which I am master.
 965
 966     Map<label> globalToLocalRegion;
 967     pointField localPoints;
 968     calcLocalRegions
 969     (
 970         globalCells,
 971         globalRegion,
 972         coupledRegionToMaster,
 973
 974         globalToLocalRegion,
 975         localPoints
 976     );
 977
 978
 979
 980     // Find distribution for regions
 981     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 982
 983     labelList regionDistribution;
 984
 985     if (isA<geomDecomp>(decomposer))
 986     {
 987         regionDistribution = decomposer.decompose(localPoints);
 988     }
 989     else
 990     {
 991         labelListList regionRegions;
 992         calcRegionRegions
 993         (
 994             globalRegion,
 995             globalToLocalRegion,
 996             coupledRegionToMaster,
 997
 998             regionRegions
 999         );
1000
1001         regionDistribution = decomposer.decompose(regionRegions, localPoints);
1002     }
1003
1004
1005
1006     // Convert region-based decomposition back to cell-based one
1007     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1008
1009     // Transfer destination processor back to all. Use global reduce for now.
1010     Map<label> regionToDist(coupledRegionToMaster.size());
1011     forAllConstIter(Map<label>, coupledRegionToMaster, iter)
1012     {
1013         label region = iter.key();
1014
1015         Map<label>::const_iterator regionFnd = globalToLocalRegion.find(region);
1016
1017         if (regionFnd != globalToLocalRegion.end())
1018         {
1019             // Master cell is local. Store my distribution.
1020             regionToDist.insert(iter.key(), regionDistribution[regionFnd()]);
1021         }
1022         else
1023         {
1024             // Master cell is not on this processor. Make sure it is overridden.
1025             regionToDist.insert(iter.key(), labelMax);
1026         }
1027     }
1028     Pstream::mapCombineGather(regionToDist, minEqOp<label>());
1029     Pstream::mapCombineScatter(regionToDist);
1030
1031
1032     // Determine destination for all cells
1033     labelList distribution(mesh_.nCells());
1034
1035     forAll(globalRegion, cellI)
1036     {
1037         Map<label>::const_iterator fndRegion =
1038             regionToDist.find(globalRegion[cellI]);
1039
1040         if (fndRegion != regionToDist.end())
1041         {
1042             distribution[cellI] = fndRegion();
1043         }
1044         else
1045         {
1046             // region is local to the processor.
1047             label localRegionI = globalToLocalRegion[globalRegion[cellI]];
1048
1049             distribution[cellI] = regionDistribution[localRegionI];
1050         }
1051     }
1052
1053     return distribution;
1054 }
1055
1056
1057 Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
1058 (
1059     const bool keepZoneFaces,
1060     const bool keepBaffles,
1061     decompositionMethod& decomposer,
1062     fvMeshDistribute& distributor
1063 )
1064 {
1065     autoPtr<mapDistributePolyMesh> map;
1066
1067     if (Pstream::parRun())
1068     {
1069         //if (debug_)
1070         //{
1071         //    const_cast<Time&>(mesh_.time())++;
1072         //}
1073
1074         // Wanted distribution
1075         labelList distribution;
1076
1077         if (keepZoneFaces || keepBaffles)
1078         {
1079             // Faces where owner and neighbour are not 'connected' so can
1080             // go to different processors.
1081             boolList blockedFace(mesh_.nFaces(), true);
1082             label nUnblocked = 0;
1083             // Pairs of baffles
1084             List<labelPair> couples;
1085
1086             if (keepZoneFaces)
1087             {
1088                 // Determine decomposition to keep/move surface zones
1089                 // on one processor. The reason is that snapping will make these
1090                 // into baffles, move and convert them back so if they were
1091                 // proc boundaries after baffling&moving the points might be no
1092                 // longer snychronised so recoupling will fail. To prevent this
1093                 // keep owner&neighbour of such a surface zone on the same
1094                 // processor.
1095
1096                 const wordList& fzNames = surfaces().faceZoneNames();
1097                 const faceZoneMesh& fZones = mesh_.faceZones();
1098
1099                 // Get faces whose owner and neighbour should stay together,
1100                 // i.e. they are not 'blocked'.
1101
1102                 forAll(fzNames, surfI)
1103                 {
1104                     if (fzNames[surfI].size())
1105                     {
1106                         // Get zone
1107                         label zoneI = fZones.findZoneID(fzNames[surfI]);
1108
1109                         const faceZone& fZone = fZones[zoneI];
1110
1111                         forAll(fZone, i)
1112                         {
1113                             if (blockedFace[fZone[i]])
1114                             {
1115                                 blockedFace[fZone[i]] = false;
1116                                 nUnblocked++;
1117                             }
1118                         }
1119                     }
1120                 }
1121
1122
1123                 // If the faceZones are not synchronised the blockedFace
1124                 // might not be synchronised. If you are sure the faceZones
1125                 // are synchronised remove below check.
1126                 syncTools::syncFaceList
1127                 (
1128                     mesh_,
1129                     blockedFace,
1130                     andEqOp<bool>(),    // combine operator
1131                     false               // separation
1132                 );
1133             }
1134             reduce(nUnblocked, sumOp<label>());
1135
1136             if (keepZoneFaces)
1137             {
1138                 Info<< "Found " << nUnblocked
1139                     << " zoned faces to keep together." << endl;
1140             }
1141
1142             if (keepBaffles)
1143             {
1144                 // Get boundary baffles that need to stay together.
1145                 couples = getDuplicateFaces   // all baffles
1146                 (
1147                     identity(mesh_.nFaces()-mesh_.nInternalFaces())
1148                    +mesh_.nInternalFaces()
1149                 );
1150             }
1151             label nCouples = returnReduce(couples.size(), sumOp<label>());
1152
1153             if (keepBaffles)
1154             {
1155                 Info<< "Found " << nCouples << " baffles to keep together."
1156                     << endl;
1157             }
1158
1159             if (nUnblocked > 0 || nCouples > 0)
1160             {
1161                 Info<< "Applying special decomposition to keep baffles"
1162                     << " and zoned faces together." << endl;
1163
1164                 distribution = decomposeCombineRegions
1165                 (
1166                     blockedFace,
1167                     couples,
1168                     decomposer
1169                 );
1170
1171                 labelList nProcCells(distributor.countCells(distribution));
1172                 Pstream::listCombineGather(nProcCells, plusEqOp<label>());
1173                 Pstream::listCombineScatter(nProcCells);
1174
1175                 Info<< "Calculated decomposition:" << endl;
1176                 forAll(nProcCells, procI)
1177                 {
1178                     Info<< "    " << procI << '\t' << nProcCells[procI] << endl;
1179                 }
1180                 Info<< endl;
1181             }
1182             else
1183             {
1184                 // Normal decomposition
1185                 distribution = decomposer.decompose(mesh_.cellCentres());
1186             }
1187         }
1188         else
1189         {
1190             // Normal decomposition
1191             distribution = decomposer.decompose(mesh_.cellCentres());
1192         }
1193
1194         if (debug)
1195         {
1196             labelList nProcCells(distributor.countCells(distribution));
1197             Pout<< "Wanted distribution:" << nProcCells << endl;
1198
1199             Pstream::listCombineGather(nProcCells, plusEqOp<label>());
1200             Pstream::listCombineScatter(nProcCells);
1201
1202             Pout<< "Wanted resulting decomposition:" << endl;
1203             forAll(nProcCells, procI)
1204             {
1205                 Pout<< "    " << procI << '\t' << nProcCells[procI] << endl;
1206             }
1207             Pout<< endl;
1208         }
1209         // Do actual sending/receiving of mesh
1210         map = distributor.distribute(distribution);
1211
1212         // Update numbering of meshRefiner
1213         distribute(map);
1214     }
1215     return map;
1216 }
1217
1218
1219 // Helper function to get intersected faces
1220 Foam::labelList Foam::meshRefinement::intersectedFaces() const
1221 {
1222     label nBoundaryFaces = 0;
1223
1224     forAll(surfaceIndex_, faceI)
1225     {
1226         if (surfaceIndex_[faceI] != -1)
1227         {
1228             nBoundaryFaces++;
1229         }
1230     }
1231
1232     labelList surfaceFaces(nBoundaryFaces);
1233     nBoundaryFaces = 0;
1234
1235     forAll(surfaceIndex_, faceI)
1236     {
1237         if (surfaceIndex_[faceI] != -1)
1238         {
1239             surfaceFaces[nBoundaryFaces++] = faceI;
1240         }
1241     }
1242     return surfaceFaces;
1243 }
1244
1245
1246 // Helper function to get points used by faces
1247 Foam::labelList Foam::meshRefinement::intersectedPoints() const
1248 {
1249     const faceList& faces = mesh_.faces();
1250
1251     // Mark all points on faces that will become baffles
1252     PackedBoolList isBoundaryPoint(mesh_.nPoints(), 0u);
1253     label nBoundaryPoints = 0;
1254
1255     forAll(surfaceIndex_, faceI)
1256     {
1257         if (surfaceIndex_[faceI] != -1)
1258         {
1259             const face& f = faces[faceI];
1260
1261             forAll(f, fp)
1262             {
1263                 if (isBoundaryPoint.set(f[fp], 1u))
1264                 {
1265                     nBoundaryPoints++;
1266                 }
1267             }
1268         }
1269     }
1270
1271     //// Insert all meshed patches.
1272     //labelList adaptPatchIDs(meshedPatches());
1273     //forAll(adaptPatchIDs, i)
1274     //{
1275     //    label patchI = adaptPatchIDs[i];
1276     //
1277     //    if (patchI != -1)
1278     //    {
1279     //        const polyPatch& pp = mesh_.boundaryMesh()[patchI];
1280     //
1281     //        label faceI = pp.start();
1282     //
1283     //        forAll(pp, i)
1284     //        {
1285     //            const face& f = faces[faceI];
1286     //
1287     //            forAll(f, fp)
1288     //            {
1289     //                if (isBoundaryPoint.set(f[fp], 1u))
1290     //                    nBoundaryPoints++;
1291     //                }
1292     //            }
1293     //            faceI++;
1294     //        }
1295     //    }
1296     //}
1297
1298
1299     // Pack
1300     labelList boundaryPoints(nBoundaryPoints);
1301     nBoundaryPoints = 0;
1302     forAll(isBoundaryPoint, pointI)
1303     {
1304         if (isBoundaryPoint.get(pointI) == 1u)
1305         {
1306             boundaryPoints[nBoundaryPoints++] = pointI;
1307         }
1308     }
1309
1310     return boundaryPoints;
1311 }
1312
1313
1314 //- Create patch from set of patches
1315 Foam::autoPtr<Foam::indirectPrimitivePatch> Foam::meshRefinement::makePatch
1316 (
1317     const polyMesh& mesh,
1318     const labelList& patchIDs
1319 )
1320 {
1321     const polyBoundaryMesh& patches = mesh.boundaryMesh();
1322
1323     // Count faces.
1324     label nFaces = 0;
1325
1326     forAll(patchIDs, i)
1327     {
1328         const polyPatch& pp = patches[patchIDs[i]];
1329
1330         nFaces += pp.size();
1331     }
1332
1333     // Collect faces.
1334     labelList addressing(nFaces);
1335     nFaces = 0;
1336
1337     forAll(patchIDs, i)
1338     {
1339         const polyPatch& pp = patches[patchIDs[i]];
1340
1341         label meshFaceI = pp.start();
1342
1343         forAll(pp, i)
1344         {
1345             addressing[nFaces++] = meshFaceI++;
1346         }
1347     }
1348
1349     return autoPtr<indirectPrimitivePatch>
1350     (
1351         new indirectPrimitivePatch
1352         (
1353             IndirectList<face>(mesh.faces(), addressing),
1354             mesh.points()
1355         )
1356     );
1357 }
1358
1359
1360 // Construct pointVectorField with correct boundary conditions
1361 Foam::tmp<Foam::pointVectorField> Foam::meshRefinement::makeDisplacementField
1362 (
1363     const pointMesh& pMesh,
1364     const labelList& adaptPatchIDs
1365 )
1366 {
1367     const polyMesh& mesh = pMesh();
1368
1369     // Construct displacement field.
1370     const pointBoundaryMesh& pointPatches = pMesh.boundary();
1371
1372     wordList patchFieldTypes
1373     (
1374         pointPatches.size(),
1375         slipPointPatchVectorField::typeName
1376     );
1377
1378     forAll(adaptPatchIDs, i)
1379     {
1380         patchFieldTypes[adaptPatchIDs[i]] =
1381             fixedValuePointPatchVectorField::typeName;
1382     }
1383
1384     forAll(pointPatches, patchI)
1385     {
1386         if (isA<globalPointPatch>(pointPatches[patchI]))
1387         {
1388             patchFieldTypes[patchI] = globalPointPatchVectorField::typeName;
1389         }
1390         else if (isA<processorPointPatch>(pointPatches[patchI]))
1391         {
1392             patchFieldTypes[patchI] = calculatedPointPatchVectorField::typeName;
1393         }
1394     }
1395
1396     tmp<pointVectorField> tfld
1397     (
1398         new pointVectorField
1399         (
1400             IOobject
1401             (
1402                 "pointDisplacement",
1403                 mesh.time().timeName(), //timeName(),
1404                 mesh,
1405                 IOobject::NO_READ,
1406                 IOobject::AUTO_WRITE
1407             ),
1408             pMesh,
1409             dimensionedVector("displacement", dimLength, vector::zero),
1410             patchFieldTypes
1411         )
1412     );
1413     return tfld;
1414 }
1415
1416
1417 void Foam::meshRefinement::checkCoupledFaceZones(const polyMesh& mesh)
1418 {
1419     const faceZoneMesh& fZones = mesh.faceZones();
1420
1421     // Check any zones are present anywhere and in same order
1422
1423     {
1424         List<wordList> zoneNames(Pstream::nProcs());
1425         zoneNames[Pstream::myProcNo()] = fZones.names();
1426         Pstream::gatherList(zoneNames);
1427         Pstream::scatterList(zoneNames);
1428         // All have same data now. Check.
1429         forAll(zoneNames, procI)
1430         {
1431             if (procI != Pstream::myProcNo())
1432             {
1433                 if (zoneNames[procI] != zoneNames[Pstream::myProcNo()])
1434                 {
1435                     FatalErrorIn
1436                     (
1437                         "meshRefinement::checkCoupledFaceZones(const polyMesh&)"
1438                     )   << "faceZones are not synchronised on processors." << nl
1439                         << "Processor " << procI << " has faceZones "
1440                         << zoneNames[procI] << nl
1441                         << "Processor " << Pstream::myProcNo()
1442                         << " has faceZones "
1443                         << zoneNames[Pstream::myProcNo()] << nl
1444                         << exit(FatalError);
1445                 }
1446             }
1447         }
1448     }
1449
1450     // Check that coupled faces are present on both sides.
1451
1452     labelList faceToZone(mesh.nFaces()-mesh.nInternalFaces(), -1);
1453
1454     forAll(fZones, zoneI)
1455     {
1456         const faceZone& fZone = fZones[zoneI];
1457
1458         forAll(fZone, i)
1459         {
1460             label bFaceI = fZone[i]-mesh.nInternalFaces();
1461
1462             if (bFaceI >= 0)
1463             {
1464                 if (faceToZone[bFaceI] == -1)
1465                 {
1466                     faceToZone[bFaceI] = zoneI;
1467                 }
1468                 else if (faceToZone[bFaceI] == zoneI)
1469                 {
1470                     FatalErrorIn
1471                     (
1472                         "meshRefinement::checkCoupledFaceZones(const polyMesh&)"
1473                     )   << "Face " << fZone[i] << " in zone "
1474                         << fZone.name()
1475                         << " is twice in zone!"
1476                         << abort(FatalError);
1477                 }
1478                 else
1479                 {
1480                     FatalErrorIn
1481                     (
1482                         "meshRefinement::checkCoupledFaceZones(const polyMesh&)"
1483                     )   << "Face " << fZone[i] << " in zone "
1484                         << fZone.name()
1485                         << " is also in zone "
1486                         << fZones[faceToZone[bFaceI]].name()
1487                         << abort(FatalError);
1488                 }
1489             }
1490         }
1491     }
1492
1493     labelList neiFaceToZone(faceToZone);
1494     syncTools::swapBoundaryFaceList(mesh, neiFaceToZone, false);
1495
1496     forAll(faceToZone, i)
1497     {
1498         if (faceToZone[i] != neiFaceToZone[i])
1499         {
1500             FatalErrorIn
1501             (
1502                 "meshRefinement::checkCoupledFaceZones(const polyMesh&)"
1503             )   << "Face " << mesh.nInternalFaces()+i
1504                 << " is in zone " << faceToZone[i]
1505                 << ", its coupled face is in zone " << neiFaceToZone[i]
1506                 << abort(FatalError);
1507         }
1508     }
1509 }
1510
1511
1512 // Adds patch if not yet there. Returns patchID.
1513 Foam::label Foam::meshRefinement::addPatch
1514 (
1515     fvMesh& mesh,
1516     const word& patchName,
1517     const word& patchType
1518 )
1519 {
1520     polyBoundaryMesh& polyPatches =
1521         const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
1522
1523     label patchI = polyPatches.findPatchID(patchName);
1524     if (patchI != -1)
1525     {
1526         if (polyPatches[patchI].type() == patchType)
1527         {
1528             // Already there
1529             return patchI;
1530         }
1531         //else
1532         //{
1533         //    FatalErrorIn
1534         //    (
1535         //        "meshRefinement::addPatch(fvMesh&, const word&, const word&)"
1536         //    )   << "Patch " << patchName << " already exists but with type "
1537         //        << patchType << nl
1538         //        << "Current patch names:" << polyPatches.names()
1539         //        << exit(FatalError);
1540         //}
1541     }
1542
1543
1544     label insertPatchI = polyPatches.size();
1545     label startFaceI = mesh.nFaces();
1546
1547     forAll(polyPatches, patchI)
1548     {
1549         const polyPatch& pp = polyPatches[patchI];
1550
1551         if (isA<processorPolyPatch>(pp))
1552         {
1553             insertPatchI = patchI;
1554             startFaceI = pp.start();
1555             break;
1556         }
1557     }
1558
1559
1560     // Below is all quite a hack. Feel free to change once there is a better
1561     // mechanism to insert and reorder patches.
1562
1563     // Clear local fields and e.g. polyMesh parallelInfo.
1564     mesh.clearOut();
1565
1566     label sz = polyPatches.size();
1567
1568     fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
1569
1570     // Add polyPatch at the end
1571     polyPatches.setSize(sz+1);
1572     polyPatches.set
1573     (
1574         sz,
1575         polyPatch::New
1576         (
1577             patchType,
1578             patchName,
1579             0,              // size
1580             startFaceI,
1581             insertPatchI,
1582             polyPatches
1583         )
1584     );
1585     fvPatches.setSize(sz+1);
1586     fvPatches.set
1587     (
1588         sz,
1589         fvPatch::New
1590         (
1591             polyPatches[sz],  // point to newly added polyPatch
1592             mesh.boundary()
1593         )
1594     );
1595
1596     addPatchFields<volScalarField>
1597     (
1598         mesh,
1599         calculatedFvPatchField<scalar>::typeName
1600     );
1601     addPatchFields<volVectorField>
1602     (
1603         mesh,
1604         calculatedFvPatchField<vector>::typeName
1605     );
1606     addPatchFields<volSphericalTensorField>
1607     (
1608         mesh,
1609         calculatedFvPatchField<sphericalTensor>::typeName
1610     );
1611     addPatchFields<volSymmTensorField>
1612     (
1613         mesh,
1614         calculatedFvPatchField<symmTensor>::typeName
1615     );
1616     addPatchFields<volTensorField>
1617     (
1618         mesh,
1619         calculatedFvPatchField<tensor>::typeName
1620     );
1621
1622     // Surface fields
1623
1624     addPatchFields<surfaceScalarField>
1625     (
1626         mesh,
1627         calculatedFvPatchField<scalar>::typeName
1628     );
1629     addPatchFields<surfaceVectorField>
1630     (
1631         mesh,
1632         calculatedFvPatchField<vector>::typeName
1633     );
1634     addPatchFields<surfaceSphericalTensorField>
1635     (
1636         mesh,
1637         calculatedFvPatchField<sphericalTensor>::typeName
1638     );
1639     addPatchFields<surfaceSymmTensorField>
1640     (
1641         mesh,
1642         calculatedFvPatchField<symmTensor>::typeName
1643     );
1644     addPatchFields<surfaceTensorField>
1645     (
1646         mesh,
1647         calculatedFvPatchField<tensor>::typeName
1648     );
1649
1650     // Create reordering list
1651     // patches before insert position stay as is
1652     labelList oldToNew(sz+1);
1653     for (label i = 0; i < insertPatchI; i++)
1654     {
1655         oldToNew[i] = i;
1656     }
1657     // patches after insert position move one up
1658     for (label i = insertPatchI; i < sz; i++)
1659     {
1660         oldToNew[i] = i+1;
1661     }
1662     // appended patch gets moved to insert position
1663     oldToNew[sz] = insertPatchI;
1664
1665     // Shuffle into place
1666     polyPatches.reorder(oldToNew);
1667     fvPatches.reorder(oldToNew);
1668
1669     reorderPatchFields<volScalarField>(mesh, oldToNew);
1670     reorderPatchFields<volVectorField>(mesh, oldToNew);
1671     reorderPatchFields<volSphericalTensorField>(mesh, oldToNew);
1672     reorderPatchFields<volSymmTensorField>(mesh, oldToNew);
1673     reorderPatchFields<volTensorField>(mesh, oldToNew);
1674     reorderPatchFields<surfaceScalarField>(mesh, oldToNew);
1675     reorderPatchFields<surfaceVectorField>(mesh, oldToNew);
1676     reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
1677     reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
1678     reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
1679
1680     return insertPatchI;
1681 }
1682
1683
1684 Foam::label Foam::meshRefinement::addMeshedPatch
1685 (
1686     const word& name,
1687     const word& type
1688 )
1689 {
1690     label meshedI = findIndex(meshedPatches_, name);
1691
1692     if (meshedI != -1)
1693     {
1694         // Already there. Get corresponding polypatch
1695         return mesh_.boundaryMesh().findPatchID(name);
1696     }
1697     else
1698     {
1699         // Add patch
1700         label patchI = addPatch(mesh_, name, type);
1701
1702         // Store
1703         label sz = meshedPatches_.size();
1704         meshedPatches_.setSize(sz+1);
1705         meshedPatches_[sz] = name;
1706
1707         return patchI;
1708     }
1709 }
1710
1711
1712 Foam::labelList Foam::meshRefinement::meshedPatches() const
1713 {
1714     labelList patchIDs(meshedPatches_.size());
1715     forAll(meshedPatches_, i)
1716     {
1717         patchIDs[i] = mesh_.boundaryMesh().findPatchID(meshedPatches_[i]);
1718
1719         if (patchIDs[i] == -1)
1720         {
1721             FatalErrorIn("meshRefinement::meshedPatches() const")
1722                 << "Problem : did not find patch " << meshedPatches_[i]
1723                 << abort(FatalError);
1724         }
1725     }
1726
1727     return patchIDs;
1728 }
1729
1730
1731 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMeshRegions
1732 (
1733     const point& keepPoint
1734 )
1735 {
1736     // Determine connected regions. regionSplit is the labelList with the
1737     // region per cell.
1738     regionSplit cellRegion(mesh_);
1739
1740     label regionI = -1;
1741
1742     label cellI = mesh_.findCell(keepPoint);
1743
1744     if (cellI != -1)
1745     {
1746         regionI = cellRegion[cellI];
1747     }
1748
1749     reduce(regionI, maxOp<label>());
1750
1751     if (regionI == -1)
1752     {
1753         FatalErrorIn
1754         (
1755             "meshRefinement::splitMeshRegions(const point&)"
1756         )   << "Point " << keepPoint
1757             << " is not inside the mesh." << nl
1758             << "Bounding box of the mesh:" << mesh_.globalData().bb()
1759             << exit(FatalError);
1760     }
1761
1762     // Subset
1763     // ~~~~~~
1764
1765     // Get cells to remove
1766     DynamicList<label> cellsToRemove(mesh_.nCells());
1767     forAll(cellRegion, cellI)
1768     {
1769         if (cellRegion[cellI] != regionI)
1770         {
1771             cellsToRemove.append(cellI);
1772         }
1773     }
1774     cellsToRemove.shrink();
1775
1776     label nCellsToKeep = mesh_.nCells() - cellsToRemove.size();
1777     reduce(nCellsToKeep, sumOp<label>());
1778
1779     Info<< "Keeping all cells in region " << regionI
1780         << " containing point " << keepPoint << endl
1781         << "Selected for keeping : "
1782         << nCellsToKeep
1783         << " cells." << endl;
1784
1785
1786     // Remove cells
1787     removeCells cellRemover(mesh_);
1788
1789     labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
1790
1791     if (exposedFaces.size())
1792     {
1793         FatalErrorIn
1794         (
1795             "meshRefinement::splitMeshRegions(const point&)"
1796         )   << "Removing non-reachable cells should only expose boundary faces"
1797             << nl
1798             << "ExposedFaces:" << exposedFaces << abort(FatalError);
1799     }
1800
1801     return doRemoveCells
1802     (
1803         cellsToRemove,
1804         exposedFaces,
1805         labelList(exposedFaces.size(),-1),  // irrelevant since 0 size.
1806         cellRemover
1807     );
1808 }
1809
1810
1811 void Foam::meshRefinement::distribute(const mapDistributePolyMesh& map)
1812 {
1813     // mesh_ already distributed; distribute my member data
1814
1815     // surfaceQueries_ ok.
1816
1817     // refinement
1818     meshCutter_.distribute(map);
1819
1820     // surfaceIndex is face data.
1821     map.distributeFaceData(surfaceIndex_);
1822
1823     // maintainedFaces are indices of faces.
1824     forAll(userFaceData_, i)
1825     {
1826         map.distributeFaceData(userFaceData_[i].second());
1827     }
1828
1829     // Redistribute surface and any fields on it.
1830     {
1831         Random rndGen(653213);
1832
1833         // Get local mesh bounding box. Single box for now.
1834         List<treeBoundBox> meshBb(1);
1835         treeBoundBox& bb = meshBb[0];
1836         bb = treeBoundBox(mesh_.points());
1837         bb = bb.extend(rndGen, 1E-4);
1838
1839         // Distribute all geometry (so refinementSurfaces and shellSurfaces)
1840         searchableSurfaces& geometry =
1841             const_cast<searchableSurfaces&>(surfaces_.geometry());
1842
1843         forAll(geometry, i)
1844         {
1845             autoPtr<mapDistribute> faceMap;
1846             autoPtr<mapDistribute> pointMap;
1847             geometry[i].distribute
1848             (
1849                 meshBb,
1850                 false,          // do not keep outside triangles
1851                 faceMap,
1852                 pointMap
1853             );
1854
1855             if (faceMap.valid())
1856             {
1857                 // (ab)use the instance() to signal current modification time
1858                 geometry[i].instance() = geometry[i].time().timeName();
1859             }
1860
1861             faceMap.clear();
1862             pointMap.clear();
1863         }
1864     }
1865 }
1866
1867
1868 void Foam::meshRefinement::updateMesh
1869 (
1870     const mapPolyMesh& map,
1871     const labelList& changedFaces
1872 )
1873 {
1874     Map<label> dummyMap(0);
1875
1876     updateMesh(map, changedFaces, dummyMap, dummyMap, dummyMap);
1877 }
1878
1879
1880 void Foam::meshRefinement::storeData
1881 (
1882     const labelList& pointsToStore,
1883     const labelList& facesToStore,
1884     const labelList& cellsToStore
1885 )
1886 {
1887     // For now only meshCutter has storable/retrievable data.
1888     meshCutter_.storeData
1889     (
1890         pointsToStore,
1891         facesToStore,
1892         cellsToStore
1893     );
1894 }
1895
1896
1897 void Foam::meshRefinement::updateMesh
1898 (
1899     const mapPolyMesh& map,
1900     const labelList& changedFaces,
1901     const Map<label>& pointsToRestore,
1902     const Map<label>& facesToRestore,
1903     const Map<label>& cellsToRestore
1904 )
1905 {
1906     // For now only meshCutter has storable/retrievable data.
1907
1908     // Update numbering of cells/vertices.
1909     meshCutter_.updateMesh
1910     (
1911         map,
1912         pointsToRestore,
1913         facesToRestore,
1914         cellsToRestore
1915     );
1916
1917     // Update surfaceIndex
1918     updateList(map.faceMap(), -1, surfaceIndex_);
1919
1920     // Update cached intersection information
1921     updateIntersections(changedFaces);
1922
1923     // Update maintained faces
1924     forAll(userFaceData_, i)
1925     {
1926         labelList& data = userFaceData_[i].second();
1927
1928         if (userFaceData_[i].first() == KEEPALL)
1929         {
1930             // extend list with face-from-face data
1931             updateList(map.faceMap(), -1, data);
1932         }
1933         else if (userFaceData_[i].first() == MASTERONLY)
1934         {
1935             // keep master only
1936             labelList newFaceData(map.faceMap().size(), -1);
1937
1938             forAll(newFaceData, faceI)
1939             {
1940                 label oldFaceI = map.faceMap()[faceI];
1941
1942                 if (oldFaceI >= 0 && map.reverseFaceMap()[oldFaceI] == faceI)
1943                 {
1944                     newFaceData[faceI] = data[oldFaceI];
1945                 }
1946             }
1947             data.transfer(newFaceData);
1948         }
1949         else
1950         {
1951             // remove any face that has been refined i.e. referenced more than
1952             // once.
1953
1954             // 1. Determine all old faces that get referenced more than once.
1955             // These get marked with -1 in reverseFaceMap
1956             labelList reverseFaceMap(map.reverseFaceMap());
1957
1958             forAll(map.faceMap(), faceI)
1959             {
1960                 label oldFaceI = map.faceMap()[faceI];
1961
1962                 if (oldFaceI >= 0)
1963                 {
1964                     if (reverseFaceMap[oldFaceI] != faceI)
1965                     {
1966                         // faceI is slave face. Mark old face.
1967                         reverseFaceMap[oldFaceI] = -1;
1968                     }
1969                 }
1970             }
1971
1972             // 2. Map only faces with intact reverseFaceMap
1973             labelList newFaceData(map.faceMap().size(), -1);
1974             forAll(newFaceData, faceI)
1975             {
1976                 label oldFaceI = map.faceMap()[faceI];
1977
1978                 if (oldFaceI >= 0)
1979                 {
1980                     if (reverseFaceMap[oldFaceI] == faceI)
1981                     {
1982                         newFaceData[faceI] = data[oldFaceI];
1983                     }
1984                 }
1985             }
1986             data.transfer(newFaceData);
1987         }
1988     }
1989 }
1990
1991
1992 bool Foam::meshRefinement::write() const
1993 {
1994     bool writeOk =
1995         mesh_.write()
1996      && meshCutter_.write()
1997      && surfaceIndex_.write();
1998
1999
2000     // Make sure that any distributed surfaces (so ones which probably have
2001     // been changed) get written as well.
2002     // Note: should ideally have some 'modified' flag to say whether it
2003     // has been changed or not.
2004     searchableSurfaces& geometry =
2005         const_cast<searchableSurfaces&>(surfaces_.geometry());
2006
2007     forAll(geometry, i)
2008     {
2009         searchableSurface& s = geometry[i];
2010
2011         // Check if instance() of surface is not constant or system.
2012         // Is good hint that surface is distributed.
2013         if
2014         (
2015             s.instance() != s.time().system()
2016          && s.instance() != s.time().caseSystem()
2017          && s.instance() != s.time().constant()
2018          && s.instance() != s.time().caseConstant()
2019         )
2020         {
2021             // Make sure it gets written to current time, not constant.
2022             s.instance() = s.time().timeName();
2023             writeOk = writeOk && s.write();
2024         }
2025     }
2026
2027     return writeOk;
2028 }
2029
2030
2031 void Foam::meshRefinement::printMeshInfo(const bool debug, const string& msg)
2032  const
2033 {
2034     const globalMeshData& pData = mesh_.globalData();
2035
2036     if (debug)
2037     {
2038         Pout<< msg.c_str()
2039             << " : cells(local):" << mesh_.nCells()
2040             << "  faces(local):" << mesh_.nFaces()
2041             << "  points(local):" << mesh_.nPoints()
2042             << endl;
2043     }
2044     else
2045     {
2046         Info<< msg.c_str()
2047             << " : cells:" << pData.nTotalCells()
2048             << "  faces:" << pData.nTotalFaces()
2049             << "  points:" << pData.nTotalPoints()
2050             << endl;
2051     }
2052
2053
2054     //if (debug)
2055     {
2056         const labelList& cellLevel = meshCutter_.cellLevel();
2057
2058         labelList nCells(gMax(cellLevel)+1, 0);
2059
2060         forAll(cellLevel, cellI)
2061         {
2062             nCells[cellLevel[cellI]]++;
2063         }
2064
2065         Pstream::listCombineGather(nCells, plusEqOp<label>());
2066         Pstream::listCombineScatter(nCells);
2067
2068         Info<< "Cells per refinement level:" << endl;
2069         forAll(nCells, levelI)
2070         {
2071             Info<< "    " << levelI << '\t' << nCells[levelI]
2072                 << endl;
2073         }
2074     }
2075 }
2076
2077
2078 //- Return either time().constant() or oldInstance
2079 Foam::word Foam::meshRefinement::timeName() const
2080 {
2081     if (overwrite_ && mesh_.time().timeIndex() == 0)
2082     {
2083         return oldInstance_;
2084     }
2085     else
2086     {
2087         return mesh_.time().timeName();
2088     }
2089 }
2090
2091
2092 void Foam::meshRefinement::dumpRefinementLevel() const
2093 {
2094     volScalarField volRefLevel
2095     (
2096         IOobject
2097         (
2098             "cellLevel",
2099             timeName(),
2100             mesh_,
2101             IOobject::NO_READ,
2102             IOobject::AUTO_WRITE,
2103             false
2104         ),
2105         mesh_,
2106         dimensionedScalar("zero", dimless, 0)
2107     );
2108
2109     const labelList& cellLevel = meshCutter_.cellLevel();
2110
2111     forAll(volRefLevel, cellI)
2112     {
2113         volRefLevel[cellI] = cellLevel[cellI];
2114     }
2115
2116     volRefLevel.write();
2117
2118
2119     const pointMesh& pMesh = pointMesh::New(mesh_);
2120
2121     pointScalarField pointRefLevel
2122     (
2123         IOobject
2124         (
2125             "pointLevel",
2126             timeName(),
2127             mesh_,
2128             IOobject::NO_READ,
2129             IOobject::NO_WRITE,
2130             false
2131         ),
2132         pMesh,
2133         dimensionedScalar("zero", dimless, 0)
2134     );
2135
2136     const labelList& pointLevel = meshCutter_.pointLevel();
2137
2138     forAll(pointRefLevel, pointI)
2139     {
2140         pointRefLevel[pointI] = pointLevel[pointI];
2141     }
2142
2143     pointRefLevel.write();
2144 }
2145
2146
2147 void Foam::meshRefinement::dumpIntersections(const fileName& prefix) const
2148 {
2149     {
2150         const pointField& cellCentres = mesh_.cellCentres();
2151
2152         OFstream str(prefix + "_edges.obj");
2153         label vertI = 0;
2154         Pout<< "meshRefinement::dumpIntersections :"
2155             << " Writing cellcentre-cellcentre intersections to file "
2156             << str.name() << endl;
2157
2158
2159         // Redo all intersections
2160         // ~~~~~~~~~~~~~~~~~~~~~~
2161
2162         // Get boundary face centre and level. Coupled aware.
2163         labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
2164         pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
2165         calcNeighbourData(neiLevel, neiCc);
2166
2167         labelList intersectionFaces(intersectedFaces());
2168
2169         // Collect segments we want to test for
2170         pointField start(intersectionFaces.size());
2171         pointField end(intersectionFaces.size());
2172
2173         forAll(intersectionFaces, i)
2174         {
2175             label faceI = intersectionFaces[i];
2176             start[i] = cellCentres[mesh_.faceOwner()[faceI]];
2177
2178             if (mesh_.isInternalFace(faceI))
2179             {
2180                 end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
2181             }
2182             else
2183             {
2184                 end[i] = neiCc[faceI-mesh_.nInternalFaces()];
2185             }
2186         }
2187
2188         // Do tests in one go
2189         labelList surfaceHit;
2190         List<pointIndexHit> surfaceHitInfo;
2191         surfaces_.findAnyIntersection
2192         (
2193             start,
2194             end,
2195             surfaceHit,
2196             surfaceHitInfo
2197         );
2198
2199         forAll(intersectionFaces, i)
2200         {
2201             if (surfaceHit[i] != -1)
2202             {
2203                 meshTools::writeOBJ(str, start[i]);
2204                 vertI++;
2205                 meshTools::writeOBJ(str, surfaceHitInfo[i].hitPoint());
2206                 vertI++;
2207                 meshTools::writeOBJ(str, end[i]);
2208                 vertI++;
2209                 str << "l " << vertI-2 << ' ' << vertI-1 << nl
2210                     << "l " << vertI-1 << ' ' << vertI << nl;
2211             }
2212         }
2213     }
2214
2215     // Convert to vtk format
2216     string cmd
2217     (
2218         "objToVTK " + prefix + "_edges.obj " + prefix + "_edges.vtk > /dev/null"
2219     );
2220     system(cmd.c_str());
2221
2222     Pout<< endl;
2223 }
2224
2225
2226 void Foam::meshRefinement::write
2227 (
2228     const label flag,
2229     const fileName& prefix
2230 ) const
2231 {
2232     if (flag & MESH)
2233     {
2234         write();
2235     }
2236     if (flag & SCALARLEVELS)
2237     {
2238         dumpRefinementLevel();
2239     }
2240     if (flag & OBJINTERSECTIONS && prefix.size())
2241     {
2242         dumpIntersections(prefix);
2243     }
2244 }
2245
2246
2247 // ************************************************************************* //