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