applications/utilities/mesh/advanced/PDRMesh/PDRMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Description
  26     Mesh and field preparation utility for PDR type simulations.
  27
  28     Reads
  29     - cellSet giving blockedCells
  30     - faceSets giving blockedFaces and the patch they should go into
  31
  32     NOTE: To avoid exposing wrong fields values faceSets should include
  33     faces contained in the blockedCells cellset.
  34
  35     - coupledFaces reads coupledFacesSet to introduces mixe-coupled baffles
  36
  37     Subsets out the blocked cells and splits the blockedFaces and updates
  38     fields.
  39
  40     The face splitting is done by duplicating the faces. No duplication of
  41     points for now (so checkMesh will show a lot of 'duplicate face' messages)
  42
  43 \*---------------------------------------------------------------------------*/
  44
  45 #include "fvMeshSubset.H"
  46 #include "argList.H"
  47 #include "cellSet.H"
  48 #include "IOobjectList.H"
  49 #include "volFields.H"
  50 #include "mapPolyMesh.H"
  51 #include "faceSet.H"
  52 #include "cellSet.H"
  53 #include "syncTools.H"
  54 #include "polyTopoChange.H"
  55 #include "polyModifyFace.H"
  56 #include "polyAddFace.H"
  57 #include "regionSplit.H"
  58 #include "Tuple2.H"
  59 #include "cyclicFvPatch.H"
  60
  61 using namespace Foam;
  62
  63 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  64
  65 void modifyOrAddFace
  66 (
  67     polyTopoChange& meshMod,
  68     const face& f,
  69     const label faceI,
  70     const label own,
  71     const bool flipFaceFlux,
  72     const label newPatchI,
  73     const label zoneID,
  74     const bool zoneFlip,
  75
  76     PackedBoolList& modifiedFace
  77 )
  78 {
  79     if (!modifiedFace[faceI])
  80     {
  81         // First usage of face. Modify.
  82         meshMod.setAction
  83         (
  84             polyModifyFace
  85             (
  86                 f,                          // modified face
  87                 faceI,                      // label of face
  88                 own,                        // owner
  89                 -1,                         // neighbour
  90                 flipFaceFlux,               // face flip
  91                 newPatchI,                  // patch for face
  92                 false,                      // remove from zone
  93                 zoneID,                     // zone for face
  94                 zoneFlip                    // face flip in zone
  95             )
  96         );
  97         modifiedFace[faceI] = 1;
  98     }
  99     else
 100     {
 101         // Second or more usage of face. Add.
 102         meshMod.setAction
 103         (
 104             polyAddFace
 105             (
 106                 f,                          // modified face
 107                 own,                        // owner
 108                 -1,                         // neighbour
 109                 -1,                         // master point
 110                 -1,                         // master edge
 111                 faceI,                      // master face
 112                 flipFaceFlux,               // face flip
 113                 newPatchI,                  // patch for face
 114                 zoneID,                     // zone for face
 115                 zoneFlip                    // face flip in zone
 116             )
 117         );
 118     }
 119 }
 120
 121
 122 template<class Type>
 123 void subsetVolFields
 124 (
 125     const fvMeshSubset& subsetter,
 126     const IOobjectList& objectsList,
 127     const label patchI,
 128     const Type& exposedValue,
 129     const word GeomVolType,
 130     PtrList<GeometricField<Type, fvPatchField, volMesh> >& subFields
 131 )
 132 {
 133     const fvMesh& baseMesh = subsetter.baseMesh();
 134
 135     label i = 0;
 136
 137     forAllConstIter(IOobjectList , objectsList, iter)
 138     {
 139         if (iter()->headerClassName() == GeomVolType)
 140         {
 141             const word fieldName = iter()->name();
 142
 143             Info<< "Subsetting field " << fieldName << endl;
 144
 145             GeometricField<Type, fvPatchField, volMesh> volField
 146             (
 147                 *iter(),
 148                 baseMesh
 149             );
 150
 151             subFields.set(i, subsetter.interpolate(volField));
 152
 153             // Explicitly set exposed faces (in patchI) to exposedValue.
 154             if (patchI >= 0)
 155             {
 156                 fvPatchField<Type>& fld =
 157                     subFields[i++].boundaryField()[patchI];
 158
 159                 label newStart = fld.patch().patch().start();
 160
 161                 label oldPatchI = subsetter.patchMap()[patchI];
 162
 163                 if (oldPatchI == -1)
 164                 {
 165                     // New patch. Reset whole value.
 166                     fld = exposedValue;
 167                 }
 168                 else
 169                 {
 170                     // Reset those faces that originate from different patch
 171                     // or internal faces.
 172                     label oldSize = volField.boundaryField()[oldPatchI].size();
 173                     label oldStart = volField.boundaryField()
 174                     [
 175                         oldPatchI
 176                     ].patch().patch().start();
 177
 178                     forAll(fld, j)
 179                     {
 180                         label oldFaceI = subsetter.faceMap()[newStart+j];
 181
 182                         if(oldFaceI < oldStart || oldFaceI >= oldStart+oldSize)
 183                         {
 184                             fld[j] = exposedValue;
 185                         }
 186                     }
 187                 }
 188             }
 189         }
 190     }
 191 }
 192
 193
 194 template<class Type>
 195 void subsetSurfaceFields
 196 (
 197     const fvMeshSubset& subsetter,
 198     const IOobjectList& objectsList,
 199     const label patchI,
 200     const Type& exposedValue,
 201     const word GeomSurfType,
 202     PtrList<GeometricField<Type, fvsPatchField, surfaceMesh> >& subFields
 203 )
 204 {
 205     const fvMesh& baseMesh = subsetter.baseMesh();
 206
 207     label i(0);
 208
 209     forAllConstIter(IOobjectList , objectsList, iter)
 210     {
 211         if (iter()->headerClassName() == GeomSurfType)
 212         {
 213             const word& fieldName = iter.key();
 214
 215             Info<< "Subsetting field " << fieldName << endl;
 216
 217             GeometricField<Type, fvsPatchField, surfaceMesh> volField
 218             (
 219                 *iter(),
 220                 baseMesh
 221             );
 222
 223             subFields.set(i, subsetter.interpolate(volField));
 224
 225
 226             // Explicitly set exposed faces (in patchI) to exposedValue.
 227             if (patchI >= 0)
 228             {
 229                 fvsPatchField<Type>& fld =
 230                     subFields[i++].boundaryField()[patchI];
 231
 232                 label newStart = fld.patch().patch().start();
 233
 234                 label oldPatchI = subsetter.patchMap()[patchI];
 235
 236                 if (oldPatchI == -1)
 237                 {
 238                     // New patch. Reset whole value.
 239                     fld = exposedValue;
 240                 }
 241                 else
 242                 {
 243                     // Reset those faces that originate from different patch
 244                     // or internal faces.
 245                     label oldSize = volField.boundaryField()[oldPatchI].size();
 246                     label oldStart = volField.boundaryField()
 247                     [
 248                         oldPatchI
 249                     ].patch().patch().start();
 250
 251                     forAll(fld, j)
 252                     {
 253                         label oldFaceI = subsetter.faceMap()[newStart+j];
 254
 255                         if(oldFaceI < oldStart || oldFaceI >= oldStart+oldSize)
 256                         {
 257                             fld[j] = exposedValue;
 258                         }
 259                     }
 260                 }
 261             }
 262         }
 263     }
 264 }
 265
 266
 267 // Faces introduced into zero-sized patches don't get a value at all.
 268 // This is hack to set them to an initial value.
 269 template<class GeoField>
 270 void initCreatedPatches
 271 (
 272     const fvMesh& mesh,
 273     const mapPolyMesh& map,
 274     const typename GeoField::value_type initValue
 275 )
 276 {
 277     HashTable<const GeoField*> fields
 278     (
 279         mesh.objectRegistry::lookupClass<GeoField>()
 280     );
 281
 282     for
 283     (
 284         typename HashTable<const GeoField*>::
 285             iterator fieldIter = fields.begin();
 286         fieldIter != fields.end();
 287         ++fieldIter
 288     )
 289     {
 290         GeoField& field = const_cast<GeoField&>(*fieldIter());
 291
 292         forAll(field.boundaryField(), patchi)
 293         {
 294             if (map.oldPatchSizes()[patchi] == 0)
 295             {
 296                 // Not mapped.
 297                 field.boundaryField()[patchi] = initValue;
 298
 299                 if (field.boundaryField()[patchi].fixesValue())
 300                 {
 301                     field.boundaryField()[patchi] == initValue;
 302                 }
 303             }
 304         }
 305     }
 306 }
 307
 308
 309 void createCoupledBaffles
 310 (
 311     fvMesh& mesh,
 312     const labelList& coupledWantedPatch,
 313     polyTopoChange& meshMod,
 314     PackedBoolList&  modifiedFace
 315 )
 316 {
 317     const faceZoneMesh& faceZones = mesh.faceZones();
 318
 319     forAll(coupledWantedPatch, faceI)
 320     {
 321         if (coupledWantedPatch[faceI] != -1)
 322         {
 323             const face& f = mesh.faces()[faceI];
 324             label zoneID = faceZones.whichZone(faceI);
 325             bool zoneFlip = false;
 326
 327             if (zoneID >= 0)
 328             {
 329                 const faceZone& fZone = faceZones[zoneID];
 330                 zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
 331             }
 332
 333             // Use owner side of face
 334             modifyOrAddFace
 335             (
 336                 meshMod,
 337                 f,                          // modified face
 338                 faceI,                      // label of face
 339                 mesh.faceOwner()[faceI],    // owner
 340                 false,                      // face flip
 341                 coupledWantedPatch[faceI],  // patch for face
 342                 zoneID,                     // zone for face
 343                 zoneFlip,                   // face flip in zone
 344                 modifiedFace                // modify or add status
 345             );
 346
 347             if (mesh.isInternalFace(faceI))
 348             {
 349                 label zoneID = faceZones.whichZone(faceI);
 350                 bool zoneFlip = false;
 351
 352                 if (zoneID >= 0)
 353                 {
 354                     const faceZone& fZone = faceZones[zoneID];
 355                     zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
 356                 }
 357                 // Use neighbour side of face
 358                 modifyOrAddFace
 359                 (
 360                     meshMod,
 361                     f.reverseFace(),            // modified face
 362                     faceI,                      // label of face
 363                     mesh.faceNeighbour()[faceI],// owner
 364                     false,                      // face flip
 365                     coupledWantedPatch[faceI],  // patch for face
 366                     zoneID,                     // zone for face
 367                     zoneFlip,                   // face flip in zone
 368                     modifiedFace                // modify or add status
 369                 );
 370             }
 371         }
 372     }
 373 }
 374
 375
 376 void createCyclicCoupledBaffles
 377 (
 378     fvMesh& mesh,
 379     const labelList& cyclicMasterPatch,
 380     const labelList& cyclicSlavePatch,
 381     polyTopoChange& meshMod,
 382     PackedBoolList&  modifiedFace
 383 )
 384 {
 385     const faceZoneMesh& faceZones = mesh.faceZones();
 386
 387     forAll(cyclicMasterPatch, faceI)
 388     {
 389         if (cyclicMasterPatch[faceI] != -1)
 390         {
 391             const face& f = mesh.faces()[faceI];
 392
 393             label zoneID = faceZones.whichZone(faceI);
 394             bool zoneFlip = false;
 395
 396             if (zoneID >= 0)
 397             {
 398                 const faceZone& fZone = faceZones[zoneID];
 399                 zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
 400             }
 401
 402             modifyOrAddFace
 403             (
 404                 meshMod,
 405                 f.reverseFace(),                    // modified face
 406                 faceI,                              // label of face
 407                 mesh.faceNeighbour()[faceI],        // owner
 408                 false,                              // face flip
 409                 cyclicMasterPatch[faceI],           // patch for face
 410                 zoneID,                             // zone for face
 411                 zoneFlip,                           // face flip in zone
 412                 modifiedFace                        // modify or add
 413             );
 414         }
 415     }
 416
 417     forAll(cyclicSlavePatch, faceI)
 418     {
 419         if (cyclicSlavePatch[faceI] != -1)
 420         {
 421             const face& f = mesh.faces()[faceI];
 422             if (mesh.isInternalFace(faceI))
 423             {
 424                 label zoneID = faceZones.whichZone(faceI);
 425                 bool zoneFlip = false;
 426
 427                 if (zoneID >= 0)
 428                 {
 429                     const faceZone& fZone = faceZones[zoneID];
 430                     zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
 431                 }
 432             // Use owner side of face
 433                 modifyOrAddFace
 434                 (
 435                     meshMod,
 436                     f,                          // modified face
 437                     faceI,                      // label of face
 438                     mesh.faceOwner()[faceI],    // owner
 439                     false,                      // face flip
 440                     cyclicSlavePatch[faceI],    // patch for face
 441                     zoneID,                     // zone for face
 442                     zoneFlip,                   // face flip in zone
 443                     modifiedFace                // modify or add status
 444                 );
 445             }
 446         }
 447     }
 448 }
 449
 450
 451 void createBaffles
 452 (
 453     fvMesh& mesh,
 454     const labelList& wantedPatch,
 455     polyTopoChange& meshMod
 456 )
 457 {
 458     const faceZoneMesh& faceZones = mesh.faceZones();
 459     Info << "faceZone:createBaffle " << faceZones << endl;
 460     forAll(wantedPatch, faceI)
 461     {
 462         if (wantedPatch[faceI] != -1)
 463         {
 464             const face& f = mesh.faces()[faceI];
 465
 466             label zoneID = faceZones.whichZone(faceI);
 467             bool zoneFlip = false;
 468
 469             if (zoneID >= 0)
 470             {
 471                 const faceZone& fZone = faceZones[zoneID];
 472                 zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
 473             }
 474
 475             meshMod.setAction
 476             (
 477                 polyModifyFace
 478                 (
 479                     f,                          // modified face
 480                     faceI,                      // label of face
 481                     mesh.faceOwner()[faceI],    // owner
 482                     -1,                         // neighbour
 483                     false,                      // face flip
 484                     wantedPatch[faceI],         // patch for face
 485                     false,                      // remove from zone
 486                     zoneID,                     // zone for face
 487                     zoneFlip                    // face flip in zone
 488                 )
 489             );
 490
 491             if (mesh.isInternalFace(faceI))
 492             {
 493                 label zoneID = faceZones.whichZone(faceI);
 494                 bool zoneFlip = false;
 495
 496                 if (zoneID >= 0)
 497                 {
 498                     const faceZone& fZone = faceZones[zoneID];
 499                     zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
 500                 }
 501
 502                 meshMod.setAction
 503                 (
 504                     polyAddFace
 505                     (
 506                         f.reverseFace(),            // modified face
 507                         mesh.faceNeighbour()[faceI],// owner
 508                         -1,                         // neighbour
 509                         -1,                         // masterPointID
 510                         -1,                         // masterEdgeID
 511                         faceI,                      // masterFaceID,
 512                         false,                      // face flip
 513                         wantedPatch[faceI],         // patch for face
 514                         zoneID,                     // zone for face
 515                         zoneFlip                    // face flip in zone
 516                     )
 517                 );
 518             }
 519         }
 520     }
 521 }
 522
 523
 524 // Wrapper around find patch. Also makes sure same patch in parallel.
 525 label findPatch(const polyBoundaryMesh& patches, const word& patchName)
 526 {
 527     label patchI = patches.findPatchID(patchName);
 528
 529     if (patchI == -1)
 530     {
 531         FatalErrorIn("findPatch(const polyBoundaryMesh&, const word&)")
 532             << "Illegal patch " << patchName
 533             << nl << "Valid patches are " << patches.names()
 534             << exit(FatalError);
 535     }
 536
 537     // Check same patch for all procs
 538     {
 539         label newPatch = patchI;
 540         reduce(newPatch, minOp<label>());
 541
 542         if (newPatch != patchI)
 543         {
 544             FatalErrorIn("findPatch(const polyBoundaryMesh&, const word&)")
 545                 << "Patch " << patchName
 546                 << " should have the same patch index on all processors." << nl
 547                 << "On my processor it has index " << patchI
 548                 << " ; on some other processor it has index " << newPatch
 549                 << exit(FatalError);
 550         }
 551     }
 552     return patchI;
 553 }
 554
 555
 556 // Main program:
 557
 558 int main(int argc, char *argv[])
 559 {
 560     #include "addOverwriteOption.H"
 561     #include "setRootCase.H"
 562     #include "createTime.H"
 563     runTime.functionObjects().off();
 564     #include "createMesh.H"
 565
 566     // Read control dictionary
 567     // ~~~~~~~~~~~~~~~~~~~~~~~
 568
 569     IOdictionary dict
 570     (
 571         IOobject
 572         (
 573             "PDRMeshDict",
 574             runTime.system(),
 575             mesh,
 576             IOobject::MUST_READ_IF_MODIFIED,
 577             IOobject::NO_WRITE
 578         )
 579     );
 580
 581     // Per faceSet the patch to put the baffles into
 582     const List<Pair<word> > setsAndPatches(dict.lookup("blockedFaces"));
 583
 584     // Per faceSet the patch to put the coupled baffles into
 585     DynamicList<FixedList<word, 3> > coupledAndPatches(10);
 586     const dictionary& functionDicts = dict.subDict("coupledFaces");
 587     forAllConstIter(dictionary, functionDicts, iter)
 588     {
 589         // safety:
 590         if (!iter().isDict())
 591         {
 592             continue;
 593         }
 594         const word& key = iter().keyword();
 595
 596         const dictionary& dict = iter().dict();
 597         const word cyclicName = dict.lookup("cyclicMasterPatchName");
 598         const word wallName = dict.lookup("wallPatchName");
 599         FixedList<word, 3> nameAndType;
 600         nameAndType[0] = key;
 601         nameAndType[1] = wallName;
 602         nameAndType[2] = cyclicName;
 603         coupledAndPatches.append(nameAndType);
 604     }
 605
 606     forAll(setsAndPatches, setI)
 607     {
 608         Info<< "Faces in faceSet " << setsAndPatches[setI][0]
 609             << " become baffles in patch " << setsAndPatches[setI][1]
 610             << endl;
 611     }
 612
 613     forAll(coupledAndPatches, setI)
 614     {
 615         Info<< "Faces in faceSet " << coupledAndPatches[setI][0]
 616             << " become coupled baffles in patch " << coupledAndPatches[setI][1]
 617             << endl;
 618     }
 619
 620     // All exposed faces that are not explicitly marked to be put into a patch
 621     const word defaultPatch(dict.lookup("defaultPatch"));
 622
 623     Info<< "Faces that get exposed become boundary faces in patch "
 624         << defaultPatch << endl;
 625
 626     const word blockedSetName(dict.lookup("blockedCells"));
 627
 628     Info<< "Reading blocked cells from cellSet " << blockedSetName
 629         << endl;
 630
 631     const bool overwrite = args.optionFound("overwrite");
 632
 633
 634     // Read faceSets, lookup patches
 635     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 636
 637     // Check that face sets don't have coincident faces
 638     labelList wantedPatch(mesh.nFaces(), -1);
 639     forAll(setsAndPatches, setI)
 640     {
 641         faceSet fSet(mesh, setsAndPatches[setI][0]);
 642
 643         label patchI = findPatch
 644         (
 645             mesh.boundaryMesh(),
 646             setsAndPatches[setI][1]
 647         );
 648
 649         forAllConstIter(faceSet, fSet, iter)
 650         {
 651             if (wantedPatch[iter.key()] != -1)
 652             {
 653                 FatalErrorIn(args.executable())
 654                     << "Face " << iter.key()
 655                     << " is in faceSet " << setsAndPatches[setI][0]
 656                     << " destined for patch " << setsAndPatches[setI][1]
 657                     << " but also in patch " << wantedPatch[iter.key()]
 658                     << exit(FatalError);
 659             }
 660             wantedPatch[iter.key()] = patchI;
 661         }
 662     }
 663
 664     // Per face the patch for coupled baffle or -1.
 665     labelList coupledWantedPatch(mesh.nFaces(), -1);
 666     labelList cyclicWantedPatch_half0(mesh.nFaces(), -1);
 667     labelList cyclicWantedPatch_half1(mesh.nFaces(), -1);
 668
 669     forAll(coupledAndPatches, setI)
 670     {
 671         const polyBoundaryMesh& patches = mesh.boundaryMesh();
 672         const label cyclicId =
 673             findPatch(patches, coupledAndPatches[setI][2]);
 674
 675         const label cyclicSlaveId = findPatch
 676         (
 677             patches,
 678             refCast<const cyclicFvPatch>
 679             (
 680                 mesh.boundary()[cyclicId]
 681             ).neighbFvPatch().name()
 682         );
 683
 684         faceSet fSet(mesh, coupledAndPatches[setI][0]);
 685         label patchI = findPatch(patches, coupledAndPatches[setI][1]);
 686
 687         forAllConstIter(faceSet, fSet, iter)
 688         {
 689             if (coupledWantedPatch[iter.key()] != -1)
 690             {
 691                 FatalErrorIn(args.executable())
 692                     << "Face " << iter.key()
 693                     << " is in faceSet " << coupledAndPatches[setI][0]
 694                     << " destined for patch " << coupledAndPatches[setI][1]
 695                     << " but also in patch " << coupledWantedPatch[iter.key()]
 696                     << exit(FatalError);
 697             }
 698                 coupledWantedPatch[iter.key()] = patchI;
 699                 cyclicWantedPatch_half0[iter.key()] = cyclicId;
 700                 cyclicWantedPatch_half1[iter.key()] = cyclicSlaveId;
 701         }
 702     }
 703
 704     // Exposed faces patch
 705     label defaultPatchI = findPatch(mesh.boundaryMesh(), defaultPatch);
 706
 707
 708     //
 709     // Removing blockedCells (like subsetMesh)
 710     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 711     //
 712
 713     // Create mesh subsetting engine
 714     fvMeshSubset subsetter(mesh);
 715
 716     {
 717
 718         cellSet blockedCells(mesh, blockedSetName);
 719
 720         // invert
 721         blockedCells.invert(mesh.nCells());
 722
 723         // Create subsetted mesh.
 724         subsetter.setLargeCellSubset(blockedCells, defaultPatchI, true);
 725     }
 726
 727
 728     // Subset wantedPatch. Note that might also include boundary faces
 729     // that have been exposed by subsetting.
 730     wantedPatch = IndirectList<label>(wantedPatch, subsetter.faceMap())();
 731
 732     coupledWantedPatch = IndirectList<label>
 733     (
 734         coupledWantedPatch,
 735         subsetter.faceMap()
 736     )();
 737
 738     cyclicWantedPatch_half0 = IndirectList<label>
 739     (
 740         cyclicWantedPatch_half0,
 741         subsetter.faceMap()
 742     )();
 743
 744     cyclicWantedPatch_half1 = IndirectList<label>
 745     (
 746         cyclicWantedPatch_half1,
 747         subsetter.faceMap()
 748     )();
 749
 750     // Read all fields in time and constant folders
 751     IOobjectList objects(mesh, runTime.timeName());
 752     IOobjectList timeObjects(IOobjectList(mesh, mesh.facesInstance()));
 753     forAllConstIter(IOobjectList, timeObjects, iter)
 754     {
 755         if
 756         (
 757             iter()->headerClassName() == volScalarField::typeName
 758          || iter()->headerClassName() == volVectorField::typeName
 759          || iter()->headerClassName() == volSphericalTensorField::typeName
 760          || iter()->headerClassName() == volTensorField::typeName
 761          || iter()->headerClassName() == volSymmTensorField::typeName
 762          || iter()->headerClassName() == surfaceScalarField::typeName
 763          || iter()->headerClassName() == surfaceVectorField::typeName
 764          || iter()->headerClassName()
 765             == surfaceSphericalTensorField::typeName
 766          || iter()->headerClassName() == surfaceSymmTensorField::typeName
 767          || iter()->headerClassName() == surfaceTensorField::typeName
 768         )
 769         {
 770             objects.add(*iter());
 771         }
 772     }
 773
 774     // Read vol fields and subset.
 775
 776     wordList scalarNames(objects.names(volScalarField::typeName));
 777     PtrList<volScalarField> scalarFlds(scalarNames.size());
 778     subsetVolFields
 779     (
 780         subsetter,
 781         objects,
 782         defaultPatchI,
 783         pTraits<scalar>::zero,
 784         volScalarField::typeName,
 785         scalarFlds
 786     );
 787
 788     wordList vectorNames(objects.names(volVectorField::typeName));
 789     PtrList<volVectorField> vectorFlds(vectorNames.size());
 790     subsetVolFields
 791     (
 792         subsetter,
 793         objects,
 794         defaultPatchI,
 795         pTraits<vector>::zero,
 796         volVectorField::typeName,
 797         vectorFlds
 798     );
 799
 800     wordList sphericalTensorNames
 801     (
 802         objects.names(volSphericalTensorField::typeName)
 803     );
 804     PtrList<volSphericalTensorField> sphericalTensorFlds
 805     (
 806         sphericalTensorNames.size()
 807     );
 808     subsetVolFields
 809     (
 810         subsetter,
 811         objects,
 812         defaultPatchI,
 813         pTraits<sphericalTensor>::zero,
 814         volSphericalTensorField::typeName,
 815         sphericalTensorFlds
 816     );
 817
 818     wordList symmTensorNames(objects.names(volSymmTensorField::typeName));
 819     PtrList<volSymmTensorField> symmTensorFlds(symmTensorNames.size());
 820     subsetVolFields
 821     (
 822         subsetter,
 823         objects,
 824         defaultPatchI,
 825         pTraits<symmTensor>::zero,
 826         volSymmTensorField::typeName,
 827         symmTensorFlds
 828     );
 829
 830     wordList tensorNames(objects.names(volTensorField::typeName));
 831     PtrList<volTensorField> tensorFlds(tensorNames.size());
 832     subsetVolFields
 833     (
 834         subsetter,
 835         objects,
 836         defaultPatchI,
 837         pTraits<tensor>::zero,
 838         volTensorField::typeName,
 839         tensorFlds
 840     );
 841
 842     // Read surface fields and subset.
 843
 844     wordList surfScalarNames(objects.names(surfaceScalarField::typeName));
 845     PtrList<surfaceScalarField> surfScalarFlds(surfScalarNames.size());
 846     subsetSurfaceFields
 847     (
 848         subsetter,
 849         objects,
 850         defaultPatchI,
 851         pTraits<scalar>::zero,
 852         surfaceScalarField::typeName,
 853         surfScalarFlds
 854     );
 855
 856     wordList surfVectorNames(objects.names(surfaceVectorField::typeName));
 857     PtrList<surfaceVectorField> surfVectorFlds(surfVectorNames.size());
 858     subsetSurfaceFields
 859     (
 860         subsetter,
 861         objects,
 862         defaultPatchI,
 863         pTraits<vector>::zero,
 864         surfaceVectorField::typeName,
 865         surfVectorFlds
 866     );
 867
 868     wordList surfSphericalTensorNames
 869     (
 870         objects.names(surfaceSphericalTensorField::typeName)
 871     );
 872     PtrList<surfaceSphericalTensorField> surfSphericalTensorFlds
 873     (
 874         surfSphericalTensorNames.size()
 875     );
 876     subsetSurfaceFields
 877     (
 878         subsetter,
 879         objects,
 880         defaultPatchI,
 881         pTraits<sphericalTensor>::zero,
 882         surfaceSphericalTensorField::typeName,
 883         surfSphericalTensorFlds
 884     );
 885
 886     wordList surfSymmTensorNames
 887     (
 888         objects.names(surfaceSymmTensorField::typeName)
 889     );
 890
 891     PtrList<surfaceSymmTensorField> surfSymmTensorFlds
 892     (
 893         surfSymmTensorNames.size()
 894     );
 895
 896     subsetSurfaceFields
 897     (
 898         subsetter,
 899         objects,
 900         defaultPatchI,
 901         pTraits<symmTensor>::zero,
 902         surfaceSymmTensorField::typeName,
 903         surfSymmTensorFlds
 904     );
 905
 906     wordList surfTensorNames(objects.names(surfaceTensorField::typeName));
 907     PtrList<surfaceTensorField> surfTensorFlds(surfTensorNames.size());
 908     subsetSurfaceFields
 909     (
 910         subsetter,
 911         objects,
 912         defaultPatchI,
 913         pTraits<tensor>::zero,
 914         surfaceTensorField::typeName,
 915         surfTensorFlds
 916     );
 917
 918     if (!overwrite)
 919     {
 920         runTime++;
 921     }
 922
 923     Info<< "Writing mesh without blockedCells to time " << runTime.value()
 924         << endl;
 925
 926     // Subsetting adds 'subset' prefix. Rename field to be like original.
 927     forAll(scalarFlds, i)
 928     {
 929         scalarFlds[i].rename(scalarNames[i]);
 930         scalarFlds[i].writeOpt() = IOobject::AUTO_WRITE;
 931     }
 932     forAll(vectorFlds, i)
 933     {
 934         vectorFlds[i].rename(vectorNames[i]);
 935         vectorFlds[i].writeOpt() = IOobject::AUTO_WRITE;
 936     }
 937     forAll(sphericalTensorFlds, i)
 938     {
 939         sphericalTensorFlds[i].rename(sphericalTensorNames[i]);
 940         sphericalTensorFlds[i].writeOpt() = IOobject::AUTO_WRITE;
 941     }
 942     forAll(symmTensorFlds, i)
 943     {
 944         symmTensorFlds[i].rename(symmTensorNames[i]);
 945         symmTensorFlds[i].writeOpt() = IOobject::AUTO_WRITE;
 946     }
 947     forAll(tensorFlds, i)
 948     {
 949         tensorFlds[i].rename(tensorNames[i]);
 950         tensorFlds[i].writeOpt() = IOobject::AUTO_WRITE;
 951     }
 952
 953     // Surface ones.
 954     forAll(surfScalarFlds, i)
 955     {
 956         surfScalarFlds[i].rename(surfScalarNames[i]);
 957         surfScalarFlds[i].writeOpt() = IOobject::AUTO_WRITE;
 958     }
 959     forAll(surfVectorFlds, i)
 960     {
 961         surfVectorFlds[i].rename(surfVectorNames[i]);
 962         surfVectorFlds[i].writeOpt() = IOobject::AUTO_WRITE;
 963     }
 964     forAll(surfSphericalTensorFlds, i)
 965     {
 966         surfSphericalTensorFlds[i].rename(surfSphericalTensorNames[i]);
 967         surfSphericalTensorFlds[i].writeOpt() = IOobject::AUTO_WRITE;
 968     }
 969     forAll(surfSymmTensorFlds, i)
 970     {
 971         surfSymmTensorFlds[i].rename(surfSymmTensorNames[i]);
 972         surfSymmTensorFlds[i].writeOpt() = IOobject::AUTO_WRITE;
 973     }
 974     forAll(surfTensorNames, i)
 975     {
 976         surfTensorFlds[i].rename(surfTensorNames[i]);
 977         surfTensorFlds[i].writeOpt() = IOobject::AUTO_WRITE;
 978     }
 979
 980     subsetter.subMesh().write();
 981
 982
 983     //
 984     // Splitting blockedFaces
 985     // ~~~~~~~~~~~~~~~~~~~~~~
 986     //
 987
 988     // Synchronize wantedPatch across coupled patches.
 989     syncTools::syncFaceList
 990     (
 991         subsetter.subMesh(),
 992         wantedPatch,
 993         maxEqOp<label>()
 994     );
 995
 996     // Synchronize coupledWantedPatch across coupled patches.
 997     syncTools::syncFaceList
 998     (
 999         subsetter.subMesh(),
1000         coupledWantedPatch,
1001         maxEqOp<label>()
1002     );
1003
1004     // Synchronize cyclicWantedPatch across coupled patches.
1005     syncTools::syncFaceList
1006     (
1007         subsetter.subMesh(),
1008         cyclicWantedPatch_half0,
1009         maxEqOp<label>()
1010     );
1011
1012     // Synchronize cyclicWantedPatch across coupled patches.
1013     syncTools::syncFaceList
1014     (
1015         subsetter.subMesh(),
1016         cyclicWantedPatch_half1,
1017         maxEqOp<label>()
1018     );
1019
1020     // Topochange container
1021     polyTopoChange meshMod(subsetter.subMesh());
1022
1023
1024     // Whether first use of face (modify) or consecutive (add)
1025     PackedBoolList modifiedFace(mesh.nFaces());
1026
1027     // Create coupled wall-side baffles
1028     createCoupledBaffles
1029     (
1030         subsetter.subMesh(),
1031         coupledWantedPatch,
1032         meshMod,
1033         modifiedFace
1034     );
1035
1036     // Create coupled master/slave cyclic baffles
1037     createCyclicCoupledBaffles
1038     (
1039         subsetter.subMesh(),
1040         cyclicWantedPatch_half0,
1041         cyclicWantedPatch_half1,
1042         meshMod,
1043         modifiedFace
1044     );
1045
1046     // Create wall baffles
1047     createBaffles
1048     (
1049         subsetter.subMesh(),
1050         wantedPatch,
1051         meshMod
1052     );
1053
1054     if (!overwrite)
1055     {
1056         runTime++;
1057     }
1058
1059     // Change the mesh. Change points directly (no inflation).
1060     autoPtr<mapPolyMesh> map = meshMod.changeMesh(subsetter.subMesh(), false);
1061
1062     // Update fields
1063     subsetter.subMesh().updateMesh(map);
1064
1065     // Fix faces that get mapped to zero-sized patches (these don't get any
1066     // value)
1067     initCreatedPatches<volScalarField>
1068     (
1069         subsetter.subMesh(),
1070         map,
1071         0.0
1072     );
1073     initCreatedPatches<volVectorField>
1074     (
1075         subsetter.subMesh(),
1076         map,
1077         vector::zero
1078     );
1079     initCreatedPatches<volSphericalTensorField>
1080     (
1081         subsetter.subMesh(),
1082         map,
1083         sphericalTensor::zero
1084     );
1085     initCreatedPatches<volSymmTensorField>
1086     (
1087         subsetter.subMesh(),
1088         map,
1089         symmTensor::zero
1090     );
1091     initCreatedPatches<volTensorField>
1092     (
1093         subsetter.subMesh(),
1094         map,
1095         tensor::zero
1096     );
1097
1098     initCreatedPatches<surfaceScalarField>
1099     (
1100         subsetter.subMesh(),
1101         map,
1102         0.0
1103     );
1104     initCreatedPatches<surfaceVectorField>
1105     (
1106         subsetter.subMesh(),
1107         map,
1108         vector::zero
1109     );
1110     initCreatedPatches<surfaceSphericalTensorField>
1111     (
1112         subsetter.subMesh(),
1113         map,
1114         sphericalTensor::zero
1115     );
1116     initCreatedPatches<surfaceSymmTensorField>
1117     (
1118         subsetter.subMesh(),
1119         map,
1120         symmTensor::zero
1121     );
1122     initCreatedPatches<surfaceTensorField>
1123     (
1124         subsetter.subMesh(),
1125         map,
1126         tensor::zero
1127     );
1128
1129
1130     // Move mesh (since morphing might not do this)
1131     if (map().hasMotionPoints())
1132     {
1133         subsetter.subMesh().movePoints(map().preMotionPoints());
1134     }
1135
1136     Info<< "Writing mesh with split blockedFaces to time " << runTime.value()
1137         << endl;
1138
1139     subsetter.subMesh().write();
1140
1141
1142     //
1143     // Removing inaccessible regions
1144     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1145     //
1146
1147     // Determine connected regions. regionSplit is the labelList with the
1148     // region per cell.
1149     regionSplit cellRegion(subsetter.subMesh());
1150
1151     if (cellRegion.nRegions() > 1)
1152     {
1153         WarningIn(args.executable())
1154             << "Removing blocked faces and cells created "
1155             << cellRegion.nRegions()
1156             << " regions that are not connected via a face." << nl
1157             << "    This is not supported in solvers." << nl
1158             << "    Use" << nl << nl
1159             << "    splitMeshRegions <root> <case> -largestOnly" << nl << nl
1160             << "    to extract a single region of the mesh." << nl
1161             << "    This mesh will be written to a new timedirectory"
1162             << " so might have to be moved back to constant/" << nl
1163             << endl;
1164
1165         word startFrom(runTime.controlDict().lookup("startFrom"));
1166
1167         if (startFrom != "latestTime")
1168         {
1169             WarningIn(args.executable())
1170                 << "To run splitMeshRegions please set your"
1171                 << " startFrom entry to latestTime" << endl;
1172         }
1173     }
1174
1175     Info << nl << "End" << endl;
1176
1177     return 0;
1178 }
1179
1180
1181 // ************************************************************************* //