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