From 307a13f08b17bdcd5f4290ea00bada8cd6a37532 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Mon, 6 Feb 2012 16:39:43 +0000
Subject: [PATCH] ENH: wallBoundedStreamLine: near-wall tracking

---
 .../functionObjects/field/Make/files          |    5 +
 .../field/wallBoundedStreamLine/controlDict   |  144 ++
 .../wallBoundedStreamLine.C                   |  887 +++++++++++
 .../wallBoundedStreamLine.H                   |  226 +++
 .../wallBoundedStreamLineFunctionObject.C     |   42 +
 .../wallBoundedStreamLineFunctionObject.H     |   55 +
 .../wallBoundedStreamLineParticle.C           | 1324 +++++++++++++++++
 .../wallBoundedStreamLineParticle.H           |  368 +++++
 .../wallBoundedStreamLineParticleCloud.C      |   64 +
 .../wallBoundedStreamLineParticleCloud.H      |   99 ++
 10 files changed, 3214 insertions(+)
 create mode 100644 src/postProcessing/functionObjects/field/wallBoundedStreamLine/controlDict
 create mode 100644 src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.C
 create mode 100644 src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.H
 create mode 100644 src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineFunctionObject.C
 create mode 100644 src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineFunctionObject.H
 create mode 100644 src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C
 create mode 100644 src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H
 create mode 100644 src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticleCloud.C
 create mode 100644 src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticleCloud.H

diff --git a/src/postProcessing/functionObjects/field/Make/files b/src/postProcessing/functionObjects/field/Make/files
index 1ea75cf3ff..ba19ff80a0 100644
--- a/src/postProcessing/functionObjects/field/Make/files
+++ b/src/postProcessing/functionObjects/field/Make/files
@@ -29,6 +29,11 @@ streamLine/streamLineParticle.C
 streamLine/streamLineParticleCloud.C
 streamLine/streamLineFunctionObject.C
 
+wallBoundedStreamLine/wallBoundedStreamLine.C
+wallBoundedStreamLine/wallBoundedStreamLineFunctionObject.C
+wallBoundedStreamLine/wallBoundedStreamLineParticle.C
+wallBoundedStreamLine/wallBoundedStreamLineParticleCloud.C
+
 surfaceInterpolateFields/surfaceInterpolateFields.C
 surfaceInterpolateFields/surfaceInterpolateFieldsFunctionObject.C
 
diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/controlDict b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/controlDict
new file mode 100644
index 0000000000..4a3f4fe29c
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/controlDict
@@ -0,0 +1,144 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.0.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     icoFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         0.5;
+
+deltaT          0.005;
+
+writeControl    timeStep;
+
+writeInterval   20;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+functions
+{
+    readFields
+    {
+        // Where to load it from (if not already in solver)
+        functionObjectLibs ("libfieldFunctionObjects.so");
+
+        type            readFields;
+        fields          (p U k);
+    }
+
+    near
+    {
+        // Where to load it from
+        functionObjectLibs ("libfieldFunctionObjects.so");
+
+        type nearWallFields;
+
+        // Output every
+        outputControl   outputTime;
+        //outputInterval  1;
+
+        // Fields to be sampled. Per field original name and mapped field to
+        // create.
+        fields
+        (
+            (U UNear)
+        );
+
+        // Patches/groups to sample (regular expressions)
+        patches (motorBike);
+
+        // Distance to sample
+        distance 0.001;
+    }
+    streamLines
+    {
+        // Where to load it from (if not already in solver)
+        functionObjectLibs ("libfieldFunctionObjects.so");
+        type            wallBoundedStreamLine;
+
+        // Output every
+        outputControl   timeStep;   //outputTime;
+        // outputInterval 10;
+
+        setFormat       vtk; //gnuplot; //xmgr; //raw; //jplot;
+
+        // Velocity field to use for tracking.
+        UName UNear;
+
+        // Interpolation method. Default is cellPoint. See sampleDict.
+        //interpolationScheme pointMVC;
+
+        // Tracked forwards (+U) or backwards (-U)
+        trackForward    true;
+
+        interpolationScheme cellPoint;
+
+        // Names of fields to sample. Should contain above velocity field!
+        fields (p U k UNear);
+
+        // Steps particles can travel before being removed
+        lifeTime        100;
+
+        // Cloud name to use
+        cloudName       particleTracks;
+
+        // Seeding method. See the sampleSets in sampleDict.
+        seedSampleSet   patchSeed;    //cloud;//triSurfaceMeshPointSet;
+
+        uniformCoeffs
+        {
+            type        uniform;
+            axis        x;  //distance;
+
+            start       (0.0035 0.0999 0.0001);
+            end         (0.0035 0.0999 0.0099);
+            nPoints     20;
+        }
+
+        cloudCoeffs
+        {
+            type        cloud;
+            axis        x;  //distance;
+            points      ((0.351516548679288 -0.0116085375585099 1.24));
+        }
+        patchSeedCoeffs
+        {
+            type        patchSeed;//patchCloud; //cloud;  //uniform;
+            patches     (motorBike);
+            axis        x;  //distance;
+            maxPoints   20000;
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.C b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.C
new file mode 100644
index 0000000000..4b6f6761f7
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.C
@@ -0,0 +1,887 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "Pstream.H"
+#include "functionObjectList.H"
+#include "wallBoundedStreamLine.H"
+#include "fvMesh.H"
+#include "wallBoundedStreamLineParticleCloud.H"
+#include "ReadFields.H"
+#include "meshSearch.H"
+#include "sampledSet.H"
+#include "globalIndex.H"
+#include "mapDistribute.H"
+#include "interpolationCellPoint.H"
+#include "PatchTools.H"
+#include "meshSearchMeshObject.H"
+#include "faceSet.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(Foam::wallBoundedStreamLine, 0);
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::indirectPrimitivePatch>
+Foam::wallBoundedStreamLine::wallPatch() const
+{
+    const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    label nFaces = 0;
+
+    forAll(patches, patchI)
+    {
+        //if (!polyPatch::constraintType(patches[patchI].type()))
+        if (isA<wallPolyPatch>(patches[patchI]))
+        {
+            nFaces += patches[patchI].size();
+        }
+    }
+
+    labelList addressing(nFaces);
+
+    nFaces = 0;
+
+    forAll(patches, patchI)
+    {
+        //if (!polyPatch::constraintType(patches[patchI].type()))
+        if (isA<wallPolyPatch>(patches[patchI]))
+        {
+            const polyPatch& pp = patches[patchI];
+
+            forAll(pp, i)
+            {
+                addressing[nFaces++] = pp.start()+i;
+            }
+        }
+    }
+
+    return autoPtr<indirectPrimitivePatch>
+    (
+        new indirectPrimitivePatch
+        (
+            IndirectList<face>
+            (
+                mesh.faces(),
+                addressing
+            ),
+            mesh.points()
+        )
+    );
+}
+
+
+Foam::tetIndices Foam::wallBoundedStreamLine::findNearestTet
+(
+    const PackedBoolList& isWallPatch,
+    const point& seedPt,
+    const label cellI
+) const
+{
+    const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
+
+    const cell& cFaces = mesh.cells()[cellI];
+
+    label minFaceI = -1;
+    label minTetPtI = -1;
+    scalar minDistSqr = sqr(GREAT);
+
+    forAll(cFaces, cFaceI)
+    {
+        label faceI = cFaces[cFaceI];
+
+        if (isWallPatch[faceI])
+        {
+            const face& f = mesh.faces()[faceI];
+            const label fp0 = mesh.tetBasePtIs()[faceI];
+            const point& basePoint = mesh.points()[f[fp0]];
+
+            label fp = f.fcIndex(fp0);
+            for (label i = 2; i < f.size(); i++)
+            {
+                const point& thisPoint = mesh.points()[f[fp]];
+                label nextFp = f.fcIndex(fp);
+                const point& nextPoint = mesh.points()[f[nextFp]];
+
+                const triPointRef tri(basePoint, thisPoint, nextPoint);
+
+                scalar d2 = magSqr(tri.centre() - seedPt);
+                if (d2 < minDistSqr)
+                {
+                    minDistSqr = d2;
+                    minFaceI = faceI;
+                    minTetPtI = i-1;
+                }
+                fp = nextFp;
+            }
+        }
+    }
+
+    // Put particle in tet
+    return tetIndices
+    (
+        cellI,
+        minFaceI,
+        minTetPtI,
+        mesh
+    );
+}
+
+
+void Foam::wallBoundedStreamLine::track()
+{
+    const Time& runTime = obr_.time();
+    const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
+
+
+    // Determine the 'wall' patches
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // These are the faces that need to be followed
+
+    autoPtr<indirectPrimitivePatch> boundaryPatch(wallPatch());
+    PackedBoolList isWallPatch(mesh.nFaces());
+    forAll(boundaryPatch().addressing(), i)
+    {
+        isWallPatch[boundaryPatch().addressing()[i]] = 1;
+    }
+
+
+
+    // Find nearest wall particle for the seedPoints
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    IDLList<wallBoundedStreamLineParticle> initialParticles;
+    wallBoundedStreamLineParticleCloud particles
+    (
+        mesh,
+        cloudName_,
+        initialParticles
+    );
+
+    {
+        // Get the seed points
+        // ~~~~~~~~~~~~~~~~~~~
+
+        const sampledSet& seedPoints = sampledSetPtr_();
+
+
+        forAll(seedPoints, i)
+        {
+            const point& seedPt = seedPoints[i];
+            label cellI = seedPoints.cells()[i];
+
+            tetIndices ids(findNearestTet(isWallPatch, seedPt, cellI));
+
+            if (ids.face() != -1 && isWallPatch[ids.face()])
+            {
+                //Pout<< "Seeding particle :" << nl
+                //    << "     seedPt:" << seedPt << nl
+                //    << "     face  :" << ids.face() << nl
+                //    << "     at    :" << mesh.faceCentres()[ids.face()] << nl
+                //    << "     cell  :" << mesh.cellCentres()[ids.cell()] << nl
+                //    << endl;
+
+                particles.addParticle
+                (
+                    new wallBoundedStreamLineParticle
+                    (
+                        mesh,
+                        ids.faceTri(mesh).centre(),
+                        ids.cell(),
+                        ids.face(),     // tetFace
+                        ids.tetPt(),
+                        -1,             // not on a mesh edge
+                        -1,             // not on a diagonal edge
+                        lifeTime_       // lifetime
+                    )
+                );
+            }
+            else
+            {
+                Pout<< type() << " : ignoring seed " << seedPt
+                    << " since not in wall cell." << endl;
+            }
+        }
+    }
+
+    label nSeeds = returnReduce(particles.size(), sumOp<label>());
+
+    Info<< type() << " : seeded " << nSeeds << " particles." << endl;
+
+
+
+    // Read or lookup fields
+    PtrList<volScalarField> vsFlds;
+    PtrList<interpolation<scalar> > vsInterp;
+    PtrList<volVectorField> vvFlds;
+    PtrList<interpolation<vector> > vvInterp;
+
+    label UIndex = -1;
+
+    if (loadFromFiles_)
+    {
+        IOobjectList allObjects(mesh, runTime.timeName());
+
+        IOobjectList objects(2*fields_.size());
+        forAll(fields_, i)
+        {
+            objects.add(*allObjects[fields_[i]]);
+        }
+
+        ReadFields(mesh, objects, vsFlds);
+        vsInterp.setSize(vsFlds.size());
+        forAll(vsFlds, i)
+        {
+            vsInterp.set
+            (
+                i,
+                interpolation<scalar>::New
+                (
+                    interpolationScheme_,
+                    vsFlds[i]
+                )
+            );
+        }
+        ReadFields(mesh, objects, vvFlds);
+        vvInterp.setSize(vvFlds.size());
+        forAll(vvFlds, i)
+        {
+            vvInterp.set
+            (
+                i,
+                interpolation<vector>::New
+                (
+                    interpolationScheme_,
+                    vvFlds[i]
+                )
+            );
+        }
+    }
+    else
+    {
+        label nScalar = 0;
+        label nVector = 0;
+
+        forAll(fields_, i)
+        {
+            if (mesh.foundObject<volScalarField>(fields_[i]))
+            {
+                nScalar++;
+            }
+            else if (mesh.foundObject<volVectorField>(fields_[i]))
+            {
+                nVector++;
+            }
+            else
+            {
+                FatalErrorIn("wallBoundedStreamLine::execute()")
+                    << "Cannot find field " << fields_[i] << endl
+                    << "Valid scalar fields are:"
+                    << mesh.names(volScalarField::typeName) << endl
+                    << "Valid vector fields are:"
+                    << mesh.names(volVectorField::typeName)
+                    << exit(FatalError);
+            }
+        }
+        vsInterp.setSize(nScalar);
+        nScalar = 0;
+        vvInterp.setSize(nVector);
+        nVector = 0;
+
+        forAll(fields_, i)
+        {
+            if (mesh.foundObject<volScalarField>(fields_[i]))
+            {
+                const volScalarField& f = mesh.lookupObject<volScalarField>
+                (
+                    fields_[i]
+                );
+                vsInterp.set
+                (
+                    nScalar++,
+                    interpolation<scalar>::New
+                    (
+                        interpolationScheme_,
+                        f
+                    )
+                );
+            }
+            else if (mesh.foundObject<volVectorField>(fields_[i]))
+            {
+                const volVectorField& f = mesh.lookupObject<volVectorField>
+                (
+                    fields_[i]
+                );
+
+                if (f.name() == UName_)
+                {
+                    UIndex = nVector;
+                }
+
+                vvInterp.set
+                (
+                    nVector++,
+                    interpolation<vector>::New
+                    (
+                        interpolationScheme_,
+                        f
+                    )
+                );
+            }
+        }
+    }
+
+    // Store the names
+    scalarNames_.setSize(vsInterp.size());
+    forAll(vsInterp, i)
+    {
+        scalarNames_[i] = vsInterp[i].psi().name();
+    }
+    vectorNames_.setSize(vvInterp.size());
+    forAll(vvInterp, i)
+    {
+        vectorNames_[i] = vvInterp[i].psi().name();
+    }
+
+    // Check that we know the index of U in the interpolators.
+
+    if (UIndex == -1)
+    {
+        FatalErrorIn("wallBoundedStreamLine::execute()")
+            << "Cannot find field to move particles with : " << UName_
+            << endl
+            << "This field has to be present in the sampled fields "
+            << fields_
+            << " and in the objectRegistry." << endl
+            << exit(FatalError);
+    }
+
+    // Sampled data
+    // ~~~~~~~~~~~~
+
+    // Size to maximum expected sizes.
+    allTracks_.clear();
+    allTracks_.setCapacity(nSeeds);
+    allScalars_.setSize(vsInterp.size());
+    forAll(allScalars_, i)
+    {
+        allScalars_[i].clear();
+        allScalars_[i].setCapacity(nSeeds);
+    }
+    allVectors_.setSize(vvInterp.size());
+    forAll(allVectors_, i)
+    {
+        allVectors_[i].clear();
+        allVectors_[i].setCapacity(nSeeds);
+    }
+
+
+    // additional particle info
+    wallBoundedStreamLineParticle::trackingData td
+    (
+        particles,
+        vsInterp,
+        vvInterp,
+        UIndex,         // index of U in vvInterp
+        trackForward_,  // track in +u direction?
+        isWallPatch,    // which faces are to follow
+
+        allTracks_,
+        allScalars_,
+        allVectors_
+    );
+
+
+    // Set very large dt. Note: cannot use GREAT since 1/GREAT is SMALL
+    // which is a trigger value for the tracking...
+    const scalar trackTime = Foam::sqrt(GREAT);
+
+    // Track
+    particles.move(td, trackTime);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::wallBoundedStreamLine::wallBoundedStreamLine
+(
+    const word& name,
+    const objectRegistry& obr,
+    const dictionary& dict,
+    const bool loadFromFiles
+)
+:
+    dict_(dict),
+    name_(name),
+    obr_(obr),
+    loadFromFiles_(loadFromFiles),
+    active_(true)
+{
+    // Only active if a fvMesh is available
+    if (isA<fvMesh>(obr_))
+    {
+        read(dict_);
+    }
+    else
+    {
+        active_ = false;
+        WarningIn
+        (
+            "wallBoundedStreamLine::wallBoundedStreamLine\n"
+            "(\n"
+                "const word&,\n"
+                "const objectRegistry&,\n"
+                "const dictionary&,\n"
+                "const bool\n"
+            ")"
+        )   << "No fvMesh available, deactivating."
+            << nl << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::wallBoundedStreamLine::~wallBoundedStreamLine()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::wallBoundedStreamLine::read(const dictionary& dict)
+{
+    if (active_)
+    {
+        //dict_ = dict;
+        dict.lookup("fields") >> fields_;
+        if (dict.found("UName"))
+        {
+            dict.lookup("UName") >> UName_;
+        }
+        else
+        {
+            UName_ = "U";
+            if (dict.found("U"))
+            {
+                IOWarningIn
+                (
+                    "wallBoundedStreamLine::read(const dictionary&)",
+                    dict
+                )   << "Using deprecated entry \"U\"."
+                    << " Please use \"UName\" instead."
+                    << endl;
+                dict.lookup("U") >> UName_;
+            }
+        }
+
+        if (findIndex(fields_, UName_) == -1)
+        {
+            FatalIOErrorIn
+            (
+                "wallBoundedStreamLine::read(const dictionary&)",
+                dict
+            )   << "Velocity field for tracking " << UName_
+                << " should be present in the list of fields " << fields_
+                << exit(FatalIOError);
+        }
+
+
+        dict.lookup("trackForward") >> trackForward_;
+        dict.lookup("lifeTime") >> lifeTime_;
+        if (lifeTime_ < 1)
+        {
+            FatalErrorIn(":wallBoundedStreamLine::read(const dictionary&)")
+                << "Illegal value " << lifeTime_ << " for lifeTime"
+                << exit(FatalError);
+        }
+
+
+        interpolationScheme_ = dict.lookupOrDefault
+        (
+            "interpolationScheme",
+            interpolationCellPoint<scalar>::typeName
+        );
+
+        //Info<< typeName << " using interpolation " << interpolationScheme_
+        //    << endl;
+
+        cloudName_ = dict.lookupOrDefault<word>
+        (
+            "cloudName",
+            "wallBoundedStreamLine"
+        );
+        dict.lookup("seedSampleSet") >> seedSet_;
+
+        const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
+
+        const dictionary& coeffsDict = dict.subDict(seedSet_ + "Coeffs");
+
+        sampledSetPtr_ = sampledSet::New
+        (
+            seedSet_,
+            mesh,
+            meshSearchMeshObject::New(mesh),
+            coeffsDict
+        );
+        coeffsDict.lookup("axis") >> sampledSetAxis_;
+
+        scalarFormatterPtr_ = writer<scalar>::New(dict.lookup("setFormat"));
+        vectorFormatterPtr_ = writer<vector>::New(dict.lookup("setFormat"));
+
+
+        // Make sure that the mesh is trackable
+        if (debug)
+        {
+            // 1. positive volume decomposition tets
+            faceSet faces(mesh, "lowQualityTetFaces", mesh.nFaces()/100+1);
+            if
+            (
+                polyMeshTetDecomposition::checkFaceTets
+                (
+                    mesh,
+                    polyMeshTetDecomposition::minTetQuality,
+                    true,
+                    &faces
+                )
+            )
+            {
+                label nFaces = returnReduce(faces.size(), sumOp<label>());
+
+                FatalErrorIn("wallBoundedStreamLine::track()")
+                    << "Found " << nFaces
+                    <<" faces with low quality or negative volume "
+                    << "decomposition tets. Writing to faceSet " << faces.name()
+                    << exit(FatalError);
+            }
+
+            // 2. all edges on a cell having two faces
+            EdgeMap<label> numFacesPerEdge;
+            forAll(mesh.cells(), cellI)
+            {
+                const cell& cFaces = mesh.cells()[cellI];
+
+                numFacesPerEdge.clear();
+
+                forAll(cFaces, cFaceI)
+                {
+                    label faceI = cFaces[cFaceI];
+                    const face& f = mesh.faces()[faceI];
+                    forAll(f, fp)
+                    {
+                        const edge e(f[fp], f.nextLabel(fp));
+                        EdgeMap<label>::iterator eFnd =
+                            numFacesPerEdge.find(e);
+                        if (eFnd != numFacesPerEdge.end())
+                        {
+                            eFnd()++;
+                        }
+                        else
+                        {
+                            numFacesPerEdge.insert(e, 1);
+                        }
+                    }
+                }
+
+                forAllConstIter(EdgeMap<label>, numFacesPerEdge, iter)
+                {
+                    if (iter() != 2)
+                    {
+                        FatalErrorIn
+                        (
+                            "wallBoundedStreamLine::read(..)"
+                        )   << "problem cell:" << cellI
+                            << abort(FatalError);
+                    }
+                }
+            }
+        }
+    }
+}
+
+
+void Foam::wallBoundedStreamLine::execute()
+{}
+
+
+void Foam::wallBoundedStreamLine::end()
+{}
+
+
+void Foam::wallBoundedStreamLine::write()
+{
+    if (active_)
+    {
+        const Time& runTime = obr_.time();
+        const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
+
+
+        // Do all injection and tracking
+        track();
+
+
+        if (Pstream::parRun())
+        {
+            // Append slave tracks to master ones
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            globalIndex globalTrackIDs(allTracks_.size());
+
+            // Construct a distribution map to pull all to the master.
+            labelListList sendMap(Pstream::nProcs());
+            labelListList recvMap(Pstream::nProcs());
+
+            if (Pstream::master())
+            {
+                // Master: receive all. My own first, then consecutive
+                // processors.
+                label trackI = 0;
+
+                forAll(recvMap, procI)
+                {
+                    labelList& fromProc = recvMap[procI];
+                    fromProc.setSize(globalTrackIDs.localSize(procI));
+                    forAll(fromProc, i)
+                    {
+                        fromProc[i] = trackI++;
+                    }
+                }
+            }
+
+            labelList& toMaster = sendMap[0];
+            toMaster.setSize(globalTrackIDs.localSize());
+            forAll(toMaster, i)
+            {
+                toMaster[i] = i;
+            }
+
+            const mapDistribute distMap
+            (
+                globalTrackIDs.size(),
+                sendMap.xfer(),
+                recvMap.xfer()
+            );
+
+
+            // Distribute the track positions. Note: use scheduled comms
+            // to prevent buffering.
+            mapDistribute::distribute
+            (
+                Pstream::scheduled,
+                distMap.schedule(),
+                distMap.constructSize(),
+                distMap.subMap(),
+                distMap.constructMap(),
+                allTracks_
+            );
+
+            // Distribute the scalars
+            forAll(allScalars_, scalarI)
+            {
+                mapDistribute::distribute
+                (
+                    Pstream::scheduled,
+                    distMap.schedule(),
+                    distMap.constructSize(),
+                    distMap.subMap(),
+                    distMap.constructMap(),
+                    allScalars_[scalarI]
+                );
+            }
+            // Distribute the vectors
+            forAll(allVectors_, vectorI)
+            {
+                mapDistribute::distribute
+                (
+                    Pstream::scheduled,
+                    distMap.schedule(),
+                    distMap.constructSize(),
+                    distMap.subMap(),
+                    distMap.constructMap(),
+                    allVectors_[vectorI]
+                );
+            }
+        }
+
+
+        label n = 0;
+        forAll(allTracks_, trackI)
+        {
+            n += allTracks_[trackI].size();
+        }
+
+        Info<< "Tracks:" << allTracks_.size()
+            << "  total samples:" << n << endl;
+
+
+        // Massage into form suitable for writers
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        if (Pstream::master() && allTracks_.size())
+        {
+            // Make output directory
+
+            fileName vtkPath
+            (
+                Pstream::parRun()
+              ? runTime.path()/".."/"sets"/name()
+              : runTime.path()/"sets"/name()
+            );
+            if (mesh.name() != fvMesh::defaultRegion)
+            {
+                vtkPath = vtkPath/mesh.name();
+            }
+            vtkPath = vtkPath/mesh.time().timeName();
+
+            mkDir(vtkPath);
+
+            // Convert track positions
+
+            PtrList<coordSet> tracks(allTracks_.size());
+            forAll(allTracks_, trackI)
+            {
+                tracks.set
+                (
+                    trackI,
+                    new coordSet
+                    (
+                        "track" + Foam::name(trackI),
+                        sampledSetAxis_                 //"xyz"
+                    )
+                );
+                tracks[trackI].transfer(allTracks_[trackI]);
+            }
+
+            // Convert scalar values
+
+            if (allScalars_.size() > 0)
+            {
+                List<List<scalarField> > scalarValues(allScalars_.size());
+
+                forAll(allScalars_, scalarI)
+                {
+                    DynamicList<scalarList>& allTrackVals =
+                        allScalars_[scalarI];
+                    scalarValues[scalarI].setSize(allTrackVals.size());
+
+                    forAll(allTrackVals, trackI)
+                    {
+                        scalarList& trackVals = allTrackVals[trackI];
+                        scalarValues[scalarI][trackI].transfer(trackVals);
+                    }
+                }
+
+                fileName vtkFile
+                (
+                    vtkPath
+                  / scalarFormatterPtr_().getFileName
+                    (
+                        tracks[0],
+                        scalarNames_
+                    )
+                );
+
+                Info<< "Writing data to " << vtkFile.path() << endl;
+
+                scalarFormatterPtr_().write
+                (
+                    true,           // writeTracks
+                    tracks,
+                    scalarNames_,
+                    scalarValues,
+                    OFstream(vtkFile)()
+                );
+            }
+
+            // Convert vector values
+
+            if (allVectors_.size() > 0)
+            {
+                List<List<vectorField> > vectorValues(allVectors_.size());
+
+                forAll(allVectors_, vectorI)
+                {
+                    DynamicList<vectorList>& allTrackVals =
+                        allVectors_[vectorI];
+                    vectorValues[vectorI].setSize(allTrackVals.size());
+
+                    forAll(allTrackVals, trackI)
+                    {
+                        vectorList& trackVals = allTrackVals[trackI];
+                        vectorValues[vectorI][trackI].transfer(trackVals);
+                    }
+                }
+
+                fileName vtkFile
+                (
+                    vtkPath
+                  / vectorFormatterPtr_().getFileName
+                    (
+                        tracks[0],
+                        vectorNames_
+                    )
+                );
+
+                //Info<< "Writing vector data to " << vtkFile << endl;
+
+                vectorFormatterPtr_().write
+                (
+                    true,           // writeTracks
+                    tracks,
+                    vectorNames_,
+                    vectorValues,
+                    OFstream(vtkFile)()
+                );
+            }
+        }
+    }
+}
+
+
+void Foam::wallBoundedStreamLine::updateMesh(const mapPolyMesh&)
+{
+    read(dict_);
+}
+
+
+void Foam::wallBoundedStreamLine::movePoints(const pointField&)
+{
+    // Moving mesh affects the search tree
+    read(dict_);
+}
+
+
+//void Foam::wallBoundedStreamLine::readUpdate
+//(const polyMesh::readUpdateState state)
+//{
+//    if (state != UNCHANGED)
+//    {
+//        read(dict_);
+//    }
+//}
+
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.H b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.H
new file mode 100644
index 0000000000..6a39d25b37
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.H
@@ -0,0 +1,226 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::wallBoundedStreamLine
+
+Description
+    Generation of streamlines. Samples along track of passive particle.
+
+SourceFiles
+    wallBoundedStreamLine.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef wallBoundedStreamLine_H
+#define wallBoundedStreamLine_H
+
+#include "volFieldsFwd.H"
+#include "pointFieldFwd.H"
+#include "Switch.H"
+#include "DynamicList.H"
+#include "scalarList.H"
+#include "vectorList.H"
+#include "polyMesh.H"
+#include "writer.H"
+#include "indirectPrimitivePatch.H"
+#include "tetIndices.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class objectRegistry;
+class dictionary;
+class mapPolyMesh;
+class meshSearch;
+class sampledSet;
+
+/*---------------------------------------------------------------------------*\
+                         Class wallBoundedStreamLine Declaration
+\*---------------------------------------------------------------------------*/
+
+class wallBoundedStreamLine
+{
+    // Private data
+
+        //- Input dictionary
+        dictionary dict_;
+
+        //- Name of this set of field averages.
+        word name_;
+
+        //- Database this class is registered to
+        const objectRegistry& obr_;
+
+        //- Load fields from files (not from objectRegistry)
+        bool loadFromFiles_;
+
+        //- On/off switch
+        bool active_;
+
+
+        //- List of fields to sample
+        wordList fields_;
+
+        //- Field to transport particle with
+        word UName_;
+
+        //- Interpolation scheme to use
+        word interpolationScheme_;
+
+        //- Whether to use +u or -u
+        bool trackForward_;
+
+        //- Maximum lifetime (= number of cells) of particle
+        label lifeTime_;
+
+        //- Optional specified name of particles
+        word cloudName_;
+
+        //- Type of seed
+        word seedSet_;
+
+        //- Names of scalar fields
+        wordList scalarNames_;
+
+        //- Names of vector fields
+        wordList vectorNames_;
+
+
+        // Demand driven
+
+            //- Mesh searching enigne
+            autoPtr<meshSearch> meshSearchPtr_;
+
+            //- Seed set engine
+            autoPtr<sampledSet> sampledSetPtr_;
+
+            //- Axis of the sampled points to output
+            word sampledSetAxis_;
+
+            //- File output writer
+            autoPtr<writer<scalar> > scalarFormatterPtr_;
+
+            autoPtr<writer<vector> > vectorFormatterPtr_;
+
+
+        // Generated data
+
+            //- All tracks. Per particle the points it passed through
+            DynamicList<List<point> > allTracks_;
+
+            //- Per scalarField, per particle, the sampled value.
+            List<DynamicList<scalarList> > allScalars_;
+
+            //- Per scalarField, per particle, the sampled value.
+            List<DynamicList<vectorList> > allVectors_;
+
+
+        //- Construct patch out of all wall patch faces
+        autoPtr<indirectPrimitivePatch> wallPatch() const;
+
+        //- Find wall tet on cell
+        tetIndices findNearestTet
+        (
+            const PackedBoolList& isWallPatch,
+            const point& seedPt,
+            const label cellI
+        ) const;
+
+        //- Do all seeding and tracking
+        void track();
+
+        //- Disallow default bitwise copy construct
+        wallBoundedStreamLine(const wallBoundedStreamLine&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const wallBoundedStreamLine&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("wallBoundedStreamLine");
+
+
+    // Constructors
+
+        //- Construct for given objectRegistry and dictionary.
+        //  Allow the possibility to load fields from files
+        wallBoundedStreamLine
+        (
+            const word& name,
+            const objectRegistry&,
+            const dictionary&,
+            const bool loadFromFiles = false
+        );
+
+
+    //- Destructor
+    virtual ~wallBoundedStreamLine();
+
+
+    // Member Functions
+
+        //- Return name of the set of field averages
+        virtual const word& name() const
+        {
+            return name_;
+        }
+
+        //- Read the field average data
+        virtual void read(const dictionary&);
+
+        //- Execute the averaging
+        virtual void execute();
+
+        //- Execute the averaging at the final time-loop, currently does nothing
+        virtual void end();
+
+        //- Calculate the field average data and write
+        virtual void write();
+
+        //- Update for changes of mesh
+        virtual void updateMesh(const mapPolyMesh&);
+
+        //- Update for mesh point-motion
+        virtual void movePoints(const pointField&);
+
+        ////- Update for changes of mesh due to readUpdate
+        //virtual void readUpdate(const polyMesh::readUpdateState state);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineFunctionObject.C b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineFunctionObject.C
new file mode 100644
index 0000000000..eaa1cc7eae
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineFunctionObject.C
@@ -0,0 +1,42 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "wallBoundedStreamLineFunctionObject.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineNamedTemplateTypeNameAndDebug(wallBoundedStreamLineFunctionObject, 0);
+
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        wallBoundedStreamLineFunctionObject,
+        dictionary
+    );
+}
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineFunctionObject.H b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineFunctionObject.H
new file mode 100644
index 0000000000..0ce95cc592
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineFunctionObject.H
@@ -0,0 +1,55 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Typedef
+    Foam::wallBoundedStreamLineFunctionObject
+
+Description
+    FunctionObject wrapper around wallBoundedStreamLines
+    to allow them to be created via
+    the functions entry within controlDict.
+
+SourceFiles
+    wallBoundedStreamLineFunctionObject.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef wallBoundedStreamLineFunctionObject_H
+#define wallBoundedStreamLineFunctionObject_H
+
+#include "wallBoundedStreamLine.H"
+#include "OutputFilterFunctionObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef OutputFilterFunctionObject<wallBoundedStreamLine>
+        wallBoundedStreamLineFunctionObject;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C
new file mode 100644
index 0000000000..cd7605a56e
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C
@@ -0,0 +1,1324 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "wallBoundedStreamLineParticle.H"
+#include "vectorFieldIOField.H"
+#include "polyMesh.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+//    defineParticleTypeNameAndDebug(wallBoundedStreamLineParticle, 0);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+//// Check position is inside tet
+//void Foam::wallBoundedStreamLineParticle::checkInside() const
+//{
+//    const tetIndices ti(currentTetIndices());
+//    const tetPointRef tpr(ti.tet(mesh_));
+//    if (!tpr.inside(position()))
+//    {
+//        FatalErrorIn("wallBoundedStreamLineParticle::checkInside(..)")
+//            << "Particle:" //<< static_cast<const particle&>(*this)
+//            << info()
+//            << "is not inside " << tpr
+//            << abort(FatalError);
+//    }
+//}
+//
+//
+//void Foam::wallBoundedStreamLineParticle::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
+//        (
+//            "wallBoundedStreamLineParticle::checkOnEdge()"
+//        )   << "Problem :"
+//            << " particle:" //<< static_cast<const particle&>(*this)
+//            << info()
+//            << "edge:" << e
+//            << " at:" << ln
+//            << " distance:" << ph.distance()
+//            << abort(FatalError);
+//    }
+//}
+//
+//
+//void Foam::wallBoundedStreamLineParticle::checkOnTriangle(const point& p)
+//const
+//{
+//    const triFace tri(currentTetIndices().faceTriIs(mesh_));
+//    pointHit ph = tri.nearestPoint(p, mesh_.points());
+//    if (ph.distance() > 1e-9)
+//    {
+//        FatalErrorIn
+//        (
+//            "wallBoundedStreamLineParticle::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::wallBoundedStreamLineParticle::currentEdge() const
+{
+    if ((meshEdgeStart_ != -1) == (diagEdge_ != -1))
+    {
+        FatalErrorIn("wallBoundedStreamLineParticle::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::wallBoundedStreamLineParticle::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
+            (
+                "wallBoundedStreamLineParticle::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
+        (
+            "wallBoundedStreamLineParticle::crossEdgeConnectedFace"
+            "(const edge&)"
+        )   << "Problem" << abort(FatalError);
+    }
+
+
+    // Check that edge (as indicated by meshEdgeStart_) is indeed one that
+    // contains the position.
+    //checkOnEdge();
+}
+
+
+void Foam::wallBoundedStreamLineParticle::crossDiagonalEdge()
+{
+    if (diagEdge_ == -1)
+    {
+        FatalErrorIn("wallBoundedStreamLineParticle::crossDiagonalEdge()")
+            << "Particle:" //<< static_cast<const particle&>(*this)
+            << info()
+            << "not on a diagonal edge" << abort(FatalError);
+    }
+    if (meshEdgeStart_ != -1)
+    {
+        FatalErrorIn("wallBoundedStreamLineParticle::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("wallBoundedStreamLineParticle::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::wallBoundedStreamLineParticle::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("wallBoundedStreamLineParticle::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("wallBoundedStreamLineParticle::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::wallBoundedStreamLineParticle::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
+            (
+                "wallBoundedStreamLineParticle::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
+    (
+        "wallBoundedStreamLineParticle::isTriAlongTrack"
+        "(const point&)"
+    )   << "Problem" << abort(FatalError);
+
+    return false;
+}
+
+
+void Foam::wallBoundedStreamLineParticle::patchInteraction
+(
+    trackingData& td,
+    const scalar trackFraction
+)
+{
+    typedef trackingData::cloudType cloudType;
+    typedef 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);
+            }
+        }
+    }
+}
+
+//- 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.
+Foam::scalar Foam::wallBoundedStreamLineParticle::trackToEdge
+(
+    trackingData& 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
+            (
+                "wallBoundedStreamLineParticle::trackToEdge"
+                "(trackingData&, 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();
+                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;
+}
+
+
+Foam::vector Foam::wallBoundedStreamLineParticle::interpolateFields
+(
+    const trackingData& td,
+    const point& position,
+    const label cellI,
+    const label faceI
+)
+{
+    if (cellI == -1)
+    {
+        FatalErrorIn("wallBoundedStreamLineParticle::interpolateFields(..)")
+            << "Cell:" << cellI << abort(FatalError);
+    }
+
+    const tetIndices ti = currentTetIndices();
+
+
+    sampledScalars_.setSize(td.vsInterp_.size());
+    forAll(td.vsInterp_, scalarI)
+    {
+        sampledScalars_[scalarI].append
+        (
+            td.vsInterp_[scalarI].interpolate
+            (
+                position,
+                ti,     //cellI,
+                faceI
+            )
+        );
+    }
+
+    sampledVectors_.setSize(td.vvInterp_.size());
+    forAll(td.vvInterp_, vectorI)
+    {
+        sampledVectors_[vectorI].append
+        (
+            td.vvInterp_[vectorI].interpolate
+            (
+                position,
+                ti,     //cellI,
+                faceI
+            )
+        );
+    }
+
+    const DynamicList<vector>& U = sampledVectors_[td.UIndex_];
+
+    return U.last();
+}
+
+
+Foam::vector Foam::wallBoundedStreamLineParticle::sample
+(
+    trackingData& td
+)
+{
+    sampledPositions_.append(position());
+    vector U = interpolateFields(td, position(), cell(), tetFace());
+
+    if (!td.trackForward_)
+    {
+        U = -U;
+    }
+
+    scalar magU = mag(U);
+
+    if (magU < SMALL)
+    {
+        // Stagnant particle. Might as well stop
+        lifeTime_ = 0;
+    }
+
+    U /= magU;
+
+    return U;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::wallBoundedStreamLineParticle::wallBoundedStreamLineParticle
+(
+    const polyMesh& mesh,
+    const vector& position,
+    const label cellI,
+    const label tetFaceI,
+    const label tetPtI,
+    const label meshEdgeStart,
+    const label diagEdge,
+    const label lifeTime
+)
+:
+    particle(mesh, position, cellI, tetFaceI, tetPtI),
+    meshEdgeStart_(meshEdgeStart),
+    diagEdge_(diagEdge),
+    lifeTime_(lifeTime)
+{
+    //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::wallBoundedStreamLineParticle::wallBoundedStreamLineParticle
+(
+    const polyMesh& mesh,
+    Istream& is,
+    bool readFields
+)
+:
+    particle(mesh, is, readFields)
+{
+    if (readFields)
+    {
+        List<scalarList> sampledScalars;
+        List<vectorList> sampledVectors;
+
+        is  >> meshEdgeStart_ >> diagEdge_ >> lifeTime_
+            >> sampledPositions_ >> sampledScalars >> sampledVectors;
+
+        sampledScalars_.setSize(sampledScalars.size());
+        forAll(sampledScalars, i)
+        {
+            sampledScalars_[i].transfer(sampledScalars[i]);
+        }
+        sampledVectors_.setSize(sampledVectors.size());
+        forAll(sampledVectors, i)
+        {
+            sampledVectors_[i].transfer(sampledVectors[i]);
+        }
+    }
+
+    // Check state of Istream
+    is.check
+    (
+        "wallBoundedStreamLineParticle::wallBoundedStreamLineParticle"
+        "(const Cloud<wallBoundedStreamLineParticle>&, Istream&, bool)"
+    );
+}
+
+
+Foam::wallBoundedStreamLineParticle::wallBoundedStreamLineParticle
+(
+    const wallBoundedStreamLineParticle& p
+)
+:
+    particle(p),
+    meshEdgeStart_(p.meshEdgeStart_),
+    diagEdge_(p.diagEdge_),
+    lifeTime_(p.lifeTime_),
+    sampledPositions_(p.sampledPositions_),
+    sampledScalars_(p.sampledScalars_),
+    sampledVectors_(p.sampledVectors_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::wallBoundedStreamLineParticle::move
+(
+    trackingData& td,
+    const scalar trackTime
+)
+{
+    wallBoundedStreamLineParticle& p = static_cast
+    <
+        wallBoundedStreamLineParticle&
+    >(*this);
+
+
+    // Check position is inside tet
+    //checkInside();
+
+    td.switchProcessor = false;
+    td.keepParticle = true;
+
+    scalar tEnd = (1.0 - stepFraction())*trackTime;
+    scalar maxDt = mesh_.bounds().mag();
+
+    while
+    (
+        td.keepParticle
+    && !td.switchProcessor
+    && lifeTime_ > 0
+    )
+    {
+        // set the lagrangian time-step
+        scalar dt = maxDt;
+
+        --lifeTime_;
+
+        // Update sampled velocity and fields if position changed
+        vector U = sample(td);
+
+        // !user parameter!
+        if (dt < SMALL)
+        {
+            // Force removal
+            lifeTime_ = 0;
+            break;
+        }
+
+
+        scalar fraction = trackToEdge(td, position() + dt*U);
+        dt *= fraction;
+
+        tEnd -= dt;
+        stepFraction() = 1.0 - tEnd/trackTime;
+
+
+        if (tEnd <= ROOTVSMALL)
+        {
+            // Force removal
+            lifeTime_ = 0;
+        }
+
+        if
+        (
+            !td.keepParticle
+        ||  td.switchProcessor
+        ||  lifeTime_ == 0
+        )
+        {
+            break;
+        }
+    }
+
+
+    if (!td.keepParticle || lifeTime_ == 0)
+    {
+        if (lifeTime_ == 0)
+        {
+            if (debug)
+            {
+                Pout<< "wallBoundedStreamLineParticle :"
+                    << " Removing stagnant particle:"
+                    << p.position()
+                    << " sampled positions:" << sampledPositions_.size()
+                    << endl;
+            }
+            td.keepParticle = false;
+        }
+        else
+        {
+            // Normal exit. Store last position and fields
+            sample(td);
+
+            if (debug)
+            {
+                Pout<< "wallBoundedStreamLineParticle : Removing particle:"
+                    << p.position()
+                    << " sampled positions:" << sampledPositions_.size()
+                    << endl;
+            }
+        }
+
+        // Transfer particle data into trackingData.
+        {
+            //td.allPositions_.append(sampledPositions_);
+            td.allPositions_.append(vectorList());
+            vectorList& top = td.allPositions_.last();
+            top.transfer(sampledPositions_);
+        }
+
+        forAll(sampledScalars_, i)
+        {
+            //td.allScalars_[i].append(sampledScalars_[i]);
+            td.allScalars_[i].append(scalarList());
+            scalarList& top = td.allScalars_[i].last();
+            top.transfer(sampledScalars_[i]);
+        }
+        forAll(sampledVectors_, i)
+        {
+            //td.allVectors_[i].append(sampledVectors_[i]);
+            td.allVectors_[i].append(vectorList());
+            vectorList& top = td.allVectors_[i].last();
+            top.transfer(sampledVectors_[i]);
+        }
+    }
+
+    return td.keepParticle;
+}
+
+
+bool Foam::wallBoundedStreamLineParticle::hitPatch
+(
+    const polyPatch&,
+    trackingData& td,
+    const label patchI,
+    const scalar trackFraction,
+    const tetIndices& tetIs
+)
+{
+    // Disable generic patch interaction
+    return false;
+}
+
+
+void Foam::wallBoundedStreamLineParticle::hitWedgePatch
+(
+    const wedgePolyPatch& pp,
+    trackingData& td
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+void Foam::wallBoundedStreamLineParticle::hitSymmetryPatch
+(
+    const symmetryPolyPatch& pp,
+    trackingData& td
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+void Foam::wallBoundedStreamLineParticle::hitCyclicPatch
+(
+    const cyclicPolyPatch& pp,
+    trackingData& td
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+void Foam::wallBoundedStreamLineParticle::hitProcessorPatch
+(
+    const processorPolyPatch& pp,
+    trackingData& 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_;
+    }
+}
+
+
+void Foam::wallBoundedStreamLineParticle::hitWallPatch
+(
+    const wallPolyPatch& wpp,
+    trackingData& td,
+    const tetIndices&
+)
+{}
+
+
+void Foam::wallBoundedStreamLineParticle::hitPatch
+(
+    const polyPatch& wpp,
+    trackingData& td
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+void Foam::wallBoundedStreamLineParticle::readFields
+(
+    Cloud<wallBoundedStreamLineParticle>& 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);
+
+    IOField<label> lifeTime
+    (
+        c.fieldIOobject("lifeTime", IOobject::MUST_READ)
+    );
+    c.checkFieldIOobject(c, lifeTime);
+
+    vectorFieldIOField sampledPositions
+    (
+        c.fieldIOobject("sampledPositions", IOobject::MUST_READ)
+    );
+    c.checkFieldIOobject(c, sampledPositions);
+
+    label i = 0;
+    forAllIter(Cloud<wallBoundedStreamLineParticle>, c, iter)
+    {
+        iter().meshEdgeStart_ = meshEdgeStart[i];
+        iter().diagEdge_ = diagEdge[i];
+        iter().lifeTime_ = lifeTime[i];
+        iter().sampledPositions_.transfer(sampledPositions[i]);
+        i++;
+    }
+}
+
+
+void Foam::wallBoundedStreamLineParticle::writeFields
+(
+    const Cloud<wallBoundedStreamLineParticle>& 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
+    );
+    IOField<label> lifeTime
+    (
+        c.fieldIOobject("lifeTime", IOobject::NO_READ),
+        np
+    );
+    vectorFieldIOField sampledPositions
+    (
+        c.fieldIOobject("sampledPositions", IOobject::NO_READ),
+        np
+    );
+
+    label i = 0;
+    forAllConstIter(Cloud<wallBoundedStreamLineParticle>, c, iter)
+    {
+        meshEdgeStart[i] = iter().meshEdgeStart_;
+        diagEdge[i] = iter().diagEdge_;
+        lifeTime[i] = iter().lifeTime_;
+        sampledPositions[i] = iter().sampledPositions_;
+        i++;
+    }
+
+    meshEdgeStart.write();
+    diagEdge.write();
+    lifeTime.write();
+    sampledPositions.write();
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const wallBoundedStreamLineParticle& p
+)
+{
+    os  << static_cast<const particle&>(p)
+        << token::SPACE << p.meshEdgeStart_
+        << token::SPACE << p.diagEdge_
+        << token::SPACE << p.lifeTime_
+        << token::SPACE << p.sampledPositions_
+        << token::SPACE << p.sampledScalars_
+        << token::SPACE << p.sampledVectors_;
+
+    // Check state of Ostream
+    os.check
+    (
+        "Ostream& operator<<(Ostream&, const wallBoundedStreamLineParticle&)"
+    );
+
+    return os;
+}
+
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const InfoProxy<wallBoundedStreamLineParticle>& ip
+)
+{
+    const wallBoundedStreamLineParticle& 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/wallBoundedStreamLineParticle.H b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H
new file mode 100644
index 0000000000..c0e19dd090
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H
@@ -0,0 +1,368 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::wallBoundedStreamLineParticle
+
+Description
+    Particle class that samples fields as it passes through. Used in streamline
+    calculation.
+
+SourceFiles
+    wallBoundedStreamLineParticle.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef wallBoundedStreamLineParticle_H
+#define wallBoundedStreamLineParticle_H
+
+#include "particle.H"
+#include "autoPtr.H"
+#include "interpolation.H"
+#include "vectorList.H"
+#include "InfoProxy.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class wallBoundedStreamLineParticleCloud;
+
+/*---------------------------------------------------------------------------*\
+                     Class wallBoundedStreamLineParticle Declaration
+\*---------------------------------------------------------------------------*/
+
+class wallBoundedStreamLineParticle
+:
+    public particle
+{
+
+public:
+
+    //- Class used to pass tracking data to the trackToFace function
+    class trackingData
+    :
+        public particle::TrackingData<Cloud<wallBoundedStreamLineParticle> >
+    {
+
+    public:
+
+
+        const PtrList<interpolation<scalar> >& vsInterp_;
+        const PtrList<interpolation<vector> >& vvInterp_;
+        const label UIndex_;
+        const bool trackForward_;
+
+        const PackedBoolList& isWallPatch_;
+
+        DynamicList<vectorList>& allPositions_;
+        List<DynamicList<scalarList> >& allScalars_;
+        List<DynamicList<vectorList> >& allVectors_;
+
+
+        // Constructors
+
+            trackingData
+            (
+                Cloud<wallBoundedStreamLineParticle>& cloud,
+                const PtrList<interpolation<scalar> >& vsInterp,
+                const PtrList<interpolation<vector> >& vvInterp,
+                const label UIndex,
+                const bool trackForward,
+                const PackedBoolList& isWallPatch,
+
+                DynamicList<List<point> >& allPositions,
+                List<DynamicList<scalarList> >& allScalars,
+                List<DynamicList<vectorList> >& allVectors
+            )
+            :
+                particle::TrackingData<Cloud<wallBoundedStreamLineParticle> >
+                (
+                    cloud
+                ),
+                vsInterp_(vsInterp),
+                vvInterp_(vvInterp),
+                UIndex_(UIndex),
+                trackForward_(trackForward),
+                isWallPatch_(isWallPatch),
+
+                allPositions_(allPositions),
+                allScalars_(allScalars),
+                allVectors_(allVectors)
+            {}
+    };
+
+
+private:
+
+    // Private 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_;
+
+        //- Lifetime of particle. Particle dies when reaches 0.
+        label lifeTime_;
+
+        //- sampled positions
+        DynamicList<point> sampledPositions_;
+
+        //- sampled scalars
+        List<DynamicList<scalar> > sampledScalars_;
+
+        //- sampled vectors
+        List<DynamicList<vector> > sampledVectors_;
+
+
+    // Private 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);
+
+        //- Do all patch interaction
+        void patchInteraction(trackingData& td, const scalar trackFraction);
+
+        //- Is current triangle in the track direction
+        bool isTriAlongTrack(const point& endPosition) const;
+
+        //- Equivalent of trackToFace
+        scalar trackToEdge
+        (
+            trackingData& td,
+            const vector& endPosition
+        );
+
+        //- Interpolate all quantities; return interpolated velocity.
+        vector interpolateFields
+        (
+            const trackingData&,
+            const point&,
+            const label cellI,
+            const label faceI
+        );
+
+        //- Interpolate all quantities, return updated velocity.
+        vector sample(trackingData& td);
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        wallBoundedStreamLineParticle
+        (
+            const polyMesh& c,
+            const vector& position,
+            const label cellI,
+            const label tetFaceI,
+            const label tetPtI,
+            const label meshEdgeStart,
+            const label diagEdge,
+            const label lifeTime
+        );
+
+        //- Construct from Istream
+        wallBoundedStreamLineParticle
+        (
+            const polyMesh& c,
+            Istream& is,
+            bool readFields = true
+        );
+
+        //- Construct copy
+        wallBoundedStreamLineParticle(const wallBoundedStreamLineParticle& p);
+
+        //- Construct and return a clone
+        autoPtr<particle> clone() const
+        {
+            return autoPtr<particle>(new wallBoundedStreamLineParticle(*this));
+        }
+
+        //- Factory class to read-construct particles used for
+        //  parallel transfer
+        class iNew
+        {
+            const polyMesh& mesh_;
+
+        public:
+
+            iNew(const polyMesh& mesh)
+            :
+                mesh_(mesh)
+            {}
+
+            autoPtr<wallBoundedStreamLineParticle> operator()
+            (
+                Istream& is
+            ) const
+            {
+                return autoPtr<wallBoundedStreamLineParticle>
+                (
+                    new wallBoundedStreamLineParticle(mesh_, is, true)
+                );
+            }
+        };
+
+
+    // Member Functions
+
+        // Tracking
+
+            //- Track all particles to their end point
+            bool move(trackingData&, const scalar trackTime);
+
+
+            //- Overridable function to handle the particle hitting a patch
+            //  Executed before other patch-hitting functions
+            bool hitPatch
+            (
+                const polyPatch&,
+                trackingData& td,
+                const label patchI,
+                const scalar trackFraction,
+                const tetIndices& tetIs
+            );
+
+            //- Overridable function to handle the particle hitting a wedge
+            void hitWedgePatch
+            (
+                const wedgePolyPatch&,
+                trackingData& td
+            );
+
+            //- Overridable function to handle the particle hitting a
+            //  symmetryPlane
+            void hitSymmetryPatch
+            (
+                const symmetryPolyPatch&,
+                trackingData& td
+            );
+
+            //- Overridable function to handle the particle hitting a cyclic
+            void hitCyclicPatch
+            (
+                const cyclicPolyPatch&,
+                trackingData& td
+            );
+
+            //- Overridable function to handle the particle hitting a
+            //- processorPatch
+            void hitProcessorPatch
+            (
+                const processorPolyPatch&,
+                trackingData& td
+            );
+
+            //- Overridable function to handle the particle hitting a wallPatch
+            void hitWallPatch
+            (
+                const wallPolyPatch&,
+                trackingData& td,
+                const tetIndices&
+            );
+
+            //- Overridable function to handle the particle hitting a polyPatch
+            void hitPatch
+            (
+                const polyPatch&,
+                trackingData& td
+            );
+
+
+        // Info
+
+            //- Return info proxy.
+            //  Used to print particle information to a stream
+            InfoProxy<wallBoundedStreamLineParticle> info() const
+            {
+                return *this;
+            }
+
+
+        // I-O
+
+            //- Read
+            static void readFields(Cloud<wallBoundedStreamLineParticle>&);
+
+            //- Write
+            static void writeFields
+            (
+                const Cloud<wallBoundedStreamLineParticle>&
+            );
+
+
+    // Ostream Operator
+
+        friend Ostream& operator<<
+        (
+            Ostream&,
+            const wallBoundedStreamLineParticle&
+        );
+
+        friend Ostream& operator<<
+        (
+            Ostream&,
+            const InfoProxy<wallBoundedStreamLineParticle>&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticleCloud.C b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticleCloud.C
new file mode 100644
index 0000000000..dbe7ce8276
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticleCloud.C
@@ -0,0 +1,64 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "wallBoundedStreamLineParticleCloud.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTemplateTypeNameAndDebug(Cloud<wallBoundedStreamLineParticle>, 0);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::wallBoundedStreamLineParticleCloud::wallBoundedStreamLineParticleCloud
+(
+    const polyMesh& mesh,
+    const word& cloudName,
+    bool readFields
+)
+:
+    Cloud<wallBoundedStreamLineParticle>(mesh, cloudName, false)
+{
+    if (readFields)
+    {
+        wallBoundedStreamLineParticle::readFields(*this);
+    }
+}
+
+
+Foam::wallBoundedStreamLineParticleCloud::wallBoundedStreamLineParticleCloud
+(
+    const polyMesh& mesh,
+    const word& cloudName,
+    const IDLList<wallBoundedStreamLineParticle>& particles
+)
+:
+    Cloud<wallBoundedStreamLineParticle>(mesh, cloudName, particles)
+{}
+
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticleCloud.H b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticleCloud.H
new file mode 100644
index 0000000000..b91e9b5d8f
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticleCloud.H
@@ -0,0 +1,99 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::wallBoundedStreamLineParticleCloud
+
+Description
+    A Cloud of streamLine particles
+
+SourceFiles
+    streamLineCloud.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef streamLineCloud_H
+#define streamLineCloud_H
+
+#include "Cloud.H"
+#include "wallBoundedStreamLineParticle.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class streamLineCloud Declaration
+\*---------------------------------------------------------------------------*/
+
+class wallBoundedStreamLineParticleCloud
+:
+    public Cloud<wallBoundedStreamLineParticle>
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        wallBoundedStreamLineParticleCloud
+        (
+            const wallBoundedStreamLineParticleCloud&
+        );
+
+        //- Disallow default bitwise assignment
+        void operator=(const wallBoundedStreamLineParticleCloud&);
+
+
+public:
+
+    //- Type of parcel the cloud was instantiated for
+    typedef wallBoundedStreamLineParticle parcelType;
+
+    // Constructors
+
+        //- Construct given mesh
+        wallBoundedStreamLineParticleCloud
+        (
+            const polyMesh&,
+            const word& cloudName = "defaultCloud",
+            bool readFields = true
+        );
+
+        //- Construct from mesh, cloud name, and a list of particles
+        wallBoundedStreamLineParticleCloud
+        (
+            const polyMesh& mesh,
+            const word& cloudName,
+            const IDLList<wallBoundedStreamLineParticle>& particles
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //