ENH: wallBoundedStreamline: survive neg-tet cells. Fixes #502.

src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,6 +35,57 @@ const std::size_t Foam::wallBoundedParticle::sizeofFields_
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
+Foam::tetIndices Foam::wallBoundedParticle::currentTetIndices() const
+{
+    // Replacement for particle::currentTetIndices that avoids error
+    // upon invalid tetBasePtIs
+
+    const faceList& pFaces = mesh_.faces();
+    const labelList& pOwner = mesh_.faceOwner();
+
+    const Foam::face& f = pFaces[tetFacei_];
+
+    bool own = (pOwner[tetFacei_] == celli_);
+
+    label faceBasePtI = mesh_.tetBasePtIs()[tetFacei_];
+    if (faceBasePtI == -1)
+    {
+        //WarningInFunction
+        //    << "No base point for face " << tetFacei_ << ", " << f
+        //    << ", produces a decomposition that has a minimum "
+        //    << "volume greater than tolerance."
+        //    << endl;
+        faceBasePtI = 0;
+    }
+
+    label facePtI = (tetPti_ + faceBasePtI) % f.size();
+    label otherFacePtI = f.fcIndex(facePtI);
+
+    label facePtAI;
+    label facePtBI;
+    if (own)
+    {
+        facePtAI = facePtI;
+        facePtBI = otherFacePtI;
+    }
+    else
+    {
+        facePtAI = otherFacePtI;
+        facePtBI = facePtI;
+    }
+
+    return tetIndices
+    (
+        celli_,
+        tetFacei_,
+        faceBasePtI,
+        facePtAI,
+        facePtBI,
+        tetPti_
+    );
+}
+
+
 Foam::edge Foam::wallBoundedParticle::currentEdge() const
 {
     if ((meshEdgeStart_ != -1) == (diagEdge_ != -1))
@@ -57,6 +108,18 @@ Foam::edge Foam::wallBoundedParticle::currentEdge() const
     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]);
@@ -66,11 +129,118 @@ Foam::edge Foam::wallBoundedParticle::currentEdge() const
 
 void Foam::wallBoundedParticle::crossEdgeConnectedFace
 (
-    const edge& meshEdge
+    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];
+
+    forAll(thisCell, cFI)
+    {
+        // Loop over all other faces of this cell and
+        // find the one that shares this edge
+
+        label fI = thisCell[cFI];
+
+        if (tetFacei == fI)
+        {
+            continue;
+        }
+
+        const Foam::face& otherFace = pFaces[fI];
+
+        label edDir = otherFace.edgeDirection(e);
+
+        if (edDir == 0)
+        {
+            continue;
+        }
+        else if (f == pFaces[fI])
+        {
+            // 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 = fI;
+
+            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 = findIndex(otherFace, 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 = findIndex(otherFace, e.end());
+            }
+
+            label tetBasePtI = mesh_.tetBasePtIs()[fI];
+
+            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
-    particle::crossEdgeConnectedFace(cell(), tetFace(), tetPt(), meshEdge);
+    crossEdgeConnectedFace(cell(), tetFace(), tetPt(), meshEdge);
 
     // Update face to be same as tracking one
     face() = tetFace();
@@ -164,14 +334,13 @@ void Foam::wallBoundedParticle::crossDiagonalEdge()
 
 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_));
-    vector n = tri.normal(mesh_.points());
-    n /= mag(n)+VSMALL;
 
     // Check which edge intersects the trajectory.
     // Project trajectory onto triangle
@@ -245,6 +414,7 @@ Foam::scalar Foam::wallBoundedParticle::trackFaceTri
 
 bool Foam::wallBoundedParticle::isTriAlongTrack
 (
+    const vector& n,
     const point& endPosition
 ) const
 {
@@ -267,10 +437,6 @@ bool Foam::wallBoundedParticle::isTriAlongTrack
 
     const vector dir = endPosition-position();
 
-    // Get normal of currentE
-    vector n = triVerts.normal(mesh_.points());
-    n /= mag(n);
-
     forAll(triVerts, i)
     {
         label j = triVerts.fcIndex(i);
@@ -382,7 +548,7 @@ Foam::Ostream& Foam::operator<<
 {
     const wallBoundedParticle& p = ip.t_;
 
-    tetPointRef tpr(p.currentTet());
+    tetPointRef tpr(p.currentTetIndices().tet(p.mesh()));
 
     os  << "    " << static_cast<const particle&>(p) << nl
         << "    tet:" << nl;

src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -123,6 +123,20 @@ protected:
         //- Construct current edge
         edge currentEdge() const;
 
+        //- Replacement for particle::currentTetIndices() that avoids bombing
+        //  out on invalid tet decomposition (tetBasePtIs = -1)
+        tetIndices currentTetIndices() const;
+
+        //- Replacement for particle::crossEdgeConnectedFace that avoids bombing
+        //  out on invalid tet decomposition (tetBasePtIs = -1)
+        void crossEdgeConnectedFace
+        (
+            const label& celli,
+            label& tetFacei,
+            label& tetPti,
+            const edge& e
+        );
+
         //- Cross mesh edge into different face on same cell
         void crossEdgeConnectedFace(const edge& meshEdge);
 
@@ -130,10 +144,10 @@ protected:
         void crossDiagonalEdge();
 
         //- Track through single triangle
-        scalar trackFaceTri(const vector& endPosition, label& minEdgei);
+        scalar trackFaceTri(const vector& n, const vector& endPosition, label&);
 
         //- Is current triangle in the track direction
-        bool isTriAlongTrack(const point& endPosition) const;
+        bool isTriAlongTrack(const vector& n, const point& endPosition) const;
 
 
     // Patch interactions

src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -175,7 +175,14 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge
             // endposition? i.e. since volume of nbr tet is positive the
             // tracking direction should be into the tet.
             tetIndices nbrTi(nbrCelli, tetFacei_, tetPti_, mesh_);
-            if ((nbrTi.faceTri(mesh_).normal() & (endPosition-position())) < 0)
+
+            bool posVol = (nbrTi.tet(mesh_).mag() > 0);
+
+            if
+            (
+                posVol
+             == ((nbrTi.faceTri(mesh_).normal() & (endPosition-position())) < 0)
+            )
             {
                 // Change into nbrCell. No need to change tetFace, tetPt.
                 //Pout<< "    crossed from cell:" << celli_
@@ -213,12 +220,22 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge
                 << abort(FatalError);
         }
 
+        const triFace tri(currentTetIndices().faceTriIs(mesh_));
+        vector n = tri.normal(mesh_.points());
+        n /= mag(n);
+
         point projectedEndPosition = endPosition;
+
+        const bool posVol = (currentTetIndices().tet(mesh_).mag() > 0);
+
+        if (!posVol)
+        {
+            // Negative tet volume. Track back by setting the end point
+            projectedEndPosition = position() - (endPosition-position());
+        }
+
         // Remove normal component
         {
-            const triFace tri(currentTetIndices().faceTriIs(mesh_));
-            vector n = tri.normal(mesh_.points());
-            n /= mag(n);
             const point& basePt = mesh_.points()[tri[0]];
             projectedEndPosition -= ((projectedEndPosition-basePt)&n)*n;
         }
@@ -234,7 +251,7 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge
         {
             // See if the current triangle has got a point on the
             // correct side of the edge.
-            doTrack = isTriAlongTrack(projectedEndPosition);
+            doTrack = isTriAlongTrack(n, projectedEndPosition);
         }
 
 
@@ -242,7 +259,7 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge
         {
             // Track across triangle. Return triangle edge crossed.
             label triEdgei = -1;
-            trackFraction = trackFaceTri(projectedEndPosition, triEdgei);
+            trackFraction = trackFaceTri(n, projectedEndPosition, triEdgei);
 
             if (triEdgei == -1)
             {

src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -120,11 +120,12 @@ Foam::vector Foam::wallBoundedStreamLineParticle::sample
     {
         // Stagnant particle. Might as well stop
         lifeTime_ = 0;
+        return vector::zero;
+    }
+    else
+    {
+        return U/magU;
     }
-
-    U /= magU;
-
-    return U;
 }
 
 
@@ -267,16 +268,6 @@ bool Foam::wallBoundedStreamLineParticle::move
             // Force removal
             lifeTime_ = 0;
         }
-
-        if
-        (
-            !td.keepParticle
-        ||  td.switchProcessor
-        ||  lifeTime_ == 0
-        )
-        {
-            break;
-        }
     }