diff --git a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C
index aec4c64304..ded2af139b 100644
--- a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C
+++ b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C
@@ -77,6 +77,9 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
+        // Store the particle positions
+        kinematicCloud.storeGlobalPositions();
+
         mesh.update();
 
         // Calculate absolute flux from the mapped surface velocity
diff --git a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C
index be0905796b..00c66823ea 100644
--- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C
+++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,6 +68,8 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
+        kinematicCloud.storeGlobalPositions();
+
         mesh.update();
 
         U.correctBoundaryConditions();
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/sprayDyMFoam.C b/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/sprayDyMFoam.C
index 16d2d524c5..40acccae9e 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/sprayDyMFoam.C
+++ b/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/sprayDyMFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -91,6 +91,9 @@ int main(int argc, char *argv[])
             // Store momentum to set rhoUf for introduced faces.
             volVectorField rhoU("rhoU", rho*U);
 
+            // Store the particle positions
+            parcels.storeGlobalPositions();
+
             // Do any mesh changes
             mesh.update();
 
diff --git a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/files b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/files
new file mode 100644
index 0000000000..430b1df24f
--- /dev/null
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/files
@@ -0,0 +1,3 @@
+uncoupledKinematicParcelDyMFoam.C
+
+EXE = $(FOAM_APPBIN)/uncoupledKinematicParcelDyMFoam
diff --git a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options
new file mode 100644
index 0000000000..2d385f9853
--- /dev/null
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options
@@ -0,0 +1,36 @@
+EXE_INC = \
+    -I.. \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude
+
+EXE_LIBS = \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -llagrangianTurbulence \
+    -lcompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lradiationModels \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -lregionModels \
+    -lsurfaceFilmModels \
+    -ldynamicMesh \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh
diff --git a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/uncoupledKinematicParcelDyMFoam.C b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/uncoupledKinematicParcelDyMFoam.C
new file mode 100644
index 0000000000..18f8e63869
--- /dev/null
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/uncoupledKinematicParcelDyMFoam.C
@@ -0,0 +1,90 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 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/>.
+
+Application
+    uncoupledKinematicParcelDyMFoam
+
+Description
+    Transient solver for the passive transport of a particle cloud.
+
+    Uses a pre- calculated velocity field to evolve the cloud.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "psiThermo.H"
+#include "turbulentFluidThermoModel.H"
+#include "basicKinematicCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addOption
+    (
+        "cloudName",
+        "name",
+        "specify alternative cloud name. default is 'kinematicCloud'"
+    );
+
+    #define NO_CONTROL
+    #include "postProcess.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+    #include "createFields.H"
+    #include "compressibleCourantNo.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        kinematicCloud.storeGlobalPositions();
+
+        mesh.update();
+
+        U.correctBoundaryConditions();
+
+        Info<< "Evolving " << kinematicCloud.name() << endl;
+        kinematicCloud.evolve();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.C b/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.C
index efb7dd2a23..98457dc2b4 100644
--- a/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.C
+++ b/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -47,13 +47,12 @@ Foam::lagrangianFieldDecomposer::lagrangianFieldDecomposer
 {
     label pi = 0;
 
-    // faceProcAddressing not required currently
-    // labelList decodedProcFaceAddressing(faceProcAddressing.size());
+    labelList decodedProcFaceAddressing(faceProcAddressing.size());
 
-    // forAll(faceProcAddressing, i)
-    // {
-    //     decodedProcFaceAddressing[i] = mag(faceProcAddressing[i]) - 1;
-    // }
+    forAll(faceProcAddressing, i)
+    {
+        decodedProcFaceAddressing[i] = mag(faceProcAddressing[i]) - 1;
+    }
 
     forAll(cellProcAddressing, procCelli)
     {
@@ -68,27 +67,28 @@ Foam::lagrangianFieldDecomposer::lagrangianFieldDecomposer
                 const indexedParticle& ppi = *iter();
                 particleIndices_[pi++] = ppi.index();
 
-                // label mappedTetFace = findIndex
-                // (
-                //     decodedProcFaceAddressing,
-                //     ppi.tetFace()
-                // );
+                label mappedTetFace = findIndex
+                (
+                    decodedProcFaceAddressing,
+                    ppi.tetFace()
+                );
 
-                // if (mappedTetFace == -1)
-                // {
-                //     FatalErrorInFunction
-                //         << "Face lookup failure." << nl
-                //         << abort(FatalError);
-                // }
+                if (mappedTetFace == -1)
+                {
+                    FatalErrorInFunction
+                        << "Face lookup failure." << nl
+                        << abort(FatalError);
+                }
 
                 positions_.append
                 (
                     new passiveParticle
                     (
                         procMesh,
-                        ppi.position(),
+                        ppi.coordinates(),
                         procCelli,
-                        false
+                        mappedTetFace,
+                        ppi.procTetPt(procMesh, procCelli, mappedTetFace)
                     )
                 );
             }
diff --git a/applications/utilities/preProcessing/mapFields/mapLagrangian.C b/applications/utilities/preProcessing/mapFields/mapLagrangian.C
index ad80ba3d07..45cd8e4bb5 100644
--- a/applications/utilities/preProcessing/mapFields/mapLagrangian.C
+++ b/applications/utilities/preProcessing/mapFields/mapLagrangian.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -98,7 +98,6 @@ void mapLagrangian(const meshToMesh0& meshToMesh0Interp)
 
     const fvMesh& meshSource = meshToMesh0Interp.fromMesh();
     const fvMesh& meshTarget = meshToMesh0Interp.toMesh();
-    const pointField& targetCc = meshTarget.cellCentres();
 
 
     fileNameList cloudDirs
@@ -182,15 +181,17 @@ void mapLagrangian(const meshToMesh0& meshToMesh0Interp)
                             new passiveParticle
                             (
                                 meshTarget,
-                                targetCc[targetCells[i]],
-                                targetCells[i]
+                                barycentric(1, 0, 0, 0),
+                                targetCells[i],
+                                meshTarget.cells()[targetCells[i]][0],
+                                1
                             )
                         );
                         passiveParticle& newP = newPtr();
 
-                        label facei = newP.track(iter().position(), td);
+                        newP.track(iter().position() - newP.position(), 0);
 
-                        if (facei < 0 && newP.cell() >= 0)
+                        if (!newP.onFace())
                         {
                             // Hit position.
                             foundCell = true;
@@ -233,11 +234,16 @@ void mapLagrangian(const meshToMesh0& meshToMesh0Interp)
                         {
                             unmappedSource.erase(sourceParticleI);
                             addParticles.append(sourceParticleI);
-                            iter().cell() = targetCell;
                             targetParcels.addParticle
                             (
-                                sourceParcels.remove(&iter())
+                                new passiveParticle
+                                (
+                                    meshTarget,
+                                    iter().position(),
+                                    targetCell
+                                )
                             );
+                            sourceParcels.remove(&iter());
                         }
                     }
                     sourceParticleI++;
diff --git a/applications/utilities/preProcessing/mapFieldsPar/mapLagrangian.C b/applications/utilities/preProcessing/mapFieldsPar/mapLagrangian.C
index 06f7edc910..2f41a6eee6 100644
--- a/applications/utilities/preProcessing/mapFieldsPar/mapLagrangian.C
+++ b/applications/utilities/preProcessing/mapFieldsPar/mapLagrangian.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -89,8 +89,6 @@ void mapLagrangian(const meshToMesh& interp)
     const polyMesh& meshTarget = interp.tgtRegion();
     const labelListList& sourceToTarget = interp.srcToTgtCellAddr();
 
-    const pointField& targetCc = meshTarget.cellCentres();
-
     fileNameList cloudDirs
     (
         readDir
@@ -173,15 +171,17 @@ void mapLagrangian(const meshToMesh& interp)
                             new passiveParticle
                             (
                                 meshTarget,
-                                targetCc[targetCells[i]],
-                                targetCells[i]
+                                barycentric(1, 0, 0, 0),
+                                targetCells[i],
+                                meshTarget.cells()[targetCells[i]][0],
+                                1
                             )
                         );
                         passiveParticle& newP = newPtr();
 
-                        label facei = newP.track(iter().position(), td);
+                        newP.track(iter().position() - newP.position(), 0);
 
-                        if (facei < 0 && newP.cell() >= 0)
+                        if (!newP.onFace())
                         {
                             // Hit position.
                             foundCell = true;
@@ -224,11 +224,16 @@ void mapLagrangian(const meshToMesh& interp)
                         {
                             unmappedSource.erase(sourceParticleI);
                             addParticles.append(sourceParticleI);
-                            iter().cell() = targetCell;
                             targetParcels.addParticle
                             (
-                                sourceParcels.remove(&iter())
+                                new passiveParticle
+                                (
+                                    meshTarget,
+                                    iter().position(),
+                                    targetCell
+                                )
                             );
+                            sourceParcels.remove(&iter());
                         }
                     }
                     sourceParticleI++;
diff --git a/etc/caseDicts/meshQualityDict b/etc/caseDicts/meshQualityDict
index e69207563a..fe5a1f321a 100644
--- a/etc/caseDicts/meshQualityDict
+++ b/etc/caseDicts/meshQualityDict
@@ -33,8 +33,7 @@ minVol 1e-13;
 
 //- Minimum quality of the tet formed by the face-centre
 //  and variable base point minimum decomposition triangles and
-//  the cell centre. This has to be a positive number for tracking
-//  to work. Set to very negative number (e.g. -1E30) to
+//  the cell centre. Set to very negative number (e.g. -1E30) to
 //  disable.
 //     <0 = inside out tet,
 //      0 = flat tet
diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index d9f2104c18..f79161fa78 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -133,6 +133,8 @@ $(spatialVectorAlgebra)/CompactSpatialTensor/compactSpatialTensor/compactSpatial
 primitives/polynomialEqns/cubicEqn/cubicEqn.C
 primitives/polynomialEqns/quadraticEqn/quadraticEqn.C
 
+primitives/Barycentric/barycentric/barycentric.C
+
 containers/HashTables/HashTable/HashTableCore.C
 containers/HashTables/StaticHashTable/StaticHashTableCore.C
 containers/Lists/SortableList/ParSortableListName.C
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
index 8437bcde5c..52b8c4c54c 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -135,6 +135,29 @@ void Foam::polyMesh::calcDirections() const
 }
 
 
+Foam::autoPtr<Foam::labelIOList> Foam::polyMesh::readTetBasePtIs() const
+{
+    IOobject io
+    (
+        "tetBasePtIs",
+        instance(),
+        meshSubDir,
+        *this,
+        IOobject::READ_IF_PRESENT,
+        IOobject::NO_WRITE
+    );
+
+    if (io.headerOk())
+    {
+        return autoPtr<labelIOList>(new labelIOList(io));
+    }
+    else
+    {
+        return autoPtr<labelIOList>(nullptr);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::polyMesh::polyMesh(const IOobject& io)
@@ -207,7 +230,7 @@ Foam::polyMesh::polyMesh(const IOobject& io)
     comm_(UPstream::worldComm),
     geometricD_(Zero),
     solutionD_(Zero),
-    tetBasePtIsPtr_(nullptr),
+    tetBasePtIsPtr_(readTetBasePtIs()),
     cellTreePtr_(nullptr),
     pointZones_
     (
@@ -401,7 +424,7 @@ Foam::polyMesh::polyMesh
     comm_(UPstream::worldComm),
     geometricD_(Zero),
     solutionD_(Zero),
-    tetBasePtIsPtr_(nullptr),
+    tetBasePtIsPtr_(readTetBasePtIs()),
     cellTreePtr_(nullptr),
     pointZones_
     (
@@ -552,7 +575,7 @@ Foam::polyMesh::polyMesh
     comm_(UPstream::worldComm),
     geometricD_(Zero),
     solutionD_(Zero),
-    tetBasePtIsPtr_(nullptr),
+    tetBasePtIsPtr_(readTetBasePtIs()),
     cellTreePtr_(nullptr),
     pointZones_
     (
@@ -815,7 +838,7 @@ Foam::label Foam::polyMesh::nSolutionD() const
 }
 
 
-const Foam::labelList& Foam::polyMesh::tetBasePtIs() const
+const Foam::labelIOList& Foam::polyMesh::tetBasePtIs() const
 {
     if (tetBasePtIsPtr_.empty())
     {
@@ -828,8 +851,17 @@ const Foam::labelList& Foam::polyMesh::tetBasePtIs() const
 
         tetBasePtIsPtr_.reset
         (
-            new labelList
+            new labelIOList
             (
+                IOobject
+                (
+                    "tetBasePtIs",
+                    instance(),
+                    meshSubDir,
+                    *this,
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::NO_WRITE
+                ),
                 polyMeshTetDecomposition::findFaceBasePts(*this)
             )
         );
@@ -1105,6 +1137,13 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
     points_.instance() = time().timeName();
     points_.eventNo() = getEvent();
 
+    if (tetBasePtIsPtr_.valid())
+    {
+        tetBasePtIsPtr_().writeOpt() = IOobject::AUTO_WRITE;
+        tetBasePtIsPtr_().instance() = time().timeName();
+        tetBasePtIsPtr_().eventNo() = getEvent();
+    }
+
     tmp<scalarField> sweptVols = primitiveMesh::movePoints
     (
         points_,
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.H b/src/OpenFOAM/meshes/polyMesh/polyMesh.H
index 1f61f73828..00ff5aa2f9 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -150,7 +150,7 @@ private:
             mutable Vector<label> solutionD_;
 
             //- Base point for face decomposition into tets
-            mutable autoPtr<labelList> tetBasePtIsPtr_;
+            mutable autoPtr<labelIOList> tetBasePtIsPtr_;
 
             //- Search tree to allow spatial cell searching
             mutable autoPtr<indexedOctree<treeDataCell>> cellTreePtr_;
@@ -208,6 +208,9 @@ private:
         //  polyhedral information
         void calcCellShapes() const;
 
+        //- Read and return the tetBasePtIs
+        autoPtr<labelIOList> readTetBasePtIs() const;
+
 
         // Helper functions for constructor from cell shapes
 
@@ -449,7 +452,7 @@ public:
             label nSolutionD() const;
 
             //- Return the tetBasePtIs
-            const labelList& tetBasePtIs() const;
+            const labelIOList& tetBasePtIs() const;
 
             //- Return the cell search tree
             const indexedOctree<treeDataCell>& cellTree() const;
@@ -605,8 +608,8 @@ public:
             //- Clear primitive data (points, faces and cells)
             void clearPrimitives();
 
-            //- Clear geometry not used for CFD (cellTree, tetBasePtIs)
-            void clearAdditionalGeom();
+            //- Clear tet base points
+            void clearTetBasePtIs();
 
             //- Clear cell tree data
             void clearCellTree();
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C b/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C
index fa1bdc8bdd..73c4537ef2 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,22 +69,6 @@ void Foam::polyMesh::clearGeom()
     geometricD_ = Zero;
     solutionD_ = Zero;
 
-    // Remove the stored tet base points
-    tetBasePtIsPtr_.clear();
-    // Remove the cell tree
-    cellTreePtr_.clear();
-}
-
-
-void Foam::polyMesh::clearAdditionalGeom()
-{
-    if (debug)
-    {
-        InfoInFunction << "Clearing additional geometric data" << endl;
-    }
-
-    // Remove the stored tet base points
-    tetBasePtIsPtr_.clear();
     // Remove the cell tree
     cellTreePtr_.clear();
 }
@@ -170,6 +154,17 @@ void Foam::polyMesh::clearOut()
 }
 
 
+void Foam::polyMesh::clearTetBasePtIs()
+{
+    if (debug)
+    {
+        InfoInFunction << "Clearing tet base points" << endl;
+    }
+
+    tetBasePtIsPtr_.clear();
+}
+
+
 void Foam::polyMesh::clearCellTree()
 {
     if (debug)
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshIO.C b/src/OpenFOAM/meshes/polyMesh/polyMeshIO.C
index 7803f094be..544ef0ecfc 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshIO.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshIO.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -59,6 +59,12 @@ void Foam::polyMesh::setInstance(const fileName& inst)
 
     cellZones_.writeOpt() = IOobject::AUTO_WRITE;
     cellZones_.instance() = inst;
+
+    if (tetBasePtIsPtr_.valid())
+    {
+        tetBasePtIsPtr_->writeOpt() = IOobject::AUTO_WRITE;
+        tetBasePtIsPtr_->instance() = inst;
+    }
 }
 
 
@@ -387,6 +393,9 @@ Foam::polyMesh::readUpdateState Foam::polyMesh::readUpdate()
             cellZones_.set(czI, newCellZones[czI].clone(cellZones_));
         }
 
+        // Re-read tet base points
+        tetBasePtIsPtr_ = readTetBasePtIs();
+
 
         if (boundaryChanged)
         {
@@ -439,6 +448,13 @@ Foam::polyMesh::readUpdateState Foam::polyMesh::readUpdate()
         points_.transfer(newPoints);
         points_.instance() = pointsInst;
 
+        // Re-read tet base points
+        autoPtr<labelIOList> newTetBasePtIsPtr = readTetBasePtIs();
+        if (newTetBasePtIsPtr.valid())
+        {
+            tetBasePtIsPtr_ = newTetBasePtIsPtr;
+        }
+
         // Calculate the geometry for the patches (transformation tensors etc.)
         boundary_.calcGeometry();
 
diff --git a/src/OpenFOAM/primitives/Barycentric/Barycentric.H b/src/OpenFOAM/primitives/Barycentric/Barycentric.H
new file mode 100644
index 0000000000..e430401cf8
--- /dev/null
+++ b/src/OpenFOAM/primitives/Barycentric/Barycentric.H
@@ -0,0 +1,117 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016-2017 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::Barycentric
+
+Description
+    Templated 3D Barycentric derived from VectorSpace. Has 4 components, one of
+    which is redundant.
+
+SourceFiles
+    BarycentricI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Barycentric_H
+#define Barycentric_H
+
+#include "VectorSpace.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class Barycentric Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Cmpt>
+class Barycentric
+:
+    public VectorSpace<Barycentric<Cmpt>, Cmpt, 4>
+{
+public:
+
+    //- Equivalent type of labels used for valid component indexing
+    typedef Barycentric<label> labelType;
+
+
+    // Member constants
+
+        //- Rank of Barycentric is 1
+        static const direction rank = 1;
+
+
+    //- Component labeling enumeration
+    enum components { A, B, C, D };
+
+
+    // Constructors
+
+        //- Construct null
+        inline Barycentric();
+
+        //- Construct initialized to zero
+        inline Barycentric(const Foam::zero);
+
+        //- Construct given four components
+        inline Barycentric
+        (
+            const Cmpt& va,
+            const Cmpt& vb,
+            const Cmpt& vc,
+            const Cmpt& vd
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            inline const Cmpt& a() const;
+            inline const Cmpt& b() const;
+            inline const Cmpt& c() const;
+            inline const Cmpt& d() const;
+
+            inline Cmpt& a();
+            inline Cmpt& b();
+            inline Cmpt& c();
+            inline Cmpt& d();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "BarycentricI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/Barycentric/BarycentricI.H b/src/OpenFOAM/primitives/Barycentric/BarycentricI.H
new file mode 100644
index 0000000000..ad649cdac5
--- /dev/null
+++ b/src/OpenFOAM/primitives/Barycentric/BarycentricI.H
@@ -0,0 +1,132 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016-2017 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
+    ANB 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Cmpt>
+inline Foam::Barycentric<Cmpt>::Barycentric()
+{}
+
+
+template<class Cmpt>
+inline Foam::Barycentric<Cmpt>::Barycentric(const Foam::zero)
+:
+    Barycentric::vsType(Zero)
+{}
+
+
+template<class Cmpt>
+inline Foam::Barycentric<Cmpt>::Barycentric
+(
+    const Cmpt& va,
+    const Cmpt& vb,
+    const Cmpt& vc,
+    const Cmpt& vd
+)
+{
+    this->v_[A] = va;
+    this->v_[B] = vb;
+    this->v_[C] = vc;
+    this->v_[D] = vd;
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Cmpt>
+inline const Cmpt& Foam::Barycentric<Cmpt>::a() const
+{
+    return this->v_[A];
+}
+
+
+template<class Cmpt>
+inline const Cmpt& Foam::Barycentric<Cmpt>::b() const
+{
+    return this->v_[B];
+}
+
+
+template<class Cmpt>
+inline const Cmpt& Foam::Barycentric<Cmpt>::c() const
+{
+    return this->v_[C];
+}
+
+
+template<class Cmpt>
+inline const Cmpt& Foam::Barycentric<Cmpt>::d() const
+{
+    return this->v_[D];
+}
+
+
+template<class Cmpt>
+inline Cmpt& Foam::Barycentric<Cmpt>::a()
+{
+    return this->v_[A];
+}
+
+
+template<class Cmpt>
+inline Cmpt& Foam::Barycentric<Cmpt>::b()
+{
+    return this->v_[B];
+}
+
+
+template<class Cmpt>
+inline Cmpt& Foam::Barycentric<Cmpt>::c()
+{
+    return this->v_[C];
+}
+
+
+template<class Cmpt>
+inline Cmpt& Foam::Barycentric<Cmpt>::d()
+{
+    return this->v_[D];
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
+
+template<class Cmpt>
+inline Cmpt operator&(const Barycentric<Cmpt>& b1, const Barycentric<Cmpt>& b2)
+{
+    return b1.a()*b2.a() + b1.b()*b2.b() + b1.c()*b2.c() + b1.d()*b2.d();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/Barycentric/BarycentricTensor.H b/src/OpenFOAM/primitives/Barycentric/BarycentricTensor.H
new file mode 100644
index 0000000000..2d83fafd8f
--- /dev/null
+++ b/src/OpenFOAM/primitives/Barycentric/BarycentricTensor.H
@@ -0,0 +1,139 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 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::BarycentricTensor
+
+Description
+    Templated 4x3 tensor derived from VectorSpace. Has 12 components. Can
+    represent a barycentric transformation as a matrix-barycentric inner-
+    product. Can alternatively represent an inverse barycentric transformation
+    as a vector-matrix inner-product.
+
+SourceFiles
+    BarycentricTensorI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef BarycentricTensor_H
+#define BarycentricTensor_H
+
+#include "Barycentric.H"
+#include "Tensor.H"
+#include "Vector.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class BarycentricTensor Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Cmpt>
+class BarycentricTensor
+:
+    public MatrixSpace<BarycentricTensor<Cmpt>, Cmpt, 4, 3>
+{
+public:
+
+    //- Equivalent type of labels used for valid component indexing
+    typedef Tensor<label> labelType;
+
+
+    // Member constants
+
+        //- Rank of BarycentricTensor is 2
+        static const direction rank = 2;
+
+
+    //- Component labeling enumeration
+    enum components { XA, XB, XC, XD, YA, YB, YC, YD, ZA, ZB, ZC, ZD };
+
+
+    // Constructors
+
+        //- Construct null
+        BarycentricTensor();
+
+        //- Construct initialised to zero
+        BarycentricTensor(const Foam::zero);
+
+        //- Construct given three barycentric components (rows)
+        BarycentricTensor
+        (
+            const Barycentric<Cmpt>& x,
+            const Barycentric<Cmpt>& y,
+            const Barycentric<Cmpt>& z
+        );
+
+        //- Construct given four vector components (columns)
+        BarycentricTensor
+        (
+            const Vector<Cmpt>& a,
+            const Vector<Cmpt>& b,
+            const Vector<Cmpt>& c,
+            const Vector<Cmpt>& d
+        );
+
+
+    // Member Functions
+
+        // Row-barycentric access
+
+            inline Barycentric<Cmpt> x() const;
+            inline Barycentric<Cmpt> y() const;
+            inline Barycentric<Cmpt> z() const;
+
+        // Column-vector access
+
+            inline Vector<Cmpt> a() const;
+            inline Vector<Cmpt> b() const;
+            inline Vector<Cmpt> c() const;
+            inline Vector<Cmpt> d() const;
+};
+
+
+template<class Cmpt>
+class typeOfTranspose<Cmpt, BarycentricTensor<Cmpt>>
+{
+public:
+
+    typedef void type;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "BarycentricTensorI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/Barycentric/BarycentricTensorI.H b/src/OpenFOAM/primitives/Barycentric/BarycentricTensorI.H
new file mode 100644
index 0000000000..84e461377f
--- /dev/null
+++ b/src/OpenFOAM/primitives/Barycentric/BarycentricTensorI.H
@@ -0,0 +1,196 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Cmpt>
+inline Foam::BarycentricTensor<Cmpt>::BarycentricTensor()
+{}
+
+
+template<class Cmpt>
+inline Foam::BarycentricTensor<Cmpt>::BarycentricTensor(const Foam::zero)
+:
+    BarycentricTensor::msType(Zero)
+{}
+
+
+template<class Cmpt>
+inline Foam::BarycentricTensor<Cmpt>::BarycentricTensor
+(
+    const Barycentric<Cmpt>& x,
+    const Barycentric<Cmpt>& y,
+    const Barycentric<Cmpt>& z
+)
+{
+    this->v_[XA] = x.a();
+    this->v_[XB] = x.b();
+    this->v_[XC] = x.c();
+    this->v_[XD] = x.d();
+
+    this->v_[YA] = y.a();
+    this->v_[YB] = y.b();
+    this->v_[YC] = y.c();
+    this->v_[YD] = y.d();
+
+    this->v_[ZA] = z.a();
+    this->v_[ZB] = z.b();
+    this->v_[ZC] = z.c();
+    this->v_[ZD] = z.d();
+}
+
+
+template<class Cmpt>
+inline Foam::BarycentricTensor<Cmpt>::BarycentricTensor
+(
+    const Vector<Cmpt>& a,
+    const Vector<Cmpt>& b,
+    const Vector<Cmpt>& c,
+    const Vector<Cmpt>& d
+)
+{
+    this->v_[XA] = a.x();
+    this->v_[XB] = b.x();
+    this->v_[XC] = c.x();
+    this->v_[XD] = d.x();
+
+    this->v_[YA] = a.y();
+    this->v_[YB] = b.y();
+    this->v_[YC] = c.y();
+    this->v_[YD] = d.y();
+
+    this->v_[ZA] = a.z();
+    this->v_[ZB] = b.z();
+    this->v_[ZC] = c.z();
+    this->v_[ZD] = d.z();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Cmpt>
+inline Foam::Barycentric<Cmpt> Foam::BarycentricTensor<Cmpt>::x() const
+{
+    return
+        Barycentric<Cmpt>
+        (
+            this->v_[XA],
+            this->v_[XB],
+            this->v_[XC],
+            this->v_[XD]
+        );
+}
+
+
+template<class Cmpt>
+inline Foam::Barycentric<Cmpt> Foam::BarycentricTensor<Cmpt>::y() const
+{
+    return
+        Barycentric<Cmpt>
+        (
+            this->v_[YA],
+            this->v_[YB],
+            this->v_[YC],
+            this->v_[YD]
+        );
+}
+
+
+template<class Cmpt>
+inline Foam::Barycentric<Cmpt> Foam::BarycentricTensor<Cmpt>::z() const
+{
+    return
+        Barycentric<Cmpt>
+        (
+            this->v_[ZA],
+            this->v_[ZB],
+            this->v_[ZC],
+            this->v_[ZD]
+        );
+}
+
+
+template<class Cmpt>
+inline Foam::Vector<Cmpt> Foam::BarycentricTensor<Cmpt>::a() const
+{
+    return Vector<Cmpt>(this->v_[XA], this->v_[YA], this->v_[ZA]);
+}
+
+
+template<class Cmpt>
+inline Foam::Vector<Cmpt> Foam::BarycentricTensor<Cmpt>::b() const
+{
+    return Vector<Cmpt>(this->v_[XB], this->v_[YB], this->v_[ZB]);
+}
+
+
+template<class Cmpt>
+inline Foam::Vector<Cmpt> Foam::BarycentricTensor<Cmpt>::c() const
+{
+    return Vector<Cmpt>(this->v_[XC], this->v_[YC], this->v_[ZC]);
+}
+
+
+template<class Cmpt>
+inline Foam::Vector<Cmpt> Foam::BarycentricTensor<Cmpt>::d() const
+{
+    return Vector<Cmpt>(this->v_[XD], this->v_[YD], this->v_[ZD]);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
+
+template<class Cmpt>
+inline Vector<Cmpt> operator&
+(
+    const BarycentricTensor<Cmpt>& T,
+    const Barycentric<Cmpt>& b
+)
+{
+    return Vector<Cmpt>(T.x() & b, T.y() & b, T.z() & b);
+}
+
+
+template<class Cmpt>
+inline Barycentric<Cmpt> operator&
+(
+    const Vector<Cmpt>& v,
+    const BarycentricTensor<Cmpt>& T
+)
+{
+    return Barycentric<Cmpt>(v & T.a(), v & T.b(), v & T.c(), v & T.d());
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/Barycentric/barycentric/barycentric.C b/src/OpenFOAM/primitives/Barycentric/barycentric/barycentric.C
new file mode 100644
index 0000000000..f42aa1a9fe
--- /dev/null
+++ b/src/OpenFOAM/primitives/Barycentric/barycentric/barycentric.C
@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 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 "barycentric.H"
+#include "Random.H"
+#include "cachedRandom.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::barycentric barycentric01
+(
+    Foam::scalar s,
+    Foam::scalar t,
+    Foam::scalar u
+)
+{
+    // Transform the random point in the unit cube to a random point in the
+    // unit tet by means of a series of reflections. See
+    // <http://vcg.isti.cnr.it/jgt/tetra.htm> for details.
+
+    if (s + t > 1)
+    {
+        s = 1 - s;
+        t = 1 - t;
+    }
+
+    if (s + t + u > 1)
+    {
+        Foam::scalar temp = u;
+
+        if (t + u > 1)
+        {
+            u = 1 - s - t;
+            t = 1 - temp;
+        }
+        else
+        {
+            u = s + t + u - 1;
+            s = 1 - t - temp;
+        }
+    }
+
+    return Foam::barycentric(1 - s - t - u, s, t, u);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::barycentric Foam::barycentric01(Random& rndGen)
+{
+    return
+        ::barycentric01
+        (
+            rndGen.scalar01(),
+            rndGen.scalar01(),
+            rndGen.scalar01()
+        );
+}
+
+
+Foam::barycentric Foam::barycentric01(cachedRandom& rndGen)
+{
+    return
+        ::barycentric01
+        (
+            rndGen.sample01<scalar>(),
+            rndGen.sample01<scalar>(),
+            rndGen.sample01<scalar>()
+        );
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/Barycentric/barycentric/barycentric.H b/src/OpenFOAM/primitives/Barycentric/barycentric/barycentric.H
new file mode 100644
index 0000000000..34cef96a1f
--- /dev/null
+++ b/src/OpenFOAM/primitives/Barycentric/barycentric/barycentric.H
@@ -0,0 +1,74 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016-2017 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::barycentric
+
+Description
+    A scalar version of the templated Barycentric
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef barycentric_H
+#define barycentric_H
+
+#include "scalar.H"
+#include "Barycentric.H"
+#include "contiguous.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class Random;
+class cachedRandom;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+typedef Barycentric<scalar> barycentric;
+
+
+//- Generate a random barycentric coordinate within the unit tetrahedron
+barycentric barycentric01(Random& rndGen);
+barycentric barycentric01(cachedRandom& rndGen);
+
+
+template<>
+inline bool contiguous<barycentric>()
+{
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/Barycentric/barycentricTensor/barycentricTensor.H b/src/OpenFOAM/primitives/Barycentric/barycentricTensor/barycentricTensor.H
new file mode 100644
index 0000000000..3bd3f06195
--- /dev/null
+++ b/src/OpenFOAM/primitives/Barycentric/barycentricTensor/barycentricTensor.H
@@ -0,0 +1,64 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016-2017 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::barycentricTensor
+
+Description
+    A scalar version of the templated BarycentricTensor
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef barycentricTensor_H
+#define barycentricTensor_H
+
+#include "scalar.H"
+#include "BarycentricTensor.H"
+#include "contiguous.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+typedef BarycentricTensor<scalar> barycentricTensor;
+
+
+template<>
+inline bool contiguous<barycentricTensor>()
+{
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.H b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.H
index db92376145..633fe479f7 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.H
+++ b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.H
@@ -92,6 +92,9 @@ public:
         //- Evaluate at x
         inline scalar value(const scalar x) const;
 
+        //- Evaluate the derivative at x
+        inline scalar derivative(const scalar x) const;
+
         //- Estimate the error of evaluation at x
         inline scalar error(const scalar x) const;
 
diff --git a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqnI.H b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqnI.H
index eb6a3283e4..da5fc3cba1 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqnI.H
+++ b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqnI.H
@@ -106,6 +106,12 @@ inline Foam::scalar Foam::cubicEqn::value(const scalar x) const
 }
 
 
+inline Foam::scalar Foam::cubicEqn::derivative(const scalar x) const
+{
+    return x*(x*3*a() + 2*b()) + c();
+}
+
+
 inline Foam::scalar Foam::cubicEqn::error(const scalar x) const
 {
     return mag(SMALL*x*(x*(x*3*a() + 2*b()) + c()));
diff --git a/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqn.H b/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqn.H
index debc52cc01..7f88336e63 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqn.H
+++ b/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqn.H
@@ -81,6 +81,9 @@ public:
         //- Evaluate at x
         inline scalar value(const scalar x) const;
 
+        //- Evaluate the derivative at x
+        inline scalar derivative(const scalar x) const;
+
         //- Estimate the error of evaluation at x
         inline scalar error(const scalar x) const;
 
diff --git a/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H b/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H
index e86ccbdd0a..17bfdf452d 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H
+++ b/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H
@@ -74,6 +74,12 @@ inline Foam::scalar Foam::linearEqn::value(const scalar x) const
 }
 
 
+inline Foam::scalar Foam::linearEqn::derivative(const scalar x) const
+{
+    return a();
+}
+
+
 inline Foam::scalar Foam::linearEqn::error(const scalar x) const
 {
     return mag(SMALL*x*a());
diff --git a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.H b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.H
index 1327dc15c9..c676227de1 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.H
+++ b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.H
@@ -84,6 +84,9 @@ public:
         //- Evaluate at x
         inline scalar value(const scalar x) const;
 
+        //- Evaluate the derivative at x
+        inline scalar derivative(const scalar x) const;
+
         //- Estimate the error of evaluation at x
         inline scalar error(const scalar x) const;
 
diff --git a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqnI.H b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqnI.H
index 4b4aacae61..7672a74d72 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqnI.H
+++ b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqnI.H
@@ -92,6 +92,12 @@ inline Foam::scalar Foam::quadraticEqn::value(const scalar x) const
 }
 
 
+inline Foam::scalar Foam::quadraticEqn::derivative(const scalar x) const
+{
+    return x*2*a() + b();
+}
+
+
 inline Foam::scalar Foam::quadraticEqn::error(const scalar x) const
 {
     return mag(SMALL*x*(x*2*a() + b()));
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
index ee586e2969..e7d6f558d1 100644
--- a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
+++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
@@ -727,7 +727,7 @@ void Foam::motionSmootherAlgo::movePoints()
 {
     // Make sure to clear out tetPtIs since used in checks (sometimes, should
     // really check)
-    mesh_.clearAdditionalGeom();
+    mesh_.clearTetBasePtIs();
     pp_.movePoints(mesh_.points());
 }
 
diff --git a/src/functionObjects/field/Make/files b/src/functionObjects/field/Make/files
index b44dc3819f..5161f46cf2 100644
--- a/src/functionObjects/field/Make/files
+++ b/src/functionObjects/field/Make/files
@@ -25,11 +25,6 @@ streamLine/streamLineBase.C
 streamLine/streamLineParticle.C
 streamLine/streamLineParticleCloud.C
 
-wallBoundedStreamLine/wallBoundedStreamLine.C
-wallBoundedStreamLine/wallBoundedStreamLineParticle.C
-wallBoundedStreamLine/wallBoundedStreamLineParticleCloud.C
-wallBoundedStreamLine/wallBoundedParticle.C
-
 surfaceInterpolate/surfaceInterpolate.C
 
 regionSizeDistribution/regionSizeDistribution.C
diff --git a/src/functionObjects/field/nearWallFields/findCellParticle.C b/src/functionObjects/field/nearWallFields/findCellParticle.C
index 6c69be4f18..e8efde1166 100644
--- a/src/functionObjects/field/nearWallFields/findCellParticle.C
+++ b/src/functionObjects/field/nearWallFields/findCellParticle.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,7 +30,7 @@ License
 Foam::findCellParticle::findCellParticle
 (
     const polyMesh& mesh,
-    const vector& position,
+    const barycentric& coordinates,
     const label celli,
     const label tetFacei,
     const label tetPti,
@@ -38,7 +38,24 @@ Foam::findCellParticle::findCellParticle
     const label data
 )
 :
-    particle(mesh, position, celli, tetFacei, tetPti),
+    particle(mesh, coordinates, celli, tetFacei, tetPti),
+    start_(position()),
+    end_(end),
+    data_(data)
+{}
+
+
+Foam::findCellParticle::findCellParticle
+(
+    const polyMesh& mesh,
+    const vector& position,
+    const label celli,
+    const point& end,
+    const label data
+)
+:
+    particle(mesh, position, celli),
+    start_(this->position()),
     end_(end),
     data_(data)
 {}
@@ -57,15 +74,15 @@ Foam::findCellParticle::findCellParticle
     {
         if (is.format() == IOstream::ASCII)
         {
-            is >> end_;
+            is >> start_ >> end_;
             data_ = readLabel(is);
         }
         else
         {
             is.read
             (
-                reinterpret_cast<char*>(&end_),
-                sizeof(end_) + sizeof(data_)
+                reinterpret_cast<char*>(&start_),
+                sizeof(start_) + sizeof(end_) + sizeof(data_)
             );
         }
     }
@@ -85,21 +102,13 @@ bool Foam::findCellParticle::move
     td.switchProcessor = false;
     td.keepParticle = true;
 
-    scalar tEnd = (1.0 - stepFraction())*maxTrackLen;
-    scalar dtMax = tEnd;
-
-    while (td.keepParticle && !td.switchProcessor && tEnd > SMALL)
+    while (td.keepParticle && !td.switchProcessor && stepFraction() < 1)
     {
-        // set the lagrangian time-step
-        scalar dt = min(dtMax, tEnd);
-
-        dt *= trackToFace(end_, td);
-
-        tEnd -= dt;
-        stepFraction() = 1.0 - tEnd/maxTrackLen;
+        const scalar f = 1 - stepFraction();
+        trackToFace(f*(end_ - start_), f, td);
     }
 
-    if (tEnd < SMALL || !td.keepParticle)
+    if (stepFraction() == 1 || !td.keepParticle)
     {
         // Hit endpoint or patch. If patch hit could do fancy stuff but just
         // to use the patch point is good enough for now.
@@ -209,6 +218,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const findCellParticle& p)
     if (os.format() == IOstream::ASCII)
     {
         os  << static_cast<const particle&>(p)
+            << token::SPACE << p.start_
             << token::SPACE << p.end_
             << token::SPACE << p.data_;
     }
@@ -217,8 +227,8 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const findCellParticle& p)
         os  << static_cast<const particle&>(p);
         os.write
         (
-            reinterpret_cast<const char*>(&p.end_),
-            sizeof(p.end_) + sizeof(p.data_)
+            reinterpret_cast<const char*>(&p.start_),
+            sizeof(p.start_) + sizeof(p.end_) + sizeof(p.data_)
         );
     }
 
diff --git a/src/functionObjects/field/nearWallFields/findCellParticle.H b/src/functionObjects/field/nearWallFields/findCellParticle.H
index d4bb541807..ccc44e9238 100644
--- a/src/functionObjects/field/nearWallFields/findCellParticle.H
+++ b/src/functionObjects/field/nearWallFields/findCellParticle.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -63,6 +63,9 @@ class findCellParticle
 {
     // Private data
 
+        //- Start point to track from
+        point start_;
+
         //- End point to track to
         point end_;
 
@@ -119,7 +122,7 @@ public:
         findCellParticle
         (
             const polyMesh& mesh,
-            const vector& position,
+            const barycentric& coordinates,
             const label celli,
             const label tetFacei,
             const label tetPti,
@@ -127,6 +130,17 @@ public:
             const label data
         );
 
+        //- Construct from a position and a cell, searching for the rest of the
+        //  required topology
+        findCellParticle
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label celli,
+            const point& end,
+            const label data
+        );
+
         //- Construct from Istream
         findCellParticle
         (
@@ -166,18 +180,42 @@ public:
 
     // Member Functions
 
+        //- Point to track from
+        const point& start() const
+        {
+            return start_;
+        }
+
+        //- Point to track from
+        point& start()
+        {
+            return start_;
+        }
+
         //- Point to track to
         const point& end() const
         {
             return end_;
         }
 
+        //- Point to track to
+        point& end()
+        {
+            return end_;
+        }
+
         //- Transported label
         label data() const
         {
             return data_;
         }
 
+        //- Transported label
+        label& data()
+        {
+            return data_;
+        }
+
 
         // Tracking
 
diff --git a/src/functionObjects/field/nearWallFields/nearWallFields.C b/src/functionObjects/field/nearWallFields/nearWallFields.C
index abf2e0d400..515f22763a 100644
--- a/src/functionObjects/field/nearWallFields/nearWallFields.C
+++ b/src/functionObjects/field/nearWallFields/nearWallFields.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -126,32 +126,19 @@ void Foam::functionObjects::nearWallFields::calcAddressing()
 
             const point end = start-distance_*nf[patchFacei];
 
-            if (tetFacei == -1)
-            {
-                WarningInFunction << "Did not find point " << start
-                    << " inside cell " << celli
-                    << ". Not seeding particle originating from face centre "
-                    << mesh_.faceCentres()[meshFacei]
-                    << " with cell centre " << mesh_.cellCentres()[celli]
-                    << endl;
-            }
-            else
-            {
-                // Add to cloud. Add originating face as passive data
-                cloud.addParticle
+            // Add a particle to the cloud with originating face as passive data
+            cloud.addParticle
+            (
+                new findCellParticle
                 (
-                    new findCellParticle
-                    (
-                        mesh_,
-                        start,
-                        celli,
-                        tetFacei,
-                        tetPti,
-                        end,
-                        globalWalls.toGlobal(nPatchFaces)    // passive data
-                    )
-                );
-            }
+                    mesh_,
+                    start,
+                    -1,
+                    end,
+                    globalWalls.toGlobal(nPatchFaces) // passive data
+                )
+            );
+
             nPatchFaces++;
         }
     }
diff --git a/src/functionObjects/field/streamLine/streamLineParticle.C b/src/functionObjects/field/streamLine/streamLineParticle.C
index 04189af71e..c587a7279c 100644
--- a/src/functionObjects/field/streamLine/streamLineParticle.C
+++ b/src/functionObjects/field/streamLine/streamLineParticle.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,35 +26,8 @@ License
 #include "streamLineParticle.H"
 #include "vectorFieldIOField.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-//    defineParticleTypeNameAndDebug(streamLineParticle, 0);
-}
-
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::scalar Foam::streamLineParticle::calcSubCycleDeltaT
-(
-    trackingData& td,
-    const scalar dt,
-    const vector& U
-) const
-{
-    particle testParticle(*this);
-
-    bool oldKeepParticle = td.keepParticle;
-    bool oldSwitchProcessor = td.switchProcessor;
-    scalar fraction = testParticle.trackToFace(position()+dt*U, td);
-    td.keepParticle = oldKeepParticle;
-    td.switchProcessor = oldSwitchProcessor;
-    // Adapt the dt to subdivide the trajectory into substeps.
-    return dt*fraction/td.nSubCycle_;
-}
-
-
 Foam::vector Foam::streamLineParticle::interpolateFields
 (
     const trackingData& td,
@@ -129,7 +102,6 @@ Foam::streamLineParticle::streamLineParticle
 {
     if (readFields)
     {
-        //if (is.format() == IOstream::ASCII)
         List<scalarList> sampledScalars;
         List<vectorList> sampledVectors;
 
@@ -169,31 +141,22 @@ Foam::streamLineParticle::streamLineParticle
 bool Foam::streamLineParticle::move
 (
     trackingData& td,
-    const scalar trackTime
+    const scalar
 )
 {
-    streamLineParticle& p = static_cast<streamLineParticle&>(*this);
-
     td.switchProcessor = false;
     td.keepParticle = true;
 
-    scalar tEnd = (1.0 - stepFraction())*trackTime;
-    scalar maxDt = mesh_.bounds().mag();
+    const scalar maxDt = mesh().bounds().mag();
 
-    while
-    (
-        td.keepParticle
-    && !td.switchProcessor
-    && lifeTime_ > 0
-    )
+    while (td.keepParticle && !td.switchProcessor && lifeTime_ > 0)
     {
-        // set the lagrangian time-step
         scalar dt = maxDt;
 
         // Cross cell in steps:
         // - at subiter 0 calculate dt to cross cell in nSubCycle steps
         // - at the last subiter do all of the remaining track
-        for (label subIter = 0; subIter < td.nSubCycle_; subIter++)
+        for (label subIter = 0; subIter < max(1, td.nSubCycle_); subIter++)
         {
             --lifeTime_;
 
@@ -219,37 +182,27 @@ bool Foam::streamLineParticle::move
 
             if (td.trackLength_ < GREAT)
             {
+                // No sub-cycling. Track a set length on each step.
                 dt = td.trackLength_;
-                //Pout<< "    subiteration " << subIter
-                //    << " : fixed length: updated dt:" << dt << endl;
             }
-            else if (subIter == 0 && td.nSubCycle_ > 1)
+            else if (subIter == 0)
             {
-                // Adapt dt to cross cell in a few steps
-                dt = calcSubCycleDeltaT(td, dt, U);
+                // Sub-cycling. Cross the cell in nSubCycle steps.
+                particle copy(*this);
+                copy.trackToFace(maxDt*U, 1);
+                dt *= (copy.stepFraction() - stepFraction())/td.nSubCycle_;
             }
             else if (subIter == td.nSubCycle_ - 1)
             {
-                // Do full step on last subcycle
+                // Sub-cycling. Track the whole cell on the last step.
                 dt = maxDt;
             }
 
-
-            scalar fraction = trackToFace(position() + dt*U, td);
-            dt *= fraction;
-
-            tEnd -= dt;
-            stepFraction() = 1.0 - tEnd/trackTime;
-
-            if (tEnd <= ROOTVSMALL)
-            {
-                // Force removal
-                lifeTime_ = 0;
-            }
+            trackToFace(dt*U, 0, td);
 
             if
             (
-                face() != -1
+                onFace()
             || !td.keepParticle
             ||  td.switchProcessor
             ||  lifeTime_ == 0
@@ -260,17 +213,16 @@ bool Foam::streamLineParticle::move
         }
     }
 
-
     if (!td.keepParticle || lifeTime_ == 0)
     {
         if (lifeTime_ == 0)
         {
+            // Failure exit. Particle stagnated or it's life ran out.
             if (debug)
             {
-                Pout<< "streamLineParticle : Removing stagnant particle:"
-                    << p.position()
-                    << " sampled positions:" << sampledPositions_.size()
-                    << endl;
+                Pout<< "streamLineParticle: Removing stagnant particle:"
+                    << position() << " sampled positions:"
+                    << sampledPositions_.size() << endl;
             }
             td.keepParticle = false;
         }
@@ -282,29 +234,25 @@ bool Foam::streamLineParticle::move
 
             if (debug)
             {
-                Pout<< "streamLineParticle : Removing particle:"
-                    << p.position()
+                Pout<< "streamLineParticle: Removing particle:" << 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]);
@@ -428,18 +376,11 @@ void Foam::streamLineParticle::readFields(Cloud<streamLineParticle>& c)
     );
     c.checkFieldIOobject(c, sampledPositions);
 
-//    vectorFieldIOField sampleVelocity
-//    (
-//        c.fieldIOobject("sampleVelocity", IOobject::MUST_READ)
-//    );
-//    c.checkFieldIOobject(c, sampleVelocity);
-
     label i = 0;
     forAllIter(Cloud<streamLineParticle>, c, iter)
     {
         iter().lifeTime_ = lifeTime[i];
         iter().sampledPositions_.transfer(sampledPositions[i]);
-//        iter().sampleVelocity_.transfer(sampleVelocity[i]);
         i++;
     }
 }
@@ -461,24 +402,17 @@ void Foam::streamLineParticle::writeFields(const Cloud<streamLineParticle>& c)
         c.fieldIOobject("sampledPositions", IOobject::NO_READ),
         np
     );
-//    vectorFieldIOField sampleVelocity
-//    (
-//        c.fieldIOobject("sampleVelocity", IOobject::NO_READ),
-//        np
-//    );
 
     label i = 0;
     forAllConstIter(Cloud<streamLineParticle>, c, iter)
     {
         lifeTime[i] = iter().lifeTime_;
         sampledPositions[i] = iter().sampledPositions_;
-//        sampleVelocity[i] = iter().sampleVelocity_;
         i++;
     }
 
     lifeTime.write();
     sampledPositions.write();
-//    sampleVelocity.write();
 }
 
 
diff --git a/src/functionObjects/field/streamLine/streamLineParticle.H b/src/functionObjects/field/streamLine/streamLineParticle.H
index bcb4fd68c2..42282e1790 100644
--- a/src/functionObjects/field/streamLine/streamLineParticle.H
+++ b/src/functionObjects/field/streamLine/streamLineParticle.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -46,16 +46,12 @@ SourceFiles
 namespace Foam
 {
 
-class streamLineParticleCloud;
-
-
 // Forward declaration of friend functions and operators
 
 class streamLineParticle;
 
 Ostream& operator<<(Ostream&, const streamLineParticle&);
 
-
 /*---------------------------------------------------------------------------*\
                      Class streamLineParticle Declaration
 \*---------------------------------------------------------------------------*/
@@ -64,32 +60,38 @@ class streamLineParticle
 :
     public particle
 {
-
 public:
 
-    //- Class used to pass tracking data to the trackToFace function
     class trackingData
     :
         public particle::TrackingData<Cloud<streamLineParticle>>
     {
-
     public:
 
+        // Public data
 
-        const PtrList<interpolation<scalar>>& vsInterp_;
-        const PtrList<interpolation<vector>>& vvInterp_;
-        const label UIndex_;
-        const bool trackForward_;
-        const label nSubCycle_;
-        const scalar trackLength_;
+            const PtrList<interpolation<scalar>>& vsInterp_;
 
-        DynamicList<vectorList>& allPositions_;
-        List<DynamicList<scalarList>>& allScalars_;
-        List<DynamicList<vectorList>>& allVectors_;
+            const PtrList<interpolation<vector>>& vvInterp_;
+
+            const label UIndex_;
+
+            const bool trackForward_;
+
+            const label nSubCycle_;
+
+            const scalar trackLength_;
+
+            DynamicList<vectorList>& allPositions_;
+
+            List<DynamicList<scalarList>>& allScalars_;
+
+            List<DynamicList<vectorList>>& allVectors_;
 
 
         // Constructors
 
+            //- Construct from components
             trackingData
             (
                 Cloud<streamLineParticle>& cloud,
@@ -99,7 +101,6 @@ public:
                 const bool trackForward,
                 const label nSubCycle,
                 const scalar trackLength,
-
                 DynamicList<List<point>>& allPositions,
                 List<DynamicList<scalarList>>& allScalars,
                 List<DynamicList<vectorList>>& allVectors
@@ -112,7 +113,6 @@ public:
                 trackForward_(trackForward),
                 nSubCycle_(nSubCycle),
                 trackLength_(trackLength),
-
                 allPositions_(allPositions),
                 allScalars_(allScalars),
                 allVectors_(allVectors)
@@ -139,22 +139,6 @@ private:
 
     // Private Member Functions
 
-        //- Estimate dt to cross from current face to next one in nSubCycle
-        //  steps.
-        scalar calcSubCycleDeltaT
-        (
-            trackingData& td,
-            const scalar dt,
-            const vector& U
-        ) const;
-
-        void constrainVelocity
-        (
-            trackingData& td,
-            const scalar dt,
-            vector& U
-        );
-
         //- Interpolate all quantities; return interpolated velocity.
         vector interpolateFields
         (
@@ -195,8 +179,7 @@ public:
             return autoPtr<particle>(new streamLineParticle(*this));
         }
 
-        //- Factory class to read-construct particles used for
-        //  parallel transfer
+        //- Factory class to read-construct particles used for parallel transfer
         class iNew
         {
             const polyMesh& mesh_;
@@ -223,8 +206,7 @@ public:
         // Tracking
 
             //- Track all particles to their end point
-            bool move(trackingData&, const scalar trackTime);
-
+            bool move(trackingData&, const scalar);
 
             //- Overridable function to handle the particle hitting a patch
             //  Executed before other patch-hitting functions
diff --git a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C
index fadaf872bc..f16a25aa50 100644
--- a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C
+++ b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -176,16 +176,7 @@ void Foam::DSMCCloud<ParcelType>::initialise
 
                     U += velocity;
 
-                    addNewParcel
-                    (
-                        p,
-                        U,
-                        Ei,
-                        celli,
-                        cellTetIs.face(),
-                        cellTetIs.tetPt(),
-                        typeId
-                    );
+                    addNewParcel(p, celli, U, Ei, typeId);
                 }
             }
         }
@@ -465,27 +456,13 @@ template<class ParcelType>
 void Foam::DSMCCloud<ParcelType>::addNewParcel
 (
     const vector& position,
+    const label celli,
     const vector& U,
     const scalar Ei,
-    const label celli,
-    const label tetFacei,
-    const label tetPti,
     const label typeId
 )
 {
-    ParcelType* pPtr = new ParcelType
-    (
-        mesh_,
-        position,
-        U,
-        Ei,
-        celli,
-        tetFacei,
-        tetPti,
-        typeId
-    );
-
-    this->addParticle(pPtr);
+    this->addParticle(new ParcelType(mesh_, position, celli, U, Ei, typeId));
 }
 
 
@@ -1110,9 +1087,7 @@ void Foam::DSMCCloud<ParcelType>::dumpParticlePositions() const
 template<class ParcelType>
 void Foam::DSMCCloud<ParcelType>::autoMap(const mapPolyMesh& mapper)
 {
-    typedef typename  ParcelType::trackingData tdType;
-    tdType td(*this);
-    Cloud<ParcelType>::template autoMap<tdType>(td, mapper);
+    Cloud<ParcelType>::autoMap(mapper);
 
     // Update the cell occupancy field
     cellOccupancy_.setSize(mesh_.nCells());
diff --git a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.H b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.H
index d7ea2be9b4..81bbc2ff64 100644
--- a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.H
+++ b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -44,6 +44,7 @@ SourceFiles
 #include "fvMesh.H"
 #include "volFields.H"
 #include "scalarIOField.H"
+#include "barycentric.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -451,11 +452,9 @@ public:
             void addNewParcel
             (
                 const vector& position,
+                const label celli,
                 const vector& U,
                 const scalar Ei,
-                const label celli,
-                const label tetFacei,
-                const label tetPti,
                 const label typeId
             );
 
diff --git a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C
index 554f3b52ba..3df37ae160 100644
--- a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C
+++ b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -41,19 +41,16 @@ bool Foam::DSMCParcel<ParcelType>::move(TrackData& td, const scalar trackTime)
     const polyMesh& mesh = td.cloud().pMesh();
     const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
 
-    scalar tEnd = (1.0 - p.stepFraction())*trackTime;
-    const scalar dtMax = tEnd;
-
     // For reduced-D cases, the velocity used to track needs to be
     // constrained, but the actual U_ of the parcel must not be
     // altered or used, as it is altered by patch interactions an
     // needs to retain its 3D value for collision purposes.
     vector Utracking = U_;
 
-    while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
+    while (td.keepParticle && !td.switchProcessor && p.stepFraction() < 1)
     {
         // Apply correction to position for reduced-D cases
-        meshTools::constrainToMeshCentre(mesh, p.position());
+        p.constrainToMeshCentre();
 
         Utracking = U_;
 
@@ -61,16 +58,10 @@ bool Foam::DSMCParcel<ParcelType>::move(TrackData& td, const scalar trackTime)
         // reduced-D cases
         meshTools::constrainDirection(mesh, mesh.solutionD(), Utracking);
 
-        // Set the Lagrangian time-step
-        scalar dt = min(dtMax, tEnd);
+        const scalar f = 1 - p.stepFraction();
+        p.trackToFace(f*trackTime*Utracking, f, td);
 
-        dt *= p.trackToFace(p.position() + dt*Utracking, td);
-
-        tEnd -= dt;
-
-        p.stepFraction() = 1.0 - tEnd/trackTime;
-
-        if (p.onBoundary() && td.keepParticle)
+        if (p.onBoundaryFace() && td.keepParticle)
         {
             if (isA<processorPolyPatch>(pbMesh[p.patch(p.face())]))
             {
diff --git a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.H b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.H
index 453cb74315..b3525cd7db 100644
--- a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.H
+++ b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -176,12 +176,24 @@ public:
         inline DSMCParcel
         (
             const polyMesh& mesh,
-            const vector& position,
-            const vector& U,
-            const scalar Ei,
+            const barycentric& coordinates,
             const label celli,
             const label tetFacei,
             const label tetPti,
+            const vector& U,
+            const scalar Ei,
+            const label typeId
+        );
+
+        //- Construct from a position and a cell, searching for the rest of the
+        //  required topology
+        inline DSMCParcel
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label celli,
+            const vector& U,
+            const scalar Ei,
             const label typeId
         );
 
diff --git a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcelI.H b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcelI.H
index 5fb82c5cda..2bff944e14 100644
--- a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcelI.H
+++ b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcelI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -55,16 +55,34 @@ template<class ParcelType>
 inline Foam::DSMCParcel<ParcelType>::DSMCParcel
 (
     const polyMesh& mesh,
-    const vector& position,
-    const vector& U,
-    const scalar Ei,
+    const barycentric& coordinates,
     const label celli,
     const label tetFacei,
     const label tetPti,
+    const vector& U,
+    const scalar Ei,
     const label typeId
 )
 :
-    ParcelType(mesh, position, celli, tetFacei, tetPti),
+    ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
+    U_(U),
+    Ei_(Ei),
+    typeId_(typeId)
+{}
+
+
+template<class ParcelType>
+inline Foam::DSMCParcel<ParcelType>::DSMCParcel
+(
+    const polyMesh& mesh,
+    const vector& position,
+    const label celli,
+    const vector& U,
+    const scalar Ei,
+    const label typeId
+)
+:
+    ParcelType(mesh, position, celli),
     U_(U),
     Ei_(Ei),
     typeId_(typeId)
diff --git a/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C b/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
index 4ed53b5cf6..3b85d70587 100644
--- a/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
+++ b/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -403,16 +403,7 @@ void Foam::FreeStream<CloudType>::inflow()
                         cloud.constProps(typeId).internalDegreesOfFreedom()
                     );
 
-                    cloud.addNewParcel
-                    (
-                        p,
-                        U,
-                        Ei,
-                        celli,
-                        globalFaceIndex,
-                        faceTetIs.tetPt(),
-                        typeId
-                    );
+                    cloud.addNewParcel(p, celli, U, Ei, typeId);
 
                     particlesInserted++;
                 }
diff --git a/src/lagrangian/basic/Cloud/Cloud.C b/src/lagrangian/basic/Cloud/Cloud.C
index 9b66f53706..f0f7da835c 100644
--- a/src/lagrangian/basic/Cloud/Cloud.C
+++ b/src/lagrangian/basic/Cloud/Cloud.C
@@ -105,7 +105,8 @@ Foam::Cloud<ParticleType>::Cloud
     polyMesh_(pMesh),
     labels_(),
     nTrackingRescues_(),
-    cellWallFacesPtr_()
+    cellWallFacesPtr_(),
+    globalPositionsPtr_()
 {
     checkPatches();
 
@@ -234,6 +235,8 @@ void Foam::Cloud<ParticleType>::move(TrackData& td, const scalar trackTime)
     // Allocate transfer buffers
     PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
 
+    // Clear the global positions as there are about to change
+    globalPositionsPtr_.clear();
 
     // While there are particles to transfer
     while (true)
@@ -393,21 +396,16 @@ void Foam::Cloud<ParticleType>::move(TrackData& td, const scalar trackTime)
 
 
 template<class ParticleType>
-template<class TrackData>
-void Foam::Cloud<ParticleType>::autoMap
-(
-    TrackData& td,
-    const mapPolyMesh& mapper
-)
+void Foam::Cloud<ParticleType>::autoMap(const mapPolyMesh& mapper)
 {
-    if (cloud::debug)
+    if (!globalPositionsPtr_.valid())
     {
-        InfoInFunction << "for lagrangian cloud " << cloud::name() << endl;
+        FatalErrorInFunction
+            << "Global positions are not available. "
+            << "Cloud::storeGlobalPositions has not been called."
+            << exit(FatalError);
     }
 
-    const labelList& reverseCellMap = mapper.reverseCellMap();
-    const labelList& reverseFaceMap = mapper.reverseFaceMap();
-
     // Reset stored data that relies on the mesh
     //    polyMesh_.clearCellTree();
     cellWallFacesPtr_.clear();
@@ -417,51 +415,13 @@ void Foam::Cloud<ParticleType>::autoMap
     // there is a comms mismatch.
     polyMesh_.tetBasePtIs();
 
+    const vectorField& positions = globalPositionsPtr_();
 
-    forAllIter(typename Cloud<ParticleType>, *this, pIter)
+    label i = 0;
+    forAllIter(typename Cloud<ParticleType>, *this, iter)
     {
-        ParticleType& p = pIter();
-
-        if (reverseCellMap[p.cell()] >= 0)
-        {
-            p.cell() = reverseCellMap[p.cell()];
-
-            if (p.face() >= 0 && reverseFaceMap[p.face()] >= 0)
-            {
-                p.face() = reverseFaceMap[p.face()];
-            }
-            else
-            {
-                p.face() = -1;
-            }
-
-            p.initCellFacePt();
-        }
-        else
-        {
-            label trackStartCell = mapper.mergedCell(p.cell());
-
-            if (trackStartCell < 0)
-            {
-                trackStartCell = 0;
-                p.cell() = 0;
-            }
-            else
-            {
-                p.cell() = trackStartCell;
-            }
-
-            vector pos = p.position();
-
-            const_cast<vector&>(p.position()) =
-                polyMesh_.cellCentres()[trackStartCell];
-
-            p.stepFraction() = 0;
-
-            p.initCellFacePt();
-
-            p.track(pos, td);
-        }
+        iter().autoMap(positions[i], mapper);
+        ++ i;
     }
 }
 
@@ -485,6 +445,27 @@ void Foam::Cloud<ParticleType>::writePositions() const
 }
 
 
+template<class ParticleType>
+void Foam::Cloud<ParticleType>::storeGlobalPositions() const
+{
+    // Store the global positions for later use by autoMap. It would be
+    // preferable not to need this. If the mapPolyMesh object passed to autoMap
+    // had a copy of the old mesh then the global positions could be recovered
+    // within autoMap, and this pre-processing would not be necessary.
+
+    globalPositionsPtr_.reset(new vectorField(this->size()));
+
+    vectorField& positions = globalPositionsPtr_();
+
+    label i = 0;
+    forAllConstIter(typename Cloud<ParticleType>, *this, iter)
+    {
+        positions[i] = iter().position();
+        ++ i;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
 
 #include "CloudIO.C"
diff --git a/src/lagrangian/basic/Cloud/Cloud.H b/src/lagrangian/basic/Cloud/Cloud.H
index 2ad874691b..27be0c3a69 100644
--- a/src/lagrangian/basic/Cloud/Cloud.H
+++ b/src/lagrangian/basic/Cloud/Cloud.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -87,6 +87,9 @@ class Cloud
         //- Does the cell have wall faces
         mutable autoPtr<PackedBoolList> cellWallFacesPtr_;
 
+        //- Temporary storage for the global particle positions
+        mutable autoPtr<vectorField> globalPositionsPtr_;
+
 
     // Private Member Functions
 
@@ -162,7 +165,7 @@ public:
                 return IDLList<ParticleType>::size();
             };
 
-            DynamicList<label>& labels()
+            DynamicList<label>& labels() const
             {
                 return labels_;
             }
@@ -256,8 +259,7 @@ public:
 
             //- Remap the cells of particles corresponding to the
             //  mesh topology change
-            template<class TrackData>
-            void autoMap(TrackData& td, const mapPolyMesh&);
+            void autoMap(const mapPolyMesh&);
 
 
         // Read
@@ -304,6 +306,10 @@ public:
             //- Write positions to \<cloudName\>_positions.obj file
             void writePositions() const;
 
+            //- Call this before a topology change. Stores the particles global
+            //  positions in the database for use during mapping.
+            void storeGlobalPositions() const;
+
 
     // Ostream Operator
 
diff --git a/src/lagrangian/basic/Cloud/CloudIO.C b/src/lagrangian/basic/Cloud/CloudIO.C
index 83369d1a4f..7e5131fd9b 100644
--- a/src/lagrangian/basic/Cloud/CloudIO.C
+++ b/src/lagrangian/basic/Cloud/CloudIO.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -132,13 +132,6 @@ void Foam::Cloud<ParticleType>::initCloud(const bool checkClass)
     // them, otherwise, if some processors have no particles then
     // there is a comms mismatch.
     polyMesh_.tetBasePtIs();
-
-    forAllIter(typename Cloud<ParticleType>, *this, pIter)
-    {
-        ParticleType& p = pIter();
-
-        p.initCellFacePt();
-    }
 }
 
 
diff --git a/src/lagrangian/basic/InteractionLists/InteractionLists.C b/src/lagrangian/basic/InteractionLists/InteractionLists.C
index 4223cd6778..b223232bd3 100644
--- a/src/lagrangian/basic/InteractionLists/InteractionLists.C
+++ b/src/lagrangian/basic/InteractionLists/InteractionLists.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -965,14 +965,7 @@ void Foam::InteractionLists<ParticleType>::prepareParticleToBeReferred
         globalTransforms.transformIndex(ciat)
     );
 
-    particle->position() = transform.invTransformPosition(particle->position());
-
-    particle->transformProperties(-transform.t());
-
-    if (transform.hasR())
-    {
-        particle->transformProperties(transform.R().T());
-    }
+    particle->prepareForInteractionListReferral(transform);
 }
 
 
@@ -1230,6 +1223,15 @@ void Foam::InteractionLists<ParticleType>::receiveReferredData
         }
     }
 
+    forAll(referredParticles_, refCelli)
+    {
+        IDLList<ParticleType>& refCell = referredParticles_[refCelli];
+        forAllIter(typename IDLList<ParticleType>, refCell, iter)
+        {
+            iter().correctAfterInteractionListReferral(ril_[refCelli][0]);
+        }
+    }
+
     fillReferredParticleCloud();
 
     wallFaceMap().receive(pBufs, referredWallData_);
diff --git a/src/lagrangian/basic/indexedParticle/indexedParticle.H b/src/lagrangian/basic/indexedParticle/indexedParticle.H
index 75ed0783e4..c619cc6694 100644
--- a/src/lagrangian/basic/indexedParticle/indexedParticle.H
+++ b/src/lagrangian/basic/indexedParticle/indexedParticle.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -65,27 +65,14 @@ public:
         indexedParticle
         (
             const polyMesh& mesh,
-            const vector& position,
+            const barycentric& coordinates,
             const label celli,
             const label tetFacei,
             const label tetPti,
             const label index = 0
         )
         :
-            particle(mesh, position, celli, tetFacei, tetPti),
-            index_(index)
-        {}
-
-        //- Construct from components, with searching for tetFace and tetPt
-        indexedParticle
-        (
-            const polyMesh& mesh,
-            const vector& position,
-            const label celli,
-            const label index = 0
-        )
-        :
-            particle(mesh, position, celli),
+            particle(mesh, coordinates, celli, tetFacei, tetPti),
             index_(index)
         {}
 
diff --git a/src/lagrangian/basic/particle/particle.C b/src/lagrangian/basic/particle/particle.C
index 88de1e7dda..ea1a1a07a3 100644
--- a/src/lagrangian/basic/particle/particle.C
+++ b/src/lagrangian/basic/particle/particle.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,41 +25,570 @@ License
 
 #include "particle.H"
 #include "transform.H"
+#include "treeDataCell.H"
+#include "cubicEqn.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 Foam::label Foam::particle::particleCount_ = 0;
 
-const Foam::scalar Foam::particle::trackingCorrectionTol = 1e-5;
-
-const Foam::scalar Foam::particle::lambdaDistanceToleranceCoeff = 1e3*SMALL;
-
-const Foam::scalar Foam::particle::minStepFractionTol = 1e5*SMALL;
-
 namespace Foam
 {
     defineTypeNameAndDebug(particle, 0);
 }
 
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::particle::tetFaceIndices
+(
+    label& baseI,
+    label& vertex1I,
+    label& vertex2I
+) const
+{
+    const Foam::face& f = mesh_.faces()[tetFacei_];
+
+    baseI = max(0, mesh_.tetBasePtIs()[tetFacei_]);
+
+    vertex1I = (baseI + tetPti_) % f.size();
+
+    vertex2I = f.fcIndex(vertex1I);
+
+    if (mesh_.faceOwner()[tetFacei_] != celli_)
+    {
+        Swap(vertex1I, vertex2I);
+    }
+}
+
+
+void Foam::particle::tetMeshIndices
+(
+    label& basei,
+    label& vertex1i,
+    label& vertex2i
+) const
+{
+    const Foam::face& f = mesh_.faces()[tetFacei_];
+
+    tetFaceIndices(basei, vertex1i, vertex2i);
+
+    basei = f[basei];
+    vertex1i = f[vertex1i];
+    vertex2i = f[vertex2i];
+}
+
+
+void Foam::particle::tetGeometry
+(
+    vector& centre,
+    vector& base,
+    vector& vertex1,
+    vector& vertex2
+) const
+{
+    label basei, vertex1i, vertex2i;
+    tetMeshIndices(basei, vertex1i, vertex2i);
+
+    centre = mesh_.cellCentres()[celli_];
+    base = mesh_.points()[basei];
+    vertex1 = mesh_.points()[vertex1i];
+    vertex2 = mesh_.points()[vertex2i];
+}
+
+
+Foam::barycentricTensor Foam::particle::tetTransform() const
+{
+    vector centre, base, vertex1, vertex2;
+    tetGeometry(centre, base, vertex1, vertex2);
+
+    return barycentricTensor(centre, base, vertex1, vertex2);
+}
+
+
+void Foam::particle::tetReverseTransform
+(
+    vector& centre,
+    scalar& detA,
+    barycentricTensor& T
+) const
+{
+    barycentricTensor A = tetTransform();
+
+    const vector ab = A.b() - A.a();
+    const vector ac = A.c() - A.a();
+    const vector ad = A.d() - A.a();
+    const vector bc = A.c() - A.b();
+    const vector bd = A.d() - A.b();
+
+    centre = A.a();
+
+    detA = ab & (ac ^ ad);
+
+    T = barycentricTensor
+    (
+        bd ^ bc,
+        ac ^ ad,
+        ad ^ ab,
+        ab ^ ac
+    );
+}
+
+
+void Foam::particle::movingTetGeometry
+(
+    const scalar fraction,
+    Pair<vector>& centre,
+    Pair<vector>& base,
+    Pair<vector>& vertex1,
+    Pair<vector>& vertex2
+) const
+{
+    label basei, vertex1i, vertex2i;
+    tetMeshIndices(basei, vertex1i, vertex2i);
+
+    const pointField& ptsOld = mesh_.oldPoints();
+    const pointField& ptsNew = mesh_.points();
+
+    // !!! <-- We would be better off using mesh_.cellCentres() here. However,
+    // we need to put a mesh_.oldCellCentres() method in for this to work. The
+    // values obtained from the mesh and those obtained from the cell do not
+    // necessarily match. See mantis #1993.
+    const vector ccOld = mesh_.cells()[celli_].centre(ptsOld, mesh_.faces());
+    const vector ccNew = mesh_.cells()[celli_].centre(ptsNew, mesh_.faces());
+
+    scalar f0 = stepFraction_, f1 = fraction;
+
+    // Old and new points and cell centres are not sub-cycled. If we are sub-
+    // cycling, then we have to account for the timestep change here by
+    // modifying the fractions that we take of the old and new geometry.
+    if (mesh_.time().subCycling())
+    {
+        const TimeState& tsNew = mesh_.time();
+        const TimeState& tsOld = mesh_.time().prevTimeState();
+
+        const scalar tFrac =
+        (
+            (tsNew.value() - tsNew.deltaTValue())
+          - (tsOld.value() - tsOld.deltaTValue())
+        )/tsOld.deltaTValue();
+
+        const scalar dtFrac = tsNew.deltaTValue()/tsOld.deltaTValue();
+
+        f0 = tFrac + dtFrac*f0;
+        f1 = dtFrac*f1;
+    }
+
+    centre[0] = ccOld + f0*(ccNew - ccOld);
+    base[0] = ptsOld[basei] + f0*(ptsNew[basei] - ptsOld[basei]);
+    vertex1[0] = ptsOld[vertex1i] + f0*(ptsNew[vertex1i] - ptsOld[vertex1i]);
+    vertex2[0] = ptsOld[vertex2i] + f0*(ptsNew[vertex2i] - ptsOld[vertex2i]);
+
+    centre[1] = f1*(ccNew - ccOld);
+    base[1] = f1*(ptsNew[basei] - ptsOld[basei]);
+    vertex1[1] = f1*(ptsNew[vertex1i] - ptsOld[vertex1i]);
+    vertex2[1] = f1*(ptsNew[vertex2i] - ptsOld[vertex2i]);
+}
+
+
+Foam::Pair<Foam::barycentricTensor> Foam::particle::movingTetTransform
+(
+    const scalar fraction
+) const
+{
+    Pair<vector> centre, base, vertex1, vertex2;
+    movingTetGeometry(fraction, centre, base, vertex1, vertex2);
+
+    return
+        Pair<barycentricTensor>
+        (
+            barycentricTensor(centre[0], base[0], vertex1[0], vertex2[0]),
+            barycentricTensor(centre[1], base[1], vertex1[1], vertex2[1])
+        );
+}
+
+
+void Foam::particle::movingTetReverseTransform
+(
+    const scalar fraction,
+    Pair<vector>& centre,
+    FixedList<scalar, 4>& detA,
+    FixedList<barycentricTensor, 3>& T
+) const
+{
+    Pair<barycentricTensor> A = movingTetTransform(fraction);
+
+    const Pair<vector> ab(A[0].b() - A[0].a(), A[1].b() - A[1].a());
+    const Pair<vector> ac(A[0].c() - A[0].a(), A[1].c() - A[1].a());
+    const Pair<vector> ad(A[0].d() - A[0].a(), A[1].d() - A[1].a());
+    const Pair<vector> bc(A[0].c() - A[0].b(), A[1].c() - A[1].b());
+    const Pair<vector> bd(A[0].d() - A[0].b(), A[1].d() - A[1].b());
+
+    centre[0] = A[0].a();
+    centre[1] = A[1].a();
+
+    detA[0] = ab[0] & (ac[0] ^ ad[0]);
+    detA[1] =
+        (ab[1] & (ac[0] ^ ad[0]))
+      + (ab[0] & (ac[1] ^ ad[0]))
+      + (ab[0] & (ac[0] ^ ad[1]));
+    detA[2] =
+        (ab[0] & (ac[1] ^ ad[1]))
+      + (ab[1] & (ac[0] ^ ad[1]))
+      + (ab[1] & (ac[1] ^ ad[0]));
+    detA[3] = ab[1] & (ac[1] ^ ad[1]);
+
+    T[0] = barycentricTensor
+    (
+        bd[0] ^ bc[0],
+        ac[0] ^ ad[0],
+        ad[0] ^ ab[0],
+        ab[0] ^ ac[0]
+    );
+    T[1] = barycentricTensor
+    (
+        (bd[0] ^ bc[1]) + (bd[1] ^ bc[0]),
+        (ac[0] ^ ad[1]) + (ac[1] ^ ad[0]),
+        (ad[0] ^ ab[1]) + (ad[1] ^ ab[0]),
+        (ab[0] ^ ac[1]) + (ab[1] ^ ac[0])
+    );
+    T[2] = barycentricTensor
+    (
+        bd[1] ^ bc[1],
+        ac[1] ^ ad[1],
+        ad[1] ^ ab[1],
+        ab[1] ^ ac[1]
+    );
+}
+
+
+void Foam::particle::reflect()
+{
+    Swap(coordinates_.c(), coordinates_.d());
+}
+
+
+void Foam::particle::rotate(const bool reverse)
+{
+    if (!reverse)
+    {
+        scalar temp = coordinates_.b();
+        coordinates_.b() = coordinates_.c();
+        coordinates_.c() = coordinates_.d();
+        coordinates_.d() = temp;
+    }
+    else
+    {
+        scalar temp = coordinates_.d();
+        coordinates_.d() = coordinates_.c();
+        coordinates_.c() = coordinates_.b();
+        coordinates_.b() = temp;
+    }
+}
+
+
+void Foam::particle::changeTet(const label tetTriI)
+{
+    const bool isOwner = mesh_.faceOwner()[tetFacei_] == celli_;
+
+    const label firstTetPtI = 1;
+    const label lastTetPtI = mesh_.faces()[tetFacei_].size() - 2;
+
+    if (tetTriI == 1)
+    {
+        changeFace(tetTriI);
+    }
+    else if (tetTriI == 2)
+    {
+        if (isOwner)
+        {
+            if (tetPti_ == lastTetPtI)
+            {
+                changeFace(tetTriI);
+            }
+            else
+            {
+                reflect();
+                tetPti_ += 1;
+            }
+        }
+        else
+        {
+            if (tetPti_ == firstTetPtI)
+            {
+                changeFace(tetTriI);
+            }
+            else
+            {
+                reflect();
+                tetPti_ -= 1;
+            }
+        }
+    }
+    else if (tetTriI == 3)
+    {
+        if (isOwner)
+        {
+            if (tetPti_ == firstTetPtI)
+            {
+                changeFace(tetTriI);
+            }
+            else
+            {
+                reflect();
+                tetPti_ -= 1;
+            }
+        }
+        else
+        {
+            if (tetPti_ == lastTetPtI)
+            {
+                changeFace(tetTriI);
+            }
+            else
+            {
+                reflect();
+                tetPti_ += 1;
+            }
+        }
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Changing tet without changing cell should only happen when the "
+            << "track is on triangle 1, 2 or 3."
+            << exit(FatalError);
+    }
+}
+
+
+void Foam::particle::changeFace(const label tetTriI)
+{
+    // Get the tet topology
+    label basei, vertex1i, vertex2i;
+    tetMeshIndices(basei, vertex1i, vertex2i);
+
+    // Get the shared edge and the pre-rotation
+    edge sharedEdge;
+    if (tetTriI == 1)
+    {
+        sharedEdge = edge(vertex1i, vertex2i);
+    }
+    else if (tetTriI == 2)
+    {
+        sharedEdge = edge(vertex2i, basei);
+    }
+    else if (tetTriI == 3)
+    {
+        sharedEdge = edge(basei, vertex1i);
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Changing face without changing cell should only happen when the"
+            << " track is on triangle 1, 2 or 3."
+            << exit(FatalError);
+    }
+
+    // Find the face in the same cell that shares the edge, and the
+    // corresponding tetrahedra point
+    tetPti_ = -1;
+    forAll(mesh_.cells()[celli_], cellFaceI)
+    {
+        const label newFaceI = mesh_.cells()[celli_][cellFaceI];
+        const class face& newFace = mesh_.faces()[newFaceI];
+
+        // Exclude the current face
+        if (tetFacei_ == newFaceI)
+        {
+            continue;
+        }
+
+        // Loop over the edges, looking for the shared one
+        label edgeComp = 0;
+        label edgeI = 0;
+        for (; edgeI < newFace.size(); ++ edgeI)
+        {
+            edgeComp = edge::compare(sharedEdge, newFace.faceEdge(edgeI));
+
+            if (edgeComp != 0)
+            {
+                break;
+            }
+        }
+
+        // If the face does not contain the edge, then move on to the next face
+        if (edgeComp == 0)
+        {
+            continue;
+        }
+
+        // Make the edge index relative to the base point
+        const label newBaseI = max(0, mesh_.tetBasePtIs()[newFaceI]);
+        edgeI = (edgeI - newBaseI + newFace.size()) % newFace.size();
+
+        // If the edge is next the base point (i.e., the index is 0 or n - 1),
+        // then we swap it for the adjacent edge. This new edge is opposite the
+        // base point, and defines the tet with the original edge in it.
+        edgeI = min(max(1, edgeI), newFace.size() - 2);
+
+        // Set the new face and tet point
+        tetFacei_ = newFaceI;
+        tetPti_ = edgeI;
+
+        // Exit the loop now that the the tet point has been found
+        break;
+    }
+
+    if (tetPti_ == -1)
+    {
+        FatalErrorInFunction
+            << "The search for an edge-connected face and tet-point failed."
+            << exit(FatalError);
+    }
+
+    // Pre-rotation puts the shared edge opposite the base of the tetrahedron
+    if (sharedEdge.otherVertex(vertex1i) == -1)
+    {
+        rotate(false);
+    }
+    else if (sharedEdge.otherVertex(vertex2i) == -1)
+    {
+        rotate(true);
+    }
+
+    // Update the new tet topology
+    tetMeshIndices(basei, vertex1i, vertex2i);
+
+    // Reflect to account for the change of triangle orientation on the new face
+    reflect();
+
+    // Post rotation puts the shared edge back in the correct location
+    if (sharedEdge.otherVertex(vertex1i) == -1)
+    {
+        rotate(true);
+    }
+    else if (sharedEdge.otherVertex(vertex2i) == -1)
+    {
+        rotate(false);
+    }
+}
+
+
+void Foam::particle::changeCell()
+{
+    // Set the cell to be the one on the other side of the face
+    const label ownerCellI = mesh_.faceOwner()[tetFacei_];
+    const bool isOwner = celli_ == ownerCellI;
+    celli_ = isOwner ? mesh_.faceNeighbour()[tetFacei_] : ownerCellI;
+
+    // Reflect to account for the change of triangle orientation in the new cell
+    reflect();
+}
+
+
+void Foam::particle::locate
+(
+    const vector& position,
+    const vector* direction,
+    const label celli,
+    const bool boundaryFail,
+    const string boundaryMsg
+)
+{
+    celli_ = celli;
+
+    // Find the cell, if it has not been given
+    if (celli_ < 0)
+    {
+        celli_ = mesh_.cellTree().findInside(position);
+    }
+    if (celli_ < 0)
+    {
+        FatalErrorInFunction
+            << "Cell not found for particle position " << position << "."
+            << exit(FatalError);
+    }
+
+    // Put the particle at the cell centre and in a random tet
+    coordinates_ = barycentric(1, 0, 0, 0);
+    tetFacei_ = mesh_.cells()[celli_][0];
+    tetPti_ = 1;
+    facei_ = -1;
+
+    // Track to the injection point
+    track(position - mesh_.cellCentres()[celli_], 0);
+    if (!onFace())
+    {
+        return;
+    }
+
+    // We hit a boundary ...
+    if (boundaryFail)
+    {
+        FatalErrorInFunction << boundaryMsg << exit(FatalError);
+    }
+    else
+    {
+        // Re-do the track, but this time do the bit tangential to the
+        // direction/patch first. This gets us as close as possible to the
+        // original path/position.
+
+        if (direction == nullptr)
+        {
+            const polyPatch& p = mesh_.boundaryMesh()[patch(facei_)];
+            direction = &p.faceNormals()[patchFace(patch(facei_), facei_)];
+        }
+
+        const vector n = *direction/mag(*direction);
+        const vector s = position - mesh_.cellCentres()[celli_];
+        const vector sN = (s & n)*n;
+        const vector sT = s - sN;
+
+        coordinates_ = barycentric(1, 0, 0, 0);
+        tetFacei_ = mesh_.cells()[celli_][0];
+        tetPti_ = 1;
+        facei_ = -1;
+
+        track(sT, 0);
+        track(sN, 0);
+
+        static label nWarnings = 0;
+        static const label maxNWarnings = 100;
+        if (nWarnings < maxNWarnings)
+        {
+            WarningInFunction << boundaryMsg << endl;
+            ++ nWarnings;
+        }
+        if (nWarnings == maxNWarnings)
+        {
+            WarningInFunction
+                << "Suppressing any further warnings about particles being "
+                << "located outside of the mesh." << endl;
+            ++ nWarnings;
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::particle::particle
 (
     const polyMesh& mesh,
-    const vector& position,
+    const barycentric& coordinates,
     const label celli,
     const label tetFacei,
     const label tetPti
 )
 :
     mesh_(mesh),
-    position_(position),
+    coordinates_(coordinates),
     celli_(celli),
-    facei_(-1),
-    stepFraction_(0.0),
     tetFacei_(tetFacei),
     tetPti_(tetPti),
+    facei_(-1),
+    stepFraction_(0.0),
     origProc_(Pstream::myProcNo()),
     origId_(getNewParticleID())
 {}
@@ -69,36 +598,39 @@ Foam::particle::particle
 (
     const polyMesh& mesh,
     const vector& position,
-    const label celli,
-    bool doCellFacePt
+    const label celli
 )
 :
     mesh_(mesh),
-    position_(position),
+    coordinates_(- VGREAT, - VGREAT, - VGREAT, - VGREAT),
     celli_(celli),
-    facei_(-1),
-    stepFraction_(0.0),
     tetFacei_(-1),
     tetPti_(-1),
+    facei_(-1),
+    stepFraction_(0.0),
     origProc_(Pstream::myProcNo()),
     origId_(getNewParticleID())
 {
-    if (doCellFacePt)
-    {
-        initCellFacePt();
-    }
+    locate
+    (
+        position,
+        nullptr,
+        celli,
+        false,
+        "Particle initialised with a location outside of the mesh."
+    );
 }
 
 
 Foam::particle::particle(const particle& p)
 :
     mesh_(p.mesh_),
-    position_(p.position_),
+    coordinates_(p.coordinates_),
     celli_(p.celli_),
-    facei_(p.facei_),
-    stepFraction_(p.stepFraction_),
     tetFacei_(p.tetFacei_),
     tetPti_(p.tetPti_),
+    facei_(p.facei_),
+    stepFraction_(p.stepFraction_),
     origProc_(p.origProc_),
     origId_(p.origId_)
 {}
@@ -107,12 +639,12 @@ Foam::particle::particle(const particle& p)
 Foam::particle::particle(const particle& p, const polyMesh& mesh)
 :
     mesh_(mesh),
-    position_(p.position_),
+    coordinates_(p.coordinates_),
     celli_(p.celli_),
-    facei_(p.facei_),
-    stepFraction_(p.stepFraction_),
     tetFacei_(p.tetFacei_),
     tetPti_(p.tetPti_),
+    facei_(p.facei_),
+    stepFraction_(p.stepFraction_),
     origProc_(p.origProc_),
     origId_(p.origId_)
 {}
@@ -120,6 +652,319 @@ Foam::particle::particle(const particle& p, const polyMesh& mesh)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::scalar Foam::particle::track
+(
+    const vector& displacement,
+    const scalar fraction
+)
+{
+    scalar f = trackToFace(displacement, fraction);
+
+    while (onInternalFace())
+    {
+        changeCell();
+
+        f *= trackToFace(f*displacement, f*fraction);
+    }
+
+    return f;
+}
+
+
+Foam::scalar Foam::particle::trackToFace
+(
+    const vector& displacement,
+    const scalar fraction
+)
+{
+    scalar f = 1;
+
+    label tetTriI = onFace() ? 0 : -1;
+
+    facei_ = -1;
+
+    while (true)
+    {
+        if (mesh_.moving())
+        {
+            f *= trackToMovingTri(f*displacement, f*fraction, tetTriI);
+        }
+        else
+        {
+            f *= trackToStationaryTri(f*displacement, f*fraction, tetTriI);
+        }
+
+        if (tetTriI == -1)
+        {
+            // The track has completed within the current tet
+            return 0;
+        }
+        else if (tetTriI == 0)
+        {
+            // The track has hit a face, so set the current face and return
+            facei_ = tetFacei_;
+            return f;
+        }
+        else
+        {
+            // Move to the next tet and continue the track
+            changeTet(tetTriI);
+        }
+    }
+}
+
+
+Foam::scalar Foam::particle::trackToStationaryTri
+(
+    const vector& displacement,
+    const scalar fraction,
+    label& tetTriI
+)
+{
+    const vector x0 = position();
+    const vector x1 = displacement;
+    const barycentric y0 = coordinates_;
+
+    if (debug)
+    {
+        Info<< "Tracking from " << x0 << " to " << x0 + x1 << endl;
+    }
+
+    // Get the tet geometry
+    vector centre;
+    scalar detA;
+    barycentricTensor T;
+    tetReverseTransform(centre, detA, T);
+
+    if (debug)
+    {
+        vector o, b, v1, v2;
+        tetGeometry(o, b, v1, v2);
+        Info<< "Tet points o=" << o << ", b=" << b
+            << ", v1=" << v1 << ", v2=" << v2 << endl
+            << "Tet determinant = " << detA << endl
+            << "Start local coordinates = " << y0 << endl;
+    }
+
+    // Calculate the local tracking displacement
+    barycentric Tx1(x1 & T);
+
+    if (debug)
+    {
+        Info<< "Local displacement = " << Tx1 << "/" << detA << endl;
+    }
+
+    // Calculate the hit fraction
+    label iH = -1;
+    scalar muH = detA == 0 ? VGREAT : mag(1/detA);
+    for (label i = 0; i < 4; ++ i)
+    {
+        if (std::isnormal(Tx1[i]) && Tx1[i] < 0)
+        {
+            scalar mu = - y0[i]/Tx1[i];
+
+            if (debug)
+            {
+                Info<< "Hit on tet face " << i << " at local coordinate "
+                    << y0 + mu*Tx1 << ", " << mu*detA*100 << "\% of the "
+                    << "way along the track" << endl;
+            }
+
+            if (0 <= mu && mu < muH)
+            {
+                iH = i;
+                muH = mu;
+            }
+        }
+    }
+
+    // Set the new coordinates
+    barycentric yH = y0 + muH*Tx1;
+
+    // Clamp to zero any negative coordinates generated by round-off error
+    for (label i = 0; i < 4; ++ i)
+    {
+        yH.replace(i, i == iH ? 0 : max(0, yH[i]));
+    }
+
+    // Re-normalise if within the tet
+    if (iH == -1)
+    {
+        yH /= cmptSum(yH);
+    }
+
+    // Set the new position and hit index
+    coordinates_ = yH;
+    tetTriI = iH;
+
+    if (debug)
+    {
+        if (iH != -1)
+        {
+            Info<< "Track hit tet face " << iH << " first" << endl;
+        }
+        else
+        {
+            Info<< "Track hit no tet faces" << endl;
+        }
+        Info<< "End local coordinates = " << yH << endl
+            << "End global coordinates = " << position() << endl
+            << "Tracking displacement = " << position() - x0 << endl
+            << muH*detA*100 << "\% of the step from " << stepFraction_ << " to "
+            << stepFraction_ + fraction << " completed" << endl << endl;
+    }
+
+    // Set the proportion of the track that has been completed
+    stepFraction_ += fraction*muH*detA;
+
+    return 1 - muH*detA;
+}
+
+
+Foam::scalar Foam::particle::trackToMovingTri
+(
+    const vector& displacement,
+    const scalar fraction,
+    label& tetTriI
+)
+{
+    const vector x0 = position();
+    const vector x1 = displacement;
+    const barycentric y0 = coordinates_;
+
+    if (debug)
+    {
+        Info<< "Tracking from " << x0 << " to " << x0 + x1 << endl;
+    }
+
+    // Get the tet geometry
+    Pair<vector> centre;
+    FixedList<scalar, 4> detA;
+    FixedList<barycentricTensor, 3> T;
+    movingTetReverseTransform(fraction, centre, detA, T);
+
+    const vector x0Rel = x0 - centre[0];
+    const vector x1Rel = x1 - centre[1];
+
+    // Form the determinant and hit equations
+    cubicEqn detAEqn(sqr(detA[0])*detA[3], detA[0]*detA[2], detA[1], 1);
+    const barycentric yC(1, 0, 0, 0);
+    const barycentric hitEqnA = ((x1Rel & T[2]) + detA[3]*yC)*sqr(detA[0]);
+    const barycentric hitEqnB =
+        ((x1Rel & T[1]) + (x0Rel & T[2]) + detA[2]*yC)*detA[0];
+    const barycentric hitEqnC = (x1Rel & T[0]) + (x0Rel & T[1]) + detA[1]*yC;
+    const barycentric hitEqnD = y0;
+    FixedList<cubicEqn, 4> hitEqn;
+    forAll(hitEqn, i)
+    {
+        hitEqn[i] = cubicEqn(hitEqnA[i], hitEqnB[i], hitEqnC[i], hitEqnD[i]);
+    }
+
+    // Calculate the hit fraction
+    label iH = -1;
+    scalar muH = detA[0] == 0 ? VGREAT : mag(1/detA[0]);
+    for (label i = 0; i < 4; ++ i)
+    {
+        const Roots<3> mu = hitEqn[i].roots();
+
+        for (label j = 0; j < 3; ++ j)
+        {
+            if (mu.type(j) == roots::real && hitEqn[i].derivative(mu[j]) < 0)
+            {
+                if (0 <= mu[j] && mu[j] < muH)
+                {
+                    iH = i;
+                    muH = mu[j];
+                }
+            }
+        }
+    }
+
+    // Set the new coordinates
+    barycentric yH
+    (
+        hitEqn[0].value(muH),
+        hitEqn[1].value(muH),
+        hitEqn[2].value(muH),
+        hitEqn[3].value(muH)
+    );
+    // !!! <-- This fails if the tet collapses onto the particle, as detA tends
+    // to zero at the hit. In this instance, we can differentiate the hit and
+    // detA polynomials to find a limiting location, but this will not be on a
+    // triangle. We will then need to do a second track through the degenerate
+    // tet to find the final hit position. This second track is over zero
+    // distance and therefore can be of the static mesh type. This has not yet
+    // been implemented.
+    const scalar detAH = detAEqn.value(muH);
+    if (!std::isnormal(detAH))
+    {
+        FatalErrorInFunction
+            << "A moving tet collapsed onto a particle. This is not supported. "
+            << "The mesh is too poor, or the motion too severe, for particle "
+            << "tracking to function." << exit(FatalError);
+    }
+    yH /= detAH;
+
+    // Clamp to zero any negative coordinates generated by round-off error
+    for (label i = 0; i < 4; ++ i)
+    {
+        yH.replace(i, i == iH ? 0 : max(0, yH[i]));
+    }
+
+    // Re-normalise if within the tet
+    if (iH == -1)
+    {
+        yH /= cmptSum(yH);
+    }
+
+    // Set the new position and hit index
+    coordinates_ = yH;
+    tetTriI = iH;
+
+    // Set the proportion of the track that has been completed
+    stepFraction_ += fraction*muH*detA[0];
+
+    if (debug)
+    {
+        if (iH != -1)
+        {
+            Info<< "Track hit tet face " << iH << " first" << endl;
+        }
+        else
+        {
+            Info<< "Track hit no tet faces" << endl;
+        }
+        Info<< "End local coordinates = " << yH << endl
+            << "End global coordinates = " << position() << endl;
+    }
+
+    return 1 - muH*detA[0];
+}
+
+
+void Foam::particle::constrainToMeshCentre()
+{
+    const Vector<label>& dirs = mesh_.geometricD();
+
+    bool isConstrained = false;
+    forAll(dirs, dirI)
+    {
+        if (dirs[dirI] == -1)
+        {
+            isConstrained = true;
+            break;
+        }
+    }
+
+    if (isConstrained)
+    {
+        vector pos = position();
+        meshTools::constrainToMeshCentre(mesh_, pos);
+        track(pos - position(), 0);
+    }
+}
+
+
 void Foam::particle::transformProperties(const tensor&)
 {}
 
@@ -134,6 +979,107 @@ Foam::scalar Foam::particle::wallImpactDistance(const vector&) const
 }
 
 
+void Foam::particle::prepareForInteractionListReferral
+(
+    const vectorTensorTransform& transform
+)
+{
+    // Get the transformed position
+    const vector pos = transform.invTransformPosition(position());
+
+    // Break the topology
+    celli_ = -1;
+    tetFacei_ = -1;
+    tetPti_ = -1;
+    facei_ = -1;
+
+    // Store the position in the barycentric data
+    coordinates_ = barycentric(1 - cmptSum(pos), pos.x(), pos.y(), pos.z());
+
+    // Transform the properties
+    transformProperties(- transform.t());
+    if (transform.hasR())
+    {
+        transformProperties(transform.R().T());
+    }
+}
+
+
+void Foam::particle::correctAfterInteractionListReferral(const label celli)
+{
+    // Get the position from the barycentric data
+    const vector pos(coordinates_.b(), coordinates_.c(), coordinates_.d());
+
+    // Create some arbitrary topology for the supplied cell
+    celli_ = celli;
+    tetFacei_ = mesh_.cells()[celli_][0];
+    tetPti_ = 1;
+    facei_ = -1;
+
+    // Get the reverse transform and directly set the coordinates from the
+    // position. This isn't likely to be correct; the particle is probably not
+    // in this tet. It will, however, generate the correct vector when the
+    // position method is called. A referred particle should never be tracked,
+    // so this approximate topology is good enough. By using the nearby cell we
+    // minimise the error associated with the incorrect topology.
+    if (mesh_.moving())
+    {
+        NotImplemented;
+    }
+    else
+    {
+        vector centre;
+        scalar detA;
+        barycentricTensor T;
+        tetReverseTransform(centre, detA, T);
+        coordinates_ = barycentric(1, 0, 0, 0) + ((pos - centre) & T/detA);
+    }
+}
+
+
+Foam::label Foam::particle::procTetPt
+(
+    const polyMesh& procMesh,
+    const label procCell,
+    const label procTetFace
+) const
+{
+    // The tet point on the procMesh differs from the current tet point if the
+    // mesh and procMesh faces are of differing orientation. The change is the
+    // same as in particle::correctAfterParallelTransfer.
+
+    if
+    (
+        (mesh_.faceOwner()[tetFacei_] == celli_)
+     == (procMesh.faceOwner()[procTetFace] == procCell)
+    )
+    {
+        return tetPti_;
+    }
+    else
+    {
+        return procMesh.faces()[procTetFace].size() - 1 - tetPti_;
+    }
+}
+
+
+void Foam::particle::autoMap
+(
+    const vector& position,
+    const mapPolyMesh& mapper
+)
+{
+    locate
+    (
+        position,
+        nullptr,
+        mapper.reverseCellMap()[celli_],
+        true,
+        "Particle mapped to a location outside of the mesh."
+    );
+}
+
+
 // * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
 
 bool Foam::operator==(const particle& pA, const particle& pB)
diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H
index d94db0dd56..ea678968ce 100644
--- a/src/lagrangian/basic/particle/particle.H
+++ b/src/lagrangian/basic/particle/particle.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -33,6 +33,8 @@ Description
 #define particle_H
 
 #include "vector.H"
+#include "barycentric.H"
+#include "barycentricTensor.H"
 #include "Cloud.H"
 #include "IDLList.H"
 #include "pointField.H"
@@ -42,6 +44,7 @@ Description
 #include "FixedList.H"
 #include "polyMeshTetDecomposition.H"
 #include "particleMacros.H"
+#include "vectorTensorTransform.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -130,25 +133,19 @@ public:
     };
 
 
-protected:
+private:
 
-    // Protected data
+    // Private data
 
         //- Reference to the polyMesh database
         const polyMesh& mesh_;
 
-        //- Position of particle
-        vector position_;
+        //- Coordinates of particle
+        barycentric coordinates_;
 
         //- Index of the cell it is in
         label celli_;
 
-        //- Face index if the particle is on a face otherwise -1
-        label facei_;
-
-        //- Fraction of time-step completed
-        scalar stepFraction_;
-
         //- Index of the face that owns the decomposed tet that the
         //  particle is in
         label tetFacei_;
@@ -158,6 +155,12 @@ protected:
         //  point.
         label tetPti_;
 
+        //- Face index if the particle is on a face otherwise -1
+        label facei_;
+
+        //- Fraction of time-step completed
+        scalar stepFraction_;
+
         //- Originating processor id
         label origProc_;
 
@@ -165,74 +168,137 @@ protected:
         label origId_;
 
 
+private:
+
     // Private Member Functions
 
-        //- Find the tet tri faces between position and tet centre
-        void findTris
-        (
-            const vector& position,
-            DynamicList<label>& faceList,
-            const tetPointRef& tet,
-            const FixedList<vector, 4>& tetAreas,
-            const FixedList<label, 4>& tetPlaneBasePtIs,
-            const scalar tol
-        ) const;
+        // Tetrahedra functions
 
-        //- Find the lambda value for the line to-from across the
-        //  given tri face, where p = from + lambda*(to - from)
-        inline scalar tetLambda
-        (
-            const vector& from,
-            const vector& to,
-            const label triI,
-            const vector& tetArea,
-            const label tetPlaneBasePtI,
-            const label celli,
-            const label tetFacei,
-            const label tetPti,
-            const scalar tol
-        ) const;
+            //- Get indices into the current face for the face-bound vertices of
+            //  the current tet.
+            void tetFaceIndices
+            (
+                label& baseI,
+                label& vertex1I,
+                label& vertex2I
+            ) const;
 
-        //- Find the lambda value for a moving tri face
-        inline scalar movingTetLambda
-        (
-            const vector& from,
-            const vector& to,
-            const label triI,
-            const vector& tetArea,
-            const label tetPlaneBasePtI,
-            const label celli,
-            const label tetFacei,
-            const label tetPti,
-            const scalar tol
-        ) const;
+            //- Get indices into the mesh points for the face-bound vertices of
+            //  the current tet.
+            void tetMeshIndices
+            (
+                label& baseI,
+                label& vertex1I,
+                label& vertex2I
+            ) const;
 
-        //- Modify the tet owner data by crossing triI
-        inline void tetNeighbour(label triI);
+            //- Get the vertices of the current tet
+            void tetGeometry
+            (
+                vector& centre,
+                vector& base,
+                vector& vertex1,
+                vector& vertex2
+            ) const;
 
-        //- Cross the from the given face across the given edge of the
-        //  given cell to find the resulting face and tetPti
-        inline void crossEdgeConnectedFace
-        (
-            const label& celli,
-            label& tetFacei,
-            label& tetPti,
-            const edge& e
-        );
+            //- Get the transformation associated with the current tet. This
+            //  will convert a barycentric position within the tet to a
+            //  cartesian position in the global coordinate system. The
+            //  conversion is x = A & y, where x is the cartesian position, y is
+            //  the barycentric position and A is the transformation tensor.
+            barycentricTensor tetTransform() const;
 
-        //- Hit wall faces in the current cell if the
-        //- wallImpactDistance is non-zero.  They may not be in
-        //- Different tets to the current.
-        template<class CloudType>
-        inline void hitWallFaces
-        (
-            const CloudType& td,
-            const vector& from,
-            const vector& to,
-            scalar& lambdaMin,
-            tetIndices& closestTetIs
-        );
+            //- Get the reverse transform associated with the current tet. The
+            //  conversion is detA*y = (x - centre) & T. The variables x, y and
+            //  centre have the same meaning as for the forward transform. T is
+            //  the transposed inverse of the forward transform tensor, A,
+            //  multiplied by its determinant, detA. This separation allows
+            //  the barycentric tracking algorithm to function on inverted or
+            //  degenerate tetrahedra.
+            void tetReverseTransform
+            (
+                vector& centre,
+                scalar& detA,
+                barycentricTensor& T
+            ) const;
 
+            //- Get the vertices of the current moving tet. Two values are
+            //  returned for each vertex. The first is a constant, and the
+            //  second is a linear coefficient of the track fraction.
+            void movingTetGeometry
+            (
+                const scalar endStepFraction,
+                Pair<vector>& centre,
+                Pair<vector>& base,
+                Pair<vector>& vertex1,
+                Pair<vector>& vertex2
+            ) const;
+
+            //- Get the transformation associated with the current, moving, tet.
+            //  This is of the same form as for the static case. As with the
+            //  moving geometry, a linear function of the tracking fraction is
+            //  returned for each component.
+            Pair<barycentricTensor> movingTetTransform
+            (
+                const scalar endStepFraction
+            ) const;
+
+            //- Get the reverse transformation associated with the current,
+            //  moving, tet. This is of the same form as for the static case. As
+            //  with the moving geometry, a function of the tracking fraction is
+            //  returned for each component. The functions are higher order than
+            //  for the forward transform; the determinant is cubic, and the
+            //  tensor is quadratic.
+            void movingTetReverseTransform
+            (
+                const scalar endStepFraction,
+                Pair<vector>& centre,
+                FixedList<scalar, 4>& detA,
+                FixedList<barycentricTensor, 3>& T
+            ) const;
+
+
+        // Transformations
+
+            //- Reflection transform. Corrects the coordinates when the particle
+            //  moves between two tets which share a base vertex, but for which
+            //  the other two non cell-centre vertices are reversed. All hits
+            //  which retain the same face behave this way, as do face hits.
+            void reflect();
+
+            //- Rotation transform. Corrects the coordinates when the particle
+            //  moves between two tets with different base vertices, but are
+            //  otherwise similarly oriented. Hits which change the face within
+            //  the cell make use of both this and the reflect transform.
+            void rotate(const bool direction);
+
+
+        // Topology changes
+
+            //- Change tet within a cell. Called after a triangle is hit.
+            void changeTet(const label tetTriI);
+
+            //- Change tet face within a cell. Called by changeTet.
+            void changeFace(const label tetTriI);
+
+            //- Change cell. Called when the particle hits an internal face.
+            void changeCell();
+
+
+        // Geometry changes
+
+            //- Locate the particle at the given position
+            void locate
+            (
+                const vector& position,
+                const vector* direction,
+                const label celli,
+                const bool boundaryFail,
+                const string boundaryMsg
+            );
+
+
+protected:
 
     // Patch interactions
 
@@ -315,8 +381,8 @@ public:
         //- String representation of properties
         DefinePropertyList
         (
-            "(Px Py Pz) celli facei stepFraction "
-            "tetFacei tetPti origProc origId"
+            "(Px Py Pz) celli tetFacei tetPti "
+            "facei stepFraction origProc origId"
         );
 
         //- String representation of property types
@@ -328,17 +394,6 @@ public:
         //- Cumulative particle counter - used to provode unique ID
         static label particleCount_;
 
-        //- Fraction of distance to tet centre to move a particle to
-        // 'rescue' it from a tracking problem
-        static const scalar trackingCorrectionTol;
-
-        //- Fraction of the cell volume to use in determining tolerance values
-        //  for the denominator and numerator of lambda
-        static const scalar lambdaDistanceToleranceCoeff;
-
-        //- Minimum stepFraction tolerance
-        static const scalar minStepFractionTol;
-
 
     // Constructors
 
@@ -346,20 +401,19 @@ public:
         particle
         (
             const polyMesh& mesh,
-            const vector& position,
+            const barycentric& coordinates,
             const label celli,
             const label tetFacei,
             const label tetPti
         );
 
-        //- Construct from components, tetFacei_ and tetPti_ are not
-        //  supplied so they will be deduced by a search
+        //- Construct from a position and a cell, searching for the rest of the
+        //  required topology
         particle
         (
             const polyMesh& mesh,
             const vector& position,
-            const label celli,
-            bool doCellFacePt = true
+            const label celli
         );
 
         //- Construct from Istream
@@ -412,27 +466,15 @@ public:
             //- Return the mesh database
             inline const polyMesh& mesh() const;
 
-            //- Return current particle position
-            inline const vector& position() const;
-
-            //- Return current particle position
-            inline vector& position();
-
-            //- Return current cell particle is in
-            inline label& cell();
+            //- Return current particle coordinates
+            inline const barycentric& coordinates() const;
 
             //- Return current cell particle is in
             inline label cell() const;
 
-            //- Return current tet face particle is in
-            inline label& tetFace();
-
             //- Return current tet face particle is in
             inline label tetFace() const;
 
-            //- Return current tet face particle is in
-            inline label& tetPt();
-
             //- Return current tet face particle is in
             inline label tetPt() const;
 
@@ -453,43 +495,32 @@ public:
             //  on oldPoints
             inline vector oldNormal() const;
 
-            //- Return current face particle is on otherwise -1
-            inline label& face();
-
             //- Return current face particle is on otherwise -1
             inline label face() const;
 
-            //- Return the impact model to be used, soft or hard (default).
-            inline bool softImpact() const;
-
             //- Return the particle current time
             inline scalar currentTime() const;
 
 
         // Check
 
-            //- Check the stored cell value (setting if necessary) and
-            //  initialise the tetFace and tetPt values
-            inline void initCellFacePt();
+            //- Is the particle on a face?
+            inline bool onFace() const;
 
-            //- Is the particle on the boundary/(or outside the domain)?
-            inline bool onBoundary() const;
+            //- Is the particle on an internal face?
+            inline bool onInternalFace() const;
 
-            //- Is this global face an internal face?
-            inline bool internalFace(const label facei) const;
+            //- Is the particle on a boundary face?
+            inline bool onBoundaryFace() const;
 
-            //- Is this global face a boundary face?
-            inline bool boundaryFace(const label facei) const;
-
-            //- Which patch is particle on
+            //- Which patch is particle on // <-- !!!
             inline label patch(const label facei) const;
 
-            //- Which face of this patch is this particle on
-            inline label patchFace
-            (
-                const label patchi,
-                const label facei
-            ) const;
+            //- Which face of this patch is this particle on // <-- !!!
+            inline label patchFace(const label patchi, const label facei) const;
+
+            //- Return current particle position
+            inline vector position() const;
 
             //- Return the fraction of time-step completed
             inline scalar& stepFraction();
@@ -512,28 +543,57 @@ public:
 
         // Track
 
-            //- Track particle to end of trajectory
-            //  or until it hits the boundary.
-            //  On entry 'stepFraction()' should be set to the fraction of the
-            //  time-step at which the tracking starts and on exit it contains
-            //  the fraction of the time-step completed.
-            //  Returns the boundary face index if the track stops at the
-            //  boundary, -1 otherwise.
-            template<class TrackData>
-            label track(const vector& endPosition, TrackData& td);
+            //- Track along the displacement for a given fraction of the overall
+            //  step. End when the track is complete, or when a boundary is hit.
+            //  On exit, stepFraction_ will have been incremented to the current
+            //  position, and facei_ will be set to the index of the boundary
+            //  face that was hit, or -1 if the track completed within a cell.
+            //  The proportion of the displacement still to be completed is
+            //  returned.
+            scalar track
+            (
+                const vector& displacement,
+                const scalar fraction
+            );
 
-            //- Track particle to a given position and returns 1.0 if the
-            //  trajectory is completed without hitting a face otherwise
-            //  stops at the face and returns the fraction of the trajectory
-            //  completed.
-            //  on entry 'stepFraction()' should be set to the fraction of the
-            //  time-step at which the tracking starts.
-            template<class TrackData>
-            scalar trackToFace(const vector& endPosition, TrackData& td);
+            //- As particle::track, but also stops on internal faces.
+            scalar trackToFace
+            (
+                const vector& displacement,
+                const scalar fraction
+            );
 
-            //- Return the index of the face to be used in the interpolation
-            //  routine
-            inline label faceInterpolation() const;
+            //- As particle::trackToFace, but also stops on tet triangles. On
+            //  exit, tetTriI is set to the index of the tet triangle that was
+            //  hit, or -1 if the end position was reached.
+            scalar trackToStationaryTri
+            (
+                const vector& displacement,
+                const scalar fraction,
+                label& tetTriI
+            );
+
+            //- As particle::trackToTri, but for moving meshes
+            scalar trackToMovingTri
+            (
+                const vector& displacement,
+                const scalar fraction,
+                label& tetTriI
+            );
+
+            //- As non-templated particle::trackToFace, but with additional
+            //  boundary handling.
+            template<class TrackData>
+            void trackToFace
+            (
+                const vector& displacement,
+                const scalar fraction,
+                TrackData& td
+            );
+
+            //- Set the constrained components of the particle position to the
+            //  mesh centre.
+            void constrainToMeshCentre();
 
 
     // Transformations
@@ -564,6 +624,38 @@ public:
         void correctAfterParallelTransfer(const label patchi, TrackData& td);
 
 
+    // Interaction list referral
+
+        //- Break the topology and store the particle position so that the
+        //  particle can be referred.
+        void prepareForInteractionListReferral
+        (
+            const vectorTensorTransform& transform
+        );
+
+        //- Correct the topology after referral. The particle may still be
+        //  outside the stored tet and therefore not track-able.
+        void correctAfterInteractionListReferral(const label celli);
+
+
+    // Decompose and reconstruct
+
+        //- Return the tet point approproate for decomposition or reconstruction
+        //  to or from the given mesh.
+        label procTetPt
+        (
+            const polyMesh& procMesh,
+            const label procCell,
+            const label procTetFace
+        ) const;
+
+
+    // Mapping
+
+        //- Map after a topology change
+        void autoMap(const vector& position, const mapPolyMesh& mapper);
+
+
     // I-O
 
         //- Read the fields associated with the owner cloud
diff --git a/src/lagrangian/basic/particle/particleI.H b/src/lagrangian/basic/particle/particleI.H
index 5df1a22e92..074bd8a5b0 100644
--- a/src/lagrangian/basic/particle/particleI.H
+++ b/src/lagrangian/basic/particle/particleI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,541 +26,6 @@ License
 #include "polyMesh.H"
 #include "Time.H"
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-inline void Foam::particle::findTris
-(
-    const vector& position,
-    DynamicList<label>& faceList,
-    const tetPointRef& tet,
-    const FixedList<vector, 4>& tetAreas,
-    const FixedList<label, 4>& tetPlaneBasePtIs,
-    const scalar tol
-) const
-{
-    faceList.clear();
-
-    const point Ct = tet.centre();
-
-    for (label i = 0; i < 4; i++)
-    {
-        scalar lambda = tetLambda
-        (
-            Ct,
-            position,
-            i,
-            tetAreas[i],
-            tetPlaneBasePtIs[i],
-            celli_,
-            tetFacei_,
-            tetPti_,
-            tol
-        );
-
-        if ((lambda > 0.0) && (lambda < 1.0))
-        {
-            faceList.append(i);
-        }
-    }
-}
-
-
-inline Foam::scalar Foam::particle::tetLambda
-(
-    const vector& from,
-    const vector& to,
-    const label triI,
-    const vector& n,
-    const label tetPlaneBasePtI,
-    const label celli,
-    const label tetFacei,
-    const label tetPti,
-    const scalar tol
-) const
-{
-    const pointField& pPts = mesh_.points();
-
-    if (mesh_.moving())
-    {
-        return movingTetLambda
-        (
-            from,
-            to,
-            triI,
-            n,
-            tetPlaneBasePtI,
-            celli,
-            tetFacei,
-            tetPti,
-            tol
-        );
-    }
-
-    const point& base = pPts[tetPlaneBasePtI];
-
-    scalar lambdaNumerator = (base - from) & n;
-    scalar lambdaDenominator = (to - from) & n;
-
-    // n carries the area of the tet faces, so the dot product with a
-    // delta-length has the units of volume.  Comparing the component of each
-    // delta-length in the direction of n times the face area to a fraction of
-    // the cell volume.
-
-    if (mag(lambdaDenominator) < tol)
-    {
-        if (mag(lambdaNumerator) < tol)
-        {
-            // Track starts on the face, and is potentially
-            // parallel to it.  +-tol/+-tol is not a good
-            // comparison, return 0.0, in anticipation of tet
-            // centre correction.
-            return 0.0;
-        }
-        else
-        {
-            if (mag((to - from)) < tol/mag(n))
-            {
-                // 'Zero' length track (compared to the tolerance, which is
-                // based on the cell volume, divided by the tet face area), not
-                // along the face, face cannot be crossed.
-                return GREAT;
-            }
-            else
-            {
-                // Trajectory is non-zero and parallel to face
-                lambdaDenominator = sign(lambdaDenominator)*SMALL;
-            }
-        }
-    }
-
-    return lambdaNumerator/lambdaDenominator;
-}
-
-
-inline Foam::scalar Foam::particle::movingTetLambda
-(
-    const vector& from,
-    const vector& to,
-    const label triI,
-    const vector& n,
-    const label tetPlaneBasePtI,
-    const label celli,
-    const label tetFacei,
-    const label tetPti,
-    const scalar tol
-) const
-{
-    const pointField& pPts = mesh_.points();
-    const pointField& oldPPts = mesh_.oldPoints();
-
-    // Base point of plane at end of motion
-    const point& b = pPts[tetPlaneBasePtI];
-
-    // n: Normal of plane at end of motion
-
-    // Base point of plane at start of timestep
-    const point& b00 = oldPPts[tetPlaneBasePtI];
-
-    // Base point of plane at start of tracking portion (cast forward by
-    // stepFraction)
-    point b0 = b00 + stepFraction_*(b - b00);
-
-    // Normal of plane at start of tracking portion
-    vector n0 = Zero;
-
-    {
-        tetIndices tetIs(celli, tetFacei, tetPti, mesh_);
-
-        // Tet at timestep start
-        tetPointRef tet00 = tetIs.oldTet(mesh_);
-
-        // Tet at timestep end
-        tetPointRef tet = tetIs.tet(mesh_);
-
-        point tet0PtA = tet00.a() + stepFraction_*(tet.a() - tet00.a());
-        point tet0PtB = tet00.b() + stepFraction_*(tet.b() - tet00.b());
-        point tet0PtC = tet00.c() + stepFraction_*(tet.c() - tet00.c());
-        point tet0PtD = tet00.d() + stepFraction_*(tet.d() - tet00.d());
-
-        // Tracking portion start tet (cast forward by stepFraction)
-        tetPointRef tet0(tet0PtA, tet0PtB, tet0PtC, tet0PtD);
-
-        switch (triI)
-        {
-            case 0:
-            {
-                n0 = tet0.Sa();
-                break;
-            }
-            case 1:
-            {
-                n0 = tet0.Sb();
-                break;
-            }
-            case 2:
-            {
-                n0 = tet0.Sc();
-                break;
-            }
-            case 3:
-            {
-                n0 = tet0.Sd();
-                break;
-            }
-            default:
-            {
-                break;
-            }
-        }
-    }
-
-    if (mag(n0) < SMALL)
-    {
-        // If the old normal is zero (for example in layer addition)
-        // then use the current normal;
-        n0 = n;
-    }
-
-    scalar lambdaNumerator = 0;
-    scalar lambdaDenominator = 0;
-
-    vector dP = to - from;
-    vector dN = n - n0;
-    vector dB = b - b0;
-    vector dS = from - b0;
-
-    if (mag(dN) > SMALL)
-    {
-        scalar a = (dP - dB) & dN;
-        scalar b = ((dP - dB) & n0) + (dS & dN);
-        scalar c = dS & n0;
-
-        if (mag(a) > SMALL)
-        {
-
-            // Solve quadratic for lambda
-            scalar discriminant = sqr(b) - 4.0*a*c;
-
-            if (discriminant < 0)
-            {
-                // Imaginary roots only - face not crossed
-                return GREAT;
-            }
-            else
-            {
-                scalar q = -0.5*(b + sign(b)*Foam::sqrt(discriminant));
-
-                if (mag(q) < VSMALL)
-                {
-                    // If q is zero, then l1 = q/a is the required
-                    // value of lambda, and is zero.
-                    return 0.0;
-                }
-
-                scalar l1 = q/a;
-                scalar l2 = c/q;
-
-                // There will be two roots, a big one and a little
-                // one, choose the little one.
-
-                if (mag(l1) < mag(l2))
-                {
-                    return l1;
-                }
-                else
-                {
-                    return l2;
-                }
-            }
-        }
-        {
-            // When a is zero, solve the first order polynomial
-            lambdaNumerator = -c;
-            lambdaDenominator = b;
-        }
-    }
-    else
-    {
-        // When n = n0 is zero, there is no plane rotation, solve the
-        // first order polynomial
-        lambdaNumerator = -(dS & n0);
-        lambdaDenominator = ((dP - dB) & n0);
-    }
-
-    if (mag(lambdaDenominator) < tol)
-    {
-        if (mag(lambdaNumerator) < tol)
-        {
-            // Track starts on the face, and is potentially
-            // parallel to it.  +-tol)/+-tol is not a good
-            // comparison, return 0.0, in anticipation of tet
-            // centre correction.
-            return 0.0;
-        }
-        else
-        {
-            if (mag((to - from)) < tol/mag(n))
-            {
-                // Zero length track, not along the face, face
-                // cannot be crossed.
-                return GREAT;
-            }
-            else
-            {
-                // Trajectory is non-zero and parallel to face
-                lambdaDenominator = sign(lambdaDenominator)*SMALL;
-            }
-        }
-    }
-
-    return lambdaNumerator/lambdaDenominator;
-}
-
-
-
-inline void Foam::particle::tetNeighbour(label triI)
-{
-    const labelList& pOwner = mesh_.faceOwner();
-    const faceList& pFaces = mesh_.faces();
-
-    bool own = (pOwner[tetFacei_] == celli_);
-
-    const Foam::face& f = pFaces[tetFacei_];
-
-    label tetBasePtI = mesh_.tetBasePtIs()[tetFacei_];
-
-    if (tetBasePtI == -1)
-    {
-        FatalErrorInFunction
-            << "No base point for face " << tetFacei_ << ", " << f
-            << ", produces a valid tet decomposition."
-            << abort(FatalError);
-    }
-
-    label facePtI = (tetPti_ + tetBasePtI) % f.size();
-    label otherFacePtI = f.fcIndex(facePtI);
-
-    switch (triI)
-    {
-        case 0:
-        {
-            // Crossing this triangle changes tet to that in the
-            // neighbour cell over tetFacei
-
-            // Modification of celli_ will happen by other indexing,
-            // tetFacei_ and tetPti don't change.
-
-            break;
-        }
-        case 1:
-        {
-            crossEdgeConnectedFace
-            (
-                celli_,
-                tetFacei_,
-                tetPti_,
-                Foam::edge(f[facePtI], f[otherFacePtI])
-            );
-
-            break;
-        }
-        case 2:
-        {
-            if (own)
-            {
-                if (tetPti_ < f.size() - 2)
-                {
-                    tetPti_ = f.fcIndex(tetPti_);
-                }
-                else
-                {
-                    crossEdgeConnectedFace
-                    (
-                        celli_,
-                        tetFacei_,
-                        tetPti_,
-                        Foam::edge(f[tetBasePtI], f[otherFacePtI])
-                    );
-                }
-            }
-            else
-            {
-                if (tetPti_ > 1)
-                {
-                    tetPti_ = f.rcIndex(tetPti_);
-                }
-                else
-                {
-                    crossEdgeConnectedFace
-                    (
-                        celli_,
-                        tetFacei_,
-                        tetPti_,
-                        Foam::edge(f[tetBasePtI], f[facePtI])
-                    );
-                }
-            }
-
-            break;
-        }
-        case 3:
-        {
-            if (own)
-            {
-                if (tetPti_ > 1)
-                {
-                    tetPti_ = f.rcIndex(tetPti_);
-                }
-                else
-                {
-                    crossEdgeConnectedFace
-                    (
-                        celli_,
-                        tetFacei_,
-                        tetPti_,
-                        Foam::edge(f[tetBasePtI], f[facePtI])
-                    );
-                }
-            }
-            else
-            {
-                if (tetPti_ < f.size() - 2)
-                {
-                    tetPti_ = f.fcIndex(tetPti_);
-                }
-                else
-                {
-                    crossEdgeConnectedFace
-                    (
-                        celli_,
-                        tetFacei_,
-                        tetPti_,
-                        Foam::edge(f[tetBasePtI], f[otherFacePtI])
-                    );
-                }
-            }
-
-            break;
-        }
-        default:
-        {
-            FatalErrorInFunction
-                << "Tet tri face index error, can only be 0..3, supplied "
-                << triI << abort(FatalError);
-
-            break;
-        }
-    }
-}
-
-
-inline void Foam::particle::crossEdgeConnectedFace
-(
-    const label& celli,
-    label& tetFacei,
-    label& tetPti,
-    const edge& e
-)
-{
-    const faceList& pFaces = mesh_.faces();
-    const cellList& pCells = mesh_.cells();
-
-    const Foam::face& f = pFaces[tetFacei];
-
-    const Foam::cell& thisCell = pCells[celli];
-
-    forAll(thisCell, cFI)
-    {
-        // Loop over all other faces of this cell and
-        // find the one that shares this edge
-
-        label fI = thisCell[cFI];
-
-        if (tetFacei == fI)
-        {
-            continue;
-        }
-
-        const Foam::face& otherFace = pFaces[fI];
-
-        label edDir = otherFace.edgeDirection(e);
-
-        if (edDir == 0)
-        {
-            continue;
-        }
-        else if (f == pFaces[fI])
-        {
-            // This is a necessary condition if using duplicate baffles
-            // (so coincident faces). We need to make sure we don't cross into
-            // the face with the same vertices since we might enter a tracking
-            // loop where it never exits. This test should be cheap
-            // for most meshes so can be left in for 'normal' meshes.
-            continue;
-        }
-        else
-        {
-            //Found edge on other face
-            tetFacei = fI;
-
-            label eIndex = -1;
-
-            if (edDir == 1)
-            {
-                // Edge is in the forward circulation of this face, so
-                // work with the start point of the edge
-                eIndex = findIndex(otherFace, e.start());
-            }
-            else
-            {
-                // edDir == -1, so the edge is in the reverse
-                // circulation of this face, so work with the end
-                // point of the edge
-                eIndex = findIndex(otherFace, e.end());
-            }
-
-            label tetBasePtI = mesh_.tetBasePtIs()[fI];
-
-            if (tetBasePtI == -1)
-            {
-                FatalErrorInFunction
-                    << "No base point for face " << fI << ", " << f
-                    << ", produces a decomposition that has a minimum "
-                    << "volume greater than tolerance."
-                    << abort(FatalError);
-            }
-
-            // Find eIndex relative to the base point on new face
-            eIndex -= tetBasePtI;
-
-            if (neg(eIndex))
-            {
-                eIndex = (eIndex + otherFace.size()) % otherFace.size();
-            }
-
-            if (eIndex == 0)
-            {
-                // The point is the base point, so this is first tet
-                // in the face circulation
-                tetPti = 1;
-            }
-            else if (eIndex == otherFace.size() - 1)
-            {
-                // The point is the last before the base point, so
-                // this is the last tet in the face circulation
-                tetPti = otherFace.size() - 2;
-            }
-            else
-            {
-                tetPti = eIndex;
-            }
-
-            break;
-        }
-    }
-}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 inline Foam::label Foam::particle::getNewParticleID() const
@@ -583,15 +48,9 @@ inline const Foam::polyMesh& Foam::particle::mesh() const
 }
 
 
-inline const Foam::vector& Foam::particle::position() const
+inline const Foam::barycentric& Foam::particle::coordinates() const
 {
-    return position_;
-}
-
-
-inline Foam::vector& Foam::particle::position()
-{
-    return position_;
+    return coordinates_;
 }
 
 
@@ -601,36 +60,18 @@ inline Foam::label Foam::particle::cell() const
 }
 
 
-inline Foam::label& Foam::particle::cell()
-{
-    return celli_;
-}
-
-
 inline Foam::label Foam::particle::tetFace() const
 {
     return tetFacei_;
 }
 
 
-inline Foam::label& Foam::particle::tetFace()
-{
-    return tetFacei_;
-}
-
-
 inline Foam::label Foam::particle::tetPt() const
 {
     return tetPti_;
 }
 
 
-inline Foam::label& Foam::particle::tetPt()
-{
-    return tetPti_;
-}
-
-
 inline Foam::tetIndices Foam::particle::currentTetIndices() const
 {
     return tetIndices(celli_, tetFacei_, tetPti_, mesh_);
@@ -661,161 +102,37 @@ inline Foam::label Foam::particle::face() const
 }
 
 
-inline Foam::label& Foam::particle::face()
+inline bool Foam::particle::onFace() const
 {
-    return facei_;
+    return facei_ >= 0;
 }
 
 
-inline void Foam::particle::initCellFacePt()
+inline bool Foam::particle::onInternalFace() const
 {
-    if (celli_ == -1)
-    {
-        mesh_.findCellFacePt
-        (
-            position_,
-            celli_,
-            tetFacei_,
-            tetPti_
-        );
+    return onFace() && mesh_.isInternalFace(facei_);
+}
 
-        if (debug && celli_ == -1)
-        {
-            WarningInFunction
-                << "cell, tetFace and tetPt search failure for position "
-                << position_
-                << endl;
-        }
+
+inline bool Foam::particle::onBoundaryFace() const
+{
+    return onFace() && !mesh_.isInternalFace(facei_);
+}
+
+
+inline Foam::vector Foam::particle::position() const
+{
+    if (mesh_.moving())
+    {
+        return movingTetTransform(0)[0] & coordinates_;
     }
     else
     {
-        mesh_.findTetFacePt(celli_, position_, tetFacei_, tetPti_);
-
-        if (tetFacei_ == -1 || tetPti_ == -1)
-        {
-            label oldCelli = celli_;
-
-            mesh_.findCellFacePt
-            (
-                position_,
-                celli_,
-                tetFacei_,
-                tetPti_
-            );
-
-            if (celli_ == -1 || tetFacei_ == -1 || tetPti_ == -1)
-            {
-                // The particle has entered this function with a cell number,
-                // but hasn't been able to find a cell to occupy.
-
-                if (!mesh_.pointInCellBB(position_, oldCelli, 0.1))
-                {
-                    // If the position is not inside the (slightly extended)
-                    // bound-box of the cell that it thought it should be in,
-                    // then this is considered an error.
-
-                    if (debug)
-                    {
-                        WarningInFunction
-                            << "position " << position_
-                            << " not in specified cell " << oldCelli
-                            << " with centre " << mesh_.cellCentres()[oldCelli]
-                            << endl;
-                    }
-
-                    celli_ = -1;
-                    tetFacei_ = -1;
-                    tetPti_ = -1;
-
-                    return;
-                }
-
-                // The position is in the (slightly extended) bound-box of the
-                // cell.  This situation may arise because the face
-                // decomposition of the cell is not the same as when the
-                // particle acquired the cell index.  For example, it has been
-                // read into a mesh that has made a different face base-point
-                // decision for a boundary face and now this particle is in a
-                // position that is not in the mesh.  Gradually move the
-                // particle towards the centre of the cell that it thought that
-                // it was in.
-
-                celli_ = oldCelli;
-
-                point newPosition = position_;
-
-                const point& cC = mesh_.cellCentres()[celli_];
-
-                label trap(1.0/trackingCorrectionTol + 1);
-
-                label iterNo = 0;
-
-                do
-                {
-                    newPosition += trackingCorrectionTol*(cC - position_);
-
-                    mesh_.findTetFacePt
-                    (
-                        celli_,
-                        newPosition,
-                        tetFacei_,
-                        tetPti_
-                    );
-
-                    iterNo++;
-
-                } while (tetFacei_ < 0  && iterNo <= trap);
-
-                if (tetFacei_ == -1)
-                {
-                    FatalErrorInFunction
-                        << "cell, tetFace and tetPt search failure at position "
-                        << position_ << abort(FatalError);
-                }
-
-                if (debug)
-                {
-                    WarningInFunction
-                        << "Particle moved from " << position_
-                        << " to " << newPosition
-                        << " in cell " << celli_
-                        << " tetFace " << tetFacei_
-                        << " tetPt " << tetPti_ << nl
-                        << "    (A fraction of "
-                        << 1.0 - mag(cC - newPosition)/mag(cC - position_)
-                        << " of the distance to the cell centre)"
-                        << " because a decomposition tetFace and tetPt "
-                        << "could not be found."
-                        << endl;
-                }
-
-                position_ = newPosition;
-            }
-
-            if (debug && celli_ != oldCelli)
-            {
-                WarningInFunction
-                    << "Particle at position " << position_
-                    << " searched for a cell, tetFace and tetPt." << nl
-                    << "    Found"
-                    << " cell " << celli_
-                    << " tetFace " << tetFacei_
-                    << " tetPt " << tetPti_ << nl
-                    << "    This is a different cell to that which was supplied"
-                    << " (" << oldCelli << ")." << nl
-                    << endl;
-            }
-        }
+        return tetTransform() & coordinates_;
     }
 }
 
 
-inline bool Foam::particle::onBoundary() const
-{
-    return facei_ != -1 && facei_ >= mesh_.nInternalFaces();
-}
-
-
 inline Foam::scalar& Foam::particle::stepFraction()
 {
     return stepFraction_;
@@ -852,29 +169,9 @@ inline Foam::label& Foam::particle::origId()
 }
 
 
-inline bool Foam::particle::softImpact() const
-{
-    return false;
-}
-
-
 inline Foam::scalar Foam::particle::currentTime() const
 {
-    return
-        mesh_.time().value()
-      + stepFraction_*mesh_.time().deltaTValue();
-}
-
-
-inline bool Foam::particle::internalFace(const label facei) const
-{
-    return mesh_.isInternalFace(facei);
-}
-
-
-bool Foam::particle::boundaryFace(const label facei) const
-{
-    return !internalFace(facei);
+    return mesh_.time().value() + stepFraction_*mesh_.time().deltaTValue();
 }
 
 
@@ -894,10 +191,4 @@ inline Foam::label Foam::particle::patchFace
 }
 
 
-inline Foam::label Foam::particle::faceInterpolation() const
-{
-    return facei_;
-}
-
-
 // ************************************************************************* //
diff --git a/src/lagrangian/basic/particle/particleIO.C b/src/lagrangian/basic/particle/particleIO.C
index b9e40bc3d9..e2842be409 100644
--- a/src/lagrangian/basic/particle/particleIO.C
+++ b/src/lagrangian/basic/particle/particleIO.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -33,12 +33,12 @@ Foam::string Foam::particle::propertyTypes_ = Foam::particle::propertyTypes();
 
 const std::size_t Foam::particle::sizeofPosition_
 (
-    offsetof(particle, facei_) - offsetof(particle, position_)
+    offsetof(particle, facei_) - offsetof(particle, coordinates_)
 );
 
 const std::size_t Foam::particle::sizeofFields
 (
-    sizeof(particle) - offsetof(particle, position_)
+    sizeof(particle) - offsetof(particle, coordinates_)
 );
 
 
@@ -47,38 +47,33 @@ const std::size_t Foam::particle::sizeofFields
 Foam::particle::particle(const polyMesh& mesh, Istream& is, bool readFields)
 :
     mesh_(mesh),
-    position_(),
+    coordinates_(),
     celli_(-1),
-    facei_(-1),
-    stepFraction_(0.0),
     tetFacei_(-1),
     tetPti_(-1),
+    facei_(-1),
+    stepFraction_(0.0),
     origProc_(Pstream::myProcNo()),
     origId_(-1)
 {
     if (is.format() == IOstream::ASCII)
     {
-        is  >> position_ >> celli_;
+        is  >> coordinates_ >> celli_ >> tetFacei_ >> tetPti_;
 
         if (readFields)
         {
-            is  >> facei_
-                >> stepFraction_
-                >> tetFacei_
-                >> tetPti_
-                >> origProc_
-                >> origId_;
+            is  >> facei_ >> stepFraction_ >> origProc_ >> origId_;
         }
     }
     else
     {
         if (readFields)
         {
-            is.read(reinterpret_cast<char*>(&position_), sizeofFields);
+            is.read(reinterpret_cast<char*>(&coordinates_), sizeofFields_);
         }
         else
         {
-            is.read(reinterpret_cast<char*>(&position_), sizeofPosition_);
+            is.read(reinterpret_cast<char*>(&coordinates_), sizeofPosition_);
         }
     }
 
@@ -91,11 +86,14 @@ void Foam::particle::writePosition(Ostream& os) const
 {
     if (os.format() == IOstream::ASCII)
     {
-        os  << position_ << token::SPACE << celli_;
+        os  << coordinates_
+            << token::SPACE << celli_
+            << token::SPACE << tetFacei_
+            << token::SPACE << tetPti_;
     }
     else
     {
-        os.write(reinterpret_cast<const char*>(&position_), sizeofPosition_);
+        os.write(reinterpret_cast<const char*>(&coordinates_), sizeofPosition_);
     }
 
     // Check state of Ostream
@@ -107,12 +105,12 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const particle& p)
 {
     if (os.format() == IOstream::ASCII)
     {
-        os  << p.position_
+        os  << p.coordinates_
             << token::SPACE << p.celli_
-            << token::SPACE << p.facei_
-            << token::SPACE << p.stepFraction_
             << token::SPACE << p.tetFacei_
             << token::SPACE << p.tetPti_
+            << token::SPACE << p.facei_
+            << token::SPACE << p.stepFraction_
             << token::SPACE << p.origProc_
             << token::SPACE << p.origId_;
     }
@@ -120,8 +118,8 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const particle& p)
     {
         os.write
         (
-            reinterpret_cast<const char*>(&p.position_),
-            particle::sizeofFields
+            reinterpret_cast<const char*>(&p.coordinates_),
+            particle::sizeofFields_
         );
     }
 
diff --git a/src/lagrangian/basic/particle/particleTemplates.C b/src/lagrangian/basic/particle/particleTemplates.C
index 95f018665c..49e8825736 100644
--- a/src/lagrangian/basic/particle/particleTemplates.C
+++ b/src/lagrangian/basic/particle/particleTemplates.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -44,7 +44,7 @@ void Foam::particle::prepareForParallelTransfer
 )
 {
     // Convert the face index to be local to the processor patch
-    facei_ = patchFace(patchi, facei_);
+    facei_ = mesh_.boundaryMesh()[patchi].whichFace(facei_);
 }
 
 
@@ -58,12 +58,6 @@ void Foam::particle::correctAfterParallelTransfer
     const coupledPolyPatch& ppp =
         refCast<const coupledPolyPatch>(mesh_.boundaryMesh()[patchi]);
 
-    celli_ = ppp.faceCells()[facei_];
-
-    // Have patch transform the position
-    ppp.transformPosition(position_, facei_);
-
-    // Transform the properties
     if (!ppp.parallel())
     {
         const tensor& T =
@@ -85,41 +79,22 @@ void Foam::particle::correctAfterParallelTransfer
         transformProperties(-s);
     }
 
-    tetFacei_ = facei_ + ppp.start();
-
-    // Faces either side of a coupled patch have matched base indices,
-    // tetPti is specified relative to the base point, already and
-    // opposite circulation directions by design, so if the vertices
-    // are:
-    // source:
-    // face    (a b c d e f)
-    // fPtI     0 1 2 3 4 5
-    //            +
-    // destination:
-    // face    (a f e d c b)
-    // fPtI     0 1 2 3 4 5
-    //                  +
-    // where a is the base point of the face are matching , and we
-    // have fPtI = 1 on the source processor face, i.e. vertex b, then
-    // this because of the face circulation direction change, vertex c
-    // is the characterising point on the destination processor face,
-    // giving the destination fPtI as:
-    //     fPtI_d = f.size() - 1 - fPtI_s = 6 - 1 - 1 = 4
-    // This relationship can be verified for other points and sizes of
-    // face.
-
+    // Set the topology
+    celli_ = ppp.faceCells()[facei_];
+    facei_ += ppp.start();
+    tetFacei_ = facei_;
+    // Faces either side of a coupled patch are numbered in opposite directions
+    // as their normals both point away from their connected cells. The tet
+    // point therefore counts in the opposite direction from the base point.
     tetPti_ = mesh_.faces()[tetFacei_].size() - 1 - tetPti_;
 
-    // Reset the face index for the next tracking operation
-    if (stepFraction_ > (1.0 - SMALL))
-    {
-        stepFraction_ = 1.0;
-        facei_ = -1;
-    }
-    else
-    {
-        facei_ += ppp.start();
-    }
+    // Reflect to account for the change of triangle orientation in the new cell
+    reflect();
+
+    // Note that the position does not need transforming explicitly. The face-
+    // triangle on the receive patch is the transformation of the one on the
+    // send patch, so whilst the barycentric coordinates remain the same, the
+    // change of triangle implicitly transforms the position.
 }
 
 
@@ -158,7 +133,6 @@ void Foam::particle::readFields(CloudType& c)
 template<class CloudType>
 void Foam::particle::writeFields(const CloudType& c)
 {
-    // Write the cloud position file
     IOPosition<CloudType> ioP(c);
     ioP.write();
 
@@ -169,7 +143,11 @@ void Foam::particle::writeFields(const CloudType& c)
         c.fieldIOobject("origProcId", IOobject::NO_READ),
         np
     );
-    IOField<label> origId(c.fieldIOobject("origId", IOobject::NO_READ), np);
+    IOField<label> origId
+    (
+        c.fieldIOobject("origId", IOobject::NO_READ),
+        np
+    );
 
     label i = 0;
     forAllConstIter(typename CloudType, c, iter)
@@ -208,387 +186,55 @@ void Foam::particle::writeObjects(const CloudType& c, objectRegistry& obr)
 
 
 template<class TrackData>
-Foam::label Foam::particle::track(const vector& endPosition, TrackData& td)
-{
-    facei_ = -1;
-
-    // Tracks to endPosition or stop on boundary
-    while (!onBoundary() && stepFraction_ < 1.0 - SMALL)
-    {
-        stepFraction_ += trackToFace(endPosition, td)*(1.0 - stepFraction_);
-    }
-
-    return facei_;
-}
-
-
-template<class TrackData>
-Foam::scalar Foam::particle::trackToFace
+void Foam::particle::trackToFace
 (
-    const vector& endPosition,
+    const vector& displacement,
+    const scalar fraction,
     TrackData& td
 )
 {
-    typedef typename TrackData::cloudType cloudType;
-    typedef typename cloudType::particleType particleType;
+    // Track
+    trackToFace(displacement, fraction);
 
-    cloudType& cloud = td.cloud();
-
-
-    const faceList& pFaces = mesh_.faces();
-    const pointField& pPts = mesh_.points();
-    const vectorField& pC = mesh_.cellCentres();
-
-    facei_ = -1;
-
-    // Pout<< "Particle " << origId_ << " " << origProc_
-    //     << " Tracking from " << position_
-    //     << " to " << endPosition
-    //     << endl;
-
-    // Pout<< "stepFraction " << stepFraction_ << nl
-    //     << "celli " << celli_ << nl
-    //     << "tetFacei " << tetFacei_ << nl
-    //     << "tetPti " << tetPti_
-    //     << endl;
-
-    scalar trackFraction = 0.0;
-
-    // Minimum tetrahedron decomposition of each cell of the mesh into
-    // using the cell centre, base point on face, and further two
-    // points on the face.  For each face of n points, there are n - 2
-    // tets generated.
-
-    // The points for each tet are organised to match those used in the
-    // tetrahedron class, supplying them in the order:
-    //     Cc, basePt, pA, pB
-    // where:
-    //   + Cc is the cell centre;
-    //   + basePt is the base point on the face;
-    //   + pA and pB are the remaining points on the face, such that
-    //     the circulation, {basePt, pA, pB} produces a positive
-    //     normal by the right-hand rule.  pA and pB are chosen from
-    //     tetPti_ do accomplish this depending if the cell owns the
-    //     face, tetPti_ is the vertex that characterises the tet, and
-    //     is the first vertex on the tet when circulating around the
-    //     face. Therefore, the same tetPti represents the same face
-    //     triangle for both the owner and neighbour cell.
-    //
-    // Each tet has its four triangles represented in the same order:
-    // 0) tri joining a tet to the tet across the face in next cell.
-    //    This is the triangle opposite Cc.
-    // 1) tri joining a tet to the tet that is in the same cell, but
-    //    belongs to the face that shares the edge of the current face
-    //    that doesn't contain basePt.  This is the triangle opposite
-    //    basePt.
-
-    // 2) tri joining a tet to the tet that is in the same cell, but
-    //    belongs to the face that shares the tet-edge (basePt - pB).
-    //    This may be on the same face, or a different one.  This is
-    //    the triangle opposite basePt.  This is the triangle opposite
-    //    pA.
-
-    // 4) tri joining a tet to the tet that is in the same cell, but
-    //    belongs to the face that shares the tet-edge (basePt - pA).
-    //    This may be on the same face, or a different one.  This is
-    //    the triangle opposite basePt.  This is the triangle opposite
-    //    pA.
-
-    // Which tri (0..3) of the tet has been crossed
-    label triI = -1;
-
-    // Determine which face was actually crossed.  lambdaMin < SMALL
-    // is considered a trigger for a tracking correction towards the
-    // current tet centre.
-    scalar lambdaMin = VGREAT;
-
-    DynamicList<label>& tris = cloud.labels();
-
-    // Tet indices that will be set by hitWallFaces if a wall face is
-    // to be hit, or are set when any wall tri of a tet is hit.
-    // Carries the description of the tet on which the cell face has
-    // been hit.  For the case of being set in hitWallFaces, this may
-    // be a different tet to the one that the particle occupies.
-    tetIndices faceHitTetIs;
-
-    // What tolerance is appropriate the minimum lambda numerator and
-    // denominator for tracking in this cell.
-    scalar lambdaDistanceTolerance =
-        lambdaDistanceToleranceCoeff*mesh_.cellVolumes()[celli_];
-
-    do
+    // If the track is complete, return
+    if (!onFace())
     {
-        if (triI != -1)
-        {
-            // Change tet ownership because a tri face has been crossed
-            tetNeighbour(triI);
-        }
-
-        const Foam::face& f = pFaces[tetFacei_];
-
-        bool own = (mesh_.faceOwner()[tetFacei_] == celli_);
-
-        label tetBasePtI = mesh_.tetBasePtIs()[tetFacei_];
-
-        label basePtI = f[tetBasePtI];
-
-        label facePtI = (tetPti_ + tetBasePtI) % f.size();
-        label otherFacePtI = f.fcIndex(facePtI);
-
-        label fPtAI = -1;
-        label fPtBI = -1;
-
-        if (own)
-        {
-            fPtAI = facePtI;
-            fPtBI = otherFacePtI;
-        }
-        else
-        {
-            fPtAI = otherFacePtI;
-            fPtBI = facePtI;
-        }
-
-        tetPointRef tet
-        (
-            pC[celli_],
-            pPts[basePtI],
-            pPts[f[fPtAI]],
-            pPts[f[fPtBI]]
-        );
-
-        if (lambdaMin < SMALL)
-        {
-            // Apply tracking correction towards tet centre
-
-            if (debug)
-            {
-                Pout<< "tracking rescue using tetCentre from " << position();
-            }
-
-            position_ += trackingCorrectionTol*(tet.centre() - position_);
-
-            if (debug)
-            {
-                Pout<< " to " << position() << " due to "
-                    << (tet.centre() - position_) << endl;
-            }
-
-            cloud.trackingRescue();
-
-            return trackFraction;
-        }
-
-        if (triI != -1 && mesh_.moving())
-        {
-            // Mesh motion requires stepFraction to be correct for
-            // each tracking portion, so trackToFace must return after
-            // every lambda calculation.
-            return trackFraction;
-        }
-
-        FixedList<vector, 4> tetAreas;
-
-        tetAreas[0] = tet.Sa();
-        tetAreas[1] = tet.Sb();
-        tetAreas[2] = tet.Sc();
-        tetAreas[3] = tet.Sd();
-
-        FixedList<label, 4> tetPlaneBasePtIs;
-
-        tetPlaneBasePtIs[0] = basePtI;
-        tetPlaneBasePtIs[1] = f[fPtAI];
-        tetPlaneBasePtIs[2] = basePtI;
-        tetPlaneBasePtIs[3] = basePtI;
-
-        findTris
-        (
-            endPosition,
-            tris,
-            tet,
-            tetAreas,
-            tetPlaneBasePtIs,
-            lambdaDistanceTolerance
-        );
-
-        // Reset variables for new track
-        triI = -1;
-        lambdaMin = VGREAT;
-
-        // Pout<< "tris " << tris << endl;
-
-        // Sets a value for lambdaMin and facei_ if a wall face is hit
-        // by the track.
-        hitWallFaces
-        (
-            cloud,
-            position_,
-            endPosition,
-            lambdaMin,
-            faceHitTetIs
-        );
-
-        // Did not hit any tet tri faces, and no wall face has been
-        // found to hit.
-        if (tris.empty() && facei_ < 0)
-        {
-            position_ = endPosition;
-
-            return 1.0;
-        }
-        else
-        {
-            // Loop over all found tris and see if any of them find a
-            // lambda value smaller than that found for a wall face.
-            forAll(tris, i)
-            {
-                label tI = tris[i];
-
-                scalar lam = tetLambda
-                (
-                    position_,
-                    endPosition,
-                    tI,
-                    tetAreas[tI],
-                    tetPlaneBasePtIs[tI],
-                    celli_,
-                    tetFacei_,
-                    tetPti_,
-                    lambdaDistanceTolerance
-                );
-
-                if (lam < lambdaMin)
-                {
-                    lambdaMin = lam;
-
-                    triI = tI;
-                }
-            }
-        }
-
-        if (triI == 0)
-        {
-            // This must be a cell face crossing
-            facei_ = tetFacei_;
-
-            // Set the faceHitTetIs to those for the current tet in case a
-            // wall interaction is required with the cell face
-            faceHitTetIs = tetIndices
-            (
-                celli_,
-                tetFacei_,
-                tetBasePtI,
-                fPtAI,
-                fPtBI,
-                tetPti_
-            );
-        }
-        else if (triI > 0)
-        {
-            // A tri was found to be crossed before a wall face was hit (if any)
-            facei_ = -1;
-        }
-
-        // Pout<< "track loop " << position_ << " " << endPosition << nl
-        //     << "    " << celli_
-        //     << "    " << facei_
-        //     << " " << tetFacei_
-        //     << " " << tetPti_
-        //     << " " << triI
-        //     << " " << lambdaMin
-        //     << " " << trackFraction
-        //     << endl;
-
-        // Pout<< "# Tracking loop tet "
-        //     << origId_ << " " << origProc_<< nl
-        //     << "# face: " << tetFacei_ << nl
-        //     << "# tetPti: " << tetPti_ << nl
-        //     << "# tetBasePtI: " << mesh_.tetBasePtIs()[tetFacei_] << nl
-        //     << "# tet.mag(): " << tet.mag() << nl
-        //     << "# tet.quality(): " << tet.quality()
-        //     << endl;
-
-        // meshTools::writeOBJ(Pout, tet.a());
-        // meshTools::writeOBJ(Pout, tet.b());
-        // meshTools::writeOBJ(Pout, tet.c());
-        // meshTools::writeOBJ(Pout, tet.d());
-
-        // Pout<< "f 1 3 2" << nl
-        //     << "f 2 3 4" << nl
-        //     << "f 1 4 3" << nl
-        //     << "f 1 2 4" << endl;
-
-        // The particle can be 'outside' the tet.  This will yield a
-        // lambda larger than 1, or smaller than 0.  For values < 0,
-        // the particle travels away from the tet and we don't move
-        // the particle, only change tet/cell.  For values larger than
-        // 1, we move the particle to endPosition before the tet/cell
-        // change.
-        if (lambdaMin > SMALL)
-        {
-            if (lambdaMin <= 1.0)
-            {
-                trackFraction += lambdaMin*(1 - trackFraction);
-                position_ += lambdaMin*(endPosition - position_);
-            }
-            else
-            {
-                position_ = endPosition;
-
-                return 1.0;
-            }
-        }
-        else
-        {
-            // Set lambdaMin to zero to force a towards-tet-centre
-            // correction.
-            lambdaMin = 0.0;
-        }
-
-    } while (facei_ < 0);
+        return;
+    }
 
+    // Hit face/patch processing
+    typedef typename TrackData::cloudType::particleType particleType;
     particleType& p = static_cast<particleType&>(*this);
     p.hitFace(td);
-
-    if (internalFace(facei_))
+    if (onInternalFace())
     {
-        // Change tet ownership because a tri face has been crossed,
-        // in general this is:
-        //     tetNeighbour(triI);
-        // but triI must be 0;
-        // No modifications are required for triI = 0, no call required to
-        //     tetNeighbour(0);
-
-        if (celli_ == mesh_.faceOwner()[facei_])
-        {
-            celli_ = mesh_.faceNeighbour()[facei_];
-        }
-        else if (celli_ == mesh_.faceNeighbour()[facei_])
-        {
-            celli_ = mesh_.faceOwner()[facei_];
-        }
-        else
-        {
-            FatalErrorInFunction
-                << "addressing failure" << abort(FatalError);
-        }
+        changeCell();
     }
     else
     {
         label origFacei = facei_;
-        label patchi = patch(facei_);
+        label patchi = mesh_.boundaryMesh().whichPatch(facei_);
 
-        // No action taken for tetPti_ for tetFacei_ here, handled by
-        // patch interaction call or later during processor transfer.
+        // No action is taken for tetPti_ for tetFacei_ here. These are handled
+        // by the patch interaction call or later during processor transfer.
+
+        const tetIndices faceHitTetIs =
+            polyMeshTetDecomposition::triangleTetIndices
+            (
+                mesh_,
+                tetFacei_,
+                celli_,
+                tetPti_
+            );
 
         if
         (
-            !p.hitPatch
+           !p.hitPatch
             (
                 mesh_.boundaryMesh()[patchi],
                 td,
                 patchi,
-                trackFraction,
+                stepFraction(),
                 faceHitTetIs
             )
         )
@@ -596,7 +242,7 @@ Foam::scalar Foam::particle::trackToFace
             // Did patch interaction model switch patches?
             if (facei_ != origFacei)
             {
-                patchi = patch(facei_);
+                patchi = mesh_.boundaryMesh().whichPatch(facei_);
             }
 
             const polyPatch& patch = mesh_.boundaryMesh()[patchi];
@@ -635,7 +281,7 @@ Foam::scalar Foam::particle::trackToFace
                 (
                     static_cast<const cyclicAMIPolyPatch&>(patch),
                     td,
-                    endPosition - position_
+                    displacement
                 );
             }
             else if (isA<processorPolyPatch>(patch))
@@ -658,290 +304,6 @@ Foam::scalar Foam::particle::trackToFace
             }
         }
     }
-
-    if (lambdaMin < SMALL)
-    {
-        // Apply tracking correction towards tet centre.
-        // Generate current tet to find centre to apply correction.
-
-        tetPointRef tet = currentTet();
-
-        if (debug)
-        {
-            Pout<< "tracking rescue for lambdaMin:" << lambdaMin
-                << "from " << position();
-        }
-
-        position_ += trackingCorrectionTol*(tet.centre() - position_);
-
-        if
-        (
-            cloud.hasWallImpactDistance()
-         && !internalFace(faceHitTetIs.face())
-         && cloud.cellHasWallFaces()[faceHitTetIs.cell()]
-        )
-        {
-            const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-            label fI = faceHitTetIs.face();
-
-            label patchi = patches.patchID()[fI - mesh_.nInternalFaces()];
-
-            if (isA<wallPolyPatch>(patches[patchi]))
-            {
-                // In the case of collision with a wall where there is
-                // a non-zero wallImpactDistance, it is possible for
-                // there to be a tracking correction required to bring
-                // the particle into the domain, but the position of
-                // the particle is further from the wall than the tet
-                // centre, in which case the normal correction can be
-                // counter-productive, i.e. pushes the particle
-                // further out of the domain.  In this case it is the
-                // position that hit the wall that is in need of a
-                // rescue correction.
-
-                triPointRef wallTri = faceHitTetIs.faceTri(mesh_);
-
-                tetPointRef wallTet = faceHitTetIs.tet(mesh_);
-
-                vector nHat = wallTri.normal();
-                nHat /= mag(nHat);
-
-                const scalar r = p.wallImpactDistance(nHat);
-
-                // Removing (approximately) the wallTri normal
-                // component of the existing correction, to avoid the
-                // situation where the existing correction in the wall
-                // normal direction is larger towards the wall than
-                // the new correction is away from it.
-                position_ +=
-                    trackingCorrectionTol
-                   *(
-                        (wallTet.centre() - (position_ + r*nHat))
-                      - (nHat & (tet.centre() - position_))*nHat
-                    );
-            }
-        }
-
-        if (debug)
-        {
-            Pout<< " to " << position() << endl;
-        }
-
-        cloud.trackingRescue();
-    }
-
-    return trackFraction;
-}
-
-
-template<class CloudType>
-void Foam::particle::hitWallFaces
-(
-    const CloudType& cloud,
-    const vector& from,
-    const vector& to,
-    scalar& lambdaMin,
-    tetIndices& closestTetIs
-)
-{
-    typedef typename CloudType::particleType particleType;
-
-    if (!(cloud.hasWallImpactDistance() && cloud.cellHasWallFaces()[celli_]))
-    {
-        return;
-    }
-
-    particleType& p = static_cast<particleType&>(*this);
-
-    const faceList& pFaces = mesh_.faces();
-
-    const Foam::cell& thisCell = mesh_.cells()[celli_];
-
-    scalar lambdaDistanceTolerance =
-        lambdaDistanceToleranceCoeff*mesh_.cellVolumes()[celli_];
-
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-    forAll(thisCell, cFI)
-    {
-        label fI = thisCell[cFI];
-
-        if (internalFace(fI))
-        {
-            continue;
-        }
-
-        label patchi = patches.patchID()[fI - mesh_.nInternalFaces()];
-
-        if (isA<wallPolyPatch>(patches[patchi]))
-        {
-            // Get the decomposition of this wall face
-
-            const List<tetIndices> faceTetIs =
-                polyMeshTetDecomposition::faceTetIndices(mesh_, fI, celli_);
-
-            const Foam::face& f = pFaces[fI];
-
-            forAll(faceTetIs, tI)
-            {
-                const tetIndices& tetIs = faceTetIs[tI];
-
-                triPointRef tri = tetIs.faceTri(mesh_);
-
-                vector n = tri.normal();
-
-                vector nHat = n/mag(n);
-
-                // Radius of particle with respect to this wall face
-                // triangle.  Assuming that the wallImpactDistance
-                // does not change as the particle or the mesh moves
-                // forward in time.
-                scalar r = p.wallImpactDistance(nHat);
-
-                vector toPlusRNHat = to + r*nHat;
-
-                // triI = 0 because it is the cell face tri of the tet
-                // we are concerned with.
-                scalar tetClambda = tetLambda
-                (
-                    tetIs.tet(mesh_).centre(),
-                    toPlusRNHat,
-                    0,
-                    n,
-                    f[tetIs.faceBasePt()],
-                    celli_,
-                    fI,
-                    tetIs.tetPt(),
-                    lambdaDistanceTolerance
-                );
-
-                if ((tetClambda <= 0.0) || (tetClambda >= 1.0))
-                {
-                    // toPlusRNHat is not on the outside of the plane of
-                    // the wall face tri, the tri cannot be hit.
-                    continue;
-                }
-
-                // Check if the actual trajectory of the near-tri
-                // points intersects the triangle.
-
-                vector fromPlusRNHat = from + r*nHat;
-
-                // triI = 0 because it is the cell face tri of the tet
-                // we are concerned with.
-                scalar lambda = tetLambda
-                (
-                    fromPlusRNHat,
-                    toPlusRNHat,
-                    0,
-                    n,
-                    f[tetIs.faceBasePt()],
-                    celli_,
-                    fI,
-                    tetIs.tetPt(),
-                    lambdaDistanceTolerance
-                );
-
-                pointHit hitInfo(Zero);
-
-                if (mesh_.moving())
-                {
-                    // For a moving mesh, the position of wall
-                    // triangle needs to be moved in time to be
-                    // consistent with the moment defined by the
-                    // current value of stepFraction and the value of
-                    // lambda just calculated.
-
-                    // Total fraction thought the timestep of the
-                    // motion, including stepFraction before the
-                    // current tracking step and the current
-                    // lambda
-                    // i.e.
-                    // let s = stepFraction, l = lambda
-                    // Motion of x in time:
-                    // |-----------------|---------|---------|
-                    // x00               x0        xi        x
-                    //
-                    // where xi is the correct value of x at the required
-                    // tracking instant.
-                    //
-                    // x0 = x00 + s*(x - x00) = s*x + (1 - s)*x00
-                    //
-                    // i.e. the motion covered by previous tracking portions
-                    // within this timestep, and
-                    //
-                    // xi = x0 + l*(x - x0)
-                    //    = l*x + (1 - l)*x0
-                    //    = l*x + (1 - l)*(s*x + (1 - s)*x00)
-                    //    = (s + l - s*l)*x + (1 - (s + l - s*l))*x00
-                    //
-                    // let m = (s + l - s*l)
-                    //
-                    // xi = m*x + (1 - m)*x00 = x00 + m*(x - x00);
-                    //
-                    // In the same form as before.
-
-                    // Clip lambda to 0.0-1.0 to ensure that sensible
-                    // positions are used for triangle intersections.
-                    scalar lam = max(0.0, min(1.0, lambda));
-
-                    scalar m = stepFraction_ + lam - (stepFraction_*lam);
-
-                    triPointRef tri00 = tetIs.oldFaceTri(mesh_);
-
-                    // Use SMALL positive tolerance to make the triangle
-                    // slightly "fat" to improve robustness.  Intersection
-                    // is calculated as the ray (from + r*nHat) -> (to +
-                    // r*nHat).
-
-                    point tPtA = tri00.a() + m*(tri.a() - tri00.a());
-                    point tPtB = tri00.b() + m*(tri.b() - tri00.b());
-                    point tPtC = tri00.c() + m*(tri.c() - tri00.c());
-
-                    triPointRef t(tPtA, tPtB, tPtC);
-
-                    // The point fromPlusRNHat + m*(to - from) is on the
-                    // plane of the triangle.  Determine the
-                    // intersection with this triangle by testing if
-                    // this point is inside or outside of the triangle.
-                    hitInfo = t.intersection
-                    (
-                        fromPlusRNHat + m*(to - from),
-                        t.normal(),
-                        intersection::FULL_RAY,
-                        SMALL
-                    );
-                }
-                else
-                {
-                    // Use SMALL positive tolerance to make the triangle
-                    // slightly "fat" to improve robustness.  Intersection
-                    // is calculated as the ray (from + r*nHat) -> (to +
-                    // r*nHat).
-                    hitInfo = tri.intersection
-                    (
-                        fromPlusRNHat,
-                        (to - from),
-                        intersection::FULL_RAY,
-                        SMALL
-                    );
-                }
-
-                if (hitInfo.hit())
-                {
-                    if (lambda < lambdaMin)
-                    {
-                        lambdaMin = lambda;
-
-                        facei_ = fI;
-
-                        closestTetIs = tetIs;
-                    }
-                }
-            }
-        }
-    }
 }
 
 
@@ -1017,22 +379,17 @@ void Foam::particle::hitCyclicPatch
     TrackData& td
 )
 {
-    facei_ = cpp.transformGlobalFace(facei_);
+    const cyclicPolyPatch& receiveCpp = cpp.neighbPatch();
+    const label receiveFacei = receiveCpp.whichFace(facei_);
 
+    // Set the topology
+    facei_ = tetFacei_ = cpp.transformGlobalFace(facei_);
     celli_ = mesh_.faceOwner()[facei_];
-
-    tetFacei_ = facei_;
-
-    // See note in correctAfterParallelTransfer for tetPti_ addressing.
+    // See note in correctAfterParallelTransfer for tetPti addressing ...
     tetPti_ = mesh_.faces()[tetFacei_].size() - 1 - tetPti_;
 
-    const cyclicPolyPatch& receiveCpp = cpp.neighbPatch();
-    label patchFacei = receiveCpp.whichFace(facei_);
-
-    // Now the particle is on the receiving side
-
-    // Have patch transform the position
-    receiveCpp.transformPosition(position_, patchFacei);
+    // Reflect to account for the change of triangle orientation in the new cell
+    reflect();
 
     // Transform the properties
     if (!receiveCpp.parallel())
@@ -1041,7 +398,7 @@ void Foam::particle::hitCyclicPatch
         (
             receiveCpp.forwardT().size() == 1
           ? receiveCpp.forwardT()[0]
-          : receiveCpp.forwardT()[patchFacei]
+          : receiveCpp.forwardT()[receiveFacei]
         );
         transformProperties(T);
     }
@@ -1051,7 +408,7 @@ void Foam::particle::hitCyclicPatch
         (
             (receiveCpp.separation().size() == 1)
           ? receiveCpp.separation()[0]
-          : receiveCpp.separation()[patchFacei]
+          : receiveCpp.separation()[receiveFacei]
         );
         transformProperties(-s);
     }
@@ -1066,38 +423,41 @@ void Foam::particle::hitCyclicAMIPatch
     const vector& direction
 )
 {
+    vector pos = position();
+
     const cyclicAMIPolyPatch& receiveCpp = cpp.neighbPatch();
+    const label sendFacei = cpp.whichFace(facei_);
+    const label receiveFacei = cpp.pointFace(sendFacei, direction, pos);
 
-    // Patch face index on sending side
-    label patchFacei = facei_ - cpp.start();
-
-    // Patch face index on receiving side - also updates position
-    patchFacei = cpp.pointFace(patchFacei, direction, position_);
-
-    if (patchFacei < 0)
+    if (receiveFacei < 0)
     {
-        // If the patch face of the particle is not known assume that
-        // the particle is lost and to be deleted
+        // If the patch face of the particle is not known assume that the
+        // particle is lost and mark it to be deleted.
         td.keepParticle = false;
-
         WarningInFunction
             << "Particle lost across " << cyclicAMIPolyPatch::typeName
             << " patches " << cpp.name() << " and " << receiveCpp.name()
-            << " at position " << position_
-            << endl;
+            << " at position " << pos << endl;
     }
 
-    // Convert face index into global numbering
-    facei_ = patchFacei + receiveCpp.start();
+    // Set the topology
+    facei_ = tetFacei_ = receiveFacei + receiveCpp.start();
 
-    celli_ = mesh_.faceOwner()[facei_];
+    // Locate the particle on the recieving side
+    locate
+    (
+        pos,
+        &direction,
+        mesh_.faceOwner()[facei_],
+        false,
+        "Particle crossed between " + cyclicAMIPolyPatch::typeName +
+        " patches " + cpp.name() + " and " + receiveCpp.name() +
+        " to a location outside of the mesh."
+    );
 
-    tetFacei_ = facei_;
-
-    // See note in correctAfterParallelTransfer for tetPti_ addressing.
-    tetPti_ = mesh_.faces()[tetFacei_].size() - 1 - tetPti_;
-
-    // Now the particle is on the receiving side
+    // The particle must remain associated with a face for the tracking to
+    // register as incomplete
+    facei_ = tetFacei_;
 
     // Transform the properties
     if (!receiveCpp.parallel())
@@ -1106,7 +466,7 @@ void Foam::particle::hitCyclicAMIPatch
         (
             receiveCpp.forwardT().size() == 1
           ? receiveCpp.forwardT()[0]
-          : receiveCpp.forwardT()[patchFacei]
+          : receiveCpp.forwardT()[receiveFacei]
         );
         transformProperties(T);
     }
@@ -1116,7 +476,7 @@ void Foam::particle::hitCyclicAMIPatch
         (
             (receiveCpp.separation().size() == 1)
           ? receiveCpp.separation()[0]
-          : receiveCpp.separation()[patchFacei]
+          : receiveCpp.separation()[receiveFacei]
         );
         transformProperties(-s);
     }
diff --git a/src/lagrangian/basic/passiveParticle/passiveParticle.H b/src/lagrangian/basic/passiveParticle/passiveParticle.H
index 0624039fdc..35ef33c0ae 100644
--- a/src/lagrangian/basic/passiveParticle/passiveParticle.H
+++ b/src/lagrangian/basic/passiveParticle/passiveParticle.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -61,26 +61,26 @@ public:
         passiveParticle
         (
             const polyMesh& mesh,
-            const vector& position,
+            const barycentric& coordinates,
             const label celli,
             const label tetFacei,
             const label tetPti
         )
         :
-            particle(mesh, position, celli, tetFacei, tetPti)
+            particle(mesh, coordinates, celli, tetFacei, tetPti)
         {}
 
-        //- Construct from components, with searching for tetFace and
-        //  tetPt unless disabled by doCellFacePt = false.
+
+        //- Construct from a position and a cell, searching for the rest of the
+        //  required topology
         passiveParticle
         (
             const polyMesh& mesh,
             const vector& position,
-            const label celli,
-            bool doCellFacePt = true
+            const label celli
         )
         :
-            particle(mesh, position, celli, doCellFacePt)
+            particle(mesh, position, celli)
         {}
 
         //- Construct from Istream
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index 0dfb75094d..6b22c85869 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -318,7 +318,7 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
     ),
     rndGen_(Pstream::myProcNo()),
     cellOccupancyPtr_(),
-    cellLengthScale_(cbrt(mesh_.V())),
+    cellLengthScale_(mag(cbrt(mesh_.V()))),
     rho_(rho),
     U_(U),
     mu_(mu),
@@ -845,18 +845,14 @@ void Foam::KinematicCloud<CloudType>::updateMesh()
 {
     updateCellOccupancy();
     injectors_.updateMesh();
-    cellLengthScale_ = cbrt(mesh_.V());
+    cellLengthScale_ = mag(cbrt(mesh_.V()));
 }
 
 
 template<class CloudType>
 void Foam::KinematicCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
 {
-    typedef typename particle::TrackingData<KinematicCloud<CloudType>> tdType;
-
-    tdType td(*this);
-
-    Cloud<parcelType>::template autoMap<tdType>(td, mapper);
+    Cloud<parcelType>::autoMap(mapper);
 
     updateMesh();
 }
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
index 6447662255..04862d198e 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -333,11 +333,7 @@ void Foam::ReactingCloud<CloudType>::evolve()
 template<class CloudType>
 void Foam::ReactingCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
 {
-    typedef typename particle::TrackingData<ReactingCloud<CloudType>> tdType;
-
-    tdType td(*this);
-
-    Cloud<parcelType>::template autoMap<tdType>(td, mapper);
+    Cloud<parcelType>::autoMap(mapper);
 
     this->updateMesh();
 }
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
index 21b03b25dc..3d9d76d60c 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -259,12 +259,7 @@ void Foam::ReactingMultiphaseCloud<CloudType>::autoMap
     const mapPolyMesh& mapper
 )
 {
-    typedef typename particle::TrackingData<ReactingMultiphaseCloud<CloudType>>
-        tdType;
-
-    tdType td(*this);
-
-    Cloud<parcelType>::template autoMap<tdType>(td, mapper);
+    Cloud<parcelType>::autoMap(mapper);
 
     this->updateMesh();
 }
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
index cde6a17b62..45f683308b 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -480,11 +480,7 @@ void Foam::ThermoCloud<CloudType>::evolve()
 template<class CloudType>
 void Foam::ThermoCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
 {
-    typedef typename particle::TrackingData<ThermoCloud<CloudType>> tdType;
-
-    tdType td(*this);
-
-    Cloud<parcelType>::template autoMap<tdType>(td, mapper);
+    Cloud<parcelType>::autoMap(mapper);
 
     this->updateMesh();
 }
diff --git a/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.H b/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.H
index e1a8c7544e..ed405f81c1 100644
--- a/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -172,22 +172,31 @@ public:
 
     // Constructors
 
-        //- Construct from owner, position, and cloud owner
+        //- Construct from mesh, coordinates and topology
         //  Other properties initialised as null
         inline CollidingParcel
         (
             const polyMesh& mesh,
-            const vector& position,
+            const barycentric& coordinates,
             const label celli,
             const label tetFacei,
             const label tetPti
         );
 
+        //- Construct from a position and a cell, searching for the rest of the
+        //  required topology. Other properties are initialised as null.
+        inline CollidingParcel
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label celli
+        );
+
         //- Construct from components
         inline CollidingParcel
         (
             const polyMesh& mesh,
-            const vector& position,
+            const barycentric& coordinates,
             const label celli,
             const label tetFacei,
             const label tetPti,
diff --git a/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcelI.H b/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcelI.H
index c40c4e5d54..c817a75bde 100644
--- a/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcelI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -63,13 +63,13 @@ template<class ParcelType>
 inline Foam::CollidingParcel<ParcelType>::CollidingParcel
 (
     const polyMesh& owner,
-    const vector& position,
+    const barycentric& coordinates,
     const label celli,
     const label tetFacei,
     const label tetPti
 )
 :
-    ParcelType(owner, position, celli, tetFacei, tetPti),
+    ParcelType(owner, coordinates, celli, tetFacei, tetPti),
     f_(Zero),
     angularMomentum_(Zero),
     torque_(Zero),
@@ -82,6 +82,22 @@ inline Foam::CollidingParcel<ParcelType>::CollidingParcel
 (
     const polyMesh& owner,
     const vector& position,
+    const label celli
+)
+:
+    ParcelType(owner, position, celli),
+    f_(Zero),
+    angularMomentum_(Zero),
+    torque_(Zero),
+    collisionRecords_()
+{}
+
+
+template<class ParcelType>
+inline Foam::CollidingParcel<ParcelType>::CollidingParcel
+(
+    const polyMesh& owner,
+    const barycentric& coordinates,
     const label celli,
     const label tetFacei,
     const label tetPti,
@@ -99,7 +115,7 @@ inline Foam::CollidingParcel<ParcelType>::CollidingParcel
     ParcelType
     (
         owner,
-        position,
+        coordinates,
         celli,
         tetFacei,
         tetPti,
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index b593d4393b..c2d11aee3c 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -268,71 +268,50 @@ bool Foam::KinematicParcel<ParcelType>::move
     const cloudSolution& solution = td.cloud().solution();
     const scalarField& cellLengthScale = td.cloud().cellLengthScale();
 
-    scalar tEnd = (1.0 - p.stepFraction())*trackTime;
-    scalar dtMax = solution.deltaTMax(trackTime);
-
-    bool tracking = true;
-    label nTrackingStalled = 0;
-
-    while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
+    while (td.keepParticle && !td.switchProcessor && p.stepFraction() < 1)
     {
         // Apply correction to position for reduced-D cases
-        meshTools::constrainToMeshCentre(mesh, p.position());
+        p.constrainToMeshCentre();
 
-        const point start(p.position());
-
-        // Set the Lagrangian time-step
-        scalar dt = min(dtMax, tEnd);
-
-        // Cache the parcel current cell as this will change if a face is hit
+        // Cache the current position, cell and step-fraction
+        const point start = p.position();
         const label celli = p.cell();
+        const scalar sfrac = p.stepFraction();
 
-        const scalar magU = mag(U_);
-        if (p.active() && tracking && (magU > ROOTVSMALL))
+        // Total displacement over the time-step
+        const vector s = trackTime*U_;
+
+        // Cell length scale
+        const scalar l = cellLengthScale[p.cell()];
+
+        // Fraction of the displacement to track in this loop. This is limited
+        // to ensure that the both the time and distance tracked is less than
+        // maxCo times the total value.
+        scalar f = 1 - p.stepFraction();
+        f = min(f, maxCo);
+        f = min(f, maxCo*l/max(SMALL*l, mag(s)));
+        if (p.active())
         {
-            const scalar d = dt*magU;
-            const scalar deltaLMax = solution.deltaLMax(cellLengthScale[celli]);
-            const scalar dCorr = min(d, deltaLMax);
-            dt *=
-                dCorr/d
-               *p.trackToFace(p.position() + dCorr*U_/magU, td);
+            // Track to the next face
+            p.trackToFace(f*s, f, td);
         }
-
-        tEnd -= dt;
-
-        const scalar newStepFraction = 1.0 - tEnd/trackTime;
-
-        if (tracking)
+        else
         {
-            if
-            (
-                mag(p.stepFraction() - newStepFraction)
-              < particle::minStepFractionTol
-            )
+            // Abandon the track, and move to the end of the sub-step. If the
+            // the mesh is moving, this will implicitly move the parcel.
+            if (mesh.moving())
             {
-                nTrackingStalled++;
-
-                if (nTrackingStalled > maxTrackAttempts)
-                {
-                    tracking = false;
-                }
-            }
-            else
-            {
-                nTrackingStalled = 0;
+                WarningInFunction
+                    << "Tracking was abandoned on a moving mesh. Parcels may "
+                    << "move unphysically as a result." << endl;
             }
+            p.stepFraction() += f;
         }
 
-        p.stepFraction() = newStepFraction;
-
-        bool calcParcel = true;
-        if (!tracking && solution.steadyState())
-        {
-            calcParcel = false;
-        }
+        const scalar dt = (p.stepFraction() - sfrac)*trackTime;
 
         // Avoid problems with extremely small timesteps
-        if ((dt > ROOTVSMALL) && calcParcel)
+        if (!td.cloud().solution().steadyState() && dt > ROOTVSMALL)
         {
             // Update cell based properties
             p.setCellValues(td, dt, celli);
@@ -345,7 +324,7 @@ bool Foam::KinematicParcel<ParcelType>::move
             p.calc(td, dt, celli);
         }
 
-        if (p.onBoundary() && td.keepParticle)
+        if (p.onBoundaryFace() && td.keepParticle)
         {
             if (isA<processorPolyPatch>(pbMesh[p.patch(p.face())]))
             {
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
index dd45972ca2..06a1e879db 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -340,22 +340,31 @@ public:
 
     // Constructors
 
-        //- Construct from owner, position, and cloud owner
+        //- Construct from mesh, coordinates and topology
         //  Other properties initialised as null
         inline KinematicParcel
         (
             const polyMesh& mesh,
-            const vector& position,
+            const barycentric& coordinates,
             const label celli,
             const label tetFacei,
             const label tetPti
         );
 
+        //- Construct from a position and a cell, searching for the rest of the
+        //  required topology. Other properties are initialised as null.
+        inline KinematicParcel
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label celli
+        );
+
         //- Construct from components
         inline KinematicParcel
         (
             const polyMesh& mesh,
-            const vector& position,
+            const barycentric& coordinates,
             const label celli,
             const label tetFacei,
             const label tetPti,
@@ -495,9 +504,6 @@ public:
 
         // Helper functions
 
-            //- Return the index of the face used in the interpolation routine
-            inline label faceInterpolation() const;
-
             //- Cell owner mass
             inline scalar massCell(const label celli) const;
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
index a0e60d35ad..749c136883 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -73,13 +73,13 @@ template<class ParcelType>
 inline Foam::KinematicParcel<ParcelType>::KinematicParcel
 (
     const polyMesh& owner,
-    const vector& position,
+    const barycentric& coordinates,
     const label celli,
     const label tetFacei,
     const label tetPti
 )
 :
-    ParcelType(owner, position, celli, tetFacei, tetPti),
+    ParcelType(owner, coordinates, celli, tetFacei, tetPti),
     active_(true),
     typeId_(-1),
     nParticle_(0),
@@ -101,6 +101,31 @@ inline Foam::KinematicParcel<ParcelType>::KinematicParcel
 (
     const polyMesh& owner,
     const vector& position,
+    const label celli
+)
+:
+    ParcelType(owner, position, celli),
+    active_(true),
+    typeId_(-1),
+    nParticle_(0),
+    d_(0.0),
+    dTarget_(0.0),
+    U_(Zero),
+    rho_(0.0),
+    age_(0.0),
+    tTurb_(0.0),
+    UTurb_(Zero),
+    rhoc_(0.0),
+    Uc_(Zero),
+    muc_(0.0)
+{}
+
+
+template<class ParcelType>
+inline Foam::KinematicParcel<ParcelType>::KinematicParcel
+(
+    const polyMesh& owner,
+    const barycentric& coordinates,
     const label celli,
     const label tetFacei,
     const label tetPti,
@@ -112,7 +137,7 @@ inline Foam::KinematicParcel<ParcelType>::KinematicParcel
     const constantProperties& constProps
 )
 :
-    ParcelType(owner, position, celli, tetFacei, tetPti),
+    ParcelType(owner, coordinates, celli, tetFacei, tetPti),
     active_(true),
     typeId_(typeId),
     nParticle_(nParticle0),
@@ -334,21 +359,6 @@ inline Foam::vector& Foam::KinematicParcel<ParcelType>::UTurb()
 }
 
 
-template<class ParcelType>
-inline Foam::label Foam::KinematicParcel<ParcelType>::faceInterpolation() const
-{
-    // Use volume-based interpolation if dealing with external faces
-    if (this->cloud().internalFace(this->face()))
-    {
-        return this->face();
-    }
-    else
-    {
-        return -1;
-    }
-}
-
-
 template<class ParcelType>
 inline Foam::scalar Foam::KinematicParcel<ParcelType>::massCell
 (
diff --git a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcel.H b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcel.H
index c932590a9c..9bae31a577 100644
--- a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -190,22 +190,31 @@ public:
 
     // Constructors
 
-        //- Construct from owner, position, and cloud owner
+        //- Construct from mesh, coordinates and topology
         //  Other properties initialised as null
         inline MPPICParcel
         (
             const polyMesh& mesh,
-            const vector& position,
+            const barycentric& coordinates,
             const label celli,
             const label tetFacei,
             const label tetPti
         );
 
+        //- Construct from a position and a cell, searching for the rest of the
+        //  required topology. Other properties are initialised as null.
+        inline MPPICParcel
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label celli
+        );
+
         //- Construct from components
         inline MPPICParcel
         (
             const polyMesh& mesh,
-            const vector& position,
+            const barycentric& coordinates,
             const label celli,
             const label tetFacei,
             const label tetPti,
diff --git a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelI.H b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelI.H
index 1af5110019..f4fd43e14b 100644
--- a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,13 +29,13 @@ template<class ParcelType>
 inline Foam::MPPICParcel<ParcelType>::MPPICParcel
 (
     const polyMesh& owner,
-    const vector& position,
+    const barycentric& coordinates,
     const label celli,
     const label tetFacei,
     const label tetPti
 )
 :
-    ParcelType(owner, position, celli, tetFacei, tetPti),
+    ParcelType(owner, coordinates, celli, tetFacei, tetPti),
     UCorrect_(Zero)
 {}
 
@@ -45,6 +45,19 @@ inline Foam::MPPICParcel<ParcelType>::MPPICParcel
 (
     const polyMesh& owner,
     const vector& position,
+    const label celli
+)
+:
+    ParcelType(owner, position, celli),
+    UCorrect_(Zero)
+{}
+
+
+template<class ParcelType>
+inline Foam::MPPICParcel<ParcelType>::MPPICParcel
+(
+    const polyMesh& owner,
+    const barycentric& coordinates,
     const label celli,
     const label tetFacei,
     const label tetPti,
@@ -60,7 +73,7 @@ inline Foam::MPPICParcel<ParcelType>::MPPICParcel
     ParcelType
     (
         owner,
-        position,
+        coordinates,
         celli,
         tetFacei,
         tetPti,
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
index c9525627d4..0aaaf25a2a 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -278,23 +278,31 @@ public:
 
     // Constructors
 
-        //- Construct from owner, position, and cloud owner
+        //- Construct from mesh, position and topology
         //  Other properties initialised as null
         inline ReactingMultiphaseParcel
         (
             const polyMesh& mesh,
-            const vector& position,
+            const barycentric& coordinates,
             const label celli,
             const label tetFacei,
             const label tetPti
         );
 
+        //- Construct from a position and a cell, searching for the rest of the
+        //  required topology. Other properties are initialised as null.
+        inline ReactingMultiphaseParcel
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label celli
+        );
 
         //- Construct from components
         inline ReactingMultiphaseParcel
         (
             const polyMesh& mesh,
-            const vector& position,
+            const barycentric& coordinates,
             const label celli,
             const label tetFacei,
             const label tetPti,
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H
index c64e29b262..f689e91a9e 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,13 +68,13 @@ template<class ParcelType>
 inline Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
 (
     const polyMesh& mesh,
-    const vector& position,
+    const barycentric& coordinates,
     const label celli,
     const label tetFacei,
     const label tetPti
 )
 :
-    ParcelType(mesh, position, celli, tetFacei, tetPti),
+    ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
     YGas_(0),
     YLiquid_(0),
     YSolid_(0),
@@ -87,6 +87,22 @@ inline Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
 (
     const polyMesh& mesh,
     const vector& position,
+    const label celli
+)
+:
+    ParcelType(mesh, position, celli),
+    YGas_(0),
+    YLiquid_(0),
+    YSolid_(0),
+    canCombust_(0)
+{}
+
+
+template<class ParcelType>
+inline Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
+(
+    const polyMesh& mesh,
+    const barycentric& coordinates,
     const label celli,
     const label tetFacei,
     const label tetPti,
@@ -108,7 +124,7 @@ inline Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
     ParcelType
     (
         mesh,
-        position,
+        coordinates,
         celli,
         tetFacei,
         tetPti,
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
index c5463d9296..6335ca9dbd 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -232,22 +232,31 @@ public:
 
     // Constructors
 
-        //- Construct from owner, position, and cloud owner
+        //- Construct from mesh, coordinates and topology
         //  Other properties initialised as null
         inline ReactingParcel
         (
             const polyMesh& mesh,
-            const vector& position,
+            const barycentric& coordinates,
             const label celli,
             const label tetFacei,
             const label tetPti
         );
 
+        //- Construct from a position and a cell, searching for the rest of the
+        //  required topology. Other properties are initialised as null.
+        inline ReactingParcel
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label celli
+        );
+
         //- Construct from components
         inline ReactingParcel
         (
             const polyMesh& mesh,
-            const vector& position,
+            const barycentric& coordinates,
             const label celli,
             const label tetFacei,
             const label tetPti,
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
index a6359c4006..08699da9b8 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -63,13 +63,13 @@ template<class ParcelType>
 inline Foam::ReactingParcel<ParcelType>::ReactingParcel
 (
     const polyMesh& mesh,
-    const vector& position,
+    const barycentric& coordinates,
     const label celli,
     const label tetFacei,
     const label tetPti
 )
 :
-    ParcelType(mesh, position, celli, tetFacei, tetPti),
+    ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
     mass0_(0.0),
     Y_(0),
     pc_(0.0)
@@ -81,6 +81,21 @@ inline Foam::ReactingParcel<ParcelType>::ReactingParcel
 (
     const polyMesh& mesh,
     const vector& position,
+    const label celli
+)
+:
+    ParcelType(mesh, position, celli),
+    mass0_(0.0),
+    Y_(0),
+    pc_(0.0)
+{}
+
+
+template<class ParcelType>
+inline Foam::ReactingParcel<ParcelType>::ReactingParcel
+(
+    const polyMesh& mesh,
+    const barycentric& coordinates,
     const label celli,
     const label tetFacei,
     const label tetPti,
@@ -99,7 +114,7 @@ inline Foam::ReactingParcel<ParcelType>::ReactingParcel
     ParcelType
     (
         mesh,
-        position,
+        coordinates,
         celli,
         tetFacei,
         tetPti,
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
index c0a6777e06..4129e25835 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -288,22 +288,31 @@ public:
 
     // Constructors
 
-        //- Construct from owner, position, and cloud owner
+        //- Construct from mesh, coordinates and topology
         //  Other properties initialised as null
         inline ThermoParcel
         (
             const polyMesh& mesh,
-            const vector& position,
+            const barycentric& coordinates,
             const label celli,
             const label tetFacei,
             const label tetPti
         );
 
+        //- Construct from a position and a cell, searching for the rest of the
+        //  required topology. Other properties are initialised as null.
+        inline ThermoParcel
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label celli
+        );
+
         //- Construct from components
         inline ThermoParcel
         (
             const polyMesh& mesh,
-            const vector& position,
+            const barycentric& coordinates,
             const label celli,
             const label tetFacei,
             const label tetPti,
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H
index 9b2f365b44..9c35a55a6f 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -74,13 +74,13 @@ template<class ParcelType>
 inline Foam::ThermoParcel<ParcelType>::ThermoParcel
 (
     const polyMesh& mesh,
-    const vector& position,
+    const barycentric& coordinates,
     const label celli,
     const label tetFacei,
     const label tetPti
 )
 :
-    ParcelType(mesh, position, celli, tetFacei, tetPti),
+    ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
     T_(0.0),
     Cp_(0.0),
     Tc_(0.0),
@@ -93,6 +93,22 @@ inline Foam::ThermoParcel<ParcelType>::ThermoParcel
 (
     const polyMesh& mesh,
     const vector& position,
+    const label celli
+)
+:
+    ParcelType(mesh, position, celli),
+    T_(0.0),
+    Cp_(0.0),
+    Tc_(0.0),
+    Cpc_(0.0)
+{}
+
+
+template<class ParcelType>
+inline Foam::ThermoParcel<ParcelType>::ThermoParcel
+(
+    const polyMesh& mesh,
+    const barycentric& coordinates,
     const label celli,
     const label tetFacei,
     const label tetPti,
@@ -110,7 +126,7 @@ inline Foam::ThermoParcel<ParcelType>::ThermoParcel
     ParcelType
     (
         mesh,
-        position,
+        coordinates,
         celli,
         tetFacei,
         tetPti,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
index 9f496a63f1..16568fcc9a 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -485,8 +485,7 @@ void Foam::InjectionModel<CloudType>::inject(TrackData& td)
                     meshTools::constrainToMeshCentre(mesh, pos);
 
                     // Create a new parcel
-                    parcelType* pPtr =
-                        new parcelType(mesh, pos, celli, tetFacei, tetPti);
+                    parcelType* pPtr = new parcelType(mesh, pos, celli);
 
                     // Check/set new parcel thermo properties
                     cloud.setParcelThermoProperties(*pPtr, dt);
@@ -610,8 +609,7 @@ void Foam::InjectionModel<CloudType>::injectSteadyState
             meshTools::constrainToMeshCentre(mesh, pos);
 
             // Create a new parcel
-            parcelType* pPtr =
-                new parcelType(mesh, pos, celli, tetFacei, tetPti);
+            parcelType* pPtr = new parcelType(mesh, pos, celli);
 
             // Check/set new parcel thermo properties
             cloud.setParcelThermoProperties(*pPtr, 0.0);
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C
index cca6266939..eb8fc4416f 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -146,16 +146,6 @@ void Foam::SurfaceFilmModel<CloudType>::inject(TrackData& td)
             {
                 const label celli = injectorCellsPatch[j];
 
-                // The position could bein any tet of the decomposed cell,
-                // so arbitrarily choose the first face of the cell as the
-                // tetFace and the first point on the face after the base
-                // point as the tetPt.  The tracking will pick the cell
-                // consistent with the motion in the first tracking step.
-                const label tetFacei = this->owner().mesh().cells()[celli][0];
-                const label tetPti = 1;
-
-//                const point& pos = this->owner().mesh().C()[celli];
-
                 const scalar offset =
                     max
                     (
@@ -166,14 +156,7 @@ void Foam::SurfaceFilmModel<CloudType>::inject(TrackData& td)
 
                 // Create a new parcel
                 parcelType* pPtr =
-                    new parcelType
-                    (
-                        this->owner().pMesh(),
-                        pos,
-                        celli,
-                        tetFacei,
-                        tetPti
-                    );
+                    new parcelType(this->owner().pMesh(), pos, celli);
 
                 // Check/set new parcel thermo properties
                 td.cloud().setParcelThermoProperties(*pPtr, 0.0);
diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
index 08237cb398..0cb1d22242 100644
--- a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
+++ b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
@@ -446,7 +446,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
         }
 
         // Perturb new parcels towards the owner cell centre
-        pPtr->position() += 0.5*rndGen_.sample01<scalar>()*(posC - posCf);
+        pPtr->track(0.5*rndGen_.sample01<scalar>()*(posC - posCf), 0);
 
         pPtr->nParticle() = npNew[i];
 
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
index a5ba420f34..f565adbfac 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -85,18 +85,10 @@ bool Foam::molecule::move(molecule::trackingData& td, const scalar trackTime)
     {
         // Leapfrog tracking part
 
-        scalar tEnd = (1.0 - stepFraction())*trackTime;
-        scalar dtMax = tEnd;
-
-        while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
+        while (td.keepParticle && !td.switchProcessor && stepFraction() < 1)
         {
-            // set the lagrangian time-step
-            scalar dt = min(dtMax, tEnd);
-
-            dt *= trackToFace(position() + dt*v_, td);
-
-            tEnd -= dt;
-            stepFraction() = 1.0 - tEnd/trackTime;
+            const scalar f = 1 - stepFraction();
+            trackToFace(f*trackTime*v_, f, td);
         }
     }
     else if (td.part() == 2)
@@ -205,7 +197,7 @@ void Foam::molecule::transformProperties(const tensor& T)
 
     rf_ = transform(T, rf_);
 
-    sitePositions_ = position_ + (T & (sitePositions_ - position_));
+    sitePositions_ = position() + (T & (sitePositions_ - position()));
 
     siteForces_ = T & siteForces_;
 }
@@ -226,7 +218,7 @@ void Foam::molecule::transformProperties(const vector& separation)
 
 void Foam::molecule::setSitePositions(const constantProperties& constProps)
 {
-    sitePositions_ = position_ + (Q_ & constProps.siteReferencePositions());
+    sitePositions_ = position() + (Q_ & constProps.siteReferencePositions());
 }
 
 
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
index 46ace15276..5fd058d473 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -238,7 +238,7 @@ public:
         inline molecule
         (
             const polyMesh& mesh,
-            const vector& position,
+            const barycentric& coordinates,
             const label celli,
             const label tetFacei,
             const label tetPti,
@@ -253,6 +253,24 @@ public:
             const label id
         );
 
+        //- Construct from a position and a cell, searching for the rest of the
+        //  required topology
+        inline molecule
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label celli,
+            const tensor& Q,
+            const vector& v,
+            const vector& a,
+            const vector& pi,
+            const vector& tau,
+            const vector& specialPosition,
+            const constantProperties& constProps,
+            const label special,
+            const label id
+        );
+
         //- Construct from Istream
         molecule
         (
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
index 5c3058d0d4..b1c9a354b4 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -220,7 +220,7 @@ inline Foam::molecule::constantProperties::constantProperties
 inline Foam::molecule::molecule
 (
     const polyMesh& mesh,
-    const vector& position,
+    const barycentric& coordinates,
     const label celli,
     const label tetFacei,
     const label tetPti,
@@ -236,7 +236,42 @@ inline Foam::molecule::molecule
 
 )
 :
-    particle(mesh, position, celli, tetFacei, tetPti),
+    particle(mesh, coordinates, celli, tetFacei, tetPti),
+    Q_(Q),
+    v_(v),
+    a_(a),
+    pi_(pi),
+    tau_(tau),
+    specialPosition_(specialPosition),
+    potentialEnergy_(0.0),
+    rf_(Zero),
+    special_(special),
+    id_(id),
+    siteForces_(constProps.nSites(), Zero),
+    sitePositions_(constProps.nSites())
+{
+    setSitePositions(constProps);
+}
+
+
+inline Foam::molecule::molecule
+(
+    const polyMesh& mesh,
+    const vector& position,
+    const label celli,
+    const tensor& Q,
+    const vector& v,
+    const vector& a,
+    const vector& pi,
+    const vector& tau,
+    const vector& specialPosition,
+    const constantProperties& constProps,
+    const label special,
+    const label id
+
+)
+:
+    particle(mesh, position, celli),
     Q_(Q),
     v_(v),
     a_(a),
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
index 5318c57867..3ff6175c44 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
@@ -746,17 +746,8 @@ void Foam::moleculeCloud::initialiseMolecules
                                 globalPosition
                             );
 
-                            label cell = -1;
-                            label tetFace = -1;
-                            label tetPt = -1;
-
-                            mesh_.findCellFacePt
-                            (
-                                globalPosition,
-                                cell,
-                                tetFace,
-                                tetPt
-                            );
+                            const label cell =
+                                mesh_.cellTree().findInside(globalPosition);
 
                             if (findIndex(zone, cell) != -1)
                             {
@@ -764,8 +755,6 @@ void Foam::moleculeCloud::initialiseMolecules
                                 (
                                     globalPosition,
                                     cell,
-                                    tetFace,
-                                    tetPt,
                                     id,
                                     tethered,
                                     temperature,
@@ -834,17 +823,11 @@ void Foam::moleculeCloud::initialiseMolecules
                                                 globalPosition
                                             );
 
-                                        label cell = -1;
-                                        label tetFace = -1;
-                                        label tetPt = -1;
-
-                                        mesh_.findCellFacePt
-                                        (
-                                            globalPosition,
-                                            cell,
-                                            tetFace,
-                                            tetPt
-                                        );
+                                        const label cell =
+                                            mesh_.cellTree().findInside
+                                            (
+                                                globalPosition
+                                            );
 
                                         if (findIndex(zone, cell) != -1)
                                         {
@@ -852,8 +835,6 @@ void Foam::moleculeCloud::initialiseMolecules
                                             (
                                                 globalPosition,
                                                 cell,
-                                                tetFace,
-                                                tetPt,
                                                 id,
                                                 tethered,
                                                 temperature,
@@ -913,17 +894,11 @@ void Foam::moleculeCloud::initialiseMolecules
                                                 globalPosition
                                             );
 
-                                        label cell = -1;
-                                        label tetFace = -1;
-                                        label tetPt = -1;
-
-                                        mesh_.findCellFacePt
-                                        (
-                                            globalPosition,
-                                            cell,
-                                            tetFace,
-                                            tetPt
-                                        );
+                                        const label cell =
+                                            mesh_.cellTree().findInside
+                                            (
+                                                globalPosition
+                                            );
 
                                         if (findIndex(zone, cell) != -1)
                                         {
@@ -931,8 +906,6 @@ void Foam::moleculeCloud::initialiseMolecules
                                             (
                                                 globalPosition,
                                                 cell,
-                                                tetFace,
-                                                tetPt,
                                                 id,
                                                 tethered,
                                                 temperature,
@@ -996,26 +969,12 @@ void Foam::moleculeCloud::createMolecule
 (
     const point& position,
     label cell,
-    label tetFace,
-    label tetPt,
     label id,
     bool tethered,
     scalar temperature,
     const vector& bulkVelocity
 )
 {
-    if (cell == -1)
-    {
-        mesh_.findCellFacePt(position, cell, tetFace, tetPt);
-    }
-
-    if (cell == -1)
-    {
-        FatalErrorInFunction
-            << "Position specified does not correspond to a mesh cell." << nl
-            << abort(FatalError);
-    }
-
     point specialPosition(Zero);
 
     label special = 0;
@@ -1068,8 +1027,6 @@ void Foam::moleculeCloud::createMolecule
             mesh_,
             position,
             cell,
-            tetFace,
-            tetPt,
             Q,
             v,
             Zero,
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
index 35a109e940..4dda921e89 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -114,8 +114,6 @@ private:
         (
             const point& position,
             label cell,
-            label tetFace,
-            label tetPt,
             label id,
             bool tethered,
             scalar temperature,
diff --git a/src/lagrangian/solidParticle/solidParticle.C b/src/lagrangian/solidParticle/solidParticle.C
index b94e99a5e1..ec189b4a14 100644
--- a/src/lagrangian/solidParticle/solidParticle.C
+++ b/src/lagrangian/solidParticle/solidParticle.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,34 +43,27 @@ bool Foam::solidParticle::move
     td.switchProcessor = false;
     td.keepParticle = true;
 
-    const polyBoundaryMesh& pbMesh = mesh_.boundaryMesh();
+    const polyBoundaryMesh& pbMesh = mesh().boundaryMesh();
 
-    scalar tEnd = (1.0 - stepFraction())*trackTime;
-    scalar dtMax = tEnd;
-
-    while (td.keepParticle && !td.switchProcessor && tEnd > SMALL)
+    while (td.keepParticle && !td.switchProcessor && stepFraction() < 1)
     {
         if (debug)
         {
-            Info<< "Time = " << mesh_.time().timeName()
+            Info<< "Time = " << mesh().time().timeName()
                 << " trackTime = " << trackTime
-                << " tEnd = " << tEnd
                 << " steptFraction() = " << stepFraction() << endl;
         }
 
-        // set the lagrangian time-step
-        scalar dt = min(dtMax, tEnd);
 
-        // remember which cell the parcel is in
-        // since this will change if a face is hit
-        label celli = cell();
+        const label celli = cell();
+        const scalar sfrac = stepFraction();
 
-        dt *= trackToFace(position() + dt*U_, td);
+        const scalar f = 1 - stepFraction();
+        trackToFace(f*trackTime*U_, f, td);
 
-        tEnd -= dt;
-        stepFraction() = 1.0 - tEnd/trackTime;
+        const scalar dt = (stepFraction() - sfrac)*trackTime;
 
-        cellPointWeight cpw(mesh_, position(), celli, face());
+        cellPointWeight cpw(mesh(), position(), celli, face());
         scalar rhoc = td.rhoInterp().interpolate(cpw);
         vector Uc = td.UInterp().interpolate(cpw);
         scalar nuc = td.nuInterp().interpolate(cpw);
@@ -90,7 +83,7 @@ bool Foam::solidParticle::move
 
         U_ = (U_ + dt*(Dc*Uc + (1.0 - rhoc/rhop)*td.g()))/(1.0 + dt*Dc);
 
-        if (onBoundary() && td.keepParticle)
+        if (onBoundaryFace() && td.keepParticle)
         {
             if (isA<processorPolyPatch>(pbMesh[patch(face())]))
             {
@@ -133,7 +126,7 @@ void Foam::solidParticle::hitWallPatch
     const tetIndices& tetIs
 )
 {
-    vector nw = tetIs.faceTri(mesh_).normal();
+    vector nw = tetIs.faceTri(mesh()).normal();
     nw /= mag(nw);
 
     scalar Un = U_ & nw;
diff --git a/src/lagrangian/solidParticle/solidParticle.H b/src/lagrangian/solidParticle/solidParticle.H
index f763abf339..bc03cd3b2e 100644
--- a/src/lagrangian/solidParticle/solidParticle.H
+++ b/src/lagrangian/solidParticle/solidParticle.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -133,7 +133,7 @@ public:
         inline solidParticle
         (
             const polyMesh& mesh,
-            const vector& position,
+            const barycentric& coordinates,
             const label celli,
             const label tetFacei,
             const label tetPti,
diff --git a/src/lagrangian/solidParticle/solidParticleI.H b/src/lagrangian/solidParticle/solidParticleI.H
index c1b6fdd418..a9a43591a3 100644
--- a/src/lagrangian/solidParticle/solidParticleI.H
+++ b/src/lagrangian/solidParticle/solidParticleI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,7 +45,7 @@ inline Foam::solidParticle::trackingData::trackingData
 inline Foam::solidParticle::solidParticle
 (
     const polyMesh& mesh,
-    const vector& position,
+    const barycentric& coordinates,
     const label celli,
     const label tetFacei,
     const label tetPti,
@@ -53,7 +53,7 @@ inline Foam::solidParticle::solidParticle
     const vector& U
 )
 :
-    particle(mesh, position, celli, tetFacei, tetPti),
+    particle(mesh, coordinates, celli, tetFacei, tetPti),
     d_(d),
     U_(U)
 {}
diff --git a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.H b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.H
index 98bf7fcc22..c5cd250bfe 100644
--- a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.H
+++ b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -219,22 +219,31 @@ public:
 
     // Constructors
 
-        //- Construct from owner, position, and cloud owner
+        //- Construct from mesh, coordinates and topology
         //  Other properties initialised as null
         inline SprayParcel
         (
             const polyMesh& mesh,
-            const vector& position,
+            const barycentric& coordinates,
             const label celli,
             const label tetFacei,
             const label tetPti
         );
 
+        //- Construct from a position and a cell, searching for the rest of the
+        //  required topology. Other properties are initialised as null.
+        inline SprayParcel
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label celli
+        );
+
         //- Construct from components
         inline SprayParcel
         (
             const polyMesh& mesh,
-            const vector& position,
+            const barycentric& coordinates,
             const label celli,
             const label tetFacei,
             const label tetPti,
diff --git a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcelI.H b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcelI.H
index d0f7391b82..44a8714e2c 100644
--- a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcelI.H
+++ b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcelI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -107,15 +107,15 @@ template<class ParcelType>
 inline Foam::SprayParcel<ParcelType>::SprayParcel
 (
     const polyMesh& mesh,
-    const vector& position,
+    const barycentric& coordinates,
     const label celli,
     const label tetFacei,
     const label tetPti
 )
 :
-    ParcelType(mesh, position, celli, tetFacei, tetPti),
+    ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
     d0_(this->d()),
-    position0_(position),
+    position0_(this->position()),
     sigma_(0.0),
     mu_(0.0),
     liquidCore_(0.0),
@@ -135,6 +135,31 @@ inline Foam::SprayParcel<ParcelType>::SprayParcel
 (
     const polyMesh& mesh,
     const vector& position,
+    const label celli
+)
+:
+    ParcelType(mesh, position, celli),
+    d0_(this->d()),
+    position0_(this->position()),
+    sigma_(0.0),
+    mu_(0.0),
+    liquidCore_(0.0),
+    KHindex_(0.0),
+    y_(0.0),
+    yDot_(0.0),
+    tc_(0.0),
+    ms_(0.0),
+    injector_(1.0),
+    tMom_(GREAT),
+    user_(0.0)
+{}
+
+
+template<class ParcelType>
+inline Foam::SprayParcel<ParcelType>::SprayParcel
+(
+    const polyMesh& mesh,
+    const barycentric& coordinates,
     const label celli,
     const label tetFacei,
     const label tetPti,
@@ -162,7 +187,7 @@ inline Foam::SprayParcel<ParcelType>::SprayParcel
     ParcelType
     (
         mesh,
-        position,
+        coordinates,
         celli,
         tetFacei,
         tetPti,
@@ -178,7 +203,7 @@ inline Foam::SprayParcel<ParcelType>::SprayParcel
         constProps
     ),
     d0_(d0),
-    position0_(position),
+    position0_(this->position()),
     sigma_(constProps.sigma0()),
     mu_(constProps.mu0()),
     liquidCore_(liquidCore),
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
index 5d69b69098..12de3f54e2 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -350,18 +350,7 @@ void Foam::meshRefinement::markFeatureCellLevel
     {
         const point& keepPoint = keepPoints[i];
 
-        label cellI = -1;
-        label tetFaceI = -1;
-        label tetPtI = -1;
-
-
-        // Force construction of search tree even if processor holds no
-        // cells
-        (void)mesh_.cellTree();
-        if (mesh_.nCells())
-        {
-            mesh_.findCellFacePt(keepPoint, cellI, tetFaceI, tetPtI);
-        }
+        const label celli = mesh_.cellTree().findInside(keepPoint);
 
         if (cellI != -1)
         {
@@ -404,10 +393,15 @@ void Foam::meshRefinement::markFeatureCellLevel
                             (
                                 mesh_,
                                 keepPoint,
+<<<<<<< HEAD
                                 cellI,
                                 tetFaceI,
                                 tetPtI,
                                 featureMesh.points()[pointI],   // endpos
+=======
+                                celli,
+                                featureMesh.points()[pointi],   // endpos
+>>>>>>> 371762757... Lagrangian: Rewrite of the particle tracking algorithm to function in
                                 featureLevel,                   // level
                                 featI,                          // featureMesh
                                 pointI,                         // end point
@@ -449,10 +443,15 @@ void Foam::meshRefinement::markFeatureCellLevel
                             (
                                 mesh_,
                                 keepPoint,
+<<<<<<< HEAD
                                 cellI,
                                 tetFaceI,
                                 tetPtI,
                                 featureMesh.points()[pointI],   // endpos
+=======
+                                celli,
+                                featureMesh.points()[pointi],   // endpos
+>>>>>>> 371762757... Lagrangian: Rewrite of the particle tracking algorithm to function in
                                 featureLevel,                   // level
                                 featI,                          // featureMesh
                                 pointI,                         // end point
@@ -545,8 +544,14 @@ void Foam::meshRefinement::markFeatureCellLevel
                 label otherPointI = e.otherVertex(pointI);
 
                 trackedParticle* tp(new trackedParticle(startTp));
+<<<<<<< HEAD
                 tp->end() = featureMesh.points()[otherPointI];
                 tp->j() = otherPointI;
+=======
+                tp->start() = tp->position();
+                tp->end() = featureMesh.points()[otherPointi];
+                tp->j() = otherPointi;
+>>>>>>> 371762757... Lagrangian: Rewrite of the particle tracking algorithm to function in
                 tp->k() = edgeI;
 
                 if (debug&meshRefinement::FEATURESEEDS)
@@ -605,8 +610,14 @@ void Foam::meshRefinement::markFeatureCellLevel
                     const edge& e = featureMesh.edges()[edgeI];
                     label otherPointI = e.otherVertex(pointI);
 
+<<<<<<< HEAD
                     tp.end() = featureMesh.points()[otherPointI];
                     tp.j() = otherPointI;
+=======
+                    tp.start() = tp.position();
+                    tp.end() = featureMesh.points()[otherPointi];
+                    tp.j() = otherPointi;
+>>>>>>> 371762757... Lagrangian: Rewrite of the particle tracking algorithm to function in
                     tp.k() = edgeI;
                     keepParticle = true;
                     break;
diff --git a/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.C b/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.C
index e4ba1c5409..71f44569be 100644
--- a/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.C
+++ b/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,7 +31,7 @@ License
 Foam::trackedParticle::trackedParticle
 (
     const polyMesh& mesh,
-    const vector& position,
+    const barycentric& coordinates,
     const label celli,
     const label tetFacei,
     const label tetPtI,
@@ -42,7 +42,30 @@ Foam::trackedParticle::trackedParticle
     const label k
 )
 :
-    particle(mesh, position, celli, tetFacei, tetPtI),
+    particle(mesh, coordinates, celli, tetFacei, tetPtI),
+    start_(position()),
+    end_(end),
+    level_(level),
+    i_(i),
+    j_(j),
+    k_(k)
+{}
+
+
+Foam::trackedParticle::trackedParticle
+(
+    const polyMesh& mesh,
+    const vector& position,
+    const label celli,
+    const point& end,
+    const label level,
+    const label i,
+    const label j,
+    const label k
+)
+:
+    particle(mesh, position, celli),
+    start_(this->position()),
     end_(end),
     level_(level),
     i_(i),
@@ -64,7 +87,7 @@ Foam::trackedParticle::trackedParticle
     {
         if (is.format() == IOstream::ASCII)
         {
-            is >> end_;
+            is >> start_ >> end_;
             level_ = readLabel(is);
             i_ = readLabel(is);
             j_ = readLabel(is);
@@ -74,8 +97,8 @@ Foam::trackedParticle::trackedParticle
         {
             is.read
             (
-                reinterpret_cast<char*>(&end_),
-                sizeof(end_) + sizeof(level_)
+                reinterpret_cast<char*>(&start_),
+                sizeof(start_) + sizeof(end_) + sizeof(level_)
               + sizeof(i_) + sizeof(j_) + sizeof(k_)
             );
         }
@@ -96,9 +119,8 @@ bool Foam::trackedParticle::move
     td.switchProcessor = false;
 
     scalar tEnd = (1.0 - stepFraction())*trackTime;
-    scalar dtMax = tEnd;
 
-    if (tEnd <= SMALL && onBoundary())
+    if (tEnd <= SMALL && onBoundaryFace())
     {
         // This is a hack to handle particles reaching their endpoint
         // on a processor boundary. If the endpoint is on a processor face
@@ -111,18 +133,13 @@ bool Foam::trackedParticle::move
     {
         td.keepParticle = true;
 
-        while (td.keepParticle && !td.switchProcessor && tEnd > SMALL)
+        while (td.keepParticle && !td.switchProcessor && stepFraction() < 1)
         {
-            // set the lagrangian time-step
-            scalar dt = min(dtMax, tEnd);
-
             // mark visited cell with max level.
             td.maxLevel_[cell()] = max(td.maxLevel_[cell()], level_);
 
-            dt *= trackToFace(end_, td);
-
-            tEnd -= dt;
-            stepFraction() = 1.0 - tEnd/trackTime;
+            const scalar f = 1 - stepFraction();
+            trackToFace(f*(end_ - start_), f, td);
         }
     }
 
@@ -248,6 +265,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const trackedParticle& p)
     if (os.format() == IOstream::ASCII)
     {
         os  << static_cast<const particle&>(p)
+            << token::SPACE << p.start_
             << token::SPACE << p.end_
             << token::SPACE << p.level_
             << token::SPACE << p.i_
@@ -259,8 +277,8 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const trackedParticle& p)
         os  << static_cast<const particle&>(p);
         os.write
         (
-            reinterpret_cast<const char*>(&p.end_),
-            sizeof(p.end_) + sizeof(p.level_)
+            reinterpret_cast<const char*>(&p.start_),
+            sizeof(p.start_) + sizeof(p.end_) + sizeof(p.level_)
           + sizeof(p.i_) + sizeof(p.j_) + sizeof(p.k_)
         );
     }
diff --git a/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.H b/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.H
index faf5519461..02a30b3a23 100644
--- a/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.H
+++ b/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -64,6 +64,9 @@ class trackedParticle
 {
     // Private data
 
+        //- Start point to track from
+        point start_;
+
         //- End point to track to
         point end_;
 
@@ -120,7 +123,7 @@ public:
         trackedParticle
         (
             const polyMesh& mesh,
-            const vector& position,
+            const barycentric& coordinates,
             const label celli,
             const label tetFacei,
             const label tetPtI,
@@ -131,6 +134,20 @@ public:
             const label k
         );
 
+        //- Construct from a position and a cell, searching for the rest of the
+        //  required topology
+        trackedParticle
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label celli,
+            const point& end,
+            const label level,
+            const label i,
+            const label j,
+            const label k
+        );
+
         //- Construct from Istream
         trackedParticle
         (
@@ -170,6 +187,12 @@ public:
 
     // Member Functions
 
+        //- Point to track from
+        point& start()
+        {
+            return start_;
+        }
+
         //- Point to track to
         point& end()
         {
diff --git a/src/parallel/reconstruct/reconstruct/reconstructLagrangianPositions.C b/src/parallel/reconstruct/reconstruct/reconstructLagrangianPositions.C
index 2136819b6c..c8e8be4c5e 100644
--- a/src/parallel/reconstruct/reconstruct/reconstructLagrangianPositions.C
+++ b/src/parallel/reconstruct/reconstruct/reconstructLagrangianPositions.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,9 +48,7 @@ void Foam::reconstructLagrangianPositions
     forAll(meshes, i)
     {
         const labelList& cellMap = cellProcAddressing[i];
-
-        // faceProcAddressing not required currently.
-        // const labelList& faceMap = faceProcAddressing[i];
+        const labelList& faceMap = faceProcAddressing[i];
 
         Cloud<passiveParticle> lpi(meshes[i], cloudName, false);
 
@@ -58,18 +56,21 @@ void Foam::reconstructLagrangianPositions
         {
             const passiveParticle& ppi = iter();
 
-            // // Inverting sign if necessary and subtracting 1 from
-            // // faceProcAddressing
-            // label mappedTetFace = mag(faceMap[ppi.tetFace()]) - 1;
+            const label mappedCell = cellMap[ppi.cell()];
+
+            // Inverting sign if necessary and subtracting 1 from
+            // faceProcAddressing
+            label mappedTetFace = mag(faceMap[ppi.tetFace()]) - 1;
 
             lagrangianPositions.append
             (
                 new passiveParticle
                 (
                     mesh,
-                    ppi.position(),
-                    cellMap[ppi.cell()],
-                    false
+                    ppi.coordinates(),
+                    mappedCell,
+                    mappedTetFace,
+                    ppi.procTetPt(mesh, mappedCell, mappedTetFace)
                 )
             );
         }
diff --git a/src/sampling/sampledSet/face/faceOnlySet.C b/src/sampling/sampledSet/face/faceOnlySet.C
index b7d217a35c..429d1c084b 100644
--- a/src/sampling/sampledSet/face/faceOnlySet.C
+++ b/src/sampling/sampledSet/face/faceOnlySet.C
@@ -56,13 +56,15 @@ bool Foam::faceOnlySet::trackToBoundary
 {
     particle::TrackingData<passiveParticleCloud> trackData(particleCloud);
 
-    const point& trackPt = singleParticle.position();
+    point trackPt = singleParticle.position();
 
     while(true)
     {
         point oldPoint = trackPt;
 
-        singleParticle.trackToFace(end_, trackData);
+        singleParticle.trackToFace(end_ - start_, 0, trackData);
+
+        trackPt = singleParticle.position();
 
         if (singleParticle.face() != -1 && mag(oldPoint - trackPt) > smallDist)
         {
@@ -78,7 +80,7 @@ bool Foam::faceOnlySet::trackToBoundary
             // End reached
             return false;
         }
-        else if (singleParticle.onBoundary())
+        else if (singleParticle.onBoundaryFace())
         {
             // Boundary reached
             return true;
diff --git a/src/sampling/sampledSet/polyLine/polyLineSet.C b/src/sampling/sampledSet/polyLine/polyLineSet.C
index 816945ec5b..358162e7a6 100644
--- a/src/sampling/sampledSet/polyLine/polyLineSet.C
+++ b/src/sampling/sampledSet/polyLine/polyLineSet.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,29 +56,16 @@ bool Foam::polyLineSet::trackToBoundary
 {
     particle::TrackingData<passiveParticleCloud> trackData(particleCloud);
 
-    // Alias
-    const point& trackPt = singleParticle.position();
-
     while (true)
     {
         // Local geometry info
         const vector offset = sampleCoords_[sampleI+1] - sampleCoords_[sampleI];
         const scalar smallDist = mag(tol*offset);
 
-        point oldPos = trackPt;
-        label facei = -1;
-        do
-        {
-            singleParticle.stepFraction() = 0;
-            singleParticle.track(sampleCoords_[sampleI+1], trackData);
-        }
-        while
-        (
-            !singleParticle.onBoundary()
-         && (mag(trackPt - oldPos) < smallDist)
-        );
+        singleParticle.track(offset, 0);
+        const point trackPt = singleParticle.position();
 
-        if (singleParticle.onBoundary())
+        if (singleParticle.onBoundaryFace())
         {
             //Info<< "trackToBoundary : reached boundary"
             //    << "  trackPt:" << trackPt << endl;
@@ -94,7 +81,7 @@ bool Foam::polyLineSet::trackToBoundary
                 //    << endl;
                 samplingPts.append(trackPt);
                 samplingCells.append(singleParticle.cell());
-                samplingFaces.append(facei);
+                samplingFaces.append(singleParticle.face());
 
                 // trackPt is at sampleI+1
                 samplingCurveDist.append(1.0*(sampleI+1));
diff --git a/src/sampling/sampledSet/sampledSet/sampledSet.C b/src/sampling/sampledSet/sampledSet/sampledSet.C
index 78694160f3..1ad1364f7e 100644
--- a/src/sampling/sampledSet/sampledSet/sampledSet.C
+++ b/src/sampling/sampledSet/sampledSet/sampledSet.C
@@ -209,17 +209,25 @@ Foam::point Foam::sampledSet::pushIn
     label tetPtI;
     mesh().findTetFacePt(celli, facePt, tetFacei, tetPtI);
 
+    // This is the tolerance that was defined as a static constant of the
+    // particle class. It is no longer used by particle, following the switch to
+    // barycentric tracking. The only place in which the tolerance is now used
+    // is here. I'm not sure what the purpose of this code is, but it probably
+    // wants removing. It is doing tet-searches for a particle position. This
+    // should almost certainly be left to the particle class.
+    const scalar trackingCorrectionTol = 1e-5;
+
     if (tetFacei == -1 || tetPtI == -1)
     {
         newPosition = facePt;
 
-        label trap(1.0/particle::trackingCorrectionTol + 1);
+        label trap(1.0/trackingCorrectionTol + 1);
 
         label iterNo = 0;
 
         do
         {
-            newPosition += particle::trackingCorrectionTol*(cC - facePt);
+            newPosition += trackingCorrectionTol*(cC - facePt);
 
             mesh().findTetFacePt
             (
diff --git a/src/sampling/sampledSet/uniform/uniformSet.C b/src/sampling/sampledSet/uniform/uniformSet.C
index c16502d214..c09403a71c 100644
--- a/src/sampling/sampledSet/uniform/uniformSet.C
+++ b/src/sampling/sampledSet/uniform/uniformSet.C
@@ -94,8 +94,7 @@ bool Foam::uniformSet::trackToBoundary
     const vector smallVec = tol*offset;
     const scalar smallDist = mag(smallVec);
 
-    // Alias
-    const point& trackPt = singleParticle.position();
+    point trackPt = singleParticle.position();
 
     particle::TrackingData<passiveParticleCloud> trackData(particleCloud);
 
@@ -153,32 +152,10 @@ bool Foam::uniformSet::trackToBoundary
                 << "  to:" << samplePt << endl;
         }
 
-        point oldPos = trackPt;
-        label facei = -1;
-        do
-        {
-            singleParticle.stepFraction() = 0;
-            singleParticle.track(samplePt, trackData);
+        singleParticle.track(samplePt - trackPt, 0);
+        trackPt = singleParticle.position();
 
-            if (debug)
-            {
-                Pout<< "Result of tracking "
-                    << "  trackPt:" << trackPt
-                    << "  trackCelli:" << singleParticle.cell()
-                    << "  trackFacei:" << singleParticle.face()
-                    << "  onBoundary:" << singleParticle.onBoundary()
-                    << "  samplePt:" << samplePt
-                    << "  smallDist:" << smallDist
-                    << endl;
-            }
-        }
-        while
-        (
-            !singleParticle.onBoundary()
-         && (mag(trackPt - oldPos) < smallDist)
-        );
-
-        if (singleParticle.onBoundary())
+        if (singleParticle.onBoundaryFace())
         {
             //Pout<< "trackToBoundary : reached boundary" << endl;
             if (mag(trackPt - samplePt) < smallDist)
@@ -188,7 +165,7 @@ bool Foam::uniformSet::trackToBoundary
                 // Reached samplePt on boundary
                 samplingPts.append(trackPt);
                 samplingCells.append(singleParticle.cell());
-                samplingFaces.append(facei);
+                samplingFaces.append(singleParticle.face());
                 samplingCurveDist.append(mag(trackPt - start_));
             }
 
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/controlDict b/tutorials/incompressible/simpleFoam/motorBike/system/controlDict
index a626f4d688..c767564e02 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/system/controlDict
+++ b/tutorials/incompressible/simpleFoam/motorBike/system/controlDict
@@ -47,7 +47,6 @@ runTimeModifiable true;
 functions
 {
     #include "streamLines"
-    #include "wallBoundedStreamLines"
     #include "cuttingPlane"
     #include "forceCoeffs"
     #include "ensightWrite"