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)" 315 "Write vertices/blocks for blockMeshDict" 321 "Upper limit on the number of files written. " 322 "Default is 10, using 0 suppresses file writing." 328 "Reconstruct and write problem triangles/edges in selected format" 335 const bool checkSelfIntersect =
args.
found(
"checkSelfIntersection");
336 const bool splitNonManifold =
args.
found(
"splitNonManifold");
337 const label outputThreshold =
340 const bool writeSets = !surfaceFormat.empty();
360 Info<<
"Reading surface from " 375 const fileName surfFilePath(surfName.path());
376 word surfFileStem(surfName.stem());
379 if (surfName.has_ext(
"gz"))
390 Info<<
"// blockMeshDict info" <<
nl <<
nl;
392 Info<<
"vertices\n(" <<
nl;
395 Info<<
" " << cornerPts[pti] <<
nl;
402 <<
" hex (0 1 2 3 4 5 6 7) (10 10 10) simpleGrading (1 1 1)\n" 406 <<
"patches\n();" <<
endl;
420 label region = surf[facei].region();
422 if (region < 0 || region >= regionSize.size())
425 <<
"Triangle " << facei <<
" vertices " << surf[facei]
426 <<
" has region " << region <<
" which is outside the range" 432 regionSize[region]++;
436 Info<<
"Region\tSize" <<
nl 437 <<
"------\t----" <<
nl;
441 << regionSize[patchi] <<
nl;
458 illegalFaces.append(facei);
462 if (illegalFaces.size())
464 Info<<
"Surface has " << illegalFaces.size()
465 <<
" illegal triangles." <<
endl;
477 subSurf.triFaceFaces(faces);
483 (surfFilePath / surfFileStem),
495 Info<<
"Wrote illegal triangles to " 498 else if (outputThreshold > 0)
504 if (i >= outputThreshold)
506 Info<<
"Suppressing further warning messages." 507 <<
" Use -outputThreshold to increase number" 508 <<
" of warnings." <<
endl;
514 Info<<
"Dumping conflicting face labels to " << str.name()
516 <<
"Paste this into the input for surfaceSubset" <<
endl;
522 Info<<
"Surface has no illegal triangles." <<
endl;
538 if (
f[0] ==
f[1] ||
f[0] ==
f[2] ||
f[1] ==
f[2])
555 labelList binCount = countBins(0, 1, 20, triQ);
557 Info<<
"Triangle quality (equilateral=1, collapsed=0):" 564 scalar dist = (1.0 - 0.0)/20.0;
568 Info<<
" " <<
min <<
" .. " <<
min+dist <<
" : " 569 << 1.0/surf.
size() * binCount[bini]
577 const label minIndex = minMaxIds.
first();
578 const label maxIndex = minMaxIds.
second();
580 Info<<
" min " << triQ[minIndex] <<
" for triangle " << minIndex
582 <<
" max " << triQ[maxIndex] <<
" for triangle " << maxIndex
587 if (triQ[minIndex] < SMALL)
590 <<
"Minimum triangle quality is " 591 << triQ[minIndex] <<
". This might give problems in" 592 <<
" self-intersection testing later on." <<
endl;
606 (surfFilePath / surfFileStem),
614 Info<<
"Wrote triangle-quality to " 617 else if (outputThreshold > 0)
623 if (triQ[facei] < 1
e-11)
625 problemFaces.append(facei);
629 if (!problemFaces.empty())
633 Info<<
"Dumping bad quality faces to " << str.name() <<
endl 634 <<
"Paste this into the input for surfaceSubset" <<
nl 654 edgeMag[edgei] = edges[edgei].mag(localPoints);
659 const label minEdgei = minMaxIds.
first();
660 const label maxEdgei = minMaxIds.
second();
662 const edge& minE = edges[minEdgei];
663 const edge& maxE = edges[maxEdgei];
667 <<
" min " << edgeMag[minEdgei] <<
" for edge " << minEdgei
668 <<
" points " << localPoints[minE[0]] << localPoints[minE[1]]
670 <<
" max " << edgeMag[maxEdgei] <<
" for edge " << maxEdgei
671 <<
" points " << localPoints[maxE[0]] << localPoints[maxE[1]]
685 scalar smallDim = 1
e-6 * bb.mag();
687 Info<<
"Checking for points less than 1e-6 of bounding box (" 688 << bb.span() <<
" metre) apart." 696 for (label i = 1; i < sortedMag.size(); i++)
698 label pti = sortedMag.indices()[i];
700 label prevPti = sortedMag.indices()[i-1];
702 if (
mag(localPoints[pti] - localPoints[prevPti]) < smallDim)
711 const edge&
e = edges[pEdges[i]];
713 if (
e[0] == prevPti ||
e[1] == prevPti)
722 if (nClose < outputThreshold)
726 Info<<
" close unconnected points " 727 << pti <<
' ' << localPoints[pti]
728 <<
" and " << prevPti <<
' ' 729 << localPoints[prevPti]
731 <<
mag(localPoints[pti] - localPoints[prevPti])
736 Info<<
" small edge between points " 737 << pti <<
' ' << localPoints[pti]
738 <<
" and " << prevPti <<
' ' 739 << localPoints[prevPti]
741 <<
mag(localPoints[pti] - localPoints[prevPti])
745 else if (nClose == outputThreshold)
747 Info<<
"Suppressing further warning messages." 748 <<
" Use -outputThreshold to increase number" 749 <<
" of warnings." <<
endl;
756 Info<<
"Found " << nClose <<
" nearby points." <<
nl 773 const labelList& myFaces = edgeFaces[edgei];
775 if (myFaces.
size() == 1)
777 problemFaces.
append(myFaces[0]);
778 openEdges.append(edgei);
784 const labelList& myFaces = edgeFaces[edgei];
786 if (myFaces.
size() > 2)
790 problemFaces.append(myFaces[myFacei]);
792 multipleEdges.append(edgei);
795 problemFaces.shrink();
797 if (openEdges.size() || multipleEdges.size())
799 Info<<
"Surface is not closed since not all edges (" 800 << edgeFaces.
size() <<
") connected to " 801 <<
"two faces:" <<
endl 802 <<
" connected to one face : " << openEdges.size() <<
endl 803 <<
" connected to >2 faces : " << multipleEdges.size() <<
endl;
805 Info<<
"Conflicting face labels:" << problemFaces.size() <<
endl;
807 if (!edgeFormat.empty() && openEdges.size())
815 Info<<
"Writing open edges to " 819 if (!edgeFormat.empty() && multipleEdges.size())
827 Info<<
"Writing multiply connected edges to " 829 writeEdgeSet(
outputName, surf, multipleEdges);
834 Info<<
"Surface is closed. All edges connected to two faces." <<
endl;
845 if (splitNonManifold)
849 if (edgeFaces[edgei].size() > 2)
851 borderEdge[edgei] =
true;
854 syncEdges(surf, borderEdge);
860 Info<<
"Number of unconnected parts : " << numZones <<
endl;
862 if (numZones > 1 && outputThreshold > 0)
864 Info<<
"Splitting surface into parts ..." <<
endl <<
endl;
881 if (numZones > outputThreshold)
883 Info<<
"Limiting number of files to " << outputThreshold
889 min(outputThreshold, numZones),
914 syncEdges(surf, borderEdge);
923 <<
"Number of zones (connected area with consistent normal) : " 924 << numNormalZones <<
endl;
926 if (numNormalZones > 1)
928 Info<<
"More than one normal orientation." <<
endl;
930 if (outputThreshold > 0)
947 if (numNormalZones > outputThreshold)
949 Info<<
"Limiting number of files to " << outputThreshold
955 min(outputThreshold, numNormalZones),
958 surfFileStem +
"_normal" 969 if (checkSelfIntersect)
971 Info<<
"Checking self-intersection." <<
endl;
978 if (outputThreshold > 0)
1002 intStreamPtr().write(hitInfo.point());
1008 if (hitInfo2.hit() && hitInfo.index() != hitInfo2.index())
1014 intStreamPtr().write(hitInfo2.point());
1062 Info<<
"Surface is not self-intersecting" <<
endl;
1066 Info<<
"Surface is self-intersecting at " << nInt
1067 <<
" locations." <<
endl;
1071 Info<<
"Writing intersection points to " 1072 << 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.
vtk::lineWriter writer(edgeCentres, edgeList::null(), fileName(aMesh.time().globalPath()/"finiteArea-edgesCentres"))
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.
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.
void write(const word &fieldName, const UList< Type > &field)
Write primitive field of CellData (Poly or Line) or PointData values.
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...
label size() const noexcept
The number of elements in table.
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 vertex labels. This can correspond to a directed graph edge or an edge 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 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...
static void addVerboseOption(const string &usage="", bool advanced=false)
Enable a 'verbose' bool option, with usage information.
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...
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)