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

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