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();
158 point& bbMin = allBb.min();
163 point& bbMax = allBb.max();
168 indexedOctree<treeDataFace> faceTree
170 treeDataFace(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().objectIndex(pHit.index());
218 cutFace[facei] =
true;
224 pt = pHit.point() + 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" <<
nl 296 <<
"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),
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 877 <<
" outside : " <<
count(*
this, OUTSIDE) <<
endl;
cellType
Equivalent to enumeration in "vtkCellType.h" (should be uint8_t)
Various (local, not parallel) searches on polyMesh; uses (demand driven) octree to search...
label trimCutCells(const label nLayers, const label meshType, const label fillType)
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Mixed uniform/non-uniform (eg, after reduction)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
List< edge > edgeList
List of edge.
constexpr char nl
The newline '\n' character (0x0a)
Ostream & endl(Ostream &os)
Add newline and flush stream.
PointIndexHit< point > pointIndexHit
A PointIndexHit with a 3D point.
void operator=(const cellClassification &)
List< labelList > labelListList
List of labelList.
Various functions to operate on Lists.
#define forAll(list, i)
Loop across all elements in list.
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
cellClassification(const polyMesh &mesh, const meshSearch &meshQuery, const triSurfaceSearch &surfQuery, const pointField &outsidePoints)
Construct from mesh and surface and point(s) on outside.
List< face > faceList
List of faces.
Helper class to search on triSurface.
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
bool returnReduceAnd(const bool value, const label comm=UPstream::worldComm)
Perform logical (and) MPI Allreduce on a copy. Uses UPstream::reduceAnd.
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
void setSize(const label n)
Alias for resize()
'Cuts' a mesh with a surface.
void operator=(const UList< label > &list)
Assignment to UList operator. Takes linear time.
label fillRegionPoints(const label meshType, const label fillType, const label maxIter)
Find regionPoints and fill all neighbours. Iterate until nothing.
errorManip< error > abort(error &err)
label size() const noexcept
The number of elements in the container.
int debug
Static debugging option.
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
OBJstream os(runTime.globalPath()/outputName)
defineTypeNameAndDebug(combustionModel, 0)
PrimitivePatch< List< face >, const pointField & > primitiveFacePatch
A PrimitivePatch with List storage for the faces, const reference for the point field.
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::HEX);labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells].reset(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
vector point
Point is a vector.
cpuTimePosix cpuTime
Selection of preferred clock mechanism for the elapsed cpu time.
label growSurface(const label meshType, const label fillType)
Sets vertex neighbours of meshType cells to fillType.
label nCells() const noexcept
Number of mesh cells.
label fillRegionEdges(const label meshType, const label fillType, const label maxIter)
Find regionEdges and fill one neighbour. Iterate until nothing.
void writeStats(Ostream &os) const
Write statistics on cell types to Ostream.
fileName search(const word &file, const fileName &directory)
Recursively search the given directory for the file.
Mesh consisting of general polyhedral cells.
List< label > labelList
A List of labels.
label fillHangingCells(const label meshType, const label fillType, const label maxIter)
Find hanging cells (cells with all points on outside) and set their.
List< bool > boolList
A List of bools.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
const volScalarField & p0