53 Foam::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;
143 void 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;
235 void Foam::shortestPathSet::sync
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();
284 bool Foam::shortestPathSet::touchesWall
286 const polyMesh&
mesh,
290 const bitSet& isLeakPoint
297 const label nextFp =
f.fcIndex(fp);
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);
376 bool 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)
780 Foam::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);
911 void 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.valid() ? isBlockedFace() : isBoundaryFace),
1000 samplingCurveDist.size()
1001 ? samplingCurveDist.last()
1013 if (!foundPath && !markLeakPath)
1021 if (nLeakFaces > nOldLeakFaces)
1026 if (isBlockedFace.valid())
1028 isBlockedFace.clear();
1030 markLeakPath =
false;
1039 if (markLeakPath && !foundPath)
1051 if (!isBlockedFace.valid())
1054 isBlockedFace.reset(
new bitSet(isBoundaryFace));
1057 markLeakPath =
true;
1069 Pout<<
"Writing new isLeakCell to " << str.
name() <<
endl;
1072 str.write(leakCcs[i]);
1082 Pout<<
"Writing new leak-path to " << str.
name() <<
endl;
1085 label samplei = oldNSamplingPts+1;
1086 samplei < samplingPts.size();
1090 Pout<<
" passing through cell "
1091 << samplingCells[samplei]
1093 <<
" distance:" << samplingCurveDist[samplei]
1100 samplingPts[samplei-1],
1101 samplingPts[samplei]
1111 const fvMesh& fm = refCast<const fvMesh>(
mesh);
1113 const_cast<fvMesh&
>(fm).setInstance(fm.time().timeName());
1119 fm.time().timeName(),
1127 forAll(isLeakCell, celli)
1129 fld[celli] = isLeakCell[celli];
1131 fld.correctBoundaryConditions();
1135 if (maxIter > 1 && iter == maxIter)
1138 <<
" leak paths" <<
nl <<
"This can cause problems when using the"
1139 <<
" paths to close leaks" <<
endl;
1144 void Foam::shortestPathSet::genSamples
1146 const bool markLeakPath,
1147 const label maxIter,
1148 const polyMesh&
mesh,
1150 const bitSet& isBlockedFace
1154 DynamicList<point> samplingPts;
1155 DynamicList<label> samplingCells;
1156 DynamicList<label> samplingFaces;
1157 DynamicList<label> samplingSegments;
1158 DynamicList<scalar> samplingCurveDist;
1165 for (label patchi : wallPatches)
1167 const polyPatch& pp = pbm[patchi];
1170 isBlockedPoint.set(pp[i]);
1175 forAll(isBlockedFace, facei)
1177 if (isBlockedFace[facei])
1179 isBlockedPoint.set(
mesh.
faces()[facei]);
1187 orEqOp<unsigned int>(),
1195 for (
const auto& pointi : isBlockedPoint)
1199 Pout<<
"Writing " << str.nVertices() <<
" points to " << str.
name()
1205 bitSet isLeakPoint(isBlockedPoint);
1211 label prevSegmenti = 0;
1212 scalar prevDistance = 0.0;
1214 for (
auto insidePoint : insidePoints_)
1218 for (
auto outsidePoint : outsidePoints_)
1220 const label nOldSamples = samplingSegments.size();
1244 label maxSegment = 0;
1245 scalar maxDistance = 0.0;
1246 for (label i = nOldSamples; i < samplingSegments.size(); ++i)
1248 samplingSegments[i] += prevSegmenti;
1249 maxSegment =
max(maxSegment, samplingSegments[i]);
1250 samplingCurveDist[i] += prevDistance;
1251 maxDistance =
max(maxDistance, samplingCurveDist[i]);
1253 prevSegmenti =
returnReduce(maxSegment, maxOp<label>());
1254 prevDistance =
returnReduce(maxDistance, maxOp<scalar>());
1260 erodeFaceSet(
mesh, isBlockedPoint, isLeakFace);
1262 leakFaces_ = isLeakFace.sortedToc();
1272 Pout<<
"Writing " << leakFaces.size() <<
" faces to " << str.
name()
1277 samplingPts.shrink();
1278 samplingCells.shrink();
1279 samplingFaces.shrink();
1280 samplingSegments.shrink();
1281 samplingCurveDist.shrink();
1286 std::move(samplingPts),
1287 std::move(samplingCells),
1288 std::move(samplingFaces),
1289 std::move(samplingSegments),
1290 std::move(samplingCurveDist)
1308 const bool markLeakPath,
1309 const label maxIter,
1318 outsidePoints_(outsidePoints)
1330 Info<<
"shortestPathSet : Writing blocked faces to "
1331 << outputDir <<
endl;
1355 setPatch.localFaces(),
1356 setPatch.meshPoints(),
1357 setPatch.meshPointMap(),
1360 uniqueMeshPointLabels,
1375 (outputDir /
"blockedFace"),
1386 setPatch.localPoints(),
1387 setPatch.localFaces(),
1388 (outputDir /
"blockedFace"),
1428 if (!pp.
coupled() && !isA<emptyPolyPatch>(pp))
1430 wallPatches.append(patchi);
1434 genSamples(markLeakPath, maxIter,
mesh, wallPatches,
bitSet());