wallBoundedParticle.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
    Copyright (C) 2017-2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 "wallBoundedParticle.H"
 
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
const std::size_t Foam::wallBoundedParticle::sizeofFields_
(
    sizeof(wallBoundedParticle) - sizeof(particle)
);
 
 
// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
Foam::edge Foam::wallBoundedParticle::currentEdge() const
{
    if ((meshEdgeStart_ != -1) == (diagEdge_ != -1))
    {
        FatalErrorInFunction
            << "Particle:"
            << info()
            << "cannot both be on a mesh edge and a face-diagonal edge."
            << " meshEdgeStart_:" << meshEdgeStart_
            << " diagEdge_:" << diagEdge_
            << abort(FatalError);
    }
 
    const Foam::face& f = mesh().faces()[tetFace()];
 
    if (meshEdgeStart_ != -1)
    {
        return edge(f[meshEdgeStart_], f.nextLabel(meshEdgeStart_));
    }
    else
    {
        label faceBasePti = mesh().tetBasePtIs()[tetFace()];
        if (faceBasePti == -1)
        {
            //FatalErrorInFunction
            //WarningInFunction
            //    << "No base point for face " << tetFace() << ", " << f
            //    << ", produces a decomposition that has a minimum "
            //    << "volume greater than tolerance."
            //    //<< abort(FatalError);
            //    << endl;
            faceBasePti = 0;
        }
 
        label diagPti = (faceBasePti + diagEdge_)%f.size();
 
        return edge(f[faceBasePti], f[diagPti]);
    }
}
 
 
void Foam::wallBoundedParticle::crossEdgeConnectedFace
(
    const label& celli,
    label& tetFacei,
    label& tetPti,
    const edge& e
)
{
    const faceList& pFaces = mesh().faces();
    const cellList& pCells = mesh().cells();
 
    const Foam::face& f = pFaces[tetFacei];
 
    const Foam::cell& thisCell = pCells[celli];
 
    for (const label facei : thisCell)
    {
        // Loop over all other faces of this cell and
        // find the one that shares this edge
 
        if (tetFacei == facei)
        {
            continue;
        }
 
        const Foam::face& otherFace = pFaces[facei];
 
        label edDir = otherFace.edgeDirection(e);
 
        if (edDir == 0)
        {
            continue;
        }
        else if (f == pFaces[facei])
        {
            // This is a necessary condition if using duplicate baffles
            // (so coincident faces). We need to make sure we don't cross into
            // the face with the same vertices since we might enter a tracking
            // loop where it never exits. This test should be cheap
            // for most meshes so can be left in for 'normal' meshes.
            continue;
        }
        else
        {
            //Found edge on other face
            tetFacei = facei;
 
            label eIndex = -1;
 
            if (edDir == 1)
            {
                // Edge is in the forward circulation of this face, so
                // work with the start point of the edge
                eIndex = otherFace.find(e.start());
            }
            else
            {
                // edDir == -1, so the edge is in the reverse
                // circulation of this face, so work with the end
                // point of the edge
                eIndex = otherFace.find(e.end());
            }
 
            label tetBasePtI = mesh().tetBasePtIs()[facei];
 
            if (tetBasePtI == -1)
            {
                //FatalErrorInFunction
                //WarningInFunction
                //    << "No base point for face " << fI << ", " << f
                //    << ", produces a decomposition that has a minimum "
                //    << "volume greater than tolerance."
                //    //<< abort(FatalError);
                //    << endl;
                tetBasePtI = 0;
            }
 
            // Find eIndex relative to the base point on new face
            eIndex -= tetBasePtI;
 
            if (neg(eIndex))
            {
                eIndex = (eIndex + otherFace.size()) % otherFace.size();
            }
 
            if (eIndex == 0)
            {
                // The point is the base point, so this is first tet
                // in the face circulation
                tetPti = 1;
            }
            else if (eIndex == otherFace.size() - 1)
            {
                // The point is the last before the base point, so
                // this is the last tet in the face circulation
                tetPti = otherFace.size() - 2;
            }
            else
            {
                tetPti = eIndex;
            }
 
            break;
        }
    }
}
 
 
void Foam::wallBoundedParticle::crossEdgeConnectedFace(const edge& meshEdge)
{
    // Update tetFace, tetPt
    crossEdgeConnectedFace(cell(), tetFace(), tetPt(), meshEdge);
 
    // Update face to be same as tracking one
    face() = tetFace();
 
    // And adapt meshEdgeStart_.
    const Foam::face& f = mesh().faces()[tetFace()];
    label fp = f.find(meshEdge[0]);
 
    if (f.nextLabel(fp) == meshEdge[1])
    {
        meshEdgeStart_ = fp;
    }
    else
    {
        label fpMin1 = f.rcIndex(fp);
 
        if (f[fpMin1] == meshEdge[1])
        {
            meshEdgeStart_ = fpMin1;
        }
        else
        {
            FatalErrorInFunction
                << "Problem :"
                << " particle:"
                << info()
                << "face:" << tetFace()
                << " verts:" << f
                << " meshEdge:" << meshEdge
                << abort(FatalError);
        }
    }
 
    diagEdge_ = -1;
 
    // Check that still on same mesh edge
    const edge eNew(f[meshEdgeStart_], f.nextLabel(meshEdgeStart_));
    if (eNew != meshEdge)
    {
        FatalErrorInFunction
            << "Problem" << abort(FatalError);
    }
}
 
 
void Foam::wallBoundedParticle::crossDiagonalEdge()
{
    if (diagEdge_ == -1)
    {
        FatalErrorInFunction
            << "Particle:"
            << info()
            << "not on a diagonal edge" << abort(FatalError);
    }
    if (meshEdgeStart_ != -1)
    {
        FatalErrorInFunction
            << "Particle:"
            << info()
            << "meshEdgeStart_:" << meshEdgeStart_ << abort(FatalError);
    }
 
    const Foam::face& f = mesh().faces()[tetFace()];
 
    // tetPti starts from 1, goes up to f.size()-2
 
    if (tetPt() == diagEdge_)
    {
        tetPt() = f.rcIndex(tetPt());
    }
    else
    {
        label nextTetPt = f.fcIndex(tetPt());
        if (diagEdge_ == nextTetPt)
        {
            tetPt() = nextTetPt;
        }
        else
        {
            FatalErrorInFunction
                << "Particle:"
                << info()
                << "tetPt:" << tetPt()
                << " diagEdge:" << diagEdge_ << abort(FatalError);
        }
    }
 
    meshEdgeStart_ = -1;
}
 
 
Foam::scalar Foam::wallBoundedParticle::trackFaceTri
(
    const vector& n,
    const vector& endPosition,
    label& minEdgei
)
{
    // Track p from position to endPosition
    const triFace tri(currentTetIndices().faceTriIs(mesh(), false));
 
    // Check which edge intersects the trajectory.
    // Project trajectory onto triangle
    minEdgei = -1;
    scalar minS = 1;        // end position
 
    edge currentE(-1, -1);
    if (meshEdgeStart_ != -1 || diagEdge_ != -1)
    {
        currentE = currentEdge();
    }
 
    // Determine path along line position+s*d to see where intersections are.
    forAll(tri, i)
    {
        label j = tri.fcIndex(i);
 
        const point& pt0 = mesh().points()[tri[i]];
        const point& pt1 = mesh().points()[tri[j]];
 
        if (edge(tri[i], tri[j]) == currentE)
        {
            // Do not check particle is on
            continue;
        }
 
        // Outwards pointing normal
        vector edgeNormal = (pt1 - pt0)^n;
        edgeNormal.normalise();
 
        // Determine whether position and end point on either side of edge.
        scalar sEnd = (endPosition - pt0)&edgeNormal;
        if (sEnd >= 0)
        {
            // endPos is outside triangle. localPosition_ should always be
            // inside.
            scalar sStart = (localPosition_ - pt0)&edgeNormal;
            if (mag(sEnd - sStart) > VSMALL)
            {
                scalar s = sStart/(sStart - sEnd);
 
                if (s >= 0 && s < minS)
                {
                    minS = s;
                    minEdgei = i;
                }
            }
        }
    }
 
    if (minEdgei != -1)
    {
        localPosition_ += minS*(endPosition - localPosition_);
    }
    else
    {
        // Did not hit any edge so tracked to the end position. Set position
        // without any calculation to avoid truncation errors.
        localPosition_ = endPosition;
        minS = 1.0;
    }
 
    // Project localPosition_ back onto plane of triangle
    const point& triPt = mesh().points()[tri[0]];
    localPosition_ -= ((localPosition_ - triPt)&n)*n;
 
    return minS;
}
 
 
bool Foam::wallBoundedParticle::isTriAlongTrack
(
    const vector& n,
    const point& endPosition
) const
{
    const triFace triVerts(currentTetIndices().faceTriIs(mesh(), false));
    const edge currentE = currentEdge();
 
    if
    (
        currentE[0] == currentE[1]
     || !triVerts.found(currentE[0])
     || !triVerts.found(currentE[1])
    )
    {
        FatalErrorInFunction
            << "Edge " << currentE << " not on triangle " << triVerts
            << info()
            << abort(FatalError);
    }
 
 
    const vector dir = endPosition - localPosition_;
 
    forAll(triVerts, i)
    {
        label j = triVerts.fcIndex(i);
        const point& pt0 = mesh().points()[triVerts[i]];
        const point& pt1 = mesh().points()[triVerts[j]];
 
        if (edge(triVerts[i], triVerts[j]) == currentE)
        {
            vector edgeNormal = (pt1-pt0)^n;
            return (dir&edgeNormal) < 0;
        }
    }
 
    FatalErrorInFunction
        << "Problem" << abort(FatalError);
 
    return false;
}
 
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
Foam::wallBoundedParticle::wallBoundedParticle
(
    const polyMesh& mesh,
    const point& position,
    const label celli,
    const label tetFacei,
    const label tetPti,
    const label meshEdgeStart,
    const label diagEdge
)
:
    particle(mesh, position, celli, tetFacei, tetPti, false),
    localPosition_(position),
    meshEdgeStart_(meshEdgeStart),
    diagEdge_(diagEdge)
{}
 
 
Foam::wallBoundedParticle::wallBoundedParticle
(
    const polyMesh& mesh,
    Istream& is,
    bool readFields,
    bool newFormat
)
:
    particle(mesh, is, readFields, newFormat)
{
    if (readFields)
    {
        if (is.format() == IOstream::ASCII)
        {
            is  >> localPosition_ >> meshEdgeStart_ >> diagEdge_;
        }
        if (!is.checkLabelSize<>() || !is.checkScalarSize<>())
        {
            // Non-native label or scalar size
 
            is.beginRawRead();
 
            readRawScalar(is, localPosition_.data(), vector::nComponents);
            readRawLabel(is, &meshEdgeStart_);
            readRawLabel(is, &diagEdge_);
 
            is.endRawRead();
        }
        else
        {
            is.read(reinterpret_cast<char*>(&localPosition_), sizeofFields_);
        }
    }
 
    is.check(FUNCTION_NAME);
}
 
 
Foam::wallBoundedParticle::wallBoundedParticle
(
    const wallBoundedParticle& p
)
:
    particle(p),
    localPosition_(p.localPosition_),
    meshEdgeStart_(p.meshEdgeStart_),
    diagEdge_(p.diagEdge_)
{}
 
 
// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
Foam::Ostream& Foam::operator<<
(
    Ostream& os,
    const wallBoundedParticle& p
)
{
    if (os.format() == IOstream::ASCII)
    {
        os  << static_cast<const particle&>(p)
        << token::SPACE << p.localPosition_
        << token::SPACE << p.meshEdgeStart_
        << token::SPACE << p.diagEdge_;
    }
    else
    {
        os  << static_cast<const particle&>(p);
        os.write
        (
            reinterpret_cast<const char*>(&p.localPosition_),
            wallBoundedParticle::sizeofFields_
        );
    }
 
    return os;
}
 
 
Foam::Ostream& Foam::operator<<
(
    Ostream& os,
    const InfoProxy<wallBoundedParticle>& ip
)
{
    const wallBoundedParticle& p = ip.t_;
 
    tetPointRef tpr(p.currentTetIndices().tet(p.mesh()));
 
    os  << "    " << static_cast<const particle&>(p) << nl
        << "    tet:" << nl;
    os  << "    ";
    meshTools::writeOBJ(os, tpr.a());
    os  << "    ";
    meshTools::writeOBJ(os, tpr.b());
    os  << "    ";
    meshTools::writeOBJ(os, tpr.c());
    os  << "    ";
    meshTools::writeOBJ(os, tpr.d());
    os  << "    l 1 2" << nl
        << "    l 1 3" << nl
        << "    l 1 4" << nl
        << "    l 2 3" << nl
        << "    l 2 4" << nl
        << "    l 3 4" << nl;
    os  << "    ";
    meshTools::writeOBJ(os, p.localPosition_);
 
    return os;
}
 
 
// ************************************************************************* //