applications/utilities/surface/surfaceCheck/surfaceCheck.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "triangle.H"
  27 #include "triSurface.H"
  28 #include "triSurfaceTools.H"
  29 #include "triSurfaceSearch.H"
  30 #include "argList.H"
  31 #include "OFstream.H"
  32 #include "surfaceIntersection.H"
  33 #include "SortableList.H"
  34 #include "PatchTools.H"
  35
  36 using namespace Foam;
  37
  38 // Does face use valid vertices?
  39 bool validTri
  40 (
  41     const bool verbose,
  42     const triSurface& surf,
  43     const label faceI
  44 )
  45 {
  46     // Simple check on indices ok.
  47
  48     const labelledTri& f = surf[faceI];
  49
  50     forAll(f, fp)
  51     {
  52         if (f[fp] < 0 || f[fp] >= surf.points().size())
  53         {
  54             WarningIn("validTri(const triSurface&, const label)")
  55                 << "triangle " << faceI << " vertices " << f
  56                 << " uses point indices outside point range 0.."
  57                 << surf.points().size()-1 << endl;
  58             return false;
  59         }
  60     }
  61
  62     if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
  63     {
  64         WarningIn("validTri(const triSurface&, const label)")
  65             << "triangle " << faceI
  66             << " uses non-unique vertices " << f
  67             << " coords:" << f.points(surf.points())
  68             << endl;
  69         return false;
  70     }
  71
  72     // duplicate triangle check
  73
  74     const labelList& fFaces = surf.faceFaces()[faceI];
  75
  76     // Check if faceNeighbours use same points as this face.
  77     // Note: discards normal information - sides of baffle are merged.
  78     forAll(fFaces, i)
  79     {
  80         label nbrFaceI = fFaces[i];
  81
  82         if (nbrFaceI <= faceI)
  83         {
  84             // lower numbered faces already checked
  85             continue;
  86         }
  87
  88         const labelledTri& nbrF = surf[nbrFaceI];
  89
  90         if
  91         (
  92             ((f[0] == nbrF[0]) || (f[0] == nbrF[1]) || (f[0] == nbrF[2]))
  93          && ((f[1] == nbrF[0]) || (f[1] == nbrF[1]) || (f[1] == nbrF[2]))
  94          && ((f[2] == nbrF[0]) || (f[2] == nbrF[1]) || (f[2] == nbrF[2]))
  95         )
  96         {
  97             WarningIn("validTri(const triSurface&, const label)")
  98                 << "triangle " << faceI << " vertices " << f
  99                 << " has the same vertices as triangle " << nbrFaceI
 100                 << " vertices " << nbrF
 101                 << " coords:" << f.points(surf.points())
 102                 << endl;
 103
 104             return false;
 105         }
 106     }
 107     return true;
 108 }
 109
 110
 111 labelList countBins
 112 (
 113     const scalar min,
 114     const scalar max,
 115     const label nBins,
 116     const scalarField& vals
 117 )
 118 {
 119     scalar dist = nBins/(max - min);
 120
 121     labelList binCount(nBins, 0);
 122
 123     forAll(vals, i)
 124     {
 125         scalar val = vals[i];
 126
 127         label index = -1;
 128
 129         if (Foam::mag(val - min) < SMALL)
 130         {
 131             index = 0;
 132         }
 133         else if (val >= max - SMALL)
 134         {
 135             index = nBins - 1;
 136         }
 137         else
 138         {
 139             index = label((val - min)*dist);
 140
 141             if ((index < 0) || (index >= nBins))
 142             {
 143                 WarningIn
 144                 (
 145                     "countBins(const scalar, const scalar, const label"
 146                     ", const scalarField&)"
 147                 )   << "value " << val << " at index " << i
 148                     << " outside range " << min << " .. " << max << endl;
 149
 150                 if (index < 0)
 151                 {
 152                     index = 0;
 153                 }
 154                 else
 155                 {
 156                     index = nBins - 1;
 157                 }
 158             }
 159         }
 160         binCount[index]++;
 161     }
 162
 163     return binCount;
 164 }
 165
 166
 167 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 168 // Main program:
 169
 170 int main(int argc, char *argv[])
 171 {
 172     argList::noParallel();
 173     argList::validArgs.append("surfaceFile");
 174     argList::addBoolOption
 175     (
 176         "checkSelfIntersection",
 177         "also check for self-intersection"
 178     );
 179     argList::addBoolOption
 180     (
 181         "verbose",
 182         "verbose operation"
 183     );
 184     argList::addBoolOption
 185     (
 186         "blockMesh",
 187         "write vertices/blocks for blockMeshDict"
 188     );
 189
 190     argList args(argc, argv);
 191
 192     const fileName surfFileName = args[1];
 193     const bool checkSelfIntersect = args.optionFound("checkSelfIntersection");
 194     const bool verbose = args.optionFound("verbose");
 195
 196     Info<< "Reading surface from " << surfFileName << " ..." << nl << endl;
 197
 198
 199     // Read
 200     // ~~~~
 201
 202     triSurface surf(surfFileName);
 203
 204
 205     Info<< "Statistics:" << endl;
 206     surf.writeStats(Info);
 207     Info<< endl;
 208
 209     // write bounding box corners
 210     if (args.optionFound("blockMesh"))
 211     {
 212         pointField cornerPts(boundBox(surf.points()).points());
 213
 214         Info<<"// blockMeshDict info" << nl
 215             <<"vertices\n(" << nl;
 216         forAll(cornerPts, ptI)
 217         {
 218             Info << "    " << cornerPts[ptI] << nl;
 219         }
 220
 221         // number of divisions needs adjustment later
 222         Info<<");\n" << nl
 223             <<"blocks\n"
 224             <<"(\n"
 225             <<"    hex (0 1 2 3 4 5 6 7) (10 10 10) simpleGrading (1 1 1)\n"
 226             <<");\n" << nl;
 227
 228         Info<<"edges\n();" << nl
 229             <<"patches\n();" << endl;
 230     }
 231
 232
 233     // Region sizes
 234     // ~~~~~~~~~~~~
 235
 236     {
 237         labelList regionSize(surf.patches().size(), 0);
 238
 239         forAll(surf, faceI)
 240         {
 241             label region = surf[faceI].region();
 242
 243             if (region < 0 || region >= regionSize.size())
 244             {
 245                 WarningIn(args.executable())
 246                     << "Triangle " << faceI << " vertices " << surf[faceI]
 247                     << " has region " << region << " which is outside the range"
 248                     << " of regions 0.." << surf.patches().size()-1
 249                     << endl;
 250             }
 251             else
 252             {
 253                 regionSize[region]++;
 254             }
 255         }
 256
 257         Info<< "Region\tSize" << nl
 258             << "------\t----" << nl;
 259         forAll(surf.patches(), patchI)
 260         {
 261             Info<< surf.patches()[patchI].name() << '\t'
 262                 << regionSize[patchI] << nl;
 263         }
 264         Info<< nl << endl;
 265     }
 266
 267
 268     // Check triangles
 269     // ~~~~~~~~~~~~~~~
 270
 271     {
 272         DynamicList<label> illegalFaces(surf.size()/100 + 1);
 273
 274         forAll(surf, faceI)
 275         {
 276             if (!validTri(verbose, surf, faceI))
 277             {
 278                 illegalFaces.append(faceI);
 279             }
 280         }
 281
 282         if (illegalFaces.size())
 283         {
 284             Info<< "Surface has " << illegalFaces.size()
 285                 << " illegal triangles." << endl;
 286
 287             OFstream str("illegalFaces");
 288             Info<< "Dumping conflicting face labels to " << str.name() << endl
 289                 << "Paste this into the input for surfaceSubset" << endl;
 290             str << illegalFaces;
 291         }
 292         else
 293         {
 294             Info<< "Surface has no illegal triangles." << endl;
 295         }
 296         Info<< endl;
 297     }
 298
 299
 300
 301     // Triangle quality
 302     // ~~~~~~~~~~~~~~~~
 303
 304     {
 305         scalarField triQ(surf.size(), 0);
 306         forAll(surf, faceI)
 307         {
 308             const labelledTri& f = surf[faceI];
 309
 310             if (f[0] == f[1] || f[0] == f[2] || f[1] == f[2])
 311             {
 312                 //WarningIn(args.executable())
 313                 //    << "Illegal triangle " << faceI << " vertices " << f
 314                 //    << " coords " << f.points(surf.points()) << endl;
 315             }
 316             else
 317             {
 318                 triPointRef tri
 319                 (
 320                     surf.points()[f[0]],
 321                     surf.points()[f[1]],
 322                     surf.points()[f[2]]
 323                 );
 324
 325                 vector ba(tri.b() - tri.a());
 326                 ba /= mag(ba) + VSMALL;
 327
 328                 vector ca(tri.c() - tri.a());
 329                 ca /= mag(ca) + VSMALL;
 330
 331                 if (mag(ba&ca) > 1-1E-3)
 332                 {
 333                     triQ[faceI] = SMALL;
 334                 }
 335                 else
 336                 {
 337                     triQ[faceI] = triPointRef
 338                     (
 339                         surf.points()[f[0]],
 340                         surf.points()[f[1]],
 341                         surf.points()[f[2]]
 342                     ).quality();
 343                 }
 344             }
 345         }
 346
 347         labelList binCount = countBins(0, 1, 20, triQ);
 348
 349         Info<< "Triangle quality (equilateral=1, collapsed=0):"
 350             << endl;
 351
 352
 353         OSstream& os = Info;
 354         os.width(4);
 355
 356         scalar dist = (1.0 - 0.0)/20.0;
 357         scalar min = 0;
 358         forAll(binCount, binI)
 359         {
 360             Info<< "    " << min << " .. " << min+dist << "  : "
 361                 << 1.0/surf.size() * binCount[binI]
 362                 << endl;
 363             min += dist;
 364         }
 365         Info<< endl;
 366
 367         label minIndex = findMin(triQ);
 368         label maxIndex = findMax(triQ);
 369
 370         Info<< "    min " << triQ[minIndex] << " for triangle " << minIndex
 371             << nl
 372             << "    max " << triQ[maxIndex] << " for triangle " << maxIndex
 373             << nl
 374             << endl;
 375
 376
 377         if (triQ[minIndex] < SMALL)
 378         {
 379             WarningIn(args.executable()) << "Minimum triangle quality is "
 380                 << triQ[minIndex] << ". This might give problems in"
 381                 << " self-intersection testing later on." << endl;
 382         }
 383
 384         // Dump for subsetting
 385         {
 386             DynamicList<label> problemFaces(surf.size()/100+1);
 387
 388             forAll(triQ, faceI)
 389             {
 390                 if (triQ[faceI] < 1E-11)
 391                 {
 392                     problemFaces.append(faceI);
 393                 }
 394             }
 395
 396             if (!problemFaces.empty())
 397             {
 398                 OFstream str("badFaces");
 399
 400                 Info<< "Dumping bad quality faces to " << str.name() << endl
 401                     << "Paste this into the input for surfaceSubset" << nl
 402                     << nl << endl;
 403
 404                 str << problemFaces;
 405             }
 406         }
 407     }
 408
 409
 410
 411     // Edges
 412     // ~~~~~
 413     {
 414         const edgeList& edges = surf.edges();
 415         const pointField& localPoints = surf.localPoints();
 416
 417         scalarField edgeMag(edges.size());
 418
 419         forAll(edges, edgeI)
 420         {
 421             edgeMag[edgeI] = edges[edgeI].mag(localPoints);
 422         }
 423
 424         label minEdgeI = findMin(edgeMag);
 425         label maxEdgeI = findMax(edgeMag);
 426
 427         const edge& minE = edges[minEdgeI];
 428         const edge& maxE = edges[maxEdgeI];
 429
 430
 431         Info<< "Edges:" << nl
 432             << "    min " << edgeMag[minEdgeI] << " for edge " << minEdgeI
 433             << " points " << localPoints[minE[0]] << localPoints[minE[1]]
 434             << nl
 435             << "    max " << edgeMag[maxEdgeI] << " for edge " << maxEdgeI
 436             << " points " << localPoints[maxE[0]] << localPoints[maxE[1]]
 437             << nl
 438             << endl;
 439     }
 440
 441
 442
 443     // Close points
 444     // ~~~~~~~~~~~~
 445     {
 446         const edgeList& edges = surf.edges();
 447         const pointField& localPoints = surf.localPoints();
 448
 449         const boundBox bb(localPoints);
 450         scalar smallDim = 1E-6 * bb.mag();
 451
 452         Info<< "Checking for points less than 1E-6 of bounding box ("
 453             << bb.span() << " meter) apart."
 454             << endl;
 455
 456         // Sort points
 457         SortableList<scalar> sortedMag(mag(localPoints));
 458
 459         label nClose = 0;
 460
 461         for (label i = 1; i < sortedMag.size(); i++)
 462         {
 463             label ptI = sortedMag.indices()[i];
 464
 465             label prevPtI = sortedMag.indices()[i-1];
 466
 467             if (mag(localPoints[ptI] - localPoints[prevPtI]) < smallDim)
 468             {
 469                 // Check if neighbours.
 470                 const labelList& pEdges = surf.pointEdges()[ptI];
 471
 472                 label edgeI = -1;
 473
 474                 forAll(pEdges, i)
 475                 {
 476                     const edge& e = edges[pEdges[i]];
 477
 478                     if (e[0] == prevPtI || e[1] == prevPtI)
 479                     {
 480                         // point1 and point0 are connected through edge.
 481                         edgeI = pEdges[i];
 482
 483                         break;
 484                     }
 485                 }
 486
 487                 nClose++;
 488
 489                 if (edgeI == -1)
 490                 {
 491                     Info<< "    close unconnected points "
 492                         << ptI << ' ' << localPoints[ptI]
 493                         << " and " << prevPtI << ' '
 494                         << localPoints[prevPtI]
 495                         << " distance:"
 496                         << mag(localPoints[ptI] - localPoints[prevPtI])
 497                         << endl;
 498                 }
 499                 else
 500                 {
 501                     Info<< "    small edge between points "
 502                         << ptI << ' ' << localPoints[ptI]
 503                         << " and " << prevPtI << ' '
 504                         << localPoints[prevPtI]
 505                         << " distance:"
 506                         << mag(localPoints[ptI] - localPoints[prevPtI])
 507                         << endl;
 508                 }
 509             }
 510         }
 511
 512         Info<< "Found " << nClose << " nearby points." << nl
 513             << endl;
 514     }
 515
 516
 517
 518     // Check manifold
 519     // ~~~~~~~~~~~~~~
 520
 521     DynamicList<label> problemFaces(surf.size()/100 + 1);
 522
 523     const labelListList& eFaces = surf.edgeFaces();
 524
 525     label nSingleEdges = 0;
 526     forAll(eFaces, edgeI)
 527     {
 528         const labelList& myFaces = eFaces[edgeI];
 529
 530         if (myFaces.size() == 1)
 531         {
 532             problemFaces.append(myFaces[0]);
 533
 534             nSingleEdges++;
 535         }
 536     }
 537
 538     label nMultEdges = 0;
 539     forAll(eFaces, edgeI)
 540     {
 541         const labelList& myFaces = eFaces[edgeI];
 542
 543         if (myFaces.size() > 2)
 544         {
 545             forAll(myFaces, myFaceI)
 546             {
 547                 problemFaces.append(myFaces[myFaceI]);
 548             }
 549
 550             nMultEdges++;
 551         }
 552     }
 553     problemFaces.shrink();
 554
 555     if ((nSingleEdges != 0) || (nMultEdges != 0))
 556     {
 557         Info<< "Surface is not closed since not all edges connected to "
 558             << "two faces:" << endl
 559             << "    connected to one face : " << nSingleEdges << endl
 560             << "    connected to >2 faces : " << nMultEdges << endl;
 561
 562         Info<< "Conflicting face labels:" << problemFaces.size() << endl;
 563
 564         OFstream str("problemFaces");
 565
 566         Info<< "Dumping conflicting face labels to " << str.name() << endl
 567             << "Paste this into the input for surfaceSubset" << endl;
 568
 569         str << problemFaces;
 570     }
 571     else
 572     {
 573         Info<< "Surface is closed. All edges connected to two faces." << endl;
 574     }
 575     Info<< endl;
 576
 577
 578
 579     // Check singly connected domain
 580     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 581
 582     labelList faceZone;
 583     label numZones = surf.markZones(boolList(surf.nEdges(), false), faceZone);
 584
 585     Info<< "Number of unconnected parts : " << numZones << endl;
 586
 587     if (numZones > 1)
 588     {
 589         Info<< "Splitting surface into parts ..." << endl << endl;
 590
 591         fileName surfFileNameBase(surfFileName.name());
 592
 593         for (label zone = 0; zone < numZones; zone++)
 594         {
 595             boolList includeMap(surf.size(), false);
 596
 597             forAll(faceZone, faceI)
 598             {
 599                 if (faceZone[faceI] == zone)
 600                 {
 601                     includeMap[faceI] = true;
 602                 }
 603             }
 604
 605             labelList pointMap;
 606             labelList faceMap;
 607
 608             triSurface subSurf
 609             (
 610                 surf.subsetMesh
 611                 (
 612                     includeMap,
 613                     pointMap,
 614                     faceMap
 615                 )
 616             );
 617
 618             fileName subFileName
 619             (
 620                 surfFileNameBase.lessExt()
 621               + "_"
 622               + name(zone)
 623               + ".obj"
 624             );
 625
 626             Info<< "writing part " << zone << " size " << subSurf.size()
 627                 << " to " << subFileName << endl;
 628
 629             subSurf.write(subFileName);
 630         }
 631
 632         return 0;
 633     }
 634
 635
 636
 637     // Check orientation
 638     // ~~~~~~~~~~~~~~~~~
 639
 640     labelHashSet borderEdge(surf.size()/1000);
 641     PatchTools::checkOrientation(surf, false, &borderEdge);
 642
 643     //
 644     // Colour all faces into zones using borderEdge
 645     //
 646     labelList normalZone;
 647     label numNormalZones = PatchTools::markZones(surf, borderEdge, normalZone);
 648
 649     Info<< endl
 650         << "Number of zones (connected area with consistent normal) : "
 651         << numNormalZones << endl;
 652
 653     if (numNormalZones > 1)
 654     {
 655         Info<< "More than one normal orientation." << endl;
 656     }
 657     Info<< endl;
 658
 659
 660
 661     // Check self-intersection
 662     // ~~~~~~~~~~~~~~~~~~~~~~~
 663
 664     if (checkSelfIntersect)
 665     {
 666         Info<< "Checking self-intersection." << endl;
 667
 668         triSurfaceSearch querySurf(surf);
 669         surfaceIntersection inter(querySurf);
 670
 671         if (inter.cutEdges().empty() && inter.cutPoints().empty())
 672         {
 673             Info<< "Surface is not self-intersecting" << endl;
 674         }
 675         else
 676         {
 677             Info<< "Surface is self-intersecting" << endl;
 678             Info<< "Writing edges of intersection to selfInter.obj" << endl;
 679
 680             OFstream intStream("selfInter.obj");
 681             forAll(inter.cutPoints(), cutPointI)
 682             {
 683                 const point& pt = inter.cutPoints()[cutPointI];
 684
 685                 intStream << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z()
 686                     << endl;
 687             }
 688             forAll(inter.cutEdges(), cutEdgeI)
 689             {
 690                 const edge& e = inter.cutEdges()[cutEdgeI];
 691
 692                 intStream << "l " << e.start()+1 << ' ' << e.end()+1 << endl;
 693             }
 694         }
 695         Info<< endl;
 696     }
 697
 698
 699     Info<< "\nEnd\n" << endl;
 700
 701     return 0;
 702 }
 703
 704
 705 // ************************************************************************* //