55 const labelList& elems,
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()];
195 vector edgeNormal(end - 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;
250void 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;
348void 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;
398bool 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)
426void 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);
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;
Various functions to operate on Lists.
void setSize(const label n)
Alias for resize()
void operator=(const UList< label > &a)
Assignment to UList operator. Takes linear time.
A HashTable to objects of type <T> with a label key.
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
A list of faces which address into the list of points.
const Map< label > & meshPointMap() const
Mesh point map.
bool checkPointManifold(const bool report=false, labelHashSet *setPtr=nullptr) const
Checks primitivePatch for faces sharing point but not edge.
const labelListList & pointFaces() const
Return point-face addressing.
const labelListList & edgeFaces() const
Return edge-face addressing.
iterator end() noexcept
Return an iterator to end traversing the UList.
label size() const noexcept
The number of elements in the UList.
void size(const label n)
Older name for setAddressableSize.
'Cuts' a mesh with a surface.
label fillRegionPoints(const label meshType, const label fillType, const label maxIter)
Find regionPoints and fill all neighbours. Iterate until nothing.
label fillHangingCells(const label meshType, const label fillType, const label maxIter)
Find hanging cells (cells with all points on outside) and set their.
label growSurface(const label meshType, const label fillType)
Sets vertex neighbours of meshType cells to fillType.
void operator=(const cellClassification &)
label trimCutCells(const label nLayers, const label meshType, const label fillType)
void writeStats(Ostream &os) const
Write statistics on cell types to Ostream.
label fillRegionEdges(const label meshType, const label fillType, const label maxIter)
Find regionEdges and fill one neighbour. Iterate until nothing.
Various (local, not parallel) searches on polyMesh; uses (demand driven) octree to search.
Mesh consisting of general polyhedral cells.
label nCells() const noexcept
Number of mesh cells.
Helper class to search on triSurface.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
const volScalarField & p0
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
OBJstream os(runTime.globalPath()/outputName)
cellType
Equivalent to enumeration in "vtkCellType.h" (should be uint8_t)
List< label > labelList
A List of labels.
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
vectorField pointField
pointField is a vectorField.
vector point
Point is a vector.
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
errorManip< error > abort(error &err)
List< bool > boolList
A List of bools.
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
List< edge > edgeList
A List of edges.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
List< face > faceList
A List of faces.
errorManipArg< error, int > exit(error &err, const int errNo=1)
cpuTimeCxx cpuTime
Selection of preferred clock mechanism for the elapsed cpu time.
fileName search(const word &file, const fileName &directory)
Recursively search the given directory for the file.
constexpr char nl
The newline '\n' character (0x0a)
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
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]
#define forAll(list, i)
Loop across all elements in list.