src/meshTools/cellClassification/cellClassification.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 "cellClassification.H"
  27 #include "triSurfaceSearch.H"
  28 #include "indexedOctree.H"
  29 #include "treeDataFace.H"
  30 #include "meshSearch.H"
  31 #include "cellInfo.H"
  32 #include "polyMesh.H"
  33 #include "MeshWave.H"
  34 #include "ListOps.H"
  35 #include "meshTools.H"
  36 #include "cpuTime.H"
  37 #include "globalMeshData.H"
  38
  39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  40
  41 defineTypeNameAndDebug(Foam::cellClassification, 0);
  42
  43
  44 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  45
  46 Foam::label Foam::cellClassification::count
  47 (
  48     const labelList& elems,
  49     const label elem
  50 )
  51 {
  52     label cnt = 0;
  53
  54     forAll(elems, elemI)
  55     {
  56         if (elems[elemI] == elem)
  57         {
  58             cnt++;
  59         }
  60     }
  61     return cnt;
  62
  63 }
  64
  65
  66 // Mark all faces that are cut by the surface. Two pass:
  67 // Pass1: mark all mesh edges that intersect surface (quick since triangle
  68 // pierce test).
  69 // Pass2: Check for all point neighbours of these faces whether any of their
  70 // faces are pierced.
  71 Foam::boolList Foam::cellClassification::markFaces
  72 (
  73     const triSurfaceSearch& search
  74 ) const
  75 {
  76     cpuTime timer;
  77
  78     boolList cutFace(mesh_.nFaces(), false);
  79
  80     label nCutFaces = 0;
  81
  82     // Intersect mesh edges with surface (is fast) and mark all faces that
  83     // use edge.
  84
  85     forAll(mesh_.edges(), edgeI)
  86     {
  87         if (debug && (edgeI % 10000 == 0))
  88         {
  89             Pout<< "Intersecting mesh edge " << edgeI << " with surface"
  90                 << endl;
  91         }
  92
  93         const edge& e = mesh_.edges()[edgeI];
  94
  95         const point& p0 = mesh_.points()[e.start()];
  96         const point& p1 = mesh_.points()[e.end()];
  97
  98         pointIndexHit pHit(search.tree().findLineAny(p0, p1));
  99
 100         if (pHit.hit())
 101         {
 102             const labelList& myFaces = mesh_.edgeFaces()[edgeI];
 103
 104             forAll(myFaces, myFaceI)
 105             {
 106                 label faceI = myFaces[myFaceI];
 107
 108                 if (!cutFace[faceI])
 109                 {
 110                     cutFace[faceI] = true;
 111
 112                     nCutFaces++;
 113                 }
 114             }
 115         }
 116     }
 117
 118     if (debug)
 119     {
 120         Pout<< "Intersected edges of mesh with surface in = "
 121             << timer.cpuTimeIncrement() << " s\n" << endl << endl;
 122     }
 123
 124     //
 125     // Construct octree on faces that have not yet been marked as cut
 126     //
 127
 128     labelList allFaces(mesh_.nFaces() - nCutFaces);
 129
 130     label allFaceI = 0;
 131
 132     forAll(cutFace, faceI)
 133     {
 134         if (!cutFace[faceI])
 135         {
 136             allFaces[allFaceI++] = faceI;
 137         }
 138     }
 139
 140     if (debug)
 141     {
 142         Pout<< "Testing " << allFaceI << " faces for piercing by surface"
 143             << endl;
 144     }
 145
 146     treeBoundBox allBb(mesh_.points());
 147
 148     // Extend domain slightly (also makes it 3D if was 2D)
 149     scalar tol = 1E-6 * allBb.avgDim();
 150
 151     point& bbMin = allBb.min();
 152     bbMin.x() -= tol;
 153     bbMin.y() -= tol;
 154     bbMin.z() -= tol;
 155
 156     point& bbMax = allBb.max();
 157     bbMax.x() += 2*tol;
 158     bbMax.y() += 2*tol;
 159     bbMax.z() += 2*tol;
 160
 161     indexedOctree<treeDataFace> faceTree
 162     (
 163         treeDataFace(false, mesh_, allFaces),
 164         allBb,      // overall search domain
 165         8,          // maxLevel
 166         10,         // leafsize
 167         3.0         // duplicity
 168     );
 169
 170     const triSurface& surf = search.surface();
 171     const edgeList& edges = surf.edges();
 172     const pointField& localPoints = surf.localPoints();
 173
 174     label nAddFaces = 0;
 175
 176     forAll(edges, edgeI)
 177     {
 178         if (debug && (edgeI % 10000 == 0))
 179         {
 180             Pout<< "Intersecting surface edge " << edgeI
 181                 << " with mesh faces" << endl;
 182         }
 183         const edge& e = edges[edgeI];
 184
 185         const point& start = localPoints[e.start()];
 186         const point& end = localPoints[e.end()];
 187
 188         vector edgeNormal(end - start);
 189         const scalar edgeMag = mag(edgeNormal);
 190         const vector smallVec = 1E-9*edgeNormal;
 191
 192         edgeNormal /= edgeMag+VSMALL;
 193
 194         // Current start of pierce test
 195         point pt = start;
 196
 197         while (true)
 198         {
 199             pointIndexHit pHit(faceTree.findLine(pt, end));
 200
 201             if (!pHit.hit())
 202             {
 203                 break;
 204             }
 205             else
 206             {
 207                 label faceI = faceTree.shapes().faceLabels()[pHit.index()];
 208
 209                 if (!cutFace[faceI])
 210                 {
 211                     cutFace[faceI] = true;
 212
 213                     nAddFaces++;
 214                 }
 215
 216                 // Restart from previous endpoint
 217                 pt = pHit.hitPoint() + smallVec;
 218
 219                 if (((pt-start) & edgeNormal) >= edgeMag)
 220                 {
 221                     break;
 222                 }
 223             }
 224         }
 225     }
 226
 227     if (debug)
 228     {
 229         Pout<< "Detected an additional " << nAddFaces << " faces cut"
 230             << endl;
 231
 232         Pout<< "Intersected edges of surface with mesh faces in = "
 233             << timer.cpuTimeIncrement() << " s\n" << endl << endl;
 234     }
 235
 236     return cutFace;
 237 }
 238
 239
 240 // Determine faces cut by surface and use to divide cells into types. See
 241 // cellInfo. All cells reachable from outsidePts are considered to be of type
 242 // 'outside'
 243 void Foam::cellClassification::markCells
 244 (
 245     const meshSearch& queryMesh,
 246     const boolList& piercedFace,
 247     const pointField& outsidePts
 248 )
 249 {
 250     // Use meshwave to partition mesh, starting from outsidePt
 251
 252
 253     // Construct null; sets type to NOTSET
 254     List<cellInfo> cellInfoList(mesh_.nCells());
 255
 256     // Mark cut cells first
 257     forAll(piercedFace, faceI)
 258     {
 259         if (piercedFace[faceI])
 260         {
 261             cellInfoList[mesh_.faceOwner()[faceI]] =
 262                 cellInfo(cellClassification::CUT);
 263
 264             if (mesh_.isInternalFace(faceI))
 265             {
 266                 cellInfoList[mesh_.faceNeighbour()[faceI]] =
 267                     cellInfo(cellClassification::CUT);
 268             }
 269         }
 270     }
 271
 272     //
 273     // Mark cells containing outside points as being outside
 274     //
 275
 276     // Coarse guess number of faces
 277     labelHashSet outsideFacesMap(outsidePts.size() * 6 * 2);
 278
 279     forAll(outsidePts, outsidePtI)
 280     {
 281         // Use linear search for points.
 282         label cellI = queryMesh.findCell(outsidePts[outsidePtI], -1, false);
 283
 284         if (returnReduce(cellI, maxOp<label>()) == -1)
 285         {
 286             FatalErrorIn
 287             (
 288                 "List<cellClassification::cType> markCells"
 289                 "(const meshSearch&, const boolList&, const pointField&)"
 290             )   << "outsidePoint " << outsidePts[outsidePtI]
 291                 << " is not inside any cell"
 292                 << nl << "It might be on a face or outside the geometry"
 293                 << exit(FatalError);
 294         }
 295
 296         if (cellI >= 0)
 297         {
 298             cellInfoList[cellI] = cellInfo(cellClassification::OUTSIDE);
 299
 300             // Mark faces of cellI
 301             const labelList& myFaces = mesh_.cells()[cellI];
 302             forAll(myFaces, myFaceI)
 303             {
 304                 outsideFacesMap.insert(myFaces[myFaceI]);
 305             }
 306         }
 307     }
 308
 309
 310     //
 311     // Mark faces to start wave from
 312     //
 313
 314     labelList changedFaces(outsideFacesMap.toc());
 315
 316     List<cellInfo> changedFacesInfo
 317     (
 318         changedFaces.size(),
 319         cellInfo(cellClassification::OUTSIDE)
 320     );
 321
 322     MeshWave<cellInfo> cellInfoCalc
 323     (
 324         mesh_,
 325         changedFaces,                       // Labels of changed faces
 326         changedFacesInfo,                   // Information on changed faces
 327         cellInfoList,                       // Information on all cells
 328         mesh_.globalData().nTotalCells()+1  // max iterations
 329     );
 330
 331     // Get information out of cellInfoList
 332     const List<cellInfo>& allInfo = cellInfoCalc.allCellInfo();
 333
 334     forAll(allInfo, cellI)
 335     {
 336         label t = allInfo[cellI].type();
 337
 338         if (t == cellClassification::NOTSET)
 339         {
 340             t = cellClassification::INSIDE;
 341         }
 342         operator[](cellI) = t;
 343     }
 344 }
 345
 346
 347 void Foam::cellClassification::classifyPoints
 348 (
 349     const label meshType,
 350     const labelList& cellType,
 351     List<pointStatus>& pointSide
 352 ) const
 353 {
 354     pointSide.setSize(mesh_.nPoints());
 355
 356     forAll(mesh_.pointCells(), pointI)
 357     {
 358         const labelList& pCells = mesh_.pointCells()[pointI];
 359
 360         pointSide[pointI] = UNSET;
 361
 362         forAll(pCells, i)
 363         {
 364             label type = cellType[pCells[i]];
 365
 366             if (type == meshType)
 367             {
 368                 if (pointSide[pointI] == UNSET)
 369                 {
 370                     pointSide[pointI] = MESH;
 371                 }
 372                 else if (pointSide[pointI] == NONMESH)
 373                 {
 374                     pointSide[pointI] = MIXED;
 375
 376                     break;
 377                 }
 378             }
 379             else
 380             {
 381                 if (pointSide[pointI] == UNSET)
 382                 {
 383                     pointSide[pointI] = NONMESH;
 384                 }
 385                 else if (pointSide[pointI] == MESH)
 386                 {
 387                     pointSide[pointI] = MIXED;
 388
 389                     break;
 390                 }
 391             }
 392         }
 393     }
 394 }
 395
 396
 397 bool Foam::cellClassification::usesMixedPointsOnly
 398 (
 399     const List<pointStatus>& pointSide,
 400     const label cellI
 401 ) const
 402 {
 403     const faceList& faces = mesh_.faces();
 404
 405     const cell& cFaces = mesh_.cells()[cellI];
 406
 407     forAll(cFaces, cFaceI)
 408     {
 409         const face& f = faces[cFaces[cFaceI]];
 410
 411         forAll(f, fp)
 412         {
 413             if (pointSide[f[fp]] != MIXED)
 414             {
 415                 return false;
 416             }
 417         }
 418     }
 419
 420     // All points are mixed.
 421     return true;
 422 }
 423
 424
 425 void Foam::cellClassification::getMeshOutside
 426 (
 427     const label meshType,
 428     faceList& outsideFaces,
 429     labelList& outsideOwner
 430 ) const
 431 {
 432     const labelList& own = mesh_.faceOwner();
 433     const labelList& nbr = mesh_.faceNeighbour();
 434     const faceList& faces = mesh_.faces();
 435
 436     outsideFaces.setSize(mesh_.nFaces());
 437     outsideOwner.setSize(mesh_.nFaces());
 438     label outsideI = 0;
 439
 440     // Get faces on interface between meshType and non-meshType
 441
 442     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
 443     {
 444         label ownType = operator[](own[faceI]);
 445         label nbrType = operator[](nbr[faceI]);
 446
 447         if (ownType == meshType && nbrType != meshType)
 448         {
 449             outsideFaces[outsideI] = faces[faceI];
 450             outsideOwner[outsideI] = own[faceI];    // meshType cell
 451             outsideI++;
 452         }
 453         else if (ownType != meshType && nbrType == meshType)
 454         {
 455             outsideFaces[outsideI] = faces[faceI];
 456             outsideOwner[outsideI] = nbr[faceI];    // meshType cell
 457             outsideI++;
 458         }
 459     }
 460
 461     // Get faces on outside of real mesh with cells of meshType.
 462
 463     for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
 464     {
 465         if (operator[](own[faceI]) == meshType)
 466         {
 467             outsideFaces[outsideI] = faces[faceI];
 468             outsideOwner[outsideI] = own[faceI];    // meshType cell
 469             outsideI++;
 470         }
 471     }
 472     outsideFaces.setSize(outsideI);
 473     outsideOwner.setSize(outsideI);
 474 }
 475
 476
 477 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 478
 479 // Construct from mesh and surface and point(s) on outside
 480 Foam::cellClassification::cellClassification
 481 (
 482     const polyMesh& mesh,
 483     const meshSearch& meshQuery,
 484     const triSurfaceSearch& surfQuery,
 485     const pointField& outsidePoints
 486 )
 487 :
 488     labelList(mesh.nCells(), cellClassification::NOTSET),
 489     mesh_(mesh)
 490 {
 491     markCells
 492     (
 493         meshQuery,
 494         markFaces(surfQuery),
 495         outsidePoints
 496     );
 497 }
 498
 499
 500 // Construct from mesh and meshType.
 501 Foam::cellClassification::cellClassification
 502 (
 503     const polyMesh& mesh,
 504     const labelList& cellType
 505 )
 506 :
 507     labelList(cellType),
 508     mesh_(mesh)
 509 {
 510     if (mesh_.nCells() != size())
 511     {
 512         FatalErrorIn
 513         (
 514             "cellClassification::cellClassification"
 515             "(const polyMesh&, const labelList&)"
 516         )   << "Number of elements of cellType argument is not equal to the"
 517             << " number of cells" << abort(FatalError);
 518     }
 519 }
 520
 521
 522 // Construct as copy
 523 Foam::cellClassification::cellClassification(const cellClassification& cType)
 524 :
 525     labelList(cType),
 526     mesh_(cType.mesh())
 527 {}
 528
 529
 530 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 531
 532
 533 // Makes cutCells further than nLayers away from meshType to
 534 // fillType. Returns number of cells changed.
 535 Foam::label Foam::cellClassification::trimCutCells
 536 (
 537     const label nLayers,
 538     const label meshType,
 539     const label fillType
 540 )
 541 {
 542     // Temporary cell type for growing.
 543     labelList newCellType(*this);
 544
 545 //    // Split types into outside and rest
 546 //    forAll(*this, cellI)
 547 //    {
 548 //        label type = operator[](cellI);
 549 //
 550 //        if (type == meshType)
 551 //        {
 552 //            newCellType[cellI] = type;
 553 //        }
 554 //        else
 555 //        {
 556 //            newCellType[cellI] = fillType;
 557 //        }
 558 //    }
 559
 560     newCellType = *this;
 561
 562
 563     // Do point-cell-point walk on newCellType out from meshType cells
 564     for (label iter = 0; iter < nLayers; iter++)
 565     {
 566         // Get status of points: visible from meshType (pointSide == MESH)
 567         // or fillType/cutCells cells (pointSide == NONMESH) or
 568         // both (pointSide == MIXED)
 569         List<pointStatus> pointSide(mesh_.nPoints());
 570         classifyPoints(meshType, newCellType, pointSide);
 571
 572         // Grow layer of outside cells
 573         forAll(pointSide, pointI)
 574         {
 575             if (pointSide[pointI] == MIXED)
 576             {
 577                 // Make cut
 578                 const labelList& pCells = mesh_.pointCells()[pointI];
 579
 580                 forAll(pCells, i)
 581                 {
 582                     label type = newCellType[pCells[i]];
 583
 584                     if (type == cellClassification::CUT)
 585                     {
 586                         // Found cut cell sharing point with
 587                         // mesh type cell. Grow.
 588                         newCellType[pCells[i]] = meshType;
 589                     }
 590                 }
 591             }
 592         }
 593     }
 594
 595     // Merge newCellType into *this:
 596     // - leave meshType cells intact
 597     // - leave nonMesh cells intact
 598     // - make cutcells fillType if they were not marked in newCellType
 599
 600     label nChanged = 0;
 601
 602     forAll(newCellType, cellI)
 603     {
 604         if (operator[](cellI) == cellClassification::CUT)
 605         {
 606             if (newCellType[cellI] != meshType)
 607             {
 608                 // Cell was cutCell but further than nLayers away from
 609                 // meshType. Convert to fillType.
 610                 operator[](cellI) = fillType;
 611                 nChanged++;
 612             }
 613         }
 614     }
 615
 616     return nChanged;
 617 }
 618
 619
 620 // Grow surface by pointNeighbours of type meshType
 621 Foam::label Foam::cellClassification::growSurface
 622 (
 623     const label meshType,
 624     const label fillType
 625 )
 626 {
 627     boolList hasMeshType(mesh_.nPoints(), false);
 628
 629     // Mark points used by meshType cells
 630
 631     forAll(mesh_.pointCells(), pointI)
 632     {
 633         const labelList& myCells = mesh_.pointCells()[pointI];
 634
 635         // Check if one of cells has meshType
 636         forAll(myCells, myCellI)
 637         {
 638             label type = operator[](myCells[myCellI]);
 639
 640             if (type == meshType)
 641             {
 642                 hasMeshType[pointI] = true;
 643
 644                 break;
 645             }
 646         }
 647     }
 648
 649     // Change neighbours of marked points
 650
 651     label nChanged = 0;
 652
 653     forAll(hasMeshType, pointI)
 654     {
 655         if (hasMeshType[pointI])
 656         {
 657             const labelList& myCells = mesh_.pointCells()[pointI];
 658
 659             forAll(myCells, myCellI)
 660             {
 661                 if (operator[](myCells[myCellI]) != meshType)
 662                 {
 663                     operator[](myCells[myCellI]) = fillType;
 664
 665                     nChanged++;
 666                 }
 667             }
 668         }
 669     }
 670     return nChanged;
 671 }
 672
 673
 674 // Check all points used by cells of meshType for use of at least one point
 675 // which is not on the outside. If all points are on the outside of the mesh
 676 // this would probably result in a flattened cell when projecting the cell
 677 // onto the surface.
 678 Foam::label Foam::cellClassification::fillHangingCells
 679 (
 680     const label meshType,
 681     const label fillType,
 682     const label maxIter
 683 )
 684 {
 685     label nTotChanged = 0;
 686
 687     for (label iter = 0; iter < maxIter; iter++)
 688     {
 689         label nChanged = 0;
 690
 691         // Get status of points: visible from meshType or non-meshType cells
 692         List<pointStatus> pointSide(mesh_.nPoints());
 693         classifyPoints(meshType, *this, pointSide);
 694
 695         // Check all cells using mixed point type for whether they use mixed
 696         // points only. Note: could probably speed this up by counting number
 697         // of mixed verts per face and mixed faces per cell or something?
 698         forAll(pointSide, pointI)
 699         {
 700             if (pointSide[pointI] == MIXED)
 701             {
 702                 const labelList& pCells = mesh_.pointCells()[pointI];
 703
 704                 forAll(pCells, i)
 705                 {
 706                     label cellI = pCells[i];
 707
 708                     if (operator[](cellI) == meshType)
 709                     {
 710                         if (usesMixedPointsOnly(pointSide, cellI))
 711                         {
 712                             operator[](cellI) = fillType;
 713
 714                             nChanged++;
 715                         }
 716                     }
 717                 }
 718             }
 719         }
 720         nTotChanged += nChanged;
 721
 722         Pout<< "removeHangingCells : changed " << nChanged
 723             << " hanging cells" << endl;
 724
 725         if (nChanged == 0)
 726         {
 727             break;
 728         }
 729     }
 730
 731     return nTotChanged;
 732 }
 733
 734
 735 Foam::label Foam::cellClassification::fillRegionEdges
 736 (
 737     const label meshType,
 738     const label fillType,
 739     const label maxIter
 740 )
 741 {
 742     label nTotChanged = 0;
 743
 744     for (label iter = 0; iter < maxIter; iter++)
 745     {
 746         // Get interface between meshType cells and non-meshType cells as a list
 747         // of faces and for each face the cell which is the meshType.
 748         faceList outsideFaces;
 749         labelList outsideOwner;
 750
 751         getMeshOutside(meshType, outsideFaces, outsideOwner);
 752
 753         // Build primitivePatch out of it and check it for problems.
 754         primitiveFacePatch fp(outsideFaces, mesh_.points());
 755
 756         const labelListList& edgeFaces = fp.edgeFaces();
 757
 758         label nChanged = 0;
 759
 760         // Check all edgeFaces for non-manifoldness
 761
 762         forAll(edgeFaces, edgeI)
 763         {
 764             const labelList& eFaces = edgeFaces[edgeI];
 765
 766             if (eFaces.size() > 2)
 767             {
 768                 // patch connected through pinched edge. Remove first face using
 769                 // edge (and not yet changed)
 770
 771                 forAll(eFaces, i)
 772                 {
 773                     label patchFaceI = eFaces[i];
 774
 775                     label ownerCell = outsideOwner[patchFaceI];
 776
 777                     if (operator[](ownerCell) == meshType)
 778                     {
 779                         operator[](ownerCell) = fillType;
 780
 781                         nChanged++;
 782
 783                         break;
 784                     }
 785                 }
 786             }
 787         }
 788
 789         nTotChanged += nChanged;
 790
 791         Pout<< "fillRegionEdges : changed " << nChanged
 792             << " cells using multiply connected edges" << endl;
 793
 794         if (nChanged == 0)
 795         {
 796             break;
 797         }
 798     }
 799
 800     return nTotChanged;
 801 }
 802
 803
 804 Foam::label Foam::cellClassification::fillRegionPoints
 805 (
 806     const label meshType,
 807     const label fillType,
 808     const label maxIter
 809 )
 810 {
 811     label nTotChanged = 0;
 812
 813     for (label iter = 0; iter < maxIter; iter++)
 814     {
 815         // Get interface between meshType cells and non-meshType cells as a list
 816         // of faces and for each face the cell which is the meshType.
 817         faceList outsideFaces;
 818         labelList outsideOwner;
 819
 820         getMeshOutside(meshType, outsideFaces, outsideOwner);
 821
 822         // Build primitivePatch out of it and check it for problems.
 823         primitiveFacePatch fp(outsideFaces, mesh_.points());
 824
 825         labelHashSet nonManifoldPoints;
 826
 827         // Check for non-manifold points.
 828         fp.checkPointManifold(false, &nonManifoldPoints);
 829
 830         const Map<label>& meshPointMap = fp.meshPointMap();
 831
 832         label nChanged = 0;
 833
 834         forAllConstIter(labelHashSet, nonManifoldPoints, iter)
 835         {
 836             // Find a face on fp using point and remove it.
 837             const label patchPointI = meshPointMap[iter.key()];
 838
 839             const labelList& pFaces = fp.pointFaces()[patchPointI];
 840
 841             // Remove any face using conflicting point. Does first face which
 842             // has not yet been done. Could be more intelligent and decide which
 843             // one would be best to remove.
 844             forAll(pFaces, i)
 845             {
 846                 const label patchFaceI = pFaces[i];
 847                 const label ownerCell  = outsideOwner[patchFaceI];
 848
 849                 if (operator[](ownerCell) == meshType)
 850                 {
 851                     operator[](ownerCell) = fillType;
 852
 853                     nChanged++;
 854                     break;
 855                 }
 856             }
 857         }
 858
 859         nTotChanged += nChanged;
 860
 861         Pout<< "fillRegionPoints : changed " << nChanged
 862             << " cells using multiply connected points" << endl;
 863
 864         if (nChanged == 0)
 865         {
 866             break;
 867         }
 868     }
 869
 870     return nTotChanged;
 871 }
 872
 873
 874 void Foam::cellClassification::writeStats(Ostream& os) const
 875 {
 876     os  << "Cells:" << size() << endl
 877         << "    notset  : " << count(*this, NOTSET) << endl
 878         << "    cut     : " << count(*this, CUT) << endl
 879         << "    inside  : " << count(*this, INSIDE) << endl
 880         << "    outside : " << count(*this, OUTSIDE) << endl;
 881 }
 882
 883
 884 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 885
 886 void Foam::cellClassification::operator=(const Foam::cellClassification& rhs)
 887 {
 888     labelList::operator=(rhs);
 889 }
 890
 891
 892 // ************************************************************************* //