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,
137 const word& surfFileStem
148 (surfFilePath / surfFileStem),
163 const label nFaceZones,
166 const word& surfFileStem
177 includeMap[facei] =
true;
186 / surfFileStem +
"_" +
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);
290 int main(
int argc,
char *argv[])
294 "Check geometric and topological quality of a surface" 302 "checkSelfIntersection",
303 "Also check for self-intersection" 308 "Split surface along non-manifold edges " 309 "(default split is fully disconnected)" 313 "Additional verbosity" 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();
363 Info<<
"Reading surface from " 378 const fileName surfFilePath(surfName.path());
379 word surfFileStem(surfName.stem());
382 if (surfName.has_ext(
"gz"))
393 Info<<
"// blockMeshDict info" <<
nl <<
nl;
395 Info<<
"vertices\n(" <<
nl;
398 Info<<
" " << cornerPts[pti] <<
nl;
405 <<
" hex (0 1 2 3 4 5 6 7) (10 10 10) simpleGrading (1 1 1)\n" 409 <<
"patches\n();" <<
endl;
423 label region = surf[facei].region();
425 if (region < 0 || region >= regionSize.size())
428 <<
"Triangle " << facei <<
" vertices " << surf[facei]
429 <<
" has region " << region <<
" which is outside the range" 435 regionSize[region]++;
439 Info<<
"Region\tSize" <<
nl 440 <<
"------\t----" <<
nl;
444 << regionSize[patchi] <<
nl;
461 illegalFaces.append(facei);
465 if (illegalFaces.size())
467 Info<<
"Surface has " << illegalFaces.size()
468 <<
" illegal triangles." <<
endl;
480 subSurf.triFaceFaces(faces);
486 (surfFilePath / surfFileStem),
498 Info<<
"Wrote illegal triangles to " 501 else if (outputThreshold > 0)
507 if (i >= outputThreshold)
509 Info<<
"Suppressing further warning messages." 510 <<
" Use -outputThreshold to increase number" 511 <<
" of warnings." <<
endl;
517 Info<<
"Dumping conflicting face labels to " << str.name()
519 <<
"Paste this into the input for surfaceSubset" <<
endl;
525 Info<<
"Surface has no illegal triangles." <<
endl;
541 if (
f[0] ==
f[1] ||
f[0] ==
f[2] ||
f[1] ==
f[2])
558 labelList binCount = countBins(0, 1, 20, triQ);
560 Info<<
"Triangle quality (equilateral=1, collapsed=0):" 567 scalar dist = (1.0 - 0.0)/20.0;
571 Info<<
" " <<
min <<
" .. " <<
min+dist <<
" : " 572 << 1.0/surf.
size() * binCount[bini]
580 const label minIndex = minMaxIds.
first();
581 const label maxIndex = minMaxIds.
second();
583 Info<<
" min " << triQ[minIndex] <<
" for triangle " << minIndex
585 <<
" max " << triQ[maxIndex] <<
" for triangle " << maxIndex
590 if (triQ[minIndex] < SMALL)
593 <<
"Minimum triangle quality is " 594 << triQ[minIndex] <<
". This might give problems in" 595 <<
" self-intersection testing later on." <<
endl;
609 (surfFilePath / surfFileStem),
617 Info<<
"Wrote triangle-quality to " 620 else if (outputThreshold > 0)
626 if (triQ[facei] < 1
e-11)
628 problemFaces.append(facei);
632 if (!problemFaces.empty())
636 Info<<
"Dumping bad quality faces to " << str.name() <<
endl 637 <<
"Paste this into the input for surfaceSubset" <<
nl 657 edgeMag[edgei] = edges[edgei].mag(localPoints);
662 const label minEdgei = minMaxIds.
first();
663 const label maxEdgei = minMaxIds.
second();
665 const edge& minE = edges[minEdgei];
666 const edge& maxE = edges[maxEdgei];
670 <<
" min " << edgeMag[minEdgei] <<
" for edge " << minEdgei
671 <<
" points " << localPoints[minE[0]] << localPoints[minE[1]]
673 <<
" max " << edgeMag[maxEdgei] <<
" for edge " << maxEdgei
674 <<
" points " << localPoints[maxE[0]] << localPoints[maxE[1]]
688 scalar smallDim = 1
e-6 * bb.mag();
690 Info<<
"Checking for points less than 1e-6 of bounding box (" 691 << bb.span() <<
" metre) apart." 699 for (label i = 1; i < sortedMag.size(); i++)
701 label pti = sortedMag.indices()[i];
703 label prevPti = sortedMag.indices()[i-1];
705 if (
mag(localPoints[pti] - localPoints[prevPti]) < smallDim)
714 const edge&
e = edges[pEdges[i]];
716 if (
e[0] == prevPti ||
e[1] == prevPti)
725 if (nClose < outputThreshold)
729 Info<<
" close unconnected points " 730 << pti <<
' ' << localPoints[pti]
731 <<
" and " << prevPti <<
' ' 732 << localPoints[prevPti]
734 <<
mag(localPoints[pti] - localPoints[prevPti])
739 Info<<
" small edge between points " 740 << pti <<
' ' << localPoints[pti]
741 <<
" and " << prevPti <<
' ' 742 << localPoints[prevPti]
744 <<
mag(localPoints[pti] - localPoints[prevPti])
748 else if (nClose == outputThreshold)
750 Info<<
"Suppressing further warning messages." 751 <<
" Use -outputThreshold to increase number" 752 <<
" of warnings." <<
endl;
759 Info<<
"Found " << nClose <<
" nearby points." <<
nl 776 const labelList& myFaces = edgeFaces[edgei];
778 if (myFaces.
size() == 1)
780 problemFaces.
append(myFaces[0]);
781 openEdges.append(edgei);
787 const labelList& myFaces = edgeFaces[edgei];
789 if (myFaces.
size() > 2)
793 problemFaces.append(myFaces[myFacei]);
795 multipleEdges.append(edgei);
798 problemFaces.shrink();
800 if (openEdges.size() || multipleEdges.size())
802 Info<<
"Surface is not closed since not all edges (" 803 << edgeFaces.
size() <<
") connected to " 804 <<
"two faces:" <<
endl 805 <<
" connected to one face : " << openEdges.size() <<
endl 806 <<
" connected to >2 faces : " << multipleEdges.size() <<
endl;
808 Info<<
"Conflicting face labels:" << problemFaces.size() <<
endl;
810 if (!edgeFormat.empty() && openEdges.size())
818 Info<<
"Writing open edges to " 822 if (!edgeFormat.empty() && multipleEdges.size())
830 Info<<
"Writing multiply connected edges to " 832 writeEdgeSet(
outputName, 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),
917 syncEdges(surf, borderEdge);
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 surfFileStem +
"_normal" 972 if (checkSelfIntersect)
974 Info<<
"Checking self-intersection." <<
endl;
981 if (outputThreshold > 0)
1005 intStreamPtr().write(hitInfo.point());
1011 if (hitInfo2.hit() && hitInfo.index() != hitInfo2.index())
1017 intStreamPtr().write(hitInfo2.point());
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;
const labelListList & pointEdges() const
Return point-edge addressing.
const T & first() const noexcept
Access the first element.
label nPoints() const
Number of points supporting patch faces.
triSurface subsetMesh(const UList< bool > &include, labelList &pointMap, labelList &faceMap) const
Return a new surface subsetted on the selected faces.
Generic output stream using a standard (STL) stream.
Field< label > labelField
Specialisation of Field<T> for label.
static void addNote(const string ¬e)
Add extra notes for the usage information.
void size(const label n)
Older name for setAddressableSize.
labelPair findMinMax(const ListType &input, label start=0)
Linear search for the index of the min/max element, similar to std::minmax_element but for lists and ...
A class for handling file names.
A surfaceWriter for VTK legacy (.vtk) or XML (.vtp) format.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
const Field< point_type > & localPoints() const
Return pointField of points in patch.
void writeStats(Ostream &os) const
Write some statistics.
void append(const T &val)
Append an element at the end of the list.
A list that is sorted upon construction or when explicitly requested with the sort() method...
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
virtual void open(const fileName &outputPath)
Open for output on specified path, using existing surface.
Output to file stream, using an OSstream.
constexpr char nl
The newline '\n' character (0x0a)
Ostream & endl(Ostream &os)
Add newline and flush stream.
static void addBoolOption(const word &optName, const string &usage="", bool advanced=false)
Add a bool option to validOptions with usage information.
static void noParallel()
Remove the parallel options.
A bounding box defined in terms of min/max extrema points.
label size() const noexcept
The number of elements in table.
This class describes the interaction of an object (often a face) and a point. It carries the info of ...
T getOrDefault(const word &optName, const T &deflt) const
Get a value from the named option if present, or return default.
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
#define forAll(list, i)
Loop across all elements in list.
void reset(T *p=nullptr) noexcept
Delete managed object and set to new given pointer.
virtual void clear()
Close any open output, remove association with a surface and expire the writer. The parallel flag rem...
Helper class to search on triSurface.
const geometricSurfacePatchList & patches() const noexcept
word outputName("finiteArea-edges.obj")
const dimensionedScalar e
Elementary charge.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Base class for mesh zones.
A class for handling words, derived from Foam::string.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Extract command arguments and options from the supplied argc and argv parameters. ...
Tree tree(triangles.begin(), triangles.end())
static void addVerboseOption(const string &usage, bool advanced=false)
Enable a 'verbose' bool option, with usage information.
static void addOption(const word &optName, const string ¶m="", const string &usage="", bool advanced=false)
Add an option to validOptions with usage information.
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
const labelListList & edgeFaces() const
Return edge-face addressing.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
A triFace with additional (region) index.
const Field< point_type > & points() const noexcept
Return reference to global points.
void append(const T &val)
Copy append an element to the end of this list.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
fileName relativePath(const fileName &input, const bool caseTag=false) const
Return the input relative to the globalPath by stripping off a leading value of the globalPath...
label nEdges() const
Number of edges in patch.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
An OFstream that keeps track of vertices and provides convenience output methods for OBJ files...
vtk::internalMeshWriter writer(topoMesh, topoCells, vtk::formatType::INLINE_ASCII, runTime.path()/"blockTopology")
void triFaceFaces(List< face > &plainFaceList) const
Create a list of faces from the triFaces.
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
OBJstream os(runTime.globalPath()/outputName)
virtual bool open(const fileName &file, bool parallel=UPstream::parRun())
Open file for writing (creates parent directory).
Mesh data needed to do the Finite Area discretisation.
virtual int width() const
Get width of output field.
T get(const label index) const
Get a value from the argument at index.
Non-pointer based hierarchical recursive searching.
#define WarningInFunction
Report a warning using Foam::Warning.
A List with indirect addressing. Like IndirectList but does not store addressing. ...
static void addArgument(const string &argName, const string &usage="")
Append a (mandatory) argument to validArgs.
triangle< point, const point & > triPointRef
A triangle using referred points.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Base class for surface writers.
Pointer management similar to std::unique_ptr, with some additional methods and type checking...
virtual fileName write()=0
Write separate surface geometry to file.
A subset of mesh faces organised as a primitive patch.
label markZones(const boolList &borderEdge, labelList &faceZone) const
(size and) fills faceZone with zone of face. Zone is area
const T & second() const noexcept
Access the second element.
Triangulated surface description with patch information.
Foam::argList args(argc, argv)
static autoPtr< surfaceWriter > New(const word &writeType)
Return a reference to the selected surfaceWriter.
bool found(const word &optName) const
Return true if the named option is found.
bool remove_ext()
Remove extension, return true if string changed.
static constexpr const zero Zero
Global zero (0)