53 Foam::label Foam::shortestPathSet::findMinFace
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];
113 label faceI = cFaces2[i];
122 if (isLeakPoint[
f[fp]])
129 if (nLeak < minLeakPoints)
131 minLeakPoints = nLeak;
143 void Foam::shortestPathSet::calculateDistance
146 const polyMesh&
mesh,
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)
170 faceDist.append(topoDistanceData<label>(0, 123));
180 topoDistanceData<label>
194 const fvMesh& fm = refCast<const fvMesh>(
mesh);
203 fm.time().timeName(),
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
299 if (isLeakPoint[
f[fp]] && isLeakPoint[
f[nextFp]])
302 if (isLeakFace.set(facei))
315 const polyMesh&
mesh,
316 const bitSet& isLeakCell
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,
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,
430 if (celli != insideCelli && celli != outsideCelli)
432 if (isLeakCell[celli])
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])
469 if (isLeakFace[facei])
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;
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);
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;
707 frontCellI =
pp.faceCells()[i];
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)
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;
955 bool markLeakPath =
false;
958 for (iter = 0; iter < maxIter; iter++)
965 const label oldNSamplingPts(samplingPts.size());
967 bool foundPath = genSingleLeakPath
1000 samplingCurveDist.size()
1001 ? samplingCurveDist.last()
1013 if (!foundPath && !markLeakPath)
1021 if (nLeakFaces > nOldLeakFaces)
1027 markLeakPath =
false;
1036 if (markLeakPath && !foundPath)
1054 markLeakPath =
true;
1066 Pout<<
"Writing new isLeakCell to " << str.
name() <<
endl;
1076 Pout<<
"Writing new leak-path to " << str.
name() <<
endl;
1079 label samplei = oldNSamplingPts+1;
1080 samplei < samplingPts.size();
1084 Pout<<
" passing through cell " 1085 << samplingCells[samplei]
1087 <<
" distance:" << samplingCurveDist[samplei]
1092 samplingPts[samplei-1],
1093 samplingPts[samplei]
1102 const fvMesh& fm = refCast<const fvMesh>(
mesh);
1104 const_cast<fvMesh&
>(fm).setInstance(fm.time().timeName());
1110 fm.time().timeName(),
1118 for (
const label celli : isLeakCell)
1120 fld[celli] = scalar(1);
1122 fld.correctBoundaryConditions();
1126 if (maxIter > 1 && iter == maxIter)
1129 <<
" leak paths" <<
nl <<
"This can cause problems when using the" 1130 <<
" paths to close leaks" <<
endl;
1135 void Foam::shortestPathSet::genSamples
1137 const bool markLeakPath,
1138 const label maxIter,
1139 const polyMesh&
mesh,
1145 DynamicList<point> samplingPts;
1146 DynamicList<label> samplingCells;
1147 DynamicList<label> samplingFaces;
1148 DynamicList<label> samplingSegments;
1149 DynamicList<scalar> samplingCurveDist;
1156 for (label patchi : wallPatches)
1158 const polyPatch&
pp =
pbm[patchi];
1161 isBlockedPoint.set(
pp[i]);
1170 isBlockedPoint.set(
mesh.
faces()[facei]);
1178 orEqOp<unsigned int>(),
1186 for (
const auto& pointi : isBlockedPoint)
1190 Pout<<
"Writing " << str.nVertices()
1191 <<
" points to " << str.
name() <<
endl;
1196 bitSet isLeakPoint(isBlockedPoint);
1202 label prevSegmenti = 0;
1203 scalar prevDistance = 0.0;
1205 for (
auto insidePoint : insidePoints_)
1209 for (
auto outsidePoint : outsidePoints_)
1211 const label nOldSamples = samplingSegments.size();
1235 label maxSegment = 0;
1236 scalar maxDistance = 0.0;
1237 for (label i = nOldSamples; i < samplingSegments.size(); ++i)
1239 samplingSegments[i] += prevSegmenti;
1240 maxSegment =
max(maxSegment, samplingSegments[i]);
1241 samplingCurveDist[i] += prevDistance;
1242 maxDistance =
max(maxDistance, samplingCurveDist[i]);
1244 prevSegmenti =
returnReduce(maxSegment, maxOp<label>());
1245 prevDistance =
returnReduce(maxDistance, maxOp<scalar>());
1251 erodeFaceSet(
mesh, isBlockedPoint, isLeakFace);
1253 leakFaces_ = isLeakFace.sortedToc();
1263 Pout<<
"Writing " << leakFaces.size()
1264 <<
" faces to " << str.
name() <<
endl;
1268 samplingPts.shrink();
1269 samplingCells.shrink();
1270 samplingFaces.shrink();
1271 samplingSegments.shrink();
1272 samplingCurveDist.shrink();
1277 std::move(samplingPts),
1278 std::move(samplingCells),
1279 std::move(samplingFaces),
1280 std::move(samplingSegments),
1281 std::move(samplingCurveDist)
1299 const bool markLeakPath,
1300 const label maxIter,
1309 outsidePoints_(outsidePoints)
1321 Info<<
"shortestPathSet : Writing blocked faces to " 1322 << outputDir <<
endl;
1346 setPatch.localFaces(),
1347 setPatch.meshPoints(),
1348 setPatch.meshPointMap(),
1351 uniqueMeshPointLabels,
1366 (outputDir /
"blockedFace"),
1377 setPatch.localPoints(),
1378 setPatch.localFaces(),
1379 (outputDir /
"blockedFace"),
1410 const label maxIter(
dict.getOrDefault<label>(
"maxIter", 50));
1411 const bool markLeakPath(
dict.getOrDefault(
"markLeakPath",
true));
1419 if (!
pp.coupled() && !isA<emptyPolyPatch>(
pp))
1421 wallPatches.append(patchi);
1425 genSamples(markLeakPath, maxIter,
mesh, wallPatches,
bitSet());
List< ReturnType > get(const UPtrList< T > &list, const AccessOp &aop)
List of values generated by applying the access operation to each list item.
shortestPathSet(const word &name, const polyMesh &mesh, const meshSearch &searchEngine, const word &axis, const bool markLeakPath, const label maxIter, const labelUList &wallPatches, const pointField &insidePoints, const pointField &outsidePoints, const boolList &blockedFace)
Construct from components. blockedFace is an optional specification.
Various (local, not parallel) searches on polyMesh; uses (demand driven) octree to search...
const polyBoundaryMesh & pbm
A class for handling file names.
List< wallPoints > allCellInfo(mesh_.nCells())
errorManipArg< error, int > exit(error &err, const int errNo=1)
vtk::lineWriter writer(edgeCentres, edgeList::null(), fileName(aMesh.time().globalPath()/"finiteArea-edgesCentres"))
void transfer(List< T > &list)
Transfer the contents of the argument List into this list and annul the argument list.
label nPoints() const noexcept
Number of mesh points.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
fileName timePath() const
Return current time path = path/timeName.
void append(const T &val)
Append an element at the end of the list.
virtual const labelList & faceNeighbour() const
Return face neighbour.
const bitSet isBlockedFace(intersectedFaces())
const mapDistribute & globalEdgeSlavesMap() const
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
constexpr char nl
The newline '\n' character (0x0a)
Write faces/points (optionally with fields) as a vtp file or a legacy vtk file.
Ostream & endl(Ostream &os)
Add newline and flush stream.
static bool & parRun() noexcept
Test if this a parallel run.
const scalarList & distance() const noexcept
Return the cumulative distance.
scalar distance(const vector &p1, const vector &p2)
labelList findIndices(const ListType &input, typename ListType::const_reference val, label start=0)
Linear search to find all occurrences of given element.
const cellList & cells() const
Ignore writing from objectRegistry::writeObject()
const dimensionSet dimless
Dimensionless.
const Time & time() const
Return the top-level database.
label nFaces() const noexcept
Number of mesh faces.
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
A PrimitivePatch with an IndirectList for the faces, const reference for the point field...
virtual const fileName & name() const override
Get the name of the output serial stream. (eg, the name of the Fstream file name) ...
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
List< labelList > labelListList
List of labelList.
Macros for easy insertion into run-time selection tables.
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
UList< label > labelUList
A UList of labels.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
virtual const pointField & points() const
Return raw points.
#define forAll(list, i)
Loop across all elements in list.
label fcIndex(const label i) const noexcept
The forward circular index. The next index in the list which returns to the first at the end of the l...
GeometricField< scalar, fvPatchField, volMesh > volScalarField
List< face > faceList
List of faces.
void write(vtk::formatter &fmt, const Type &val, const label n=1)
Component-wise write of a value (N times)
A list of faces which address into the list of points.
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.
const fileName & pointsInstance() const
Return the current instance directory for points.
bool mkDir(const fileName &pathName, mode_t mode=0777)
Make a directory and return an error if it could not be created.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
Holds list of sampling points which is filled at construction time. Various implementations of this b...
fileName globalPath() const
Return global path for the case = rootPath/globalCaseName. Same as TimePaths::globalPath() ...
const polyBoundaryMesh & boundaryMesh() const noexcept
Return boundary mesh.
A class for handling words, derived from Foam::string.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
const Time & time() const noexcept
Return time registry.
static void combineReduce(T &value, const CombineOp &cop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce inplace (cf. MPI Allreduce) applying cop to inplace combine value from different processors...
virtual const labelList & faceOwner() const
Return face owner.
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
const labelListList & edgeFaces() const
Return edge-face addressing.
const globalMeshData & globalData() const
Return parallel info (demand-driven)
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
label nInternalFaces() const noexcept
Number of internal faces.
virtual const faceList & faces() const
Return raw faces.
FaceCellWave< wallPoints, wallPoints::trackData > wallDistCalc(mesh_, changedFaces, faceDist, allFaceInfo, allCellInfo, 0, td)
const vectorField & cellCentres() const
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO...
const polyMesh & mesh() const noexcept
int debug
Static debugging option.
defineTypeNameAndDebug(combustionModel, 0)
static bool clean(std::string &str)
Cleanup filename string, possibly applies other transformations such as changing the path separator e...
label nTotalCells() const noexcept
Total global number of mesh cells.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< ' ';}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< ' ';}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< ' ';}gmvFile<< nl;for(const word &name :lagrangianScalarNames){ IOField< scalar > fld(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
const indirectPrimitivePatch & coupledPatch() const
Return patch of all coupled faces.
const vectorField & faceCentres() const
Calculates points shared by more than two processor patches or cyclic patches.
PrimitivePatch< UIndirectList< face >, const pointField & > uindirectPrimitivePatch
A PrimitivePatch with UIndirectList for the faces, const reference for the point field.
bool empty() const noexcept
True if the list is empty (ie, size() is zero).
vector point
Point is a vector.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
#define WarningInFunction
Report a warning using Foam::Warning.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
label nCells() const noexcept
Number of mesh cells.
static bool master(const label communicator=worldComm)
True if process corresponds to the master rank in the communicator.
const polyBoundaryMesh & patches
static word outputPrefix
Directory prefix.
label findCell(const point &p, const cellDecomposition=CELL_TETS) const
Find cell enclosing this location and return index.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Pointer management similar to std::unique_ptr, with some additional methods and type checking...
void reset()
Clear all bits but do not adjust the addressable size.
Mesh consisting of general polyhedral cells.
List< label > labelList
A List of labels.
A patch is a list of labels that address the faces in the global face list.
label nBoundaryFaces() const noexcept
Number of boundary faces (== nFaces - nInternalFaces)
bool returnReduceOr(const bool value, const label comm=UPstream::worldComm)
Perform logical (or) MPI Allreduce on a copy. Uses UPstream::reduceOr.
List< bool > boolList
A List of bools.
A List with indirect addressing.
List< wallPoints > allFaceInfo(mesh_.nFaces())
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
label size() const noexcept
Number of entries.
void reduce(T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce inplace (cf. MPI Allreduce) using linear/tree communication schedule.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
addToRunTimeSelectionTable(functionObject, pointHistory, dictionary)
static constexpr const zero Zero
Global zero (0)