src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.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 "OFstream.H"
  35 #include "cellSet.H"
  36 #include "searchableSurfaces.H"
  37 #include "polyMeshGeometry.H"
  38 #include "IOmanip.H"
  39 #include "unitConversion.H"
  40
  41 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  42
  43 void Foam::meshRefinement::markBoundaryFace
  44 (
  45     const label faceI,
  46     boolList& isBoundaryFace,
  47     boolList& isBoundaryEdge,
  48     boolList& isBoundaryPoint
  49 ) const
  50 {
  51     isBoundaryFace[faceI] = true;
  52
  53     const labelList& fEdges = mesh_.faceEdges(faceI);
  54
  55     forAll(fEdges, fp)
  56     {
  57         isBoundaryEdge[fEdges[fp]] = true;
  58     }
  59
  60     const face& f = mesh_.faces()[faceI];
  61
  62     forAll(f, fp)
  63     {
  64         isBoundaryPoint[f[fp]] = true;
  65     }
  66 }
  67
  68
  69 void Foam::meshRefinement::findNearest
  70 (
  71     const labelList& meshFaces,
  72     List<pointIndexHit>& nearestInfo,
  73     labelList& nearestSurface,
  74     labelList& nearestRegion,
  75     vectorField& nearestNormal
  76 ) const
  77 {
  78     pointField fc(meshFaces.size());
  79     forAll(meshFaces, i)
  80     {
  81         fc[i] = mesh_.faceCentres()[meshFaces[i]];
  82     }
  83
  84     const labelList allSurfaces(identity(surfaces_.surfaces().size()));
  85
  86     surfaces_.findNearest
  87     (
  88         allSurfaces,
  89         fc,
  90         scalarField(fc.size(), sqr(GREAT)),    // sqr of attraction
  91         nearestSurface,
  92         nearestInfo
  93     );
  94
  95     // Do normal testing per surface.
  96     nearestNormal.setSize(nearestInfo.size());
  97     nearestRegion.setSize(nearestInfo.size());
  98
  99     forAll(allSurfaces, surfI)
 100     {
 101         DynamicList<pointIndexHit> localHits;
 102
 103         forAll(nearestSurface, i)
 104         {
 105             if (nearestSurface[i] == surfI)
 106             {
 107                 localHits.append(nearestInfo[i]);
 108             }
 109         }
 110
 111         label geomI = surfaces_.surfaces()[surfI];
 112
 113         pointField localNormals;
 114         surfaces_.geometry()[geomI].getNormal(localHits, localNormals);
 115
 116         labelList localRegion;
 117         surfaces_.geometry()[geomI].getRegion(localHits, localRegion);
 118
 119         label localI = 0;
 120         forAll(nearestSurface, i)
 121         {
 122             if (nearestSurface[i] == surfI)
 123             {
 124                 nearestNormal[i] = localNormals[localI];
 125                 nearestRegion[i] = localRegion[localI];
 126                 localI++;
 127             }
 128         }
 129     }
 130 }
 131
 132
 133 Foam::Map<Foam::label> Foam::meshRefinement::findEdgeConnectedProblemCells
 134 (
 135     const scalarField& perpendicularAngle,
 136     const labelList& globalToPatch
 137 ) const
 138 {
 139     // Construct addressing engine from all patches added for meshing.
 140     autoPtr<indirectPrimitivePatch> ppPtr
 141     (
 142         meshRefinement::makePatch
 143         (
 144             mesh_,
 145             meshedPatches()
 146         )
 147     );
 148     const indirectPrimitivePatch& pp = ppPtr();
 149
 150
 151     // 1. Collect faces to test
 152     // ~~~~~~~~~~~~~~~~~~~~~~~~
 153
 154     DynamicList<label> candidateFaces(pp.size()/20);
 155
 156     const labelListList& edgeFaces = pp.edgeFaces();
 157
 158     const labelList& cellLevel = meshCutter_.cellLevel();
 159
 160     forAll(edgeFaces, edgeI)
 161     {
 162         const labelList& eFaces = edgeFaces[edgeI];
 163
 164         if (eFaces.size() == 2)
 165         {
 166             label face0 = pp.addressing()[eFaces[0]];
 167             label face1 = pp.addressing()[eFaces[1]];
 168
 169             label cell0 = mesh_.faceOwner()[face0];
 170             label cell1 = mesh_.faceOwner()[face1];
 171
 172             if (cellLevel[cell0] > cellLevel[cell1])
 173             {
 174                 // cell0 smaller.
 175                 const vector& n0 = pp.faceNormals()[eFaces[0]];
 176                 const vector& n1 = pp.faceNormals()[eFaces[1]];
 177
 178                 if (mag(n0 & n1) < 0.1)
 179                 {
 180                     candidateFaces.append(face0);
 181                 }
 182             }
 183             else if (cellLevel[cell1] > cellLevel[cell0])
 184             {
 185                 // cell1 smaller.
 186                 const vector& n0 = pp.faceNormals()[eFaces[0]];
 187                 const vector& n1 = pp.faceNormals()[eFaces[1]];
 188
 189                 if (mag(n0 & n1) < 0.1)
 190                 {
 191                     candidateFaces.append(face1);
 192                 }
 193             }
 194         }
 195     }
 196     candidateFaces.shrink();
 197
 198     Info<< "Testing " << returnReduce(candidateFaces.size(), sumOp<label>())
 199         << " faces on edge-connected cells of differing level."
 200         << endl;
 201
 202     if (debug)
 203     {
 204         faceSet fSet(mesh_, "edgeConnectedFaces", candidateFaces);
 205         fSet.instance() = timeName();
 206         Pout<< "Writing " << fSet.size()
 207             << " with problematic topology to faceSet "
 208             << fSet.objectPath() << endl;
 209         fSet.write();
 210     }
 211
 212
 213     // 2. Find nearest surface on candidate faces
 214     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 215
 216     List<pointIndexHit> nearestInfo;
 217     labelList nearestSurface;
 218     labelList nearestRegion;
 219     vectorField nearestNormal;
 220     findNearest
 221     (
 222         candidateFaces,
 223         nearestInfo,
 224         nearestSurface,
 225         nearestRegion,
 226         nearestNormal
 227     );
 228
 229
 230     // 3. Test angle to surface
 231     // ~~~~~~~~~~~~~~~~~~~~~~~~
 232
 233     Map<label> candidateCells(candidateFaces.size());
 234
 235     faceSet perpFaces(mesh_, "perpendicularFaces", pp.size()/100);
 236
 237     forAll(candidateFaces, i)
 238     {
 239         label faceI = candidateFaces[i];
 240
 241         vector n = mesh_.faceAreas()[faceI];
 242         n /= mag(n);
 243
 244         label region = surfaces_.globalRegion
 245         (
 246             nearestSurface[i],
 247             nearestRegion[i]
 248         );
 249
 250         scalar angle = degToRad(perpendicularAngle[region]);
 251
 252         if (angle >= 0)
 253         {
 254             if (mag(n & nearestNormal[i]) < Foam::sin(angle))
 255             {
 256                 perpFaces.insert(faceI);
 257                 candidateCells.insert
 258                 (
 259                     mesh_.faceOwner()[faceI],
 260                     globalToPatch[region]
 261                 );
 262             }
 263         }
 264     }
 265
 266     if (debug)
 267     {
 268         perpFaces.instance() = timeName();
 269         Pout<< "Writing " << perpFaces.size()
 270             << " faces that are perpendicular to the surface to set "
 271             << perpFaces.objectPath() << endl;
 272         perpFaces.write();
 273     }
 274     return candidateCells;
 275 }
 276
 277
 278 // Check if moving face to new points causes a 'collapsed' face.
 279 // Uses new point position only for the face, not the neighbouring
 280 // cell centres
 281 bool Foam::meshRefinement::isCollapsedFace
 282 (
 283     const pointField& points,
 284     const pointField& neiCc,
 285     const scalar minFaceArea,
 286     const scalar maxNonOrtho,
 287     const label faceI
 288 ) const
 289 {
 290     vector s = mesh_.faces()[faceI].normal(points);
 291     scalar magS = mag(s);
 292
 293     // Check face area
 294     if (magS < minFaceArea)
 295     {
 296         return true;
 297     }
 298
 299     // Check orthogonality
 300     const point& ownCc = mesh_.cellCentres()[mesh_.faceOwner()[faceI]];
 301
 302     if (mesh_.isInternalFace(faceI))
 303     {
 304         label nei = mesh_.faceNeighbour()[faceI];
 305         vector d = ownCc - mesh_.cellCentres()[nei];
 306
 307         scalar dDotS = (d & s)/(mag(d)*magS + VSMALL);
 308
 309         if (dDotS < maxNonOrtho)
 310         {
 311             return true;
 312         }
 313         else
 314         {
 315             return false;
 316         }
 317     }
 318     else
 319     {
 320         label patchI = mesh_.boundaryMesh().whichPatch(faceI);
 321
 322         if (mesh_.boundaryMesh()[patchI].coupled())
 323         {
 324             vector d = ownCc - neiCc[faceI-mesh_.nInternalFaces()];
 325
 326             scalar dDotS = (d & s)/(mag(d)*magS + VSMALL);
 327
 328             if (dDotS < maxNonOrtho)
 329             {
 330                 return true;
 331             }
 332             else
 333             {
 334                 return false;
 335             }
 336         }
 337         else
 338         {
 339             // Collapsing normal boundary face does not cause problems with
 340             // non-orthogonality
 341             return true;
 342         }
 343     }
 344 }
 345
 346
 347 // Check if moving cell to new points causes it to collapse.
 348 bool Foam::meshRefinement::isCollapsedCell
 349 (
 350     const pointField& points,
 351     const scalar volFraction,
 352     const label cellI
 353 ) const
 354 {
 355     scalar vol = mesh_.cells()[cellI].mag(points, mesh_.faces());
 356
 357     if (vol/mesh_.cellVolumes()[cellI] < volFraction)
 358     {
 359         return true;
 360     }
 361     else
 362     {
 363         return false;
 364     }
 365 }
 366
 367
 368 // Returns list with for every internal face -1 or the patch they should
 369 // be baffled into. Gets run after createBaffles so all the unzoned surface
 370 // intersections have already been turned into baffles. (Note: zoned surfaces
 371 // are not baffled at this stage)
 372 // Used to remove cells by baffling all their faces and have the
 373 // splitMeshRegions chuck away non used regions.
 374 Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
 375 (
 376     const dictionary& motionDict,
 377     const bool removeEdgeConnectedCells,
 378     const scalarField& perpendicularAngle,
 379     const labelList& globalToPatch
 380 ) const
 381 {
 382     const labelList& cellLevel = meshCutter_.cellLevel();
 383     const labelList& pointLevel = meshCutter_.pointLevel();
 384     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 385
 386     // Per internal face (boundary faces not used) the patch that the
 387     // baffle should get (or -1)
 388     labelList facePatch(mesh_.nFaces(), -1);
 389
 390     // Mark all points and edges on baffle patches (so not on any inlets,
 391     // outlets etc.)
 392     boolList isBoundaryPoint(mesh_.nPoints(), false);
 393     boolList isBoundaryEdge(mesh_.nEdges(), false);
 394     boolList isBoundaryFace(mesh_.nFaces(), false);
 395
 396     // Fill boundary data. All elements on meshed patches get marked.
 397     // Get the labels of added patches.
 398     labelList adaptPatchIDs(meshedPatches());
 399
 400     forAll(adaptPatchIDs, i)
 401     {
 402         label patchI = adaptPatchIDs[i];
 403
 404         const polyPatch& pp = patches[patchI];
 405
 406         label faceI = pp.start();
 407
 408         forAll(pp, j)
 409         {
 410             markBoundaryFace
 411             (
 412                 faceI,
 413                 isBoundaryFace,
 414                 isBoundaryEdge,
 415                 isBoundaryPoint
 416             );
 417
 418             faceI++;
 419         }
 420     }
 421
 422     // Swap neighbouring cell centres and cell level
 423     labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
 424     pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
 425     calcNeighbourData(neiLevel, neiCc);
 426
 427
 428     // Count of faces marked for baffling
 429     label nBaffleFaces = 0;
 430     PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
 431
 432     // Count of faces not baffled since would not cause a collapse
 433     label nPrevented = 0;
 434
 435     if (removeEdgeConnectedCells && max(perpendicularAngle) >= 0)
 436     {
 437         Info<< "markFacesOnProblemCells :"
 438             << " Checking for edge-connected cells of highly differing sizes."
 439             << endl;
 440
 441         // Pick up the cells that need to be removed and (a guess for)
 442         // the patch they should be patched with.
 443         Map<label> problemCells
 444         (
 445             findEdgeConnectedProblemCells
 446             (
 447                 perpendicularAngle,
 448                 globalToPatch
 449             )
 450         );
 451
 452         // Baffle all faces of cells that need to be removed
 453         forAllConstIter(Map<label>, problemCells, iter)
 454         {
 455             const cell& cFaces = mesh_.cells()[iter.key()];
 456
 457             forAll(cFaces, i)
 458             {
 459                 label faceI = cFaces[i];
 460
 461                 if (facePatch[faceI] == -1 && mesh_.isInternalFace(faceI))
 462                 {
 463                     facePatch[faceI] = getBafflePatch(facePatch, faceI);
 464                     nBaffleFaces++;
 465
 466                     // Mark face as a 'boundary'
 467                     markBoundaryFace
 468                     (
 469                         faceI,
 470                         isBoundaryFace,
 471                         isBoundaryEdge,
 472                         isBoundaryPoint
 473                     );
 474                 }
 475             }
 476         }
 477         Info<< "markFacesOnProblemCells : Marked "
 478             << returnReduce(nBaffleFaces, sumOp<label>())
 479             << " additional internal faces to be converted into baffles"
 480             << " due to "
 481             << returnReduce(problemCells.size(), sumOp<label>())
 482             << " cells edge-connected to lower level cells." << endl;
 483
 484         if (debug)
 485         {
 486             cellSet problemCellSet(mesh_, "problemCells", problemCells.toc());
 487             problemCellSet.instance() = timeName();
 488             Pout<< "Writing " << problemCellSet.size()
 489                 << " cells that are edge connected to coarser cell to set "
 490                 << problemCellSet.objectPath() << endl;
 491             problemCellSet.write();
 492         }
 493     }
 494
 495     syncTools::syncPointList
 496     (
 497         mesh_,
 498         isBoundaryPoint,
 499         orEqOp<bool>(),
 500         false               // null value
 501     );
 502
 503     syncTools::syncEdgeList
 504     (
 505         mesh_,
 506         isBoundaryEdge,
 507         orEqOp<bool>(),
 508         false               // null value
 509     );
 510
 511     syncTools::syncFaceList
 512     (
 513         mesh_,
 514         isBoundaryFace,
 515         orEqOp<bool>()
 516     );
 517
 518
 519     // See if checking for collapse
 520     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 521
 522     // Collapse checking parameters
 523     const scalar volFraction =
 524         motionDict.lookupOrDefault<scalar>("minVolCollapseRatio", -1);
 525
 526     const bool checkCollapse = (volFraction > 0);
 527     scalar minArea = -1;
 528     scalar maxNonOrtho = -1;
 529
 530
 531     // Find nearest (non-baffle) surface
 532     pointField newPoints;
 533
 534     if (checkCollapse)
 535     {
 536         minArea = readScalar(motionDict.lookup("minArea"));
 537         maxNonOrtho = readScalar(motionDict.lookup("maxNonOrtho"));
 538
 539         Info<< "markFacesOnProblemCells :"
 540             << " Deleting all-anchor surface cells only if"
 541             << "snapping them violates mesh quality constraints:" << nl
 542             << "    snapped/original cell volume < " << volFraction << nl
 543             << "    face area                    < " << minArea << nl
 544             << "    non-orthogonality            > " << maxNonOrtho << nl
 545             << endl;
 546
 547         // Construct addressing engine.
 548         autoPtr<indirectPrimitivePatch> ppPtr
 549         (
 550             meshRefinement::makePatch
 551             (
 552                 mesh_,
 553                 adaptPatchIDs
 554             )
 555         );
 556         const indirectPrimitivePatch& pp = ppPtr();
 557         const pointField& localPoints = pp.localPoints();
 558         const labelList& meshPoints = pp.meshPoints();
 559
 560         List<pointIndexHit> hitInfo;
 561         labelList hitSurface;
 562         surfaces_.findNearest
 563         (
 564             surfaces_.getUnnamedSurfaces(),
 565             localPoints,
 566             scalarField(localPoints.size(), sqr(GREAT)),    // sqr of attraction
 567             hitSurface,
 568             hitInfo
 569         );
 570
 571         // Start of from current points
 572         newPoints = mesh_.points();
 573
 574         forAll(hitInfo, i)
 575         {
 576             if (hitInfo[i].hit())
 577             {
 578                 newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
 579             }
 580         }
 581     }
 582
 583
 584
 585     // For each cell count the number of anchor points that are on
 586     // the boundary:
 587     // 8 : check the number of (baffle) boundary faces. If 3 or more block
 588     //     off the cell since the cell would get squeezed down to a diamond
 589     //     (probably; if the 3 or more faces are unrefined (only use the
 590     //      anchor points))
 591     // 7 : store. Used to check later on whether there are points with
 592     //     3 or more of these cells. (note that on a flat surface a boundary
 593     //     point will only have 4 cells connected to it)
 594
 595
 596     // Does cell have exactly 7 of its 8 anchor points on the boundary?
 597     PackedBoolList hasSevenBoundaryAnchorPoints(mesh_.nCells());
 598     // If so what is the remaining non-boundary anchor point?
 599     labelHashSet nonBoundaryAnchors(mesh_.nCells()/10000);
 600
 601     // On-the-fly addressing storage.
 602     DynamicList<label> dynFEdges;
 603     DynamicList<label> dynCPoints;
 604
 605     forAll(cellLevel, cellI)
 606     {
 607         const labelList& cPoints = mesh_.cellPoints(cellI, dynCPoints);
 608
 609         // Get number of anchor points (pointLevel <= cellLevel)
 610
 611         label nBoundaryAnchors = 0;
 612         label nNonAnchorBoundary = 0;
 613         label nonBoundaryAnchor = -1;
 614
 615         forAll(cPoints, i)
 616         {
 617             label pointI = cPoints[i];
 618
 619             if (pointLevel[pointI] <= cellLevel[cellI])
 620             {
 621                 // Anchor point
 622                 if (isBoundaryPoint[pointI])
 623                 {
 624                     nBoundaryAnchors++;
 625                 }
 626                 else
 627                 {
 628                     // Anchor point which is not on the surface
 629                     nonBoundaryAnchor = pointI;
 630                 }
 631             }
 632             else if (isBoundaryPoint[pointI])
 633             {
 634                 nNonAnchorBoundary++;
 635             }
 636         }
 637
 638         if (nBoundaryAnchors == 8)
 639         {
 640             const cell& cFaces = mesh_.cells()[cellI];
 641
 642             // Count boundary faces.
 643             label nBfaces = 0;
 644
 645             forAll(cFaces, cFaceI)
 646             {
 647                 if (isBoundaryFace[cFaces[cFaceI]])
 648                 {
 649                     nBfaces++;
 650                 }
 651             }
 652
 653             // If nBfaces > 1 make all non-boundary non-baffle faces baffles.
 654             // We assume that this situation is where there is a single
 655             // cell sticking out which would get flattened.
 656
 657             // Eugene: delete cell no matter what.
 658             //if (nBfaces > 1)
 659             {
 660                 if
 661                 (
 662                     checkCollapse
 663                 && !isCollapsedCell(newPoints, volFraction, cellI)
 664                 )
 665                 {
 666                     nPrevented++;
 667                     //Pout<< "Preventing collapse of 8 anchor point cell "
 668                     //    << cellI << " at " << mesh_.cellCentres()[cellI]
 669                     //    << " since new volume "
 670                     //    << mesh_.cells()[cellI].mag(newPoints, mesh_.faces())
 671                     //    << " old volume " << mesh_.cellVolumes()[cellI]
 672                     //    << endl;
 673                 }
 674                 else
 675                 {
 676                     // Block all faces of cell
 677                     forAll(cFaces, cf)
 678                     {
 679                         label faceI = cFaces[cf];
 680
 681                         if
 682                         (
 683                             facePatch[faceI] == -1
 684                          && mesh_.isInternalFace(faceI)
 685                         )
 686                         {
 687                             facePatch[faceI] = getBafflePatch(facePatch, faceI);
 688                             nBaffleFaces++;
 689
 690                             // Mark face as a 'boundary'
 691                             markBoundaryFace
 692                             (
 693                                 faceI,
 694                                 isBoundaryFace,
 695                                 isBoundaryEdge,
 696                                 isBoundaryPoint
 697                             );
 698                         }
 699                     }
 700                 }
 701             }
 702         }
 703         else if (nBoundaryAnchors == 7)
 704         {
 705             // Mark the cell. Store the (single!) non-boundary anchor point.
 706             hasSevenBoundaryAnchorPoints.set(cellI, 1u);
 707             nonBoundaryAnchors.insert(nonBoundaryAnchor);
 708         }
 709     }
 710
 711
 712     // Loop over all points. If a point is connected to 4 or more cells
 713     // with 7 anchor points on the boundary set those cell's non-boundary faces
 714     // to baffles
 715
 716     DynamicList<label> dynPCells;
 717
 718     forAllConstIter(labelHashSet, nonBoundaryAnchors, iter)
 719     {
 720         label pointI = iter.key();
 721
 722         const labelList& pCells = mesh_.pointCells(pointI, dynPCells);
 723
 724         // Count number of 'hasSevenBoundaryAnchorPoints' cells.
 725         label n = 0;
 726
 727         forAll(pCells, i)
 728         {
 729             if (hasSevenBoundaryAnchorPoints.get(pCells[i]) == 1u)
 730             {
 731                 n++;
 732             }
 733         }
 734
 735         if (n > 3)
 736         {
 737             // Point in danger of being what? Remove all 7-cells.
 738             forAll(pCells, i)
 739             {
 740                 label cellI = pCells[i];
 741
 742                 if (hasSevenBoundaryAnchorPoints.get(cellI) == 1u)
 743                 {
 744                     if
 745                     (
 746                         checkCollapse
 747                     && !isCollapsedCell(newPoints, volFraction, cellI)
 748                     )
 749                     {
 750                         nPrevented++;
 751                         //Pout<< "Preventing collapse of 7 anchor cell "
 752                         //    << cellI
 753                         //    << " at " << mesh_.cellCentres()[cellI]
 754                         //    << " since new volume "
 755                         //    << mesh_.cells()[cellI].mag
 756                         //        (newPoints, mesh_.faces())
 757                         //    << " old volume " << mesh_.cellVolumes()[cellI]
 758                         //    << endl;
 759                     }
 760                     else
 761                     {
 762                         const cell& cFaces = mesh_.cells()[cellI];
 763
 764                         forAll(cFaces, cf)
 765                         {
 766                             label faceI = cFaces[cf];
 767
 768                             if
 769                             (
 770                                 facePatch[faceI] == -1
 771                              && mesh_.isInternalFace(faceI)
 772                             )
 773                             {
 774                                 facePatch[faceI] = getBafflePatch
 775                                 (
 776                                     facePatch,
 777                                     faceI
 778                                 );
 779                                 nBaffleFaces++;
 780
 781                                 // Mark face as a 'boundary'
 782                                 markBoundaryFace
 783                                 (
 784                                     faceI,
 785                                     isBoundaryFace,
 786                                     isBoundaryEdge,
 787                                     isBoundaryPoint
 788                                 );
 789                             }
 790                         }
 791                     }
 792                 }
 793             }
 794         }
 795     }
 796
 797
 798     // Sync all. (note that pointdata and facedata not used anymore but sync
 799     // anyway)
 800
 801     syncTools::syncPointList
 802     (
 803         mesh_,
 804         isBoundaryPoint,
 805         orEqOp<bool>(),
 806         false               // null value
 807     );
 808
 809     syncTools::syncEdgeList
 810     (
 811         mesh_,
 812         isBoundaryEdge,
 813         orEqOp<bool>(),
 814         false               // null value
 815     );
 816
 817     syncTools::syncFaceList
 818     (
 819         mesh_,
 820         isBoundaryFace,
 821         orEqOp<bool>()
 822     );
 823
 824
 825     // Find faces with all edges on the boundary and make them baffles
 826     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
 827     {
 828         if (facePatch[faceI] == -1)
 829         {
 830             const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
 831             label nFaceBoundaryEdges = 0;
 832
 833             forAll(fEdges, fe)
 834             {
 835                 if (isBoundaryEdge[fEdges[fe]])
 836                 {
 837                     nFaceBoundaryEdges++;
 838                 }
 839             }
 840
 841             if (nFaceBoundaryEdges == fEdges.size())
 842             {
 843                 if
 844                 (
 845                     checkCollapse
 846                 && !isCollapsedFace
 847                     (
 848                         newPoints,
 849                         neiCc,
 850                         minArea,
 851                         maxNonOrtho,
 852                         faceI
 853                     )
 854                 )
 855                 {
 856                     nPrevented++;
 857                     //Pout<< "Preventing collapse of face " << faceI
 858                     //    << " with all boundary edges "
 859                     //    << " at " << mesh_.faceCentres()[faceI]
 860                     //    << endl;
 861                 }
 862                 else
 863                 {
 864                     facePatch[faceI] = getBafflePatch(facePatch, faceI);
 865                     nBaffleFaces++;
 866
 867                     // Do NOT update boundary data since this would grow blocked
 868                     // faces across gaps.
 869                 }
 870             }
 871         }
 872     }
 873
 874     forAll(patches, patchI)
 875     {
 876         const polyPatch& pp = patches[patchI];
 877
 878         if (pp.coupled())
 879         {
 880             label faceI = pp.start();
 881
 882             forAll(pp, i)
 883             {
 884                 if (facePatch[faceI] == -1)
 885                 {
 886                     const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
 887                     label nFaceBoundaryEdges = 0;
 888
 889                     forAll(fEdges, fe)
 890                     {
 891                         if (isBoundaryEdge[fEdges[fe]])
 892                         {
 893                             nFaceBoundaryEdges++;
 894                         }
 895                     }
 896
 897                     if (nFaceBoundaryEdges == fEdges.size())
 898                     {
 899                         if
 900                         (
 901                             checkCollapse
 902                         && !isCollapsedFace
 903                             (
 904                                 newPoints,
 905                                 neiCc,
 906                                 minArea,
 907                                 maxNonOrtho,
 908                                 faceI
 909                             )
 910                         )
 911                         {
 912                             nPrevented++;
 913                             //Pout<< "Preventing collapse of coupled face "
 914                             //    << faceI
 915                             //    << " with all boundary edges "
 916                             //    << " at " << mesh_.faceCentres()[faceI]
 917                             //    << endl;
 918                         }
 919                         else
 920                         {
 921                             facePatch[faceI] = getBafflePatch(facePatch, faceI);
 922                             if (isMasterFace[faceI])
 923                             {
 924                                 nBaffleFaces++;
 925                             }
 926
 927                             // Do NOT update boundary data since this would grow
 928                             // blocked faces across gaps.
 929                         }
 930                     }
 931                 }
 932
 933                 faceI++;
 934             }
 935         }
 936     }
 937
 938     Info<< "markFacesOnProblemCells : marked "
 939         << returnReduce(nBaffleFaces, sumOp<label>())
 940         << " additional internal faces to be converted into baffles."
 941         << endl;
 942
 943     if (checkCollapse)
 944     {
 945         Info<< "markFacesOnProblemCells : prevented "
 946             << returnReduce(nPrevented, sumOp<label>())
 947             << " internal faces fom getting converted into baffles."
 948             << endl;
 949     }
 950
 951     return facePatch;
 952 }
 953
 954
 955 //// Mark faces to be baffled to prevent snapping problems. Does
 956 //// test to find nearest surface and checks which faces would get squashed.
 957 //Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
 958 //(
 959 //    const dictionary& motionDict
 960 //) const
 961 //{
 962 //    // Construct addressing engine.
 963 //    autoPtr<indirectPrimitivePatch> ppPtr
 964 //    (
 965 //        meshRefinement::makePatch
 966 //        (
 967 //            mesh_,
 968 //            meshedPatches()
 969 //        )
 970 //    );
 971 //    const indirectPrimitivePatch& pp = ppPtr();
 972 //    const pointField& localPoints = pp.localPoints();
 973 //    const labelList& meshPoints = pp.meshPoints();
 974 //
 975 //    // Find nearest (non-baffle) surface
 976 //    pointField newPoints(mesh_.points());
 977 //    {
 978 //        List<pointIndexHit> hitInfo;
 979 //        labelList hitSurface;
 980 //        surfaces_.findNearest
 981 //        (
 982 //            surfaces_.getUnnamedSurfaces(),
 983 //            localPoints,
 984 //            scalarField(localPoints.size(), sqr(GREAT)),// sqr of attraction
 985 //            hitSurface,
 986 //            hitInfo
 987 //        );
 988 //
 989 //        forAll(hitInfo, i)
 990 //        {
 991 //            if (hitInfo[i].hit())
 992 //            {
 993 //                //label pointI = meshPoints[i];
 994 //                //Pout<< "   " << pointI << " moved from "
 995 //                //    << mesh_.points()[pointI] << " by "
 996 //                //    << mag(hitInfo[i].hitPoint()-mesh_.points()[pointI])
 997 //                //    << endl;
 998 //                newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
 999 //            }
1000 //        }
1001 //    }
1002 //
1003 //    // Per face (internal or coupled!) the patch that the
1004 //    // baffle should get (or -1).
1005 //    labelList facePatch(mesh_.nFaces(), -1);
1006 //    // Count of baffled faces
1007 //    label nBaffleFaces = 0;
1008 //
1009 //    {
1010 //        pointField oldPoints(mesh_.points());
1011 //        mesh_.movePoints(newPoints);
1012 //        faceSet wrongFaces(mesh_, "wrongFaces", 100);
1013 //        {
1014 //            //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
1015 //
1016 //            // Just check the errors from squashing
1017 //            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1018 //
1019 //            const labelList allFaces(identity(mesh_.nFaces()));
1020 //            label nWrongFaces = 0;
1021 //
1022 //            scalar minArea(readScalar(motionDict.lookup("minArea")));
1023 //            if (minArea > -SMALL)
1024 //            {
1025 //                polyMeshGeometry::checkFaceArea
1026 //                (
1027 //                    false,
1028 //                    minArea,
1029 //                    mesh_,
1030 //                    mesh_.faceAreas(),
1031 //                    allFaces,
1032 //                    &wrongFaces
1033 //                );
1034 //
1035 //                label nNewWrongFaces = returnReduce
1036 //                (
1037 //                    wrongFaces.size(),
1038 //                    sumOp<label>()
1039 //                );
1040 //
1041 //                Info<< "    faces with area < "
1042 //                    << setw(5) << minArea
1043 //                    << " m^2                            : "
1044 //                    << nNewWrongFaces-nWrongFaces << endl;
1045 //
1046 //                nWrongFaces = nNewWrongFaces;
1047 //            }
1048 //
1049 ////            scalar minDet(readScalar(motionDict.lookup("minDeterminant")));
1050 //            scalar minDet = 0.01;
1051 //            if (minDet > -1)
1052 //            {
1053 //                polyMeshGeometry::checkCellDeterminant
1054 //                (
1055 //                    false,
1056 //                    minDet,
1057 //                    mesh_,
1058 //                    mesh_.faceAreas(),
1059 //                    allFaces,
1060 //                    polyMeshGeometry::affectedCells(mesh_, allFaces),
1061 //                    &wrongFaces
1062 //                );
1063 //
1064 //                label nNewWrongFaces = returnReduce
1065 //                (
1066 //                    wrongFaces.size(),
1067 //                    sumOp<label>()
1068 //                );
1069 //
1070 //                Info<< "    faces on cells with determinant < "
1071 //                    << setw(5) << minDet << "                : "
1072 //                    << nNewWrongFaces-nWrongFaces << endl;
1073 //
1074 //                nWrongFaces = nNewWrongFaces;
1075 //            }
1076 //        }
1077 //
1078 //
1079 //        forAllConstIter(faceSet, wrongFaces, iter)
1080 //        {
1081 //            label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
1082 //
1083 //            if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
1084 //            {
1085 //                facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
1086 //                nBaffleFaces++;
1087 //
1088 //                //Pout<< "    " << iter.key()
1089 //                //    //<< " on patch " << mesh_.boundaryMesh()[patchI].name()
1090 //                //    << " is destined for patch " << facePatch[iter.key()]
1091 //                //    << endl;
1092 //            }
1093 //        }
1094 //        // Restore points.
1095 //        mesh_.movePoints(oldPoints);
1096 //    }
1097 //
1098 //
1099 //    Info<< "markFacesOnProblemCellsGeometric : marked "
1100 //        << returnReduce(nBaffleFaces, sumOp<label>())
1101 //        << " additional internal and coupled faces"
1102 //        << " to be converted into baffles." << endl;
1103 //
1104 //    syncTools::syncFaceList
1105 //    (
1106 //        mesh_,
1107 //        facePatch,
1108 //        maxEqOp<label>()
1109 //    );
1110 //
1111 //    return facePatch;
1112 //}
1113
1114
1115 // ************************************************************************* //