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