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 \*---------------------------------------------------------------------------*/
  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
 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
 203         forAll(surf, faceI)
 204         {
 205             label region = surf[faceI].region();
 206
 207             if (region < 0 || region >= regionSize.size())
 208             {
 209                 WarningIn(args.executable())
 210                     << "Triangle " << faceI << " vertices " << surf[faceI]
 211                     << " has region " << region << " which is outside the range"
 212                     << " of regions 0.." << surf.patches().size()-1
 213                     << endl;
 214             }
 215             else
 216             {
 217                 regionSize[region]++;
 218             }
 219         }
 220
 221         Pout<< "Region\tSize" << nl
 222             << "------\t----" << nl;
 223         forAll(surf.patches(), patchI)
 224         {
 225             Pout<< surf.patches()[patchI].name() << '\t'
 226                 << regionSize[patchI] << nl;
 227         }
 228         Pout<< nl << endl;
 229     }
 230
 231
 232     // Check triangles
 233     // ~~~~~~~~~~~~~~~
 234
 235     {
 236         DynamicList<label> illegalFaces(surf.size()/100 + 1);
 237
 238         forAll(surf, faceI)
 239         {
 240             if (!validTri(verbose, surf, faceI))
 241             {
 242                 illegalFaces.append(faceI);
 243             }
 244         }
 245
 246         if (illegalFaces.size())
 247         {
 248             Pout<< "Surface has " << illegalFaces.size()
 249                 << " illegal triangles." << endl;
 250
 251             OFstream str("illegalFaces");
 252             Pout<< "Dumping conflicting face labels to " << str.name() << endl
 253                 << "Paste this into the input for surfaceSubset" << endl;
 254             str << illegalFaces;
 255         }
 256         else
 257         {
 258             Pout<< "Surface has no illegal triangles." << endl;
 259         }
 260         Pout<< endl;
 261     }
 262
 263
 264
 265     // Triangle quality
 266     // ~~~~~~~~~~~~~~~~
 267
 268     {
 269         scalarField triQ(surf.size(), 0);
 270         forAll(surf, faceI)
 271         {
 272             const labelledTri& f = surf[faceI];
 273
 274             if (f[0] == f[1] || f[0] == f[2] || f[1] == f[2])
 275             {
 276                 //WarningIn(args.executable())
 277                 //    << "Illegal triangle " << faceI << " vertices " << f
 278                 //    << " coords " << f.points(surf.points()) << endl;
 279             }
 280             else
 281             {
 282                 triPointRef tri
 283                 (
 284                     surf.points()[f[0]],
 285                     surf.points()[f[1]],
 286                     surf.points()[f[2]]
 287                 );
 288
 289                 vector ba(tri.b() - tri.a());
 290                 ba /= mag(ba) + VSMALL;
 291
 292                 vector ca(tri.c() - tri.a());
 293                 ca /= mag(ca) + VSMALL;
 294
 295                 if (mag(ba&ca) > 1-1E-3)
 296                 {
 297                     triQ[faceI] = SMALL;
 298                 }
 299                 else
 300                 {
 301                     triQ[faceI] = triPointRef
 302                     (
 303                         surf.points()[f[0]],
 304                         surf.points()[f[1]],
 305                         surf.points()[f[2]]
 306                     ).quality();
 307                 }
 308             }
 309         }
 310
 311         labelList binCount = countBins(0, 1, 20, triQ);
 312
 313         Pout<< "Triangle quality (equilateral=1, collapsed=0):"
 314             << endl;
 315
 316
 317         OSstream& os = Pout;
 318         os.width(4);
 319
 320         scalar dist = (1.0 - 0.0)/20.0;
 321         scalar min = 0;
 322         forAll(binCount, binI)
 323         {
 324             Pout<< "    " << min << " .. " << min+dist << "  : "
 325                 << 1.0/surf.size() * binCount[binI]
 326                 << endl;
 327             min += dist;
 328         }
 329         Pout<< endl;
 330
 331         label minIndex = findMin(triQ);
 332         label maxIndex = findMax(triQ);
 333
 334         Pout<< "    min " << triQ[minIndex] << " for triangle " << minIndex
 335             << nl
 336             << "    max " << triQ[maxIndex] << " for triangle " << maxIndex
 337             << nl
 338             << endl;
 339
 340
 341         if (triQ[minIndex] < SMALL)
 342         {
 343             WarningIn(args.executable()) << "Minimum triangle quality is "
 344                 << triQ[minIndex] << ". This might give problems in"
 345                 << " self-intersection testing later on." << endl;
 346         }
 347
 348         // Dump for subsetting
 349         {
 350             DynamicList<label> problemFaces(surf.size()/100+1);
 351
 352             forAll(triQ, faceI)
 353             {
 354                 if (triQ[faceI] < 1E-11)
 355                 {
 356                     problemFaces.append(faceI);
 357                 }
 358             }
 359             OFstream str("badFaces");
 360
 361             Pout<< "Dumping bad quality faces to " << str.name() << endl
 362                 << "Paste this into the input for surfaceSubset" << nl
 363                 << nl << endl;
 364
 365             str << problemFaces;
 366         }
 367     }
 368
 369
 370
 371     // Edges
 372     // ~~~~~
 373     {
 374         const edgeList& edges = surf.edges();
 375         const pointField& localPoints = surf.localPoints();
 376
 377         scalarField edgeMag(edges.size());
 378
 379         forAll(edges, edgeI)
 380         {
 381             edgeMag[edgeI] = edges[edgeI].mag(localPoints);
 382         }
 383
 384         label minEdgeI = findMin(edgeMag);
 385         label maxEdgeI = findMax(edgeMag);
 386
 387         const edge& minE = edges[minEdgeI];
 388         const edge& maxE = edges[maxEdgeI];
 389
 390
 391         Pout<< "Edges:" << nl
 392             << "    min " << edgeMag[minEdgeI] << " for edge " << minEdgeI
 393             << " points " << localPoints[minE[0]] << localPoints[minE[1]]
 394             << nl
 395             << "    max " << edgeMag[maxEdgeI] << " for edge " << maxEdgeI
 396             << " points " << localPoints[maxE[0]] << localPoints[maxE[1]]
 397             << nl
 398             << endl;
 399     }
 400
 401
 402
 403     // Close points
 404     // ~~~~~~~~~~~~
 405     {
 406         const edgeList& edges = surf.edges();
 407         const pointField& localPoints = surf.localPoints();
 408
 409         const boundBox bb(localPoints);
 410         scalar smallDim = 1E-6 * bb.mag();
 411
 412         Pout<< "Checking for points less than 1E-6 of bounding box ("
 413             << bb.span() << " meter) apart."
 414             << endl;
 415
 416         // Sort points
 417         SortableList<scalar> sortedMag(mag(localPoints));
 418
 419         label nClose = 0;
 420
 421         for (label i = 1; i < sortedMag.size(); i++)
 422         {
 423             label ptI = sortedMag.indices()[i];
 424
 425             label prevPtI = sortedMag.indices()[i-1];
 426
 427             if (mag(localPoints[ptI] - localPoints[prevPtI]) < smallDim)
 428             {
 429                 // Check if neighbours.
 430                 const labelList& pEdges = surf.pointEdges()[ptI];
 431
 432                 label edgeI = -1;
 433
 434                 forAll(pEdges, i)
 435                 {
 436                     const edge& e = edges[pEdges[i]];
 437
 438                     if (e[0] == prevPtI || e[1] == prevPtI)
 439                     {
 440                         // point1 and point0 are connected through edge.
 441                         edgeI = pEdges[i];
 442
 443                         break;
 444                     }
 445                 }
 446
 447                 nClose++;
 448
 449                 if (edgeI == -1)
 450                 {
 451                     Pout<< "    close unconnected points "
 452                         << ptI << ' ' << localPoints[ptI]
 453                         << " and " << prevPtI << ' '
 454                         << localPoints[prevPtI]
 455                         << " distance:"
 456                         << mag(localPoints[ptI] - localPoints[prevPtI])
 457                         << endl;
 458                 }
 459                 else
 460                 {
 461                     Pout<< "    small edge between points "
 462                         << ptI << ' ' << localPoints[ptI]
 463                         << " and " << prevPtI << ' '
 464                         << localPoints[prevPtI]
 465                         << " distance:"
 466                         << mag(localPoints[ptI] - localPoints[prevPtI])
 467                         << endl;
 468                 }
 469             }
 470         }
 471
 472         Pout<< "Found " << nClose << " nearby points." << nl
 473             << endl;
 474     }
 475
 476
 477
 478     // Check manifold
 479     // ~~~~~~~~~~~~~~
 480
 481     DynamicList<label> problemFaces(surf.size()/100 + 1);
 482
 483     const labelListList& eFaces = surf.edgeFaces();
 484
 485     label nSingleEdges = 0;
 486     forAll(eFaces, edgeI)
 487     {
 488         const labelList& myFaces = eFaces[edgeI];
 489
 490         if (myFaces.size() == 1)
 491         {
 492             problemFaces.append(myFaces[0]);
 493
 494             nSingleEdges++;
 495         }
 496     }
 497
 498     label nMultEdges = 0;
 499     forAll(eFaces, edgeI)
 500     {
 501         const labelList& myFaces = eFaces[edgeI];
 502
 503         if (myFaces.size() > 2)
 504         {
 505             forAll(myFaces, myFaceI)
 506             {
 507                 problemFaces.append(myFaces[myFaceI]);
 508             }
 509
 510             nMultEdges++;
 511         }
 512     }
 513     problemFaces.shrink();
 514
 515     if ((nSingleEdges != 0) || (nMultEdges != 0))
 516     {
 517         Pout<< "Surface is not closed since not all edges connected to "
 518             << "two faces:" << endl
 519             << "    connected to one face : " << nSingleEdges << endl
 520             << "    connected to >2 faces : " << nMultEdges << endl;
 521
 522         Pout<< "Conflicting face labels:" << problemFaces.size() << endl;
 523
 524         OFstream str("problemFaces");
 525
 526         Pout<< "Dumping conflicting face labels to " << str.name() << endl
 527             << "Paste this into the input for surfaceSubset" << endl;
 528
 529         str << problemFaces;
 530     }
 531     else
 532     {
 533         Pout<< "Surface is closed. All edges connected to two faces." << endl;
 534     }
 535     Pout<< endl;
 536
 537
 538
 539     // Check singly connected domain
 540     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 541
 542     labelList faceZone;
 543     label numZones = surf.markZones(boolList(surf.nEdges(), false), faceZone);
 544
 545     Pout<< "Number of unconnected parts : " << numZones << endl;
 546
 547     if (numZones > 1)
 548     {
 549         Pout<< "Splitting surface into parts ..." << endl << endl;
 550
 551         fileName surfFileNameBase(surfFileName.name());
 552
 553         for(label zone = 0; zone < numZones; zone++)
 554         {
 555             boolList includeMap(surf.size(), false);
 556
 557             forAll(faceZone, faceI)
 558             {
 559                 if (faceZone[faceI] == zone)
 560                 {
 561                     includeMap[faceI] = true;
 562                 }
 563             }
 564
 565             labelList pointMap;
 566             labelList faceMap;
 567
 568             triSurface subSurf
 569             (
 570                 surf.subsetMesh
 571                 (
 572                     includeMap,
 573                     pointMap,
 574                     faceMap
 575                 )
 576             );
 577
 578             fileName subFileName
 579             (
 580                 surfFileNameBase.lessExt()
 581               + "_"
 582               + name(zone)
 583               + ".ftr"
 584             );
 585
 586             Pout<< "writing part " << zone << " size " << subSurf.size()
 587                 << " to " << subFileName << endl;
 588
 589             subSurf.write(subFileName);
 590         }
 591
 592         return 0;
 593     }
 594
 595
 596
 597     // Check orientation
 598     // ~~~~~~~~~~~~~~~~~
 599
 600     labelHashSet borderEdge(surf.size()/1000);
 601     PatchTools::checkOrientation(surf, false, &borderEdge);
 602
 603     //
 604     // Colour all faces into zones using borderEdge
 605     //
 606     labelList normalZone;
 607     label numNormalZones = PatchTools::markZones(surf, borderEdge, normalZone);
 608
 609     Pout<< endl
 610         << "Number of zones (connected area with consistent normal) : "
 611         << numNormalZones << endl;
 612
 613     if (numNormalZones > 1)
 614     {
 615         Pout<< "More than one normal orientation." << endl;
 616     }
 617     Pout<< endl;
 618
 619
 620
 621     // Check self-intersection
 622     // ~~~~~~~~~~~~~~~~~~~~~~~
 623
 624     if (checkSelfIntersection)
 625     {
 626         Pout<< "Checking self-intersection." << endl;
 627
 628         triSurfaceSearch querySurf(surf);
 629         surfaceIntersection inter(querySurf);
 630
 631         if (inter.cutEdges().empty() && inter.cutPoints().empty())
 632         {
 633             Pout<< "Surface is not self-intersecting" << endl;
 634         }
 635         else
 636         {
 637             Pout<< "Surface is self-intersecting" << endl;
 638             Pout<< "Writing edges of intersection to selfInter.obj" << endl;
 639
 640             OFstream intStream("selfInter.obj");
 641             forAll(inter.cutPoints(), cutPointI)
 642             {
 643                 const point& pt = inter.cutPoints()[cutPointI];
 644
 645                 intStream << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z()
 646                     << endl;
 647             }
 648             forAll(inter.cutEdges(), cutEdgeI)
 649             {
 650                 const edge& e = inter.cutEdges()[cutEdgeI];
 651
 652                 intStream << "l " << e.start()+1 << ' ' << e.end()+1 << endl;
 653             }
 654         }
 655         Pout<< endl;
 656     }
 657
 658
 659     Pout<< "End\n" << endl;
 660
 661     return 0;
 662 }
 663
 664
 665 // ************************************************************************* //