applications/utilities/surface/surfaceCheck/surfaceCheck.C

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