53Foam::label Foam::shortestPathSet::findMinFace
57 const List<topoDistanceData<label>>& allFaceInfo,
58 const bitSet& isLeakPoint,
59 const bool distanceMode,
63 const cell& cFaces2 =
mesh.
cells()[cellI];
72 label faceI = cFaces2[i];
73 const topoDistanceData<label>& info = allFaceInfo[faceI];
74 if (info.distance() < minDist)
76 minDist = info.distance();
80 else if (info.distance() == minDist)
92 scalar minDist2 = ROOTVGREAT;
95 label faceI = cFaces2[i];
96 if (allFaceInfo[faceI].
distance() == minDist)
113 label faceI = cFaces2[i];
114 if (allFaceInfo[faceI].
distance() == minDist)
122 if (isLeakPoint[
f[fp]])
129 if (nLeak < minLeakPoints)
131 minLeakPoints = nLeak;
143void Foam::shortestPathSet::calculateDistance
146 const polyMesh&
mesh,
149 List<topoDistanceData<label>>& allFaceInfo,
150 List<topoDistanceData<label>>& allCellInfo
153 int dummyTrackData = 0;
156 DynamicList<topoDistanceData<label>> faceDist;
157 DynamicList<label> cFaces1;
162 faceDist.reserve(cFaces.size());
163 cFaces1.reserve(cFaces.size());
165 for (label facei : cFaces)
167 if (!allFaceInfo[facei].valid(dummyTrackData))
170 faceDist.append(topoDistanceData<label>(0, 123));
180 topoDistanceData<label>
194 const fvMesh& fm = refCast<const fvMesh>(
mesh);
203 fm.time().timeName(),
211 forAll(allCellInfo, celli)
213 fld[celli] = allCellInfo[celli].distance();
218 SubList<topoDistanceData<label>>
p(pp.patchSlice(allFaceInfo));
222 pfld[i] = 1.0*
p[i].distance();
224 fld.boundaryFieldRef()[patchi] == pfld;
228 Pout<<
"Writing distance field for initial cell " << cellI
229 <<
" to " <<
fld.objectPath() <<
endl;
237 const polyMesh&
mesh,
242 bool& findMinDistance
249 orEqOp<unsigned int>(),
256 orEqOp<unsigned int>()
260 typedef Tuple2<label, Tuple2<point, bool>> ProcData;
265 Tuple2<point, bool>(origin, findMinDistance)
270 [](ProcData&
x,
const ProcData&
y)
278 origin = searchData.second().first();
279 findMinDistance = searchData.second().second();
284bool Foam::shortestPathSet::touchesWall
286 const polyMesh&
mesh,
290 const bitSet& isLeakPoint
299 if (isLeakPoint[
f[fp]] && isLeakPoint[
f[nextFp]])
302 if (isLeakFace.set(facei))
315 const polyMesh&
mesh,
316 const bitSet& isLeakCell
330 for (
const polyPatch& pp :
patches)
338 for (
const label celli : faceCells)
340 nbrLeakCell[bFacei] = isLeakCell[celli];
359 if (isLeakCell[own[facei]] && isLeakCell[nbr[facei]])
361 isLeakFace.set(facei);
369 isLeakFace.set(facei);
376bool Foam::shortestPathSet::genSingleLeakPath
378 const bool markLeakPath,
380 const polyMesh&
mesh,
381 const bitSet& isBlockedFace,
382 const point& insidePoint,
383 const label insideCelli,
384 const point& outsidePoint,
385 const label outsideCelli,
388 const scalar trackLength,
389 DynamicList<point>& samplingPts,
390 DynamicList<label>& samplingCells,
391 DynamicList<label>& samplingFaces,
392 DynamicList<label>& samplingSegments,
393 DynamicList<scalar>& samplingCurveDist,
401 List<topoDistanceData<label>>& allFaceInfo,
402 List<topoDistanceData<label>>& allCellInfo
410 allFaceInfo = topoDistanceData<label>();
412 allCellInfo = topoDistanceData<label>();
415 forAll(isBlockedFace, facei)
417 if (isBlockedFace[facei])
419 allFaceInfo[facei] = maxData;
430 if (celli != insideCelli && celli != outsideCelli)
432 if (isLeakCell[celli])
434 allCellInfo[celli] = maxData;
446 bitSet isLeakCellWithout(isLeakCell);
447 if (insideCelli != -1)
449 isLeakCellWithout.unset(insideCelli);
451 if (outsideCelli != -1)
453 isLeakCellWithout.unset(outsideCelli);
455 const bitSet isPathFace(pathFaces(
mesh, isLeakCellWithout));
458 if (isPathFace[facei])
460 allFaceInfo[facei] = maxData;
469 if (isLeakFace[facei])
471 allFaceInfo[facei] = maxData;
478 calculateDistance(iter,
mesh, insideCelli, allFaceInfo, allCellInfo);
486 bool targetFound =
false;
487 if (outsideCelli != -1)
490 targetFound = allCellInfo[outsideCelli].valid(dummyTrackData);
491 if (iter == 0 && !targetFound)
494 <<
"Point " << outsidePoint
495 <<
" not reachable by walk from " << insidePoint
496 <<
". Probably mesh has island/regions."
497 <<
" Skipped route detection." <<
endl;
500 reduce(targetFound, orOp<bool>());
518 label frontCellI = outsideCelli;
519 point origin(outsidePoint);
520 bool findMinDistance =
true;
527 label frontFaceI = -1;
530 if (frontCellI != -1)
549 frontFaceI = findMinFace
561 bitSet isNewLeakPoint(isLeakPoint);
564 if (isBlockedFace.size() && isBlockedFace[frontFaceI])
575 if (nbrCellI == frontCellI)
580 if (nbrCellI == insideCelli)
587 frontCellI = nbrCellI;
590 frontFaceI = findMinFace
600 const topoDistanceData<label>& cInfo = allCellInfo[frontCellI];
601 const topoDistanceData<label>& fInfo = allFaceInfo[frontFaceI];
603 if (fInfo.distance() <= cInfo.distance())
607 samplingCells.append(frontCellI);
608 samplingFaces.append(-1);
609 samplingSegments.append(iter);
610 samplingCurveDist.append
612 trackLength+cInfo.distance()
614 isLeakCell.set(frontCellI);
630 isNewLeakPoint.set(
mesh.
faces()[frontFaceI]);
632 findMinDistance =
false;
636 isLeakPoint.
transfer(isNewLeakPoint);
660 if (frontFaceI != -1)
667 if (frontCellI != -1)
669 minCellDistance = allCellInfo[frontCellI].distance();
671 reduce(minCellDistance, minOp<label>());
690 const label oldFrontFaceI = frontFaceI;
695 const polyPatch& pp = pbm[patchI];
698 label faceI = pp.start()+i;
703 && (isBlockedFace.empty() || !isBlockedFace[faceI])
707 frontCellI = pp.faceCells()[i];
715 && allCellInfo[frontCellI].
distance() < minCellDistance
718 const topoDistanceData<label>& cInfo = allCellInfo[frontCellI];
721 samplingCells.append(frontCellI);
722 samplingFaces.append(-1);
723 samplingSegments.append(iter);
724 samplingCurveDist.append
726 trackLength+cInfo.distance()
728 isLeakCell.set(frontCellI);
744 isLeakPoint.set(
mesh.
faces()[frontFaceI]);
746 findMinDistance =
false;
753 if (insideCelli != -1 && frontCellI == insideCelli)
780Foam::label Foam::shortestPathSet::erodeFaceSet
782 const polyMesh&
mesh,
783 const bitSet& isBlockedPoint,
794 <<
" isBlockedPoint:" << isBlockedPoint.size()
795 <<
" isLeakFace:" << isLeakFace.size()
804 label nTotalEroded = 0;
808 bitSet newIsLeakFace(isLeakFace);
812 const labelList meshFaceIDs(isLeakFace.toc());
815 UIndirectList<face>(
mesh.
faces(), meshFaceIDs),
824 nEdgeFaces[edgei] = edgeFaces[edgei].size();
830 bitSet sameEdgeOrientation;
843 UIndirectList<label>(coupledNEdgeFaces, coupledEdges) =
844 UIndirectList<label>(nEdgeFaces, patchEdges);
850 globalData.globalEdgeSlaves(),
851 globalData.globalEdgeTransformedSlaves(),
857 UIndirectList<label>(nEdgeFaces, patchEdges) =
858 UIndirectList<label>(coupledNEdgeFaces, coupledEdges);
864 if (nEdgeFaces[edgei] == 1)
866 const edge&
e = pp.edges()[edgei];
867 const label mp0 = pp.meshPoints()[
e[0]];
868 const label mp1 = pp.meshPoints()[
e[1]];
870 if (!isBlockedPoint[mp0] || !isBlockedPoint[mp1])
873 const label patchFacei = edgeFaces[edgei][0];
874 const label meshFacei = pp.addressing()[patchFacei];
880 if (newIsLeakFace.unset(meshFacei))
888 reduce(nEroded, sumOp<label>());
889 nTotalEroded += nEroded;
901 orEqOp<unsigned int>()
904 isLeakFace = std::move(newIsLeakFace);
911void Foam::shortestPathSet::genSamples
913 const bool addLeakPath,
915 const polyMesh&
mesh,
916 const bitSet& isBoundaryFace,
917 const point& insidePoint,
918 const label insideCelli,
919 const point& outsidePoint,
921 DynamicList<point>& samplingPts,
922 DynamicList<label>& samplingCells,
923 DynamicList<label>& samplingFaces,
924 DynamicList<label>& samplingSegments,
925 DynamicList<scalar>& samplingCurveDist,
944 scalar trackLength = 0;
946 List<topoDistanceData<label>> allFaceInfo(
mesh.
nFaces());
947 List<topoDistanceData<label>> allCellInfo(
mesh.
nCells());
952 autoPtr<bitSet> isBlockedFace;
955 bool markLeakPath =
false;
958 for (iter = 0; iter < maxIter; iter++)
965 const label oldNSamplingPts(samplingPts.size());
967 bool foundPath = genSingleLeakPath
972 (isBlockedFace ? *isBlockedFace : isBoundaryFace),
1000 samplingCurveDist.size()
1001 ? samplingCurveDist.last()
1013 if (!foundPath && !markLeakPath)
1021 if (nLeakFaces > nOldLeakFaces)
1026 isBlockedFace.reset(
nullptr);
1027 markLeakPath =
false;
1036 if (markLeakPath && !foundPath)
1051 isBlockedFace.reset(
new bitSet(isBoundaryFace));
1054 markLeakPath =
true;
1066 Pout<<
"Writing new isLeakCell to " << str.
name() <<
endl;
1069 str.write(leakCcs[i]);
1079 Pout<<
"Writing new leak-path to " << str.
name() <<
endl;
1082 label samplei = oldNSamplingPts+1;
1083 samplei < samplingPts.size();
1087 Pout<<
" passing through cell "
1088 << samplingCells[samplei]
1090 <<
" distance:" << samplingCurveDist[samplei]
1097 samplingPts[samplei-1],
1098 samplingPts[samplei]
1108 const fvMesh& fm = refCast<const fvMesh>(
mesh);
1110 const_cast<fvMesh&
>(fm).setInstance(fm.time().timeName());
1116 fm.time().timeName(),
1124 forAll(isLeakCell, celli)
1126 fld[celli] = isLeakCell[celli];
1128 fld.correctBoundaryConditions();
1132 if (maxIter > 1 && iter == maxIter)
1135 <<
" leak paths" <<
nl <<
"This can cause problems when using the"
1136 <<
" paths to close leaks" <<
endl;
1141void Foam::shortestPathSet::genSamples
1143 const bool markLeakPath,
1144 const label maxIter,
1145 const polyMesh&
mesh,
1147 const bitSet& isBlockedFace
1151 DynamicList<point> samplingPts;
1152 DynamicList<label> samplingCells;
1153 DynamicList<label> samplingFaces;
1154 DynamicList<label> samplingSegments;
1155 DynamicList<scalar> samplingCurveDist;
1162 for (label patchi : wallPatches)
1164 const polyPatch& pp = pbm[patchi];
1167 isBlockedPoint.set(pp[i]);
1172 forAll(isBlockedFace, facei)
1174 if (isBlockedFace[facei])
1176 isBlockedPoint.set(
mesh.
faces()[facei]);
1184 orEqOp<unsigned int>(),
1192 for (
const auto& pointi : isBlockedPoint)
1196 Pout<<
"Writing " << str.nVertices() <<
" points to " << str.
name()
1202 bitSet isLeakPoint(isBlockedPoint);
1208 label prevSegmenti = 0;
1209 scalar prevDistance = 0.0;
1211 for (
auto insidePoint : insidePoints_)
1215 for (
auto outsidePoint : outsidePoints_)
1217 const label nOldSamples = samplingSegments.size();
1241 label maxSegment = 0;
1242 scalar maxDistance = 0.0;
1243 for (label i = nOldSamples; i < samplingSegments.size(); ++i)
1245 samplingSegments[i] += prevSegmenti;
1246 maxSegment =
max(maxSegment, samplingSegments[i]);
1247 samplingCurveDist[i] += prevDistance;
1248 maxDistance =
max(maxDistance, samplingCurveDist[i]);
1250 prevSegmenti =
returnReduce(maxSegment, maxOp<label>());
1251 prevDistance =
returnReduce(maxDistance, maxOp<scalar>());
1257 erodeFaceSet(
mesh, isBlockedPoint, isLeakFace);
1259 leakFaces_ = isLeakFace.sortedToc();
1269 Pout<<
"Writing " << leakFaces.size() <<
" faces to " << str.
name()
1274 samplingPts.shrink();
1275 samplingCells.shrink();
1276 samplingFaces.shrink();
1277 samplingSegments.shrink();
1278 samplingCurveDist.shrink();
1283 std::move(samplingPts),
1284 std::move(samplingCells),
1285 std::move(samplingFaces),
1286 std::move(samplingSegments),
1287 std::move(samplingCurveDist)
1305 const bool markLeakPath,
1306 const label maxIter,
1315 outsidePoints_(outsidePoints)
1327 Info<<
"shortestPathSet : Writing blocked faces to "
1328 << outputDir <<
endl;
1357 uniqueMeshPointLabels,
1372 (outputDir /
"blockedFace"),
1385 (outputDir /
"blockedFace"),
1425 if (!pp.
coupled() && !isA<emptyPolyPatch>(pp))
1427 wallPatches.
append(patchi);
1431 genSamples(markLeakPath, maxIter,
mesh, wallPatches,
bitSet());
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
vtk::internalMeshWriter writer(topoMesh, topoCells, vtk::formatType::INLINE_ASCII, runTime.path()/"blockTopology")
Info<< nl<< "Wrote faMesh in vtk format: "<< writer.output().name()<< nl;}{ vtk::lineWriter writer(aMesh.points(), aMesh.edges(), fileName(aMesh.mesh().time().globalPath()/"finiteArea-edges"));writer.writeGeometry();writer.beginCellData(4);writer.writeProcIDs();{ Field< scalar > fld(faMeshTools::flattenEdgeField(aMesh.magLe(), true))
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
void append(const T &val)
Copy append an element to the end of this list.
A List with indirect addressing.
void transfer(List< T > &list)
void append(const T &val)
Append an element at the end of the list.
virtual const fileName & name() const
Get the name of the stream.
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
A list of faces which address into the list of points.
const Map< label > & meshPointMap() const
Mesh point map.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
const Field< point_type > & localPoints() const
Return pointField of points in patch.
const List< face_type > & localFaces() const
Return patch faces addressing into local point list.
static void combineAllGather(const List< commsStruct > &comms, T &value, const CombineOp &cop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
fileName globalPath() const
Return global path for the case.
fileName timePath() const
Return current time path.
label fcIndex(const label i) const noexcept
static bool & parRun() noexcept
Test if this a parallel run.
label size() const noexcept
The number of elements in the list.
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
const scalarList & distance() const noexcept
Return the cumulative distance.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
A class for handling file names.
static bool clean(std::string &str)
static word outputPrefix
Directory prefix.
void sync()
Do all: synchronise all IOFields and objectRegistry.
const Time & time() const
Return the top-level database.
const mapDistribute & globalEdgeSlavesMap() const
const indirectPrimitivePatch & coupledPatch() const
Return patch of all coupled faces.
label nTotalCells() const noexcept
Return total number of cells in decomposed mesh.
Calculates points shared by more than two processor patches or cyclic patches.
Various (local, not parallel) searches on polyMesh; uses (demand driven) octree to search.
const Time & time() const noexcept
Return time registry.
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
Mesh consisting of general polyhedral cells.
virtual const faceList & faces() const
Return raw faces.
label findCell(const point &p, const cellDecomposition=CELL_TETS) const
Find cell enclosing this location and return index.
virtual const labelList & faceOwner() const
Return face owner.
const globalMeshData & globalData() const
Return parallel info.
const fileName & pointsInstance() const
Return the current instance directory for points.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
virtual const labelList & faceNeighbour() const
Return face neighbour.
virtual const pointField & points() const
Return raw points.
A patch is a list of labels that address the faces in the global face list.
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
label nBoundaryFaces() const noexcept
Number of boundary faces (== nFaces - nInternalFaces)
const vectorField & faceCentres() const
label nInternalFaces() const noexcept
Number of internal faces.
const vectorField & cellCentres() const
label nPoints() const noexcept
Number of mesh points.
label nCells() const noexcept
Number of mesh cells.
label nFaces() const noexcept
Number of mesh faces.
const cellList & cells() const
Holds list of sampling points which is filled at construction time. Various implementations of this b...
const polyMesh & mesh() const noexcept
Finds shortest path (in terms of cell centres) to walk on mesh from any point in insidePoints to any ...
splitCell * master() const
Write faces/points (optionally with fields) as a vtp file or a legacy vtk file.
A class for handling words, derived from Foam::string.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
const polyBoundaryMesh & patches
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
#define WarningInFunction
Report a warning using Foam::Warning.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
scalar distance(const vector &p1, const vector &p2)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
const dimensionSet dimless
Dimensionless.
List< label > labelList
A List of labels.
bool mkDir(const fileName &pathName, mode_t mode=0777)
Make a directory and return an error if it could not be created.
line< point, const point & > linePointRef
A line using referred points.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
vector point
Point is a vector.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
PrimitivePatch< UIndirectList< face >, const pointField & > uindirectPrimitivePatch
A PrimitivePatch with UIndirectList for the faces, const reference for the point field.
List< labelList > labelListList
A List of labelList.
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
List< bool > boolList
A List of bools.
static constexpr const zero Zero
Global zero (0)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
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.
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
A PrimitivePatch with an IndirectList for the faces, const reference for the point field.
labelList findIndices(const ListType &input, typename ListType::const_reference val, label start=0)
Linear search to find all occurrences of given element.
errorManipArg< error, int > exit(error &err, const int errNo=1)
UList< label > labelUList
A UList of labels.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
constexpr char nl
The newline '\n' character (0x0a)
#define forAll(list, i)
Loop across all elements in list.