primitiveMeshCheck.C 42.9 KB
Newer Older
1 2 3 4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
OpenFOAM bot's avatar
OpenFOAM bot committed
5
    \\  /    A nd           | www.openfoam.com
6
     \\/     M anipulation  |
OpenFOAM bot's avatar
OpenFOAM bot committed
7
-------------------------------------------------------------------------------
OpenFOAM bot's avatar
OpenFOAM bot committed
8
    Copyright (C) 2011-2016 OpenFOAM Foundation
Mark Olesen's avatar
Mark Olesen committed
9
    Copyright (C) 2019-2020 OpenCFD Ltd.
10 11 12 13
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

mattijs's avatar
mattijs committed
14 15 16 17
    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
18 19 20 21 22 23 24

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
mattijs's avatar
mattijs committed
25
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
26 27 28 29 30 31

\*---------------------------------------------------------------------------*/

#include "primitiveMesh.H"
#include "pyramidPointFaceRef.H"
#include "ListOps.H"
32
#include "unitConversion.H"
33
#include "SortableList.H"
Mark Olesen's avatar
Mark Olesen committed
34
#include "edgeHashes.H"
35
#include "primitiveMeshTools.H"
36 37 38

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

39 40 41 42
Foam::scalar Foam::primitiveMesh::closedThreshold_  = 1.0e-6;
Foam::scalar Foam::primitiveMesh::aspectThreshold_  = 1000;
Foam::scalar Foam::primitiveMesh::nonOrthThreshold_ = 70;    // deg
Foam::scalar Foam::primitiveMesh::skewThreshold_    = 4;
mattijs's avatar
mattijs committed
43
Foam::scalar Foam::primitiveMesh::planarCosAngle_   = 1.0e-6;
44 45 46 47


// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

48 49 50 51
bool Foam::primitiveMesh::checkClosedBoundary
(
    const vectorField& areas,
    const bool report,
52
    const bitSet& internalOrCoupledFaces
53
) const
54
{
Mark Olesen's avatar
Mark Olesen committed
55
    DebugInFunction << "Checking if boundary is closed" << endl;
56 57 58 59

    // Loop through all boundary faces and sum up the face area vectors.
    // For a closed boundary, this should be zero in all vector components

60 61
    Vector<solveScalar> sumClosed(Zero);
    solveScalar sumMagClosedBoundary = 0;
62

63
    for (label facei = nInternalFaces(); facei < areas.size(); facei++)
64
    {
65
        if (!internalOrCoupledFaces.size() || !internalOrCoupledFaces[facei])
66
        {
67
            sumClosed += Vector<solveScalar>(areas[facei]);
68
            sumMagClosedBoundary += mag(areas[facei]);
69
        }
70 71
    }

72 73
    reduce(sumClosed, sumOp<Vector<solveScalar>>());
    reduce(sumMagClosedBoundary, sumOp<solveScalar>());
74

75
    Vector<solveScalar> openness = sumClosed/(sumMagClosedBoundary + VSMALL);
76 77 78 79 80 81 82 83 84 85 86 87 88

    if (cmptMax(cmptMag(openness)) > closedThreshold_)
    {
        if (debug || report)
        {
            Info<< " ***Boundary openness " << openness
                << " possible hole in boundary description."
                << endl;
        }

        return true;
    }

89 90 91 92
    if (debug || report)
    {
        Info<< "    Boundary openness " << openness << " OK."
            << endl;
93
    }
94 95

    return false;
96 97 98
}


99
bool Foam::primitiveMesh::checkClosedCells
100
(
101 102
    const vectorField& faceAreas,
    const scalarField& cellVolumes,
103 104
    const bool report,
    labelHashSet* setPtr,
105 106
    labelHashSet* aspectSetPtr,
    const Vector<label>& meshD
107 108
) const
{
Mark Olesen's avatar
Mark Olesen committed
109
    DebugInFunction << "Checking if cells are closed" << endl;
110 111 112 113 114 115

    // Check that all cells labels are valid
    const cellList& c = cells();

    label nErrorClosed = 0;

116
    forAll(c, cI)
117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142
    {
        const cell& curCell = c[cI];

        if (min(curCell) < 0 || max(curCell) > nFaces())
        {
            if (setPtr)
            {
                setPtr->insert(cI);
            }

            nErrorClosed++;
        }
    }

    if (nErrorClosed > 0)
    {
        if (debug || report)
        {
            Info<< " ***Cells with invalid face labels found, number of cells "
                << nErrorClosed << endl;
        }

        return true;
    }


143 144 145 146 147 148 149 150 151 152 153
    scalarField openness;
    scalarField aspectRatio;
    primitiveMeshTools::cellClosedness
    (
        *this,
        meshD,
        faceAreas,
        cellVolumes,
        openness,
        aspectRatio
    );
154 155

    label nOpen = 0;
156
    scalar maxOpennessCell = max(openness);
157
    label nAspect = 0;
158
    scalar maxAspectRatio = max(aspectRatio);
159

160
    // Check the sums
161
    forAll(openness, celli)
162
    {
163
        if (openness[celli] > closedThreshold_)
164 165 166
        {
            if (setPtr)
            {
167
                setPtr->insert(celli);
168 169 170 171 172
            }

            nOpen++;
        }

173
        if (aspectRatio[celli] > aspectThreshold_)
174 175 176
        {
            if (aspectSetPtr)
            {
177
                aspectSetPtr->insert(celli);
178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226
            }

            nAspect++;
        }
    }

    reduce(nOpen, sumOp<label>());
    reduce(maxOpennessCell, maxOp<scalar>());

    reduce(nAspect, sumOp<label>());
    reduce(maxAspectRatio, maxOp<scalar>());


    if (nOpen > 0)
    {
        if (debug || report)
        {
            Info<< " ***Open cells found, max cell openness: "
                << maxOpennessCell << ", number of open cells " << nOpen
                << endl;
        }

        return true;
    }

    if (nAspect > 0)
    {
        if (debug || report)
        {
            Info<< " ***High aspect ratio cells found, Max aspect ratio: "
                << maxAspectRatio
                << ", number of cells " << nAspect
                << endl;
        }

        return true;
    }

    if (debug || report)
    {
        Info<< "    Max cell openness = " << maxOpennessCell << " OK." << nl
            << "    Max aspect ratio = " << maxAspectRatio << " OK."
            << endl;
    }

    return false;
}


227
bool Foam::primitiveMesh::checkFaceAreas
228
(
229
    const vectorField& faceAreas,
230
    const bool report,
231
    const bool detailedReport,
232 233 234
    labelHashSet* setPtr
) const
{
Mark Olesen's avatar
Mark Olesen committed
235
    DebugInFunction << "Checking face area magnitudes" << endl;
236

237
    const scalarField magFaceAreas(mag(faceAreas));
238 239 240 241

    scalar minArea = GREAT;
    scalar maxArea = -GREAT;

242
    forAll(magFaceAreas, facei)
243
    {
244
        if (magFaceAreas[facei] < VSMALL)
245 246 247
        {
            if (setPtr)
            {
248
                setPtr->insert(facei);
249
            }
250 251
            if (detailedReport)
            {
252
                if (isInternalFace(facei))
253 254
                {
                    Pout<< "Zero or negative face area detected for "
255 256 257 258
                        << "internal face "<< facei << " between cells "
                        << faceOwner()[facei] << " and "
                        << faceNeighbour()[facei]
                        << ".  Face area magnitude = " << magFaceAreas[facei]
259 260 261 262 263
                        << endl;
                }
                else
                {
                    Pout<< "Zero or negative face area detected for "
264 265 266
                        << "boundary face " << facei << " next to cell "
                        << faceOwner()[facei] << ".  Face area magnitude = "
                        << magFaceAreas[facei] << endl;
267 268
                }
            }
269 270
        }

271 272
        minArea = min(minArea, magFaceAreas[facei]);
        maxArea = max(maxArea, magFaceAreas[facei]);
273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288
    }

    reduce(minArea, minOp<scalar>());
    reduce(maxArea, maxOp<scalar>());

    if (minArea < VSMALL)
    {
        if (debug || report)
        {
            Info<< " ***Zero or negative face area detected.  "
                "Minimum area: " << minArea << endl;
        }

        return true;
    }

289 290 291 292 293
    if (debug || report)
    {
        Info<< "    Minimum face area = " << minArea
            << ". Maximum face area = " << maxArea
            << ".  Face area magnitudes OK." << endl;
294
    }
295 296

    return false;
297 298 299
}


300
bool Foam::primitiveMesh::checkCellVolumes
301
(
302
    const scalarField& vols,
303
    const bool report,
304
    const bool detailedReport,
305 306 307
    labelHashSet* setPtr
) const
{
Mark Olesen's avatar
Mark Olesen committed
308
    DebugInFunction << "Checking cell volumes" << endl;
309 310 311 312 313 314

    scalar minVolume = GREAT;
    scalar maxVolume = -GREAT;

    label nNegVolCells = 0;

315
    forAll(vols, celli)
316
    {
317
        if (vols[celli] < VSMALL)
318 319 320
        {
            if (setPtr)
            {
321
                setPtr->insert(celli);
322
            }
323 324 325
            if (detailedReport)
            {
                Pout<< "Zero or negative cell volume detected for cell "
326
                    << celli << ".  Volume = " << vols[celli] << endl;
327
            }
328 329 330 331

            nNegVolCells++;
        }

332 333
        minVolume = min(minVolume, vols[celli]);
        maxVolume = max(maxVolume, vols[celli]);
334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352
    }

    reduce(minVolume, minOp<scalar>());
    reduce(maxVolume, maxOp<scalar>());
    reduce(nNegVolCells, sumOp<label>());

    if (minVolume < VSMALL)
    {
        if (debug || report)
        {
            Info<< " ***Zero or negative cell volume detected.  "
                << "Minimum negative volume: " << minVolume
                << ", Number of negative volume cells: " << nNegVolCells
                << endl;
        }

        return true;
    }

353 354 355 356 357 358
    if (debug || report)
    {
        Info<< "    Min volume = " << minVolume
            << ". Max volume = " << maxVolume
            << ".  Total volume = " << gSum(vols)
            << ".  Cell volumes OK." << endl;
359
    }
360 361

    return false;
362 363 364
}


365
bool Foam::primitiveMesh::checkFaceOrthogonality
366
(
367 368
    const vectorField& fAreas,
    const vectorField& cellCtrs,
369 370 371 372
    const bool report,
    labelHashSet* setPtr
) const
{
Mark Olesen's avatar
Mark Olesen committed
373
    DebugInFunction << "Checking mesh non-orthogonality" << endl;
374

375 376 377 378 379 380 381
    tmp<scalarField> tortho = primitiveMeshTools::faceOrthogonality
    (
        *this,
        fAreas,
        cellCtrs
    );
    const scalarField& ortho = tortho();
382 383 384

    // Severe nonorthogonality threshold
    const scalar severeNonorthogonalityThreshold =
385
        ::cos(degToRad(nonOrthThreshold_));
386

387
    scalar minDDotS = min(ortho);
388

389
    scalar sumDDotS = sum(ortho);
390 391 392 393 394 395

    label severeNonOrth = 0;

    label errorNonOrth = 0;


396
    forAll(ortho, facei)
397
    {
398
        if (ortho[facei] < severeNonorthogonalityThreshold)
399
        {
400
            if (ortho[facei] > SMALL)
401 402 403
            {
                if (setPtr)
                {
404
                    setPtr->insert(facei);
405 406 407 408 409 410 411 412
                }

                severeNonOrth++;
            }
            else
            {
                if (setPtr)
                {
413
                    setPtr->insert(facei);
414 415 416 417 418 419 420 421 422 423 424 425 426 427
                }

                errorNonOrth++;
            }
        }
    }

    reduce(minDDotS, minOp<scalar>());
    reduce(sumDDotS, sumOp<scalar>());
    reduce(severeNonOrth, sumOp<label>());
    reduce(errorNonOrth, sumOp<label>());

    if (debug || report)
    {
428
        label neiSize = ortho.size();
429 430 431 432 433 434 435
        reduce(neiSize, sumOp<label>());

        if (neiSize > 0)
        {
            if (debug || report)
            {
                Info<< "    Mesh non-orthogonality Max: "
436 437
                    << radToDeg(::acos(minDDotS))
                    << " average: " << radToDeg(::acos(sumDDotS/neiSize))
438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459
                    << endl;
            }
        }

        if (severeNonOrth > 0)
        {
            Info<< "   *Number of severely non-orthogonal faces: "
                << severeNonOrth << "." << endl;
        }
    }

    if (errorNonOrth > 0)
    {
        if (debug || report)
        {
            Info<< " ***Number of non-orthogonality errors: "
                << errorNonOrth << "." << endl;
        }

        return true;
    }

460 461 462
    if (debug || report)
    {
        Info<< "    Non-orthogonality check OK." << endl;
463
    }
464 465

    return false;
466 467 468
}


469
bool Foam::primitiveMesh::checkFacePyramids
470
(
471 472
    const pointField& points,
    const vectorField& ctrs,
473
    const bool report,
474
    const bool detailedReport,
475 476 477 478
    const scalar minPyrVol,
    labelHashSet* setPtr
) const
{
Mark Olesen's avatar
Mark Olesen committed
479
    DebugInFunction << "Checking face orientation" << endl;
480 481 482 483 484

    const labelList& own = faceOwner();
    const labelList& nei = faceNeighbour();
    const faceList& f = faces();

485 486 487 488 489 490 491 492 493 494 495 496

    scalarField ownPyrVol;
    scalarField neiPyrVol;
    primitiveMeshTools::facePyramidVolume
    (
        *this,
        points,
        ctrs,
        ownPyrVol,
        neiPyrVol
    );

497 498 499

    label nErrorPyrs = 0;

500
    forAll(ownPyrVol, facei)
501
    {
502
        if (ownPyrVol[facei] < minPyrVol)
503 504 505
        {
            if (setPtr)
            {
506
                setPtr->insert(facei);
507
            }
508 509
            if (detailedReport)
            {
510 511 512
                Pout<< "Negative pyramid volume: " << ownPyrVol[facei]
                    << " for face " << facei << " " << f[facei]
                    << "  and owner cell: " << own[facei] << endl
513
                    << "Owner cell vertex labels: "
514
                    << cells()[own[facei]].labels(faces())
515 516
                    << endl;
            }
517 518 519 520

            nErrorPyrs++;
        }

521
        if (isInternalFace(facei))
522
        {
523
            if (neiPyrVol[facei] < minPyrVol)
524 525 526
            {
                if (setPtr)
                {
527
                    setPtr->insert(facei);
528
                }
529 530
                if (detailedReport)
                {
531 532 533
                    Pout<< "Negative pyramid volume: " << neiPyrVol[facei]
                        << " for face " << facei << " " << f[facei]
                        << "  and neighbour cell: " << nei[facei] << nl
534
                        << "Neighbour cell vertex labels: "
535
                        << cells()[nei[facei]].labels(faces())
536 537
                        << endl;
                }
538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556
                nErrorPyrs++;
            }
        }
    }

    reduce(nErrorPyrs, sumOp<label>());

    if (nErrorPyrs > 0)
    {
        if (debug || report)
        {
            Info<< " ***Error in face pyramids: "
                << nErrorPyrs << " faces are incorrectly oriented."
                << endl;
        }

        return true;
    }

557 558 559
    if (debug || report)
    {
        Info<< "    Face pyramids OK." << endl;
560
    }
561 562

    return false;
563 564 565
}


566
bool Foam::primitiveMesh::checkFaceSkewness
567
(
568 569 570 571
    const pointField& points,
    const vectorField& fCtrs,
    const vectorField& fAreas,
    const vectorField& cellCtrs,
572 573 574 575
    const bool report,
    labelHashSet* setPtr
) const
{
Mark Olesen's avatar
Mark Olesen committed
576
    DebugInFunction << "Checking face skewness" << endl;
577 578 579 580

    // Warn if the skew correction vector is more than skewWarning times
    // larger than the face area vector

581 582 583 584 585 586 587 588 589 590 591
    tmp<scalarField> tskewness = primitiveMeshTools::faceSkewness
    (
        *this,
        points,
        fCtrs,
        fAreas,
        cellCtrs
    );
    const scalarField& skewness = tskewness();

    scalar maxSkew = max(skewness);
592 593
    label nWarnSkew = 0;

594
    forAll(skewness, facei)
595 596 597
    {
        // Check if the skewness vector is greater than the PN vector.
        // This does not cause trouble but is a good indication of a poor mesh.
598
        if (skewness[facei] > skewThreshold_)
599 600 601
        {
            if (setPtr)
            {
602
                setPtr->insert(facei);
603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624
            }

            nWarnSkew++;
        }
    }

    reduce(maxSkew, maxOp<scalar>());
    reduce(nWarnSkew, sumOp<label>());

    if (nWarnSkew > 0)
    {
        if (debug || report)
        {
            Info<< " ***Max skewness = " << maxSkew
                << ", " << nWarnSkew << " highly skew faces detected"
                   " which may impair the quality of the results"
                << endl;
        }

        return true;
    }

625 626 627
    if (debug || report)
    {
        Info<< "    Max skewness = " << maxSkew << " OK." << endl;
628
    }
629 630

    return false;
631 632 633
}


634
bool Foam::primitiveMesh::checkFaceAngles
635
(
636 637
    const pointField& points,
    const vectorField& faceAreas,
638 639 640 641 642
    const bool report,
    const scalar maxDeg,
    labelHashSet* setPtr
) const
{
Mark Olesen's avatar
Mark Olesen committed
643
    DebugInFunction << "Checking face angles" << endl;
644 645 646

    if (maxDeg < -SMALL || maxDeg > 180+SMALL)
    {
647 648
        FatalErrorInFunction
            << "maxDeg should be [0..180] but is now " << maxDeg
649 650 651
            << exit(FatalError);
    }

652
    const scalar maxSin = Foam::sin(degToRad(maxDeg));
653 654


655 656 657 658 659 660 661 662
    tmp<scalarField> tfaceAngles = primitiveMeshTools::faceConcavity
    (
        maxSin,
        *this,
        points,
        faceAreas
    );
    const scalarField& faceAngles = tfaceAngles();
663

664
    scalar maxEdgeSin = max(faceAngles);
665

666
    label nConcave = 0;
667

668
    forAll(faceAngles, facei)
669
    {
670
        if (faceAngles[facei] > SMALL)
671
        {
672
            nConcave++;
673

674
            if (setPtr)
675
            {
676
                setPtr->insert(facei);
677 678 679 680 681 682 683 684 685 686
            }
        }
    }

    reduce(nConcave, sumOp<label>());
    reduce(maxEdgeSin, maxOp<scalar>());

    if (nConcave > 0)
    {
        scalar maxConcaveDegr =
687
            radToDeg(Foam::asin(Foam::min(1.0, maxEdgeSin)));
688 689 690 691 692 693 694 695 696 697 698 699

        if (debug || report)
        {
            Info<< "   *There are " << nConcave
                << " faces with concave angles between consecutive"
                << " edges. Max concave angle = " << maxConcaveDegr
                << " degrees." << endl;
        }

        return true;
    }

700 701 702
    if (debug || report)
    {
        Info<< "    All angles in faces OK." << endl;
703
    }
704 705

    return false;
706 707 708
}


709
bool Foam::primitiveMesh::checkFaceFlatness
710
(
711 712 713
    const pointField& points,
    const vectorField& faceCentres,
    const vectorField& faceAreas,
714 715 716 717 718
    const bool report,
    const scalar warnFlatness,
    labelHashSet* setPtr
) const
{
Mark Olesen's avatar
Mark Olesen committed
719
    DebugInFunction << "Checking face flatness" << endl;
720 721 722

    if (warnFlatness < 0 || warnFlatness > 1)
    {
723
        FatalErrorInFunction
Mark Olesen's avatar
Mark Olesen committed
724
            << "warnFlatness should be [0..1] but is " << warnFlatness
725 726 727 728 729
            << exit(FatalError);
    }

    const faceList& fcs = faces();

730 731 732 733 734 735 736 737
    tmp<scalarField> tfaceFlatness = primitiveMeshTools::faceFlatness
    (
        *this,
        points,
        faceCentres,
        faceAreas
    );
    const scalarField& faceFlatness = tfaceFlatness();
738

739
    scalarField magAreas(mag(faceAreas));
740 741 742 743

    scalar minFlatness = GREAT;
    scalar sumFlatness = 0;
    label nSummed = 0;
744
    label nWarped = 0;
745

746
    forAll(faceFlatness, facei)
747
    {
748
        if (fcs[facei].size() > 3 && magAreas[facei] > VSMALL)
749
        {
750
            sumFlatness += faceFlatness[facei];
751 752
            nSummed++;

753
            minFlatness = min(minFlatness, faceFlatness[facei]);
754

755
            if (faceFlatness[facei] < warnFlatness)
756 757 758 759 760
            {
                nWarped++;

                if (setPtr)
                {
761
                    setPtr->insert(facei);
762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777
                }
            }
        }
    }


    reduce(nWarped, sumOp<label>());
    reduce(minFlatness, minOp<scalar>());

    reduce(nSummed, sumOp<label>());
    reduce(sumFlatness, sumOp<scalar>());

    if (debug || report)
    {
        if (nSummed > 0)
        {
laurence's avatar
laurence committed
778 779 780
            Info<< "    Face flatness (1 = flat, 0 = butterfly) : min = "
                << minFlatness << "  average = " << sumFlatness / nSummed
                << endl;
781 782 783 784
        }
    }


785
    if (nWarped > 0)
786 787 788 789 790 791 792 793 794 795 796 797 798 799
    {
        if (debug || report)
        {
            Info<< "   *There are " << nWarped
                << " faces with ratio between projected and actual area < "
                << warnFlatness << endl;

            Info<< "    Minimum ratio (minimum flatness, maximum warpage) = "
                << minFlatness << endl;
        }

        return true;
    }

800 801 802
    if (debug || report)
    {
        Info<< "    All face flatness OK." << endl;
803
    }
804 805

    return false;
806 807 808
}


809
bool Foam::primitiveMesh::checkConcaveCells
810
(
811 812
    const vectorField& fAreas,
    const pointField& fCentres,
813 814 815 816
    const bool report,
    labelHashSet* setPtr
) const
{
Mark Olesen's avatar
Mark Olesen committed
817
    DebugInFunction << "Checking for concave cells" << endl;
818

819 820
    const cellList& c = cells();
    const labelList& fOwner = faceOwner();
821

822
    label nConcaveCells = 0;
823

824
    forAll(c, celli)
825
    {
826
        const cell& cFaces = c[celli];
827

828
        bool concave = false;
829

830 831 832 833 834 835
        forAll(cFaces, i)
        {
            if (concave)
            {
                break;
            }
836

837
            label fI = cFaces[i];
838

839
            const point& fC = fCentres[fI];
840

841
            vector fN = fAreas[fI];
842

843
            fN /= max(mag(fN), VSMALL);
844

845 846
            // Flip normal if required so that it is always pointing out of
            // the cell
847
            if (fOwner[fI] != celli)
848 849 850
            {
                fN *= -1;
            }
851

852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873
            // Is the centre of any other face of the cell on the
            // wrong side of the plane of this face?

            forAll(cFaces, j)
            {
                if (j != i)
                {
                    label fJ = cFaces[j];

                    const point& pt = fCentres[fJ];

                    // If the cell is concave, the point will be on the
                    // positive normal side of the plane of f, defined by
                    // its centre and normal, and the angle between (pt -
                    // fC) and fN will be less than 90 degrees, so the dot
                    // product will be positive.

                    vector pC = (pt - fC);

                    pC /= max(mag(pC), VSMALL);

                    if ((pC & fN) > -planarCosAngle_)
874
                    {
875 876 877 878 879
                        // Concave or planar face

                        concave = true;

                        if (setPtr)
880
                        {
881
                            setPtr->insert(celli);
882
                        }
883 884 885 886

                        nConcaveCells++;

                        break;
887 888 889 890 891 892
                    }
                }
            }
        }
    }

893
    reduce(nConcaveCells, sumOp<label>());
894

895
    if (nConcaveCells > 0)
896 897 898
    {
        if (debug || report)
        {
899 900
            Info<< " ***Concave cells (using face planes) found,"
                << " number of cells: " << nConcaveCells << endl;
901 902 903 904
        }

        return true;
    }
905

906 907 908
    if (debug || report)
    {
        Info<< "    Concave cell check OK." << endl;
909
    }
910 911

    return false;
912 913 914
}


915
bool Foam::primitiveMesh::checkUpperTriangular
916 917 918 919 920
(
    const bool report,
    labelHashSet* setPtr
) const
{
Mark Olesen's avatar
Mark Olesen committed
921
    DebugInFunction << "Checking face ordering" << endl;
922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937

    // Check whether internal faces are ordered in the upper triangular order
    const labelList& own = faceOwner();
    const labelList& nei = faceNeighbour();

    const cellList& c = cells();

    label internal = nInternalFaces();

    // Has error occurred?
    bool error = false;
    // Have multiple faces been detected?
    label nMultipleCells = false;

    // Loop through faceCells once more and make sure that for internal cell
    // the first label is smaller
938
    for (label facei = 0; facei < internal; facei++)
939
    {
940
        if (own[facei] >= nei[facei])
941 942 943 944 945
        {
            error  = true;

            if (setPtr)
            {
946
                setPtr->insert(facei);
947 948 949 950 951 952 953 954
            }
        }
    }

    // Loop through all cells. For each cell, find the face that is internal
    // and add it to the check list (upper triangular order).
    // Once the list is completed, check it against the faceCell list

955
    forAll(c, celli)
956
    {
957
        const labelList& curFaces = c[celli];
958 959 960 961 962 963

        // Neighbouring cells
        SortableList<label> nbr(curFaces.size());

        forAll(curFaces, i)
        {
964
            label facei = curFaces[i];
965

966
            if (facei >= nInternalFaces())
967 968 969 970 971 972
            {
                // Sort last
                nbr[i] = labelMax;
            }
            else
            {
973
                label nbrCelli = nei[facei];
974

975
                if (nbrCelli == celli)
976
                {
977
                    nbrCelli = own[facei];
978 979
                }

980
                if (celli < nbrCelli)
981
                {
982
                    // celli is master
983
                    nbr[i] = nbrCelli;
984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063
                }
                else
                {
                    // nbrCell is master. Let it handle this face.
                    nbr[i] = labelMax;
                }
            }
        }

        nbr.sort();

        // Now nbr holds the cellCells in incremental order. Check:
        // - neighbouring cells appear only once. Since nbr is sorted this
        //   is simple check on consecutive elements
        // - faces indexed in same order as nbr are incrementing as well.

        label prevCell = nbr[0];
        label prevFace = curFaces[nbr.indices()[0]];

        bool hasMultipleFaces = false;

        for (label i = 1; i < nbr.size(); i++)
        {
            label thisCell = nbr[i];
            label thisFace = curFaces[nbr.indices()[i]];

            if (thisCell == labelMax)
            {
                break;
            }

            if (thisCell == prevCell)
            {
                hasMultipleFaces = true;

                if (setPtr)
                {
                    setPtr->insert(prevFace);
                    setPtr->insert(thisFace);
                }
            }
            else if (thisFace < prevFace)
            {
                error = true;

                if (setPtr)
                {
                    setPtr->insert(thisFace);
                }
            }

            prevCell = thisCell;
            prevFace = thisFace;
        }

        if (hasMultipleFaces)
        {
            nMultipleCells++;
        }
    }

    reduce(error, orOp<bool>());
    reduce(nMultipleCells, sumOp<label>());

    if ((debug || report) && nMultipleCells > 0)
    {
        Info<< "  <<Found " << nMultipleCells
            << " neighbouring cells with multiple inbetween faces." << endl;
    }

    if (error)
    {
        if (debug || report)
        {
            Info<< " ***Faces not in upper triangular order." << endl;
        }

        return true;
    }

1064 1065 1066
    if (debug || report)
    {
        Info<< "    Upper triangular ordering OK." << endl;
1067
    }
1068 1069

    return false;
1070 1071 1072
}


1073
bool Foam::primitiveMesh::checkCellsZipUp
1074 1075 1076 1077 1078
(
    const bool report,
    labelHashSet* setPtr
) const
{
Mark Olesen's avatar
Mark Olesen committed
1079
    DebugInFunction << "Checking topological cell openness" << endl;
1080 1081 1082 1083 1084 1085

    label nOpenCells = 0;

    const faceList& f = faces();
    const cellList& c = cells();

1086
    forAll(c, celli)
1087
    {
1088
        const labelList& curFaces = c[celli];
1089

1090
        const edgeList cellEdges = c[celli].edges(f);
1091

1092
        labelList edgeUsage(cellEdges.size(), Zero);
1093

1094
        forAll(curFaces, facei)
1095
        {
1096
            edgeList curFaceEdges = f[curFaces[facei]].edges();
1097

1098
            forAll(curFaceEdges, faceEdgeI)
1099 1100 1101
            {
                const edge& curEdge = curFaceEdges[faceEdgeI];

Mark Olesen's avatar