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