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