surfaceCheck.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9 Copyright (C) 2016-2021 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27Application
28 surfaceCheck
29
30Group
31 grpSurfaceUtilities
32
33Description
34 Check geometric and topological quality of a surface.
35
36Usage
37 \b surfaceCheck [OPTION] surfaceFile
38
39 Options:
40 - \par -checkSelfIntersection
41 Check for self-intersection.
42
43 - \par -splitNonManifold
44 Split surface along non-manifold edges.
45
46 - \par -verbose
47 Extra verbosity.
48
49 - \par -blockMesh
50 Write vertices/blocks for tight-fitting 1 cell blockMeshDict.
51
52 - \par -outputThreshold <num files>
53 Upper limit on the number of files written.
54 Prevent surfaces with many disconnected parts from writing many files.
55 The default is 10. A value of 0 suppresses file writing.
56
57\*---------------------------------------------------------------------------*/
58
59#include "triangle.H"
60#include "edgeHashes.H"
61#include "triSurface.H"
62#include "triSurfaceSearch.H"
63#include "triSurfaceTools.H"
64#include "argList.H"
65#include "OFstream.H"
66#include "OBJstream.H"
67#include "SortableList.H"
68#include "PatchTools.H"
69#include "vtkSurfaceWriter.H"
70#include "functionObject.H"
71#include "DynamicField.H"
72#include "edgeMesh.H"
73
74using namespace Foam;
75
76labelList countBins
77(
78 const scalar min,
79 const scalar max,
80 const label nBins,
81 const scalarField& vals
82)
83{
84 scalar dist = nBins/(max - min);
85
86 labelList binCount(nBins, Zero);
87
88 forAll(vals, i)
89 {
90 scalar val = vals[i];
91
92 label index = -1;
93
94 if (Foam::mag(val - min) < SMALL)
95 {
96 index = 0;
97 }
98 else if (val >= max - SMALL)
99 {
100 index = nBins - 1;
101 }
102 else
103 {
104 index = label((val - min)*dist);
105
106 if ((index < 0) || (index >= nBins))
107 {
109 << "value " << val << " at index " << i
110 << " outside range " << min << " .. " << max << endl;
111
112 if (index < 0)
113 {
114 index = 0;
115 }
116 else
117 {
118 index = nBins - 1;
119 }
120 }
121 }
122 binCount[index]++;
123 }
124
125 return binCount;
126}
127
128
129
130void writeZoning
131(
133 const triSurface& surf,
134 const labelList& faceZone,
135 const word& fieldName,
136 const fileName& surfFilePath,
137 const fileName& surfFileNameBase
138)
139{
140 // Transcribe faces
141 faceList faces;
142 surf.triFaceFaces(faces);
143
144 writer.open
145 (
146 surf.points(),
147 faces,
148 (surfFilePath / surfFileNameBase),
149 false // serial - already merged
150 );
151
152 fileName outputName = writer.write(fieldName, labelField(faceZone));
153
154 writer.clear();
155
156 Info<< "Wrote zoning to " << outputName << nl << endl;
157}
158
159
160void writeParts
161(
162 const triSurface& surf,
163 const label nFaceZones,
164 const labelList& faceZone,
165 const fileName& surfFilePath,
166 const fileName& surfFileNameBase
167)
168{
169 for (label zone = 0; zone < nFaceZones; zone++)
170 {
171 boolList includeMap(surf.size(), false);
172
173 forAll(faceZone, facei)
174 {
175 if (faceZone[facei] == zone)
176 {
177 includeMap[facei] = true;
178 }
179 }
180
181 triSurface subSurf(surf.subsetMesh(includeMap));
182
183 fileName subName
184 (
185 surfFilePath
186 / surfFileNameBase + "_" + name(zone) + ".obj"
187 );
188
189 Info<< "writing part " << zone << " size " << subSurf.size()
190 << " to " << subName << endl;
191
192 subSurf.write(subName);
193 }
194}
195
196
197void syncEdges(const triSurface& p, labelHashSet& markedEdges)
198{
199 // See comment below about having duplicate edges
200
201 const edgeList& edges = p.edges();
202 edgeHashSet edgeSet(2*markedEdges.size());
203
204 for (const label edgei : markedEdges)
205 {
206 edgeSet.insert(edges[edgei]);
207 }
208
209 forAll(edges, edgei)
210 {
211 if (edgeSet.found(edges[edgei]))
212 {
213 markedEdges.insert(edgei);
214 }
215 }
216}
217
218
219void syncEdges(const triSurface& p, boolList& isMarkedEdge)
220{
221 // See comment below about having duplicate edges
222
223 const edgeList& edges = p.edges();
224
225 label n = 0;
226 forAll(isMarkedEdge, edgei)
227 {
228 if (isMarkedEdge[edgei])
229 {
230 n++;
231 }
232 }
233
234 edgeHashSet edgeSet(2*n);
235
236 forAll(isMarkedEdge, edgei)
237 {
238 if (isMarkedEdge[edgei])
239 {
240 edgeSet.insert(edges[edgei]);
241 }
242 }
243
244 forAll(edges, edgei)
245 {
246 if (edgeSet.found(edges[edgei]))
247 {
248 isMarkedEdge[edgei] = true;
249 }
250 }
251}
252
253
254void writeEdgeSet
255(
256 const word& setName,
257 const triSurface& surf,
258 const labelUList& edgeSet
259)
260{
261 // Get compact edge mesh
262 labelList pointToCompact(surf.nPoints(), -1);
263 DynamicField<point> compactPoints(edgeSet.size());
264 DynamicList<edge> compactEdges(edgeSet.size());
265 for (label edgei : edgeSet)
266 {
267 const edge& e = surf.edges()[edgei];
268 edge compactEdge(-1, -1);
269 forAll(e, ep)
270 {
271 label& compacti = pointToCompact[e[ep]];
272 if (compacti == -1)
273 {
274 compacti = compactPoints.size();
275 label pointi = surf.meshPoints()[e[ep]];
276 compactPoints.append(surf.points()[pointi]);
277 }
278 compactEdge[ep] = compacti;
279 }
280 compactEdges.append(compactEdge);
281 }
282
283 edgeMesh eMesh(std::move(compactPoints), std::move(compactEdges));
284 eMesh.write(setName);
285}
286
287
288// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
289
290int main(int argc, char *argv[])
291{
292 argList::addNote
293 (
294 "Check geometric and topological quality of a surface"
295 );
296
297 argList::noParallel();
298 argList::addArgument("input", "The input surface file");
299
300 argList::addBoolOption
301 (
302 "checkSelfIntersection",
303 "Also check for self-intersection"
304 );
305 argList::addBoolOption
306 (
307 "splitNonManifold",
308 "Split surface along non-manifold edges "
309 "(default split is fully disconnected)"
310 );
311 argList::addVerboseOption
312 (
313 "Additional verbosity"
314 );
315 argList::addBoolOption
316 (
317 "blockMesh",
318 "Write vertices/blocks for blockMeshDict"
319 );
320 argList::addOption
321 (
322 "outputThreshold",
323 "number",
324 "Upper limit on the number of files written. "
325 "Default is 10, using 0 suppresses file writing."
326 );
327 argList::addOption
328 (
329 "writeSets",
330 "surfaceFormat",
331 "Reconstruct and write problem triangles/edges in selected format"
332 );
333
334
335 argList args(argc, argv);
336
337 const auto surfFileName = args.get<fileName>(1);
338 const bool checkSelfIntersect = args.found("checkSelfIntersection");
339 const bool splitNonManifold = args.found("splitNonManifold");
340 const label outputThreshold =
341 args.getOrDefault<label>("outputThreshold", 10);
342 const word surfaceFormat = args.getOrDefault<word>("writeSets", "");
343 const bool writeSets = !surfaceFormat.empty();
344
345 autoPtr<surfaceWriter> surfWriter;
346 word edgeFormat;
347 if (writeSets)
348 {
349 surfWriter = surfaceWriter::New(surfaceFormat);
350
351 // Option1: hard-coded format
352 edgeFormat = "obj";
355 //edgeFormat = surfaceFormat;
356 //if (!edgeMesh::canWriteType(edgeFormat))
357 //{
358 // edgeFormat = "obj";
359 //}
360 }
361
362
363 Info<< "Reading surface from " << surfFileName << " ..." << nl << endl;
364
365 // Read
366 // ~~~~
367
368 triSurface surf(surfFileName);
369
370
371 Info<< "Statistics:" << endl;
372 surf.writeStats(Info);
373 Info<< endl;
374
375
376 // Determine path and extension
377 fileName surfFileNameBase(surfFileName.name());
378 const word fileType = surfFileNameBase.ext();
379 // Strip extension
380 surfFileNameBase = surfFileNameBase.lessExt();
381 // If extension was .gz strip original extension
382 if (fileType == "gz")
383 {
384 surfFileNameBase = surfFileNameBase.lessExt();
385 }
386 const fileName surfFilePath(surfFileName.path());
387
388
389 // write bounding box corners
390 if (args.found("blockMesh"))
391 {
392 pointField cornerPts(boundBox(surf.points(), false).points());
393
394 Info<< "// blockMeshDict info" << nl << nl;
395
396 Info<< "vertices\n(" << nl;
397 forAll(cornerPts, pti)
398 {
399 Info<< " " << cornerPts[pti] << nl;
400 }
401
402 // number of divisions needs adjustment later
403 Info<< ");\n" << nl
404 << "blocks\n"
405 << "(\n"
406 << " hex (0 1 2 3 4 5 6 7) (10 10 10) simpleGrading (1 1 1)\n"
407 << ");\n" << nl;
408
409 Info<< "edges\n();" << nl
410 << "patches\n();" << endl;
411
412 Info<< nl << "// end blockMeshDict info" << nl << endl;
413 }
414
415
416 // Region sizes
417 // ~~~~~~~~~~~~
418
419 {
420 labelList regionSize(surf.patches().size(), Zero);
421
422 forAll(surf, facei)
423 {
424 label region = surf[facei].region();
425
426 if (region < 0 || region >= regionSize.size())
427 {
429 << "Triangle " << facei << " vertices " << surf[facei]
430 << " has region " << region << " which is outside the range"
431 << " of regions 0.." << surf.patches().size()-1
432 << endl;
433 }
434 else
435 {
436 regionSize[region]++;
437 }
438 }
439
440 Info<< "Region\tSize" << nl
441 << "------\t----" << nl;
442 forAll(surf.patches(), patchi)
443 {
444 Info<< surf.patches()[patchi].name() << '\t'
445 << regionSize[patchi] << nl;
446 }
447 Info<< nl << endl;
448 }
449
450
451 // Check triangles
452 // ~~~~~~~~~~~~~~~
453
454 {
455 DynamicList<label> illegalFaces(surf.size()/100 + 1);
456
457 forAll(surf, facei)
458 {
459 // Check silently
460 if (!triSurfaceTools::validTri(surf, facei, false))
461 {
462 illegalFaces.append(facei);
463 }
464 }
465
466 if (illegalFaces.size())
467 {
468 Info<< "Surface has " << illegalFaces.size()
469 << " illegal triangles." << endl;
470
471 if (surfWriter)
472 {
473 boolList isIllegalFace(surf.size(), false);
474 UIndirectList<bool>(isIllegalFace, illegalFaces) = true;
475
476 triSurface subSurf(surf.subsetMesh(isIllegalFace));
477
478
479 // Transcribe faces
480 faceList faces;
481 subSurf.triFaceFaces(faces);
482
483 surfWriter->open
484 (
485 subSurf.points(),
486 faces,
487 (surfFilePath / surfFileNameBase),
488 false // serial - already merged
489 );
490
491 fileName outputName = surfWriter->write
492 (
493 "illegal",
494 scalarField(subSurf.size(), Zero)
495 );
496
497 surfWriter->clear();
498
499 Info<< "Wrote illegal triangles to "
500 << outputName << nl << endl;
501 }
502 else if (outputThreshold > 0)
503 {
504 forAll(illegalFaces, i)
505 {
506 // Force warning message
507 triSurfaceTools::validTri(surf, illegalFaces[i], true);
508 if (i >= outputThreshold)
509 {
510 Info<< "Suppressing further warning messages."
511 << " Use -outputThreshold to increase number"
512 << " of warnings." << endl;
513 break;
514 }
515 }
516
517 OFstream str("illegalFaces");
518 Info<< "Dumping conflicting face labels to " << str.name()
519 << endl
520 << "Paste this into the input for surfaceSubset" << endl;
521 str << illegalFaces;
522 }
523 }
524 else
525 {
526 Info<< "Surface has no illegal triangles." << endl;
527 }
528 Info<< endl;
529 }
530
531
532
533 // Triangle quality
534 // ~~~~~~~~~~~~~~~~
535
536 {
537 scalarField triQ(surf.size(), Zero);
538 forAll(surf, facei)
539 {
540 const labelledTri& f = surf[facei];
541
542 if (f[0] == f[1] || f[0] == f[2] || f[1] == f[2])
543 {
544 //WarningInFunction
545 // << "Illegal triangle " << facei << " vertices " << f
546 // << " coords " << f.points(surf.points()) << endl;
547 }
548 else
549 {
550 triQ[facei] = triPointRef
551 (
552 surf.points()[f[0]],
553 surf.points()[f[1]],
554 surf.points()[f[2]]
555 ).quality();
556 }
557 }
558
559 labelList binCount = countBins(0, 1, 20, triQ);
560
561 Info<< "Triangle quality (equilateral=1, collapsed=0):"
562 << endl;
563
564
565 OSstream& os = Info;
566 os.width(4);
567
568 scalar dist = (1.0 - 0.0)/20.0;
569 scalar min = 0;
570 forAll(binCount, bini)
571 {
572 Info<< " " << min << " .. " << min+dist << " : "
573 << 1.0/surf.size() * binCount[bini]
574 << endl;
575 min += dist;
576 }
577 Info<< endl;
578
579 labelPair minMaxIds = findMinMax(triQ);
580
581 const label minIndex = minMaxIds.first();
582 const label maxIndex = minMaxIds.second();
583
584 Info<< " min " << triQ[minIndex] << " for triangle " << minIndex
585 << nl
586 << " max " << triQ[maxIndex] << " for triangle " << maxIndex
587 << nl
588 << endl;
589
590
591 if (triQ[minIndex] < SMALL)
592 {
594 << "Minimum triangle quality is "
595 << triQ[minIndex] << ". This might give problems in"
596 << " self-intersection testing later on." << endl;
597 }
598
599 // Dump for subsetting
600 if (surfWriter)
601 {
602 // Transcribe faces
603 faceList faces(surf.size());
604 surf.triFaceFaces(faces);
605
606 surfWriter->open
607 (
608 surf.points(),
609 faces,
610 (surfFilePath / surfFileNameBase),
611 false // serial - already merged
612 );
613
614 fileName outputName = surfWriter->write("quality", triQ);
615
616 surfWriter->clear();
617
618 Info<< "Wrote triangle-quality to "
619 << outputName << nl << endl;
620 }
621 else if (outputThreshold > 0)
622 {
623 DynamicList<label> problemFaces(surf.size()/100+1);
624
625 forAll(triQ, facei)
626 {
627 if (triQ[facei] < 1e-11)
628 {
629 problemFaces.append(facei);
630 }
631 }
632
633 if (!problemFaces.empty())
634 {
635 OFstream str("badFaces");
636
637 Info<< "Dumping bad quality faces to " << str.name() << endl
638 << "Paste this into the input for surfaceSubset" << nl
639 << nl << endl;
640
641 str << problemFaces;
642 }
643 }
644 }
645
646
647
648 // Edges
649 // ~~~~~
650 {
651 const edgeList& edges = surf.edges();
652 const pointField& localPoints = surf.localPoints();
653
654 scalarField edgeMag(edges.size());
655
656 forAll(edges, edgei)
657 {
658 edgeMag[edgei] = edges[edgei].mag(localPoints);
659 }
660
661 labelPair minMaxIds = findMinMax(edgeMag);
662
663 const label minEdgei = minMaxIds.first();
664 const label maxEdgei = minMaxIds.second();
665
666 const edge& minE = edges[minEdgei];
667 const edge& maxE = edges[maxEdgei];
668
669
670 Info<< "Edges:" << nl
671 << " min " << edgeMag[minEdgei] << " for edge " << minEdgei
672 << " points " << localPoints[minE[0]] << localPoints[minE[1]]
673 << nl
674 << " max " << edgeMag[maxEdgei] << " for edge " << maxEdgei
675 << " points " << localPoints[maxE[0]] << localPoints[maxE[1]]
676 << nl
677 << endl;
678 }
679
680
681
682 // Close points
683 // ~~~~~~~~~~~~
684 {
685 const edgeList& edges = surf.edges();
686 const pointField& localPoints = surf.localPoints();
687
688 const boundBox bb(localPoints);
689 scalar smallDim = 1e-6 * bb.mag();
690
691 Info<< "Checking for points less than 1e-6 of bounding box ("
692 << bb.span() << " metre) apart."
693 << endl;
694
695 // Sort points
696 SortableList<scalar> sortedMag(mag(localPoints));
697
698 label nClose = 0;
699
700 for (label i = 1; i < sortedMag.size(); i++)
701 {
702 label pti = sortedMag.indices()[i];
703
704 label prevPti = sortedMag.indices()[i-1];
705
706 if (mag(localPoints[pti] - localPoints[prevPti]) < smallDim)
707 {
708 // Check if neighbours.
709 const labelList& pEdges = surf.pointEdges()[pti];
710
711 label edgei = -1;
712
713 forAll(pEdges, i)
714 {
715 const edge& e = edges[pEdges[i]];
716
717 if (e[0] == prevPti || e[1] == prevPti)
718 {
719 // point1 and point0 are connected through edge.
720 edgei = pEdges[i];
721
722 break;
723 }
724 }
725
726 if (nClose < outputThreshold)
727 {
728 if (edgei == -1)
729 {
730 Info<< " close unconnected points "
731 << pti << ' ' << localPoints[pti]
732 << " and " << prevPti << ' '
733 << localPoints[prevPti]
734 << " distance:"
735 << mag(localPoints[pti] - localPoints[prevPti])
736 << endl;
737 }
738 else
739 {
740 Info<< " small edge between points "
741 << pti << ' ' << localPoints[pti]
742 << " and " << prevPti << ' '
743 << localPoints[prevPti]
744 << " distance:"
745 << mag(localPoints[pti] - localPoints[prevPti])
746 << endl;
747 }
748 }
749 else if (nClose == outputThreshold)
750 {
751 Info<< "Suppressing further warning messages."
752 << " Use -outputThreshold to increase number"
753 << " of warnings." << endl;
754 }
755
756 nClose++;
757 }
758 }
759
760 Info<< "Found " << nClose << " nearby points." << nl
761 << endl;
762 }
763
764
765
766 // Check manifold
767 // ~~~~~~~~~~~~~~
768
769 DynamicList<label> problemFaces(surf.size()/100 + 1);
770 DynamicList<label> openEdges(surf.nEdges()/100 + 1);
771 DynamicList<label> multipleEdges(surf.nEdges()/100 + 1);
772
773 const labelListList& edgeFaces = surf.edgeFaces();
774
775 forAll(edgeFaces, edgei)
776 {
777 const labelList& myFaces = edgeFaces[edgei];
778
779 if (myFaces.size() == 1)
780 {
781 problemFaces.append(myFaces[0]);
782 openEdges.append(edgei);
783 }
784 }
785
786 forAll(edgeFaces, edgei)
787 {
788 const labelList& myFaces = edgeFaces[edgei];
789
790 if (myFaces.size() > 2)
791 {
792 forAll(myFaces, myFacei)
793 {
794 problemFaces.append(myFaces[myFacei]);
795 }
796 multipleEdges.append(edgei);
797 }
798 }
799 problemFaces.shrink();
800
801 if (openEdges.size() || multipleEdges.size())
802 {
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;
808
809 Info<< "Conflicting face labels:" << problemFaces.size() << endl;
810
811 if (!edgeFormat.empty() && openEdges.size())
812 {
813 const fileName openName
814 (
815 surfFileName.lessExt()
816 + "_open."
817 + edgeFormat
818 );
819 Info<< "Writing open edges to " << openName << " ..." << endl;
820 writeEdgeSet(openName, surf, openEdges);
821 }
822 if (!edgeFormat.empty() && multipleEdges.size())
823 {
824 const fileName multName
825 (
826 surfFileName.lessExt()
827 + "_multiply."
828 + edgeFormat
829 );
830 Info<< "Writing multiply connected edges to "
831 << multName << " ..." << endl;
832 writeEdgeSet(multName, surf, multipleEdges);
833 }
834 }
835 else
836 {
837 Info<< "Surface is closed. All edges connected to two faces." << endl;
838 }
839 Info<< endl;
840
841
842
843 // Check singly connected domain
844 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
845
846 {
847 boolList borderEdge(surf.nEdges(), false);
848 if (splitNonManifold)
849 {
850 forAll(edgeFaces, edgei)
851 {
852 if (edgeFaces[edgei].size() > 2)
853 {
854 borderEdge[edgei] = true;
855 }
856 }
857 syncEdges(surf, borderEdge);
858 }
859
861 label numZones = surf.markZones(borderEdge, faceZone);
862
863 Info<< "Number of unconnected parts : " << numZones << endl;
864
865 if (numZones > 1 && outputThreshold > 0)
866 {
867 Info<< "Splitting surface into parts ..." << endl << endl;
868
869 if (!surfWriter)
870 {
871 surfWriter.reset(new surfaceWriters::vtkWriter());
872 }
873
874 writeZoning
875 (
876 *surfWriter,
877 surf,
878 faceZone,
879 "zone",
880 surfFilePath,
881 surfFileNameBase
882 );
883
884 if (numZones > outputThreshold)
885 {
886 Info<< "Limiting number of files to " << outputThreshold
887 << endl;
888 }
889 writeParts
890 (
891 surf,
892 min(outputThreshold, numZones),
893 faceZone,
894 surfFilePath,
895 surfFileNameBase
896 );
897 }
898 }
899
900
901
902 // Check orientation
903 // ~~~~~~~~~~~~~~~~~
904
905 labelHashSet borderEdge(surf.size()/1000);
906 PatchTools::checkOrientation(surf, false, &borderEdge);
907
908 // Bit strange: if a triangle has two same vertices (illegal!) it will
909 // still have three distinct edges (two of which have the same vertices).
910 // In this case the faceEdges addressing is not symmetric, i.e. a
911 // neighbouring, valid, triangle will have correct addressing so 3 distinct
912 // edges so it will miss one of those two identical edges.
913 // - we don't want to fix this in PrimitivePatch since it is too specific
914 // - instead just make sure we mark all identical edges consistently
915 // when we use them for marking.
916
917 syncEdges(surf, borderEdge);
918
919 //
920 // Colour all faces into zones using borderEdge
921 //
922 labelList normalZone;
923 label numNormalZones = PatchTools::markZones(surf, borderEdge, normalZone);
924
925 Info<< endl
926 << "Number of zones (connected area with consistent normal) : "
927 << numNormalZones << endl;
928
929 if (numNormalZones > 1)
930 {
931 Info<< "More than one normal orientation." << endl;
932
933 if (outputThreshold > 0)
934 {
935 if (!surfWriter)
936 {
937 surfWriter.reset(new surfaceWriters::vtkWriter());
938 }
939
940 writeZoning
941 (
942 *surfWriter,
943 surf,
944 normalZone,
945 "normal",
946 surfFilePath,
947 surfFileNameBase
948 );
949
950 if (numNormalZones > outputThreshold)
951 {
952 Info<< "Limiting number of files to " << outputThreshold
953 << endl;
954 }
955 writeParts
956 (
957 surf,
958 min(outputThreshold, numNormalZones),
959 normalZone,
960 surfFilePath,
961 surfFileNameBase + "_normal"
962 );
963 }
964 }
965 Info<< endl;
966
967
968
969 // Check self-intersection
970 // ~~~~~~~~~~~~~~~~~~~~~~~
971
972 if (checkSelfIntersect)
973 {
974 Info<< "Checking self-intersection." << endl;
975
976 triSurfaceSearch querySurf(surf);
977
978 const indexedOctree<treeDataTriSurface>& tree = querySurf.tree();
979
980 autoPtr<OBJstream> intStreamPtr;
981 if (outputThreshold > 0)
982 {
983 intStreamPtr.reset(new OBJstream("selfInterPoints.obj"));
984 }
985
986 label nInt = 0;
987
988 forAll(surf.edges(), edgei)
989 {
990 const edge& e = surf.edges()[edgei];
991 const point& start = surf.points()[surf.meshPoints()[e[0]]];
992 const point& end = surf.points()[surf.meshPoints()[e[1]]];
993
994 // Exclude hits of connected triangles
995 treeDataTriSurface::findSelfIntersectOp exclOp(tree, edgei);
996
997 pointIndexHit hitInfo(tree.findLineAny(start, end, exclOp));
998
999 if (hitInfo.hit())
1000 {
1001 nInt++;
1002
1003 if (intStreamPtr)
1004 {
1005 intStreamPtr().write(hitInfo.hitPoint());
1006 }
1007
1008 // Try and find from other side.
1009 pointIndexHit hitInfo2(tree.findLineAny(end, start, exclOp));
1010
1011 if (hitInfo2.hit() && hitInfo.index() != hitInfo2.index())
1012 {
1013 nInt++;
1014
1015 if (intStreamPtr)
1016 {
1017 intStreamPtr().write(hitInfo2.hitPoint());
1018 }
1019 }
1020 }
1021 }
1022
1024 //{
1025 // const pointField& localPoints = surf.localPoints();
1026 //
1027 // const boundBox bb(localPoints);
1028 // scalar smallDim = 1e-6 * bb.mag();
1029 // scalar smallDimSqr = Foam::sqr(smallDim);
1030 //
1031 // const pointField& faceCentres = surf.faceCentres();
1032 // forAll(faceCentres, faceI)
1033 // {
1034 // const point& fc = faceCentres[faceI];
1035 // pointIndexHit hitInfo
1036 // (
1037 // tree.findNearest
1038 // (
1039 // fc,
1040 // smallDimSqr,
1041 // findSelfNearOp(tree, faceI)
1042 // )
1043 // );
1044 //
1045 // if (hitInfo.hit() && intStreamPtr)
1046 // {
1047 // intStreamPtr().write(hitInfo.hitPoint());
1048 //
1049 // label nearFaceI = hitInfo.index();
1050 // triPointRef nearTri(surf[nearFaceI].tri(surf.points()));
1051 // triStreamPtr().write
1052 // (
1053 // surf[faceI].tri(surf.points()),
1054 // false
1055 // );
1056 // triStreamPtr().write(nearTri, false);
1057 // nInt++;
1058 // }
1059 // }
1060 //}
1061
1062
1063 if (nInt == 0)
1064 {
1065 Info<< "Surface is not self-intersecting" << endl;
1066 }
1067 else
1068 {
1069 Info<< "Surface is self-intersecting at " << nInt
1070 << " locations." << endl;
1071
1072 if (intStreamPtr)
1073 {
1074 Info<< "Writing intersection points to "
1075 << intStreamPtr().name() << endl;
1076 }
1077 }
1078 Info<< endl;
1079 }
1080
1081
1082 Info<< "\nEnd\n" << endl;
1083
1084 return 0;
1085}
1086
1087
1088// ************************************************************************* //
label n
Y[inertIndex] max(0.0)
vtk::internalMeshWriter writer(topoMesh, topoCells, vtk::formatType::INLINE_ASCII, runTime.path()/"blockTopology")
Dynamically sized Field.
Definition: DynamicField.H:65
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
T & first() noexcept
The first element of the list, position [0].
Definition: FixedListI.H:207
label size() const noexcept
The number of elements in table.
Definition: HashTableI.H:52
void append(const T &val)
Append an element at the end of the list.
Definition: ListI.H:175
OFstream that keeps track of vertices.
Definition: OBJstream.H:61
Output to file stream, using an OSstream.
Definition: OFstream.H:57
Generic output stream using a standard (STL) stream.
Definition: OSstream.H:57
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.
Definition: PairI.H:120
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:66
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.
Definition: SortableList.H:63
A List with indirect addressing. Like IndirectList but does not store addressing.
Definition: IndirectList.H:79
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Extract command arguments and options from the supplied argc and argv parameters.
Definition: argList.H:124
T get(const label index) const
Get a value from the argument at index.
Definition: argListI.H:278
bool found(const word &optName) const
Return true if the named option is found.
Definition: argListI.H:178
T getOrDefault(const word &optName, const T &deflt) const
Get a value from the named option if present, or return default.
Definition: argListI.H:307
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
void clear() noexcept
Same as reset(nullptr)
Definition: autoPtr.H:176
void reset(autoPtr< T > &&other) noexcept
Delete managed object and set to new given pointer.
Definition: autoPtrI.H:117
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:64
tmp< pointField > points() const
Corner points in an order corresponding to a 'hex' cell.
Definition: boundBox.C:118
Mesh data needed to do the Finite Area discretisation.
Definition: edgeFaMesh.H:56
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:66
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:67
A class for handling file names.
Definition: fileName.H:76
fileName lessExt() const
Return file name without extension (part before last .)
Definition: fileNameI.H:230
word ext() const
Return file name extension (part after last .)
Definition: fileNameI.H:218
Non-pointer based hierarchical recursive searching.
Definition: indexedOctree.H:74
pointIndexHit findLineAny(const point &start, const point &end) const
Find any intersection of line between start and end.
A triFace with additional (region) index.
Definition: labelledTri.H:60
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.
Definition: triSurface.H:79
label markZones(const boolList &borderEdge, labelList &faceZone) const
(size and) fills faceZone with zone of face. Zone is area
Definition: triSurface.C:796
void triFaceFaces(List< face > &plainFaceList) const
Create a list of faces from the triFaces.
Definition: triSurface.C:723
triSurface subsetMesh(const UList< bool > &include, labelList &pointMap, labelList &faceMap) const
Return a new surface subsetted on the selected faces.
Definition: triSurface.C:879
const geometricSurfacePatchList & patches() const noexcept
Definition: triSurface.H:399
void writeStats(Ostream &os) const
Write some statistics.
Definition: triSurfaceIO.C:353
A triangle primitive used to calculate face normals and swept volumes.
Definition: triangle.H:80
scalar quality() const
Return quality: Ratio of triangle and circum-circle.
Definition: triangleI.H:183
A class for handling words, derived from Foam::string.
Definition: word.H:68
Base class for mesh zones.
Definition: zone.H:67
volScalarField & p
word outputName("finiteArea-edges.obj")
OBJstream os(runTime.globalPath()/outputName)
#define WarningInFunction
Report a warning using Foam::Warning.
Namespace for OpenFOAM.
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.
Definition: Ostream.H:372
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.
Definition: hashSets.C:33
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Definition: stdFoam.H:158
labelList f(nPoints)
Foam::argList args(argc, argv)
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333