applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 Application
  25     splitMeshRegions
  26
  27 Description
  28     Splits mesh into multiple regions.
  29
  30     Each region is defined as a domain whose cells can all be reached by
  31     cell-face-cell walking without crossing
  32     - boundary faces
  33     - additional faces from faceset (-blockedFaces faceSet).
  34     - any face inbetween differing cellZones (-cellZones)
  35
  36     Output is:
  37     - volScalarField with regions as different scalars (-detectOnly)
  38             or
  39     - mesh with multiple regions and directMapped patches. These patches
  40       either cover the whole interface between two region (default) or
  41       only part according to faceZones (-useFaceZones)
  42             or
  43     - mesh with cells put into cellZones (-makeCellZones)
  44
  45     Note:
  46     - cellZonesOnly does not do a walk and uses the cellZones only. Use
  47     this if you don't mind having disconnected domains in a single region.
  48     This option requires all cells to be in one (and one only) cellZone.
  49
  50     - cellZonesFileOnly behaves like -cellZonesOnly but reads the cellZones
  51     from the specified file. This allows one to explicitly specify the region
  52     distribution and still have multiple cellZones per region.
  53
  54     - useCellZonesOnly does not do a walk and uses the cellZones only. Use
  55     this if you don't mind having disconnected domains in a single region.
  56     This option requires all cells to be in one (and one only) cellZone.
  57
  58
  59     - Should work in parallel.
  60     cellZones can differ on either side of processor boundaries in which case
  61     the faces get moved from processor patch to directMapped patch. Not
  62     very well tested.
  63
  64     - If a cell zone gets split into more than one region it can detect
  65     the largest matching region (-sloppyCellZones). This will accept any
  66     region that covers more than 50% of the zone. It has to be a subset
  67     so cannot have any cells in any other zone.
  68
  69     - writes maps like decomposePar back to original mesh:
  70         - pointRegionAddressing : for every point in this region the point in
  71         the original mesh
  72         - cellRegionAddressing  :   ,,      cell                ,,  cell    ,,
  73         - faceRegionAddressing  :   ,,      face                ,,  face in
  74         the original mesh + 'turning index'. For a face in the same orientation
  75         this is the original facelabel+1, for a turned face this is -facelabel-1
  76         - boundaryRegionAddressing : for every patch in this region the
  77         patch in the original mesh (or -1 if added patch)
  78 \*---------------------------------------------------------------------------*/
  79
  80 #include "SortableList.H"
  81 #include "argList.H"
  82 #include "regionSplit.H"
  83 #include "fvMeshSubset.H"
  84 #include "IOobjectList.H"
  85 #include "volFields.H"
  86 #include "faceSet.H"
  87 #include "cellSet.H"
  88 #include "polyTopoChange.H"
  89 #include "removeCells.H"
  90 #include "EdgeMap.H"
  91 #include "syncTools.H"
  92 #include "ReadFields.H"
  93 #include "directMappedWallPolyPatch.H"
  94 #include "zeroGradientFvPatchFields.H"
  95
  96 using namespace Foam;
  97
  98 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  99
 100 template<class GeoField>
 101 void addPatchFields(fvMesh& mesh, const word& patchFieldType)
 102 {
 103     HashTable<const GeoField*> flds
 104     (
 105         mesh.objectRegistry::lookupClass<GeoField>()
 106     );
 107
 108     forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
 109     {
 110         const GeoField& fld = *iter();
 111
 112         typename GeoField::GeometricBoundaryField& bfld =
 113             const_cast<typename GeoField::GeometricBoundaryField&>
 114             (
 115                 fld.boundaryField()
 116             );
 117
 118         label sz = bfld.size();
 119         bfld.setSize(sz+1);
 120         bfld.set
 121         (
 122             sz,
 123             GeoField::PatchFieldType::New
 124             (
 125                 patchFieldType,
 126                 mesh.boundary()[sz],
 127                 fld.dimensionedInternalField()
 128             )
 129         );
 130     }
 131 }
 132
 133
 134 // Remove last patch field
 135 template<class GeoField>
 136 void trimPatchFields(fvMesh& mesh, const label nPatches)
 137 {
 138     HashTable<const GeoField*> flds
 139     (
 140         mesh.objectRegistry::lookupClass<GeoField>()
 141     );
 142
 143     forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
 144     {
 145         const GeoField& fld = *iter();
 146
 147         const_cast<typename GeoField::GeometricBoundaryField&>
 148         (
 149             fld.boundaryField()
 150         ).setSize(nPatches);
 151     }
 152 }
 153
 154
 155 // Reorder patch field
 156 template<class GeoField>
 157 void reorderPatchFields(fvMesh& mesh, const labelList& oldToNew)
 158 {
 159     HashTable<const GeoField*> flds
 160     (
 161         mesh.objectRegistry::lookupClass<GeoField>()
 162     );
 163
 164     forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
 165     {
 166         const GeoField& fld = *iter();
 167
 168         typename GeoField::GeometricBoundaryField& bfld =
 169             const_cast<typename GeoField::GeometricBoundaryField&>
 170             (
 171                 fld.boundaryField()
 172             );
 173
 174         bfld.reorder(oldToNew);
 175     }
 176 }
 177
 178
 179 // Adds patch if not yet there. Returns patchID.
 180 label addPatch(fvMesh& mesh, const polyPatch& patch)
 181 {
 182     polyBoundaryMesh& polyPatches =
 183         const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
 184
 185     label patchI = polyPatches.findPatchID(patch.name());
 186     if (patchI != -1)
 187     {
 188         if (polyPatches[patchI].type() == patch.type())
 189         {
 190             // Already there
 191             return patchI;
 192         }
 193         else
 194         {
 195             FatalErrorIn("addPatch(fvMesh&, const polyPatch*)")
 196                 << "Already have patch " << patch.name()
 197                 << " but of type " << patch.type()
 198                 << exit(FatalError);
 199         }
 200     }
 201
 202
 203     label insertPatchI = polyPatches.size();
 204     label startFaceI = mesh.nFaces();
 205
 206     forAll(polyPatches, patchI)
 207     {
 208         const polyPatch& pp = polyPatches[patchI];
 209
 210         if (isA<processorPolyPatch>(pp))
 211         {
 212             insertPatchI = patchI;
 213             startFaceI = pp.start();
 214             break;
 215         }
 216     }
 217
 218
 219     // Below is all quite a hack. Feel free to change once there is a better
 220     // mechanism to insert and reorder patches.
 221
 222     // Clear local fields and e.g. polyMesh parallelInfo.
 223     mesh.clearOut();
 224
 225     label sz = polyPatches.size();
 226
 227     fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
 228
 229     // Add polyPatch at the end
 230     polyPatches.setSize(sz+1);
 231     polyPatches.set
 232     (
 233         sz,
 234         patch.clone
 235         (
 236             polyPatches,
 237             insertPatchI,   //index
 238             0,              //size
 239             startFaceI      //start
 240         )
 241     );
 242     fvPatches.setSize(sz+1);
 243     fvPatches.set
 244     (
 245         sz,
 246         fvPatch::New
 247         (
 248             polyPatches[sz],  // point to newly added polyPatch
 249             mesh.boundary()
 250         )
 251     );
 252
 253     addPatchFields<volScalarField>
 254     (
 255         mesh,
 256         calculatedFvPatchField<scalar>::typeName
 257     );
 258     addPatchFields<volVectorField>
 259     (
 260         mesh,
 261         calculatedFvPatchField<vector>::typeName
 262     );
 263     addPatchFields<volSphericalTensorField>
 264     (
 265         mesh,
 266         calculatedFvPatchField<sphericalTensor>::typeName
 267     );
 268     addPatchFields<volSymmTensorField>
 269     (
 270         mesh,
 271         calculatedFvPatchField<symmTensor>::typeName
 272     );
 273     addPatchFields<volTensorField>
 274     (
 275         mesh,
 276         calculatedFvPatchField<tensor>::typeName
 277     );
 278
 279     // Surface fields
 280
 281     addPatchFields<surfaceScalarField>
 282     (
 283         mesh,
 284         calculatedFvPatchField<scalar>::typeName
 285     );
 286     addPatchFields<surfaceVectorField>
 287     (
 288         mesh,
 289         calculatedFvPatchField<vector>::typeName
 290     );
 291     addPatchFields<surfaceSphericalTensorField>
 292     (
 293         mesh,
 294         calculatedFvPatchField<sphericalTensor>::typeName
 295     );
 296     addPatchFields<surfaceSymmTensorField>
 297     (
 298         mesh,
 299         calculatedFvPatchField<symmTensor>::typeName
 300     );
 301     addPatchFields<surfaceTensorField>
 302     (
 303         mesh,
 304         calculatedFvPatchField<tensor>::typeName
 305     );
 306
 307     // Create reordering list
 308     // patches before insert position stay as is
 309     labelList oldToNew(sz+1);
 310     for (label i = 0; i < insertPatchI; i++)
 311     {
 312         oldToNew[i] = i;
 313     }
 314     // patches after insert position move one up
 315     for (label i = insertPatchI; i < sz; i++)
 316     {
 317         oldToNew[i] = i+1;
 318     }
 319     // appended patch gets moved to insert position
 320     oldToNew[sz] = insertPatchI;
 321
 322     // Shuffle into place
 323     polyPatches.reorder(oldToNew);
 324     fvPatches.reorder(oldToNew);
 325
 326     reorderPatchFields<volScalarField>(mesh, oldToNew);
 327     reorderPatchFields<volVectorField>(mesh, oldToNew);
 328     reorderPatchFields<volSphericalTensorField>(mesh, oldToNew);
 329     reorderPatchFields<volSymmTensorField>(mesh, oldToNew);
 330     reorderPatchFields<volTensorField>(mesh, oldToNew);
 331     reorderPatchFields<surfaceScalarField>(mesh, oldToNew);
 332     reorderPatchFields<surfaceVectorField>(mesh, oldToNew);
 333     reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
 334     reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
 335     reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
 336
 337     return insertPatchI;
 338 }
 339
 340
 341 // Reorder and delete patches.
 342 void reorderPatches
 343 (
 344     fvMesh& mesh,
 345     const labelList& oldToNew,
 346     const label nNewPatches
 347 )
 348 {
 349     polyBoundaryMesh& polyPatches =
 350         const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
 351     fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
 352
 353     // Shuffle into place
 354     polyPatches.reorder(oldToNew);
 355     fvPatches.reorder(oldToNew);
 356
 357     reorderPatchFields<volScalarField>(mesh, oldToNew);
 358     reorderPatchFields<volVectorField>(mesh, oldToNew);
 359     reorderPatchFields<volSphericalTensorField>(mesh, oldToNew);
 360     reorderPatchFields<volSymmTensorField>(mesh, oldToNew);
 361     reorderPatchFields<volTensorField>(mesh, oldToNew);
 362     reorderPatchFields<surfaceScalarField>(mesh, oldToNew);
 363     reorderPatchFields<surfaceVectorField>(mesh, oldToNew);
 364     reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
 365     reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
 366     reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
 367
 368     // Remove last.
 369     polyPatches.setSize(nNewPatches);
 370     fvPatches.setSize(nNewPatches);
 371     trimPatchFields<volScalarField>(mesh, nNewPatches);
 372     trimPatchFields<volVectorField>(mesh, nNewPatches);
 373     trimPatchFields<volSphericalTensorField>(mesh, nNewPatches);
 374     trimPatchFields<volSymmTensorField>(mesh, nNewPatches);
 375     trimPatchFields<volTensorField>(mesh, nNewPatches);
 376     trimPatchFields<surfaceScalarField>(mesh, nNewPatches);
 377     trimPatchFields<surfaceVectorField>(mesh, nNewPatches);
 378     trimPatchFields<surfaceSphericalTensorField>(mesh, nNewPatches);
 379     trimPatchFields<surfaceSymmTensorField>(mesh, nNewPatches);
 380     trimPatchFields<surfaceTensorField>(mesh, nNewPatches);
 381 }
 382
 383
 384 template<class GeoField>
 385 void subsetVolFields
 386 (
 387     const fvMesh& mesh,
 388     const fvMesh& subMesh,
 389     const labelList& cellMap,
 390     const labelList& faceMap,
 391     const labelHashSet& addedPatches
 392 )
 393 {
 394     const labelList patchMap(identity(mesh.boundaryMesh().size()));
 395
 396     HashTable<const GeoField*> fields
 397     (
 398         mesh.objectRegistry::lookupClass<GeoField>()
 399     );
 400     forAllConstIter(typename HashTable<const GeoField*>, fields, iter)
 401     {
 402         const GeoField& fld = *iter();
 403
 404         Info<< "Mapping field " << fld.name() << endl;
 405
 406         tmp<GeoField> tSubFld
 407         (
 408             fvMeshSubset::interpolate
 409             (
 410                 fld,
 411                 subMesh,
 412                 patchMap,
 413                 cellMap,
 414                 faceMap
 415             )
 416         );
 417
 418         // Hack: set value to 0 for introduced patches (since don't
 419         //       get initialised.
 420         forAll(tSubFld().boundaryField(), patchI)
 421         {
 422             if (addedPatches.found(patchI))
 423             {
 424                 tSubFld().boundaryField()[patchI] ==
 425                     pTraits<typename GeoField::value_type>::zero;
 426             }
 427         }
 428
 429         // Store on subMesh
 430         GeoField* subFld = tSubFld.ptr();
 431         subFld->rename(fld.name());
 432         subFld->writeOpt() = IOobject::AUTO_WRITE;
 433         subFld->store();
 434     }
 435 }
 436
 437
 438 template<class GeoField>
 439 void subsetSurfaceFields
 440 (
 441     const fvMesh& mesh,
 442     const fvMesh& subMesh,
 443     const labelList& faceMap,
 444     const labelHashSet& addedPatches
 445 )
 446 {
 447     const labelList patchMap(identity(mesh.boundaryMesh().size()));
 448
 449     HashTable<const GeoField*> fields
 450     (
 451         mesh.objectRegistry::lookupClass<GeoField>()
 452     );
 453     forAllConstIter(typename HashTable<const GeoField*>, fields, iter)
 454     {
 455         const GeoField& fld = *iter();
 456
 457         Info<< "Mapping field " << fld.name() << endl;
 458
 459         tmp<GeoField> tSubFld
 460         (
 461             fvMeshSubset::interpolate
 462             (
 463                 fld,
 464                 subMesh,
 465                 patchMap,
 466                 faceMap
 467             )
 468         );
 469
 470         // Hack: set value to 0 for introduced patches (since don't
 471         //       get initialised.
 472         forAll(tSubFld().boundaryField(), patchI)
 473         {
 474             if (addedPatches.found(patchI))
 475             {
 476                 tSubFld().boundaryField()[patchI] ==
 477                     pTraits<typename GeoField::value_type>::zero;
 478             }
 479         }
 480
 481         // Store on subMesh
 482         GeoField* subFld = tSubFld.ptr();
 483         subFld->rename(fld.name());
 484         subFld->writeOpt() = IOobject::AUTO_WRITE;
 485         subFld->store();
 486     }
 487 }
 488
 489 // Select all cells not in the region
 490 labelList getNonRegionCells(const labelList& cellRegion, const label regionI)
 491 {
 492     DynamicList<label> nonRegionCells(cellRegion.size());
 493     forAll(cellRegion, cellI)
 494     {
 495         if (cellRegion[cellI] != regionI)
 496         {
 497             nonRegionCells.append(cellI);
 498         }
 499     }
 500     return nonRegionCells.shrink();
 501 }
 502
 503
 504 void addToInterface
 505 (
 506     const polyMesh& mesh,
 507     const label zoneID,
 508     const label ownRegion,
 509     const label neiRegion,
 510     EdgeMap<Map<label> >& regionsToSize
 511 )
 512 {
 513     edge interface
 514     (
 515         min(ownRegion, neiRegion),
 516         max(ownRegion, neiRegion)
 517     );
 518
 519     EdgeMap<Map<label> >::iterator iter = regionsToSize.find
 520     (
 521         interface
 522     );
 523
 524     if (iter != regionsToSize.end())
 525     {
 526         // Check if zone present
 527         Map<label>::iterator zoneFnd = iter().find(zoneID);
 528         if (zoneFnd != iter().end())
 529         {
 530             // Found zone. Increment count.
 531             zoneFnd()++;
 532         }
 533         else
 534         {
 535             // New or no zone. Insert with count 1.
 536             iter().insert(zoneID, 1);
 537         }
 538     }
 539     else
 540     {
 541         // Create new interface of size 1.
 542         Map<label> zoneToSize;
 543         zoneToSize.insert(zoneID, 1);
 544         regionsToSize.insert(interface, zoneToSize);
 545     }
 546 }
 547
 548
 549 // Get region-region interface name and sizes.
 550 // Returns interfaces as straight list for looping in identical order.
 551 void getInterfaceSizes
 552 (
 553     const polyMesh& mesh,
 554     const bool useFaceZones,
 555     const labelList& cellRegion,
 556     const wordList& regionNames,
 557
 558     edgeList& interfaces,
 559     List<Pair<word> >& interfaceNames,
 560     labelList& interfaceSizes,
 561     labelList& faceToInterface
 562 )
 563 {
 564     // From region-region to faceZone (or -1) to number of faces.
 565
 566     EdgeMap<Map<label> > regionsToSize;
 567
 568
 569     // Internal faces
 570     // ~~~~~~~~~~~~~~
 571
 572     forAll(mesh.faceNeighbour(), faceI)
 573     {
 574         label ownRegion = cellRegion[mesh.faceOwner()[faceI]];
 575         label neiRegion = cellRegion[mesh.faceNeighbour()[faceI]];
 576
 577         if (ownRegion != neiRegion)
 578         {
 579             addToInterface
 580             (
 581                 mesh,
 582                 (useFaceZones ? mesh.faceZones().whichZone(faceI) : -1),
 583                 ownRegion,
 584                 neiRegion,
 585                 regionsToSize
 586             );
 587         }
 588     }
 589
 590     // Boundary faces
 591     // ~~~~~~~~~~~~~~
 592
 593     // Neighbour cellRegion.
 594     labelList coupledRegion(mesh.nFaces()-mesh.nInternalFaces());
 595
 596     forAll(coupledRegion, i)
 597     {
 598         label cellI = mesh.faceOwner()[i+mesh.nInternalFaces()];
 599         coupledRegion[i] = cellRegion[cellI];
 600     }
 601     syncTools::swapBoundaryFaceList(mesh, coupledRegion);
 602
 603     forAll(coupledRegion, i)
 604     {
 605         label faceI = i+mesh.nInternalFaces();
 606         label ownRegion = cellRegion[mesh.faceOwner()[faceI]];
 607         label neiRegion = coupledRegion[i];
 608
 609         if (ownRegion != neiRegion)
 610         {
 611             addToInterface
 612             (
 613                 mesh,
 614                 (useFaceZones ? mesh.faceZones().whichZone(faceI) : -1),
 615                 ownRegion,
 616                 neiRegion,
 617                 regionsToSize
 618             );
 619         }
 620     }
 621
 622
 623     if (Pstream::parRun())
 624     {
 625         if (Pstream::master())
 626         {
 627             // Receive and add to my sizes
 628             for
 629             (
 630                 int slave=Pstream::firstSlave();
 631                 slave<=Pstream::lastSlave();
 632                 slave++
 633             )
 634             {
 635                 IPstream fromSlave(Pstream::blocking, slave);
 636
 637                 EdgeMap<Map<label> > slaveSizes(fromSlave);
 638
 639                 forAllConstIter(EdgeMap<Map<label> >, slaveSizes, slaveIter)
 640                 {
 641                     EdgeMap<Map<label> >::iterator masterIter =
 642                         regionsToSize.find(slaveIter.key());
 643
 644                     if (masterIter != regionsToSize.end())
 645                     {
 646                         // Same inter-region
 647                         const Map<label>& slaveInfo = slaveIter();
 648                         Map<label>& masterInfo = masterIter();
 649
 650                         forAllConstIter(Map<label>, slaveInfo, iter)
 651                         {
 652                             label zoneID = iter.key();
 653                             label slaveSize = iter();
 654
 655                             Map<label>::iterator zoneFnd = masterInfo.find
 656                             (
 657                                 zoneID
 658                             );
 659                             if (zoneFnd != masterInfo.end())
 660                             {
 661                                 zoneFnd() += slaveSize;
 662                             }
 663                             else
 664                             {
 665                                 masterInfo.insert(zoneID, slaveSize);
 666                             }
 667                         }
 668                     }
 669                     else
 670                     {
 671                         regionsToSize.insert(slaveIter.key(), slaveIter());
 672                     }
 673                 }
 674             }
 675         }
 676         else
 677         {
 678             // Send to master
 679             {
 680                 OPstream toMaster(Pstream::blocking, Pstream::masterNo());
 681                 toMaster << regionsToSize;
 682             }
 683         }
 684     }
 685
 686     // Rework
 687
 688     Pstream::scatter(regionsToSize);
 689
 690
 691
 692     // Now we have the global sizes of all inter-regions.
 693     // Invert this on master and distribute.
 694     label nInterfaces = 0;
 695     forAllConstIter(EdgeMap<Map<label> >, regionsToSize, iter)
 696     {
 697         const Map<label>& info = iter();
 698         nInterfaces += info.size();
 699     }
 700
 701     interfaces.setSize(nInterfaces);
 702     interfaceNames.setSize(nInterfaces);
 703     interfaceSizes.setSize(nInterfaces);
 704     EdgeMap<Map<label> > regionsToInterface(nInterfaces);
 705
 706     nInterfaces = 0;
 707     forAllConstIter(EdgeMap<Map<label> >, regionsToSize, iter)
 708     {
 709         const edge& e = iter.key();
 710         const word& name0 = regionNames[e[0]];
 711         const word& name1 = regionNames[e[1]];
 712
 713         const Map<label>& info = iter();
 714         forAllConstIter(Map<label>, info, infoIter)
 715         {
 716             interfaces[nInterfaces] = iter.key();
 717             label zoneID = infoIter.key();
 718             if (zoneID == -1)
 719             {
 720                 interfaceNames[nInterfaces] = Pair<word>
 721                 (
 722                     name0 + "_to_" + name1,
 723                     name1 + "_to_" + name0
 724                 );
 725             }
 726             else
 727             {
 728                 const word& zoneName = mesh.faceZones()[zoneID].name();
 729                 interfaceNames[nInterfaces] = Pair<word>
 730                 (
 731                     zoneName + "_" + name0 + "_to_" + name1,
 732                     zoneName + "_" + name1 + "_to_" + name0
 733                 );
 734             }
 735             interfaceSizes[nInterfaces] = infoIter();
 736
 737             Map<label> zoneAndInterface;
 738             zoneAndInterface.insert(zoneID, nInterfaces);
 739             regionsToInterface.insert(e, zoneAndInterface);
 740
 741             nInterfaces++;
 742         }
 743     }
 744
 745
 746     // Now all processor have consistent interface information
 747
 748     Pstream::scatter(interfaces);
 749     Pstream::scatter(interfaceNames);
 750     Pstream::scatter(interfaceSizes);
 751     Pstream::scatter(regionsToInterface);
 752
 753     // Mark all inter-region faces.
 754     faceToInterface.setSize(mesh.nFaces(), -1);
 755
 756     forAll(mesh.faceNeighbour(), faceI)
 757     {
 758         label ownRegion = cellRegion[mesh.faceOwner()[faceI]];
 759         label neiRegion = cellRegion[mesh.faceNeighbour()[faceI]];
 760
 761         if (ownRegion != neiRegion)
 762         {
 763             label zoneID = -1;
 764             if (useFaceZones)
 765             {
 766                 zoneID = mesh.faceZones().whichZone(faceI);
 767             }
 768
 769             edge interface
 770             (
 771                 min(ownRegion, neiRegion),
 772                 max(ownRegion, neiRegion)
 773             );
 774
 775             faceToInterface[faceI] = regionsToInterface[interface][zoneID];
 776         }
 777     }
 778     forAll(coupledRegion, i)
 779     {
 780         label faceI = i+mesh.nInternalFaces();
 781         label ownRegion = cellRegion[mesh.faceOwner()[faceI]];
 782         label neiRegion = coupledRegion[i];
 783
 784         if (ownRegion != neiRegion)
 785         {
 786             label zoneID = -1;
 787             if (useFaceZones)
 788             {
 789                 zoneID = mesh.faceZones().whichZone(faceI);
 790             }
 791
 792             edge interface
 793             (
 794                 min(ownRegion, neiRegion),
 795                 max(ownRegion, neiRegion)
 796             );
 797
 798             faceToInterface[faceI] = regionsToInterface[interface][zoneID];
 799         }
 800     }
 801 }
 802
 803
 804 // Create mesh for region.
 805 autoPtr<mapPolyMesh> createRegionMesh
 806 (
 807     const fvMesh& mesh,
 808     // Region info
 809     const labelList& cellRegion,
 810     const label regionI,
 811     const word& regionName,
 812     // Interface info
 813     const labelList& interfacePatches,
 814     const labelList& faceToInterface,
 815
 816     autoPtr<fvMesh>& newMesh
 817 )
 818 {
 819     // Create dummy system/fv*
 820     {
 821         IOobject io
 822         (
 823             "fvSchemes",
 824             mesh.time().system(),
 825             regionName,
 826             mesh,
 827             IOobject::NO_READ,
 828             IOobject::NO_WRITE,
 829             false
 830         );
 831
 832         Info<< "Testing:" << io.objectPath() << endl;
 833
 834         if (!io.headerOk())
 835         // if (!exists(io.objectPath()))
 836         {
 837             Info<< "Writing dummy " << regionName/io.name() << endl;
 838             dictionary dummyDict;
 839             dictionary divDict;
 840             dummyDict.add("divSchemes", divDict);
 841             dictionary gradDict;
 842             dummyDict.add("gradSchemes", gradDict);
 843             dictionary laplDict;
 844             dummyDict.add("laplacianSchemes", laplDict);
 845
 846             IOdictionary(io, dummyDict).regIOobject::write();
 847         }
 848     }
 849     {
 850         IOobject io
 851         (
 852             "fvSolution",
 853             mesh.time().system(),
 854             regionName,
 855             mesh,
 856             IOobject::NO_READ,
 857             IOobject::NO_WRITE,
 858             false
 859         );
 860
 861         if (!io.headerOk())
 862         //if (!exists(io.objectPath()))
 863         {
 864             Info<< "Writing dummy " << regionName/io.name() << endl;
 865             dictionary dummyDict;
 866             IOdictionary(io, dummyDict).regIOobject::write();
 867         }
 868     }
 869
 870
 871     // Neighbour cellRegion.
 872     labelList coupledRegion(mesh.nFaces()-mesh.nInternalFaces());
 873
 874     forAll(coupledRegion, i)
 875     {
 876         label cellI = mesh.faceOwner()[i+mesh.nInternalFaces()];
 877         coupledRegion[i] = cellRegion[cellI];
 878     }
 879     syncTools::swapBoundaryFaceList(mesh, coupledRegion);
 880
 881
 882     // Topology change container. Start off from existing mesh.
 883     polyTopoChange meshMod(mesh);
 884
 885     // Cell remover engine
 886     removeCells cellRemover(mesh);
 887
 888     // Select all but region cells
 889     labelList cellsToRemove(getNonRegionCells(cellRegion, regionI));
 890
 891     // Find out which faces will get exposed. Note that this
 892     // gets faces in mesh face order. So both regions will get same
 893     // face in same order (important!)
 894     labelList exposedFaces = cellRemover.getExposedFaces(cellsToRemove);
 895
 896     labelList exposedPatchIDs(exposedFaces.size());
 897     forAll(exposedFaces, i)
 898     {
 899         label faceI = exposedFaces[i];
 900         label interfaceI = faceToInterface[faceI];
 901
 902         label ownRegion = cellRegion[mesh.faceOwner()[faceI]];
 903         label neiRegion = -1;
 904
 905         if (mesh.isInternalFace(faceI))
 906         {
 907             neiRegion = cellRegion[mesh.faceNeighbour()[faceI]];
 908         }
 909         else
 910         {
 911             neiRegion = coupledRegion[faceI-mesh.nInternalFaces()];
 912         }
 913
 914
 915         // Check which side is being kept - determines which of the two
 916         // patches will be used.
 917
 918         label otherRegion = -1;
 919
 920         if (ownRegion == regionI && neiRegion != regionI)
 921         {
 922             otherRegion = neiRegion;
 923         }
 924         else if (ownRegion != regionI && neiRegion == regionI)
 925         {
 926             otherRegion = ownRegion;
 927         }
 928         else
 929         {
 930             FatalErrorIn("createRegionMesh(..)")
 931                 << "Exposed face:" << faceI
 932                 << " fc:" << mesh.faceCentres()[faceI]
 933                 << " has owner region " << ownRegion
 934                 << " and neighbour region " << neiRegion
 935                 << " when handling region:" << regionI
 936                 << exit(FatalError);
 937         }
 938
 939         // Find the patch.
 940         if (regionI < otherRegion)
 941         {
 942             exposedPatchIDs[i] = interfacePatches[interfaceI];
 943         }
 944         else
 945         {
 946             exposedPatchIDs[i] = interfacePatches[interfaceI]+1;
 947         }
 948     }
 949
 950     // Remove faces
 951     cellRemover.setRefinement
 952     (
 953         cellsToRemove,
 954         exposedFaces,
 955         exposedPatchIDs,
 956         meshMod
 957     );
 958
 959     autoPtr<mapPolyMesh> map = meshMod.makeMesh
 960     (
 961         newMesh,
 962         IOobject
 963         (
 964             regionName,
 965             mesh.time().timeName(),
 966             mesh.time(),
 967             IOobject::NO_READ,
 968             IOobject::AUTO_WRITE
 969         ),
 970         mesh
 971     );
 972
 973     return map;
 974 }
 975
 976
 977 void createAndWriteRegion
 978 (
 979     const fvMesh& mesh,
 980     const labelList& cellRegion,
 981     const wordList& regionNames,
 982     const labelList& faceToInterface,
 983     const labelList& interfacePatches,
 984     const label regionI,
 985     const word& newMeshInstance
 986 )
 987 {
 988     Info<< "Creating mesh for region " << regionI
 989         << ' ' << regionNames[regionI] << endl;
 990
 991     autoPtr<fvMesh> newMesh;
 992     autoPtr<mapPolyMesh> map = createRegionMesh
 993     (
 994         mesh,
 995         cellRegion,
 996         regionI,
 997         regionNames[regionI],
 998         interfacePatches,
 999         faceToInterface,
1000         newMesh
1001     );
1002
1003
1004     // Make map of all added patches
1005     labelHashSet addedPatches(2*interfacePatches.size());
1006     forAll(interfacePatches, interfaceI)
1007     {
1008         addedPatches.insert(interfacePatches[interfaceI]);
1009         addedPatches.insert(interfacePatches[interfaceI]+1);
1010     }
1011
1012     Info<< "Mapping fields" << endl;
1013
1014     // Map existing fields
1015     newMesh().updateMesh(map());
1016
1017     // Add subsetted fields
1018     subsetVolFields<volScalarField>
1019     (
1020         mesh,
1021         newMesh(),
1022         map().cellMap(),
1023         map().faceMap(),
1024         addedPatches
1025     );
1026     subsetVolFields<volVectorField>
1027     (
1028         mesh,
1029         newMesh(),
1030         map().cellMap(),
1031         map().faceMap(),
1032         addedPatches
1033     );
1034     subsetVolFields<volSphericalTensorField>
1035     (
1036         mesh,
1037         newMesh(),
1038         map().cellMap(),
1039         map().faceMap(),
1040         addedPatches
1041     );
1042     subsetVolFields<volSymmTensorField>
1043     (
1044         mesh,
1045         newMesh(),
1046         map().cellMap(),
1047         map().faceMap(),
1048         addedPatches
1049     );
1050     subsetVolFields<volTensorField>
1051     (
1052         mesh,
1053         newMesh(),
1054         map().cellMap(),
1055         map().faceMap(),
1056         addedPatches
1057     );
1058
1059     subsetSurfaceFields<surfaceScalarField>
1060     (
1061         mesh,
1062         newMesh(),
1063         map().faceMap(),
1064         addedPatches
1065     );
1066     subsetSurfaceFields<surfaceVectorField>
1067     (
1068         mesh,
1069         newMesh(),
1070         map().faceMap(),
1071         addedPatches
1072     );
1073     subsetSurfaceFields<surfaceSphericalTensorField>
1074     (
1075         mesh,
1076         newMesh(),
1077         map().faceMap(),
1078         addedPatches
1079     );
1080     subsetSurfaceFields<surfaceSymmTensorField>
1081     (
1082         mesh,
1083         newMesh(),
1084         map().faceMap(),
1085         addedPatches
1086     );
1087     subsetSurfaceFields<surfaceTensorField>
1088     (
1089         mesh,
1090         newMesh(),
1091         map().faceMap(),
1092         addedPatches
1093     );
1094
1095
1096     const polyBoundaryMesh& newPatches = newMesh().boundaryMesh();
1097     newPatches.checkParallelSync(true);
1098
1099     // Delete empty patches
1100     // ~~~~~~~~~~~~~~~~~~~~
1101
1102     // Create reordering list to move patches-to-be-deleted to end
1103     labelList oldToNew(newPatches.size(), -1);
1104     label newI = 0;
1105
1106     Info<< "Deleting empty patches" << endl;
1107
1108     // Assumes all non-proc boundaries are on all processors!
1109     forAll(newPatches, patchI)
1110     {
1111         const polyPatch& pp = newPatches[patchI];
1112
1113         if (!isA<processorPolyPatch>(pp))
1114         {
1115             if (returnReduce(pp.size(), sumOp<label>()) > 0)
1116             {
1117                 oldToNew[patchI] = newI++;
1118             }
1119         }
1120     }
1121
1122     // Same for processor patches (but need no reduction)
1123     forAll(newPatches, patchI)
1124     {
1125         const polyPatch& pp = newPatches[patchI];
1126
1127         if (isA<processorPolyPatch>(pp) && pp.size())
1128         {
1129             oldToNew[patchI] = newI++;
1130         }
1131     }
1132
1133     const label nNewPatches = newI;
1134
1135     // Move all deleteable patches to the end
1136     forAll(oldToNew, patchI)
1137     {
1138         if (oldToNew[patchI] == -1)
1139         {
1140             oldToNew[patchI] = newI++;
1141         }
1142     }
1143
1144     reorderPatches(newMesh(), oldToNew, nNewPatches);
1145
1146
1147     Info<< "Writing new mesh" << endl;
1148
1149     newMesh().setInstance(newMeshInstance);
1150     newMesh().write();
1151
1152     // Write addressing files like decomposePar
1153     Info<< "Writing addressing to base mesh" << endl;
1154
1155     labelIOList pointProcAddressing
1156     (
1157         IOobject
1158         (
1159             "pointRegionAddressing",
1160             newMesh().facesInstance(),
1161             newMesh().meshSubDir,
1162             newMesh(),
1163             IOobject::NO_READ,
1164             IOobject::NO_WRITE,
1165             false
1166         ),
1167         map().pointMap()
1168     );
1169     Info<< "Writing map " << pointProcAddressing.name()
1170         << " from region" << regionI
1171         << " points back to base mesh." << endl;
1172     pointProcAddressing.write();
1173
1174     labelIOList faceProcAddressing
1175     (
1176         IOobject
1177         (
1178             "faceRegionAddressing",
1179             newMesh().facesInstance(),
1180             newMesh().meshSubDir,
1181             newMesh(),
1182             IOobject::NO_READ,
1183             IOobject::NO_WRITE,
1184             false
1185         ),
1186         newMesh().nFaces()
1187     );
1188     forAll(faceProcAddressing, faceI)
1189     {
1190         // face + turning index. (see decomposePar)
1191         // Is the face pointing in the same direction?
1192         label oldFaceI = map().faceMap()[faceI];
1193
1194         if
1195         (
1196             map().cellMap()[newMesh().faceOwner()[faceI]]
1197          == mesh.faceOwner()[oldFaceI]
1198         )
1199         {
1200             faceProcAddressing[faceI] = oldFaceI+1;
1201         }
1202         else
1203         {
1204             faceProcAddressing[faceI] = -(oldFaceI+1);
1205         }
1206     }
1207     Info<< "Writing map " << faceProcAddressing.name()
1208         << " from region" << regionI
1209         << " faces back to base mesh." << endl;
1210     faceProcAddressing.write();
1211
1212     labelIOList cellProcAddressing
1213     (
1214         IOobject
1215         (
1216             "cellRegionAddressing",
1217             newMesh().facesInstance(),
1218             newMesh().meshSubDir,
1219             newMesh(),
1220             IOobject::NO_READ,
1221             IOobject::NO_WRITE,
1222             false
1223         ),
1224         map().cellMap()
1225     );
1226     Info<< "Writing map " <<cellProcAddressing.name()
1227         << " from region" << regionI
1228         << " cells back to base mesh." << endl;
1229     cellProcAddressing.write();
1230
1231     labelIOList boundaryProcAddressing
1232     (
1233         IOobject
1234         (
1235             "boundaryRegionAddressing",
1236             newMesh().facesInstance(),
1237             newMesh().meshSubDir,
1238             newMesh(),
1239             IOobject::NO_READ,
1240             IOobject::NO_WRITE,
1241             false
1242         ),
1243         labelList(nNewPatches, -1)
1244     );
1245     forAll(oldToNew, i)
1246     {
1247         if (!addedPatches.found(i))
1248         {
1249             label newI = oldToNew[i];
1250             if (newI >= 0 && newI < nNewPatches)
1251             {
1252                 boundaryProcAddressing[oldToNew[i]] = i;
1253             }
1254         }
1255     }
1256     Info<< "Writing map " << boundaryProcAddressing.name()
1257         << " from region" << regionI
1258         << " boundary back to base mesh." << endl;
1259     boundaryProcAddressing.write();
1260 }
1261
1262
1263 // Create for every region-region interface with non-zero size two patches.
1264 // First one is for minimumregion to maximumregion.
1265 // Note that patches get created in same order on all processors (if parallel)
1266 // since looping over synchronised 'interfaces'.
1267 labelList addRegionPatches
1268 (
1269     fvMesh& mesh,
1270     const wordList& regionNames,
1271     const edgeList& interfaces,
1272     const List<Pair<word> >& interfaceNames
1273 )
1274 {
1275     Info<< nl << "Adding patches" << nl << endl;
1276
1277     labelList interfacePatches(interfaces.size());
1278
1279     forAll(interfaces, interI)
1280     {
1281         const edge& e = interfaces[interI];
1282         const Pair<word>& names = interfaceNames[interI];
1283
1284         //Info<< "For interface " << interI
1285         //    << " between regions " << e
1286         //    << " trying to add patches " << names << endl;
1287
1288
1289         directMappedWallPolyPatch patch1
1290         (
1291             names[0],
1292             0,                  // overridden
1293             0,                  // overridden
1294             0,                  // overridden
1295             regionNames[e[1]],  // sampleRegion
1296             directMappedPatchBase::NEARESTPATCHFACE,
1297             names[1],           // samplePatch
1298             point::zero,        // offset
1299             mesh.boundaryMesh()
1300         );
1301
1302         interfacePatches[interI] = addPatch(mesh, patch1);
1303
1304         directMappedWallPolyPatch patch2
1305         (
1306             names[1],
1307             0,
1308             0,
1309             0,
1310             regionNames[e[0]],  // sampleRegion
1311             directMappedPatchBase::NEARESTPATCHFACE,
1312             names[0],
1313             point::zero,        // offset
1314             mesh.boundaryMesh()
1315         );
1316         addPatch(mesh, patch2);
1317
1318         Info<< "For interface between region " << regionNames[e[0]]
1319             << " and " << regionNames[e[1]] << " added patches" << endl
1320             << "    " << interfacePatches[interI]
1321             << "\t" << mesh.boundaryMesh()[interfacePatches[interI]].name()
1322             << endl
1323             << "    " << interfacePatches[interI]+1
1324             << "\t" << mesh.boundaryMesh()[interfacePatches[interI]+1].name()
1325             << endl;
1326     }
1327     return interfacePatches;
1328 }
1329
1330
1331 // Find region that covers most of cell zone
1332 label findCorrespondingRegion
1333 (
1334     const labelList& existingZoneID,    // per cell the (unique) zoneID
1335     const labelList& cellRegion,
1336     const label nCellRegions,
1337     const label zoneI,
1338     const label minOverlapSize
1339 )
1340 {
1341     // Per region the number of cells in zoneI
1342     labelList cellsInZone(nCellRegions, 0);
1343
1344     forAll(cellRegion, cellI)
1345     {
1346         if (existingZoneID[cellI] == zoneI)
1347         {
1348             cellsInZone[cellRegion[cellI]]++;
1349         }
1350     }
1351
1352     Pstream::listCombineGather(cellsInZone, plusEqOp<label>());
1353     Pstream::listCombineScatter(cellsInZone);
1354
1355     // Pick region with largest overlap of zoneI
1356     label regionI = findMax(cellsInZone);
1357
1358
1359     if (cellsInZone[regionI] < minOverlapSize)
1360     {
1361         // Region covers too little of zone. Not good enough.
1362         regionI = -1;
1363     }
1364     else
1365     {
1366         // Check that region contains no cells that aren't in cellZone.
1367         forAll(cellRegion, cellI)
1368         {
1369             if (cellRegion[cellI] == regionI && existingZoneID[cellI] != zoneI)
1370             {
1371                 // cellI in regionI but not in zoneI
1372                 regionI = -1;
1373                 break;
1374             }
1375         }
1376         // If one in error, all should be in error. Note that branch gets taken
1377         // on all procs.
1378         reduce(regionI, minOp<label>());
1379     }
1380
1381     return regionI;
1382 }
1383
1384
1385 // Get zone per cell
1386 // - non-unique zoning
1387 // - coupled zones
1388 void getZoneID
1389 (
1390     const polyMesh& mesh,
1391     const cellZoneMesh& cellZones,
1392     labelList& zoneID,
1393     labelList& neiZoneID
1394 )
1395 {
1396     // Existing zoneID
1397     zoneID.setSize(mesh.nCells());
1398     zoneID = -1;
1399
1400     forAll(cellZones, zoneI)
1401     {
1402         const cellZone& cz = cellZones[zoneI];
1403
1404         forAll(cz, i)
1405         {
1406             label cellI = cz[i];
1407             if (zoneID[cellI] == -1)
1408             {
1409                 zoneID[cellI] = zoneI;
1410             }
1411             else
1412             {
1413                 FatalErrorIn("getZoneID(..)")
1414                     << "Cell " << cellI << " with cell centre "
1415                     << mesh.cellCentres()[cellI]
1416                     << " is multiple zones. This is not allowed." << endl
1417                     << "It is in zone " << cellZones[zoneID[cellI]].name()
1418                     << " and in zone " << cellZones[zoneI].name()
1419                     << exit(FatalError);
1420             }
1421         }
1422     }
1423
1424     // Neighbour zoneID.
1425     neiZoneID.setSize(mesh.nFaces()-mesh.nInternalFaces());
1426
1427     forAll(neiZoneID, i)
1428     {
1429         neiZoneID[i] = zoneID[mesh.faceOwner()[i+mesh.nInternalFaces()]];
1430     }
1431     syncTools::swapBoundaryFaceList(mesh, neiZoneID);
1432 }
1433
1434
1435 void matchRegions
1436 (
1437     const bool sloppyCellZones,
1438     const polyMesh& mesh,
1439
1440     const label nCellRegions,
1441     const labelList& cellRegion,
1442
1443     labelList& regionToZone,
1444     wordList& regionNames,
1445     labelList& zoneToRegion
1446 )
1447 {
1448     const cellZoneMesh& cellZones = mesh.cellZones();
1449
1450     regionToZone.setSize(nCellRegions, -1);
1451     regionNames.setSize(nCellRegions);
1452     zoneToRegion.setSize(cellZones.size(), -1);
1453
1454     // Get current per cell zoneID
1455     labelList zoneID(mesh.nCells(), -1);
1456     labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
1457     getZoneID(mesh, cellZones, zoneID, neiZoneID);
1458
1459     // Sizes per cellzone
1460     labelList zoneSizes(cellZones.size(), 0);
1461     {
1462         List<wordList> zoneNames(Pstream::nProcs());
1463         zoneNames[Pstream::myProcNo()] = cellZones.names();
1464         Pstream::gatherList(zoneNames);
1465         Pstream::scatterList(zoneNames);
1466
1467         forAll(zoneNames, procI)
1468         {
1469             if (zoneNames[procI] != zoneNames[0])
1470             {
1471                 FatalErrorIn("matchRegions(..)")
1472                     << "cellZones not synchronised across processors." << endl
1473                     << "Master has cellZones " << zoneNames[0] << endl
1474                     << "Processor " << procI
1475                     << " has cellZones " << zoneNames[procI]
1476                     << exit(FatalError);
1477             }
1478         }
1479
1480         forAll(cellZones, zoneI)
1481         {
1482             zoneSizes[zoneI] = returnReduce
1483             (
1484                 cellZones[zoneI].size(),
1485                 sumOp<label>()
1486             );
1487         }
1488     }
1489
1490
1491     if (sloppyCellZones)
1492     {
1493         Info<< "Trying to match regions to existing cell zones;"
1494             << " region can be subset of cell zone." << nl << endl;
1495
1496         forAll(cellZones, zoneI)
1497         {
1498             label regionI = findCorrespondingRegion
1499             (
1500                 zoneID,
1501                 cellRegion,
1502                 nCellRegions,
1503                 zoneI,
1504                 label(0.5*zoneSizes[zoneI]) // minimum overlap
1505             );
1506
1507             if (regionI != -1)
1508             {
1509                 Info<< "Sloppily matched region " << regionI
1510                     //<< " size " << regionSizes[regionI]
1511                     << " to zone " << zoneI << " size " << zoneSizes[zoneI]
1512                     << endl;
1513                 zoneToRegion[zoneI] = regionI;
1514                 regionToZone[regionI] = zoneI;
1515                 regionNames[regionI] = cellZones[zoneI].name();
1516             }
1517         }
1518     }
1519     else
1520     {
1521         Info<< "Trying to match regions to existing cell zones." << nl << endl;
1522
1523         forAll(cellZones, zoneI)
1524         {
1525             label regionI = findCorrespondingRegion
1526             (
1527                 zoneID,
1528                 cellRegion,
1529                 nCellRegions,
1530                 zoneI,
1531                 1               // minimum overlap
1532             );
1533
1534             if (regionI != -1)
1535             {
1536                 zoneToRegion[zoneI] = regionI;
1537                 regionToZone[regionI] = zoneI;
1538                 regionNames[regionI] = cellZones[zoneI].name();
1539             }
1540         }
1541     }
1542     // Allocate region names for unmatched regions.
1543     forAll(regionToZone, regionI)
1544     {
1545         if (regionToZone[regionI] == -1)
1546         {
1547             regionNames[regionI] = "domain" + Foam::name(regionI);
1548         }
1549     }
1550 }
1551
1552
1553 void writeCellToRegion(const fvMesh& mesh, const labelList& cellRegion)
1554 {
1555     // Write to manual decomposition option
1556     {
1557         labelIOList cellToRegion
1558         (
1559             IOobject
1560             (
1561                 "cellToRegion",
1562                 mesh.facesInstance(),
1563                 mesh,
1564                 IOobject::NO_READ,
1565                 IOobject::NO_WRITE,
1566                 false
1567             ),
1568             cellRegion
1569         );
1570         cellToRegion.write();
1571
1572         Info<< "Writing region per cell file (for manual decomposition) to "
1573             << cellToRegion.objectPath() << nl << endl;
1574     }
1575     // Write for postprocessing
1576     {
1577         volScalarField cellToRegion
1578         (
1579             IOobject
1580             (
1581                 "cellToRegion",
1582                 mesh.time().timeName(),
1583                 mesh,
1584                 IOobject::NO_READ,
1585                 IOobject::NO_WRITE,
1586                 false
1587             ),
1588             mesh,
1589             dimensionedScalar("zero", dimless, 0),
1590             zeroGradientFvPatchScalarField::typeName
1591         );
1592         forAll(cellRegion, cellI)
1593         {
1594             cellToRegion[cellI] = cellRegion[cellI];
1595         }
1596         cellToRegion.write();
1597
1598         Info<< "Writing region per cell as volScalarField to "
1599             << cellToRegion.objectPath() << nl << endl;
1600     }
1601 }
1602
1603
1604
1605 // Main program:
1606
1607 int main(int argc, char *argv[])
1608 {
1609     argList::addNote
1610     (
1611         "splits mesh into multiple regions (detected by walking across faces)"
1612     );
1613     #include "addOverwriteOption.H"
1614     argList::addBoolOption
1615     (
1616         "cellZones",
1617         "additionally split cellZones off into separate regions"
1618     );
1619     argList::addBoolOption
1620     (
1621         "cellZonesOnly",
1622         "use cellZones only to split mesh into regions; do not use walking"
1623     );
1624     argList::addOption
1625     (
1626         "cellZonesFileOnly",
1627         "file",
1628         "like -cellZonesOnly, but use specified file"
1629     );
1630     argList::addOption
1631     (
1632         "blockedFaces",
1633         "faceSet",
1634         "specify additional region boundaries that walking does not cross"
1635     );
1636     argList::addBoolOption
1637     (
1638         "makeCellZones",
1639         "place cells into cellZones instead of splitting mesh"
1640     );
1641     argList::addBoolOption
1642     (
1643         "largestOnly",
1644         "only write largest region"
1645     );
1646     argList::addOption
1647     (
1648         "insidePoint",
1649         "point",
1650         "only write region containing point"
1651     );
1652     argList::addBoolOption
1653     (
1654         "detectOnly",
1655         "do not write mesh"
1656     );
1657     argList::addBoolOption
1658     (
1659         "sloppyCellZones",
1660         "try to match heuristically regions to existing cell zones"
1661     );
1662     argList::addBoolOption
1663     (
1664         "useFaceZones",
1665         "use faceZones to patch inter-region faces instead of single patch"
1666     );
1667
1668     #include "setRootCase.H"
1669     #include "createTime.H"
1670     runTime.functionObjects().off();
1671     #include "createMesh.H"
1672     const word oldInstance = mesh.pointsInstance();
1673
1674     word blockedFacesName;
1675     if (args.optionReadIfPresent("blockedFaces", blockedFacesName))
1676     {
1677         Info<< "Reading blocked internal faces from faceSet "
1678             << blockedFacesName << nl << endl;
1679     }
1680
1681     const bool makeCellZones    = args.optionFound("makeCellZones");
1682     const bool largestOnly      = args.optionFound("largestOnly");
1683     const bool insidePoint      = args.optionFound("insidePoint");
1684     const bool useCellZones     = args.optionFound("cellZones");
1685     const bool useCellZonesOnly = args.optionFound("cellZonesOnly");
1686     const bool useCellZonesFile = args.optionFound("cellZonesFileOnly");
1687     const bool overwrite        = args.optionFound("overwrite");
1688     const bool detectOnly       = args.optionFound("detectOnly");
1689     const bool sloppyCellZones  = args.optionFound("sloppyCellZones");
1690     const bool useFaceZones     = args.optionFound("useFaceZones");
1691
1692     if
1693     (
1694         (useCellZonesOnly || useCellZonesFile)
1695      && (useCellZones || blockedFacesName.size())
1696     )
1697     {
1698         FatalErrorIn(args.executable())
1699             << "You cannot specify both -cellZonesOnly or -cellZonesFileOnly"
1700             << " (which specify complete split)"
1701             << " in combination with -blockedFaces or -cellZones"
1702             << " (which imply a split based on topology)"
1703             << exit(FatalError);
1704     }
1705
1706
1707     if (useFaceZones)
1708     {
1709         Info<< "Using current faceZones to divide inter-region interfaces"
1710             << " into multiple patches."
1711             << nl << endl;
1712     }
1713     else
1714     {
1715         Info<< "Creating single patch per inter-region interface."
1716             << nl << endl;
1717     }
1718
1719
1720
1721     if (insidePoint && largestOnly)
1722     {
1723         FatalErrorIn(args.executable())
1724             << "You cannot specify both -largestOnly"
1725             << " (keep region with most cells)"
1726             << " and -insidePoint (keep region containing point)"
1727             << exit(FatalError);
1728     }
1729
1730
1731     const cellZoneMesh& cellZones = mesh.cellZones();
1732
1733     // Existing zoneID
1734     labelList zoneID(mesh.nCells(), -1);
1735     // Neighbour zoneID.
1736     labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
1737     getZoneID(mesh, cellZones, zoneID, neiZoneID);
1738
1739
1740
1741     // Determine per cell the region it belongs to
1742     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1743
1744     // cellRegion is the labelList with the region per cell.
1745     labelList cellRegion;
1746     // Region per zone
1747     labelList regionToZone;
1748     // Name of region
1749     wordList regionNames;
1750     // Zone to region
1751     labelList zoneToRegion;
1752
1753     label nCellRegions = 0;
1754     if (useCellZonesOnly)
1755     {
1756         Info<< "Using current cellZones to split mesh into regions."
1757             << " This requires all"
1758             << " cells to be in one and only one cellZone." << nl << endl;
1759
1760         label unzonedCellI = findIndex(zoneID, -1);
1761         if (unzonedCellI != -1)
1762         {
1763             FatalErrorIn(args.executable())
1764                 << "For the cellZonesOnly option all cells "
1765                 << "have to be in a cellZone." << endl
1766                 << "Cell " << unzonedCellI
1767                 << " at" << mesh.cellCentres()[unzonedCellI]
1768                 << " is not in a cellZone. There might be more unzoned cells."
1769                 << exit(FatalError);
1770         }
1771         cellRegion = zoneID;
1772         nCellRegions = gMax(cellRegion)+1;
1773         regionToZone.setSize(nCellRegions);
1774         regionNames.setSize(nCellRegions);
1775         zoneToRegion.setSize(cellZones.size(), -1);
1776         for (label regionI = 0; regionI < nCellRegions; regionI++)
1777         {
1778             regionToZone[regionI] = regionI;
1779             zoneToRegion[regionI] = regionI;
1780             regionNames[regionI] = cellZones[regionI].name();
1781         }
1782     }
1783     else if (useCellZonesFile)
1784     {
1785         const word zoneFile = args.option("cellZonesFileOnly");
1786         Info<< "Reading split from cellZones file " << zoneFile << endl
1787             << "This requires all"
1788             << " cells to be in one and only one cellZone." << nl << endl;
1789
1790         cellZoneMesh newCellZones
1791         (
1792             IOobject
1793             (
1794                 zoneFile,
1795                 mesh.facesInstance(),
1796                 polyMesh::meshSubDir,
1797                 mesh,
1798                 IOobject::MUST_READ,
1799                 IOobject::NO_WRITE,
1800                 false
1801             ),
1802             mesh
1803         );
1804
1805         labelList newZoneID(mesh.nCells(), -1);
1806         labelList newNeiZoneID(mesh.nFaces()-mesh.nInternalFaces());
1807         getZoneID(mesh, newCellZones, newZoneID, newNeiZoneID);
1808
1809         label unzonedCellI = findIndex(newZoneID, -1);
1810         if (unzonedCellI != -1)
1811         {
1812             FatalErrorIn(args.executable())
1813                 << "For the cellZonesFileOnly option all cells "
1814                 << "have to be in a cellZone." << endl
1815                 << "Cell " << unzonedCellI
1816                 << " at" << mesh.cellCentres()[unzonedCellI]
1817                 << " is not in a cellZone. There might be more unzoned cells."
1818                 << exit(FatalError);
1819         }
1820         cellRegion = newZoneID;
1821         nCellRegions = gMax(cellRegion)+1;
1822         zoneToRegion.setSize(newCellZones.size(), -1);
1823         regionToZone.setSize(nCellRegions);
1824         regionNames.setSize(nCellRegions);
1825         for (label regionI = 0; regionI < nCellRegions; regionI++)
1826         {
1827             regionToZone[regionI] = regionI;
1828             zoneToRegion[regionI] = regionI;
1829             regionNames[regionI] = newCellZones[regionI].name();
1830         }
1831     }
1832     else
1833     {
1834         // Determine connected regions
1835         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
1836
1837         // Mark additional faces that are blocked
1838         boolList blockedFace;
1839
1840         // Read from faceSet
1841         if (blockedFacesName.size())
1842         {
1843             faceSet blockedFaceSet(mesh, blockedFacesName);
1844             Info<< "Read "
1845                 << returnReduce(blockedFaceSet.size(), sumOp<label>())
1846                 << " blocked faces from set " << blockedFacesName << nl << endl;
1847
1848             blockedFace.setSize(mesh.nFaces(), false);
1849
1850             forAllConstIter(faceSet, blockedFaceSet, iter)
1851             {
1852                 blockedFace[iter.key()] = true;
1853             }
1854         }
1855
1856         // Imply from differing cellZones
1857         if (useCellZones)
1858         {
1859             blockedFace.setSize(mesh.nFaces(), false);
1860
1861             for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
1862             {
1863                 label own = mesh.faceOwner()[faceI];
1864                 label nei = mesh.faceNeighbour()[faceI];
1865
1866                 if (zoneID[own] != zoneID[nei])
1867                 {
1868                     blockedFace[faceI] = true;
1869                 }
1870             }
1871
1872             // Different cellZones on either side of processor patch.
1873             forAll(neiZoneID, i)
1874             {
1875                 label faceI = i+mesh.nInternalFaces();
1876
1877                 if (zoneID[mesh.faceOwner()[faceI]] != neiZoneID[i])
1878                 {
1879                     blockedFace[faceI] = true;
1880                 }
1881             }
1882         }
1883
1884         // Do a topological walk to determine regions
1885         regionSplit regions(mesh, blockedFace);
1886         nCellRegions = regions.nRegions();
1887         cellRegion.transfer(regions);
1888
1889         // Make up region names. If possible match them to existing zones.
1890         matchRegions
1891         (
1892             sloppyCellZones,
1893             mesh,
1894             nCellRegions,
1895             cellRegion,
1896
1897             regionToZone,
1898             regionNames,
1899             zoneToRegion
1900         );
1901     }
1902
1903     Info<< endl << "Number of regions:" << nCellRegions << nl << endl;
1904
1905
1906     // Write decomposition to file
1907     writeCellToRegion(mesh, cellRegion);
1908
1909
1910
1911     // Sizes per region
1912     // ~~~~~~~~~~~~~~~~
1913
1914     labelList regionSizes(nCellRegions, 0);
1915
1916     forAll(cellRegion, cellI)
1917     {
1918         regionSizes[cellRegion[cellI]]++;
1919     }
1920     forAll(regionSizes, regionI)
1921     {
1922         reduce(regionSizes[regionI], sumOp<label>());
1923     }
1924
1925     Info<< "Region\tCells" << nl
1926         << "------\t-----" << endl;
1927
1928     forAll(regionSizes, regionI)
1929     {
1930         Info<< regionI << '\t' << regionSizes[regionI] << nl;
1931     }
1932     Info<< endl;
1933
1934
1935
1936     // Print region to zone
1937     Info<< "Region\tZone\tName" << nl
1938         << "------\t----\t----" << endl;
1939     forAll(regionToZone, regionI)
1940     {
1941         Info<< regionI << '\t' << regionToZone[regionI] << '\t'
1942             << regionNames[regionI] << nl;
1943     }
1944     Info<< endl;
1945
1946
1947
1948     // Since we're going to mess with patches and zones make sure all
1949     // is synchronised
1950     mesh.boundaryMesh().checkParallelSync(true);
1951     mesh.faceZones().checkParallelSync(true);
1952
1953
1954     // Interfaces
1955     // ----------
1956     // per interface:
1957     // - the two regions on either side
1958     // - the name
1959     // - the (global) size
1960     edgeList interfaces;
1961     List<Pair<word> > interfaceNames;
1962     labelList interfaceSizes;
1963     // per face the interface
1964     labelList faceToInterface;
1965
1966     getInterfaceSizes
1967     (
1968         mesh,
1969         useFaceZones,
1970         cellRegion,
1971         regionNames,
1972
1973         interfaces,
1974         interfaceNames,
1975         interfaceSizes,
1976         faceToInterface
1977     );
1978
1979     Info<< "Sizes of interfaces between regions:" << nl << nl
1980         << "Interface\tRegion\tRegion\tFaces" << nl
1981         << "---------\t------\t------\t-----" << endl;
1982
1983     forAll(interfaces, interI)
1984     {
1985         const edge& e = interfaces[interI];
1986
1987         Info<< interI
1988             << "\t\t" << e[0] << '\t' << e[1]
1989             << '\t' << interfaceSizes[interI] << nl;
1990     }
1991     Info<< endl;
1992
1993
1994     if (detectOnly)
1995     {
1996         return 0;
1997     }
1998
1999
2000     // Read objects in time directory
2001     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2002
2003     IOobjectList objects(mesh, runTime.timeName());
2004
2005     // Read vol fields.
2006
2007     PtrList<volScalarField> vsFlds;
2008     ReadFields(mesh, objects, vsFlds);
2009
2010     PtrList<volVectorField> vvFlds;
2011     ReadFields(mesh, objects, vvFlds);
2012
2013     PtrList<volSphericalTensorField> vstFlds;
2014     ReadFields(mesh, objects, vstFlds);
2015
2016     PtrList<volSymmTensorField> vsymtFlds;
2017     ReadFields(mesh, objects, vsymtFlds);
2018
2019     PtrList<volTensorField> vtFlds;
2020     ReadFields(mesh, objects, vtFlds);
2021
2022     // Read surface fields.
2023
2024     PtrList<surfaceScalarField> ssFlds;
2025     ReadFields(mesh, objects, ssFlds);
2026
2027     PtrList<surfaceVectorField> svFlds;
2028     ReadFields(mesh, objects, svFlds);
2029
2030     PtrList<surfaceSphericalTensorField> sstFlds;
2031     ReadFields(mesh, objects, sstFlds);
2032
2033     PtrList<surfaceSymmTensorField> ssymtFlds;
2034     ReadFields(mesh, objects, ssymtFlds);
2035
2036     PtrList<surfaceTensorField> stFlds;
2037     ReadFields(mesh, objects, stFlds);
2038
2039     Info<< endl;
2040
2041
2042     // Remove any demand-driven fields ('S', 'V' etc)
2043     mesh.clearOut();
2044
2045
2046     if (nCellRegions == 1)
2047     {
2048         Info<< "Only one region. Doing nothing." << endl;
2049     }
2050     else if (makeCellZones)
2051     {
2052         Info<< "Putting cells into cellZones instead of splitting mesh."
2053             << endl;
2054
2055         // Check if region overlaps with existing zone. If so keep.
2056
2057         for (label regionI = 0; regionI < nCellRegions; regionI++)
2058         {
2059             label zoneI = regionToZone[regionI];
2060
2061             if (zoneI != -1)
2062             {
2063                 Info<< "    Region " << regionI << " : corresponds to existing"
2064                     << " cellZone "
2065                     << zoneI << ' ' << cellZones[zoneI].name() << endl;
2066             }
2067             else
2068             {
2069                 // Create new cellZone.
2070                 labelList regionCells = findIndices(cellRegion, regionI);
2071
2072                 word zoneName = "region" + Foam::name(regionI);
2073
2074                 zoneI = cellZones.findZoneID(zoneName);
2075
2076                 if (zoneI == -1)
2077                 {
2078                     zoneI = cellZones.size();
2079                     mesh.cellZones().setSize(zoneI+1);
2080                     mesh.cellZones().set
2081                     (
2082                         zoneI,
2083                         new cellZone
2084                         (
2085                             zoneName,           //name
2086                             regionCells,        //addressing
2087                             zoneI,              //index
2088                             cellZones           //cellZoneMesh
2089                         )
2090                     );
2091                 }
2092                 else
2093                 {
2094                     mesh.cellZones()[zoneI].clearAddressing();
2095                     mesh.cellZones()[zoneI] = regionCells;
2096                 }
2097                 Info<< "    Region " << regionI << " : created new cellZone "
2098                     << zoneI << ' ' << cellZones[zoneI].name() << endl;
2099             }
2100         }
2101         mesh.cellZones().writeOpt() = IOobject::AUTO_WRITE;
2102
2103         if (!overwrite)
2104         {
2105             runTime++;
2106             mesh.setInstance(runTime.timeName());
2107         }
2108         else
2109         {
2110             mesh.setInstance(oldInstance);
2111         }
2112
2113         Info<< "Writing cellZones as new mesh to time " << runTime.timeName()
2114             << nl << endl;
2115
2116         mesh.write();
2117
2118
2119         // Write cellSets for convenience
2120         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2121
2122         Info<< "Writing cellSets corresponding to cellZones." << nl << endl;
2123
2124         forAll(cellZones, zoneI)
2125         {
2126             const cellZone& cz = cellZones[zoneI];
2127
2128             cellSet(mesh, cz.name(), cz).write();
2129         }
2130     }
2131     else
2132     {
2133         // Add patches for interfaces
2134         // ~~~~~~~~~~~~~~~~~~~~~~~~~~
2135
2136         // Add all possible patches. Empty ones get filtered later on.
2137         Info<< nl << "Adding patches" << nl << endl;
2138
2139         labelList interfacePatches
2140         (
2141             addRegionPatches
2142             (
2143                 mesh,
2144                 regionNames,
2145                 interfaces,
2146                 interfaceNames
2147             )
2148         );
2149
2150
2151         if (!overwrite)
2152         {
2153             runTime++;
2154         }
2155
2156
2157         // Create regions
2158         // ~~~~~~~~~~~~~~
2159
2160         if (insidePoint)
2161         {
2162             const point insidePoint = args.optionRead<point>("insidePoint");
2163
2164             label regionI = -1;
2165
2166             label cellI = mesh.findCell(insidePoint);
2167
2168             Info<< nl << "Found point " << insidePoint << " in cell " << cellI
2169                 << endl;
2170
2171             if (cellI != -1)
2172             {
2173                 regionI = cellRegion[cellI];
2174             }
2175
2176             reduce(regionI, maxOp<label>());
2177
2178             Info<< nl
2179                 << "Subsetting region " << regionI
2180                 << " containing point " << insidePoint << endl;
2181
2182             if (regionI == -1)
2183             {
2184                 FatalErrorIn(args.executable())
2185                     << "Point " << insidePoint
2186                     << " is not inside the mesh." << nl
2187                     << "Bounding box of the mesh:" << mesh.bounds()
2188                     << exit(FatalError);
2189             }
2190
2191             createAndWriteRegion
2192             (
2193                 mesh,
2194                 cellRegion,
2195                 regionNames,
2196                 faceToInterface,
2197                 interfacePatches,
2198                 regionI,
2199                 (overwrite ? oldInstance : runTime.timeName())
2200             );
2201         }
2202         else if (largestOnly)
2203         {
2204             label regionI = findMax(regionSizes);
2205
2206             Info<< nl
2207                 << "Subsetting region " << regionI
2208                 << " of size " << regionSizes[regionI] << endl;
2209
2210             createAndWriteRegion
2211             (
2212                 mesh,
2213                 cellRegion,
2214                 regionNames,
2215                 faceToInterface,
2216                 interfacePatches,
2217                 regionI,
2218                 (overwrite ? oldInstance : runTime.timeName())
2219             );
2220         }
2221         else
2222         {
2223             // Split all
2224             for (label regionI = 0; regionI < nCellRegions; regionI++)
2225             {
2226                 Info<< nl
2227                     << "Region " << regionI << nl
2228                     << "-------- " << endl;
2229
2230                 createAndWriteRegion
2231                 (
2232                     mesh,
2233                     cellRegion,
2234                     regionNames,
2235                     faceToInterface,
2236                     interfacePatches,
2237                     regionI,
2238                     (overwrite ? oldInstance : runTime.timeName())
2239                 );
2240             }
2241         }
2242     }
2243
2244     Info<< "End\n" << endl;
2245
2246     return 0;
2247 }
2248
2249
2250 // ************************************************************************* //