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