tetOverlapVolumeTemplates.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.
 
    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.
 
    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
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
\*---------------------------------------------------------------------------*/
 
#include "tetOverlapVolume.H"
#include "primitiveMesh.H"
 
// * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * * //
 
template<class tetPointsOp>
void Foam::tetOverlapVolume::tetTetOverlap
(
    const tetPoints& tetA,
    const tetPoints& tetB,
    tetPointsOp& insideOp
)
{
    static tetPointRef::tetIntersectionList insideTets;
    label nInside = 0;
    static tetPointRef::tetIntersectionList cutInsideTets;
    label nCutInside = 0;
 
    tetPointRef::storeOp inside(insideTets, nInside);
    tetPointRef::storeOp cutInside(cutInsideTets, nCutInside);
    tetPointRef::dummyOp outside;
 
    if (tetA.tet().mag() < minTetVolume_ || tetB.tet().mag() < minTetVolume_)
    {
        return;
    }
 
    // face0
    plane pl0(tetB[1], tetB[3], tetB[2]);
    tetA.tet().sliceWithPlane(pl0, cutInside, outside);
    if (nCutInside == 0)
    {
        return;
    }
 
    // face1
    plane pl1(tetB[0], tetB[2], tetB[3]);
    nInside = 0;
    for (label i = 0; i < nCutInside; i++)
    {
        const tetPointRef t = cutInsideTets[i].tet();
        t.sliceWithPlane(pl1, inside, outside);
    }
    if (nInside == 0)
    {
        return;
    }
 
    // face2
    plane pl2(tetB[0], tetB[3], tetB[1]);
    nCutInside = 0;
    for (label i = 0; i < nInside; i++)
    {
        const tetPointRef t = insideTets[i].tet();
        t.sliceWithPlane(pl2, cutInside, outside);
    }
    if (nCutInside == 0)
    {
        return;
    }
 
    // face3
    plane pl3(tetB[0], tetB[1], tetB[2]);
    for (label i = 0; i < nCutInside; i++)
    {
        const tetPointRef t = cutInsideTets[i].tet();
        t.sliceWithPlane(pl3, insideOp, outside);
    }
}
 
 
template<class tetsOp>
void Foam::tetOverlapVolume::cellCellOverlapMinDecomp
(
    const primitiveMesh& meshA,
    const label cellAI,
 
    const primitiveMesh& meshB,
    const label cellBI,
    const treeBoundBox& cellBbB,
    tetsOp& combineTetsOp
)
{
    const cell& cFacesA = meshA.cells()[cellAI];
    const point& ccA = meshA.cellCentres()[cellAI];
 
    const cell& cFacesB = meshB.cells()[cellBI];
    const point& ccB = meshB.cellCentres()[cellBI];
 
    forAll(cFacesA, cFA)
    {
        label faceAI = cFacesA[cFA];
 
        const face& fA = meshA.faces()[faceAI];
        const treeBoundBox pyrA = pyrBb(meshA.points(), fA, ccA);
        if (!pyrA.overlaps(cellBbB))
        {
            continue;
        }
 
        bool ownA = (meshA.faceOwner()[faceAI] == cellAI);
 
        label tetBasePtAI = 0;
 
        const point& tetBasePtA = meshA.points()[fA[tetBasePtAI]];
 
        for (label tetPtI = 1; tetPtI < fA.size() - 1; tetPtI++)
        {
            label facePtAI = (tetPtI + tetBasePtAI) % fA.size();
            label otherFacePtAI = fA.fcIndex(facePtAI);
 
            label pt0I = -1;
            label pt1I = -1;
 
            if (ownA)
            {
                pt0I = fA[facePtAI];
                pt1I = fA[otherFacePtAI];
            }
            else
            {
                pt0I = fA[otherFacePtAI];
                pt1I = fA[facePtAI];
            }
 
            const tetPoints tetA
            (
                ccA,
                tetBasePtA,
                meshA.points()[pt0I],
                meshA.points()[pt1I]
            );
            const treeBoundBox tetABb(tetA.bounds());
 
            // Loop over tets of cellB
            forAll(cFacesB, cFB)
            {
                label faceBI = cFacesB[cFB];
 
                const face& fB = meshB.faces()[faceBI];
                const treeBoundBox pyrB = pyrBb(meshB.points(), fB, ccB);
                if (!pyrB.overlaps(pyrA))
                {
                    continue;
                }
 
                bool ownB = (meshB.faceOwner()[faceBI] == cellBI);
 
                label tetBasePtBI = 0;
 
                const point& tetBasePtB = meshB.points()[fB[tetBasePtBI]];
 
                for (label tetPtI = 1; tetPtI < fB.size() - 1; tetPtI++)
                {
                    label facePtBI = (tetPtI + tetBasePtBI) % fB.size();
                    label otherFacePtBI = fB.fcIndex(facePtBI);
 
                    label pt0I = -1;
                    label pt1I = -1;
 
                    if (ownB)
                    {
                        pt0I = fB[facePtBI];
                        pt1I = fB[otherFacePtBI];
                    }
                    else
                    {
                        pt0I = fB[otherFacePtBI];
                        pt1I = fB[facePtBI];
                    }
 
                    const tetPoints tetB
                    (
                        ccB,
                        tetBasePtB,
                        meshB.points()[pt0I],
                        meshB.points()[pt1I]
                    );
                    if (!tetB.bounds().overlaps(tetABb))
                    {
                        continue;
                    }
 
                    if (combineTetsOp(tetA, tetB))
                    {
                        return;
                    }
                }
            }
        }
    }
}
 
 
// ************************************************************************* //