67 bool validTri(
const bool verbose,
const triSurface& surf,
const label faceI)
76 || (f[1] < 0) || (f[1] >= surf.
points().
size())
77 || (f[2] < 0) || (f[2] >= surf.
points().
size())
80 WarningIn(
"validTri(const triSurface&, const label)")
81 <<
"triangle " << faceI <<
" vertices " << f
82 <<
" uses point indices outside point range 0.."
88 if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
90 WarningIn(
"validTri(const triSurface&, const label)")
91 <<
"triangle " << faceI
92 <<
" uses non-unique vertices " << f
106 label nbrFaceI = fFaces[i];
108 if (nbrFaceI <= faceI)
118 ((f[0] == nbrF[0]) || (f[0] == nbrF[1]) || (f[0] == nbrF[2]))
119 && ((f[1] == nbrF[0]) || (f[1] == nbrF[1]) || (f[1] == nbrF[2]))
120 && ((f[2] == nbrF[0]) || (f[2] == nbrF[1]) || (f[2] == nbrF[2]))
123 WarningIn(
"validTri(const triSurface&, const label)")
124 <<
"triangle " << faceI <<
" vertices " << f
125 <<
" has the same vertices as triangle " << nbrFaceI
126 <<
" vertices " << nbrF
145 scalar
dist = nBins/(max -
min);
151 scalar val = vals[i];
159 else if (val >= max - SMALL)
165 index = label((val - min)*dist);
167 if ((index < 0) || (index >= nBins))
171 "countBins(const scalar, const scalar, const label"
172 ", const scalarField&)"
173 ) <<
"value " << val <<
" at index " << i
174 <<
" outside range " << min <<
" .. " << max <<
endl;
196 int main(
int argc,
char *argv[])
206 bool checkSelfIntersection =
args.
optionFound(
"checkSelfIntersection");
210 Pout<<
"Reading surface from " << surfFileName <<
" ..." <<
nl <<
endl;
232 label region = surf[faceI].region();
234 if (region < 0 || region >= regionSize.size())
237 <<
"Triangle " << faceI <<
" vertices " << surf[faceI]
238 <<
" has region " << region <<
" which is outside the range"
244 regionSize[region]++;
248 Pout<<
"Region\tSize" <<
nl
249 <<
"------\t----" <<
nl;
253 << regionSize[patchI] <<
nl;
267 if (!validTri(verbose, surf, faceI))
269 illegalFaces.append(faceI);
273 if (illegalFaces.size())
275 Pout<<
"Surface has " << illegalFaces.size()
276 <<
" illegal triangles." <<
endl;
279 Pout<<
"Dumping conflicting face labels to " << str.
name() << endl
280 <<
"Paste this into the input for surfaceSubset" <<
endl;
285 Pout<<
"Surface has no illegal triangles." <<
endl;
301 if (f[0] == f[1] || f[0] == f[2] || f[1] == f[2])
316 vector ba(tri.b() - tri.a());
317 ba /=
mag(ba) + VSMALL;
319 vector ca(tri.c() - tri.a());
320 ca /=
mag(ca) + VSMALL;
322 if (
mag(ba&ca) > 1-1
E-3)
338 labelList binCount = countBins(0, 1, 20, triQ);
340 Pout<<
"Triangle quality (equilateral=1, collapsed=0):"
347 scalar dist = (1.0 - 0.0)/20.0;
351 Pout<<
" " << min <<
" .. " << min+dist <<
" : "
352 << 1.0/surf.
size() * binCount[binI]
358 label minIndex =
findMin(triQ);
359 label maxIndex =
findMax(triQ);
361 Pout<<
" min " << triQ[minIndex] <<
" for triangle " << minIndex
363 <<
" max " << triQ[maxIndex] <<
" for triangle " << maxIndex
368 if (triQ[minIndex] < SMALL)
371 << triQ[minIndex] <<
". This might give problems in"
372 <<
" self-intersection testing later on." << endl;
381 if (triQ[faceI] < 1
E-11)
383 problemFaces.append(faceI);
388 Pout<<
"Dumping bad quality faces to " << str.
name() << endl
389 <<
"Paste this into the input for surfaceSubset" << nl
408 edgeMag[edgeI] = edges[edgeI].mag(localPoints);
411 label minEdgeI =
findMin(edgeMag);
412 label maxEdgeI =
findMax(edgeMag);
414 const edge& minE = edges[minEdgeI];
415 const edge& maxE = edges[maxEdgeI];
418 Pout<<
"Edges:" << nl
419 <<
" min " << edgeMag[minEdgeI] <<
" for edge " << minEdgeI
420 <<
" points " << localPoints[minE[0]] << localPoints[minE[1]]
422 <<
" max " << edgeMag[maxEdgeI] <<
" for edge " << maxEdgeI
423 <<
" points " << localPoints[maxE[0]] << localPoints[maxE[1]]
437 scalar smallDim = 1
E-6 * bb.mag();
439 Pout<<
"Checking for points less than 1E-6 of bounding box ("
440 << bb.span() <<
" meter) apart."
448 for (label i = 1; i < sortedMag.size(); i++)
450 label ptI = sortedMag.indices()[i];
452 label prevPtI = sortedMag.indices()[i-1];
454 if (
mag(localPoints[ptI] - localPoints[prevPtI]) < smallDim)
463 const edge&
e = edges[pEdges[i]];
465 if (e[0] == prevPtI || e[1] == prevPtI)
478 Pout<<
" close unconnected points "
479 << ptI <<
' ' << localPoints[ptI]
480 <<
" and " << prevPtI <<
' '
481 << localPoints[prevPtI]
483 <<
mag(localPoints[ptI] - localPoints[prevPtI])
488 Pout<<
" small edge between points "
489 << ptI <<
' ' << localPoints[ptI]
490 <<
" and " << prevPtI <<
' '
491 << localPoints[prevPtI]
493 <<
mag(localPoints[ptI] - localPoints[prevPtI])
499 Pout<<
"Found " << nClose <<
" nearby points." << nl
512 label nSingleEdges = 0;
515 const labelList& myFaces = eFaces[edgeI];
517 if (myFaces.
size() == 1)
519 problemFaces.
append(myFaces[0]);
525 label nMultEdges = 0;
528 const labelList& myFaces = eFaces[edgeI];
530 if (myFaces.
size() > 2)
534 problemFaces.append(myFaces[myFaceI]);
540 problemFaces.shrink();
542 if ((nSingleEdges != 0) || (nMultEdges != 0))
544 Pout<<
"Surface is not closed since not all edges connected to "
545 <<
"two faces:" << endl
546 <<
" connected to one face : " << nSingleEdges << endl
547 <<
" connected to >2 faces : " << nMultEdges <<
endl;
549 Pout<<
"Conflicting face labels:" << problemFaces.size() <<
endl;
553 Pout<<
"Dumping conflicting face labels to " << str.
name() << endl
554 <<
"Paste this into the input for surfaceSubset" <<
endl;
560 Pout<<
"Surface is closed. All edges connected to two faces." <<
endl;
572 Pout<<
"Number of unconnected parts : " << numZones <<
endl;
576 Pout<<
"Splitting surface into parts ..." << endl <<
endl;
578 fileName surfFileNameBase(surfFileName.name());
580 for(label zone = 0; zone < numZones; zone++)
586 if (faceZone[faceI] == zone)
588 includeMap[faceI] =
true;
607 surfFileNameBase.lessExt()
613 Pout<<
"writing part " << zone <<
" size " << subSurf.size()
614 <<
" to " << subFileName <<
endl;
616 subSurf.
write(subFileName);
637 <<
"Number of zones (connected area with consistent normal) : "
638 << numNormalZones <<
endl;
640 if (numNormalZones > 1)
642 Pout<<
"More than one normal orientation." <<
endl;
651 if (checkSelfIntersection)
653 Pout<<
"Checking self-intersection." <<
endl;
658 if (inter.cutEdges().empty() && inter.cutPoints().empty())
660 Pout<<
"Surface is not self-intersecting" <<
endl;
664 Pout<<
"Surface is self-intersecting" <<
endl;
665 Pout<<
"Writing edges of intersection to selfInter.obj" <<
endl;
667 OFstream intStream(
"selfInter.obj");
668 forAll(inter.cutPoints(), cutPointI)
670 const point& pt = inter.cutPoints()[cutPointI];
672 intStream <<
"v " << pt.
x() <<
' ' << pt.
y() <<
' ' << pt.
z()
675 forAll(inter.cutEdges(), cutEdgeI)
677 const edge&
e = inter.cutEdges()[cutEdgeI];
679 intStream <<
"l " << e.
start()+1 <<
' ' << e.
end()+1 <<
endl;