FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

cellClassification.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2011 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 Description
00025 
00026 \*---------------------------------------------------------------------------*/
00027 
00028 #include "cellClassification.H"
00029 #include <meshTools/triSurfaceSearch.H>
00030 #include <meshTools/indexedOctree.H>
00031 #include <meshTools/treeDataFace.H>
00032 #include <meshTools/meshSearch.H>
00033 #include "cellInfo.H"
00034 #include <OpenFOAM/polyMesh.H>
00035 #include <OpenFOAM/MeshWave.H>
00036 #include <OpenFOAM/ListOps.H>
00037 #include <meshTools/meshTools.H>
00038 
00039 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00040 
00041 namespace Foam
00042 {
00043     defineTypeNameAndDebug(cellClassification, 0);
00044 }
00045 
00046 
00047 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00048 
00049 Foam::label Foam::cellClassification::count
00050 (
00051     const labelList& elems,
00052     const label elem
00053 )
00054 {
00055     label cnt = 0;
00056 
00057     forAll(elems, elemI)
00058     {
00059         if (elems[elemI] == elem)
00060         {
00061             cnt++;
00062         }
00063     }
00064     return cnt;
00065 
00066 }
00067 
00068 
00069 // Mark all faces that are cut by the surface. Two pass:
00070 // Pass1: mark all mesh edges that intersect surface (quick since triangle
00071 // pierce test).
00072 // Pass2: Check for all point neighbours of these faces whether any of their
00073 // faces are pierced.
00074 Foam::boolList Foam::cellClassification::markFaces
00075 (
00076     const triSurfaceSearch& search
00077 ) const
00078 {
00079     cpuTime timer;
00080 
00081     boolList cutFace(mesh_.nFaces(), false);
00082 
00083     label nCutFaces = 0;
00084 
00085     // Intersect mesh edges with surface (is fast) and mark all faces that
00086     // use edge.
00087 
00088     forAll(mesh_.edges(), edgeI)
00089     {
00090         if (debug && (edgeI % 10000 == 0))
00091         {
00092             Pout<< "Intersecting mesh edge " << edgeI << " with surface"
00093                 << endl;
00094         }
00095 
00096         const edge& e = mesh_.edges()[edgeI];
00097 
00098         const point& p0 = mesh_.points()[e.start()];
00099         const point& p1 = mesh_.points()[e.end()];
00100 
00101         pointIndexHit pHit(search.tree().findLineAny(p0, p1));
00102 
00103         if (pHit.hit())
00104         {
00105             const labelList& myFaces = mesh_.edgeFaces()[edgeI];
00106 
00107             forAll(myFaces, myFaceI)
00108             {
00109                 label faceI = myFaces[myFaceI];
00110 
00111                 if (!cutFace[faceI])
00112                 {
00113                     cutFace[faceI] = true;
00114 
00115                     nCutFaces++;
00116                 }
00117             }
00118         }
00119     }
00120 
00121     if (debug)
00122     {
00123         Pout<< "Intersected edges of mesh with surface in = "
00124             << timer.cpuTimeIncrement() << " s\n" << endl << endl;
00125     }
00126 
00127     //
00128     // Construct octree on faces that have not yet been marked as cut
00129     //
00130 
00131     labelList allFaces(mesh_.nFaces() - nCutFaces);
00132 
00133     label allFaceI = 0;
00134 
00135     forAll(cutFace, faceI)
00136     {
00137         if (!cutFace[faceI])
00138         {
00139             allFaces[allFaceI++] = faceI;
00140         }
00141     }
00142 
00143     if (debug)
00144     {
00145         Pout<< "Testing " << allFaceI << " faces for piercing by surface"
00146             << endl;
00147     }
00148 
00149     treeBoundBox allBb(mesh_.points());
00150 
00151     // Extend domain slightly (also makes it 3D if was 2D)
00152     scalar tol = 1E-6 * allBb.avgDim();
00153 
00154     point& bbMin = allBb.min();
00155     bbMin.x() -= tol;
00156     bbMin.y() -= tol;
00157     bbMin.z() -= tol;
00158 
00159     point& bbMax = allBb.max();
00160     bbMax.x() += 2*tol;
00161     bbMax.y() += 2*tol;
00162     bbMax.z() += 2*tol;
00163 
00164     indexedOctree<treeDataFace> faceTree
00165     (
00166         treeDataFace(false, mesh_, allFaces),
00167         allBb,      // overall search domain
00168         8,          // maxLevel
00169         10,         // leafsize
00170         3.0         // duplicity
00171     );
00172 
00173     const triSurface& surf = search.surface();
00174     const edgeList& edges = surf.edges();
00175     const pointField& localPoints = surf.localPoints();
00176 
00177     label nAddFaces = 0;
00178 
00179     forAll(edges, edgeI)
00180     {
00181         if (debug && (edgeI % 10000 == 0))
00182         {
00183             Pout<< "Intersecting surface edge " << edgeI
00184                 << " with mesh faces" << endl;
00185         }
00186         const edge& e = edges[edgeI];
00187 
00188         const point& start = localPoints[e.start()];
00189         const point& end = localPoints[e.end()];
00190 
00191         vector edgeNormal(end - start);
00192         const scalar edgeMag = mag(edgeNormal);
00193         const vector smallVec = 1E-9*edgeNormal;
00194 
00195         edgeNormal /= edgeMag+VSMALL;
00196 
00197         // Current start of pierce test
00198         point pt = start;
00199 
00200         while (true)
00201         {
00202             pointIndexHit pHit(faceTree.findLine(pt, end));
00203 
00204             if (!pHit.hit())
00205             {
00206                 break;
00207             }
00208             else
00209             {
00210                 label faceI = faceTree.shapes().faceLabels()[pHit.index()];
00211 
00212                 if (!cutFace[faceI])
00213                 {
00214                     cutFace[faceI] = true;
00215 
00216                     nAddFaces++;
00217                 }
00218 
00219                 // Restart from previous endpoint
00220                 pt = pHit.hitPoint() + smallVec;
00221 
00222                 if (((pt-start) & edgeNormal) >= edgeMag)
00223                 {
00224                     break;
00225                 }
00226             }
00227         }
00228     }
00229 
00230     if (debug)
00231     {
00232         Pout<< "Detected an additional " << nAddFaces << " faces cut"
00233             << endl;
00234 
00235         Pout<< "Intersected edges of surface with mesh faces in = "
00236             << timer.cpuTimeIncrement() << " s\n" << endl << endl;
00237     }
00238 
00239     return cutFace;
00240 }
00241 
00242 
00243 // Determine faces cut by surface and use to divide cells into types. See
00244 // cellInfo. All cells reachable from outsidePts are considered to be of type
00245 // 'outside'
00246 void Foam::cellClassification::markCells
00247 (
00248     const meshSearch& queryMesh,
00249     const boolList& piercedFace,
00250     const pointField& outsidePts
00251 )
00252 {
00253     // Use meshwave to partition mesh, starting from outsidePt
00254 
00255 
00256     // Construct null; sets type to NOTSET
00257     List<cellInfo> cellInfoList(mesh_.nCells());
00258 
00259     // Mark cut cells first
00260     forAll(piercedFace, faceI)
00261     {
00262         if (piercedFace[faceI])
00263         {
00264             cellInfoList[mesh_.faceOwner()[faceI]] =
00265                 cellInfo(cellClassification::CUT);
00266 
00267             if (mesh_.isInternalFace(faceI))
00268             {
00269                 cellInfoList[mesh_.faceNeighbour()[faceI]] =
00270                     cellInfo(cellClassification::CUT);
00271             }
00272         }
00273     }
00274 
00275     //
00276     // Mark cells containing outside points as being outside
00277     //
00278 
00279     // Coarse guess number of faces
00280     labelHashSet outsideFacesMap(outsidePts.size() * 6 * 2);
00281 
00282     forAll(outsidePts, outsidePtI)
00283     {
00284         // Use linear search for points.
00285         label cellI = queryMesh.findCell(outsidePts[outsidePtI], -1, false);
00286 
00287         if (returnReduce(cellI, maxOp<label>()) == -1)
00288         {
00289             FatalErrorIn
00290             (
00291                 "List<cellClassification::cType> markCells"
00292                 "(const meshSearch&, const boolList&, const pointField&)"
00293             )   << "outsidePoint " << outsidePts[outsidePtI]
00294                 << " is not inside any cell"
00295                 << nl << "It might be on a face or outside the geometry"
00296                 << exit(FatalError);
00297         }
00298 
00299         if (cellI >= 0)
00300         {
00301             cellInfoList[cellI] = cellInfo(cellClassification::OUTSIDE);
00302 
00303             // Mark faces of cellI
00304             const labelList& myFaces = mesh_.cells()[cellI];
00305             forAll(myFaces, myFaceI)
00306             {
00307                 outsideFacesMap.insert(myFaces[myFaceI]);
00308             }
00309         }
00310     }
00311 
00312 
00313     //
00314     // Mark faces to start wave from
00315     //
00316 
00317     labelList changedFaces(outsideFacesMap.toc());
00318 
00319     List<cellInfo> changedFacesInfo
00320     (
00321         changedFaces.size(),
00322         cellInfo(cellClassification::OUTSIDE)
00323     );
00324 
00325     MeshWave<cellInfo> cellInfoCalc
00326     (
00327         mesh_,
00328         changedFaces,                       // Labels of changed faces
00329         changedFacesInfo,                   // Information on changed faces
00330         cellInfoList,                       // Information on all cells
00331         mesh_.globalData().nTotalCells()    // max iterations
00332     );
00333 
00334     // Get information out of cellInfoList
00335     const List<cellInfo>& allInfo = cellInfoCalc.allCellInfo();
00336 
00337     forAll(allInfo, cellI)
00338     {
00339         label t = allInfo[cellI].type();
00340 
00341         if (t == cellClassification::NOTSET)
00342         {
00343             t = cellClassification::INSIDE;
00344         }
00345         operator[](cellI) = t;
00346     }
00347 }
00348 
00349 
00350 void Foam::cellClassification::classifyPoints
00351 (
00352     const label meshType,
00353     const labelList& cellType,
00354     List<pointStatus>& pointSide
00355 ) const
00356 {
00357     pointSide.setSize(mesh_.nPoints());
00358 
00359     forAll(mesh_.pointCells(), pointI)
00360     {
00361         const labelList& pCells = mesh_.pointCells()[pointI];
00362 
00363         pointSide[pointI] = UNSET;
00364 
00365         forAll(pCells, i)
00366         {
00367             label type = cellType[pCells[i]];
00368 
00369             if (type == meshType)
00370             {
00371                 if (pointSide[pointI] == UNSET)
00372                 {
00373                     pointSide[pointI] = MESH;
00374                 }
00375                 else if (pointSide[pointI] == NONMESH)
00376                 {
00377                     pointSide[pointI] = MIXED;
00378 
00379                     break;
00380                 }
00381             }
00382             else
00383             {
00384                 if (pointSide[pointI] == UNSET)
00385                 {
00386                     pointSide[pointI] = NONMESH;
00387                 }
00388                 else if (pointSide[pointI] == MESH)
00389                 {
00390                     pointSide[pointI] = MIXED;
00391 
00392                     break;
00393                 }
00394             }
00395         }
00396     }
00397 }
00398 
00399 
00400 bool Foam::cellClassification::usesMixedPointsOnly
00401 (
00402     const List<pointStatus>& pointSide,
00403     const label cellI
00404 ) const
00405 {
00406     const faceList& faces = mesh_.faces();
00407 
00408     const cell& cFaces = mesh_.cells()[cellI];
00409 
00410     forAll(cFaces, cFaceI)
00411     {
00412         const face& f = faces[cFaces[cFaceI]];
00413 
00414         forAll(f, fp)
00415         {
00416             if (pointSide[f[fp]] != MIXED)
00417             {
00418                 return false;
00419             }
00420         }
00421     }
00422 
00423     // All points are mixed.
00424     return true;
00425 }
00426 
00427 
00428 void Foam::cellClassification::getMeshOutside
00429 (
00430     const label meshType,
00431     faceList& outsideFaces,
00432     labelList& outsideOwner
00433 ) const
00434 {
00435     const labelList& own = mesh_.faceOwner();
00436     const labelList& nbr = mesh_.faceNeighbour();
00437     const faceList& faces = mesh_.faces();
00438 
00439     outsideFaces.setSize(mesh_.nFaces());
00440     outsideOwner.setSize(mesh_.nFaces());
00441     label outsideI = 0;
00442 
00443     // Get faces on interface between meshType and non-meshType
00444 
00445     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
00446     {
00447         label ownType = operator[](own[faceI]);
00448         label nbrType = operator[](nbr[faceI]);
00449 
00450         if (ownType == meshType && nbrType != meshType)
00451         {
00452             outsideFaces[outsideI] = faces[faceI];
00453             outsideOwner[outsideI] = own[faceI];    // meshType cell
00454             outsideI++;
00455         }
00456         else if (ownType != meshType && nbrType == meshType)
00457         {
00458             outsideFaces[outsideI] = faces[faceI];
00459             outsideOwner[outsideI] = nbr[faceI];    // meshType cell
00460             outsideI++;
00461         }
00462     }
00463 
00464     // Get faces on outside of real mesh with cells of meshType.
00465 
00466     for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
00467     {
00468         if (operator[](own[faceI]) == meshType)
00469         {
00470             outsideFaces[outsideI] = faces[faceI];
00471             outsideOwner[outsideI] = own[faceI];    // meshType cell
00472             outsideI++;
00473         }
00474     }
00475     outsideFaces.setSize(outsideI);
00476     outsideOwner.setSize(outsideI);
00477 }
00478 
00479 
00480 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00481 
00482 // Construct from mesh and surface and point(s) on outside
00483 Foam::cellClassification::cellClassification
00484 (
00485     const polyMesh& mesh,
00486     const meshSearch& meshQuery,
00487     const triSurfaceSearch& surfQuery,
00488     const pointField& outsidePoints
00489 )
00490 :
00491     labelList(mesh.nCells(), cellClassification::NOTSET),
00492     mesh_(mesh)
00493 {
00494     markCells
00495     (
00496         meshQuery,
00497         markFaces(surfQuery),
00498         outsidePoints
00499     );
00500 }
00501 
00502 
00503 // Construct from mesh and meshType.
00504 Foam::cellClassification::cellClassification
00505 (
00506     const polyMesh& mesh,
00507     const labelList& cellType
00508 )
00509 :
00510     labelList(cellType),
00511     mesh_(mesh)
00512 {
00513     if (mesh_.nCells() != size())
00514     {
00515         FatalErrorIn
00516         (
00517             "cellClassification::cellClassification"
00518             "(const polyMesh&, const labelList&)"
00519         )   << "Number of elements of cellType argument is not equal to the"
00520             << " number of cells" << abort(FatalError);
00521     }
00522 }
00523 
00524 
00525 // Construct as copy
00526 Foam::cellClassification::cellClassification(const cellClassification& cType)
00527 :
00528     labelList(cType),
00529     mesh_(cType.mesh())
00530 {}
00531 
00532 
00533 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00534 
00535 
00536 // Makes cutCells further than nLayers away from meshType to
00537 // fillType. Returns number of cells changed.
00538 Foam::label Foam::cellClassification::trimCutCells
00539 (
00540     const label nLayers,
00541     const label meshType,
00542     const label fillType
00543 )
00544 {
00545     // Temporary cell type for growing.
00546     labelList newCellType(*this);
00547 
00548 //    // Split types into outside and rest
00549 //    forAll(*this, cellI)
00550 //    {
00551 //        label type = operator[](cellI);
00552 //
00553 //        if (type == meshType)
00554 //        {
00555 //            newCellType[cellI] = type;
00556 //        }
00557 //        else
00558 //        {
00559 //            newCellType[cellI] = fillType;
00560 //        }
00561 //    }
00562 
00563     newCellType = *this;
00564 
00565 
00566     // Do point-cell-point walk on newCellType out from meshType cells
00567     for (label iter = 0; iter < nLayers; iter++)
00568     {
00569         // Get status of points: visible from meshType (pointSide == MESH)
00570         // or fillType/cutCells cells (pointSide == NONMESH) or
00571         // both (pointSide == MIXED)
00572         List<pointStatus> pointSide(mesh_.nPoints());
00573         classifyPoints(meshType, newCellType, pointSide);
00574 
00575         // Grow layer of outside cells
00576         forAll(pointSide, pointI)
00577         {
00578             if (pointSide[pointI] == MIXED)
00579             {
00580                 // Make cut
00581                 const labelList& pCells = mesh_.pointCells()[pointI];
00582 
00583                 forAll(pCells, i)
00584                 {
00585                     label type = newCellType[pCells[i]];
00586 
00587                     if (type == cellClassification::CUT)
00588                     {
00589                         // Found cut cell sharing point with
00590                         // mesh type cell. Grow.
00591                         newCellType[pCells[i]] = meshType;
00592                     }
00593                 }
00594             }
00595         }
00596     }
00597 
00598     // Merge newCellType into *this:
00599     // - leave meshType cells intact
00600     // - leave nonMesh cells intact
00601     // - make cutcells fillType if they were not marked in newCellType
00602 
00603     label nChanged = 0;
00604 
00605     forAll(newCellType, cellI)
00606     {
00607         if (operator[](cellI) == cellClassification::CUT)
00608         {
00609             if (newCellType[cellI] != meshType)
00610             {
00611                 // Cell was cutCell but further than nLayers away from
00612                 // meshType. Convert to fillType.
00613                 operator[](cellI) = fillType;
00614                 nChanged++;
00615             }
00616         }
00617     }
00618 
00619     return nChanged;
00620 }
00621 
00622 
00623 // Grow surface by pointNeighbours of type meshType
00624 Foam::label Foam::cellClassification::growSurface
00625 (
00626     const label meshType,
00627     const label fillType
00628 )
00629 {
00630     boolList hasMeshType(mesh_.nPoints(), false);
00631 
00632     // Mark points used by meshType cells
00633 
00634     forAll(mesh_.pointCells(), pointI)
00635     {
00636         const labelList& myCells = mesh_.pointCells()[pointI];
00637 
00638         // Check if one of cells has meshType
00639         forAll(myCells, myCellI)
00640         {
00641             label type = operator[](myCells[myCellI]);
00642 
00643             if (type == meshType)
00644             {
00645                 hasMeshType[pointI] = true;
00646 
00647                 break;
00648             }
00649         }
00650     }
00651 
00652     // Change neighbours of marked points
00653 
00654     label nChanged = 0;
00655 
00656     forAll(hasMeshType, pointI)
00657     {
00658         if (hasMeshType[pointI])
00659         {
00660             const labelList& myCells = mesh_.pointCells()[pointI];
00661 
00662             forAll(myCells, myCellI)
00663             {
00664                 if (operator[](myCells[myCellI]) != meshType)
00665                 {
00666                     operator[](myCells[myCellI]) = fillType;
00667 
00668                     nChanged++;
00669                 }
00670             }
00671         }
00672     }
00673     return nChanged;
00674 }
00675 
00676 
00677 // Check all points used by cells of meshType for use of at least one point
00678 // which is not on the outside. If all points are on the outside of the mesh
00679 // this would probably result in a flattened cell when projecting the cell
00680 // onto the surface.
00681 Foam::label Foam::cellClassification::fillHangingCells
00682 (
00683     const label meshType,
00684     const label fillType,
00685     const label maxIter
00686 )
00687 {
00688     label nTotChanged = 0;
00689 
00690     for(label iter = 0; iter < maxIter; iter++)
00691     {
00692         label nChanged = 0;
00693 
00694         // Get status of points: visible from meshType or non-meshType cells
00695         List<pointStatus> pointSide(mesh_.nPoints());
00696         classifyPoints(meshType, *this, pointSide);
00697 
00698         // Check all cells using mixed point type for whether they use mixed
00699         // points only. Note: could probably speed this up by counting number
00700         // of mixed verts per face and mixed faces per cell or something?
00701         forAll(pointSide, pointI)
00702         {
00703             if (pointSide[pointI] == MIXED)
00704             {
00705                 const labelList& pCells = mesh_.pointCells()[pointI];
00706 
00707                 forAll(pCells, i)
00708                 {
00709                     label cellI = pCells[i];
00710 
00711                     if (operator[](cellI) == meshType)
00712                     {
00713                         if (usesMixedPointsOnly(pointSide, cellI))
00714                         {
00715                             operator[](cellI) = fillType;
00716 
00717                             nChanged++;
00718                         }
00719                     }
00720                 }
00721             }
00722         }
00723         nTotChanged += nChanged;
00724 
00725         Pout<< "removeHangingCells : changed " << nChanged
00726             << " hanging cells" << endl;
00727 
00728         if (nChanged == 0)
00729         {
00730             break;
00731         }
00732     }
00733 
00734     return nTotChanged;
00735 }
00736 
00737 
00738 Foam::label Foam::cellClassification::fillRegionEdges
00739 (
00740     const label meshType,
00741     const label fillType,
00742     const label maxIter
00743 )
00744 {
00745     label nTotChanged = 0;
00746 
00747     for(label iter = 0; iter < maxIter; iter++)
00748     {
00749         // Get interface between meshType cells and non-meshType cells as a list
00750         // of faces and for each face the cell which is the meshType.
00751         faceList outsideFaces;
00752         labelList outsideOwner;
00753 
00754         getMeshOutside(meshType, outsideFaces, outsideOwner);
00755 
00756         // Build primitivePatch out of it and check it for problems.
00757         primitiveFacePatch fp(outsideFaces, mesh_.points());
00758 
00759         const labelListList& edgeFaces = fp.edgeFaces();
00760 
00761         label nChanged = 0;
00762 
00763         // Check all edgeFaces for non-manifoldness
00764 
00765         forAll(edgeFaces, edgeI)
00766         {
00767             const labelList& eFaces = edgeFaces[edgeI];
00768 
00769             if (eFaces.size() > 2)
00770             {
00771                 // patch connected through pinched edge. Remove first face using
00772                 // edge (and not yet changed)
00773 
00774                 forAll(eFaces, i)
00775                 {
00776                     label patchFaceI = eFaces[i];
00777 
00778                     label ownerCell = outsideOwner[patchFaceI];
00779 
00780                     if (operator[](ownerCell) == meshType)
00781                     {
00782                         operator[](ownerCell) = fillType;
00783 
00784                         nChanged++;
00785 
00786                         break;
00787                     }
00788                 }
00789             }
00790         }
00791 
00792         nTotChanged += nChanged;
00793 
00794         Pout<< "fillRegionEdges : changed " << nChanged
00795             << " cells using multiply connected edges" << endl;
00796 
00797         if (nChanged == 0)
00798         {
00799             break;
00800         }
00801     }
00802 
00803     return nTotChanged;
00804 }
00805 
00806 
00807 Foam::label Foam::cellClassification::fillRegionPoints
00808 (
00809     const label meshType,
00810     const label fillType,
00811     const label maxIter
00812 )
00813 {
00814     label nTotChanged = 0;
00815 
00816     for(label iter = 0; iter < maxIter; iter++)
00817     {
00818         // Get interface between meshType cells and non-meshType cells as a list
00819         // of faces and for each face the cell which is the meshType.
00820         faceList outsideFaces;
00821         labelList outsideOwner;
00822 
00823         getMeshOutside(meshType, outsideFaces, outsideOwner);
00824 
00825         // Build primitivePatch out of it and check it for problems.
00826         primitiveFacePatch fp(outsideFaces, mesh_.points());
00827 
00828         labelHashSet nonManifoldPoints;
00829 
00830         // Check for non-manifold points.
00831         fp.checkPointManifold(false, &nonManifoldPoints);
00832 
00833         const Map<label>& meshPointMap = fp.meshPointMap();
00834 
00835         label nChanged = 0;
00836 
00837         for
00838         (
00839             labelHashSet::const_iterator iter = nonManifoldPoints.begin();
00840             iter != nonManifoldPoints.end();
00841             ++iter
00842         )
00843         {
00844             // Find a face on fp using point and remove it.
00845             label patchPointI = meshPointMap[iter.key()];
00846 
00847             const labelList& pFaces = fp.pointFaces()[patchPointI];
00848 
00849             // Remove any face using conflicting point. Does first face which
00850             // has not yet been done. Could be more intelligent and decide which
00851             // one would be best to remove.
00852             forAll(pFaces, i)
00853             {
00854                 label patchFaceI = pFaces[i];
00855 
00856                 label ownerCell = outsideOwner[patchFaceI];
00857 
00858                 if (operator[](ownerCell) == meshType)
00859                 {
00860                     operator[](ownerCell) = fillType;
00861 
00862                     nChanged++;
00863 
00864                     break;
00865                 }
00866             }
00867         }
00868 
00869         nTotChanged += nChanged;
00870 
00871         Pout<< "fillRegionPoints : changed " << nChanged
00872             << " cells using multiply connected points" << endl;
00873 
00874         if (nChanged == 0)
00875         {
00876             break;
00877         }
00878     }
00879 
00880     return nTotChanged;
00881 }
00882 
00883 
00884 void Foam::cellClassification::writeStats(Ostream& os) const
00885 {
00886     os  << "Cells:" << size() << endl
00887         << "    notset  : " << count(*this, NOTSET) << endl
00888         << "    cut     : " << count(*this, CUT) << endl
00889         << "    inside  : " << count(*this, INSIDE) << endl
00890         << "    outside : " << count(*this, OUTSIDE) << endl;
00891 }
00892 
00893 
00894 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
00895 
00896 void Foam::cellClassification::operator=(const Foam::cellClassification& rhs)
00897 {
00898     labelList::operator=(rhs);
00899 }
00900 
00901 
00902 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines