src/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*----------------------------------------------------------------------------*/
  26
  27 #include "meshRefinement.H"
  28 #include "fvMesh.H"
  29 #include "syncTools.H"
  30 #include "Time.H"
  31 #include "refinementSurfaces.H"
  32 #include "pointSet.H"
  33 #include "faceSet.H"
  34 #include "indirectPrimitivePatch.H"
  35 #include "OFstream.H"
  36 #include "cellSet.H"
  37 #include "searchableSurfaces.H"
  38 #include "polyMeshGeometry.H"
  39 #include "IOmanip.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         Pout<< "Writing " << fSet.size()
 206             << " with problematic topology to faceSet "
 207             << fSet.objectPath() << endl;
 208         fSet.write();
 209     }
 210
 211
 212     // 2. Find nearest surface on candidate faces
 213     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 214
 215     List<pointIndexHit> nearestInfo;
 216     labelList nearestSurface;
 217     labelList nearestRegion;
 218     vectorField nearestNormal;
 219     findNearest
 220     (
 221         candidateFaces,
 222         nearestInfo,
 223         nearestSurface,
 224         nearestRegion,
 225         nearestNormal
 226     );
 227
 228
 229     // 3. Test angle to surface
 230     // ~~~~~~~~~~~~~~~~~~~~~~~~
 231
 232     Map<label> candidateCells(candidateFaces.size());
 233
 234     faceSet perpFaces(mesh_, "perpendicularFaces", pp.size()/100);
 235
 236     forAll(candidateFaces, i)
 237     {
 238         label faceI = candidateFaces[i];
 239
 240         vector n = mesh_.faceAreas()[faceI];
 241         n /= mag(n);
 242
 243         label region = surfaces_.globalRegion
 244         (
 245             nearestSurface[i],
 246             nearestRegion[i]
 247         );
 248
 249         scalar angle =
 250             perpendicularAngle[region]
 251           / 180.0
 252           * mathematicalConstant::pi;
 253
 254         if (angle >= 0)
 255         {
 256             if (mag(n & nearestNormal[i]) < Foam::sin(angle))
 257             {
 258                 perpFaces.insert(faceI);
 259                 candidateCells.insert
 260                 (
 261                     mesh_.faceOwner()[faceI],
 262                     globalToPatch[region]
 263                 );
 264             }
 265         }
 266     }
 267
 268     if (debug)
 269     {
 270         Pout<< "Writing " << perpFaces.size()
 271             << " faces that are perpendicular to the surface to set "
 272             << perpFaces.objectPath() << endl;
 273         perpFaces.write();
 274     }
 275     return candidateCells;
 276 }
 277
 278
 279 // Check if moving face to new points causes a 'collapsed' face.
 280 // Uses new point position only for the face, not the neighbouring
 281 // cell centres
 282 bool Foam::meshRefinement::isCollapsedFace
 283 (
 284     const pointField& points,
 285     const pointField& neiCc,
 286     const scalar minFaceArea,
 287     const scalar maxNonOrtho,
 288     const label faceI
 289 ) const
 290 {
 291     vector s = mesh_.faces()[faceI].normal(points);
 292     scalar magS = mag(s);
 293
 294     // Check face area
 295     if (magS < minFaceArea)
 296     {
 297         return true;
 298     }
 299
 300     // Check orthogonality
 301     const point& ownCc = mesh_.cellCentres()[mesh_.faceOwner()[faceI]];
 302
 303     if (mesh_.isInternalFace(faceI))
 304     {
 305         label nei = mesh_.faceNeighbour()[faceI];
 306         vector d = ownCc - mesh_.cellCentres()[nei];
 307
 308         scalar dDotS = (d & s)/(mag(d)*magS + VSMALL);
 309
 310         if (dDotS < maxNonOrtho)
 311         {
 312             return true;
 313         }
 314         else
 315         {
 316             return false;
 317         }
 318     }
 319     else
 320     {
 321         label patchI = mesh_.boundaryMesh().whichPatch(faceI);
 322
 323         if (mesh_.boundaryMesh()[patchI].coupled())
 324         {
 325             vector d = ownCc - neiCc[faceI-mesh_.nInternalFaces()];
 326
 327             scalar dDotS = (d & s)/(mag(d)*magS + VSMALL);
 328
 329             if (dDotS < maxNonOrtho)
 330             {
 331                 return true;
 332             }
 333             else
 334             {
 335                 return false;
 336             }
 337         }
 338         else
 339         {
 340             // Collapsing normal boundary face does not cause problems with
 341             // non-orthogonality
 342             return true;
 343         }
 344     }
 345 }
 346
 347
 348 // Check if moving cell to new points causes it to collapse.
 349 bool Foam::meshRefinement::isCollapsedCell
 350 (
 351     const pointField& points,
 352     const scalar volFraction,
 353     const label cellI
 354 ) const
 355 {
 356     scalar vol = mesh_.cells()[cellI].mag(points, mesh_.faces());
 357
 358     if (vol/mesh_.cellVolumes()[cellI] < volFraction)
 359     {
 360         return true;
 361     }
 362     else
 363     {
 364         return false;
 365     }
 366 }
 367
 368
 369 // Returns list with for every internal face -1 or the patch they should
 370 // be baffled into. Gets run after createBaffles so all the unzoned surface
 371 // intersections have already been turned into baffles. (Note: zoned surfaces
 372 // are not baffled at this stage)
 373 // Used to remove cells by baffling all their faces and have the
 374 // splitMeshRegions chuck away non used regions.
 375 Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
 376 (
 377     const dictionary& motionDict,
 378     const bool removeEdgeConnectedCells,
 379     const scalarField& perpendicularAngle,
 380     const labelList& globalToPatch
 381 ) const
 382 {
 383     const labelList& cellLevel = meshCutter_.cellLevel();
 384     const labelList& pointLevel = meshCutter_.pointLevel();
 385     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 386
 387     // Per internal face (boundary faces not used) the patch that the
 388     // baffle should get (or -1)
 389     labelList facePatch(mesh_.nFaces(), -1);
 390
 391     // Mark all points and edges on baffle patches (so not on any inlets,
 392     // outlets etc.)
 393     boolList isBoundaryPoint(mesh_.nPoints(), false);
 394     boolList isBoundaryEdge(mesh_.nEdges(), false);
 395     boolList isBoundaryFace(mesh_.nFaces(), false);
 396
 397     // Fill boundary data. All elements on meshed patches get marked.
 398     // Get the labels of added patches.
 399     labelList adaptPatchIDs(meshedPatches());
 400
 401     forAll(adaptPatchIDs, i)
 402     {
 403         label patchI = adaptPatchIDs[i];
 404
 405         const polyPatch& pp = patches[patchI];
 406
 407         label faceI = pp.start();
 408
 409         forAll(pp, j)
 410         {
 411             markBoundaryFace
 412             (
 413                 faceI,
 414                 isBoundaryFace,
 415                 isBoundaryEdge,
 416                 isBoundaryPoint
 417             );
 418
 419             faceI++;
 420         }
 421     }
 422
 423     // Swap neighbouring cell centres and cell level
 424     labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
 425     pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
 426     calcNeighbourData(neiLevel, neiCc);
 427
 428
 429     // Count of faces marked for baffling
 430     label nBaffleFaces = 0;
 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             Pout<< "Writing " << problemCellSet.size()
 488                 << " cells that are edge connected to coarser cell to set "
 489                 << problemCellSet.objectPath() << endl;
 490             problemCellSet.write();
 491         }
 492     }
 493
 494     syncTools::syncPointList
 495     (
 496         mesh_,
 497         isBoundaryPoint,
 498         orEqOp<bool>(),
 499         false,              // null value
 500         false               // no separation
 501     );
 502
 503     syncTools::syncEdgeList
 504     (
 505         mesh_,
 506         isBoundaryEdge,
 507         orEqOp<bool>(),
 508         false,              // null value
 509         false               // no separation
 510     );
 511
 512     syncTools::syncFaceList
 513     (
 514         mesh_,
 515         isBoundaryFace,
 516         orEqOp<bool>(),
 517         false               // no separation
 518     );
 519
 520
 521     // See if checking for collapse
 522     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 523
 524     // Collapse checking parameters
 525     scalar volFraction = -1;
 526     if (motionDict.found("minVolCollapseRatio"))
 527     {
 528         volFraction = readScalar(motionDict.lookup("minVolCollapseRatio"));
 529     }
 530     const bool checkCollapse = (volFraction > 0);
 531     scalar minArea = -1;
 532     scalar maxNonOrtho = -1;
 533
 534
 535     // Find nearest (non-baffle) surface
 536     pointField newPoints;
 537
 538     if (checkCollapse)
 539     {
 540         minArea = readScalar(motionDict.lookup("minArea"));
 541         maxNonOrtho = readScalar(motionDict.lookup("maxNonOrtho"));
 542
 543         Info<< "markFacesOnProblemCells :"
 544             << " Deleting all-anchor surface cells only if"
 545             << "snapping them violates mesh quality constraints:" << nl
 546             << "    snapped/original cell volume < " << volFraction << nl
 547             << "    face area                    < " << minArea << nl
 548             << "    non-orthogonality            > " << maxNonOrtho << nl
 549             << endl;
 550
 551         // Construct addressing engine.
 552         autoPtr<indirectPrimitivePatch> ppPtr
 553         (
 554             meshRefinement::makePatch
 555             (
 556                 mesh_,
 557                 adaptPatchIDs
 558             )
 559         );
 560         const indirectPrimitivePatch& pp = ppPtr();
 561         const pointField& localPoints = pp.localPoints();
 562         const labelList& meshPoints = pp.meshPoints();
 563
 564         List<pointIndexHit> hitInfo;
 565         labelList hitSurface;
 566         surfaces_.findNearest
 567         (
 568             surfaces_.getUnnamedSurfaces(),
 569             localPoints,
 570             scalarField(localPoints.size(), sqr(GREAT)),    // sqr of attraction
 571             hitSurface,
 572             hitInfo
 573         );
 574
 575         // Start of from current points
 576         newPoints = mesh_.points();
 577
 578         forAll(hitInfo, i)
 579         {
 580             if (hitInfo[i].hit())
 581             {
 582                 newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
 583             }
 584         }
 585     }
 586
 587
 588
 589     // For each cell count the number of anchor points that are on
 590     // the boundary:
 591     // 8 : check the number of (baffle) boundary faces. If 3 or more block
 592     //     off the cell since the cell would get squeezed down to a diamond
 593     //     (probably; if the 3 or more faces are unrefined (only use the
 594     //      anchor points))
 595     // 7 : store. Used to check later on whether there are points with
 596     //     3 or more of these cells. (note that on a flat surface a boundary
 597     //     point will only have 4 cells connected to it)
 598
 599
 600     // Does cell have exactly 7 of its 8 anchor points on the boundary?
 601     PackedBoolList hasSevenBoundaryAnchorPoints(mesh_.nCells());
 602     // If so what is the remaining non-boundary anchor point?
 603     labelHashSet nonBoundaryAnchors(mesh_.nCells()/10000);
 604
 605     // On-the-fly addressing storage.
 606     DynamicList<label> dynFEdges;
 607     DynamicList<label> dynCPoints;
 608
 609     forAll(cellLevel, cellI)
 610     {
 611         const labelList& cPoints = mesh_.cellPoints(cellI, dynCPoints);
 612
 613         // Get number of anchor points (pointLevel <= cellLevel)
 614
 615         label nBoundaryAnchors = 0;
 616         label nNonAnchorBoundary = 0;
 617         label nonBoundaryAnchor = -1;
 618
 619         forAll(cPoints, i)
 620         {
 621             label pointI = cPoints[i];
 622
 623             if (pointLevel[pointI] <= cellLevel[cellI])
 624             {
 625                 // Anchor point
 626                 if (isBoundaryPoint[pointI])
 627                 {
 628                     nBoundaryAnchors++;
 629                 }
 630                 else
 631                 {
 632                     // Anchor point which is not on the surface
 633                     nonBoundaryAnchor = pointI;
 634                 }
 635             }
 636             else if (isBoundaryPoint[pointI])
 637             {
 638                 nNonAnchorBoundary++;
 639             }
 640         }
 641
 642         if (nBoundaryAnchors == 8)
 643         {
 644             const cell& cFaces = mesh_.cells()[cellI];
 645
 646             // Count boundary faces.
 647             label nBfaces = 0;
 648
 649             forAll(cFaces, cFaceI)
 650             {
 651                 if (isBoundaryFace[cFaces[cFaceI]])
 652                 {
 653                     nBfaces++;
 654                 }
 655             }
 656
 657             // If nBfaces > 1 make all non-boundary non-baffle faces baffles.
 658             // We assume that this situation is where there is a single
 659             // cell sticking out which would get flattened.
 660
 661             // Eugene: delete cell no matter what.
 662             //if (nBfaces > 1)
 663             {
 664                 if
 665                 (
 666                     checkCollapse
 667                 && !isCollapsedCell(newPoints, volFraction, cellI)
 668                 )
 669                 {
 670                     nPrevented++;
 671                     //Pout<< "Preventing collapse of 8 anchor point cell "
 672                     //    << cellI << " at " << mesh_.cellCentres()[cellI]
 673                     //    << " since new volume "
 674                     //    << mesh_.cells()[cellI].mag(newPoints, mesh_.faces())
 675                     //    << " old volume " << mesh_.cellVolumes()[cellI]
 676                     //    << endl;
 677                 }
 678                 else
 679                 {
 680                     // Block all faces of cell
 681                     forAll(cFaces, cf)
 682                     {
 683                         label faceI = cFaces[cf];
 684
 685                         if
 686                         (
 687                             facePatch[faceI] == -1
 688                          && mesh_.isInternalFace(faceI)
 689                         )
 690                         {
 691                             facePatch[faceI] = getBafflePatch(facePatch, faceI);
 692                             nBaffleFaces++;
 693
 694                             // Mark face as a 'boundary'
 695                             markBoundaryFace
 696                             (
 697                                 faceI,
 698                                 isBoundaryFace,
 699                                 isBoundaryEdge,
 700                                 isBoundaryPoint
 701                             );
 702                         }
 703                     }
 704                 }
 705             }
 706         }
 707         else if (nBoundaryAnchors == 7)
 708         {
 709             // Mark the cell. Store the (single!) non-boundary anchor point.
 710             hasSevenBoundaryAnchorPoints.set(cellI, 1u);
 711             nonBoundaryAnchors.insert(nonBoundaryAnchor);
 712         }
 713     }
 714
 715
 716     // Loop over all points. If a point is connected to 4 or more cells
 717     // with 7 anchor points on the boundary set those cell's non-boundary faces
 718     // to baffles
 719
 720     DynamicList<label> dynPCells;
 721
 722     forAllConstIter(labelHashSet, nonBoundaryAnchors, iter)
 723     {
 724         label pointI = iter.key();
 725
 726         const labelList& pCells = mesh_.pointCells(pointI, dynPCells);
 727
 728         // Count number of 'hasSevenBoundaryAnchorPoints' cells.
 729         label n = 0;
 730
 731         forAll(pCells, i)
 732         {
 733             if (hasSevenBoundaryAnchorPoints.get(pCells[i]) == 1u)
 734             {
 735                 n++;
 736             }
 737         }
 738
 739         if (n > 3)
 740         {
 741             // Point in danger of being what? Remove all 7-cells.
 742             forAll(pCells, i)
 743             {
 744                 label cellI = pCells[i];
 745
 746                 if (hasSevenBoundaryAnchorPoints.get(cellI) == 1u)
 747                 {
 748                     if
 749                     (
 750                         checkCollapse
 751                     && !isCollapsedCell(newPoints, volFraction, cellI)
 752                     )
 753                     {
 754                         nPrevented++;
 755                         //Pout<< "Preventing collapse of 7 anchor cell "
 756                         //    << cellI
 757                         //    << " at " << mesh_.cellCentres()[cellI]
 758                         //    << " since new volume "
 759                         //    << mesh_.cells()[cellI].mag
 760                         //        (newPoints, mesh_.faces())
 761                         //    << " old volume " << mesh_.cellVolumes()[cellI]
 762                         //    << endl;
 763                     }
 764                     else
 765                     {
 766                         const cell& cFaces = mesh_.cells()[cellI];
 767
 768                         forAll(cFaces, cf)
 769                         {
 770                             label faceI = cFaces[cf];
 771
 772                             if
 773                             (
 774                                 facePatch[faceI] == -1
 775                              && mesh_.isInternalFace(faceI)
 776                             )
 777                             {
 778                                 facePatch[faceI] = getBafflePatch
 779                                 (
 780                                     facePatch,
 781                                     faceI
 782                                 );
 783                                 nBaffleFaces++;
 784
 785                                 // Mark face as a 'boundary'
 786                                 markBoundaryFace
 787                                 (
 788                                     faceI,
 789                                     isBoundaryFace,
 790                                     isBoundaryEdge,
 791                                     isBoundaryPoint
 792                                 );
 793                             }
 794                         }
 795                     }
 796                 }
 797             }
 798         }
 799     }
 800
 801
 802     // Sync all. (note that pointdata and facedata not used anymore but sync
 803     // anyway)
 804
 805     syncTools::syncPointList
 806     (
 807         mesh_,
 808         isBoundaryPoint,
 809         orEqOp<bool>(),
 810         false,              // null value
 811         false               // no separation
 812     );
 813
 814     syncTools::syncEdgeList
 815     (
 816         mesh_,
 817         isBoundaryEdge,
 818         orEqOp<bool>(),
 819         false,              // null value
 820         false               // no separation
 821     );
 822
 823     syncTools::syncFaceList
 824     (
 825         mesh_,
 826         isBoundaryFace,
 827         orEqOp<bool>(),
 828         false               // no separation
 829     );
 830
 831
 832     // Find faces with all edges on the boundary and make them baffles
 833     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
 834     {
 835         if (facePatch[faceI] == -1)
 836         {
 837             const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
 838             label nFaceBoundaryEdges = 0;
 839
 840             forAll(fEdges, fe)
 841             {
 842                 if (isBoundaryEdge[fEdges[fe]])
 843                 {
 844                     nFaceBoundaryEdges++;
 845                 }
 846             }
 847
 848             if (nFaceBoundaryEdges == fEdges.size())
 849             {
 850                 if
 851                 (
 852                     checkCollapse
 853                 && !isCollapsedFace
 854                     (
 855                         newPoints,
 856                         neiCc,
 857                         minArea,
 858                         maxNonOrtho,
 859                         faceI
 860                     )
 861                 )
 862                 {
 863                     nPrevented++;
 864                     //Pout<< "Preventing collapse of face " << faceI
 865                     //    << " with all boundary edges "
 866                     //    << " at " << mesh_.faceCentres()[faceI]
 867                     //    << endl;
 868                 }
 869                 else
 870                 {
 871                     facePatch[faceI] = getBafflePatch(facePatch, faceI);
 872                     nBaffleFaces++;
 873
 874                     // Do NOT update boundary data since this would grow blocked
 875                     // faces across gaps.
 876                 }
 877             }
 878         }
 879     }
 880
 881     forAll(patches, patchI)
 882     {
 883         const polyPatch& pp = patches[patchI];
 884
 885         if (pp.coupled())
 886         {
 887             label faceI = pp.start();
 888
 889             forAll(pp, i)
 890             {
 891                 if (facePatch[faceI] == -1)
 892                 {
 893                     const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
 894                     label nFaceBoundaryEdges = 0;
 895
 896                     forAll(fEdges, fe)
 897                     {
 898                         if (isBoundaryEdge[fEdges[fe]])
 899                         {
 900                             nFaceBoundaryEdges++;
 901                         }
 902                     }
 903
 904                     if (nFaceBoundaryEdges == fEdges.size())
 905                     {
 906                         if
 907                         (
 908                             checkCollapse
 909                         && !isCollapsedFace
 910                             (
 911                                 newPoints,
 912                                 neiCc,
 913                                 minArea,
 914                                 maxNonOrtho,
 915                                 faceI
 916                             )
 917                         )
 918                         {
 919                             nPrevented++;
 920                             //Pout<< "Preventing collapse of coupled face "
 921                             //    << faceI
 922                             //    << " with all boundary edges "
 923                             //    << " at " << mesh_.faceCentres()[faceI]
 924                             //    << endl;
 925                         }
 926                         else
 927                         {
 928                             facePatch[faceI] = getBafflePatch(facePatch, faceI);
 929                             nBaffleFaces++;
 930
 931                             // Do NOT update boundary data since this would grow
 932                             // blocked faces across gaps.
 933                         }
 934                     }
 935                 }
 936
 937                 faceI++;
 938             }
 939         }
 940     }
 941
 942     Info<< "markFacesOnProblemCells : marked "
 943         << returnReduce(nBaffleFaces, sumOp<label>())
 944         << " additional internal faces to be converted into baffles."
 945         << endl;
 946
 947     if (checkCollapse)
 948     {
 949         Info<< "markFacesOnProblemCells : prevented "
 950             << returnReduce(nPrevented, sumOp<label>())
 951             << " internal faces fom getting converted into baffles."
 952             << endl;
 953     }
 954
 955     return facePatch;
 956 }
 957
 958
 959 //// Mark faces to be baffled to prevent snapping problems. Does
 960 //// test to find nearest surface and checks which faces would get squashed.
 961 //Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
 962 //(
 963 //    const dictionary& motionDict
 964 //) const
 965 //{
 966 //    // Construct addressing engine.
 967 //    autoPtr<indirectPrimitivePatch> ppPtr
 968 //    (
 969 //        meshRefinement::makePatch
 970 //        (
 971 //            mesh_,
 972 //            meshedPatches()
 973 //        )
 974 //    );
 975 //    const indirectPrimitivePatch& pp = ppPtr();
 976 //    const pointField& localPoints = pp.localPoints();
 977 //    const labelList& meshPoints = pp.meshPoints();
 978 //
 979 //    // Find nearest (non-baffle) surface
 980 //    pointField newPoints(mesh_.points());
 981 //    {
 982 //        List<pointIndexHit> hitInfo;
 983 //        labelList hitSurface;
 984 //        surfaces_.findNearest
 985 //        (
 986 //            surfaces_.getUnnamedSurfaces(),
 987 //            localPoints,
 988 //            scalarField(localPoints.size(), sqr(GREAT)),// sqr of attraction
 989 //            hitSurface,
 990 //            hitInfo
 991 //        );
 992 //
 993 //        forAll(hitInfo, i)
 994 //        {
 995 //            if (hitInfo[i].hit())
 996 //            {
 997 //                //label pointI = meshPoints[i];
 998 //                //Pout<< "   " << pointI << " moved from "
 999 //                //    << mesh_.points()[pointI] << " by "
1000 //                //    << mag(hitInfo[i].hitPoint()-mesh_.points()[pointI])
1001 //                //    << endl;
1002 //                newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
1003 //            }
1004 //        }
1005 //    }
1006 //
1007 //    // Per face (internal or coupled!) the patch that the
1008 //    // baffle should get (or -1).
1009 //    labelList facePatch(mesh_.nFaces(), -1);
1010 //    // Count of baffled faces
1011 //    label nBaffleFaces = 0;
1012 //
1013 //    {
1014 //        pointField oldPoints(mesh_.points());
1015 //        mesh_.movePoints(newPoints);
1016 //        faceSet wrongFaces(mesh_, "wrongFaces", 100);
1017 //        {
1018 //            //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
1019 //
1020 //            // Just check the errors from squashing
1021 //            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1022 //
1023 //            const labelList allFaces(identity(mesh_.nFaces()));
1024 //            label nWrongFaces = 0;
1025 //
1026 //            scalar minArea(readScalar(motionDict.lookup("minArea")));
1027 //            if (minArea > -SMALL)
1028 //            {
1029 //                polyMeshGeometry::checkFaceArea
1030 //                (
1031 //                    false,
1032 //                    minArea,
1033 //                    mesh_,
1034 //                    mesh_.faceAreas(),
1035 //                    allFaces,
1036 //                    &wrongFaces
1037 //                );
1038 //
1039 //                label nNewWrongFaces = returnReduce
1040 //                (
1041 //                    wrongFaces.size(),
1042 //                    sumOp<label>()
1043 //                );
1044 //
1045 //                Info<< "    faces with area < "
1046 //                    << setw(5) << minArea
1047 //                    << " m^2                            : "
1048 //                    << nNewWrongFaces-nWrongFaces << endl;
1049 //
1050 //                nWrongFaces = nNewWrongFaces;
1051 //            }
1052 //
1053 ////            scalar minDet(readScalar(motionDict.lookup("minDeterminant")));
1054 //            scalar minDet = 0.01;
1055 //            if (minDet > -1)
1056 //            {
1057 //                polyMeshGeometry::checkCellDeterminant
1058 //                (
1059 //                    false,
1060 //                    minDet,
1061 //                    mesh_,
1062 //                    mesh_.faceAreas(),
1063 //                    allFaces,
1064 //                    polyMeshGeometry::affectedCells(mesh_, allFaces),
1065 //                    &wrongFaces
1066 //                );
1067 //
1068 //                label nNewWrongFaces = returnReduce
1069 //                (
1070 //                    wrongFaces.size(),
1071 //                    sumOp<label>()
1072 //                );
1073 //
1074 //                Info<< "    faces on cells with determinant < "
1075 //                    << setw(5) << minDet << "                : "
1076 //                    << nNewWrongFaces-nWrongFaces << endl;
1077 //
1078 //                nWrongFaces = nNewWrongFaces;
1079 //            }
1080 //        }
1081 //
1082 //
1083 //        forAllConstIter(faceSet, wrongFaces, iter)
1084 //        {
1085 //            label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
1086 //
1087 //            if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
1088 //            {
1089 //                facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
1090 //                nBaffleFaces++;
1091 //
1092 //                //Pout<< "    " << iter.key()
1093 //                //    //<< " on patch " << mesh_.boundaryMesh()[patchI].name()
1094 //                //    << " is destined for patch " << facePatch[iter.key()]
1095 //                //    << endl;
1096 //            }
1097 //        }
1098 //        // Restore points.
1099 //        mesh_.movePoints(oldPoints);
1100 //    }
1101 //
1102 //
1103 //    Info<< "markFacesOnProblemCellsGeometric : marked "
1104 //        << returnReduce(nBaffleFaces, sumOp<label>())
1105 //        << " additional internal and coupled faces"
1106 //        << " to be converted into baffles." << endl;
1107 //
1108 //    syncTools::syncFaceList
1109 //    (
1110 //        mesh_,
1111 //        facePatch,
1112 //        maxEqOp<label>(),
1113 //        false               // no separation
1114 //    );
1115 //
1116 //    return facePatch;
1117 //}
1118
1119
1120 // ************************************************************************* //