applications/utilities/surface/surfaceCheck/surfaceCheck.C

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