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