meshTools.C 18.7 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
OpenFOAM bot's avatar
OpenFOAM bot committed
6 7
     \\/     M anipulation  |
-------------------------------------------------------------------------------
OpenFOAM bot's avatar
OpenFOAM bot committed
8 9
    Copyright (C) 2011-2016 OpenFOAM Foundation
    Copyright (C) 2015-2017 OpenCFD Ltd.
10 11 12 13
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

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
25
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
26 27 28 29

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

#include "meshTools.H"
mattijs's avatar
mattijs committed
30
#include "polyMesh.H"
31
#include "hexMatcher.H"
32
#include "treeBoundBox.H"
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

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

bool Foam::meshTools::visNormal
(
    const vector& n,
    const vectorField& faceNormals,
    const labelList& faceLabels
)
{
    forAll(faceLabels, i)
    {
        if ((faceNormals[faceLabels[i]] & n) < SMALL)
        {
            // Found normal in different direction from n.
            return false;
        }
    }
51

52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
    return true;
}


Foam::vectorField Foam::meshTools::calcBoxPointNormals(const primitivePatch& pp)
{
    vectorField octantNormal(8);
    octantNormal[mXmYmZ] = vector(-1, -1, -1);
    octantNormal[pXmYmZ] = vector(1, -1, -1);
    octantNormal[mXpYmZ] = vector(-1, 1, -1);
    octantNormal[pXpYmZ] = vector(1, 1, -1);

    octantNormal[mXmYpZ] = vector(-1, -1, 1);
    octantNormal[pXmYpZ] = vector(1, -1, 1);
    octantNormal[mXpYpZ] = vector(-1, 1, 1);
    octantNormal[pXpYpZ] = vector(1, 1, 1);
68

69 70 71 72 73 74 75 76 77
    octantNormal /= mag(octantNormal);


    vectorField pn(pp.nPoints());

    const vectorField& faceNormals = pp.faceNormals();
    const vectorField& pointNormals = pp.pointNormals();
    const labelListList& pointFaces = pp.pointFaces();

78
    forAll(pointFaces, pointi)
79
    {
80
        const labelList& pFaces = pointFaces[pointi];
81

82
        if (visNormal(pointNormals[pointi], faceNormals, pFaces))
83
        {
84
            pn[pointi] = pointNormals[pointi];
85 86 87
        }
        else
        {
88 89
            WarningInFunction
                << "Average point normal not visible for point:"
90
                << pp.meshPoints()[pointi] << endl;
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 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 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

            label visOctant =
                mXmYmZMask
              | pXmYmZMask
              | mXpYmZMask
              | pXpYmZMask
              | mXmYpZMask
              | pXmYpZMask
              | mXpYpZMask
              | pXpYpZMask;

            forAll(pFaces, i)
            {
                const vector& n = faceNormals[pFaces[i]];

                if (n.x() > SMALL)
                {
                    // All -x octants become invisible
                    visOctant &= ~mXmYmZMask;
                    visOctant &= ~mXmYpZMask;
                    visOctant &= ~mXpYmZMask;
                    visOctant &= ~mXpYpZMask;
                }
                else if (n.x() < -SMALL)
                {
                    // All +x octants become invisible
                    visOctant &= ~pXmYmZMask;
                    visOctant &= ~pXmYpZMask;
                    visOctant &= ~pXpYmZMask;
                    visOctant &= ~pXpYpZMask;
                }

                if (n.y() > SMALL)
                {
                    visOctant &= ~mXmYmZMask;
                    visOctant &= ~mXmYpZMask;
                    visOctant &= ~pXmYmZMask;
                    visOctant &= ~pXmYpZMask;
                }
                else if (n.y() < -SMALL)
                {
                    visOctant &= ~mXpYmZMask;
                    visOctant &= ~mXpYpZMask;
                    visOctant &= ~pXpYmZMask;
                    visOctant &= ~pXpYpZMask;
                }

                if (n.z() > SMALL)
                {
                    visOctant &= ~mXmYmZMask;
                    visOctant &= ~mXpYmZMask;
                    visOctant &= ~pXmYmZMask;
                    visOctant &= ~pXpYmZMask;
                }
                else if (n.z() < -SMALL)
                {
                    visOctant &= ~mXmYpZMask;
                    visOctant &= ~mXpYpZMask;
                    visOctant &= ~pXmYpZMask;
                    visOctant &= ~pXpYpZMask;
                }
            }

            label visI = -1;

            label mask = 1;

            for (label octant = 0; octant < 8; octant++)
            {
                if (visOctant & mask)
                {
                    // Take first visible octant found. Note: could use
                    // combination of visible octants here.
                    visI = octant;

                    break;
                }
                mask <<= 1;
            }

            if (visI != -1)
            {
                // Take a vector in this octant.
174
                pn[pointi] = octantNormal[visI];
175 176 177
            }
            else
            {
178
                pn[pointi] = Zero;
179

180
                WarningInFunction
181 182 183
                    << "No visible octant for point:" << pp.meshPoints()[pointi]
                    << " cooord:" << pp.points()[pp.meshPoints()[pointi]] << nl
                    << "Normal set to " << pn[pointi] << endl;
184 185 186 187 188 189 190 191 192 193 194 195 196 197
            }
        }
    }

    return pn;
}


Foam::vector Foam::meshTools::normEdgeVec
(
    const primitiveMesh& mesh,
    const label edgeI
)
{
198
    return mesh.edges()[edgeI].unitVec(mesh.points());
199 200 201 202 203 204 205 206 207 208 209 210 211
}


void Foam::meshTools::writeOBJ
(
    Ostream& os,
    const point& pt
)
{
    os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
}


212 213 214 215 216 217 218 219 220 221 222 223 224 225 226
void Foam::meshTools::writeOBJ
(
    Ostream& os,
    const UList<point>& pts
)
{
    forAll(pts, i)
    {
        const point& pt = pts[i];
        os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
    }

}


227 228 229 230
void Foam::meshTools::writeOBJ
(
    Ostream& os,
    const triad& t,
231
    const point& origin
232 233 234 235
)
{
    forAll(t, dirI)
    {
236
        writeOBJ(os, origin, origin + t[dirI]);
237 238 239 240
    }
}


241 242 243 244 245 246 247 248
void Foam::meshTools::writeOBJ
(
    Ostream& os,
    const point& p1,
    const point& p2,
    label& count
)
{
249 250 251
    os << "v " << p1.x() << ' ' << p1.y() << ' ' << p1.z() << nl;
    os << "v " << p2.x() << ' ' << p2.y() << ' ' << p2.z() << nl;
    os << "l " << (count + 1) << " " << (count + 2) << endl;
laurence's avatar
laurence committed
252

253 254 255 256 257 258 259 260 261 262 263
    count += 2;
}


void Foam::meshTools::writeOBJ
(
    Ostream& os,
    const point& p1,
    const point& p2
)
{
264 265 266 267 268 269 270 271 272 273 274 275 276 277 278
    os  << "v " << p1.x() << ' ' << p1.y() << ' ' << p1.z() << nl;
    os  << "vn "
        << (p2.x() - p1.x()) << ' '
        << (p2.y() - p1.y()) << ' '
        << (p2.z() - p1.z()) << endl;
}


void Foam::meshTools::writeOBJ
(
    Ostream& os,
    const treeBoundBox& bb
)
{
    writeOBJ(os, bb.points());
279

280
    for (const edge& e : treeBoundBox::edges)
281 282 283
    {
        os << "l " << (e[0] + 1) <<  ' ' << (e[1] + 1) << nl;
    }
284 285 286
}


287 288 289 290 291
void Foam::meshTools::writeOBJ
(
    Ostream& os,
    const cellList& cells,
    const faceList& faces,
292
    const UList<point>& points,
293 294 295 296 297
    const labelList& cellLabels
)
{
    labelHashSet usedFaces(4*cellLabels.size());

298
    for (const label celli : cellLabels)
299
    {
300
        usedFaces.insert(cells[celli]);
301 302 303 304 305 306 307 308 309
    }

    writeOBJ(os, faces, points, usedFaces.toc());
}


bool Foam::meshTools::edgeOnCell
(
    const primitiveMesh& mesh,
310
    const label celli,
311 312 313
    const label edgeI
)
{
314
    return mesh.edgeCells(edgeI).found(celli);
315
}
316 317 318 319 320


bool Foam::meshTools::edgeOnFace
(
    const primitiveMesh& mesh,
321
    const label facei,
322 323 324
    const label edgeI
)
{
325
    return mesh.faceEdges(facei).found(edgeI);
326
}
327 328 329 330 331


bool Foam::meshTools::faceOnCell
(
    const primitiveMesh& mesh,
332 333
    const label celli,
    const label facei
334 335
)
{
336
    if (mesh.isInternalFace(facei))
337 338 339
    {
        if
        (
340 341
            (mesh.faceOwner()[facei] == celli)
         || (mesh.faceNeighbour()[facei] == celli)
342 343 344 345 346 347 348
        )
        {
            return true;
        }
    }
    else
    {
349
        if (mesh.faceOwner()[facei] == celli)
350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367
        {
            return true;
        }
    }
    return false;
}


Foam::label Foam::meshTools::findEdge
(
    const edgeList& edges,
    const labelList& candidates,
    const label v0,
    const label v1
)
{
    forAll(candidates, i)
    {
368
        const label edgeI = candidates[i];
369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430

        const edge& e = edges[edgeI];

        if ((e[0] == v0 && e[1] == v1) || (e[0] == v1 && e[1] == v0))
        {
            return edgeI;
        }
    }
    return -1;
}


Foam::label Foam::meshTools::findEdge
(
    const primitiveMesh& mesh,
    const label v0,
    const label v1
)
{
    const edgeList& edges = mesh.edges();

    const labelList& v0Edges = mesh.pointEdges()[v0];

    forAll(v0Edges, i)
    {
        label edgeI = v0Edges[i];

        const edge& e = edges[edgeI];

        if ((e.start() == v1) || (e.end() == v1))
        {
            return edgeI;
        }
    }
    return -1;
}


Foam::label Foam::meshTools::getSharedEdge
(
    const primitiveMesh& mesh,
    const label f0,
    const label f1
)
{
    const labelList& f0Edges = mesh.faceEdges()[f0];
    const labelList& f1Edges = mesh.faceEdges()[f1];

    forAll(f0Edges, f0EdgeI)
    {
        label edge0 = f0Edges[f0EdgeI];

        forAll(f1Edges, f1EdgeI)
        {
            label edge1 = f1Edges[f1EdgeI];

            if (edge0 == edge1)
            {
                return edge0;
            }
        }
    }
431 432
    FatalErrorInFunction
        << "Faces " << f0 << " and " << f1 << " do not share an edge"
433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448
        << abort(FatalError);

    return -1;

}


Foam::label Foam::meshTools::getSharedFace
(
    const primitiveMesh& mesh,
    const label cell0I,
    const label cell1I
)
{
    const cell& cFaces = mesh.cells()[cell0I];

449
    forAll(cFaces, cFacei)
450
    {
451
        label facei = cFaces[cFacei];
452 453 454

        if
        (
455
            mesh.isInternalFace(facei)
456
         && (
457 458
                mesh.faceOwner()[facei] == cell1I
             || mesh.faceNeighbour()[facei] == cell1I
459 460 461
            )
        )
        {
462
            return facei;
463 464 465 466
        }
    }


467 468
    FatalErrorInFunction
        << "No common face for"
469 470 471 472 473 474 475 476 477 478 479 480
        << "  cell0I:" << cell0I << "  faces:" << cFaces
        << "  cell1I:" << cell1I << "  faces:"
        << mesh.cells()[cell1I]
        << abort(FatalError);

    return -1;
}


void Foam::meshTools::getEdgeFaces
(
    const primitiveMesh& mesh,
481
    const label celli,
482 483 484 485 486
    const label edgeI,
    label& face0,
    label& face1
)
{
mattijs's avatar
mattijs committed
487
    const labelList& eFaces = mesh.edgeFaces(edgeI);
488 489 490 491

    face0 = -1;
    face1 = -1;

492
    forAll(eFaces, eFacei)
493
    {
494
        label facei = eFaces[eFacei];
495

496
        if (faceOnCell(mesh, celli, facei))
497 498 499
        {
            if (face0 == -1)
            {
500
                face0 = facei;
501 502 503
            }
            else
            {
504
                face1 = facei;
505 506 507 508 509 510 511 512

                return;
            }
        }
    }

    if ((face0 == -1) || (face1 == -1))
    {
513 514
        FatalErrorInFunction
            << "Can not find faces using edge " << mesh.edges()[edgeI]
515
            << " on cell " << celli << abort(FatalError);
516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542
    }
}


Foam::label Foam::meshTools::otherEdge
(
    const primitiveMesh& mesh,
    const labelList& edgeLabels,
    const label thisEdgeI,
    const label thisVertI
)
{
    forAll(edgeLabels, edgeLabelI)
    {
        label edgeI = edgeLabels[edgeLabelI];

        if (edgeI != thisEdgeI)
        {
            const edge& e = mesh.edges()[edgeI];

            if ((e.start() == thisVertI) || (e.end() == thisVertI))
            {
                return edgeI;
            }
        }
    }

543 544
    FatalErrorInFunction
        << "Can not find edge in "
545
        << UIndirectList<edge>(mesh.edges(), edgeLabels)
546 547 548 549 550 551 552 553 554 555 556
        << " connected to edge "
        << thisEdgeI << " with vertices " << mesh.edges()[thisEdgeI]
        << " on side " << thisVertI << abort(FatalError);

    return -1;
}


Foam::label Foam::meshTools::otherFace
(
    const primitiveMesh& mesh,
557 558
    const label celli,
    const label facei,
559 560 561 562 563 564
    const label edgeI
)
{
    label face0;
    label face1;

565
    getEdgeFaces(mesh, celli, edgeI, face0, face1);
566

567
    if (face0 == facei)
568 569 570 571 572 573 574 575 576 577 578 579 580
    {
        return face1;
    }
    else
    {
        return face0;
    }
}


Foam::label Foam::meshTools::otherCell
(
    const primitiveMesh& mesh,
581
    const label otherCelli,
582
    const label facei
583 584
)
{
585
    if (!mesh.isInternalFace(facei))
586
    {
587
        FatalErrorInFunction
588
            << "Face " << facei << " is not internal"
589 590 591
            << abort(FatalError);
    }

592
    label newCelli = mesh.faceOwner()[facei];
593

594
    if (newCelli == otherCelli)
595
    {
596
        newCelli = mesh.faceNeighbour()[facei];
597
    }
598
    return newCelli;
599 600 601 602 603 604
}


Foam::label Foam::meshTools::walkFace
(
    const primitiveMesh& mesh,
605
    const label facei,
606 607 608 609 610
    const label startEdgeI,
    const label startVertI,
    const label nEdges
)
{
611
    const labelList& fEdges = mesh.faceEdges(facei);
612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627

    label edgeI = startEdgeI;

    label vertI = startVertI;

    for (label iter = 0; iter < nEdges; iter++)
    {
        edgeI = otherEdge(mesh, fEdges, edgeI, vertI);

        vertI = mesh.edges()[edgeI].otherVertex(vertI);
    }

    return edgeI;
}


628 629 630 631 632
void Foam::meshTools::constrainToMeshCentre
(
    const polyMesh& mesh,
    point& pt
)
mattijs's avatar
mattijs committed
633 634
{
    const Vector<label>& dirs = mesh.geometricD();
635

mattijs's avatar
mattijs committed
636 637 638 639 640 641 642
    const point& min = mesh.bounds().min();
    const point& max = mesh.bounds().max();

    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
    {
        if (dirs[cmpt] == -1)
        {
643
            pt[cmpt] = 0.5*(min[cmpt] + max[cmpt]);
mattijs's avatar
mattijs committed
644 645 646 647 648 649 650 651 652 653 654 655
        }
    }
}


void Foam::meshTools::constrainToMeshCentre
(
    const polyMesh& mesh,
    pointField& pts
)
{
    const Vector<label>& dirs = mesh.geometricD();
656

mattijs's avatar
mattijs committed
657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677
    const point& min = mesh.bounds().min();
    const point& max = mesh.bounds().max();

    bool isConstrained = false;
    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
    {
        if (dirs[cmpt] == -1)
        {
            isConstrained = true;
            break;
        }
    }

    if (isConstrained)
    {
        forAll(pts, i)
        {
            for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
            {
                if (dirs[cmpt] == -1)
                {
678
                    pts[i][cmpt] = 0.5*(min[cmpt] + max[cmpt]);
mattijs's avatar
mattijs committed
679 680 681 682 683 684 685
                }
            }
        }
    }
}


686 687 688 689 690 691
void Foam::meshTools::constrainDirection
(
    const polyMesh& mesh,
    const Vector<label>& dirs,
    vector& d
)
mattijs's avatar
mattijs committed
692 693 694 695 696 697 698 699 700 701 702
{
    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
    {
        if (dirs[cmpt] == -1)
        {
            d[cmpt] = 0.0;
        }
    }
}


703 704 705 706 707 708
void Foam::meshTools::constrainDirection
(
    const polyMesh& mesh,
    const Vector<label>& dirs,
    vectorField& d
)
mattijs's avatar
mattijs committed
709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735
{
    bool isConstrained = false;
    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
    {
        if (dirs[cmpt] == -1)
        {
            isConstrained = true;
            break;
        }
    }

    if (isConstrained)
    {
        forAll(d, i)
        {
            for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
            {
                if (dirs[cmpt] == -1)
                {
                    d[i][cmpt] = 0.0;
                }
            }
        }
    }
}


736 737 738
void Foam::meshTools::getParallelEdges
(
    const primitiveMesh& mesh,
739
    const label celli,
740 741 742 743 744 745 746
    const label e0,
    label& e1,
    label& e2,
    label& e3
)
{
    // Go to any face using e0
747
    label facei = meshTools::otherFace(mesh, celli, -1, e0);
748 749

    // Opposite edge on face
750
    e1 = meshTools::walkFace(mesh, facei, e0, mesh.edges()[e0].end(), 2);
751

752
    facei = meshTools::otherFace(mesh, celli, facei, e1);
753

754
    e2 = meshTools::walkFace(mesh, facei, e1, mesh.edges()[e1].end(), 2);
755

756
    facei = meshTools::otherFace(mesh, celli, facei, e2);
757

758
    e3 = meshTools::walkFace(mesh, facei, e2, mesh.edges()[e2].end(), 2);
759 760 761 762 763 764
}


Foam::vector Foam::meshTools::edgeToCutDir
(
    const primitiveMesh& mesh,
765
    const label celli,
766 767 768
    const label startEdgeI
)
{
769
    if (!hexMatcher().isA(mesh, celli))
770
    {
771
        FatalErrorInFunction
772
            << "Not a hex : cell:" << celli << abort(FatalError);
773 774 775 776 777 778 779
    }


    vector avgVec(normEdgeVec(mesh, startEdgeI));

    label edgeI = startEdgeI;

780
    label facei = -1;
781 782 783 784

    for (label i = 0; i < 3; i++)
    {
        // Step to next face, next edge
785
        facei = meshTools::otherFace(mesh, celli, facei, edgeI);
786 787 788 789 790 791 792 793 794 795 796 797 798 799

        vector eVec(normEdgeVec(mesh, edgeI));

        if ((eVec & avgVec) > 0)
        {
            avgVec += eVec;
        }
        else
        {
            avgVec -= eVec;
        }

        label vertI = mesh.edges()[edgeI].end();

800
        edgeI = meshTools::walkFace(mesh, facei, edgeI, vertI, 2);
801 802
    }

803
    avgVec.normalise();
804 805 806 807 808 809 810 811

    return avgVec;
}


Foam::label Foam::meshTools::cutDirToEdge
(
    const primitiveMesh& mesh,
812
    const label celli,
813 814 815
    const vector& cutDir
)
{
816
    if (!hexMatcher().isA(mesh, celli))
817
    {
818
        FatalErrorInFunction
819
            << "Not a hex : cell:" << celli << abort(FatalError);
820 821
    }

822
    const labelList& cEdges = mesh.cellEdges()[celli];
823 824 825 826 827 828 829 830 831 832 833 834 835 836

    labelHashSet doneEdges(2*cEdges.size());

    scalar maxCos = -GREAT;
    label maxEdgeI = -1;

    for (label i = 0; i < 4; i++)
    {
        forAll(cEdges, cEdgeI)
        {
            label e0 = cEdges[cEdgeI];

            if (!doneEdges.found(e0))
            {
837
                vector avgDir(edgeToCutDir(mesh, celli, e0));
838 839 840 841 842 843 844 845 846 847 848

                scalar cosAngle = mag(avgDir & cutDir);

                if (cosAngle > maxCos)
                {
                    maxCos = cosAngle;
                    maxEdgeI = e0;
                }

                // Mark off edges in cEdges.
                label e1, e2, e3;
849
                getParallelEdges(mesh, celli, e0, e1, e2, e3);
850 851 852 853 854 855

                doneEdges.insert(e0);
                doneEdges.insert(e1);
                doneEdges.insert(e2);
                doneEdges.insert(e3);
            }
856
        }
857 858 859 860 861 862
    }

    forAll(cEdges, cEdgeI)
    {
        if (!doneEdges.found(cEdges[cEdgeI]))
        {
863
            FatalErrorInFunction
864
                << "Cell:" << celli << " edges:" << cEdges << endl
865 866 867 868 869 870 871
                << "Edge:" << cEdges[cEdgeI] << " not yet handled"
                << abort(FatalError);
        }
    }

    if (maxEdgeI == -1)
    {
872 873
        FatalErrorInFunction
            << "Problem : did not find edge aligned with " << cutDir
874
            << " on cell " << celli << abort(FatalError);
875 876 877 878 879 880 881
    }

    return maxEdgeI;
}


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