From 46edbe2bb2c7afa54270536f730be9a6bb7c0b1e Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 14 Feb 2012 12:21:47 +0000
Subject: [PATCH] ENH: wallBoundedParticle: generic particle class for tracking
 along boundary

---
 .../wallBoundedParticle.C                     | 659 ++++++++++++++++++
 .../wallBoundedParticle.H                     | 345 +++++++++
 .../wallBoundedParticleTemplates.C            | 543 +++++++++++++++
 3 files changed, 1547 insertions(+)
 create mode 100644 src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C
 create mode 100644 src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H
 create mode 100644 src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C

diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C
new file mode 100644
index 0000000000..e73f4e72f4
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C
@@ -0,0 +1,659 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2012 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 "wallBoundedParticle.H"
+#include "vectorFieldIOField.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+//    defineParticleTypeNameAndDebug(wallBoundedParticle, 0);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+//// Check position is inside tet
+//void Foam::wallBoundedParticle::checkInside() const
+//{
+//    const tetIndices ti(currentTetIndices());
+//    const tetPointRef tpr(ti.tet(mesh_));
+//    if (!tpr.inside(position()))
+//    {
+//        FatalErrorIn("wallBoundedParticle::checkInside(..)")
+//            << "Particle:" //<< static_cast<const particle&>(*this)
+//            << info()
+//            << "is not inside " << tpr
+//            << abort(FatalError);
+//    }
+//}
+//
+//
+//void Foam::wallBoundedParticle::checkOnEdge() const
+//{
+//    // Check that edge (as indicated by meshEdgeStart_, diagEdge_) is
+//    // indeed one that contains the position.
+//    const edge e = currentEdge();
+//
+//    linePointRef ln(e.line(mesh_.points()));
+//
+//    pointHit ph(ln.nearestDist(position()));
+//
+//    if (ph.distance() > 1E-6)
+//    {
+//        FatalErrorIn
+//        (
+//            "wallBoundedParticle::checkOnEdge()"
+//        )   << "Problem :"
+//            << " particle:" //<< static_cast<const particle&>(*this)
+//            << info()
+//            << "edge:" << e
+//            << " at:" << ln
+//            << " distance:" << ph.distance()
+//            << abort(FatalError);
+//    }
+//}
+//
+//
+//void Foam::wallBoundedParticle::checkOnTriangle(const point& p)
+//const
+//{
+//    const triFace tri(currentTetIndices().faceTriIs(mesh_));
+//    pointHit ph = tri.nearestPoint(p, mesh_.points());
+//    if (ph.distance() > 1e-9)
+//    {
+//        FatalErrorIn
+//        (
+//            "wallBoundedParticle::checkOnTriangle(const point&)"
+//        )   << "Problem :"
+//            << " particle:" //<< static_cast<const particle&>(*this)
+//            << info()
+//            << "point:" << p
+//            << " distance:" << ph.distance()
+//            << abort(FatalError);
+//    }
+//}
+
+
+// Construct the edge the particle is on (according to meshEdgeStart_,
+// diagEdge_)
+Foam::edge Foam::wallBoundedParticle::currentEdge() const
+{
+    if ((meshEdgeStart_ != -1) == (diagEdge_ != -1))
+    {
+        FatalErrorIn("wallBoundedParticle::currentEdge() const")
+            << "Particle:" //<< static_cast<const particle&>(*this)
+            << 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()];
+        label diagPtI = (faceBasePtI+diagEdge_)%f.size();
+        return edge(f[faceBasePtI], f[diagPtI]);
+    }
+}
+
+
+void Foam::wallBoundedParticle::crossEdgeConnectedFace
+(
+    const edge& meshEdge
+)
+{
+    //label oldFaceI = tetFace();
+
+    // Update tetFace, tetPt
+    particle::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 = findIndex(f, meshEdge[0]);
+
+    if (f.nextLabel(fp) == meshEdge[1])
+    {
+        meshEdgeStart_ = fp;
+    }
+    else
+    {
+        label fpMin1 = f.rcIndex(fp);
+
+        if (f[fpMin1] == meshEdge[1])
+        {
+            meshEdgeStart_ = fpMin1;
+        }
+        else
+        {
+            FatalErrorIn
+            (
+                "wallBoundedParticle::crossEdgeConnectedFace"
+                "(const edge&)"
+            )   << "Problem :"
+                << " particle:" //<< static_cast<const particle&>(*this)
+                << info()
+                << "face:" << tetFace()
+                << " verts:" << f
+                << " meshEdge:" << meshEdge
+                << abort(FatalError);
+        }
+    }
+
+    diagEdge_ = -1;
+
+    //Pout<< "    crossed meshEdge "
+    //    << meshEdge.line(mesh().points())
+    //    << " from face:" << oldFaceI
+    //    << " to face:" << tetFace() << endl;
+
+
+    // Check that still on same mesh edge
+
+    const edge eNew(f[meshEdgeStart_], f.nextLabel(meshEdgeStart_));
+    if (eNew != meshEdge)
+    {
+        FatalErrorIn
+        (
+            "wallBoundedParticle::crossEdgeConnectedFace"
+            "(const edge&)"
+        )   << "Problem" << abort(FatalError);
+    }
+
+
+    // Check that edge (as indicated by meshEdgeStart_) is indeed one that
+    // contains the position.
+    //checkOnEdge();
+}
+
+
+void Foam::wallBoundedParticle::crossDiagonalEdge()
+{
+    if (diagEdge_ == -1)
+    {
+        FatalErrorIn("wallBoundedParticle::crossDiagonalEdge()")
+            << "Particle:" //<< static_cast<const particle&>(*this)
+            << info()
+            << "not on a diagonal edge" << abort(FatalError);
+    }
+    if (meshEdgeStart_ != -1)
+    {
+        FatalErrorIn("wallBoundedParticle::crossDiagonalEdge()")
+            << "Particle:" //<< static_cast<const particle&>(*this)
+            << info()
+            << "meshEdgeStart_:" << meshEdgeStart_ << abort(FatalError);
+    }
+
+    //label oldTetPt = tetPt();
+
+    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
+        {
+            FatalErrorIn("wallBoundedParticle::crossDiagonalEdge()")
+                << "Particle:" //<< static_cast<const particle&>(*this)
+                << info()
+                << "tetPt:" << tetPt()
+                << " diagEdge:" << diagEdge_ << abort(FatalError);
+        }
+    }
+
+    meshEdgeStart_ = -1;
+
+    //Pout<< "    crossed diagonalEdge "
+    //    << currentEdge().line(mesh().points())
+    //    << " from tetPt:" << oldTetPt
+    //    << " to tetPt:" << tetPt() << endl;
+}
+
+
+//- Track through a single triangle.
+// Gets passed tet+triangle the particle is in. Updates position() but nothing
+// else. Returns the triangle edge the particle is now on.
+Foam::scalar Foam::wallBoundedParticle::trackFaceTri
+(
+    const vector& endPosition,
+    label& minEdgeI
+)
+{
+    // Track p from position to endPosition
+    const triFace tri(currentTetIndices().faceTriIs(mesh_));
+    vector n = tri.normal(mesh_.points());
+    //if (mag(n) < sqr(SMALL))
+    //{
+    //    FatalErrorIn("wallBoundedParticle::trackFaceTri(..)")
+    //        << "Small triangle." //<< static_cast<const particle&>(*this)
+    //        << info()
+    //        << "n:" << n
+    //        << abort(FatalError);
+    //}
+    n /= mag(n)+VSMALL;
+
+    // Check which edge intersects the trajectory.
+    // Project trajectory onto triangle
+    minEdgeI = -1;
+    scalar minS = 1;        // end position
+
+    //const point oldPosition(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;
+
+        //if (mag(edgeNormal) < SMALL)
+        //{
+        //    FatalErrorIn("wallBoundedParticle::trackFaceTri(..)")
+        //        << "Edge not perpendicular to triangle."
+        //        //<< static_cast<const particle&>(*this)
+        //        << info()
+        //        << "triangle n:" << n
+        //        << " edgeNormal:" << edgeNormal
+        //        << " on tri:" << tri
+        //        << " at:" << pt0
+        //        << " at:" << pt1
+        //        << abort(FatalError);
+        //}
+
+
+        edgeNormal /= mag(edgeNormal)+VSMALL;
+
+        // Determine whether position and end point on either side of edge.
+        scalar sEnd = (endPosition-pt0)&edgeNormal;
+        if (sEnd >= 0)
+        {
+            // endPos is outside triangle. position() should always be
+            // inside.
+            scalar sStart = (position()-pt0)&edgeNormal;
+            if (mag(sEnd-sStart) > VSMALL)
+            {
+                scalar s = sStart/(sStart-sEnd);
+
+                if (s >= 0 && s < minS)
+                {
+                    minS = s;
+                    minEdgeI = i;
+                }
+            }
+        }
+    }
+
+    if (minEdgeI != -1)
+    {
+        position() += minS*(endPosition-position());
+    }
+    else
+    {
+        // Did not hit any edge so tracked to the end position. Set position
+        // without any calculation to avoid truncation errors.
+        position() = endPosition;
+        minS = 1.0;
+    }
+
+    // Project position() back onto plane of triangle
+    const point& triPt = mesh_.points()[tri[0]];
+    position() -= ((position()-triPt)&n)*n;
+
+
+    //Pout<< "    tracked from:" << oldPosition << " to:" << position()
+    //    << " projectedEnd:" << endPosition
+    //    << " at s:" << minS << endl;
+    //if (minEdgeI != -1)
+    //{
+    //    Pout<< "    on edge:" << minEdgeI
+    //        << " on edge:"
+    //        << mesh_.points()[tri[minEdgeI]]
+    //        << mesh_.points()[tri[tri.fcIndex(minEdgeI)]]
+    //        << endl;
+    //}
+
+    return minS;
+}
+
+
+// See if the current triangle has got a point on the
+// correct side of the edge.
+bool Foam::wallBoundedParticle::isTriAlongTrack
+(
+    const point& endPosition
+) const
+{
+    const triFace triVerts(currentTetIndices().faceTriIs(mesh_));
+    const edge currentE = currentEdge();
+
+    //if (debug)
+    {
+        if
+        (
+            currentE[0] == currentE[1]
+         || findIndex(triVerts, currentE[0]) == -1
+         || findIndex(triVerts, currentE[1]) == -1
+        )
+        {
+            FatalErrorIn
+            (
+                "wallBoundedParticle::isTriAlongTrack"
+                "(const point&)"
+            )   << "Edge " << currentE << " not on triangle " << triVerts
+                << info()
+                << abort(FatalError);
+        }
+    }
+
+
+    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);
+        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;
+        }
+    }
+
+    FatalErrorIn
+    (
+        "wallBoundedParticle::isTriAlongTrack"
+        "(const point&)"
+    )   << "Problem" << abort(FatalError);
+
+    return false;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::wallBoundedParticle::wallBoundedParticle
+(
+    const polyMesh& mesh,
+    const vector& position,
+    const label cellI,
+    const label tetFaceI,
+    const label tetPtI,
+    const label meshEdgeStart,
+    const label diagEdge
+)
+:
+    particle(mesh, position, cellI, tetFaceI, tetPtI),
+    meshEdgeStart_(meshEdgeStart),
+    diagEdge_(diagEdge)
+{
+    //checkInside();
+
+    //if (meshEdgeStart_ != -1 || diagEdge_ != -1)
+    //{
+    //    checkOnEdge();
+    //}
+
+    // Unfortunately have no access to trackdata so cannot check if particle
+    // is on a wallPatch or has an mesh edge set (either of which is
+    // a requirement).
+}
+
+
+Foam::wallBoundedParticle::wallBoundedParticle
+(
+    const polyMesh& mesh,
+    Istream& is,
+    bool readFields
+)
+:
+    particle(mesh, is, readFields)
+{
+    if (readFields)
+    {
+        is  >> meshEdgeStart_ >> diagEdge_;
+    }
+
+    // Check state of Istream
+    is.check
+    (
+        "wallBoundedParticle::wallBoundedParticle"
+        "(const Cloud<wallBoundedParticle>&, Istream&, bool)"
+    );
+}
+
+
+Foam::wallBoundedParticle::wallBoundedParticle
+(
+    const wallBoundedParticle& p
+)
+:
+    particle(p),
+    meshEdgeStart_(p.meshEdgeStart_),
+    diagEdge_(p.diagEdge_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+//void Foam::wallBoundedParticle::readFields
+//(
+//    Cloud<wallBoundedParticle>& c
+//)
+//{
+//    if (!c.size())
+//    {
+//        return;
+//    }
+//
+//    particle::readFields(c);
+//
+//    IOField<label> meshEdgeStart
+//    (
+//        c.fieldIOobject("meshEdgeStart", IOobject::MUST_READ)
+//    );
+//
+//    IOField<label> diagEdge
+//    (
+//        c.fieldIOobject("diagEdge_", IOobject::MUST_READ)
+//    );
+//    c.checkFieldIOobject(c, diagEdge);
+//
+//    label i = 0;
+//    forAllIter(Cloud<wallBoundedParticle>, c, iter)
+//    {
+//        iter().meshEdgeStart_ = meshEdgeStart[i];
+//        iter().diagEdge_ = diagEdge[i];
+//        i++;
+//    }
+//}
+//
+//
+//void Foam::wallBoundedParticle::writeFields
+//(
+//    const Cloud<wallBoundedParticle>& c
+//)
+//{
+//    particle::writeFields(c);
+//
+//    label np =  c.size();
+//
+//    IOField<label> meshEdgeStart
+//    (
+//        c.fieldIOobject("meshEdgeStart", IOobject::NO_READ),
+//        np
+//    );
+//    IOField<label> diagEdge
+//    (
+//        c.fieldIOobject("diagEdge", IOobject::NO_READ),
+//        np
+//    );
+//
+//    label i = 0;
+//    forAllConstIter(Cloud<wallBoundedParticle>, c, iter)
+//    {
+//        meshEdgeStart[i] = iter().meshEdgeStart_;
+//        diagEdge[i] = iter().diagEdge_;
+//        i++;
+//    }
+//
+//    meshEdgeStart.write();
+//    diagEdge.write();
+//}
+
+void Foam::wallBoundedParticle::write(Ostream& os, bool writeFields) const
+{
+    const particle& p = static_cast<const particle&>(*this);
+
+    if (os.format() == IOstream::ASCII)
+    {
+        // Write base particle
+        p.write(os, writeFields);
+
+        if (writeFields)
+        {
+            // Write the additional entries
+            os  << token::SPACE << meshEdgeStart_
+                << token::SPACE << diagEdge_;
+        }
+    }
+    else
+    {
+        // Write base particle
+        p.write(os, writeFields);
+
+        // Write additional entries
+        if (writeFields)
+        {
+            os.write
+            (
+                reinterpret_cast<const char*>(&meshEdgeStart_),
+                sizeof(meshEdgeStart_)
+              + sizeof(diagEdge_)
+            );
+        }
+    }
+
+    // Check state of Ostream
+    os.check("wallBoundedParticle::write(Ostream& os, bool) const");
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const wallBoundedParticle& p
+)
+{
+    // Write all data
+    p.write(os, true);
+
+    return os;
+}
+
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const InfoProxy<wallBoundedParticle>& ip
+)
+{
+    const wallBoundedParticle& p = ip.t_;
+
+    tetPointRef tpr(p.currentTet());
+
+    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.position());
+
+    return os;
+}
+
+
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H
new file mode 100644
index 0000000000..df93215be1
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H
@@ -0,0 +1,345 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2012 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/>.
+
+Class
+    Foam::wallBoundedParticle
+
+Description
+    Particle class that tracks on triangles of boundary faces. Use
+    trackToEdge similar to trackToFace on particle.
+
+SourceFiles
+    wallBoundedParticle.C
+    wallBoundedParticleTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef wallBoundedParticle_H
+#define wallBoundedParticle_H
+
+#include "particle.H"
+#include "autoPtr.H"
+#include "InfoProxy.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class wallBoundedParticle Declaration
+\*---------------------------------------------------------------------------*/
+
+class wallBoundedParticle
+:
+    public particle
+{
+
+public:
+
+    //- Class used to pass tracking data to the trackToFace function
+    template<class CloudType>
+    class TrackingData
+    :
+        public particle::TrackingData<CloudType>
+    {
+
+    public:
+
+        const PackedBoolList& isWallPatch_;
+
+        // Constructors
+
+            inline TrackingData
+            (
+                CloudType& cloud,
+                const PackedBoolList& isWallPatch
+            )
+            :
+                particle::TrackingData<CloudType>
+                (
+                    cloud
+                ),
+                isWallPatch_(isWallPatch)
+            {}
+    };
+
+
+protected:
+
+    // Protected data
+
+        //- Particle is on mesh edge:
+        //      const face& f = mesh.faces()[tetFace()]
+        //      const edge e(f[meshEdgeStart_], f.nextLabel(meshEdgeStart_));
+        //  Note that this real edge
+        //  is also one of the edges of the face-triangle (from
+        //  tetFace()+tetPt()).
+        label meshEdgeStart_;
+
+        //- Particle is on diagonal edge:
+        //      const face& f = mesh.faces()[tetFace()]
+        //      label faceBasePtI = mesh.tetBasePtIs()[faceI];
+        //      label diagPtI = (faceBasePtI+diagEdge_)%f.size();
+        //      const edge e(f[faceBasePtI], f[diagPtI]);
+        label diagEdge_;
+
+
+    // Protected Member Functions
+
+        //- Construct current edge
+        edge currentEdge() const;
+
+        //- Check if inside current tet
+        //void checkInside() const;
+
+        //- Check if on current edge
+        //void checkOnEdge() const;
+
+        //- Check if point on triangle
+        //void checkOnTriangle(const point&) const;
+
+        //- Cross mesh edge into different face on same cell
+        void crossEdgeConnectedFace(const edge& meshEdge);
+
+        //- Cross diagonal edge into different triangle on same face,cell
+        void crossDiagonalEdge();
+
+        //- Track through single triangle
+        scalar trackFaceTri(const vector& endPosition, label& minEdgeI);
+
+        //- Is current triangle in the track direction
+        bool isTriAlongTrack(const point& endPosition) const;
+
+
+    // Patch interactions
+
+        //- Do all patch interaction
+        template<class TrackData>
+        void patchInteraction(TrackData& td, const scalar trackFraction);
+
+        //- Overridable function to handle the particle hitting a patch
+        //  Executed before other patch-hitting functions
+        template<class TrackData>
+        bool hitPatch
+        (
+            const polyPatch&,
+            TrackData& td,
+            const label patchI,
+            const scalar trackFraction,
+            const tetIndices& tetIs
+        );
+
+        //- Overridable function to handle the particle hitting a wedge
+        template<class TrackData>
+        void hitWedgePatch
+        (
+            const wedgePolyPatch&,
+            TrackData& td
+        );
+
+        //- Overridable function to handle the particle hitting a
+        //  symmetryPlane
+        template<class TrackData>
+        void hitSymmetryPatch
+        (
+            const symmetryPolyPatch&,
+            TrackData& td
+        );
+
+        //- Overridable function to handle the particle hitting a cyclic
+        template<class TrackData>
+        void hitCyclicPatch
+        (
+            const cyclicPolyPatch&,
+            TrackData& td
+        );
+
+        //- Overridable function to handle the particle hitting a
+        //- processorPatch
+        template<class TrackData>
+        void hitProcessorPatch
+        (
+            const processorPolyPatch&,
+            TrackData& td
+        );
+
+        //- Overridable function to handle the particle hitting a wallPatch
+        template<class TrackData>
+        void hitWallPatch
+        (
+            const wallPolyPatch&,
+            TrackData& td,
+            const tetIndices&
+        );
+
+        //- Overridable function to handle the particle hitting a polyPatch
+        template<class TrackData>
+        void hitPatch
+        (
+            const polyPatch&,
+            TrackData& td
+        );
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        wallBoundedParticle
+        (
+            const polyMesh& c,
+            const vector& position,
+            const label cellI,
+            const label tetFaceI,
+            const label tetPtI,
+            const label meshEdgeStart,
+            const label diagEdge
+        );
+
+        //- Construct from Istream
+        wallBoundedParticle
+        (
+            const polyMesh& c,
+            Istream& is,
+            bool readFields = true
+        );
+
+        //- Construct copy
+        wallBoundedParticle(const wallBoundedParticle& p);
+
+        //- Construct and return a clone
+        autoPtr<particle> clone() const
+        {
+            return autoPtr<particle>(new wallBoundedParticle(*this));
+        }
+
+        //- Factory class to read-construct particles used for
+        //  parallel transfer
+        class iNew
+        {
+            const polyMesh& mesh_;
+
+        public:
+
+            iNew(const polyMesh& mesh)
+            :
+                mesh_(mesh)
+            {}
+
+            autoPtr<wallBoundedParticle> operator()
+            (
+                Istream& is
+            ) const
+            {
+                return autoPtr<wallBoundedParticle>
+                (
+                    new wallBoundedParticle(mesh_, is, true)
+                );
+            }
+        };
+
+
+    // Member Functions
+
+        // Access
+
+            //- -1 or label of mesh edge
+            inline label meshEdgeStart() const
+            {
+                return meshEdgeStart_;
+            }
+
+            //- -1 or diagonal edge
+            inline label diagEdge() const
+            {
+                return diagEdge_;
+            }
+
+
+        // Track
+
+            //- Equivalent of trackToFace
+            template<class TrackData>
+            scalar trackToEdge
+            (
+                TrackData& td,
+                const vector& endPosition
+            );
+
+
+        // Info
+
+            //- Return info proxy.
+            //  Used to print particle information to a stream
+            inline InfoProxy<wallBoundedParticle> info() const
+            {
+                return *this;
+            }
+
+
+        // I-O
+
+            //- Read
+            template<class CloudType>
+            static void readFields(CloudType&);
+
+            //- Write
+            template<class CloudType>
+            static void writeFields(const CloudType&);
+
+            //- Write the particle data
+            void write(Ostream& os, bool writeFields) const;
+
+
+    // Ostream Operator
+
+        friend Ostream& operator<<
+        (
+            Ostream&,
+            const wallBoundedParticle&
+        );
+
+        friend Ostream& operator<<
+        (
+            Ostream&,
+            const InfoProxy<wallBoundedParticle>&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "wallBoundedParticleTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C
new file mode 100644
index 0000000000..f79d964c7b
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C
@@ -0,0 +1,543 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2012 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 "wallBoundedParticle.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class TrackData>
+void Foam::wallBoundedParticle::patchInteraction
+(
+    TrackData& td,
+    const scalar trackFraction
+)
+{
+//    typedef TrackData::CloudType cloudType;
+    typedef typename TrackData::cloudType::particleType particleType;
+
+    particleType& p = static_cast<particleType&>(*this);
+    p.hitFace(td);
+
+    if (!internalFace(faceI_))
+    {
+        label origFaceI = faceI_;
+        label patchI = patch(faceI_);
+
+        // No action taken for tetPtI_ for tetFaceI_ here, handled by
+        // patch interaction call or later during processor transfer.
+
+
+        // Dummy tet indices. What to do here?
+        tetIndices faceHitTetIs;
+
+        if
+        (
+            !p.hitPatch
+            (
+                mesh_.boundaryMesh()[patchI],
+                td,
+                patchI,
+                trackFraction,
+                faceHitTetIs
+            )
+        )
+        {
+            // Did patch interaction model switch patches?
+            // Note: recalculate meshEdgeStart_, diagEdge_!
+            if (faceI_ != origFaceI)
+            {
+                patchI = patch(faceI_);
+            }
+
+            const polyPatch& patch = mesh_.boundaryMesh()[patchI];
+
+            if (isA<wedgePolyPatch>(patch))
+            {
+                p.hitWedgePatch
+                (
+                    static_cast<const wedgePolyPatch&>(patch), td
+                );
+            }
+            else if (isA<symmetryPolyPatch>(patch))
+            {
+                p.hitSymmetryPatch
+                (
+                    static_cast<const symmetryPolyPatch&>(patch), td
+                );
+            }
+            else if (isA<cyclicPolyPatch>(patch))
+            {
+                p.hitCyclicPatch
+                (
+                    static_cast<const cyclicPolyPatch&>(patch), td
+                );
+            }
+            else if (isA<processorPolyPatch>(patch))
+            {
+                p.hitProcessorPatch
+                (
+                    static_cast<const processorPolyPatch&>(patch), td
+                );
+            }
+            else if (isA<wallPolyPatch>(patch))
+            {
+                p.hitWallPatch
+                (
+                    static_cast<const wallPolyPatch&>(patch), td, faceHitTetIs
+                );
+            }
+            else
+            {
+                p.hitPatch(patch, td);
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+//- Track particle to a given position and returns 1.0 if the
+//  trajectory is completed without hitting a face otherwise
+//  stops at the face and returns the fraction of the trajectory
+//  completed.
+//  on entry 'stepFraction()' should be set to the fraction of the
+//  time-step at which the tracking starts.
+template<class TrackData>
+Foam::scalar Foam::wallBoundedParticle::trackToEdge
+(
+    TrackData& td,
+    const vector& endPosition
+)
+{
+    // Are we on a track face? If not we do a topological walk.
+
+    // Particle:
+    // - cell_              always set
+    // - tetFace_, tetPt_   always set (these identify tet particle is in)
+    // - optionally meshEdgeStart_ or  diagEdge_ set (edge particle is on)
+
+    //checkInside();
+    //checkOnTriangle(position());
+    //if (meshEdgeStart_ != -1 || diagEdge_ != -1)
+    //{
+    //    checkOnEdge();
+    //}
+
+    scalar trackFraction = 0.0;
+
+    if (!td.isWallPatch_[tetFace()])
+    {
+        // Don't track across face. Just walk in cell. Particle is on
+        // mesh edge (as indicated by meshEdgeStart_).
+        const edge meshEdge(currentEdge());
+
+        // If internal face check whether to go to neighbour cell or just
+        // check to the other internal tet on the edge.
+        if (mesh_.isInternalFace(tetFace()))
+        {
+            label nbrCellI =
+            (
+                cellI_ == mesh_.faceOwner()[faceI_]
+              ? mesh_.faceNeighbour()[faceI_]
+              : mesh_.faceOwner()[faceI_]
+            );
+            // Check angle to nbrCell tet. Is it in the direction of the
+            // 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)
+            {
+                // Change into nbrCell. No need to change tetFace, tetPt.
+                //Pout<< "    crossed from cell:" << cellI_
+                //    << " into " << nbrCellI << endl;
+                cellI_ = nbrCellI;
+                patchInteraction(td, trackFraction);
+            }
+            else
+            {
+                // Walk to other face on edge. Changes tetFace, tetPt but not
+                // cell.
+                crossEdgeConnectedFace(meshEdge);
+                patchInteraction(td, trackFraction);
+            }
+        }
+        else
+        {
+            // Walk to other face on edge. This might give loop since
+            // particle should have been removed?
+            crossEdgeConnectedFace(meshEdge);
+            patchInteraction(td, trackFraction);
+        }
+    }
+    else
+    {
+        // We're inside a tet on the wall. Check if the current tet is
+        // the one to cross. If not we cross into the neighbouring triangle.
+
+        if (mesh_.isInternalFace(tetFace()))
+        {
+            FatalErrorIn
+            (
+                "wallBoundedParticle::trackToEdge"
+                "(TrackData&, const vector&)"
+            )   << "Can only track on boundary faces."
+                << " Face:" << tetFace()
+                << " at:" << mesh_.faceCentres()[tetFace()]
+                << abort(FatalError);
+        }
+
+        point projectedEndPosition = endPosition;
+        // 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;
+        }
+
+
+        bool doTrack = false;
+        if (meshEdgeStart_ == -1 && diagEdge_ == -1)
+        {
+            // We're starting and not yet on an edge.
+            doTrack = true;
+        }
+        else
+        {
+            // See if the current triangle has got a point on the
+            // correct side of the edge.
+            doTrack = isTriAlongTrack(projectedEndPosition);
+        }
+
+
+        if (doTrack)
+        {
+            // Track across triangle. Return triangle edge crossed.
+            label triEdgeI = -1;
+            trackFraction = trackFaceTri(projectedEndPosition, triEdgeI);
+
+            if (triEdgeI == -1)
+            {
+                // Reached endpoint
+                //checkInside();
+                diagEdge_ = -1;
+                meshEdgeStart_ = -1;
+                return trackFraction;
+            }
+
+            const tetIndices ti(currentTetIndices());
+
+            // Triangle (faceTriIs) gets constructed from
+            //    f[faceBasePtI_],
+            //    f[facePtAI_],
+            //    f[facePtBI_]
+            //
+            // So edge indices are:
+            // 0 : edge between faceBasePtI_ and facePtAI_
+            // 1 : edge between facePtAI_ and facePtBI_ (is always a real edge)
+            // 2 : edge between facePtBI_ and faceBasePtI_
+
+            const Foam::face& f = mesh_.faces()[ti.face()];
+            const label fp0 = ti.faceBasePt();
+
+            if (triEdgeI == 0)
+            {
+                if (ti.facePtA() == f.fcIndex(fp0))
+                {
+                    //Pout<< "Real edge." << endl;
+                    diagEdge_ = -1;
+                    meshEdgeStart_ = fp0;
+                    //checkOnEdge();
+                    crossEdgeConnectedFace(currentEdge());
+                    patchInteraction(td, trackFraction);
+                }
+                else if (ti.facePtA() == f.rcIndex(fp0))
+                {
+                    //Note: should not happen since boundary face so owner
+                    //Pout<< "Real edge." << endl;
+                    FatalErrorIn("shold not happend") << info()
+                        << abort(FatalError);
+
+                    diagEdge_ = -1;
+                    meshEdgeStart_ = f.rcIndex(fp0);
+                    //checkOnEdge();
+                    crossEdgeConnectedFace(currentEdge());
+                    patchInteraction(td, trackFraction);
+                }
+                else
+                {
+                    // Get index of triangle on other side of edge.
+                    diagEdge_ = ti.facePtA()-fp0;
+                    if (diagEdge_ < 0)
+                    {
+                        diagEdge_ += f.size();
+                    }
+                    meshEdgeStart_ = -1;
+                    //checkOnEdge();
+                    crossDiagonalEdge();
+                }
+            }
+            else if (triEdgeI == 1)
+            {
+                //Pout<< "Real edge." << endl;
+                diagEdge_ = -1;
+                meshEdgeStart_ = ti.facePtA();
+                //checkOnEdge();
+                crossEdgeConnectedFace(currentEdge());
+                patchInteraction(td, trackFraction);
+            }
+            else // if (triEdgeI == 2)
+            {
+                if (ti.facePtB() == f.rcIndex(fp0))
+                {
+                    //Pout<< "Real edge." << endl;
+                    diagEdge_ = -1;
+                    meshEdgeStart_ = ti.facePtB();
+                    //checkOnEdge();
+                    crossEdgeConnectedFace(currentEdge());
+                    patchInteraction(td, trackFraction);
+                }
+                else if (ti.facePtB() == f.fcIndex(fp0))
+                {
+                    //Note: should not happen since boundary face so owner
+                    //Pout<< "Real edge." << endl;
+                    FatalErrorIn("shold not happend") << info()
+                        << abort(FatalError);
+
+                    diagEdge_ = -1;
+                    meshEdgeStart_ = fp0;
+                    //checkOnEdge();
+                    crossEdgeConnectedFace(currentEdge());
+                    patchInteraction(td, trackFraction);
+                }
+                else
+                {
+                    //Pout<< "Triangle edge." << endl;
+                    // Get index of triangle on other side of edge.
+                    diagEdge_ = ti.facePtB()-fp0;
+                    if (diagEdge_ < 0)
+                    {
+                        diagEdge_ += f.size();
+                    }
+                    meshEdgeStart_ = -1;
+                    //checkOnEdge();
+                    crossDiagonalEdge();
+                }
+            }
+        }
+        else
+        {
+            // Current tet is not the right one. Check the neighbour tet.
+
+            if (meshEdgeStart_ != -1)
+            {
+                // Particle is on mesh edge so change into other face on cell
+                crossEdgeConnectedFace(currentEdge());
+                //checkOnEdge();
+                patchInteraction(td, trackFraction);
+            }
+            else
+            {
+                // Particle is on diagonal edge so change into the other
+                // triangle.
+                crossDiagonalEdge();
+                //checkOnEdge();
+            }
+        }
+    }
+
+    //checkInside();
+
+    return trackFraction;
+}
+
+
+template<class TrackData>
+bool Foam::wallBoundedParticle::hitPatch
+(
+    const polyPatch&,
+    TrackData& td,
+    const label patchI,
+    const scalar trackFraction,
+    const tetIndices& tetIs
+)
+{
+    // Disable generic patch interaction
+    return false;
+}
+
+
+template<class TrackData>
+void Foam::wallBoundedParticle::hitWedgePatch
+(
+    const wedgePolyPatch& pp,
+    TrackData& td
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+template<class TrackData>
+void Foam::wallBoundedParticle::hitSymmetryPatch
+(
+    const symmetryPolyPatch& pp,
+    TrackData& td
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+template<class TrackData>
+void Foam::wallBoundedParticle::hitCyclicPatch
+(
+    const cyclicPolyPatch& pp,
+    TrackData& td
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+template<class TrackData>
+void Foam::wallBoundedParticle::hitProcessorPatch
+(
+    const processorPolyPatch& pp,
+    TrackData& td
+)
+{
+    // Switch particle
+    td.switchProcessor = true;
+
+    // Adapt edgeStart_ for other side.
+    // E.g. if edgeStart_ is 1 then the edge is between vertex 1 and 2 so
+    // on the other side between 2 and 3 so edgeStart_ should be
+    // f.size()-edgeStart_-1.
+
+    const Foam::face& f = mesh_.faces()[face()];
+
+    if (meshEdgeStart_ != -1)
+    {
+        meshEdgeStart_ = f.size()-meshEdgeStart_-1;
+    }
+    else
+    {
+        // diagEdge_ is relative to faceBasePt
+        diagEdge_ = f.size()-diagEdge_;
+    }
+}
+
+
+template<class TrackData>
+void Foam::wallBoundedParticle::hitWallPatch
+(
+    const wallPolyPatch& wpp,
+    TrackData& td,
+    const tetIndices&
+)
+{}
+
+
+template<class TrackData>
+void Foam::wallBoundedParticle::hitPatch
+(
+    const polyPatch& wpp,
+    TrackData& td
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+template<class CloudType>
+void Foam::wallBoundedParticle::readFields(CloudType& c)
+{
+    if (!c.size())
+    {
+        return;
+    }
+
+    particle::readFields(c);
+
+    IOField<label> meshEdgeStart
+    (
+        c.fieldIOobject("meshEdgeStart", IOobject::MUST_READ)
+    );
+
+    IOField<label> diagEdge
+    (
+        c.fieldIOobject("diagEdge_", IOobject::MUST_READ)
+    );
+    c.checkFieldIOobject(c, diagEdge);
+
+    label i = 0;
+    forAllIter(typename CloudType, c, iter)
+    {
+        iter().meshEdgeStart_ = meshEdgeStart[i];
+        iter().diagEdge_ = diagEdge[i];
+        i++;
+    }
+}
+
+
+template<class CloudType>
+void Foam::wallBoundedParticle::writeFields(const CloudType& c)
+{
+    particle::writeFields(c);
+
+    label np =  c.size();
+
+    IOField<label> meshEdgeStart
+    (
+        c.fieldIOobject("meshEdgeStart", IOobject::NO_READ),
+        np
+    );
+    IOField<label> diagEdge
+    (
+        c.fieldIOobject("diagEdge", IOobject::NO_READ),
+        np
+    );
+
+    label i = 0;
+    forAllConstIter(typename CloudType, c, iter)
+    {
+        meshEdgeStart[i] = iter().meshEdgeStart_;
+        diagEdge[i] = iter().diagEdge_;
+        i++;
+    }
+
+    meshEdgeStart.write();
+    diagEdge.write();
+}
+
+
+// ************************************************************************* //