53 Foam::label Foam::cellClassification::count
63 if (elems[elemI] == elem)
80 const triSurfaceSearch&
search
85 boolList cutFace(mesh_.nFaces(),
false);
92 forAll(mesh_.edges(), edgeI)
94 if (
debug && (edgeI % 10000 == 0))
96 Pout<<
"Intersecting mesh edge " << edgeI <<
" with surface"
100 const edge&
e = mesh_.edges()[edgeI];
102 const point&
p0 = mesh_.points()[
e.start()];
103 const point& p1 = mesh_.points()[
e.end()];
109 const labelList& myFaces = mesh_.edgeFaces()[edgeI];
113 label facei = myFaces[myFacei];
117 cutFace[facei] =
true;
127 Pout<<
"Intersected edges of mesh with surface in = "
128 << timer.cpuTimeIncrement() <<
" s\n" <<
endl <<
endl;
135 labelList allFaces(mesh_.nFaces() - nCutFaces);
143 allFaces[allFacei++] = facei;
149 Pout<<
"Testing " << allFacei <<
" faces for piercing by surface"
153 treeBoundBox allBb(mesh_.points());
156 scalar tol = 1
e-6 * allBb.avgDim();
168 indexedOctree<treeDataFace> faceTree
170 treeDataFace(
false, mesh_, allFaces),
177 const triSurface& surf =
search.surface();
178 const edgeList& edges = surf.edges();
179 const pointField& localPoints = surf.localPoints();
185 if (
debug && (edgeI % 10000 == 0))
187 Pout<<
"Intersecting surface edge " << edgeI
188 <<
" with mesh faces" <<
endl;
190 const edge&
e = edges[edgeI];
192 const point& start = localPoints[
e.start()];
196 const scalar edgeMag =
mag(edgeNormal);
197 const vector smallVec = 1
e-9*edgeNormal;
199 edgeNormal /= edgeMag+VSMALL;
214 label facei = faceTree.shapes().faceLabels()[pHit.index()];
218 cutFace[facei] =
true;
224 pt = pHit.hitPoint() + smallVec;
226 if (((pt-start) & edgeNormal) >= edgeMag)
236 Pout<<
"Detected an additional " << nAddFaces <<
" faces cut"
239 Pout<<
"Intersected edges of surface with mesh faces in = "
240 << timer.cpuTimeIncrement() <<
" s\n" <<
endl <<
endl;
250 void Foam::cellClassification::markCells
252 const meshSearch& queryMesh,
261 List<cellInfo> cellInfoList(mesh_.nCells());
264 forAll(piercedFace, facei)
266 if (piercedFace[facei])
268 cellInfoList[mesh_.faceOwner()[facei]] =
271 if (mesh_.isInternalFace(facei))
273 cellInfoList[mesh_.faceNeighbour()[facei]] =
284 labelHashSet outsideFacesMap(outsidePts.size() * 6 * 2);
286 forAll(outsidePts, outsidePtI)
289 label celli = queryMesh.findCell(outsidePts[outsidePtI], -1,
false);
294 <<
"outsidePoint " << outsidePts[outsidePtI]
295 <<
" is not inside any cell"
296 <<
nl <<
"It might be on a face or outside the geometry"
305 const labelList& myFaces = mesh_.cells()[celli];
306 outsideFacesMap.insert(myFaces);
315 labelList changedFaces(outsideFacesMap.toc());
317 List<cellInfo> changedFacesInfo
323 MeshWave<cellInfo> cellInfoCalc
329 mesh_.globalData().nTotalCells()+1
333 const List<cellInfo>& allInfo = cellInfoCalc.allCellInfo();
337 label t = allInfo[celli].type();
343 operator[](celli) = t;
348 void Foam::cellClassification::classifyPoints
350 const label meshType,
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;
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];
410 const face&
f = faces[cFaces[cFacei]];
414 if (pointSide[
f[fp]] != MIXED)
426 void Foam::cellClassification::getMeshOutside
428 const label meshType,
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());
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];
454 else if (ownType != meshType && nbrType == meshType)
456 outsideFaces[outsideI] = faces[facei];
457 outsideOwner[outsideI] = nbr[facei];
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];
473 outsideFaces.setSize(outsideI);
474 outsideOwner.setSize(outsideI);
495 markFaces(surfQuery),
511 if (mesh_.nCells() != size())
514 <<
"Number of elements of cellType argument is not equal to the"
536 const label meshType,
562 for (label iter = 0; iter < nLayers; iter++)
568 classifyPoints(meshType, newCellType, pointSide);
573 if (pointSide[pointi] == MIXED)
576 const labelList& pCells = mesh_.pointCells()[pointi];
580 label
type = newCellType[pCells[i]];
586 newCellType[pCells[i]] = meshType;
600 forAll(newCellType, celli)
604 if (newCellType[celli] != meshType)
608 operator[](celli) = fillType;
621 const label meshType,
625 boolList hasMeshType(mesh_.nPoints(),
false);
629 forAll(mesh_.pointCells(), pointi)
631 const labelList& myCells = mesh_.pointCells()[pointi];
636 label
type = operator[](myCells[myCelli]);
638 if (
type == meshType)
640 hasMeshType[pointi] =
true;
651 forAll(hasMeshType, pointi)
653 if (hasMeshType[pointi])
655 const labelList& myCells = mesh_.pointCells()[pointi];
659 if (
operator[](myCells[myCelli]) != meshType)
661 operator[](myCells[myCelli]) = fillType;
678 const label meshType,
679 const label fillType,
683 label nTotChanged = 0;
685 for (label iter = 0; iter < maxIter; iter++)
691 classifyPoints(meshType, *
this, pointSide);
698 if (pointSide[pointi] == MIXED)
700 const labelList& pCells = mesh_.pointCells()[pointi];
704 label celli = pCells[i];
706 if (
operator[](celli) == meshType)
708 if (usesMixedPointsOnly(pointSide, celli))
710 operator[](celli) = fillType;
718 nTotChanged += nChanged;
720 Pout<<
"removeHangingCells : changed " << nChanged
721 <<
" hanging cells" <<
endl;
735 const label meshType,
736 const label fillType,
740 label nTotChanged = 0;
742 for (label iter = 0; iter < maxIter; iter++)
749 getMeshOutside(meshType, outsideFaces, outsideOwner);
762 const labelList& eFaces = edgeFaces[edgeI];
764 if (eFaces.size() > 2)
771 label patchFacei = eFaces[i];
773 label ownerCell = outsideOwner[patchFacei];
775 if (
operator[](ownerCell) == meshType)
777 operator[](ownerCell) = fillType;
787 nTotChanged += nChanged;
789 Pout<<
"fillRegionEdges : changed " << nChanged
790 <<
" cells using multiply connected edges" <<
endl;
804 const label meshType,
805 const label fillType,
809 label nTotChanged = 0;
811 for (label iter = 0; iter < maxIter; ++iter)
818 getMeshOutside(meshType, outsideFaces, outsideOwner);
826 fp.checkPointManifold(
false, &nonManifoldPoints);
828 const Map<label>& meshPointMap = fp.meshPointMap();
832 for (
const label nonManPti : nonManifoldPoints)
835 const label patchPointi = meshPointMap[nonManPti];
842 for (
const label patchFacei :
pFaces)
844 const label ownerCell = outsideOwner[patchFacei];
846 if (
operator[](ownerCell) == meshType)
848 operator[](ownerCell) = fillType;
856 nTotChanged += nChanged;
858 Pout<<
"fillRegionPoints : changed " << nChanged
859 <<
" cells using multiply connected points" <<
endl;
873 os <<
"Cells:" << size() <<
endl
874 <<
" notset : " <<
count(*
this, NOTSET) <<
endl
875 <<
" cut : " <<
count(*
this, CUT) <<
endl
876 <<
" inside : " <<
count(*
this, INSIDE) <<
endl
877 <<
" outside : " <<
count(*
this, OUTSIDE) <<
endl;