49 Foam::label Foam::cellClassification::count
59 if (elems[elemI] == elem)
76 const triSurfaceSearch& search
81 boolList cutFace(mesh_.nFaces(),
false);
88 forAll(mesh_.edges(), edgeI)
90 if (debug && (edgeI % 10000 == 0))
92 Pout<<
"Intersecting mesh edge " << edgeI <<
" with surface"
96 const edge&
e = mesh_.edges()[edgeI];
98 const point& p0 = mesh_.points()[e.start()];
99 const point& p1 = mesh_.points()[e.end()];
105 const labelList& myFaces = mesh_.edgeFaces()[edgeI];
109 label faceI = myFaces[myFaceI];
113 cutFace[faceI] =
true;
123 Pout<<
"Intersected edges of mesh with surface in = "
124 << timer.cpuTimeIncrement() <<
" s\n" <<
endl <<
endl;
131 labelList allFaces(mesh_.nFaces() - nCutFaces);
139 allFaces[allFaceI++] = faceI;
145 Pout<<
"Testing " << allFaceI <<
" faces for piercing by surface"
149 treeBoundBox allBb(mesh_.points());
152 scalar tol = 1
E-6 * allBb.avgDim();
164 indexedOctree<treeDataFace> faceTree
166 treeDataFace(
false, mesh_, allFaces),
173 const triSurface& surf = search.surface();
174 const edgeList& edges = surf.edges();
175 const pointField& localPoints = surf.localPoints();
181 if (debug && (edgeI % 10000 == 0))
183 Pout<<
"Intersecting surface edge " << edgeI
184 <<
" with mesh faces" <<
endl;
186 const edge& e = edges[edgeI];
188 const point& start = localPoints[e.start()];
189 const point& end = localPoints[e.end()];
191 vector edgeNormal(end - start);
192 const scalar edgeMag =
mag(edgeNormal);
193 const vector smallVec = 1
E-9*edgeNormal;
195 edgeNormal /= edgeMag+VSMALL;
210 label faceI = faceTree.shapes().faceLabels()[pHit.index()];
214 cutFace[faceI] =
true;
220 pt = pHit.hitPoint() + smallVec;
222 if (((pt-start) & edgeNormal) >= edgeMag)
232 Pout<<
"Detected an additional " << nAddFaces <<
" faces cut"
235 Pout<<
"Intersected edges of surface with mesh faces in = "
236 << timer.cpuTimeIncrement() <<
" s\n" << endl <<
endl;
246 void Foam::cellClassification::markCells
248 const meshSearch& queryMesh,
257 List<cellInfo> cellInfoList(mesh_.nCells());
260 forAll(piercedFace, faceI)
262 if (piercedFace[faceI])
264 cellInfoList[mesh_.faceOwner()[faceI]] =
267 if (mesh_.isInternalFace(faceI))
269 cellInfoList[mesh_.faceNeighbour()[faceI]] =
280 labelHashSet outsideFacesMap(outsidePts.size() * 6 * 2);
282 forAll(outsidePts, outsidePtI)
285 label cellI = queryMesh.findCell(outsidePts[outsidePtI], -1,
false);
291 "List<cellClassification::cType> markCells"
292 "(const meshSearch&, const boolList&, const pointField&)"
293 ) <<
"outsidePoint " << outsidePts[outsidePtI]
294 <<
" is not inside any cell"
295 <<
nl <<
"It might be on a face or outside the geometry"
304 const labelList& myFaces = mesh_.cells()[cellI];
307 outsideFacesMap.insert(myFaces[myFaceI]);
317 labelList changedFaces(outsideFacesMap.toc());
319 List<cellInfo> changedFacesInfo
325 MeshWave<cellInfo> cellInfoCalc
331 mesh_.globalData().nTotalCells()
335 const List<cellInfo>& allInfo = cellInfoCalc.allCellInfo();
339 label t = allInfo[cellI].type();
345 operator[](cellI) = t;
350 void Foam::cellClassification::classifyPoints
352 const label meshType,
354 List<pointStatus>& pointSide
357 pointSide.setSize(mesh_.nPoints());
359 forAll(mesh_.pointCells(), pointI)
361 const labelList& pCells = mesh_.pointCells()[pointI];
363 pointSide[pointI] = UNSET;
367 label
type = cellType[pCells[i]];
369 if (type == meshType)
371 if (pointSide[pointI] == UNSET)
373 pointSide[pointI] = MESH;
375 else if (pointSide[pointI] == NONMESH)
377 pointSide[pointI] = MIXED;
384 if (pointSide[pointI] == UNSET)
386 pointSide[pointI] = NONMESH;
388 else if (pointSide[pointI] == MESH)
390 pointSide[pointI] = MIXED;
400 bool Foam::cellClassification::usesMixedPointsOnly
402 const List<pointStatus>& pointSide,
406 const faceList& faces = mesh_.faces();
408 const cell& cFaces = mesh_.cells()[cellI];
412 const face&
f = faces[cFaces[cFaceI]];
416 if (pointSide[f[fp]] != MIXED)
428 void Foam::cellClassification::getMeshOutside
430 const label meshType,
435 const labelList& own = mesh_.faceOwner();
436 const labelList& nbr = mesh_.faceNeighbour();
437 const faceList& faces = mesh_.faces();
439 outsideFaces.
setSize(mesh_.nFaces());
440 outsideOwner.setSize(mesh_.nFaces());
445 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
447 label ownType = operator[](own[faceI]);
448 label nbrType = operator[](nbr[faceI]);
450 if (ownType == meshType && nbrType != meshType)
452 outsideFaces[outsideI] = faces[faceI];
453 outsideOwner[outsideI] = own[faceI];
456 else if (ownType != meshType && nbrType == meshType)
458 outsideFaces[outsideI] = faces[faceI];
459 outsideOwner[outsideI] = nbr[faceI];
466 for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
468 if (
operator[](own[faceI]) == meshType)
470 outsideFaces[outsideI] = faces[faceI];
471 outsideOwner[outsideI] = own[faceI];
475 outsideFaces.setSize(outsideI);
476 outsideOwner.setSize(outsideI);
497 markFaces(surfQuery),
513 if (mesh_.nCells() != size())
517 "cellClassification::cellClassification"
518 "(const polyMesh&, const labelList&)"
519 ) <<
"Number of elements of cellType argument is not equal to the"
541 const label meshType,
567 for (label iter = 0; iter < nLayers; iter++)
573 classifyPoints(meshType, newCellType, pointSide);
578 if (pointSide[pointI] == MIXED)
581 const labelList& pCells = mesh_.pointCells()[pointI];
585 label type = newCellType[pCells[i]];
591 newCellType[pCells[i]] = meshType;
605 forAll(newCellType, cellI)
609 if (newCellType[cellI] != meshType)
613 operator[](cellI) = fillType;
626 const label meshType,
630 boolList hasMeshType(mesh_.nPoints(),
false);
634 forAll(mesh_.pointCells(), pointI)
636 const labelList& myCells = mesh_.pointCells()[pointI];
641 label type = operator[](myCells[myCellI]);
643 if (type == meshType)
645 hasMeshType[pointI] =
true;
656 forAll(hasMeshType, pointI)
658 if (hasMeshType[pointI])
660 const labelList& myCells = mesh_.pointCells()[pointI];
664 if (
operator[](myCells[myCellI]) != meshType)
666 operator[](myCells[myCellI]) = fillType;
683 const label meshType,
684 const label fillType,
688 label nTotChanged = 0;
690 for(label iter = 0; iter < maxIter; iter++)
696 classifyPoints(meshType, *
this, pointSide);
703 if (pointSide[pointI] == MIXED)
705 const labelList& pCells = mesh_.pointCells()[pointI];
709 label cellI = pCells[i];
711 if (
operator[](cellI) == meshType)
713 if (usesMixedPointsOnly(pointSide, cellI))
715 operator[](cellI) = fillType;
723 nTotChanged += nChanged;
725 Pout<<
"removeHangingCells : changed " << nChanged
726 <<
" hanging cells" <<
endl;
740 const label meshType,
741 const label fillType,
745 label nTotChanged = 0;
747 for(label iter = 0; iter < maxIter; iter++)
754 getMeshOutside(meshType, outsideFaces, outsideOwner);
767 const labelList& eFaces = edgeFaces[edgeI];
769 if (eFaces.
size() > 2)
776 label patchFaceI = eFaces[i];
778 label ownerCell = outsideOwner[patchFaceI];
780 if (
operator[](ownerCell) == meshType)
782 operator[](ownerCell) = fillType;
792 nTotChanged += nChanged;
794 Pout<<
"fillRegionEdges : changed " << nChanged
795 <<
" cells using multiply connected edges" <<
endl;
809 const label meshType,
810 const label fillType,
814 label nTotChanged = 0;
816 for(label iter = 0; iter < maxIter; iter++)
823 getMeshOutside(meshType, outsideFaces, outsideOwner);
831 fp.checkPointManifold(
false, &nonManifoldPoints);
833 const Map<label>& meshPointMap = fp.meshPointMap();
840 iter != nonManifoldPoints.end();
845 label patchPointI = meshPointMap[iter.key()];
854 label patchFaceI = pFaces[i];
856 label ownerCell = outsideOwner[patchFaceI];
858 if (
operator[](ownerCell) == meshType)
860 operator[](ownerCell) = fillType;
869 nTotChanged += nChanged;
871 Pout<<
"fillRegionPoints : changed " << nChanged
872 <<
" cells using multiply connected points" <<
endl;
886 os <<
"Cells:" << size() <<
endl
887 <<
" notset : " << count(*
this, NOTSET) <<
endl
888 <<
" cut : " << count(*
this, CUT) <<
endl
889 <<
" inside : " << count(*
this, INSIDE) <<
endl
890 <<
" outside : " << count(*
this, OUTSIDE) <<
endl;