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