polyMeshGenChecksGeometry.C 59.5 KB
Newer Older
franjo's avatar
franjo committed
1 2
/*---------------------------------------------------------------------------*\
  =========                 |
Tomislav Lugaric's avatar
Tomislav Lugaric committed
3
  \\      /  F ield         | cfMesh: A library for mesh generation
franjo's avatar
franjo committed
4
   \\    /   O peration     |
5 6
    \\  /    A nd           | Copyright held by the origina author
     \\/     M anipulation  |
franjo's avatar
franjo committed
7 8
-------------------------------------------------------------------------------
License
Tomislav Lugaric's avatar
Tomislav Lugaric committed
9
    This file is part of cfMesh.
franjo's avatar
franjo committed
10

Tomislav Lugaric's avatar
Tomislav Lugaric committed
11
    cfMesh is free software; you can redistribute it and/or modify it
franjo's avatar
franjo committed
12
    under the terms of the GNU General Public License as published by the
Tomislav Lugaric's avatar
Tomislav Lugaric committed
13
    Free Software Foundation; either version 3 of the License, or (at your
franjo's avatar
franjo committed
14 15
    option) any later version.

Tomislav Lugaric's avatar
Tomislav Lugaric committed
16
    cfMesh is distributed in the hope that it will be useful, but WITHOUT
franjo's avatar
franjo committed
17 18 19 20 21
    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
Tomislav Lugaric's avatar
Tomislav Lugaric committed
22
    along with cfMesh.  If not, see <http://www.gnu.org/licenses/>.
franjo's avatar
franjo committed
23 24 25 26 27 28 29 30

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

#include "polyMeshGenChecks.H"
#include "polyMeshGenAddressing.H"
#include "pyramidPointFaceRef.H"
#include "tetrahedron.H"

31
# ifdef USE_OMP
franjo's avatar
franjo committed
32
#include <omp.h>
33
# endif
franjo's avatar
franjo committed
34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace polyMeshGenChecks
{

bool checkClosedBoundary(const polyMeshGen& mesh, const bool report)
{
    // Loop through all boundary faces and sum up the face area vectors.
    // For a closed boundary, this should be zero in all vector components
    vector sumClosed(vector::zero);
    scalar sumMagClosedBoundary = 0;

    const vectorField& areas = mesh.addressingData().faceAreas();

    for(label faceI = mesh.nInternalFaces();faceI<areas.size();++faceI)
    {
        sumClosed += areas[faceI];
        sumMagClosedBoundary += mag(areas[faceI]);
    }

    // Check the openness in the maximum direction (this is APPROXIMATE!)
    scalar maxOpen = max
    (
        sumClosed.component(vector::X),
        max
        (
            sumClosed.component(vector::Y),
            sumClosed.component(vector::Z)
        )
    );

    reduce(sumClosed, sumOp<vector>());
    reduce(maxOpen, maxOp<scalar>());

    if( maxOpen > SMALL*max(1.0, sumMagClosedBoundary) )
    {
        SeriousErrorIn
        (
            "bool checkClosedBoundary(const polyMeshGen&, const bool report)"
        )   << "Possible hole in boundary description" << endl;

        Info<< "Boundary openness in x-direction = "
            << sumClosed.component(vector::X) << endl;

        Info<< "Boundary openness in y-direction = "
            << sumClosed.component(vector::Y) << endl;

        Info<< "Boundary openness in z-direction = "
            << sumClosed.component(vector::Z) << endl;

        return true;
    }
    else
    {
        if( report )
        {
            Info<< "Boundary openness in x-direction = "
                << sumClosed.component(vector::X) << endl;

            Info<< "Boundary openness in y-direction = "
                << sumClosed.component(vector::Y) << endl;

            Info<< "Boundary openness in z-direction = "
                << sumClosed.component(vector::Z) << endl;

            Info<< "Boundary closed (OK)." << endl;
        }

        return false;
    }
}

bool checkClosedCells
(
    const polyMeshGen& mesh,
    const bool report,
    const scalar aspectWarn,
    labelHashSet* setPtr
)
{
    // Check that all cells labels are valid
    const cellListPMG& cells = mesh.cells();
122
    const label nFaces = mesh.faces().size();
franjo's avatar
franjo committed
123 124 125

    label nErrorClosed = 0;

126
    # ifdef USE_OMP
franjo's avatar
franjo committed
127
    # pragma omp parallel for schedule(guided) reduction(+ : nErrorClosed)
128
    # endif
franjo's avatar
franjo committed
129 130 131 132 133 134 135 136 137 138 139 140 141 142 143
    forAll(cells, cI)
    {
        const cell& curCell = cells[cI];

        if( min(curCell) < 0 || max(curCell) > nFaces )
        {
            WarningIn
            (
                "bool checkClosedCells("
                "const polyMeshGen&, const bool, const scalar, labelHashSet*)"
            )   << "Cell " << cI << " contains face labels out of range: "
                << curCell << " Max face index = " << nFaces << endl;

            if( setPtr )
            {
144
                # ifdef USE_OMP
franjo's avatar
franjo committed
145
                # pragma omp critical
146
                # endif
franjo's avatar
franjo committed
147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
                setPtr->insert(cI);
            }

            ++nErrorClosed;
        }
    }

    if( nErrorClosed > 0 )
    {
        SeriousErrorIn
        (
            "bool checkClosedCells("
            "const polyMeshGen&, const bool, const scalar, labelHashSet*)"
        )  << nErrorClosed << " cells with invalid face labels found"
            << endl;

        return true;
    }

    // Loop through cell faces and sum up the face area vectors for each cell.
    // This should be zero in all vector components
    vectorField sumClosed(cells.size(), vector::zero);

    scalarField sumMagClosed(cells.size(), 0.0);

    const labelList& own = mesh.owner();
    const labelList& nei = mesh.neighbour();
    const vectorField& areas = mesh.addressingData().faceAreas();

    forAll(own, faceI)
    {
        // Add to owner
        sumClosed[own[faceI]] += areas[faceI];
        sumMagClosed[own[faceI]] += mag(areas[faceI]);
181

182 183
        if( nei[faceI] == -1 )
            continue;
184

franjo's avatar
franjo committed
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 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339
        // Subtract from neighbour
        sumClosed[nei[faceI]] -= areas[faceI];
        sumMagClosed[nei[faceI]] += mag(areas[faceI]);
    }

    label nOpen = 0;
    scalar maxOpenCell = 0;

    label nAspect = 0;
    scalar maxAspectRatio = 0;

    const scalarField& vols = mesh.addressingData().cellVolumes();

    // Check the sums
    forAll(sumClosed, cellI)
    {
        scalar maxOpen = max
        (
            sumClosed[cellI].component(vector::X),
            max
            (
                sumClosed[cellI].component(vector::Y),
                sumClosed[cellI].component(vector::Z)
            )
        );

        maxOpenCell = max(maxOpenCell, maxOpen);

        if( maxOpen > SMALL*max(1.0, sumMagClosed[cellI]) )
        {
            if( report )
            {
                Pout<< "Cell " << cellI << " is not closed. "
                    << "Face area vectors sum up to " << mag(sumClosed[cellI])
                    << " directionwise " << sumClosed[cellI] << " or "
                    << mag(sumClosed[cellI])
                       /(mag(sumMagClosed[cellI]) + VSMALL)
                    << " of the area of all faces of the cell. " << endl
                    << "    Area magnitudes sum up to "
                    << sumMagClosed[cellI] << endl;
            }

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

            ++nOpen;
        }

        scalar aspectRatio =
            1.0/6.0*sumMagClosed[cellI]/pow(vols[cellI], 2.0/3.0);

        maxAspectRatio = max(maxAspectRatio, aspectRatio);

        if( aspectRatio > aspectWarn )
        {
            if( report )
            {
                Pout<< "High aspect ratio for cell " << cellI
                    << ": " << aspectRatio << endl;
            }

            ++nAspect;
        }
    }

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

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

    if( nOpen > 0 )
    {
        SeriousErrorIn
        (
            "bool checkClosedCells("
            "const polyMeshGen&, const bool, const scalar, labelHashSet*)"
        )   << nOpen << " open cells found. Max cell openness: "
            << maxOpenCell << endl;

        return true;
    }

    if( nAspect > 0 )
    {
        SeriousErrorIn
        (
            "bool checkClosedCells("
            "const polyMeshGen&, const bool, const scalar, labelHashSet*)"
        )   << nAspect << " high aspect ratio cells found.  "
            << "Max aspect ratio: " << maxAspectRatio
            << endl;

        return true;
    }

    if( report )
    {
        Info<< "Max cell openness = " << maxOpenCell
            << "  Max aspect ratio = " << maxAspectRatio
            << ".  All cells OK.\n" << endl;
    }

    return false;
}

bool checkCellVolumes
(
    const polyMeshGen& mesh,
    const bool report,
    labelHashSet* setPtr
)
{
    const scalarField& vols = mesh.addressingData().cellVolumes();

    scalar minVolume = GREAT;
    scalar maxVolume = -GREAT;

    label nNegVolCells = 0;

    forAll(vols, cellI)
    {
        if( vols[cellI] < VSMALL )
        {
            if( report )
                SeriousErrorIn
                (
                    "bool checkCellVolumes("
                    "const polyMeshGen&, const bool, labelHashSet*)"
                )   << "Zero or negative cell volume detected for cell "
                    << cellI << ".  Volume = " << vols[cellI] << endl;

            if( setPtr )
                setPtr->insert(cellI);

            ++nNegVolCells;
        }

        minVolume = min(minVolume, vols[cellI]);
        maxVolume = max(maxVolume, vols[cellI]);
    }

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

    if( minVolume < VSMALL )
    {
        SeriousErrorIn
        (
            "bool checkCellVolumes("
            "const polyMeshGen&, const bool, labelHashSet*)"
        )   << "Zero or negative cell volume detected.  "
340
            << "Minimum negative volume: "
franjo's avatar
franjo committed
341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369
            << minVolume << ".\nNumber of negative volume cells: "
            << nNegVolCells << ".  This mesh is invalid"
            << endl;

        return true;
    }
    else
    {
        if( report )
        {
            Info<< "Min volume = " << minVolume
                << ". Max volume = " << maxVolume
                << ".  Total volume = " << sum(vols)
                << ".  Cell volumes OK.\n" << endl;
        }

        return false;
    }
}

bool checkFaceAreas
(
    const polyMeshGen& mesh,
    const bool report,
    const scalar minFaceArea,
    labelHashSet* setPtr,
    const boolList* changedFacePtr
)
{
370
    const vectorField& faceAreas = mesh.addressingData().faceAreas();
franjo's avatar
franjo committed
371 372 373 374 375 376 377

    const labelList& own = mesh.owner();
    const labelList& nei = mesh.neighbour();

    scalar minArea = VGREAT;
    scalar maxArea = -VGREAT;

378
    # ifdef USE_OMP
franjo's avatar
franjo committed
379
    # pragma omp parallel if( own.size() > 100 )
380
    # endif
franjo's avatar
franjo committed
381 382
    {
        scalar localMaxArea(-VGREAT), localMinArea(VGREAT);
383

384
        # ifdef USE_OMP
franjo's avatar
franjo committed
385
        # pragma omp for schedule(guided)
386
        # endif
387
        forAll(faceAreas, faceI)
franjo's avatar
franjo committed
388 389 390
        {
            if( changedFacePtr && !changedFacePtr->operator[](faceI) )
                continue;
391

392 393 394
            const scalar magFaceArea = mag(faceAreas[faceI]);

            if( magFaceArea < minFaceArea )
franjo's avatar
franjo committed
395 396 397 398 399 400 401 402 403
            {
                if( report )
                {
                    if( nei[faceI] != -1 )
                    {
                        Pout<< "Zero or negative face area detected for "
                            << "internal face " << faceI << " between cells "
                            << own[faceI] << " and " << nei[faceI]
                            << ".  Face area magnitude = "
404
                            << magFaceArea << endl;
franjo's avatar
franjo committed
405 406 407 408 409 410
                    }
                    else
                    {
                        Pout<< "Zero or negative face area detected for "
                            << "boundary face " << faceI << " next to cell "
                            << own[faceI] << ".  Face area magnitude = "
411
                            << magFaceArea << endl;
franjo's avatar
franjo committed
412 413
                    }
                }
414

franjo's avatar
franjo committed
415 416
                if( setPtr )
                {
417
                    # ifdef USE_OMP
franjo's avatar
franjo committed
418
                    # pragma omp critical
419
                    # endif
franjo's avatar
franjo committed
420 421 422
                    setPtr->insert(faceI);
                }
            }
423

424 425
            localMinArea = Foam::min(localMinArea, magFaceArea);
            localMaxArea = Foam::max(localMaxArea, magFaceArea);
franjo's avatar
franjo committed
426
        }
427

428
        # ifdef USE_OMP
franjo's avatar
franjo committed
429
        # pragma omp critical
430
        # endif
franjo's avatar
franjo committed
431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446
        {
            minArea = Foam::min(minArea, localMinArea);
            maxArea = Foam::max(maxArea, localMaxArea);
        }
    }

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

    if( minArea < VSMALL )
    {
        SeriousErrorIn
        (
            "bool checkFaceAreas("
            "const polyMeshGen&, const bool, const scalar,"
            " labelHashSet*, const boolList*)"
447
        )   << "Zero or negative face area detected.  Minimum negative area: "
franjo's avatar
franjo committed
448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478
            << minArea << ". This mesh is invalid"
            << endl;

        return true;
    }
    else
    {
        if( report )
        {
            Info<< "Minumum face area = " << minArea
                << ". Maximum face area = " << maxArea
                << ".  Face area magnitudes OK.\n" << endl;
        }

        return false;
    }
}

bool checkCellPartTetrahedra
(
    const polyMeshGen& mesh,
    const bool report,
    const scalar minPartTet,
    labelHashSet* setPtr,
    const boolList* changedFacePtr
)
{
    const pointFieldPMG& points = mesh.points();
    const faceListPMG& faces = mesh.faces();
    const labelList& owner = mesh.owner();
    const labelList& neighbour = mesh.neighbour();
479

franjo's avatar
franjo committed
480 481 482 483 484
    const vectorField& fCentres = mesh.addressingData().faceCentres();
    const vectorField& cCentres = mesh.addressingData().cellCentres();

    label nNegVolCells = 0;

485
    # ifdef USE_OMP
franjo's avatar
franjo committed
486 487
    # pragma omp parallel for if( owner.size() > 100 ) \
    schedule(guided) reduction(+ : nNegVolCells)
488
    # endif
franjo's avatar
franjo committed
489 490 491 492
    forAll(owner, faceI)
    {
        if( changedFacePtr && !(*changedFacePtr)[faceI] )
            continue;
493

franjo's avatar
franjo committed
494
        const face& f = faces[faceI];
495

franjo's avatar
franjo committed
496
        bool badFace(false);
497

franjo's avatar
franjo committed
498 499 500 501 502 503 504 505 506
        forAll(f, eI)
        {
            const tetrahedron<point, point> tetOwn
            (
                fCentres[faceI],
                points[f.nextLabel(eI)],
                points[f[eI]],
                cCentres[owner[faceI]]
            );
507

franjo's avatar
franjo committed
508 509 510 511
            if( tetOwn.mag() < minPartTet )
            {
                if( report )
                {
512
                    # ifdef USE_OMP
franjo's avatar
franjo committed
513
                    # pragma omp critical
514
                    # endif
franjo's avatar
franjo committed
515 516 517
                    Pout<< "Zero or negative cell volume detected for cell "
                        << owner[faceI] << "." << endl;
                }
518

franjo's avatar
franjo committed
519 520
                badFace = true;
            }
521

franjo's avatar
franjo committed
522 523
            if( neighbour[faceI] < 0 )
                continue;
524

franjo's avatar
franjo committed
525 526 527 528 529 530 531
            const tetrahedron<point, point> tetNei
            (
                fCentres[faceI],
                points[f[eI]],
                points[f.nextLabel(eI)],
                cCentres[neighbour[faceI]]
            );
532

franjo's avatar
franjo committed
533 534 535 536
            if( tetNei.mag() < minPartTet )
            {
                if( report )
                {
537
                    # ifdef USE_OMP
franjo's avatar
franjo committed
538
                    # pragma omp critical
539
                    # endif
franjo's avatar
franjo committed
540 541 542
                    Pout<< "Zero or negative cell volume detected for cell "
                        << neighbour[faceI] << "." << endl;
                }
543

franjo's avatar
franjo committed
544 545 546
                badFace = true;
            }
        }
547

franjo's avatar
franjo committed
548 549 550 551
        if( badFace )
        {
            if( setPtr )
            {
552
                # ifdef USE_OMP
franjo's avatar
franjo committed
553
                # pragma omp critical
554
                # endif
franjo's avatar
franjo committed
555 556
                setPtr->insert(faceI);
            }
557

franjo's avatar
franjo committed
558 559 560
            ++nNegVolCells;
        }
    }
561

franjo's avatar
franjo committed
562 563 564
    if( setPtr )
    {
        //- ensure that faces are selected at both sides
565
        const PtrList<processorBoundaryPatch>& procBnd = mesh.procBoundaries();
franjo's avatar
franjo committed
566 567 568 569
        forAll(procBnd, patchI)
        {
            const label start = procBnd[patchI].patchStart();
            const label size = procBnd[patchI].patchSize();
570

571
            labelLongList sendData;
franjo's avatar
franjo committed
572 573 574 575 576
            for(label faceI=0;faceI<size;++faceI)
            {
                if( setPtr->found(faceI+start) )
                    sendData.append(faceI);
            }
577

franjo's avatar
franjo committed
578 579 580 581 582 583
            OPstream toOtherProc
            (
                Pstream::blocking,
                procBnd[patchI].neiProcNo(),
                sendData.byteSize()
            );
584

franjo's avatar
franjo committed
585 586
            toOtherProc << sendData;
        }
587

franjo's avatar
franjo committed
588 589 590
        forAll(procBnd, patchI)
        {
            labelList receivedData;
591

franjo's avatar
franjo committed
592 593 594 595 596
            IPstream fromOtherProc
            (
                Pstream::blocking,
                procBnd[patchI].neiProcNo()
            );
597

franjo's avatar
franjo committed
598
            fromOtherProc >> receivedData;
599

franjo's avatar
franjo committed
600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628
            const label start = procBnd[patchI].patchStart();
            forAll(receivedData, i)
                setPtr->insert(start+receivedData[i]);
        }
    }

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

    if( nNegVolCells != 0 )
    {
        WarningIn
        (
            "bool checkCellPartTetrahedra("
            "const polyMeshGen&, const bool, const scalar,"
            " labelHashSet*, const boolList*)"
        )   << nNegVolCells << " zero or negative part tetrahedra detected."
            << endl;

        return true;
    }
    else
    {
        if( report )
            Info<< "Part cell tetrahedra OK.\n" << endl;

        return false;
    }
}

629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727
void checkFaceDotProduct
(
    const polyMeshGen& mesh,
    scalarField& faceDotProduct,
    const boolList* changedFacePtr
)
{
    //- for all internal faces check theat the d dot S product is positive
    const vectorField& centres = mesh.addressingData().cellCentres();
    const vectorField& areas = mesh.addressingData().faceAreas();

    const labelList& own = mesh.owner();
    const labelList& nei = mesh.neighbour();
    const label nInternalFaces = mesh.nInternalFaces();

    faceDotProduct.setSize(own.size());
    faceDotProduct = 1.0;

    # ifdef USE_OMP
    # pragma omp parallel if( nInternalFaces > 1000 )
    # endif
    {
        # ifdef USE_OMP
        # pragma omp for schedule(dynamic, 100)
        # endif
        for(label faceI=0;faceI<nInternalFaces;++faceI)
        {
            if( changedFacePtr && !(*changedFacePtr)[faceI] )
                continue;

            const vector d = centres[nei[faceI]] - centres[own[faceI]];
            const vector& s = areas[faceI];

            faceDotProduct[faceI] = (d & s)/(mag(d)*mag(s) + VSMALL);
        }
    }

    if( Pstream::parRun() )
    {
        const PtrList<processorBoundaryPatch>& procBoundaries =
            mesh.procBoundaries();

        forAll(procBoundaries, patchI)
        {
            const label start = procBoundaries[patchI].patchStart();

            vectorField cCentres(procBoundaries[patchI].patchSize());
            forAll(cCentres, faceI)
                cCentres[faceI] = centres[own[start+faceI]];

            OPstream toOtherProc
            (
                Pstream::blocking,
                procBoundaries[patchI].neiProcNo(),
                cCentres.byteSize()
            );

            toOtherProc << cCentres;
        }

        forAll(procBoundaries, patchI)
        {
            vectorField otherCentres;
            IPstream fromOtherProc
            (
                Pstream::blocking,
                procBoundaries[patchI].neiProcNo()
            );

            fromOtherProc >> otherCentres;

            //- calculate skewness at processor faces
            const label start = procBoundaries[patchI].patchStart();
            # ifdef USE_OMP
            # pragma omp parallel
            # endif
            {
                # ifdef USE_OMP
                # pragma omp for schedule(dynamic, 100)
                # endif
                forAll(otherCentres, fI)
                {
                    const label faceI = start + fI;

                    if( changedFacePtr && !(*changedFacePtr)[faceI] )
                        continue;

                    const point& cOwn = centres[own[faceI]];
                    const point& cNei = otherCentres[fI];
                    const vector d = cNei - cOwn;
                    const vector& s = areas[faceI];

                    faceDotProduct[faceI] = (d & s)/(mag(d)*mag(s) + VSMALL);
                }
            }
        }
    }
}

franjo's avatar
franjo committed
728 729 730 731 732 733 734 735 736
bool checkFaceDotProduct
(
    const polyMeshGen& mesh,
    const bool report,
    const scalar nonOrthWarn,
    labelHashSet* setPtr,
    const boolList* changedFacePtr
)
{
737 738 739
    //- for all internal faces check theat the d dot S product is positive
    scalarField faceDotProduct;
    checkFaceDotProduct(mesh, faceDotProduct, changedFacePtr);
franjo's avatar
franjo committed
740 741 742

    const labelList& own = mesh.owner();
    const labelList& nei = mesh.neighbour();
743
    const label nInternalFaces = mesh.nInternalFaces();
franjo's avatar
franjo committed
744 745 746

    // Severe nonorthogonality threshold
    const scalar severeNonorthogonalityThreshold =
747
        ::cos(nonOrthWarn/180.0*M_PI);
franjo's avatar
franjo committed
748 749 750 751 752 753 754 755

    scalar minDDotS = VGREAT;
    scalar sumDDotS = 0.0;

    label severeNonOrth = 0;
    label errorNonOrth = 0;
    label counter = 0;

756
    # ifdef USE_OMP
franjo's avatar
franjo committed
757 758
    # pragma omp parallel if( nInternalFaces > 1000 ) \
    reduction(+ : severeNonOrth, errorNonOrth, sumDDotS)
759
    # endif
franjo's avatar
franjo committed
760 761
    {
        scalar localMinDDotS(VGREAT);
762
        # ifdef USE_OMP
franjo's avatar
franjo committed
763
        # pragma omp for schedule(guided)
764
        # endif
franjo's avatar
franjo committed
765 766
        for(label faceI=0;faceI<nInternalFaces;++faceI)
        {
767 768 769
            if( changedFacePtr && !(*changedFacePtr)[faceI] )
                continue;

770
            const scalar dDotS = faceDotProduct[faceI];
771

franjo's avatar
franjo committed
772 773 774 775 776 777 778
            if( dDotS < severeNonorthogonalityThreshold )
            {
                if( dDotS > SMALL )
                {
                    if( report )
                    {
                        // Severe non-orthogonality but mesh still OK
779
                        # ifdef USE_OMP
franjo's avatar
franjo committed
780
                        # pragma omp critical
781
                        # endif
franjo's avatar
franjo committed
782 783 784 785
                        Pout<< "Severe non-orthogonality for face " << faceI
                            << " between cells " << own[faceI]
                            << " and " << nei[faceI]
                            << ": Angle = "
786
                            << ::acos(dDotS)/M_PI*180.0
franjo's avatar
franjo committed
787 788
                            << " deg." << endl;
                    }
789

franjo's avatar
franjo committed
790 791
                    if( setPtr )
                    {
792
                        # ifdef USE_OMP
franjo's avatar
franjo committed
793
                        # pragma omp critical
794
                        # endif
franjo's avatar
franjo committed
795 796
                        setPtr->insert(faceI);
                    }
797

franjo's avatar
franjo committed
798 799 800 801 802
                    ++severeNonOrth;
                }
                else
                {
                    ++errorNonOrth;
803

franjo's avatar
franjo committed
804 805
                    if( setPtr )
                    {
806
                        # ifdef USE_OMP
franjo's avatar
franjo committed
807
                        # pragma omp critical
808
                        # endif
franjo's avatar
franjo committed
809 810 811 812
                        setPtr->insert(faceI);
                    }
                }
            }
813

franjo's avatar
franjo committed
814 815 816
            localMinDDotS = Foam::min(dDotS, localMinDDotS);
            sumDDotS += dDotS;
        }
817

818
        # ifdef USE_OMP
franjo's avatar
franjo committed
819
        # pragma omp critical
820
        # endif
franjo's avatar
franjo committed
821 822
        minDDotS = Foam::min(minDDotS, localMinDDotS);
    }
823

franjo's avatar
franjo committed
824
    counter += nInternalFaces;
825

franjo's avatar
franjo committed
826 827
    if( Pstream::parRun() )
    {
828
        const PtrList<processorBoundaryPatch>& procBoundaries =
franjo's avatar
franjo committed
829
            mesh.procBoundaries();
830

831
        const label start = procBoundaries[0].patchStart();
832

833 834 835 836
        # ifdef USE_OMP
        # pragma omp parallel \
        reduction(+ : severeNonOrth, errorNonOrth, sumDDotS, counter)
        # endif
franjo's avatar
franjo committed
837
        {
838
            scalar localMinDDotS(VGREAT);
839

840
            # ifdef USE_OMP
841
            # pragma omp for schedule(dynamic, 100)
842
            # endif
843
            for(label faceI=start;faceI<own.size();++faceI)
franjo's avatar
franjo committed
844
            {
845 846
                if( changedFacePtr && !(*changedFacePtr)[start+faceI] )
                            continue;
847

848
                const scalar dDotS = faceDotProduct[faceI];
849

850 851 852
                if( dDotS < severeNonorthogonalityThreshold )
                {
                    if( dDotS > SMALL )
franjo's avatar
franjo committed
853
                    {
854
                        if( report )
franjo's avatar
franjo committed
855
                        {
856 857 858 859
                            // Severe non-orthogonality but mesh still OK
                            # ifdef USE_OMP
                            # pragma omp critical
                            # endif
franjo's avatar
franjo committed
860
                            {
861 862 863 864 865 866 867 868 869 870
                                const scalar angle
                                (
                                    Foam::acos(dDotS) /
                                    M_PI * 180.0
                                );
                                Pout<< "Severe non-orthogonality for face "
                                    << start+faceI
                                    << ": Angle = "
                                    << angle
                                    << " deg." << endl;
franjo's avatar
franjo committed
871 872
                            }
                        }
873

874 875 876 877 878 879
                        if( setPtr )
                        {
                            # ifdef USE_OMP
                            # pragma omp critical
                            # endif
                            setPtr->insert(start+faceI);
franjo's avatar
franjo committed
880
                        }
881 882

                        ++severeNonOrth;
franjo's avatar
franjo committed
883
                    }
884 885 886
                    else
                    {
                        ++errorNonOrth;
887

888 889 890 891 892 893 894 895
                        if( setPtr )
                        {
                            # ifdef USE_OMP
                            # pragma omp critical
                            # endif
                            setPtr->insert(start+faceI);
                        }
                    }
franjo's avatar
franjo committed
896
                }
897

898 899 900 901
                localMinDDotS = Foam::min(dDotS, localMinDDotS);
                sumDDotS += 0.5 * dDotS;

                ++counter;
franjo's avatar
franjo committed
902
            }
903

904 905 906 907
            # ifdef USE_OMP
            # pragma omp critical
            # endif
            minDDotS = Foam::min(minDDotS, localMinDDotS);
franjo's avatar
franjo committed
908 909
        }
    }
910

franjo's avatar
franjo committed
911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933
    reduce(minDDotS, minOp<scalar>());
    reduce(sumDDotS, sumOp<scalar>());
    reduce(severeNonOrth, sumOp<label>());
    reduce(errorNonOrth, sumOp<label>());

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

    // Only report if there are some internal faces
    if( counter > 0 )
    {
        if( minDDotS < severeNonorthogonalityThreshold )
        {
            Info<< "Number of non-orthogonality errors: " << errorNonOrth
                << ". Number of severely non-orthogonal faces: "
                << severeNonOrth  << "." << endl;
        }
    }

    if( report )
    {
        if( counter > 0 )
        {
            Info<< "Mesh non-orthogonality Max: "
934
                << ::acos(minDDotS)/M_PI*180.0
franjo's avatar
franjo committed
935
                << " average: " <<
936
                   ::acos(sumDDotS/counter)/M_PI*180.0
franjo's avatar
franjo committed
937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980
                << endl;
        }
    }

    if( errorNonOrth > 0 )
    {
        WarningIn
        (
            "checkFaceDotProduct("
            "const polyMeshGen&, const bool, const scalar,"
            " labelHashSet*, const boolList*)"
        )   << "Error in non-orthogonality detected" << endl;

        return true;
    }
    else
    {
        if( report )
            Info<< "Non-orthogonality check OK.\n" << endl;

        return false;
    }
}

bool checkFacePyramids
(
    const polyMeshGen& mesh,
    const bool report,
    const scalar minPyrVol,
    labelHashSet* setPtr,
    const boolList* changedFacePtr
)
{
    // check whether face area vector points to the cell with higher label
    const vectorField& ctrs = mesh.addressingData().cellCentres();

    const labelList& owner = mesh.owner();
    const labelList& neighbour = mesh.neighbour();

    const faceListPMG& faces = mesh.faces();
    const pointFieldPMG& points = mesh.points();

    label nErrorPyrs = 0;

981
    # ifdef USE_OMP
franjo's avatar
franjo committed
982
    # pragma omp parallel for schedule(guided) reduction(+ : nErrorPyrs)
983
    # endif
franjo's avatar
franjo committed
984 985
    forAll(faces, faceI)
    {
986 987
        if( changedFacePtr && !(*changedFacePtr)[faceI] )
            continue;
988

franjo's avatar
franjo committed
989 990
        // Create the owner pyramid - it will have negative volume
        const scalar pyrVol = pyramidPointFaceRef
991 992 993 994
        (
            faces[faceI],
            ctrs[owner[faceI]]
        ).mag(points);
995

franjo's avatar
franjo committed
996 997 998 999 1000 1001
        bool badFace(false);

        if( pyrVol > -minPyrVol )
        {
            if( report )
            {
1002
                # ifdef USE_OMP
franjo's avatar
franjo committed
1003
                # pragma omp critical
1004
                # endif
franjo's avatar
franjo committed
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
                Pout<< "bool checkFacePyramids("
                    << "const bool, const scalar, labelHashSet*) : "
                    << "face " << faceI << " points the wrong way. " << endl
                    << "Pyramid volume: " << -pyrVol
                    << " Face " << faces[faceI] << " area: "
                    << faces[faceI].mag(points)
                    << " Owner cell: " << owner[faceI] << endl
                    << "Owner cell vertex labels: "
                    << mesh.cells()[owner[faceI]].labels(faces)
                    << endl;
            }

            badFace = true;
        }

        if( neighbour[faceI] != -1 )
        {
            // Create the neighbour pyramid - it will have positive volume
            const scalar pyrVol =
                pyramidPointFaceRef
                (
                    faces[faceI],
                    ctrs[neighbour[faceI]]
                ).mag(points);

            if( pyrVol < minPyrVol )
            {
                if( report )
                {
1034
                    # ifdef USE_OMP
franjo's avatar
franjo committed
1035
                    # pragma omp critical
1036
                    # endif
franjo's avatar
franjo committed
1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051
                    Pout<< "bool checkFacePyramids("
                        << "const bool, const scalar, labelHashSet*) : "
                        << "face " << faceI << " points the wrong way. " << endl
                        << "Pyramid volume: " << -pyrVol
                        << " Face " << faces[faceI] << " area: "
                        << faces[faceI].mag(points)
                        << " Neighbour cell: " << neighbour[faceI] << endl
                        << "Neighbour cell vertex labels: "
                        << mesh.cells()[neighbour[faceI]].labels(faces)
                        << endl;
                }

                badFace = true;
            }
        }
1052

franjo's avatar
franjo committed
1053 1054 1055 1056
        if( badFace )
        {
            if( setPtr )
            {
1057
                # ifdef USE_OMP
franjo's avatar
franjo committed
1058
                # pragma omp critical
1059
                # endif
franjo's avatar
franjo committed
1060 1061
                setPtr->insert(faceI);
            }
1062

franjo's avatar
franjo committed
1063 1064 1065 1066 1067
            ++nErrorPyrs;
        }
    }

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

franjo's avatar
franjo committed
1069 1070 1071
    if( setPtr )
    {
        //- make sure that processor faces are marked on both sides
1072
        const PtrList<processorBoundaryPatch>& procBoundaries =
franjo's avatar
franjo committed
1073
            mesh.procBoundaries();
1074

franjo's avatar
franjo committed
1075 1076 1077 1078 1079
        //- send and receive data where needed
        forAll(procBoundaries, patchI)
        {
            const label start = procBoundaries[patchI].patchStart();
            const label size = procBoundaries[patchI].patchSize();
1080

1081
            labelLongList markedFaces;
franjo's avatar
franjo committed
1082 1083 1084 1085 1086
            for(label faceI=0;faceI<size;++faceI)
            {
                if( setPtr->found(start+faceI) )
                    markedFaces.append(faceI);
            }
1087

franjo's avatar
franjo committed
1088 1089 1090 1091 1092 1093
            OPstream toOtherProc
            (
                Pstream::blocking,
                procBoundaries[patchI].neiProcNo(),
                markedFaces.byteSize()
            );
1094

franjo's avatar
franjo committed
1095 1096
            toOtherProc << markedFaces;
        }
1097

franjo's avatar
franjo committed
1098 1099 1100 1101 1102 1103 1104 1105 1106
        forAll(procBoundaries, patchI)
        {
            labelList receivedData;
            IPstream fromOtheProc
            (
                Pstream::blocking,
                procBoundaries[patchI].neiProcNo()
            );
            fromOtheProc >> receivedData;
1107

franjo's avatar
franjo committed
1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136
            const label start = procBoundaries[patchI].patchStart();
            forAll(receivedData, i)
                setPtr->insert(start+receivedData[i]);
        }
    }

    if( nErrorPyrs > 0 )
    {
        if( Pstream::master() )
            WarningIn
            (
                "bool checkFacePyramids("
                "const polyMeshGen&, const bool, const scalar,"
                " labelHashSet*, const boolList*)"
            )   << "Error in face pyramids: " << nErrorPyrs
                << " faces pointing the wrong way!"
                << endl;

        return true;
    }
    else
    {
        if( report )
            Info<< "Face pyramids OK.\n" << endl;

        return false;
    }
}

1137
void checkFaceSkewness
franjo's avatar
franjo committed
1138 1139
(
    const polyMeshGen& mesh,
1140
    scalarField& faceSkewness,
franjo's avatar
franjo committed
1141 1142 1143 1144 1145 1146 1147
    const boolList* changedFacePtr
)
{
    //- Warn if the skew correction vector is more than skewWarning times
    //- larger than the face area vector
    const labelList& own = mesh.owner();
    const labelList& nei = mesh.neighbour();
1148
    const label nInternalFaces = mesh.nInternalFaces();
1149

franjo's avatar
franjo committed
1150 1151 1152
    const vectorField& centres = mesh.addressingData().cellCentres();
    const vectorField& fCentres = mesh.addressingData().faceCentres();

1153 1154
    faceSkewness.setSize(own.size());
    faceSkewness = 0.0;
franjo's avatar
franjo committed
1155 1156

    //- check internal faces
1157
    # ifdef USE_OMP
1158
    # pragma omp parallel for schedule(dynamic, 100)
1159
    # endif
1160
    for(label faceI=0;faceI<nInternalFaces;++faceI)
franjo's avatar
franjo committed
1161
    {
1162 1163
        if( changedFacePtr && !changedFacePtr->operator[](faceI) )
            continue;
1164

1165 1166
        const scalar dOwn = mag(fCentres[faceI] - centres[own[faceI]]);
        const scalar dNei = mag(fCentres[faceI] - centres[nei[faceI]]);
1167

1168 1169 1170
        const point faceIntersection =
            centres[own[faceI]]*dNei/(dOwn+dNei)
          + centres[nei[faceI]]*dOwn/(dOwn+dNei);
1171

1172 1173 1174
        faceSkewness[faceI] =
            mag(fCentres[faceI] - faceIntersection)
            /(mag(centres[nei[faceI]] - centres[own[faceI]]) + VSMALL);
franjo's avatar
franjo committed
1175
    }
1176

franjo's avatar
franjo committed
1177 1178 1179
    if( Pstream::parRun() )
    {
        //- check parallel boundaries
1180
        const PtrList<processorBoundaryPatch>& procBoundaries =
franjo's avatar
franjo committed
1181
            mesh.procBoundaries();
1182

franjo's avatar
franjo committed
1183 1184 1185
        forAll(procBoundaries, patchI)
        {
            const label start = procBoundaries[patchI].patchStart();
1186

franjo's avatar
franjo committed
1187 1188 1189
            vectorField cCentres(procBoundaries[patchI].patchSize());
            forAll(cCentres, faceI)
                cCentres[faceI] = centres[own[start+faceI]];
1190

franjo's avatar
franjo committed
1191 1192 1193 1194 1195 1196
            OPstream toOtherProc
            (
                Pstream::blocking,
                procBoundaries[patchI].neiProcNo(),
                cCentres.byteSize()
            );
1197

franjo's avatar
franjo committed
1198 1199
            toOtherProc << cCentres;
        }
1200

franjo's avatar
franjo committed
1201 1202 1203
        forAll(procBoundaries, patchI)
        {
            const label start = procBoundaries[patchI].patchStart();
1204

franjo's avatar
franjo committed
1205 1206 1207 1208 1209 1210
            vectorField otherCentres;
            IPstream fromOtheProc
            (
                Pstream::blocking,
                procBoundaries[patchI].neiProcNo()
            );
1211

franjo's avatar
franjo committed
1212
            fromOtheProc >> otherCentres;
1213

franjo's avatar
franjo committed
1214
            //- calculate skewness at processor faces
1215
            # ifdef USE_OMP
1216
            # pragma omp parallel
1217
            # endif
franjo's avatar
franjo committed
1218
            {
1219
                # ifdef USE_OMP
1220
                # pragma omp for schedule(dynamic, 100)
1221
                # endif
1222
                forAll(otherCentres, fI)
franjo's avatar
franjo committed
1223
                {
1224 1225
                    const label faceI = start + fI;