1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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.
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
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
27 \*---------------------------------------------------------------------------*/
29 #include "cellClassification.H"
30 #include <meshTools/triSurfaceSearch.H>
31 #include <meshTools/indexedOctree.H>
32 #include <meshTools/treeDataFace.H>
33 #include <meshTools/meshSearch.H>
35 #include <OpenFOAM/polyMesh.H>
36 #include <OpenFOAM/MeshWave.H>
37 #include <OpenFOAM/ListOps.H>
38 #include <meshTools/meshTools.H>
40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
44 defineTypeNameAndDebug(cellClassification, 0);
48 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
50 Foam::label Foam::cellClassification::count
52 const labelList& elems,
60 if (elems[elemI] == elem)
70 // Mark all faces that are cut by the surface. Two pass:
71 // Pass1: mark all mesh edges that intersect surface (quick since triangle
73 // Pass2: Check for all point neighbours of these faces whether any of their
75 Foam::boolList Foam::cellClassification::markFaces
77 const triSurfaceSearch& search
82 boolList cutFace(mesh_.nFaces(), false);
86 // Intersect mesh edges with surface (is fast) and mark all faces that
89 forAll(mesh_.edges(), edgeI)
91 if (debug && (edgeI % 10000 == 0))
93 Pout<< "Intersecting mesh edge " << edgeI << " with surface"
97 const edge& e = mesh_.edges()[edgeI];
99 const point& p0 = mesh_.points()[e.start()];
100 const point& p1 = mesh_.points()[e.end()];
102 pointIndexHit pHit(search.tree().findLineAny(p0, p1));
106 const labelList& myFaces = mesh_.edgeFaces()[edgeI];
108 forAll(myFaces, myFaceI)
110 label faceI = myFaces[myFaceI];
114 cutFace[faceI] = true;
124 Pout<< "Intersected edges of mesh with surface in = "
125 << timer.cpuTimeIncrement() << " s\n" << endl << endl;
129 // Construct octree on faces that have not yet been marked as cut
132 labelList allFaces(mesh_.nFaces() - nCutFaces);
136 forAll(cutFace, faceI)
140 allFaces[allFaceI++] = faceI;
146 Pout<< "Testing " << allFaceI << " faces for piercing by surface"
150 treeBoundBox allBb(mesh_.points());
152 // Extend domain slightly (also makes it 3D if was 2D)
153 scalar tol = 1E-6 * allBb.avgDim();
155 point& bbMin = allBb.min();
160 point& bbMax = allBb.max();
165 indexedOctree<treeDataFace> faceTree
167 treeDataFace(false, mesh_, allFaces),
168 allBb, // overall search domain
174 const triSurface& surf = search.surface();
175 const edgeList& edges = surf.edges();
176 const pointField& localPoints = surf.localPoints();
182 if (debug && (edgeI % 10000 == 0))
184 Pout<< "Intersecting surface edge " << edgeI
185 << " with mesh faces" << endl;
187 const edge& e = edges[edgeI];
189 const point& start = localPoints[e.start()];
190 const point& end = localPoints[e.end()];
192 vector edgeNormal(end - start);
193 const scalar edgeMag = mag(edgeNormal);
194 const vector smallVec = 1E-9*edgeNormal;
196 edgeNormal /= edgeMag+VSMALL;
198 // Current start of pierce test
203 pointIndexHit pHit(faceTree.findLine(pt, end));
211 label faceI = faceTree.shapes().faceLabels()[pHit.index()];
215 cutFace[faceI] = true;
220 // Restart from previous endpoint
221 pt = pHit.hitPoint() + smallVec;
223 if (((pt-start) & edgeNormal) >= edgeMag)
233 Pout<< "Detected an additional " << nAddFaces << " faces cut"
236 Pout<< "Intersected edges of surface with mesh faces in = "
237 << timer.cpuTimeIncrement() << " s\n" << endl << endl;
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
247 void Foam::cellClassification::markCells
249 const meshSearch& queryMesh,
250 const boolList& piercedFace,
251 const pointField& outsidePts
254 // Use meshwave to partition mesh, starting from outsidePt
257 // Construct null; sets type to NOTSET
258 List<cellInfo> cellInfoList(mesh_.nCells());
260 // Mark cut cells first
261 forAll(piercedFace, faceI)
263 if (piercedFace[faceI])
265 cellInfoList[mesh_.faceOwner()[faceI]] =
266 cellInfo(cellClassification::CUT);
268 if (mesh_.isInternalFace(faceI))
270 cellInfoList[mesh_.faceNeighbour()[faceI]] =
271 cellInfo(cellClassification::CUT);
277 // Mark cells containing outside points as being outside
280 // Coarse guess number of faces
281 labelHashSet outsideFacesMap(outsidePts.size() * 6 * 2);
283 forAll(outsidePts, outsidePtI)
285 // Use linear search for points.
286 label cellI = queryMesh.findCell(outsidePts[outsidePtI], -1, false);
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"
300 cellInfoList[cellI] = cellInfo(cellClassification::OUTSIDE);
302 // Mark faces of cellI
303 const labelList& myFaces = mesh_.cells()[cellI];
304 forAll(myFaces, myFaceI)
306 outsideFacesMap.insert(myFaces[myFaceI]);
312 // Mark faces to start wave from
315 labelList changedFaces(outsideFacesMap.toc());
317 List<cellInfo> changedFacesInfo
320 cellInfo(cellClassification::OUTSIDE)
323 MeshWave<cellInfo> cellInfoCalc
326 changedFaces, // Labels of changed faces
327 changedFacesInfo, // Information on changed faces
328 cellInfoList, // Information on all cells
329 mesh_.globalData().nTotalCells() // max iterations
332 // Get information out of cellInfoList
333 const List<cellInfo>& allInfo = cellInfoCalc.allCellInfo();
335 forAll(allInfo, cellI)
337 label t = allInfo[cellI].type();
339 if (t == cellClassification::NOTSET)
341 t = cellClassification::INSIDE;
343 operator[](cellI) = t;
348 void Foam::cellClassification::classifyPoints
350 const label meshType,
351 const labelList& cellType,
352 List<pointStatus>& pointSide
355 pointSide.setSize(mesh_.nPoints());
357 forAll(mesh_.pointCells(), pointI)
359 const labelList& pCells = mesh_.pointCells()[pointI];
361 pointSide[pointI] = UNSET;
365 label type = cellType[pCells[i]];
367 if (type == meshType)
369 if (pointSide[pointI] == UNSET)
371 pointSide[pointI] = MESH;
373 else if (pointSide[pointI] == NONMESH)
375 pointSide[pointI] = MIXED;
382 if (pointSide[pointI] == UNSET)
384 pointSide[pointI] = NONMESH;
386 else if (pointSide[pointI] == MESH)
388 pointSide[pointI] = MIXED;
398 bool Foam::cellClassification::usesMixedPointsOnly
400 const List<pointStatus>& pointSide,
404 const faceList& faces = mesh_.faces();
406 const cell& cFaces = mesh_.cells()[cellI];
408 forAll(cFaces, cFaceI)
410 const face& f = faces[cFaces[cFaceI]];
414 if (pointSide[f[fp]] != MIXED)
421 // All points are mixed.
426 void Foam::cellClassification::getMeshOutside
428 const label meshType,
429 faceList& outsideFaces,
430 labelList& outsideOwner
433 const labelList& own = mesh_.faceOwner();
434 const labelList& nbr = mesh_.faceNeighbour();
435 const faceList& faces = mesh_.faces();
437 outsideFaces.setSize(mesh_.nFaces());
438 outsideOwner.setSize(mesh_.nFaces());
441 // Get faces on interface between meshType and non-meshType
443 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
445 label ownType = operator[](own[faceI]);
446 label nbrType = operator[](nbr[faceI]);
448 if (ownType == meshType && nbrType != meshType)
450 outsideFaces[outsideI] = faces[faceI];
451 outsideOwner[outsideI] = own[faceI]; // meshType cell
454 else if (ownType != meshType && nbrType == meshType)
456 outsideFaces[outsideI] = faces[faceI];
457 outsideOwner[outsideI] = nbr[faceI]; // meshType cell
462 // Get faces on outside of real mesh with cells of meshType.
464 for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
466 if (operator[](own[faceI]) == meshType)
468 outsideFaces[outsideI] = faces[faceI];
469 outsideOwner[outsideI] = own[faceI]; // meshType cell
473 outsideFaces.setSize(outsideI);
474 outsideOwner.setSize(outsideI);
478 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
480 // Construct from mesh and surface and point(s) on outside
481 Foam::cellClassification::cellClassification
483 const polyMesh& mesh,
484 const meshSearch& meshQuery,
485 const triSurfaceSearch& surfQuery,
486 const pointField& outsidePoints
489 labelList(mesh.nCells(), cellClassification::NOTSET),
495 markFaces(surfQuery),
501 // Construct from mesh and meshType.
502 Foam::cellClassification::cellClassification
504 const polyMesh& mesh,
505 const labelList& cellType
511 if (mesh_.nCells() != size())
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);
524 Foam::cellClassification::cellClassification(const cellClassification& cType)
531 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
534 // Makes cutCells further than nLayers away from meshType to
535 // fillType. Returns number of cells changed.
536 Foam::label Foam::cellClassification::trimCutCells
539 const label meshType,
543 // Temporary cell type for growing.
544 labelList newCellType(*this);
546 // // Split types into outside and rest
547 // forAll(*this, cellI)
549 // label type = operator[](cellI);
551 // if (type == meshType)
553 // newCellType[cellI] = type;
557 // newCellType[cellI] = fillType;
564 // Do point-cell-point walk on newCellType out from meshType cells
565 for (label iter = 0; iter < nLayers; iter++)
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);
573 // Grow layer of outside cells
574 forAll(pointSide, pointI)
576 if (pointSide[pointI] == MIXED)
579 const labelList& pCells = mesh_.pointCells()[pointI];
583 label type = newCellType[pCells[i]];
585 if (type == cellClassification::CUT)
587 // Found cut cell sharing point with
588 // mesh type cell. Grow.
589 newCellType[pCells[i]] = meshType;
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
603 forAll(newCellType, cellI)
605 if (operator[](cellI) == cellClassification::CUT)
607 if (newCellType[cellI] != meshType)
609 // Cell was cutCell but further than nLayers away from
610 // meshType. Convert to fillType.
611 operator[](cellI) = fillType;
621 // Grow surface by pointNeighbours of type meshType
622 Foam::label Foam::cellClassification::growSurface
624 const label meshType,
628 boolList hasMeshType(mesh_.nPoints(), false);
630 // Mark points used by meshType cells
632 forAll(mesh_.pointCells(), pointI)
634 const labelList& myCells = mesh_.pointCells()[pointI];
636 // Check if one of cells has meshType
637 forAll(myCells, myCellI)
639 label type = operator[](myCells[myCellI]);
641 if (type == meshType)
643 hasMeshType[pointI] = true;
650 // Change neighbours of marked points
654 forAll(hasMeshType, pointI)
656 if (hasMeshType[pointI])
658 const labelList& myCells = mesh_.pointCells()[pointI];
660 forAll(myCells, myCellI)
662 if (operator[](myCells[myCellI]) != meshType)
664 operator[](myCells[myCellI]) = fillType;
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
679 Foam::label Foam::cellClassification::fillHangingCells
681 const label meshType,
682 const label fillType,
686 label nTotChanged = 0;
688 for(label iter = 0; iter < maxIter; iter++)
692 // Get status of points: visible from meshType or non-meshType cells
693 List<pointStatus> pointSide(mesh_.nPoints());
694 classifyPoints(meshType, *this, pointSide);
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)
701 if (pointSide[pointI] == MIXED)
703 const labelList& pCells = mesh_.pointCells()[pointI];
707 label cellI = pCells[i];
709 if (operator[](cellI) == meshType)
711 if (usesMixedPointsOnly(pointSide, cellI))
713 operator[](cellI) = fillType;
721 nTotChanged += nChanged;
723 Pout<< "removeHangingCells : changed " << nChanged
724 << " hanging cells" << endl;
736 Foam::label Foam::cellClassification::fillRegionEdges
738 const label meshType,
739 const label fillType,
743 label nTotChanged = 0;
745 for(label iter = 0; iter < maxIter; iter++)
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;
752 getMeshOutside(meshType, outsideFaces, outsideOwner);
754 // Build primitivePatch out of it and check it for problems.
755 primitiveFacePatch fp(outsideFaces, mesh_.points());
757 const labelListList& edgeFaces = fp.edgeFaces();
761 // Check all edgeFaces for non-manifoldness
763 forAll(edgeFaces, edgeI)
765 const labelList& eFaces = edgeFaces[edgeI];
767 if (eFaces.size() > 2)
769 // patch connected through pinched edge. Remove first face using
770 // edge (and not yet changed)
774 label patchFaceI = eFaces[i];
776 label ownerCell = outsideOwner[patchFaceI];
778 if (operator[](ownerCell) == meshType)
780 operator[](ownerCell) = fillType;
790 nTotChanged += nChanged;
792 Pout<< "fillRegionEdges : changed " << nChanged
793 << " cells using multiply connected edges" << endl;
805 Foam::label Foam::cellClassification::fillRegionPoints
807 const label meshType,
808 const label fillType,
812 label nTotChanged = 0;
814 for(label iter = 0; iter < maxIter; iter++)
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;
821 getMeshOutside(meshType, outsideFaces, outsideOwner);
823 // Build primitivePatch out of it and check it for problems.
824 primitiveFacePatch fp(outsideFaces, mesh_.points());
826 labelHashSet nonManifoldPoints;
828 // Check for non-manifold points.
829 fp.checkPointManifold(false, &nonManifoldPoints);
831 const Map<label>& meshPointMap = fp.meshPointMap();
837 labelHashSet::const_iterator iter = nonManifoldPoints.begin();
838 iter != nonManifoldPoints.end();
842 // Find a face on fp using point and remove it.
843 label patchPointI = meshPointMap[iter.key()];
845 const labelList& pFaces = fp.pointFaces()[patchPointI];
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.
852 label patchFaceI = pFaces[i];
854 label ownerCell = outsideOwner[patchFaceI];
856 if (operator[](ownerCell) == meshType)
858 operator[](ownerCell) = fillType;
867 nTotChanged += nChanged;
869 Pout<< "fillRegionPoints : changed " << nChanged
870 << " cells using multiply connected points" << endl;
882 void Foam::cellClassification::writeStats(Ostream& os) const
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;
892 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
894 void Foam::cellClassification::operator=(const Foam::cellClassification& rhs)
896 labelList::operator=(rhs);
900 // ************************ vim: set sw=4 sts=4 et: ************************ //