src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*----------------------------------------------------------------------------*/
  26
  27 #include "meshRefinement.H"
  28 #include "fvMesh.H"
  29 #include "syncTools.H"
  30 #include "Time.H"
  31 #include "refinementSurfaces.H"
  32 #include "pointSet.H"
  33 #include "faceSet.H"
  34 #include "indirectPrimitivePatch.H"
  35 #include "polyTopoChange.H"
  36 #include "meshTools.H"
  37 #include "polyModifyFace.H"
  38 #include "polyModifyCell.H"
  39 #include "polyAddFace.H"
  40 #include "polyRemoveFace.H"
  41 #include "polyAddPoint.H"
  42 #include "localPointRegion.H"
  43 #include "duplicatePoints.H"
  44 #include "OFstream.H"
  45 #include "regionSplit.H"
  46 #include "removeCells.H"
  47
  48 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  49
  50 // Repatches external face or creates baffle for internal face
  51 // with user specified patches (might be different for both sides).
  52 // Returns label of added face.
  53 Foam::label Foam::meshRefinement::createBaffle
  54 (
  55     const label faceI,
  56     const label ownPatch,
  57     const label neiPatch,
  58     polyTopoChange& meshMod
  59 ) const
  60 {
  61     const face& f = mesh_.faces()[faceI];
  62     label zoneID = mesh_.faceZones().whichZone(faceI);
  63     bool zoneFlip = false;
  64
  65     if (zoneID >= 0)
  66     {
  67         const faceZone& fZone = mesh_.faceZones()[zoneID];
  68         zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
  69     }
  70
  71     meshMod.setAction
  72     (
  73         polyModifyFace
  74         (
  75             f,                          // modified face
  76             faceI,                      // label of face
  77             mesh_.faceOwner()[faceI],   // owner
  78             -1,                         // neighbour
  79             false,                      // face flip
  80             ownPatch,                   // patch for face
  81             false,                      // remove from zone
  82             zoneID,                     // zone for face
  83             zoneFlip                    // face flip in zone
  84         )
  85     );
  86
  87
  88     label dupFaceI = -1;
  89
  90     if (mesh_.isInternalFace(faceI))
  91     {
  92         if (neiPatch == -1)
  93         {
  94             FatalErrorIn
  95             (
  96                 "meshRefinement::createBaffle"
  97                 "(const label, const label, const label, polyTopoChange&)"
  98             )   << "No neighbour patch for internal face " << faceI
  99                 << " fc:" << mesh_.faceCentres()[faceI]
 100                 << " ownPatch:" << ownPatch << abort(FatalError);
 101         }
 102
 103         bool reverseFlip = false;
 104         if (zoneID >= 0)
 105         {
 106             reverseFlip = !zoneFlip;
 107         }
 108
 109         dupFaceI = meshMod.setAction
 110         (
 111             polyAddFace
 112             (
 113                 f.reverseFace(),            // modified face
 114                 mesh_.faceNeighbour()[faceI],// owner
 115                 -1,                         // neighbour
 116                 -1,                         // masterPointID
 117                 -1,                         // masterEdgeID
 118                 faceI,                      // masterFaceID,
 119                 true,                       // face flip
 120                 neiPatch,                   // patch for face
 121                 zoneID,                     // zone for face
 122                 reverseFlip                 // face flip in zone
 123             )
 124         );
 125     }
 126     return dupFaceI;
 127 }
 128
 129
 130 // Get an estimate for the patch the internal face should get. Bit heuristic.
 131 Foam::label Foam::meshRefinement::getBafflePatch
 132 (
 133     const labelList& facePatch,
 134     const label faceI
 135 ) const
 136 {
 137     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 138
 139     // Loop over face points
 140     // for each point check all faces patch IDs
 141     // as soon as an ID >= 0 is found, break and assign that ID
 142     // to the current face.
 143     // Check first for real patch (so proper surface intersection and then
 144     // in facePatch array for patches to block off faces
 145
 146     forAll(mesh_.faces()[faceI], fp)
 147     {
 148         label pointI = mesh_.faces()[faceI][fp];
 149
 150         forAll(mesh_.pointFaces()[pointI], pf)
 151         {
 152             label pFaceI = mesh_.pointFaces()[pointI][pf];
 153
 154             label patchI = patches.whichPatch(pFaceI);
 155
 156             if (patchI != -1 && !patches[patchI].coupled())
 157             {
 158                 return patchI;
 159             }
 160             else if (facePatch[pFaceI] != -1)
 161             {
 162                 return facePatch[pFaceI];
 163             }
 164         }
 165     }
 166
 167     // Loop over owner and neighbour cells, looking for the first face with a
 168     // valid patch number
 169     const cell& ownFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
 170
 171     forAll(ownFaces, i)
 172     {
 173         label cFaceI = ownFaces[i];
 174
 175         label patchI = patches.whichPatch(cFaceI);
 176
 177         if (patchI != -1 && !patches[patchI].coupled())
 178         {
 179             return patchI;
 180         }
 181         else if (facePatch[cFaceI] != -1)
 182         {
 183             return facePatch[cFaceI];
 184         }
 185     }
 186
 187     if (mesh_.isInternalFace(faceI))
 188     {
 189         const cell& neiFaces = mesh_.cells()[mesh_.faceNeighbour()[faceI]];
 190
 191         forAll(neiFaces, i)
 192         {
 193             label cFaceI = neiFaces[i];
 194
 195             label patchI = patches.whichPatch(cFaceI);
 196
 197             if (patchI != -1 && !patches[patchI].coupled())
 198             {
 199                 return patchI;
 200             }
 201             else if (facePatch[cFaceI] != -1)
 202             {
 203                 return facePatch[cFaceI];
 204             }
 205         }
 206     }
 207
 208     WarningIn
 209     (
 210         "meshRefinement::getBafflePatch(const labelList&, const label)"
 211     )   << "Could not find boundary face neighbouring internal face "
 212         << faceI << " with face centre " << mesh_.faceCentres()[faceI]
 213         << nl
 214         << "Using arbitrary patch " << 0 << " instead." << endl;
 215
 216     return 0;
 217 }
 218
 219
 220 // Determine patches for baffles.
 221 void Foam::meshRefinement::getBafflePatches
 222 (
 223     const labelList& globalToPatch,
 224     const labelList& neiLevel,
 225     const pointField& neiCc,
 226
 227     labelList& ownPatch,
 228     labelList& neiPatch
 229 ) const
 230 {
 231     autoPtr<OFstream> str;
 232     label vertI = 0;
 233     if (debug&OBJINTERSECTIONS)
 234     {
 235         str.reset
 236         (
 237             new OFstream
 238             (
 239                 mesh_.time().path()/timeName()/"intersections.obj"
 240             )
 241         );
 242
 243         Pout<< "getBafflePatches : Writing surface intersections to file "
 244             << str().name() << nl << endl;
 245     }
 246
 247     const pointField& cellCentres = mesh_.cellCentres();
 248
 249     // Build list of surfaces that are not to be baffled.
 250     const wordList& cellZoneNames = surfaces_.cellZoneNames();
 251
 252     labelList surfacesToBaffle(cellZoneNames.size());
 253     label baffleI = 0;
 254     forAll(cellZoneNames, surfI)
 255     {
 256         if (cellZoneNames[surfI].size())
 257         {
 258             if (debug)
 259             {
 260                 Pout<< "getBafflePatches : Not baffling surface "
 261                     << surfaces_.names()[surfI] << endl;
 262             }
 263         }
 264         else
 265         {
 266             surfacesToBaffle[baffleI++] = surfI;
 267         }
 268     }
 269     surfacesToBaffle.setSize(baffleI);
 270
 271     ownPatch.setSize(mesh_.nFaces());
 272     ownPatch = -1;
 273     neiPatch.setSize(mesh_.nFaces());
 274     neiPatch = -1;
 275
 276
 277     // Collect candidate faces
 278     // ~~~~~~~~~~~~~~~~~~~~~~~
 279
 280     labelList testFaces(intersectedFaces());
 281
 282     // Collect segments
 283     // ~~~~~~~~~~~~~~~~
 284
 285     pointField start(testFaces.size());
 286     pointField end(testFaces.size());
 287
 288     forAll(testFaces, i)
 289     {
 290         label faceI = testFaces[i];
 291
 292         label own = mesh_.faceOwner()[faceI];
 293
 294         if (mesh_.isInternalFace(faceI))
 295         {
 296             start[i] = cellCentres[own];
 297             end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
 298         }
 299         else
 300         {
 301             start[i] = cellCentres[own];
 302             end[i] = neiCc[faceI-mesh_.nInternalFaces()];
 303         }
 304     }
 305
 306
 307     // Do test for intersections
 308     // ~~~~~~~~~~~~~~~~~~~~~~~~~
 309     labelList surface1;
 310     List<pointIndexHit> hit1;
 311     labelList region1;
 312     labelList surface2;
 313     List<pointIndexHit> hit2;
 314     labelList region2;
 315     surfaces_.findNearestIntersection
 316     (
 317         surfacesToBaffle,
 318         start,
 319         end,
 320
 321         surface1,
 322         hit1,
 323         region1,
 324
 325         surface2,
 326         hit2,
 327         region2
 328     );
 329
 330     forAll(testFaces, i)
 331     {
 332         label faceI = testFaces[i];
 333
 334         if (hit1[i].hit() && hit2[i].hit())
 335         {
 336             if (str.valid())
 337             {
 338                 meshTools::writeOBJ(str(), start[i]);
 339                 vertI++;
 340                 meshTools::writeOBJ(str(), hit1[i].rawPoint());
 341                 vertI++;
 342                 meshTools::writeOBJ(str(), hit2[i].rawPoint());
 343                 vertI++;
 344                 meshTools::writeOBJ(str(), end[i]);
 345                 vertI++;
 346                 str()<< "l " << vertI-3 << ' ' << vertI-2 << nl;
 347                 str()<< "l " << vertI-2 << ' ' << vertI-1 << nl;
 348                 str()<< "l " << vertI-1 << ' ' << vertI << nl;
 349             }
 350
 351             // Pick up the patches
 352             ownPatch[faceI] = globalToPatch
 353             [
 354                 surfaces_.globalRegion(surface1[i], region1[i])
 355             ];
 356             neiPatch[faceI] = globalToPatch
 357             [
 358                 surfaces_.globalRegion(surface2[i], region2[i])
 359             ];
 360
 361             if (ownPatch[faceI] == -1 || neiPatch[faceI] == -1)
 362             {
 363                 FatalErrorIn("getBafflePatches(..)")
 364                     << "problem." << abort(FatalError);
 365             }
 366         }
 367     }
 368
 369     // No need to parallel sync since intersection data (surfaceIndex_ etc.)
 370     // already guaranteed to be synced...
 371     // However:
 372     // - owncc and neicc are reversed on different procs so might pick
 373     //   up different regions reversed? No problem. Neighbour on one processor
 374     //   might not be owner on the other processor but the neighbour is
 375     //   not used when creating baffles from proc faces.
 376     // - tolerances issues occasionally crop up.
 377     syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>(), false);
 378     syncTools::syncFaceList(mesh_, neiPatch, maxEqOp<label>(), false);
 379 }
 380
 381
 382 // Create baffle for every face where ownPatch != -1
 383 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
 384 (
 385     const labelList& ownPatch,
 386     const labelList& neiPatch
 387 )
 388 {
 389     if
 390     (
 391         ownPatch.size() != mesh_.nFaces()
 392      || neiPatch.size() != mesh_.nFaces()
 393     )
 394     {
 395         FatalErrorIn
 396         (
 397             "meshRefinement::createBaffles"
 398             "(const labelList&, const labelList&)"
 399         )   << "Illegal size :"
 400             << " ownPatch:" << ownPatch.size()
 401             << " neiPatch:" << neiPatch.size()
 402             << ". Should be number of faces:" << mesh_.nFaces()
 403             << abort(FatalError);
 404     }
 405
 406     if (debug)
 407     {
 408         labelList syncedOwnPatch(ownPatch);
 409         syncTools::syncFaceList(mesh_, syncedOwnPatch, maxEqOp<label>(), false);
 410         labelList syncedNeiPatch(neiPatch);
 411         syncTools::syncFaceList(mesh_, syncedNeiPatch, maxEqOp<label>(), false);
 412
 413         forAll(syncedOwnPatch, faceI)
 414         {
 415             if
 416             (
 417                 (ownPatch[faceI] == -1 && syncedOwnPatch[faceI] != -1)
 418              || (neiPatch[faceI] == -1 && syncedNeiPatch[faceI] != -1)
 419             )
 420             {
 421                 FatalErrorIn
 422                 (
 423                     "meshRefinement::createBaffles"
 424                     "(const labelList&, const labelList&)"
 425                 )   << "Non synchronised at face:" << faceI
 426                     << " on patch:" << mesh_.boundaryMesh().whichPatch(faceI)
 427                     << " fc:" << mesh_.faceCentres()[faceI] << endl
 428                     << "ownPatch:" << ownPatch[faceI]
 429                     << " syncedOwnPatch:" << syncedOwnPatch[faceI]
 430                     << " neiPatch:" << neiPatch[faceI]
 431                     << " syncedNeiPatch:" << syncedNeiPatch[faceI]
 432                     << abort(FatalError);
 433             }
 434         }
 435     }
 436
 437     // Topochange container
 438     polyTopoChange meshMod(mesh_);
 439
 440     label nBaffles = 0;
 441
 442     forAll(ownPatch, faceI)
 443     {
 444         if (ownPatch[faceI] != -1)
 445         {
 446             // Create baffle or repatch face. Return label of inserted baffle
 447             // face.
 448             createBaffle
 449             (
 450                 faceI,
 451                 ownPatch[faceI],   // owner side patch
 452                 neiPatch[faceI],   // neighbour side patch
 453                 meshMod
 454             );
 455             nBaffles++;
 456         }
 457     }
 458
 459     // Change the mesh (no inflation, parallel sync)
 460     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
 461
 462     // Update fields
 463     mesh_.updateMesh(map);
 464
 465     // Move mesh if in inflation mode
 466     if (map().hasMotionPoints())
 467     {
 468         mesh_.movePoints(map().preMotionPoints());
 469     }
 470     else
 471     {
 472         // Delete mesh volumes.
 473         mesh_.clearOut();
 474     }
 475
 476     if (overwrite())
 477     {
 478         mesh_.setInstance(oldInstance());
 479     }
 480
 481     //- Redo the intersections on the newly create baffle faces. Note that
 482     //  this changes also the cell centre positions.
 483     faceSet baffledFacesSet(mesh_, "baffledFacesSet", 2*nBaffles);
 484
 485     const labelList& reverseFaceMap = map().reverseFaceMap();
 486     const labelList& faceMap = map().faceMap();
 487
 488     // Pick up owner side of baffle
 489     forAll(ownPatch, oldFaceI)
 490     {
 491         label faceI = reverseFaceMap[oldFaceI];
 492
 493         if (ownPatch[oldFaceI] != -1 && faceI >= 0)
 494         {
 495             const cell& ownFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
 496
 497             forAll(ownFaces, i)
 498             {
 499                 baffledFacesSet.insert(ownFaces[i]);
 500             }
 501         }
 502     }
 503     // Pick up neighbour side of baffle (added faces)
 504     forAll(faceMap, faceI)
 505     {
 506         label oldFaceI = faceMap[faceI];
 507
 508         if (oldFaceI >= 0 && reverseFaceMap[oldFaceI] != faceI)
 509         {
 510             const cell& ownFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
 511
 512             forAll(ownFaces, i)
 513             {
 514                 baffledFacesSet.insert(ownFaces[i]);
 515             }
 516         }
 517     }
 518     baffledFacesSet.sync(mesh_);
 519
 520     updateMesh(map, baffledFacesSet.toc());
 521
 522     return map;
 523 }
 524
 525
 526 // Return a list of coupled face pairs, i.e. faces that use the same vertices.
 527 // (this information is recalculated instead of maintained since would be too
 528 // hard across splitMeshRegions).
 529 Foam::List<Foam::labelPair> Foam::meshRefinement::getDuplicateFaces
 530 (
 531     const labelList& testFaces
 532 ) const
 533 {
 534     labelList duplicateFace
 535     (
 536         localPointRegion::findDuplicateFaces
 537         (
 538             mesh_,
 539             testFaces
 540         )
 541     );
 542
 543     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 544
 545     // Convert into list of coupled face pairs (mesh face labels).
 546     List<labelPair> duplicateFaces(testFaces.size());
 547     label dupI = 0;
 548
 549     forAll(duplicateFace, i)
 550     {
 551         label otherFaceI = duplicateFace[i];
 552
 553         if (otherFaceI != -1 && i < otherFaceI)
 554         {
 555             label meshFace0 = testFaces[i];
 556             label patch0 = patches.whichPatch(meshFace0);
 557             label meshFace1 = testFaces[otherFaceI];
 558             label patch1 = patches.whichPatch(meshFace1);
 559
 560             if
 561             (
 562                 (patch0 != -1 && isA<processorPolyPatch>(patches[patch0]))
 563              || (patch1 != -1 && isA<processorPolyPatch>(patches[patch1]))
 564             )
 565             {
 566                 FatalErrorIn
 567                 (
 568                     "meshRefinement::getDuplicateFaces"
 569                     "(const bool, const labelList&)"
 570                 )   << "One of two duplicate faces is on"
 571                     << " processorPolyPatch."
 572                     << "This is not allowed." << nl
 573                     << "Face:" << meshFace0
 574                     << " is on patch:" << patches[patch0].name()
 575                     << nl
 576                     << "Face:" << meshFace1
 577                     << " is on patch:" << patches[patch1].name()
 578                     << abort(FatalError);
 579             }
 580
 581             duplicateFaces[dupI++] = labelPair(meshFace0, meshFace1);
 582         }
 583     }
 584     duplicateFaces.setSize(dupI);
 585
 586     Info<< "getDuplicateFaces : found " << returnReduce(dupI, sumOp<label>())
 587         << " pairs of duplicate faces." << nl << endl;
 588
 589
 590     if (debug)
 591     {
 592         faceSet duplicateFaceSet(mesh_, "duplicateFaces", 2*dupI);
 593
 594         forAll(duplicateFaces, i)
 595         {
 596             duplicateFaceSet.insert(duplicateFaces[i][0]);
 597             duplicateFaceSet.insert(duplicateFaces[i][1]);
 598         }
 599         Pout<< "Writing duplicate faces (baffles) to faceSet "
 600             << duplicateFaceSet.name() << nl << endl;
 601         duplicateFaceSet.write();
 602     }
 603
 604     return duplicateFaces;
 605 }
 606
 607
 608 // Extract those baffles (duplicate) faces that are on the edge of a baffle
 609 // region. These are candidates for merging.
 610 // Done by counting the number of baffles faces per mesh edge. If edge
 611 // has 2 boundary faces and both are baffle faces it is the edge of a baffle
 612 // region.
 613 Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
 614 (
 615     const List<labelPair>& couples
 616 ) const
 617 {
 618     // Construct addressing engine for all duplicate faces (only one
 619     // for each pair)
 620
 621     // Duplicate faces in mesh labels (first face of each pair only)
 622     // (reused later on to mark off filtered couples. see below)
 623     labelList duplicateFaces(couples.size());
 624
 625
 626     // All duplicate faces on edge of the patch are to be merged.
 627     // So we count for all edges of duplicate faces how many duplicate
 628     // faces use them.
 629     labelList nBafflesPerEdge(mesh_.nEdges(), 0);
 630
 631
 632
 633     // Count number of boundary faces per edge
 634     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 635
 636     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 637
 638     forAll(patches, patchI)
 639     {
 640         const polyPatch& pp = patches[patchI];
 641
 642         // Count number of boundary faces. Discard coupled boundary faces.
 643         if (!pp.coupled())
 644         {
 645             label faceI = pp.start();
 646
 647             forAll(pp, i)
 648             {
 649                 const labelList& fEdges = mesh_.faceEdges(faceI);
 650
 651                 forAll(fEdges, fEdgeI)
 652                 {
 653                     nBafflesPerEdge[fEdges[fEdgeI]]++;
 654                 }
 655                 faceI++;
 656             }
 657         }
 658     }
 659
 660
 661     DynamicList<label> fe0;
 662     DynamicList<label> fe1;
 663
 664
 665     // Count number of duplicate boundary faces per edge
 666     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 667
 668     forAll(couples, i)
 669     {
 670         const labelList& fEdges0 = mesh_.faceEdges(couples[i].first(), fe0);
 671
 672         forAll(fEdges0, fEdgeI)
 673         {
 674             nBafflesPerEdge[fEdges0[fEdgeI]] += 1000000;
 675         }
 676
 677         const labelList& fEdges1 = mesh_.faceEdges(couples[i].second(), fe1);
 678
 679         forAll(fEdges1, fEdgeI)
 680         {
 681             nBafflesPerEdge[fEdges1[fEdgeI]] += 1000000;
 682         }
 683     }
 684
 685     // Add nBaffles on shared edges
 686     syncTools::syncEdgeList
 687     (
 688         mesh_,
 689         nBafflesPerEdge,
 690         plusEqOp<label>(),  // in-place add
 691         0,                  // initial value
 692         false               // no separation
 693     );
 694
 695
 696     // Baffles which are not next to other boundaries and baffles will have
 697     // value 2*1000000+2*1
 698
 699     List<labelPair> filteredCouples(couples.size());
 700     label filterI = 0;
 701
 702     forAll(couples, i)
 703     {
 704         const labelPair& couple = couples[i];
 705
 706         if
 707         (
 708             patches.whichPatch(couple.first())
 709          == patches.whichPatch(couple.second())
 710         )
 711         {
 712             const labelList& fEdges = mesh_.faceEdges(couples[i].first());
 713
 714             forAll(fEdges, fEdgeI)
 715             {
 716                 label edgeI = fEdges[fEdgeI];
 717
 718                 if (nBafflesPerEdge[edgeI] == 2*1000000+2*1)
 719                 {
 720                     filteredCouples[filterI++] = couples[i];
 721                     break;
 722                 }
 723             }
 724         }
 725     }
 726     filteredCouples.setSize(filterI);
 727
 728     //Info<< "filterDuplicateFaces : from "
 729     //    << returnReduce(couples.size(), sumOp<label>())
 730     //    << " down to "
 731     //    << returnReduce(filteredCouples.size(), sumOp<label>())
 732     //    << " baffles." << nl << endl;
 733
 734     return filteredCouples;
 735 }
 736
 737
 738 // Merge baffles. Gets pairs of faces.
 739 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeBaffles
 740 (
 741     const List<labelPair>& couples
 742 )
 743 {
 744     // Mesh change engine
 745     polyTopoChange meshMod(mesh_);
 746
 747     const faceList& faces = mesh_.faces();
 748     const labelList& faceOwner = mesh_.faceOwner();
 749     const faceZoneMesh& faceZones = mesh_.faceZones();
 750
 751     forAll(couples, i)
 752     {
 753         label face0 = couples[i].first();
 754         label face1 = couples[i].second();
 755
 756         // face1 < 0 signals a coupled face that has been converted to baffle.
 757
 758         label own0 = faceOwner[face0];
 759         label own1 = faceOwner[face1];
 760
 761         if (face1 < 0 || own0 < own1)
 762         {
 763             // Use face0 as the new internal face.
 764             label zoneID = faceZones.whichZone(face0);
 765             bool zoneFlip = false;
 766
 767             if (zoneID >= 0)
 768             {
 769                 const faceZone& fZone = faceZones[zoneID];
 770                 zoneFlip = fZone.flipMap()[fZone.whichFace(face0)];
 771             }
 772
 773             label nei = (face1 < 0 ? -1 : own1);
 774
 775             meshMod.setAction(polyRemoveFace(face1));
 776             meshMod.setAction
 777             (
 778                 polyModifyFace
 779                 (
 780                     faces[face0],           // modified face
 781                     face0,                  // label of face being modified
 782                     own0,                   // owner
 783                     nei,                    // neighbour
 784                     false,                  // face flip
 785                     -1,                     // patch for face
 786                     false,                  // remove from zone
 787                     zoneID,                 // zone for face
 788                     zoneFlip                // face flip in zone
 789                 )
 790             );
 791         }
 792         else
 793         {
 794             // Use face1 as the new internal face.
 795             label zoneID = faceZones.whichZone(face1);
 796             bool zoneFlip = false;
 797
 798             if (zoneID >= 0)
 799             {
 800                 const faceZone& fZone = faceZones[zoneID];
 801                 zoneFlip = fZone.flipMap()[fZone.whichFace(face1)];
 802             }
 803
 804             meshMod.setAction(polyRemoveFace(face0));
 805             meshMod.setAction
 806             (
 807                 polyModifyFace
 808                 (
 809                     faces[face1],           // modified face
 810                     face1,                  // label of face being modified
 811                     own1,                   // owner
 812                     own0,                   // neighbour
 813                     false,                  // face flip
 814                     -1,                     // patch for face
 815                     false,                  // remove from zone
 816                     zoneID,                 // zone for face
 817                     zoneFlip                // face flip in zone
 818                 )
 819             );
 820         }
 821     }
 822
 823     // Change the mesh (no inflation)
 824     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
 825
 826     // Update fields
 827     mesh_.updateMesh(map);
 828
 829     // Move mesh (since morphing does not do this)
 830     if (map().hasMotionPoints())
 831     {
 832         mesh_.movePoints(map().preMotionPoints());
 833     }
 834     else
 835     {
 836         // Delete mesh volumes.
 837         mesh_.clearOut();
 838     }
 839
 840     if (overwrite())
 841     {
 842         mesh_.setInstance(oldInstance());
 843     }
 844
 845     // Update intersections. Recalculate intersections on merged faces since
 846     // this seems to give problems? Note: should not be nessecary since
 847     // baffles preserve intersections from when they were created.
 848     labelList newExposedFaces(2*couples.size());
 849     label newI = 0;
 850
 851     forAll(couples, i)
 852     {
 853         label newFace0 = map().reverseFaceMap()[couples[i].first()];
 854         if (newFace0 != -1)
 855         {
 856             newExposedFaces[newI++] = newFace0;
 857         }
 858
 859         label newFace1 = map().reverseFaceMap()[couples[i].second()];
 860         if (newFace1 != -1)
 861         {
 862             newExposedFaces[newI++] = newFace1;
 863         }
 864     }
 865     newExposedFaces.setSize(newI);
 866     updateMesh(map, newExposedFaces);
 867
 868     return map;
 869 }
 870
 871
 872 // Finds region per cell for cells inside closed named surfaces
 873 void Foam::meshRefinement::findCellZoneGeometric
 874 (
 875     const labelList& closedNamedSurfaces,   // indices of closed surfaces
 876     labelList& namedSurfaceIndex,           // per face index of named surface
 877     const labelList& surfaceToCellZone,     // cell zone index per surface
 878
 879     labelList& cellToZone
 880 ) const
 881 {
 882     const pointField& cellCentres = mesh_.cellCentres();
 883     const labelList& faceOwner = mesh_.faceOwner();
 884     const labelList& faceNeighbour = mesh_.faceNeighbour();
 885
 886     // Check if cell centre is inside
 887     labelList insideSurfaces;
 888     surfaces_.findInside
 889     (
 890         closedNamedSurfaces,
 891         cellCentres,
 892         insideSurfaces
 893     );
 894
 895     forAll(insideSurfaces, cellI)
 896     {
 897         if (cellToZone[cellI] == -2)
 898         {
 899             label surfI = insideSurfaces[cellI];
 900
 901             if (surfI != -1)
 902             {
 903                 cellToZone[cellI] = surfaceToCellZone[surfI];
 904             }
 905         }
 906     }
 907
 908
 909     // Some cells with cell centres close to surface might have
 910     // had been put into wrong surface. Recheck with perturbed cell centre.
 911
 912
 913     // 1. Collect points
 914
 915     // Count points to test.
 916     label nCandidates = 0;
 917     forAll(namedSurfaceIndex, faceI)
 918     {
 919         label surfI = namedSurfaceIndex[faceI];
 920
 921         if (surfI != -1)
 922         {
 923             if (mesh_.isInternalFace(faceI))
 924             {
 925                 nCandidates += 2;
 926             }
 927             else
 928             {
 929                 nCandidates += 1;
 930             }
 931         }
 932     }
 933
 934     // Collect points.
 935     pointField candidatePoints(nCandidates);
 936     nCandidates = 0;
 937     forAll(namedSurfaceIndex, faceI)
 938     {
 939         label surfI = namedSurfaceIndex[faceI];
 940
 941         if (surfI != -1)
 942         {
 943             label own = faceOwner[faceI];
 944             const point& ownCc = cellCentres[own];
 945
 946             if (mesh_.isInternalFace(faceI))
 947             {
 948                 label nei = faceNeighbour[faceI];
 949                 const point& neiCc = cellCentres[nei];
 950                 // Perturbed cc
 951                 const vector d = 1E-4*(neiCc - ownCc);
 952                 candidatePoints[nCandidates++] = ownCc-d;
 953                 candidatePoints[nCandidates++] = neiCc+d;
 954             }
 955             else
 956             {
 957                 const point& neiFc = mesh_.faceCentres()[faceI];
 958                 // Perturbed cc
 959                 const vector d = 1E-4*(neiFc - ownCc);
 960                 candidatePoints[nCandidates++] = ownCc-d;
 961             }
 962         }
 963     }
 964
 965
 966     // 2. Test points for inside
 967
 968     surfaces_.findInside
 969     (
 970         closedNamedSurfaces,
 971         candidatePoints,
 972         insideSurfaces
 973     );
 974
 975
 976     // 3. Update zone information
 977
 978     nCandidates = 0;
 979     forAll(namedSurfaceIndex, faceI)
 980     {
 981         label surfI = namedSurfaceIndex[faceI];
 982
 983         if (surfI != -1)
 984         {
 985             label own = faceOwner[faceI];
 986
 987             if (mesh_.isInternalFace(faceI))
 988             {
 989                 label ownSurfI = insideSurfaces[nCandidates++];
 990                 if (ownSurfI != -1)
 991                 {
 992                     cellToZone[own] = surfaceToCellZone[ownSurfI];
 993                 }
 994
 995                 label neiSurfI = insideSurfaces[nCandidates++];
 996                 if (neiSurfI != -1)
 997                 {
 998                     label nei = faceNeighbour[faceI];
 999
1000                     cellToZone[nei] = surfaceToCellZone[neiSurfI];
1001                 }
1002             }
1003             else
1004             {
1005                 label ownSurfI = insideSurfaces[nCandidates++];
1006                 if (ownSurfI != -1)
1007                 {
1008                     cellToZone[own] = surfaceToCellZone[ownSurfI];
1009                 }
1010             }
1011         }
1012     }
1013
1014
1015     // Adapt the namedSurfaceIndex
1016     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
1017     // for if any cells were not completely covered.
1018
1019     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1020     {
1021         label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
1022         label neiZone = cellToZone[mesh_.faceNeighbour()[faceI]];
1023
1024         if (namedSurfaceIndex[faceI] == -1 && (ownZone != neiZone))
1025         {
1026             // Give face the zone of max cell zone
1027             namedSurfaceIndex[faceI] = findIndex
1028             (
1029                 surfaceToCellZone,
1030                 max(ownZone, neiZone)
1031             );
1032         }
1033     }
1034
1035     labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces());
1036     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1037
1038     forAll(patches, patchI)
1039     {
1040         const polyPatch& pp = patches[patchI];
1041
1042         if (pp.coupled())
1043         {
1044             forAll(pp, i)
1045             {
1046                 label faceI = pp.start()+i;
1047                 label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
1048                 neiCellZone[faceI-mesh_.nInternalFaces()] = ownZone;
1049             }
1050         }
1051     }
1052     syncTools::swapBoundaryFaceList(mesh_, neiCellZone, false);
1053
1054     forAll(patches, patchI)
1055     {
1056         const polyPatch& pp = patches[patchI];
1057
1058         if (pp.coupled())
1059         {
1060             forAll(pp, i)
1061             {
1062                 label faceI = pp.start()+i;
1063                 label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
1064                 label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
1065
1066                 if (namedSurfaceIndex[faceI] == -1 && (ownZone != neiZone))
1067                 {
1068                     // Give face the max cell zone
1069                     namedSurfaceIndex[faceI] = findIndex
1070                     (
1071                         surfaceToCellZone,
1072                         max(ownZone, neiZone)
1073                     );
1074                 }
1075             }
1076         }
1077     }
1078
1079     // Sync
1080     syncTools::syncFaceList
1081     (
1082         mesh_,
1083         namedSurfaceIndex,
1084         maxEqOp<label>(),
1085         false
1086     );
1087 }
1088
1089
1090 bool Foam::meshRefinement::calcRegionToZone
1091 (
1092     const label surfZoneI,
1093     const label ownRegion,
1094     const label neiRegion,
1095
1096     labelList& regionToCellZone
1097 ) const
1098 {
1099     bool changed = false;
1100
1101     // Check whether inbetween different regions
1102     if (ownRegion != neiRegion)
1103     {
1104         // Jump. Change one of the sides to my type.
1105
1106         // 1. Interface between my type and unset region.
1107         // Set region to keepRegion
1108
1109         if (regionToCellZone[ownRegion] == -2)
1110         {
1111             if (regionToCellZone[neiRegion] == surfZoneI)
1112             {
1113                 // Face between unset and my region. Put unset
1114                 // region into keepRegion
1115                 regionToCellZone[ownRegion] = -1;
1116                 changed = true;
1117             }
1118             else if (regionToCellZone[neiRegion] != -2)
1119             {
1120                 // Face between unset and other region.
1121                 // Put unset region into my region
1122                 regionToCellZone[ownRegion] = surfZoneI;
1123                 changed = true;
1124             }
1125         }
1126         else if (regionToCellZone[neiRegion] == -2)
1127         {
1128             if (regionToCellZone[ownRegion] == surfZoneI)
1129             {
1130                 // Face between unset and my region. Put unset
1131                 // region into keepRegion
1132                 regionToCellZone[neiRegion] = -1;
1133                 changed = true;
1134             }
1135             else if (regionToCellZone[ownRegion] != -2)
1136             {
1137                 // Face between unset and other region.
1138                 // Put unset region into my region
1139                 regionToCellZone[neiRegion] = surfZoneI;
1140                 changed = true;
1141             }
1142         }
1143     }
1144     return changed;
1145 }
1146
1147
1148 // Finds region per cell. Assumes:
1149 // - region containing keepPoint does not go into a cellZone
1150 // - all other regions can be found by crossing faces marked in
1151 //   namedSurfaceIndex.
1152 void Foam::meshRefinement::findCellZoneTopo
1153 (
1154     const point& keepPoint,
1155     const labelList& namedSurfaceIndex,
1156     const labelList& surfaceToCellZone,
1157     labelList& cellToZone
1158 ) const
1159 {
1160     // Analyse regions. Reuse regionsplit
1161     boolList blockedFace(mesh_.nFaces());
1162
1163     forAll(namedSurfaceIndex, faceI)
1164     {
1165         if (namedSurfaceIndex[faceI] == -1)
1166         {
1167             blockedFace[faceI] = false;
1168         }
1169         else
1170         {
1171             blockedFace[faceI] = true;
1172         }
1173     }
1174     // No need to sync since namedSurfaceIndex already is synced
1175
1176     // Set region per cell based on walking
1177     regionSplit cellRegion(mesh_, blockedFace);
1178     blockedFace.clear();
1179
1180     // Per mesh region the zone the cell should be put in.
1181     // -2   : not analysed yet
1182     // -1   : keepPoint region. Not put into any cellzone.
1183     // >= 0 : index of cellZone
1184     labelList regionToCellZone(cellRegion.nRegions(), -2);
1185
1186     // See which cells already are set in the cellToZone (from geometric
1187     // searching) and use these to take over their zones.
1188     // Note: could be improved to count number of cells per region.
1189     forAll(cellToZone, cellI)
1190     {
1191         if (cellToZone[cellI] != -2)
1192         {
1193             regionToCellZone[cellRegion[cellI]] = cellToZone[cellI];
1194         }
1195     }
1196
1197
1198
1199     // Find the region containing the keepPoint
1200     label keepRegionI = -1;
1201
1202     label cellI = mesh_.findCell(keepPoint);
1203
1204     if (cellI != -1)
1205     {
1206         keepRegionI = cellRegion[cellI];
1207     }
1208     reduce(keepRegionI, maxOp<label>());
1209
1210     Info<< "Found point " << keepPoint << " in cell " << cellI
1211         << " in global region " << keepRegionI
1212         << " out of " << cellRegion.nRegions() << " regions." << endl;
1213
1214     if (keepRegionI == -1)
1215     {
1216         FatalErrorIn
1217         (
1218             "meshRefinement::findCellZoneTopo"
1219             "(const point&, const labelList&, const labelList&, labelList&)"
1220         )   << "Point " << keepPoint
1221             << " is not inside the mesh." << nl
1222             << "Bounding box of the mesh:" << mesh_.globalData().bb()
1223             << exit(FatalError);
1224     }
1225
1226     // Mark default region with zone -1.
1227     if (regionToCellZone[keepRegionI] == -2)
1228     {
1229         regionToCellZone[keepRegionI] = -1;
1230     }
1231
1232
1233     // Find correspondence between cell zone and surface
1234     // by changing cell zone every time we cross a surface.
1235     while (true)
1236     {
1237         // Synchronise regionToCellZone.
1238         // Note:
1239         // - region numbers are identical on all processors
1240         // - keepRegion is identical ,,
1241         // - cellZones are identical ,,
1242         // This done at top of loop to account for geometric matching
1243         // not being synchronised.
1244         Pstream::listCombineGather(regionToCellZone, maxEqOp<label>());
1245         Pstream::listCombineScatter(regionToCellZone);
1246
1247
1248         bool changed = false;
1249
1250         // Internal faces
1251
1252         for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1253         {
1254             label surfI = namedSurfaceIndex[faceI];
1255
1256             if (surfI != -1)
1257             {
1258                 // Calculate region to zone from cellRegions on either side
1259                 // of internal face.
1260                 bool changedCell = calcRegionToZone
1261                 (
1262                     surfaceToCellZone[surfI],
1263                     cellRegion[mesh_.faceOwner()[faceI]],
1264                     cellRegion[mesh_.faceNeighbour()[faceI]],
1265                     regionToCellZone
1266                 );
1267
1268                 changed = changed | changedCell;
1269             }
1270         }
1271
1272         // Boundary faces
1273
1274         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1275
1276         // Get coupled neighbour cellRegion
1277         labelList neiCellRegion(mesh_.nFaces()-mesh_.nInternalFaces());
1278         forAll(patches, patchI)
1279         {
1280             const polyPatch& pp = patches[patchI];
1281
1282             if (pp.coupled())
1283             {
1284                 forAll(pp, i)
1285                 {
1286                     label faceI = pp.start()+i;
1287                     neiCellRegion[faceI-mesh_.nInternalFaces()] =
1288                         cellRegion[mesh_.faceOwner()[faceI]];
1289                 }
1290             }
1291         }
1292         syncTools::swapBoundaryFaceList(mesh_, neiCellRegion, false);
1293
1294         // Calculate region to zone from cellRegions on either side of coupled
1295         // face.
1296         forAll(patches, patchI)
1297         {
1298             const polyPatch& pp = patches[patchI];
1299
1300             if (pp.coupled())
1301             {
1302                 forAll(pp, i)
1303                 {
1304                     label faceI = pp.start()+i;
1305
1306                     label surfI = namedSurfaceIndex[faceI];
1307
1308                     if (surfI != -1)
1309                     {
1310                         bool changedCell = calcRegionToZone
1311                         (
1312                             surfaceToCellZone[surfI],
1313                             cellRegion[mesh_.faceOwner()[faceI]],
1314                             neiCellRegion[faceI-mesh_.nInternalFaces()],
1315                             regionToCellZone
1316                         );
1317
1318                         changed = changed | changedCell;
1319                     }
1320                 }
1321             }
1322         }
1323
1324         if (!returnReduce(changed, orOp<bool>()))
1325         {
1326             break;
1327         }
1328     }
1329
1330
1331     forAll(regionToCellZone, regionI)
1332     {
1333         label zoneI = regionToCellZone[regionI];
1334
1335         if (zoneI ==  -2)
1336         {
1337             FatalErrorIn
1338             (
1339                 "meshRefinement::findCellZoneTopo"
1340                 "(const point&, const labelList&, const labelList&, labelList&)"
1341             )   << "For region " << regionI << " haven't set cell zone."
1342                 << exit(FatalError);
1343         }
1344     }
1345
1346     if (debug)
1347     {
1348         forAll(regionToCellZone, regionI)
1349         {
1350             Pout<< "Region " << regionI
1351                 << " becomes cellZone:" << regionToCellZone[regionI]
1352                 << endl;
1353         }
1354     }
1355
1356     // Rework into cellToZone
1357     forAll(cellToZone, cellI)
1358     {
1359         cellToZone[cellI] = regionToCellZone[cellRegion[cellI]];
1360     }
1361 }
1362
1363
1364 // Make namedSurfaceIndex consistent with cellToZone - clear out any blocked
1365 // faces inbetween same cell zone.
1366 void Foam::meshRefinement::makeConsistentFaceIndex
1367 (
1368     const labelList& cellToZone,
1369     labelList& namedSurfaceIndex
1370 ) const
1371 {
1372     const labelList& faceOwner = mesh_.faceOwner();
1373     const labelList& faceNeighbour = mesh_.faceNeighbour();
1374
1375     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1376     {
1377         label ownZone = cellToZone[faceOwner[faceI]];
1378         label neiZone = cellToZone[faceNeighbour[faceI]];
1379
1380         if (ownZone == neiZone && namedSurfaceIndex[faceI] != -1)
1381         {
1382             namedSurfaceIndex[faceI] = -1;
1383         }
1384         else if (ownZone != neiZone && namedSurfaceIndex[faceI] == -1)
1385         {
1386             FatalErrorIn("meshRefinement::zonify()")
1387                 << "Different cell zones on either side of face " << faceI
1388                 << " at " << mesh_.faceCentres()[faceI]
1389                 << " but face not marked with a surface."
1390                 << abort(FatalError);
1391         }
1392     }
1393
1394     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1395
1396     // Get coupled neighbour cellZone
1397     labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces());
1398     forAll(patches, patchI)
1399     {
1400         const polyPatch& pp = patches[patchI];
1401
1402         if (pp.coupled())
1403         {
1404             forAll(pp, i)
1405             {
1406                 label faceI = pp.start()+i;
1407                 neiCellZone[faceI-mesh_.nInternalFaces()] =
1408                     cellToZone[mesh_.faceOwner()[faceI]];
1409             }
1410         }
1411     }
1412     syncTools::swapBoundaryFaceList(mesh_, neiCellZone, false);
1413
1414     // Use coupled cellZone to do check
1415     forAll(patches, patchI)
1416     {
1417         const polyPatch& pp = patches[patchI];
1418
1419         if (pp.coupled())
1420         {
1421             forAll(pp, i)
1422             {
1423                 label faceI = pp.start()+i;
1424
1425                 label ownZone = cellToZone[faceOwner[faceI]];
1426                 label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
1427
1428                 if (ownZone == neiZone && namedSurfaceIndex[faceI] != -1)
1429                 {
1430                     namedSurfaceIndex[faceI] = -1;
1431                 }
1432                 else if (ownZone != neiZone && namedSurfaceIndex[faceI] == -1)
1433                 {
1434                     FatalErrorIn("meshRefinement::zonify()")
1435                         << "Different cell zones on either side of face "
1436                         << faceI << " at " << mesh_.faceCentres()[faceI]
1437                         << " but face not marked with a surface."
1438                         << abort(FatalError);
1439                 }
1440             }
1441         }
1442     }
1443 }
1444
1445
1446 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
1447
1448 // Split off unreachable areas of mesh.
1449 void Foam::meshRefinement::baffleAndSplitMesh
1450 (
1451     const bool handleSnapProblems,
1452     const bool removeEdgeConnectedCells,
1453     const scalarField& perpendicularAngle,
1454     const bool mergeFreeStanding,
1455     const dictionary& motionDict,
1456     Time& runTime,
1457     const labelList& globalToPatch,
1458     const point& keepPoint
1459 )
1460 {
1461     // Introduce baffles
1462     // ~~~~~~~~~~~~~~~~~
1463
1464     // Split the mesh along internal faces wherever there is a pierce between
1465     // two cell centres.
1466
1467     Info<< "Introducing baffles for "
1468         << returnReduce(countHits(), sumOp<label>())
1469         << " faces that are intersected by the surface." << nl << endl;
1470
1471     // Swap neighbouring cell centres and cell level
1472     labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
1473     pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
1474     calcNeighbourData(neiLevel, neiCc);
1475
1476     labelList ownPatch, neiPatch;
1477     getBafflePatches
1478     (
1479         globalToPatch,
1480         neiLevel,
1481         neiCc,
1482
1483         ownPatch,
1484         neiPatch
1485     );
1486
1487     if (debug)
1488     {
1489         runTime++;
1490     }
1491
1492     createBaffles(ownPatch, neiPatch);
1493
1494     if (debug)
1495     {
1496         // Debug:test all is still synced across proc patches
1497         checkData();
1498     }
1499
1500     Info<< "Created baffles in = "
1501         << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
1502
1503     printMeshInfo(debug, "After introducing baffles");
1504
1505     if (debug)
1506     {
1507         Pout<< "Writing baffled mesh to time " << timeName()
1508             << endl;
1509         write(debug, runTime.path()/"baffles");
1510         Pout<< "Dumped debug data in = "
1511             << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
1512     }
1513
1514
1515     // Introduce baffles to delete problem cells
1516     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1517
1518     // Create some additional baffles where we want surface cells removed.
1519
1520     if (handleSnapProblems)
1521     {
1522         Info<< nl
1523             << "Introducing baffles to block off problem cells" << nl
1524             << "----------------------------------------------" << nl
1525             << endl;
1526
1527         labelList facePatch
1528         (
1529             markFacesOnProblemCells
1530             (
1531                 motionDict,
1532                 removeEdgeConnectedCells,
1533                 perpendicularAngle,
1534                 globalToPatch
1535             )
1536             //markFacesOnProblemCellsGeometric(motionDict)
1537         );
1538         Info<< "Analyzed problem cells in = "
1539             << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
1540
1541         if (debug)
1542         {
1543             faceSet problemTopo(mesh_, "problemFacesTopo", 100);
1544
1545             const labelList facePatchTopo
1546             (
1547                 markFacesOnProblemCells
1548                 (
1549                     motionDict,
1550                     removeEdgeConnectedCells,
1551                     perpendicularAngle,
1552                     globalToPatch
1553                 )
1554             );
1555             forAll(facePatchTopo, faceI)
1556             {
1557                 if (facePatchTopo[faceI] != -1)
1558                 {
1559                     problemTopo.insert(faceI);
1560                 }
1561             }
1562             Pout<< "Dumping " << problemTopo.size()
1563                 << " problem faces to " << problemTopo.objectPath() << endl;
1564             problemTopo.write();
1565         }
1566
1567         Info<< "Introducing baffles to delete problem cells." << nl << endl;
1568
1569         if (debug)
1570         {
1571             runTime++;
1572         }
1573
1574         // Create baffles with same owner and neighbour for now.
1575         createBaffles(facePatch, facePatch);
1576
1577         if (debug)
1578         {
1579             // Debug:test all is still synced across proc patches
1580             checkData();
1581         }
1582         Info<< "Created baffles in = "
1583             << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
1584
1585         printMeshInfo(debug, "After introducing baffles");
1586
1587         if (debug)
1588         {
1589             Pout<< "Writing extra baffled mesh to time "
1590                 << timeName() << endl;
1591             write(debug, runTime.path()/"extraBaffles");
1592             Pout<< "Dumped debug data in = "
1593                 << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
1594         }
1595     }
1596
1597
1598     // Select part of mesh
1599     // ~~~~~~~~~~~~~~~~~~~
1600
1601     Info<< nl
1602         << "Remove unreachable sections of mesh" << nl
1603         << "-----------------------------------" << nl
1604         << endl;
1605
1606     if (debug)
1607     {
1608         runTime++;
1609     }
1610
1611     splitMeshRegions(keepPoint);
1612
1613     if (debug)
1614     {
1615         // Debug:test all is still synced across proc patches
1616         checkData();
1617     }
1618     Info<< "Split mesh in = "
1619         << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
1620
1621     printMeshInfo(debug, "After subsetting");
1622
1623     if (debug)
1624     {
1625         Pout<< "Writing subsetted mesh to time " << timeName()
1626             << endl;
1627         write(debug, runTime.path()/timeName());
1628         Pout<< "Dumped debug data in = "
1629             << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
1630     }
1631
1632
1633     // Merge baffles
1634     // ~~~~~~~~~~~~~
1635
1636     if (mergeFreeStanding)
1637     {
1638         Info<< nl
1639             << "Merge free-standing baffles" << nl
1640             << "---------------------------" << nl
1641             << endl;
1642
1643         if (debug)
1644         {
1645             runTime++;
1646         }
1647
1648         // List of pairs of freestanding baffle faces.
1649         List<labelPair> couples
1650         (
1651             filterDuplicateFaces    // filter out freestanding baffles
1652             (
1653                 getDuplicateFaces   // get all baffles
1654                 (
1655                     identity(mesh_.nFaces()-mesh_.nInternalFaces())
1656                    +mesh_.nInternalFaces()
1657                 )
1658             )
1659         );
1660
1661         label nCouples = couples.size();
1662         reduce(nCouples, sumOp<label>());
1663
1664         Info<< "Detected free-standing baffles : " << nCouples << endl;
1665
1666         if (nCouples > 0)
1667         {
1668             // Actually merge baffles. Note: not exactly parallellized. Should
1669             // convert baffle faces into processor faces if they resulted
1670             // from them.
1671             mergeBaffles(couples);
1672
1673             if (debug)
1674             {
1675                 // Debug:test all is still synced across proc patches
1676                 checkData();
1677             }
1678         }
1679         Info<< "Merged free-standing baffles in = "
1680             << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
1681     }
1682 }
1683
1684
1685 // Split off (with optional buffer layers) unreachable areas of mesh.
1686 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
1687 (
1688     const label nBufferLayers,
1689     const labelList& globalToPatch,
1690     const point& keepPoint
1691 )
1692 {
1693     // Determine patches to put intersections into
1694     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1695
1696     // Swap neighbouring cell centres and cell level
1697     labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
1698     pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
1699     calcNeighbourData(neiLevel, neiCc);
1700
1701     labelList ownPatch, neiPatch;
1702     getBafflePatches
1703     (
1704         globalToPatch,
1705         neiLevel,
1706         neiCc,
1707
1708         ownPatch,
1709         neiPatch
1710     );
1711
1712     // Analyse regions. Reuse regionsplit
1713     boolList blockedFace(mesh_.nFaces(), false);
1714
1715     forAll(ownPatch, faceI)
1716     {
1717         if (ownPatch[faceI] != -1 || neiPatch[faceI] != -1)
1718         {
1719             blockedFace[faceI] = true;
1720         }
1721     }
1722     syncTools::syncFaceList(mesh_, blockedFace, orEqOp<bool>(), false);
1723
1724     // Set region per cell based on walking
1725     regionSplit cellRegion(mesh_, blockedFace);
1726     blockedFace.clear();
1727
1728     // Find the region containing the keepPoint
1729     label keepRegionI = -1;
1730
1731     label cellI = mesh_.findCell(keepPoint);
1732
1733     if (cellI != -1)
1734     {
1735         keepRegionI = cellRegion[cellI];
1736     }
1737     reduce(keepRegionI, maxOp<label>());
1738
1739     Info<< "Found point " << keepPoint << " in cell " << cellI
1740         << " in global region " << keepRegionI
1741         << " out of " << cellRegion.nRegions() << " regions." << endl;
1742
1743     if (keepRegionI == -1)
1744     {
1745         FatalErrorIn
1746         (
1747             "meshRefinement::splitMesh"
1748             "(const label, const labelList&, const point&)"
1749         )   << "Point " << keepPoint
1750             << " is not inside the mesh." << nl
1751             << "Bounding box of the mesh:" << mesh_.globalData().bb()
1752             << exit(FatalError);
1753     }
1754
1755
1756     // Walk out nBufferlayers from region split
1757     // (modifies cellRegion, ownPatch)
1758     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1759     // Takes over face patch onto points and then back to faces and cells
1760     // (so cell-face-point walk)
1761
1762     const labelList& faceOwner = mesh_.faceOwner();
1763     const labelList& faceNeighbour = mesh_.faceNeighbour();
1764
1765     // Patch for exposed faces for lack of anything sensible.
1766     label defaultPatch = 0;
1767     if (globalToPatch.size())
1768     {
1769         defaultPatch = globalToPatch[0];
1770     }
1771
1772     for (label i = 0; i < nBufferLayers; i++)
1773     {
1774         // 1. From cells (via faces) to points
1775
1776         labelList pointBaffle(mesh_.nPoints(), -1);
1777
1778         forAll(faceNeighbour, faceI)
1779         {
1780             const face& f = mesh_.faces()[faceI];
1781
1782             label ownRegion = cellRegion[faceOwner[faceI]];
1783             label neiRegion = cellRegion[faceNeighbour[faceI]];
1784
1785             if (ownRegion == keepRegionI && neiRegion != keepRegionI)
1786             {
1787                 // Note max(..) since possibly regionSplit might have split
1788                 // off extra unreachable parts of mesh. Note: or can this only
1789                 // happen for boundary faces?
1790                 forAll(f, fp)
1791                 {
1792                     pointBaffle[f[fp]] = max(defaultPatch, ownPatch[faceI]);
1793                 }
1794             }
1795             else if (ownRegion != keepRegionI && neiRegion == keepRegionI)
1796             {
1797                 label newPatchI = neiPatch[faceI];
1798                 if (newPatchI == -1)
1799                 {
1800                     newPatchI = max(defaultPatch, ownPatch[faceI]);
1801                 }
1802                 forAll(f, fp)
1803                 {
1804                     pointBaffle[f[fp]] = newPatchI;
1805                 }
1806             }
1807         }
1808         for
1809         (
1810             label faceI = mesh_.nInternalFaces();
1811             faceI < mesh_.nFaces();
1812             faceI++
1813         )
1814         {
1815             const face& f = mesh_.faces()[faceI];
1816
1817             label ownRegion = cellRegion[faceOwner[faceI]];
1818
1819             if (ownRegion == keepRegionI)
1820             {
1821                 forAll(f, fp)
1822                 {
1823                     pointBaffle[f[fp]] = max(defaultPatch, ownPatch[faceI]);
1824                 }
1825             }
1826         }
1827
1828         // Sync
1829         syncTools::syncPointList
1830         (
1831             mesh_,
1832             pointBaffle,
1833             maxEqOp<label>(),
1834             -1,                 // null value
1835             false               // no separation
1836         );
1837
1838
1839         // 2. From points back to faces
1840
1841         const labelListList& pointFaces = mesh_.pointFaces();
1842
1843         forAll(pointFaces, pointI)
1844         {
1845             if (pointBaffle[pointI] != -1)
1846             {
1847                 const labelList& pFaces = pointFaces[pointI];
1848
1849                 forAll(pFaces, pFaceI)
1850                 {
1851                     label faceI = pFaces[pFaceI];
1852
1853                     if (ownPatch[faceI] == -1)
1854                     {
1855                         ownPatch[faceI] = pointBaffle[pointI];
1856                     }
1857                 }
1858             }
1859         }
1860         syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>(), false);
1861
1862
1863         // 3. From faces to cells (cellRegion) and back to faces (ownPatch)
1864
1865         labelList newOwnPatch(ownPatch);
1866
1867         forAll(ownPatch, faceI)
1868         {
1869             if (ownPatch[faceI] != -1)
1870             {
1871                 label own = faceOwner[faceI];
1872
1873                 if (cellRegion[own] != keepRegionI)
1874                 {
1875                     cellRegion[own] = keepRegionI;
1876
1877                     const cell& ownFaces = mesh_.cells()[own];
1878                     forAll(ownFaces, j)
1879                     {
1880                         if (ownPatch[ownFaces[j]] == -1)
1881                         {
1882                             newOwnPatch[ownFaces[j]] = ownPatch[faceI];
1883                         }
1884                     }
1885                 }
1886                 if (mesh_.isInternalFace(faceI))
1887                 {
1888                     label nei = faceNeighbour[faceI];
1889
1890                     if (cellRegion[nei] != keepRegionI)
1891                     {
1892                         cellRegion[nei] = keepRegionI;
1893
1894                         const cell& neiFaces = mesh_.cells()[nei];
1895                         forAll(neiFaces, j)
1896                         {
1897                             if (ownPatch[neiFaces[j]] == -1)
1898                             {
1899                                 newOwnPatch[neiFaces[j]] = ownPatch[faceI];
1900                             }
1901                         }
1902                     }
1903                 }
1904             }
1905         }
1906
1907         ownPatch.transfer(newOwnPatch);
1908
1909         syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>(), false);
1910     }
1911
1912
1913
1914     // Subset
1915     // ~~~~~~
1916
1917     // Get cells to remove
1918     DynamicList<label> cellsToRemove(mesh_.nCells());
1919     forAll(cellRegion, cellI)
1920     {
1921         if (cellRegion[cellI] != keepRegionI)
1922         {
1923             cellsToRemove.append(cellI);
1924         }
1925     }
1926     cellsToRemove.shrink();
1927
1928     label nCellsToKeep = mesh_.nCells() - cellsToRemove.size();
1929     reduce(nCellsToKeep, sumOp<label>());
1930
1931     Info<< "Keeping all cells in region " << keepRegionI
1932         << " containing point " << keepPoint << endl
1933         << "Selected for keeping : " << nCellsToKeep
1934         << " cells." << endl;
1935
1936
1937     // Remove cells
1938     removeCells cellRemover(mesh_);
1939
1940     // Pick up patches for exposed faces
1941     labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
1942     labelList exposedPatches(exposedFaces.size());
1943
1944     forAll(exposedFaces, i)
1945     {
1946         label faceI = exposedFaces[i];
1947
1948         if (ownPatch[faceI] != -1)
1949         {
1950             exposedPatches[i] = ownPatch[faceI];
1951         }
1952         else
1953         {
1954             WarningIn("meshRefinement::splitMesh(..)")
1955                 << "For exposed face " << faceI
1956                 << " fc:" << mesh_.faceCentres()[faceI]
1957                 << " found no patch." << endl
1958                 << "    Taking patch " << defaultPatch
1959                 << " instead." << endl;
1960             exposedPatches[i] = defaultPatch;
1961         }
1962     }
1963
1964     return doRemoveCells
1965     (
1966         cellsToRemove,
1967         exposedFaces,
1968         exposedPatches,
1969         cellRemover
1970     );
1971 }
1972
1973
1974 // Find boundary points that connect to more than one cell region and
1975 // split them.
1976 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints()
1977 {
1978     // Topochange container
1979     polyTopoChange meshMod(mesh_);
1980
1981
1982     // Analyse which points need to be duplicated
1983     localPointRegion regionSide(mesh_);
1984
1985     label nNonManifPoints = returnReduce
1986     (
1987         regionSide.meshPointMap().size(),
1988         sumOp<label>()
1989     );
1990
1991     Info<< "dupNonManifoldPoints : Found : " << nNonManifPoints
1992         << " non-manifold points (out of "
1993         << mesh_.globalData().nTotalPoints()
1994         << ')' << endl;
1995
1996     // Topo change engine
1997     duplicatePoints pointDuplicator(mesh_);
1998
1999     // Insert changes into meshMod
2000     pointDuplicator.setRefinement(regionSide, meshMod);
2001
2002     // Change the mesh (no inflation, parallel sync)
2003     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
2004
2005     // Update fields
2006     mesh_.updateMesh(map);
2007
2008     // Move mesh if in inflation mode
2009     if (map().hasMotionPoints())
2010     {
2011         mesh_.movePoints(map().preMotionPoints());
2012     }
2013     else
2014     {
2015         // Delete mesh volumes.
2016         mesh_.clearOut();
2017     }
2018
2019     if (overwrite())
2020     {
2021         mesh_.setInstance(oldInstance());
2022     }
2023
2024     // Update intersections. Is mapping only (no faces created, positions stay
2025     // same) so no need to recalculate intersections.
2026     updateMesh(map, labelList(0));
2027
2028     return map;
2029 }
2030
2031
2032 // Zoning
2033 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
2034 (
2035     const point& keepPoint
2036 )
2037 {
2038     const wordList& cellZoneNames = surfaces_.cellZoneNames();
2039     const wordList& faceZoneNames = surfaces_.faceZoneNames();
2040
2041     labelList namedSurfaces(surfaces_.getNamedSurfaces());
2042
2043     boolList isNamedSurface(cellZoneNames.size(), false);
2044
2045     forAll(namedSurfaces, i)
2046     {
2047         label surfI = namedSurfaces[i];
2048
2049         isNamedSurface[surfI] = true;
2050
2051         Info<< "Surface : " << surfaces_.names()[surfI] << nl
2052             << "    faceZone : " << faceZoneNames[surfI] << nl
2053             << "    cellZone : " << cellZoneNames[surfI] << endl;
2054     }
2055
2056
2057     // Add zones to mesh
2058
2059     labelList surfaceToFaceZone(faceZoneNames.size(), -1);
2060     {
2061         faceZoneMesh& faceZones = mesh_.faceZones();
2062
2063         forAll(namedSurfaces, i)
2064         {
2065             label surfI = namedSurfaces[i];
2066
2067             label zoneI = faceZones.findZoneID(faceZoneNames[surfI]);
2068
2069             if (zoneI == -1)
2070             {
2071                 zoneI = faceZones.size();
2072                 faceZones.setSize(zoneI+1);
2073                 faceZones.set
2074                 (
2075                     zoneI,
2076                     new faceZone
2077                     (
2078                         faceZoneNames[surfI],   //name
2079                         labelList(0),           //addressing
2080                         boolList(0),            //flipmap
2081                         zoneI,                  //index
2082                         faceZones               //faceZoneMesh
2083                     )
2084                 );
2085             }
2086
2087             if (debug)
2088             {
2089                 Pout<< "Faces on " << surfaces_.names()[surfI]
2090                     << " will go into faceZone " << zoneI << endl;
2091             }
2092             surfaceToFaceZone[surfI] = zoneI;
2093         }
2094
2095         // Check they are synced
2096         List<wordList> allFaceZones(Pstream::nProcs());
2097         allFaceZones[Pstream::myProcNo()] = faceZones.names();
2098         Pstream::gatherList(allFaceZones);
2099         Pstream::scatterList(allFaceZones);
2100
2101         for (label procI = 1; procI < allFaceZones.size(); procI++)
2102         {
2103             if (allFaceZones[procI] != allFaceZones[0])
2104             {
2105                 FatalErrorIn
2106                 (
2107                     "meshRefinement::zonify"
2108                     "(const label, const point&)"
2109                 )   << "Zones not synchronised among processors." << nl
2110                     << " Processor0 has faceZones:" << allFaceZones[0]
2111                     << " , processor" << procI
2112                     << " has faceZones:" << allFaceZones[procI]
2113                     << exit(FatalError);
2114             }
2115         }
2116     }
2117
2118     labelList surfaceToCellZone(cellZoneNames.size(), -1);
2119     {
2120         cellZoneMesh& cellZones = mesh_.cellZones();
2121
2122         forAll(namedSurfaces, i)
2123         {
2124             label surfI = namedSurfaces[i];
2125
2126             label zoneI = cellZones.findZoneID(cellZoneNames[surfI]);
2127
2128             if (zoneI == -1)
2129             {
2130                 zoneI = cellZones.size();
2131                 cellZones.setSize(zoneI+1);
2132                 cellZones.set
2133                 (
2134                     zoneI,
2135                     new cellZone
2136                     (
2137                         cellZoneNames[surfI],   //name
2138                         labelList(0),           //addressing
2139                         zoneI,                  //index
2140                         cellZones               //cellZoneMesh
2141                     )
2142                 );
2143             }
2144
2145             if (debug)
2146             {
2147                 Pout<< "Cells inside " << surfaces_.names()[surfI]
2148                     << " will go into cellZone " << zoneI << endl;
2149             }
2150             surfaceToCellZone[surfI] = zoneI;
2151         }
2152
2153         // Check they are synced
2154         List<wordList> allCellZones(Pstream::nProcs());
2155         allCellZones[Pstream::myProcNo()] = cellZones.names();
2156         Pstream::gatherList(allCellZones);
2157         Pstream::scatterList(allCellZones);
2158
2159         for (label procI = 1; procI < allCellZones.size(); procI++)
2160         {
2161             if (allCellZones[procI] != allCellZones[0])
2162             {
2163                 FatalErrorIn
2164                 (
2165                     "meshRefinement::zonify"
2166                     "(const label, const point&)"
2167                 )   << "Zones not synchronised among processors." << nl
2168                     << " Processor0 has cellZones:" << allCellZones[0]
2169                     << " , processor" << procI
2170                     << " has cellZones:" << allCellZones[procI]
2171                     << exit(FatalError);
2172             }
2173         }
2174     }
2175
2176
2177
2178     const pointField& cellCentres = mesh_.cellCentres();
2179     const labelList& faceOwner = mesh_.faceOwner();
2180     const labelList& faceNeighbour = mesh_.faceNeighbour();
2181     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
2182
2183
2184     // Mark faces intersecting zoned surfaces
2185     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2186
2187
2188     // Like surfaceIndex_ but only for named surfaces.
2189     labelList namedSurfaceIndex(mesh_.nFaces(), -1);
2190
2191     {
2192         // Statistics: number of faces per faceZone
2193         labelList nSurfFaces(faceZoneNames.size(), 0);
2194
2195         // Swap neighbouring cell centres and cell level
2196         labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
2197         pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
2198         calcNeighbourData(neiLevel, neiCc);
2199
2200         // Note: for all internal faces? internal + coupled?
2201         // Since zonify is run after baffling the surfaceIndex_ on baffles is
2202         // not synchronised across both baffle faces. Fortunately we don't
2203         // do zonify baffle faces anyway (they are normal boundary faces).
2204
2205         // Collect candidate faces
2206         // ~~~~~~~~~~~~~~~~~~~~~~~
2207
2208         labelList testFaces(intersectedFaces());
2209
2210         // Collect segments
2211         // ~~~~~~~~~~~~~~~~
2212
2213         pointField start(testFaces.size());
2214         pointField end(testFaces.size());
2215
2216         forAll(testFaces, i)
2217         {
2218             label faceI = testFaces[i];
2219
2220             if (mesh_.isInternalFace(faceI))
2221             {
2222                 start[i] = cellCentres[faceOwner[faceI]];
2223                 end[i] = cellCentres[faceNeighbour[faceI]];
2224             }
2225             else
2226             {
2227                 start[i] = cellCentres[faceOwner[faceI]];
2228                 end[i] = neiCc[faceI-mesh_.nInternalFaces()];
2229             }
2230         }
2231
2232
2233         // Do test for intersections
2234         // ~~~~~~~~~~~~~~~~~~~~~~~~~
2235         // Note that we intersect all intersected faces again. Could reuse
2236         // the information already in surfaceIndex_.
2237
2238         labelList surface1;
2239         labelList surface2;
2240         {
2241             List<pointIndexHit> hit1;
2242             labelList region1;
2243             List<pointIndexHit> hit2;
2244             labelList region2;
2245             surfaces_.findNearestIntersection
2246             (
2247                 namedSurfaces,
2248                 start,
2249                 end,
2250
2251                 surface1,
2252                 hit1,
2253                 region1,
2254                 surface2,
2255                 hit2,
2256                 region2
2257             );
2258         }
2259
2260         forAll(testFaces, i)
2261         {
2262             label faceI = testFaces[i];
2263
2264             if (surface1[i] != -1)
2265             {
2266                 // If both hit should probably choose nearest. For later.
2267                 namedSurfaceIndex[faceI] = surface1[i];
2268                 nSurfFaces[surface1[i]]++;
2269             }
2270             else if (surface2[i] != -1)
2271             {
2272                 namedSurfaceIndex[faceI] = surface2[i];
2273                 nSurfFaces[surface2[i]]++;
2274             }
2275         }
2276
2277
2278         // surfaceIndex migh have different surfaces on both sides if
2279         // there happen to be a (obviously thin) surface with different
2280         // regions between the cell centres. If one is on a named surface
2281         // and the other is not this might give problems so sync.
2282         syncTools::syncFaceList
2283         (
2284             mesh_,
2285             namedSurfaceIndex,
2286             maxEqOp<label>(),
2287             false
2288         );
2289
2290         // Print a bit
2291         if (debug)
2292         {
2293             forAll(nSurfFaces, surfI)
2294             {
2295                 Pout<< "Surface:"
2296                     << surfaces_.names()[surfI]
2297                     << "  nZoneFaces:" << nSurfFaces[surfI] << nl;
2298             }
2299             Pout<< endl;
2300         }
2301     }
2302
2303
2304     // Put the cells into the correct zone
2305     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2306
2307     // Closed surfaces with cellZone specified.
2308     labelList closedNamedSurfaces(surfaces_.getClosedNamedSurfaces());
2309
2310     // Zone per cell:
2311     // -2 : unset
2312     // -1 : not in any zone
2313     // >=0: zoneID
2314     labelList cellToZone(mesh_.nCells(), -2);
2315
2316
2317     // Set using geometric test
2318     // ~~~~~~~~~~~~~~~~~~~~~~~~
2319
2320     if (closedNamedSurfaces.size())
2321     {
2322         findCellZoneGeometric
2323         (
2324             closedNamedSurfaces,   // indices of closed surfaces
2325             namedSurfaceIndex,     // per face index of named surface
2326             surfaceToCellZone,     // cell zone index per surface
2327             cellToZone
2328         );
2329     }
2330
2331     // Set using walking
2332     // ~~~~~~~~~~~~~~~~~
2333
2334     //if (returnReduce(nSet, sumOp<label>()) < mesh_.globalData().nTotalCells())
2335     {
2336         // Topological walk
2337         findCellZoneTopo
2338         (
2339             keepPoint,
2340             namedSurfaceIndex,
2341             surfaceToCellZone,
2342             cellToZone
2343         );
2344     }
2345
2346
2347     //// Make sure namedSurfaceIndex is unset inbetween same cell cell zones.
2348     //makeConsistentFaceIndex(cellToZone, namedSurfaceIndex);
2349
2350
2351     // Topochange container
2352     polyTopoChange meshMod(mesh_);
2353
2354
2355     // Put the faces into the correct zone
2356     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2357
2358     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
2359     {
2360         label surfI = namedSurfaceIndex[faceI];
2361
2362         if (surfI != -1)
2363         {
2364             // Orient face zone to have slave cells in max cell zone.
2365             label ownZone = cellToZone[faceOwner[faceI]];
2366             label neiZone = cellToZone[faceNeighbour[faceI]];
2367
2368             bool flip;
2369             if (ownZone == max(ownZone, neiZone))
2370             {
2371                 flip = false;
2372             }
2373             else
2374             {
2375                 flip = true;
2376             }
2377
2378             meshMod.setAction
2379             (
2380                 polyModifyFace
2381                 (
2382                     mesh_.faces()[faceI],           // modified face
2383                     faceI,                          // label of face
2384                     faceOwner[faceI],               // owner
2385                     faceNeighbour[faceI],           // neighbour
2386                     false,                          // face flip
2387                     -1,                             // patch for face
2388                     false,                          // remove from zone
2389                     surfaceToFaceZone[surfI],       // zone for face
2390                     flip                            // face flip in zone
2391                 )
2392             );
2393         }
2394     }
2395
2396     // Get coupled neighbour cellZone. Set to -1 on non-coupled patches.
2397     labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
2398     forAll(patches, patchI)
2399     {
2400         const polyPatch& pp = patches[patchI];
2401
2402         if (pp.coupled())
2403         {
2404             forAll(pp, i)
2405             {
2406                 label faceI = pp.start()+i;
2407                 neiCellZone[faceI-mesh_.nInternalFaces()] =
2408                     cellToZone[mesh_.faceOwner()[faceI]];
2409             }
2410         }
2411     }
2412     syncTools::swapBoundaryFaceList(mesh_, neiCellZone, false);
2413
2414     // Set owner as no-flip
2415     forAll(patches, patchI)
2416     {
2417         const polyPatch& pp = patches[patchI];
2418
2419         label faceI = pp.start();
2420
2421         forAll(pp, i)
2422         {
2423             label surfI = namedSurfaceIndex[faceI];
2424
2425             if (surfI != -1)
2426             {
2427                 label ownZone = cellToZone[faceOwner[faceI]];
2428                 label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
2429
2430                 bool flip;
2431                 if (ownZone == max(ownZone, neiZone))
2432                 {
2433                     flip = false;
2434                 }
2435                 else
2436                 {
2437                     flip = true;
2438                 }
2439
2440                 meshMod.setAction
2441                 (
2442                     polyModifyFace
2443                     (
2444                         mesh_.faces()[faceI],           // modified face
2445                         faceI,                          // label of face
2446                         faceOwner[faceI],               // owner
2447                         -1,                             // neighbour
2448                         false,                          // face flip
2449                         patchI,                         // patch for face
2450                         false,                          // remove from zone
2451                         surfaceToFaceZone[surfI],       // zone for face
2452                         flip                            // face flip in zone
2453                     )
2454                 );
2455             }
2456             faceI++;
2457         }
2458     }
2459
2460
2461     // Put the cells into the correct zone
2462     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2463
2464     forAll(cellToZone, cellI)
2465     {
2466         label zoneI = cellToZone[cellI];
2467
2468         if (zoneI >= 0)
2469         {
2470             meshMod.setAction
2471             (
2472                 polyModifyCell
2473                 (
2474                     cellI,
2475                     false,          // removeFromZone
2476                     zoneI
2477                 )
2478             );
2479         }
2480     }
2481
2482     // Change the mesh (no inflation, parallel sync)
2483     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
2484
2485     // Update fields
2486     mesh_.updateMesh(map);
2487
2488     // Move mesh if in inflation mode
2489     if (map().hasMotionPoints())
2490     {
2491         mesh_.movePoints(map().preMotionPoints());
2492     }
2493     else
2494     {
2495         // Delete mesh volumes.
2496         mesh_.clearOut();
2497     }
2498
2499     if (overwrite())
2500     {
2501         mesh_.setInstance(oldInstance());
2502     }
2503
2504     // None of the faces has changed, only the zones. Still...
2505     updateMesh(map, labelList());
2506
2507     return map;
2508 }
2509
2510
2511 // ************************************************************************* //