src/meshTools/cellClassification/cellClassification.C

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