84 scalar dist = nBins/(
max -
min);
98 else if (val >=
max - SMALL)
104 index = label((val - min)*dist);
106 if ((index < 0) || (index >= nBins))
109 <<
"value " << val <<
" at index " << i
110 <<
" outside range " <<
min <<
" .. " <<
max <<
endl;
135 const word& fieldName,
148 (surfFilePath / surfFileNameBase),
163 const label nFaceZones,
177 includeMap[facei] =
true;
186 / surfFileNameBase +
"_" +
name(
zone) +
".obj"
189 Info<<
"writing part " <<
zone <<
" size " << subSurf.
size()
190 <<
" to " << subName <<
endl;
192 subSurf.
write(subName);
204 for (
const label edgei : markedEdges)
206 edgeSet.insert(edges[edgei]);
211 if (edgeSet.found(edges[edgei]))
213 markedEdges.insert(edgei);
226 forAll(isMarkedEdge, edgei)
228 if (isMarkedEdge[edgei])
236 forAll(isMarkedEdge, edgei)
238 if (isMarkedEdge[edgei])
240 edgeSet.insert(edges[edgei]);
246 if (edgeSet.found(edges[edgei]))
248 isMarkedEdge[edgei] =
true;
265 for (label edgei : edgeSet)
268 edge compactEdge(-1, -1);
271 label& compacti = pointToCompact[
e[ep]];
274 compacti = compactPoints.size();
278 compactEdge[ep] = compacti;
280 compactEdges.append(compactEdge);
283 edgeMesh eMesh(std::move(compactPoints), std::move(compactEdges));
284 eMesh.write(setName);
290int main(
int argc,
char *argv[])
294 "Check geometric and topological quality of a surface"
297 argList::noParallel();
298 argList::addArgument(
"input",
"The input surface file");
300 argList::addBoolOption
302 "checkSelfIntersection",
303 "Also check for self-intersection"
305 argList::addBoolOption
308 "Split surface along non-manifold edges "
309 "(default split is fully disconnected)"
311 argList::addVerboseOption
313 "Additional verbosity"
315 argList::addBoolOption
318 "Write vertices/blocks for blockMeshDict"
324 "Upper limit on the number of files written. "
325 "Default is 10, using 0 suppresses file writing."
331 "Reconstruct and write problem triangles/edges in selected format"
338 const bool checkSelfIntersect =
args.
found(
"checkSelfIntersection");
339 const bool splitNonManifold =
args.
found(
"splitNonManifold");
340 const label outputThreshold =
343 const bool writeSets = !surfaceFormat.empty();
349 surfWriter = surfaceWriter::New(surfaceFormat);
363 Info<<
"Reading surface from " << surfFileName <<
" ..." <<
nl <<
endl;
377 fileName surfFileNameBase(surfFileName.name());
378 const word fileType = surfFileNameBase.
ext();
380 surfFileNameBase = surfFileNameBase.
lessExt();
382 if (fileType ==
"gz")
384 surfFileNameBase = surfFileNameBase.
lessExt();
386 const fileName surfFilePath(surfFileName.path());
394 Info<<
"// blockMeshDict info" <<
nl <<
nl;
396 Info<<
"vertices\n(" <<
nl;
399 Info<<
" " << cornerPts[pti] <<
nl;
406 <<
" hex (0 1 2 3 4 5 6 7) (10 10 10) simpleGrading (1 1 1)\n"
410 <<
"patches\n();" <<
endl;
424 label region = surf[facei].region();
426 if (region < 0 || region >= regionSize.size())
429 <<
"Triangle " << facei <<
" vertices " << surf[facei]
430 <<
" has region " << region <<
" which is outside the range"
436 regionSize[region]++;
440 Info<<
"Region\tSize" <<
nl
441 <<
"------\t----" <<
nl;
445 << regionSize[patchi] <<
nl;
460 if (!triSurfaceTools::validTri(surf, facei,
false))
462 illegalFaces.append(facei);
466 if (illegalFaces.size())
468 Info<<
"Surface has " << illegalFaces.size()
469 <<
" illegal triangles." <<
endl;
481 subSurf.triFaceFaces(faces);
487 (surfFilePath / surfFileNameBase),
499 Info<<
"Wrote illegal triangles to "
502 else if (outputThreshold > 0)
507 triSurfaceTools::validTri(surf, illegalFaces[i],
true);
508 if (i >= outputThreshold)
510 Info<<
"Suppressing further warning messages."
511 <<
" Use -outputThreshold to increase number"
512 <<
" of warnings." <<
endl;
518 Info<<
"Dumping conflicting face labels to " << str.name()
520 <<
"Paste this into the input for surfaceSubset" <<
endl;
526 Info<<
"Surface has no illegal triangles." <<
endl;
542 if (
f[0] ==
f[1] ||
f[0] ==
f[2] ||
f[1] ==
f[2])
559 labelList binCount = countBins(0, 1, 20, triQ);
561 Info<<
"Triangle quality (equilateral=1, collapsed=0):"
568 scalar dist = (1.0 - 0.0)/20.0;
572 Info<<
" " <<
min <<
" .. " <<
min+dist <<
" : "
573 << 1.0/surf.
size() * binCount[bini]
581 const label minIndex = minMaxIds.
first();
582 const label maxIndex = minMaxIds.
second();
584 Info<<
" min " << triQ[minIndex] <<
" for triangle " << minIndex
586 <<
" max " << triQ[maxIndex] <<
" for triangle " << maxIndex
591 if (triQ[minIndex] < SMALL)
594 <<
"Minimum triangle quality is "
595 << triQ[minIndex] <<
". This might give problems in"
596 <<
" self-intersection testing later on." <<
endl;
610 (surfFilePath / surfFileNameBase),
618 Info<<
"Wrote triangle-quality to "
621 else if (outputThreshold > 0)
627 if (triQ[facei] < 1
e-11)
629 problemFaces.append(facei);
633 if (!problemFaces.empty())
637 Info<<
"Dumping bad quality faces to " << str.name() <<
endl
638 <<
"Paste this into the input for surfaceSubset" <<
nl
658 edgeMag[edgei] = edges[edgei].mag(localPoints);
663 const label minEdgei = minMaxIds.
first();
664 const label maxEdgei = minMaxIds.
second();
666 const edge& minE = edges[minEdgei];
667 const edge& maxE = edges[maxEdgei];
671 <<
" min " << edgeMag[minEdgei] <<
" for edge " << minEdgei
672 <<
" points " << localPoints[minE[0]] << localPoints[minE[1]]
674 <<
" max " << edgeMag[maxEdgei] <<
" for edge " << maxEdgei
675 <<
" points " << localPoints[maxE[0]] << localPoints[maxE[1]]
689 scalar smallDim = 1
e-6 * bb.mag();
691 Info<<
"Checking for points less than 1e-6 of bounding box ("
692 << bb.span() <<
" metre) apart."
700 for (label i = 1; i < sortedMag.size(); i++)
702 label pti = sortedMag.indices()[i];
704 label prevPti = sortedMag.indices()[i-1];
706 if (
mag(localPoints[pti] - localPoints[prevPti]) < smallDim)
715 const edge&
e = edges[pEdges[i]];
717 if (
e[0] == prevPti ||
e[1] == prevPti)
726 if (nClose < outputThreshold)
730 Info<<
" close unconnected points "
731 << pti <<
' ' << localPoints[pti]
732 <<
" and " << prevPti <<
' '
733 << localPoints[prevPti]
735 <<
mag(localPoints[pti] - localPoints[prevPti])
740 Info<<
" small edge between points "
741 << pti <<
' ' << localPoints[pti]
742 <<
" and " << prevPti <<
' '
743 << localPoints[prevPti]
745 <<
mag(localPoints[pti] - localPoints[prevPti])
749 else if (nClose == outputThreshold)
751 Info<<
"Suppressing further warning messages."
752 <<
" Use -outputThreshold to increase number"
753 <<
" of warnings." <<
endl;
760 Info<<
"Found " << nClose <<
" nearby points." <<
nl
777 const labelList& myFaces = edgeFaces[edgei];
779 if (myFaces.
size() == 1)
781 problemFaces.
append(myFaces[0]);
782 openEdges.append(edgei);
788 const labelList& myFaces = edgeFaces[edgei];
790 if (myFaces.
size() > 2)
794 problemFaces.append(myFaces[myFacei]);
796 multipleEdges.append(edgei);
799 problemFaces.shrink();
801 if (openEdges.size() || multipleEdges.size())
803 Info<<
"Surface is not closed since not all edges ("
804 << edgeFaces.
size() <<
") connected to "
805 <<
"two faces:" <<
endl
806 <<
" connected to one face : " << openEdges.size() <<
endl
807 <<
" connected to >2 faces : " << multipleEdges.size() <<
endl;
809 Info<<
"Conflicting face labels:" << problemFaces.size() <<
endl;
811 if (!edgeFormat.empty() && openEdges.size())
815 surfFileName.lessExt()
819 Info<<
"Writing open edges to " << openName <<
" ..." <<
endl;
820 writeEdgeSet(openName, surf, openEdges);
822 if (!edgeFormat.empty() && multipleEdges.size())
826 surfFileName.lessExt()
830 Info<<
"Writing multiply connected edges to "
831 << multName <<
" ..." <<
endl;
832 writeEdgeSet(multName, surf, multipleEdges);
837 Info<<
"Surface is closed. All edges connected to two faces." <<
endl;
848 if (splitNonManifold)
852 if (edgeFaces[edgei].size() > 2)
854 borderEdge[edgei] =
true;
857 syncEdges(surf, borderEdge);
863 Info<<
"Number of unconnected parts : " << numZones <<
endl;
865 if (numZones > 1 && outputThreshold > 0)
867 Info<<
"Splitting surface into parts ..." <<
endl <<
endl;
884 if (numZones > outputThreshold)
886 Info<<
"Limiting number of files to " << outputThreshold
892 min(outputThreshold, numZones),
906 PatchTools::checkOrientation(surf,
false, &borderEdge);
917 syncEdges(surf, borderEdge);
923 label numNormalZones = PatchTools::markZones(surf, borderEdge, normalZone);
926 <<
"Number of zones (connected area with consistent normal) : "
927 << numNormalZones <<
endl;
929 if (numNormalZones > 1)
931 Info<<
"More than one normal orientation." <<
endl;
933 if (outputThreshold > 0)
950 if (numNormalZones > outputThreshold)
952 Info<<
"Limiting number of files to " << outputThreshold
958 min(outputThreshold, numNormalZones),
961 surfFileNameBase +
"_normal"
972 if (checkSelfIntersect)
974 Info<<
"Checking self-intersection." <<
endl;
981 if (outputThreshold > 0)
1005 intStreamPtr().write(hitInfo.hitPoint());
1011 if (hitInfo2.hit() && hitInfo.index() != hitInfo2.index())
1017 intStreamPtr().write(hitInfo2.hitPoint());
1065 Info<<
"Surface is not self-intersecting" <<
endl;
1069 Info<<
"Surface is self-intersecting at " << nInt
1070 <<
" locations." <<
endl;
1074 Info<<
"Writing intersection points to "
1075 << intStreamPtr().name() <<
endl;
vtk::internalMeshWriter writer(topoMesh, topoCells, vtk::formatType::INLINE_ASCII, runTime.path()/"blockTopology")
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
T & first() noexcept
The first element of the list, position [0].
label size() const noexcept
The number of elements in table.
void append(const T &val)
Append an element at the end of the list.
OFstream that keeps track of vertices.
Output to file stream, using an OSstream.
Generic output stream using a standard (STL) stream.
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
const T & second() const noexcept
Return second element, which is also the last element.
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
label nEdges() const
Number of edges in patch.
label nPoints() const
Number of points supporting patch faces.
const labelListList & pointEdges() const
Return point-edge addressing.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
const Field< point_type > & localPoints() const
Return pointField of points in patch.
const Field< point_type > & points() const noexcept
Return reference to global points.
const labelListList & edgeFaces() const
Return edge-face addressing.
A list that is sorted upon construction or when explicitly requested with the sort() method.
A List with indirect addressing. Like IndirectList but does not store addressing.
void size(const label n)
Older name for setAddressableSize.
Extract command arguments and options from the supplied argc and argv parameters.
T get(const label index) const
Get a value from the argument at index.
bool found(const word &optName) const
Return true if the named option is found.
T getOrDefault(const word &optName, const T &deflt) const
Get a value from the named option if present, or return default.
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
void clear() noexcept
Same as reset(nullptr)
void reset(autoPtr< T > &&other) noexcept
Delete managed object and set to new given pointer.
A bounding box defined in terms of min/max extrema points.
tmp< pointField > points() const
Corner points in an order corresponding to a 'hex' cell.
Mesh data needed to do the Finite Area discretisation.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
A subset of mesh faces organised as a primitive patch.
A class for handling file names.
fileName lessExt() const
Return file name without extension (part before last .)
word ext() const
Return file name extension (part after last .)
Non-pointer based hierarchical recursive searching.
pointIndexHit findLineAny(const point &start, const point &end) const
Find any intersection of line between start and end.
A triFace with additional (region) index.
Base class for surface writers.
A surfaceWriter for VTK legacy (.vtk) or XML (.vtp) format.
Helper class to search on triSurface.
Triangulated surface description with patch information.
label markZones(const boolList &borderEdge, labelList &faceZone) const
(size and) fills faceZone with zone of face. Zone is area
void triFaceFaces(List< face > &plainFaceList) const
Create a list of faces from the triFaces.
triSurface subsetMesh(const UList< bool > &include, labelList &pointMap, labelList &faceMap) const
Return a new surface subsetted on the selected faces.
const geometricSurfacePatchList & patches() const noexcept
void writeStats(Ostream &os) const
Write some statistics.
A triangle primitive used to calculate face normals and swept volumes.
scalar quality() const
Return quality: Ratio of triangle and circum-circle.
A class for handling words, derived from Foam::string.
Base class for mesh zones.
word outputName("finiteArea-edges.obj")
OBJstream os(runTime.globalPath()/outputName)
#define WarningInFunction
Report a warning using Foam::Warning.
messageStream Info
Information stream (stdout output on master, null elsewhere)
labelPair findMinMax(const ListType &input, label start=0)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
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.
Foam::argList args(argc, argv)
#define forAll(list, i)
Loop across all elements in list.