applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C

   1 #include "checkGeometry.H"
   2 #include "polyMesh.H"
   3 #include "cellSet.H"
   4 #include "faceSet.H"
   5 #include "pointSet.H"
   6 #include "EdgeMap.H"
   7 #include "wedgePolyPatch.H"
   8 #include "unitConversion.H"
   9 #include "polyMeshTetDecomposition.H"
  10
  11
  12 // Find wedge with opposite orientation. Note: does not actually check that
  13 // it is opposite, only that it has opposite normal and same axis
  14 Foam::label Foam::findOppositeWedge
  15 (
  16     const polyMesh& mesh,
  17     const wedgePolyPatch& wpp
  18 )
  19 {
  20     const polyBoundaryMesh& patches = mesh.boundaryMesh();
  21
  22     scalar wppCosAngle = wpp.centreNormal()&wpp.patchNormal();
  23
  24     forAll(patches, patchI)
  25     {
  26         if
  27         (
  28             patchI != wpp.index()
  29          && patches[patchI].size()
  30          && isA<wedgePolyPatch>(patches[patchI])
  31         )
  32         {
  33             const wedgePolyPatch& pp = refCast<const wedgePolyPatch>
  34             (
  35                 patches[patchI]
  36             );
  37
  38             // Calculate (cos of) angle to wpp (not pp!) centre normal
  39             scalar ppCosAngle = wpp.centreNormal()&pp.patchNormal();
  40
  41             if
  42             (
  43                 pp.size() == wpp.size()
  44              && mag(pp.axis() & wpp.axis()) >= (1-1E-3)
  45              && mag(ppCosAngle - wppCosAngle) >= 1E-3
  46             )
  47             {
  48                 return patchI;
  49             }
  50         }
  51     }
  52     return -1;
  53 }
  54
  55
  56 bool Foam::checkWedges
  57 (
  58     const polyMesh& mesh,
  59     const bool report,
  60     const Vector<label>& directions,
  61     labelHashSet* setPtr
  62 )
  63 {
  64     // To mark edges without calculating edge addressing
  65     EdgeMap<label> edgesInError;
  66
  67     const pointField& p = mesh.points();
  68     const faceList& fcs = mesh.faces();
  69
  70
  71     const polyBoundaryMesh& patches = mesh.boundaryMesh();
  72     forAll(patches, patchI)
  73     {
  74         if (patches[patchI].size() && isA<wedgePolyPatch>(patches[patchI]))
  75         {
  76             const wedgePolyPatch& pp = refCast<const wedgePolyPatch>
  77             (
  78                 patches[patchI]
  79             );
  80
  81             scalar wedgeAngle = acos(pp.centreNormal()&pp.patchNormal());
  82
  83             if (report)
  84             {
  85                 Info<< "    Wedge " << pp.name() << " with angle "
  86                     << radToDeg(wedgeAngle) << " degrees"
  87                     << endl;
  88             }
  89
  90             // Find opposite
  91             label oppositePatchI = findOppositeWedge(mesh, pp);
  92
  93             if (oppositePatchI == -1)
  94             {
  95                 if (report)
  96                 {
  97                     Info<< " ***Cannot find opposite wedge for wedge "
  98                         << pp.name() << endl;
  99                 }
 100                 return true;
 101             }
 102
 103             const wedgePolyPatch& opp = refCast<const wedgePolyPatch>
 104             (
 105                 patches[oppositePatchI]
 106             );
 107
 108
 109             if (mag(opp.axis() & pp.axis()) < (1-1E-3))
 110             {
 111                 if (report)
 112                 {
 113                     Info<< " ***Wedges do not have the same axis."
 114                         << " Encountered " << pp.axis()
 115                         << " on patch " << pp.name()
 116                         << " which differs from " << opp.axis()
 117                         << " on opposite wedge patch" << opp.axis()
 118                         << endl;
 119                 }
 120                 return true;
 121             }
 122
 123
 124
 125             // Mark edges on wedgePatches
 126             forAll(pp, i)
 127             {
 128                 const face& f = pp[i];
 129                 forAll(f, fp)
 130                 {
 131                     label p0 = f[fp];
 132                     label p1 = f.nextLabel(fp);
 133                     edgesInError.insert(edge(p0, p1), -1);  // non-error value
 134                 }
 135             }
 136
 137
 138             // Check that wedge patch is flat
 139             const point& p0 = p[pp.meshPoints()[0]];
 140             forAll(pp.meshPoints(), i)
 141             {
 142                 const point& pt = p[pp.meshPoints()[i]];
 143                 scalar d = mag((pt-p0) & pp.patchNormal());
 144
 145                 if (d > sqrt(SMALL))
 146                 {
 147                     if (report)
 148                     {
 149                         Info<< " ***Wedge patch " << pp.name() << " not planar."
 150                             << " Point " << pt << " is not in patch plane by "
 151                             << d << " meter."
 152                             << endl;
 153                     }
 154                     return true;
 155                 }
 156             }
 157         }
 158     }
 159
 160
 161
 162     // Check all non-wedge faces
 163     label nEdgesInError = 0;
 164
 165     forAll(fcs, faceI)
 166     {
 167         const face& f = fcs[faceI];
 168
 169         forAll(f, fp)
 170         {
 171             label p0 = f[fp];
 172             label p1 = f.nextLabel(fp);
 173             if (p0 < p1)
 174             {
 175                 vector d(p[p1]-p[p0]);
 176                 scalar magD = mag(d);
 177
 178                 if (magD > ROOTVSMALL)
 179                 {
 180                     d /= magD;
 181
 182                     // Check how many empty directions are used by the edge.
 183                     label nEmptyDirs = 0;
 184                     label nNonEmptyDirs = 0;
 185                     for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
 186                     {
 187                         if (mag(d[cmpt]) > 1e-6)
 188                         {
 189                             if (directions[cmpt] == 0)
 190                             {
 191                                 nEmptyDirs++;
 192                             }
 193                             else
 194                             {
 195                                 nNonEmptyDirs++;
 196                             }
 197                         }
 198                     }
 199
 200                     if (nEmptyDirs == 0)
 201                     {
 202                         // Purely in ok directions.
 203                     }
 204                     else if (nEmptyDirs == 1)
 205                     {
 206                         // Ok if purely in empty directions.
 207                         if (nNonEmptyDirs > 0)
 208                         {
 209                             if (edgesInError.insert(edge(p0, p1), faceI))
 210                             {
 211                                 nEdgesInError++;
 212                             }
 213                         }
 214                     }
 215                     else if (nEmptyDirs > 1)
 216                     {
 217                         // Always an error
 218                         if (edgesInError.insert(edge(p0, p1), faceI))
 219                         {
 220                             nEdgesInError++;
 221                         }
 222                     }
 223                 }
 224             }
 225         }
 226     }
 227
 228     label nErrorEdges = returnReduce(nEdgesInError, sumOp<label>());
 229
 230     if (nErrorEdges > 0)
 231     {
 232         if (report)
 233         {
 234             Info<< " ***Number of edges not aligned with or perpendicular to "
 235                 << "non-empty directions: " << nErrorEdges << endl;
 236         }
 237
 238         if (setPtr)
 239         {
 240             setPtr->resize(2*nEdgesInError);
 241             forAllConstIter(EdgeMap<label>, edgesInError, iter)
 242             {
 243                 if (iter() >= 0)
 244                 {
 245                     setPtr->insert(iter.key()[0]);
 246                     setPtr->insert(iter.key()[1]);
 247                 }
 248             }
 249         }
 250
 251         return true;
 252     }
 253     else
 254     {
 255         if (report)
 256         {
 257             Info<< "    All edges aligned with or perpendicular to "
 258                 << "non-empty directions." << endl;
 259         }
 260         return false;
 261     }
 262 }
 263
 264
 265 bool Foam::checkCoupledPoints
 266 (
 267     const polyMesh& mesh,
 268     const bool report,
 269     labelHashSet* setPtr
 270 )
 271 {
 272     const pointField& p = mesh.points();
 273     const faceList& fcs = mesh.faces();
 274     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 275
 276     // Zero'th point on coupled faces
 277     pointField nbrZeroPoint(fcs.size()-mesh.nInternalFaces(), vector::max);
 278
 279     // Exchange zero point
 280     forAll(patches, patchI)
 281     {
 282         if (patches[patchI].coupled())
 283         {
 284             const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
 285             (
 286                 patches[patchI]
 287             );
 288
 289             forAll(cpp, i)
 290             {
 291                 label bFaceI = cpp.start()+i-mesh.nInternalFaces();
 292                 const point& p0 = p[cpp[i][0]];
 293                 nbrZeroPoint[bFaceI] = p0;
 294             }
 295         }
 296     }
 297     syncTools::swapBoundaryFacePositions(mesh, nbrZeroPoint);
 298
 299     // Compare to local ones. Use same tolerance as for matching
 300     label nErrorFaces = 0;
 301     scalar avgMismatch = 0;
 302     label nCoupledFaces = 0;
 303
 304     forAll(patches, patchI)
 305     {
 306         if (patches[patchI].coupled())
 307         {
 308             const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
 309             (
 310                 patches[patchI]
 311             );
 312
 313             if (cpp.owner())
 314             {
 315                 scalarField smallDist
 316                 (
 317                     cpp.calcFaceTol
 318                     (
 319                         cpp.matchTolerance(),
 320                         cpp,
 321                         cpp.points(),
 322                         cpp.faceCentres()
 323                     )
 324                 );
 325
 326                 forAll(cpp, i)
 327                 {
 328                     label bFaceI = cpp.start()+i-mesh.nInternalFaces();
 329                     const point& p0 = p[cpp[i][0]];
 330
 331                     scalar d = mag(p0 - nbrZeroPoint[bFaceI]);
 332
 333                     if (d > smallDist[i])
 334                     {
 335                         if (setPtr)
 336                         {
 337                             setPtr->insert(cpp.start()+i);
 338                         }
 339                         nErrorFaces++;
 340                     }
 341                     avgMismatch += d;
 342                     nCoupledFaces++;
 343                 }
 344             }
 345         }
 346     }
 347
 348     reduce(nErrorFaces, sumOp<label>());
 349     reduce(avgMismatch, maxOp<scalar>());
 350     reduce(nCoupledFaces, sumOp<label>());
 351
 352     if (nCoupledFaces > 0)
 353     {
 354         avgMismatch /= nCoupledFaces;
 355     }
 356
 357     if (nErrorFaces > 0)
 358     {
 359         if (report)
 360         {
 361             Info<< "  **Error in coupled point location: "
 362                 << nErrorFaces
 363                 << " faces have their 0th vertex not opposite"
 364                 << " their coupled equivalent. Average mismatch "
 365                 << avgMismatch << "."
 366                 << endl;
 367         }
 368
 369         return true;
 370     }
 371     else
 372     {
 373         if (report)
 374         {
 375             Info<< "    Coupled point location match (average "
 376                 << avgMismatch << ") OK." << endl;
 377         }
 378
 379         return false;
 380     }
 381 }
 382
 383
 384 Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
 385 {
 386     label noFailedChecks = 0;
 387
 388     Info<< "\nChecking geometry..." << endl;
 389
 390     // Get a small relative length from the bounding box
 391     const boundBox& globalBb = mesh.bounds();
 392
 393     Info<< "    Overall domain bounding box "
 394         << globalBb.min() << " " << globalBb.max() << endl;
 395
 396
 397     // Min length
 398     scalar minDistSqr = magSqr(1e-6 * globalBb.span());
 399
 400     // Non-empty directions
 401     const Vector<label> validDirs = (mesh.geometricD() + Vector<label>::one)/2;
 402     Info<< "    Mesh (non-empty, non-wedge) directions " << validDirs << endl;
 403
 404     const Vector<label> solDirs = (mesh.solutionD() + Vector<label>::one)/2;
 405     Info<< "    Mesh (non-empty) directions " << solDirs << endl;
 406
 407     if (mesh.nGeometricD() < 3)
 408     {
 409         pointSet nonAlignedPoints(mesh, "nonAlignedEdges", mesh.nPoints()/100);
 410
 411         if
 412         (
 413             (
 414                 validDirs != solDirs
 415              && checkWedges(mesh, true, validDirs, &nonAlignedPoints)
 416             )
 417          || (
 418                 validDirs == solDirs
 419              && mesh.checkEdgeAlignment(true, validDirs, &nonAlignedPoints)
 420             )
 421         )
 422         {
 423             noFailedChecks++;
 424             label nNonAligned = returnReduce
 425             (
 426                 nonAlignedPoints.size(),
 427                 sumOp<label>()
 428             );
 429
 430             if (nNonAligned > 0)
 431             {
 432                 Info<< "  <<Writing " << nNonAligned
 433                     << " points on non-aligned edges to set "
 434                     << nonAlignedPoints.name() << endl;
 435                 nonAlignedPoints.instance() = mesh.pointsInstance();
 436                 nonAlignedPoints.write();
 437             }
 438         }
 439     }
 440
 441     if (mesh.checkClosedBoundary(true)) noFailedChecks++;
 442
 443     {
 444         cellSet cells(mesh, "nonClosedCells", mesh.nCells()/100+1);
 445         cellSet aspectCells(mesh, "highAspectRatioCells", mesh.nCells()/100+1);
 446         if
 447         (
 448             mesh.checkClosedCells
 449             (
 450                 true,
 451                 &cells,
 452                 &aspectCells,
 453                 mesh.geometricD()
 454             )
 455         )
 456         {
 457             noFailedChecks++;
 458
 459             label nNonClosed = returnReduce(cells.size(), sumOp<label>());
 460
 461             if (nNonClosed > 0)
 462             {
 463                 Info<< "  <<Writing " << nNonClosed
 464                     << " non closed cells to set " << cells.name() << endl;
 465                 cells.instance() = mesh.pointsInstance();
 466                 cells.write();
 467             }
 468         }
 469
 470         label nHighAspect = returnReduce(aspectCells.size(), sumOp<label>());
 471
 472         if (nHighAspect > 0)
 473         {
 474             Info<< "  <<Writing " << nHighAspect
 475                 << " cells with high aspect ratio to set "
 476                 << aspectCells.name() << endl;
 477             aspectCells.instance() = mesh.pointsInstance();
 478             aspectCells.write();
 479         }
 480     }
 481
 482     {
 483         faceSet faces(mesh, "zeroAreaFaces", mesh.nFaces()/100+1);
 484         if (mesh.checkFaceAreas(true, &faces))
 485         {
 486             noFailedChecks++;
 487
 488             label nFaces = returnReduce(faces.size(), sumOp<label>());
 489
 490             if (nFaces > 0)
 491             {
 492                 Info<< "  <<Writing " << nFaces
 493                     << " zero area faces to set " << faces.name() << endl;
 494                 faces.instance() = mesh.pointsInstance();
 495                 faces.write();
 496             }
 497         }
 498     }
 499
 500     {
 501         cellSet cells(mesh, "zeroVolumeCells", mesh.nCells()/100+1);
 502         if (mesh.checkCellVolumes(true, &cells))
 503         {
 504             noFailedChecks++;
 505
 506             label nCells = returnReduce(cells.size(), sumOp<label>());
 507
 508             if (nCells > 0)
 509             {
 510                 Info<< "  <<Writing " << nCells
 511                     << " zero volume cells to set " << cells.name() << endl;
 512                 cells.instance() = mesh.pointsInstance();
 513                 cells.write();
 514             }
 515         }
 516     }
 517
 518     {
 519         faceSet faces(mesh, "nonOrthoFaces", mesh.nFaces()/100+1);
 520         if (mesh.checkFaceOrthogonality(true, &faces))
 521         {
 522             noFailedChecks++;
 523         }
 524
 525         label nFaces = returnReduce(faces.size(), sumOp<label>());
 526
 527         if (nFaces > 0)
 528         {
 529             Info<< "  <<Writing " << nFaces
 530                 << " non-orthogonal faces to set " << faces.name() << endl;
 531             faces.instance() = mesh.pointsInstance();
 532             faces.write();
 533         }
 534     }
 535
 536     {
 537         faceSet faces(mesh, "wrongOrientedFaces", mesh.nFaces()/100 + 1);
 538         if (mesh.checkFacePyramids(true, -SMALL, &faces))
 539         {
 540             noFailedChecks++;
 541
 542             label nFaces = returnReduce(faces.size(), sumOp<label>());
 543
 544             if (nFaces > 0)
 545             {
 546                 Info<< "  <<Writing " << nFaces
 547                     << " faces with incorrect orientation to set "
 548                     << faces.name() << endl;
 549                 faces.instance() = mesh.pointsInstance();
 550                 faces.write();
 551             }
 552         }
 553     }
 554
 555     {
 556         faceSet faces(mesh, "skewFaces", mesh.nFaces()/100+1);
 557         if (mesh.checkFaceSkewness(true, &faces))
 558         {
 559             noFailedChecks++;
 560
 561             label nFaces = returnReduce(faces.size(), sumOp<label>());
 562
 563             if (nFaces > 0)
 564             {
 565                 Info<< "  <<Writing " << nFaces
 566                     << " skew faces to set " << faces.name() << endl;
 567                 faces.instance() = mesh.pointsInstance();
 568                 faces.write();
 569             }
 570         }
 571     }
 572
 573     {
 574         faceSet faces(mesh, "coupledFaces", mesh.nFaces()/100 + 1);
 575         if (checkCoupledPoints(mesh, true, &faces))
 576         {
 577             noFailedChecks++;
 578
 579             label nFaces = returnReduce(faces.size(), sumOp<label>());
 580
 581             if (nFaces > 0)
 582             {
 583                 Info<< "  <<Writing " << nFaces
 584                     << " faces with incorrectly matched 0th vertex to set "
 585                     << faces.name() << endl;
 586                 faces.instance() = mesh.pointsInstance();
 587                 faces.write();
 588             }
 589         }
 590     }
 591
 592     if (allGeometry)
 593     {
 594         faceSet faces(mesh, "lowQualityTetFaces", mesh.nFaces()/100+1);
 595         if
 596         (
 597             polyMeshTetDecomposition::checkFaceTets
 598             (
 599                 mesh,
 600                 polyMeshTetDecomposition::minTetQuality,
 601                 true,
 602                 &faces
 603             )
 604         )
 605         {
 606             noFailedChecks++;
 607
 608             label nFaces = returnReduce(faces.size(), sumOp<label>());
 609
 610             if (nFaces > 0)
 611             {
 612                 Info<< "  <<Writing " << nFaces
 613                     << " faces with low quality or negative volume "
 614                     << "decomposition tets to set " << faces.name() << endl;
 615                 faces.instance() = mesh.pointsInstance();
 616                 faces.write();
 617             }
 618         }
 619     }
 620
 621     if (allGeometry)
 622     {
 623         // Note use of nPoints since don't want edge construction.
 624         pointSet points(mesh, "shortEdges", mesh.nPoints()/1000 + 1);
 625         if (mesh.checkEdgeLength(true, minDistSqr, &points))
 626         {
 627             //noFailedChecks++;
 628
 629             label nPoints = returnReduce(points.size(), sumOp<label>());
 630
 631             if (nPoints > 0)
 632             {
 633                 Info<< "  <<Writing " << nPoints
 634                     << " points on short edges to set " << points.name()
 635                     << endl;
 636                 points.instance() = mesh.pointsInstance();
 637                 points.write();
 638             }
 639         }
 640
 641         label nEdgeClose = returnReduce(points.size(), sumOp<label>());
 642
 643         if (mesh.checkPointNearness(false, minDistSqr, &points))
 644         {
 645             //noFailedChecks++;
 646
 647             label nPoints = returnReduce(points.size(), sumOp<label>());
 648
 649             if (nPoints > nEdgeClose)
 650             {
 651                 pointSet nearPoints(mesh, "nearPoints", points);
 652                 Info<< "  <<Writing " << nPoints
 653                     << " near (closer than " << Foam::sqrt(minDistSqr)
 654                     << " apart) points to set " << nearPoints.name() << endl;
 655                 nearPoints.instance() = mesh.pointsInstance();
 656                 nearPoints.write();
 657             }
 658         }
 659     }
 660
 661     if (allGeometry)
 662     {
 663         faceSet faces(mesh, "concaveFaces", mesh.nFaces()/100 + 1);
 664         if (mesh.checkFaceAngles(true, 10, &faces))
 665         {
 666             //noFailedChecks++;
 667
 668             label nFaces = returnReduce(faces.size(), sumOp<label>());
 669
 670             if (nFaces > 0)
 671             {
 672                 Info<< "  <<Writing " << nFaces
 673                     << " faces with concave angles to set " << faces.name()
 674                     << endl;
 675                 faces.instance() = mesh.pointsInstance();
 676                 faces.write();
 677             }
 678         }
 679     }
 680
 681     if (allGeometry)
 682     {
 683         faceSet faces(mesh, "warpedFaces", mesh.nFaces()/100 + 1);
 684         if (mesh.checkFaceFlatness(true, 0.8, &faces))
 685         {
 686             //noFailedChecks++;
 687
 688             label nFaces = returnReduce(faces.size(), sumOp<label>());
 689
 690             if (nFaces > 0)
 691             {
 692                 Info<< "  <<Writing " << nFaces
 693                     << " warped faces to set " << faces.name() << endl;
 694                 faces.instance() = mesh.pointsInstance();
 695                 faces.write();
 696             }
 697         }
 698     }
 699
 700     if (allGeometry)
 701     {
 702         cellSet cells(mesh, "underdeterminedCells", mesh.nCells()/100);
 703         if (mesh.checkCellDeterminant(true, &cells, mesh.geometricD()))
 704         {
 705             noFailedChecks++;
 706
 707             label nCells = returnReduce(cells.size(), sumOp<label>());
 708
 709             Info<< "  <<Writing " << nCells
 710                 << " under-determined cells to set " << cells.name() << endl;
 711             cells.instance() = mesh.pointsInstance();
 712             cells.write();
 713         }
 714     }
 715
 716     if (allGeometry)
 717     {
 718         cellSet cells(mesh, "concaveCells", mesh.nCells()/100);
 719         if (mesh.checkConcaveCells(true, &cells))
 720         {
 721             noFailedChecks++;
 722
 723             label nCells = returnReduce(cells.size(), sumOp<label>());
 724
 725             Info<< "  <<Writing " << nCells
 726                 << " concave cells to set " << cells.name() << endl;
 727             cells.instance() = mesh.pointsInstance();
 728             cells.write();
 729         }
 730     }
 731
 732     return noFailedChecks;
 733 }