diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C b/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
index aa9fd58b70..fdddddfbfe 100644
--- a/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
+++ b/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
@@ -270,7 +270,7 @@ void PDRkEpsilon::correct()
     }
 
     tmp<volTensorField> tgradU = fvc::grad(U_);
-    volScalarField G = 2*mut_*(tgradU() && dev(symm(tgradU())));
+    volScalarField G("RASModel::G", mut_*(tgradU() && dev(twoSymm(tgradU()))));
     tgradU.clear();
 
     // Update espsilon and G at the wall
diff --git a/applications/utilities/postProcessing/dataConversion/foamToTecplot360/foamToTecplot360.C b/applications/utilities/postProcessing/dataConversion/foamToTecplot360/foamToTecplot360.C
index 2d11a6538d..36d57be60e 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToTecplot360/foamToTecplot360.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToTecplot360/foamToTecplot360.C
@@ -1019,6 +1019,7 @@ int main(int argc, char *argv[])
             }
         }
         writer.writeEnd();
+        Info<< endl;
 
 
 
@@ -1067,7 +1068,7 @@ int main(int argc, char *argv[])
             writer.writeInit
             (
                 runTime.caseName(),
-                cellVarNames,
+                allVarNames,
                 patchFileName,
                 DataFileType_Full
             );
@@ -1163,6 +1164,8 @@ int main(int argc, char *argv[])
                         << nl << endl;
                 }
             }
+            writer.writeEnd();
+            Info<< endl;
         }
 
 
diff --git a/src/OpenFOAM/db/IOstreams/Pstreams/IPstream.C b/src/OpenFOAM/db/IOstreams/Pstreams/IPstream.C
index 8842768098..087028c3d2 100644
--- a/src/OpenFOAM/db/IOstreams/Pstreams/IPstream.C
+++ b/src/OpenFOAM/db/IOstreams/Pstreams/IPstream.C
@@ -38,7 +38,17 @@ Foam::IPstream::IPstream
 )
 :
     Pstream(commsType, bufSize),
-    UIPstream(commsType, fromProcNo, buf_, externalBufPosition_),
+    UIPstream
+    (
+        commsType,
+        fromProcNo,
+        buf_,
+        externalBufPosition_,
+        UPstream::msgType(),        // tag
+        false,                      // do not clear buf_ if at end
+        format,
+        version
+    ),
     externalBufPosition_(0)
 {}
 
diff --git a/src/OpenFOAM/db/IOstreams/Pstreams/UIPstream.C b/src/OpenFOAM/db/IOstreams/Pstreams/UIPstream.C
index f4c37633f3..b644942c73 100644
--- a/src/OpenFOAM/db/IOstreams/Pstreams/UIPstream.C
+++ b/src/OpenFOAM/db/IOstreams/Pstreams/UIPstream.C
@@ -77,6 +77,21 @@ inline void Foam::UIPstream::readFromBuffer
 }
 
 
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::UIPstream::~UIPstream()
+{
+    if (clearAtEnd_ && eof())
+    {
+        if (debug)
+        {
+            Pout<< "UIPstream::~UIPstream() : clearing externalBuf_ of size "
+                << externalBuf_.size() << endl;
+        }
+        externalBuf_.clearStorage();
+    }
+}
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::Istream& Foam::UIPstream::read(token& t)
diff --git a/src/OpenFOAM/db/IOstreams/Pstreams/UIPstream.H b/src/OpenFOAM/db/IOstreams/Pstreams/UIPstream.H
index 30e4e96119..b4f084baed 100644
--- a/src/OpenFOAM/db/IOstreams/Pstreams/UIPstream.H
+++ b/src/OpenFOAM/db/IOstreams/Pstreams/UIPstream.H
@@ -67,6 +67,8 @@ class UIPstream
 
         const int tag_;
 
+        const bool clearAtEnd_;
+
         label messageSize_;
 
 
@@ -96,6 +98,7 @@ public:
             DynamicList<char>& externalBuf,
             label& externalBufPosition,
             const int tag = UPstream::msgType(),
+            const bool clearAtEnd = false,   // destroy externalBuf if at end
             streamFormat format=BINARY,
             versionNumber version=currentVersion
         );
@@ -104,6 +107,11 @@ public:
         UIPstream(const int fromProcNo, PstreamBuffers&);
 
 
+    // Destructor
+
+        ~UIPstream();
+
+
     // Member functions
 
         // Inquiry
diff --git a/src/Pstream/dummy/UIPread.C b/src/Pstream/dummy/UIPread.C
index da9613af37..c21a7e77b5 100644
--- a/src/Pstream/dummy/UIPread.C
+++ b/src/Pstream/dummy/UIPread.C
@@ -38,6 +38,7 @@ Foam::UIPstream::UIPstream
     DynamicList<char>& externalBuf,
     label& externalBufPosition,
     const int tag,
+    const bool clearAtEnd,
     streamFormat format,
     versionNumber version
 )
@@ -48,6 +49,7 @@ Foam::UIPstream::UIPstream
     externalBuf_(externalBuf),
     externalBufPosition_(externalBufPosition),
     tag_(tag),
+    clearAtEnd_(clearAtEnd),
     messageSize_(0)
 {
     notImplemented
@@ -59,6 +61,7 @@ Foam::UIPstream::UIPstream
             "DynamicList<char>&,"
             "label&,"
             "const int tag,"
+            "const bool,"
             "streamFormat, versionNumber"
         ")"
     );
diff --git a/src/Pstream/mpi/UIPread.C b/src/Pstream/mpi/UIPread.C
index 9a70f456b9..28c3fe3e55 100644
--- a/src/Pstream/mpi/UIPread.C
+++ b/src/Pstream/mpi/UIPread.C
@@ -42,6 +42,7 @@ Foam::UIPstream::UIPstream
     DynamicList<char>& externalBuf,
     label& externalBufPosition,
     const int tag,
+    const bool clearAtEnd,
     streamFormat format,
     versionNumber version
 )
@@ -52,6 +53,7 @@ Foam::UIPstream::UIPstream
     externalBuf_(externalBuf),
     externalBufPosition_(externalBufPosition),
     tag_(tag),
+    clearAtEnd_(clearAtEnd),
     messageSize_(0)
 {
     setOpened();
@@ -125,6 +127,7 @@ Foam::UIPstream::UIPstream(const int fromProcNo, PstreamBuffers& buffers)
     externalBuf_(buffers.recvBuf_[fromProcNo]),
     externalBufPosition_(buffers.recvBufPos_[fromProcNo]),
     tag_(buffers.tag_),
+    clearAtEnd_(true),
     messageSize_(0)
 {
     if (commsType() != UPstream::scheduled && !buffers.finishedSendsCalled_)
diff --git a/src/fvMotionSolver/Make/files b/src/fvMotionSolver/Make/files
index f229ae322e..44046b57ca 100644
--- a/src/fvMotionSolver/Make/files
+++ b/src/fvMotionSolver/Make/files
@@ -1,6 +1,7 @@
 fvMotionSolvers/fvMotionSolver/fvMotionSolver.C
 fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C
 fvMotionSolvers/displacement/displacementFvMotionSolver/displacementFvMotionSolver.C
+fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionFvMotionSolver.C
 fvMotionSolvers/displacement/interpolation/displacementInterpolationFvMotionSolver.C
 fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C
 fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionFvMotionSolver.C
new file mode 100644
index 0000000000..f86bd0439e
--- /dev/null
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionFvMotionSolver.C
@@ -0,0 +1,567 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "displacementLayeredMotionFvMotionSolver.H"
+#include "addToRunTimeSelectionTable.H"
+#include "pointEdgeStructuredWalk.H"
+#include "pointFields.H"
+#include "PointEdgeWave.H"
+#include "syncTools.H"
+#include "interpolationTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(displacementLayeredMotionFvMotionSolver, 0);
+
+    addToRunTimeSelectionTable
+    (
+        fvMotionSolver,
+        displacementLayeredMotionFvMotionSolver,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::displacementLayeredMotionFvMotionSolver::calcZoneMask
+(
+    const label cellZoneI, 
+    PackedBoolList& isZonePoint,
+    PackedBoolList& isZoneEdge
+) const
+{
+    if (cellZoneI == -1)
+    {
+        isZonePoint.setSize(mesh().nPoints());
+        isZonePoint = 1;
+
+        isZoneEdge.setSize(mesh().nEdges());
+        isZoneEdge = 1;
+    }
+    else
+    {
+        const cellZone& cz = mesh().cellZones()[cellZoneI];
+
+        label nPoints = 0;
+        forAll(cz, i)
+        {
+            const labelList& cPoints = mesh().cellPoints(cz[i]);
+            forAll(cPoints, cPointI)
+            {
+                if (!isZonePoint[cPoints[cPointI]])
+                {
+                    isZonePoint[cPoints[cPointI]] = 1;
+                    nPoints++;
+                }
+            }
+        }
+        syncTools::syncPointList
+        (
+            mesh(),
+            isZonePoint,
+            orEqOp<unsigned int>(),
+            0
+        );
+
+
+        // Mark edge inside cellZone
+        label nEdges = 0;
+        forAll(cz, i)
+        {
+            const labelList& cEdges = mesh().cellEdges(cz[i]);
+            forAll(cEdges, cEdgeI)
+            {
+                if (!isZoneEdge[cEdges[cEdgeI]])
+                {
+                    isZoneEdge[cEdges[cEdgeI]] = 1;
+                    nEdges++;
+                }
+            }
+        }
+        syncTools::syncEdgeList
+        (
+            mesh(),
+            isZoneEdge,
+            orEqOp<unsigned int>(),
+            0
+        );
+
+
+        Info<< "On cellZone " << cz.name()
+            << " marked " << returnReduce(nPoints, sumOp<label>())
+            << " points and " << returnReduce(nEdges, sumOp<label>())
+            << " edges." << endl;
+    }
+}
+
+
+// Find distance to starting point
+void Foam::displacementLayeredMotionFvMotionSolver::walkStructured
+(
+    const label cellZoneI,
+    const PackedBoolList& isZonePoint,
+    const PackedBoolList& isZoneEdge,
+    const labelList& seedPoints,
+    const vectorField& seedData,
+    scalarField& distance,
+    vectorField& data
+) const
+{
+    const pointField& points = mesh().points();
+
+    List<pointEdgeStructuredWalk> seedInfo(seedPoints.size());
+
+    forAll(seedPoints, i)
+    {
+        seedInfo[i] = pointEdgeStructuredWalk
+        (
+            true,
+            points[seedPoints[i]],
+            0.0,
+            seedData[i]
+        );
+    }
+
+    // Current info on points
+    List<pointEdgeStructuredWalk> allPointInfo(mesh().nPoints());
+    // Mark points inside cellZone
+    forAll(isZonePoint, pointI)
+    {
+        if (isZonePoint[pointI])
+        {
+            allPointInfo[pointI].inZone() = true;
+        }
+    }
+
+    // Current info on edges
+    List<pointEdgeStructuredWalk> allEdgeInfo(mesh().nEdges());
+    // Mark edges inside cellZone
+    forAll(isZoneEdge, edgeI)
+    {
+        if (isZoneEdge[edgeI])
+        {
+            allEdgeInfo[edgeI].inZone() = true;
+        }
+    }
+
+    // Walk
+    PointEdgeWave<pointEdgeStructuredWalk> wallCalc
+    (
+        mesh(),
+        seedPoints,
+        seedInfo,
+
+        allPointInfo,
+        allEdgeInfo,
+        mesh().globalData().nTotalPoints()  // max iterations
+    );
+
+    // Extract distance and passive data
+    forAll(allPointInfo, pointI)
+    {
+        if (isZonePoint[pointI])
+        {
+            distance[pointI] = allPointInfo[pointI].dist();
+            data[pointI] = allPointInfo[pointI].data();
+        }
+    }
+}
+
+
+// Evaluate faceZone patch
+Foam::tmp<Foam::vectorField>
+Foam::displacementLayeredMotionFvMotionSolver::faceZoneEvaluate
+(
+    const faceZone& fz,
+    const labelList& meshPoints,
+    const dictionary& dict,
+    const PtrList<pointVectorField>& patchDisp,
+    const label patchI
+) const
+{
+    tmp<vectorField> tfld(new vectorField(meshPoints.size()));
+    vectorField& fld = tfld();
+
+    const word& type = dict.lookup("type");
+
+    if (type == "fixedValue")
+    {
+        fld = vectorField("value", dict, meshPoints.size());
+    }
+    else if (type == "timeVaryingUniformFixedValue")
+    {
+        interpolationTable<vector> timeSeries(dict);
+
+        fld = timeSeries(mesh().time().timeOutputValue());
+    }
+    else if (type == "slip")
+    {
+        if ((patchI % 2) != 1)
+        {
+            FatalIOErrorIn
+            (
+                "displacementLayeredMotionFvMotionSolver::faceZoneEvaluate(..)",
+                *this
+            )   << "slip can only be used on second faceZonePatch of pair."
+                << "FaceZone:" << fz.name()
+                << exit(FatalIOError);
+        }
+        // Use field set by previous bc
+        fld = vectorField(patchDisp[patchI-1], meshPoints);
+    }
+    else if (type == "follow")
+    {
+        // Only on boundary faces - follow boundary conditions
+        fld = vectorField(pointDisplacement_, meshPoints);
+    }
+    else
+    {
+        FatalIOErrorIn
+        (
+            "displacementLayeredMotionFvMotionSolver::faceZoneEvaluate(..)",
+            *this
+        )   << "Unknown faceZonePatch type " << type << " for faceZone "
+            << fz.name() << exit(FatalIOError);
+    }
+    return tfld;
+}
+
+
+void Foam::displacementLayeredMotionFvMotionSolver::cellZoneSolve
+(
+    const label cellZoneI,
+    const dictionary& zoneDict
+)
+{
+    PackedBoolList isZonePoint(mesh().nPoints());
+    PackedBoolList isZoneEdge(mesh().nEdges());
+    calcZoneMask(cellZoneI, isZonePoint, isZoneEdge);
+
+    const dictionary& patchesDict = zoneDict.subDict("boundaryField");
+
+    if (patchesDict.size() != 2)
+    {
+        FatalIOErrorIn
+        (
+            "displacementLayeredMotionFvMotionSolver::"
+            "correctBoundaryConditions(..)",
+            *this
+        )   << "Can only handle 2 faceZones (= patches) per cellZone. "
+            << " cellZone:" << cellZoneI
+            << " patches:" << patchesDict.toc()
+            << exit(FatalIOError);
+    }
+
+    PtrList<scalarField> patchDist(patchesDict.size());
+    PtrList<pointVectorField> patchDisp(patchesDict.size());
+
+    // Allocate the fields
+    label patchI = 0;
+    forAllConstIter(dictionary, patchesDict, patchIter)
+    {
+        const word& faceZoneName = patchIter().keyword();
+        label zoneI = mesh().faceZones().findZoneID(faceZoneName);
+        if (zoneI == -1)
+        {
+            FatalIOErrorIn
+            (
+                "displacementLayeredMotionFvMotionSolver::"
+                "correctBoundaryConditions(..)",
+                *this
+            )   << "Cannot find faceZone " << faceZoneName
+                << endl << "Valid zones are " << mesh().faceZones().names()
+                << exit(FatalIOError);
+        }
+
+        // Determine the points of the faceZone within the cellZone
+        const faceZone& fz = mesh().faceZones()[zoneI];
+
+        patchDist.set(patchI, new scalarField(mesh().nPoints()));
+        patchDisp.set
+        (
+            patchI,
+            new pointVectorField
+            (
+                IOobject
+                (
+                    mesh().cellZones()[cellZoneI].name() + "_" + fz.name(),
+                    mesh().time().timeName(),
+                    mesh(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                pointDisplacement_  // to inherit the boundary conditions
+            )
+        );
+
+        patchI++;
+    }
+
+
+
+    // 'correctBoundaryConditions'
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Loops over all the faceZones and walks their boundary values
+
+    // Make sure we can pick up bc values from field
+    pointDisplacement_.correctBoundaryConditions();
+
+    patchI = 0;
+    forAllConstIter(dictionary, patchesDict, patchIter)
+    {
+        const word& faceZoneName = patchIter().keyword();
+        const dictionary& faceZoneDict = patchIter().dict();
+
+        // Determine the points of the faceZone within the cellZone
+        label zoneI = mesh().faceZones().findZoneID(faceZoneName);
+        const faceZone& fz = mesh().faceZones()[zoneI];
+        const labelList& fzMeshPoints = fz().meshPoints();
+        DynamicList<label> meshPoints(fzMeshPoints.size());
+        forAll(fzMeshPoints, i)
+        {
+            if (isZonePoint[fzMeshPoints[i]])
+            {
+                meshPoints.append(fzMeshPoints[i]);
+            }
+        }
+
+        // Get initial value for all the faceZone points
+        tmp<vectorField> tseed = faceZoneEvaluate
+        (
+            fz,
+            meshPoints,
+            faceZoneDict,
+            patchDisp,
+            patchI
+        );
+
+
+Info<< "For cellZone:" << cellZoneI
+    << " for faceZone:" << fz.name() << " nPoints:" << tseed().size()
+    << " have patchField:"
+    << " max:" << gMax(tseed())
+    << " min:" << gMin(tseed())
+    << " avg:" << gAverage(tseed())
+    << endl;
+
+        // Set distance and transported value
+        walkStructured
+        (
+            cellZoneI,
+            isZonePoint,
+            isZoneEdge,
+
+            meshPoints,
+            tseed,
+            patchDist[patchI],
+            patchDisp[patchI]
+        );
+
+        // Implement real bc.
+        patchDisp[patchI].correctBoundaryConditions();
+
+
+//Info<< "Writing displacement for faceZone " << fz.name()
+//    << " to " << patchDisp[patchI].name() << endl;
+//patchDisp[patchI].write();
+
+//        // Copy into pointDisplacement for other fields to use
+//        forAll(isZonePoint, pointI)
+//        {
+//            if (isZonePoint[pointI])
+//            {
+//                pointDisplacement_[pointI] = patchDisp[patchI][pointI];
+//            }
+//        }
+//        pointDisplacement_.correctBoundaryConditions();
+
+
+        patchI++;
+    }
+
+
+    // Solve
+    // ~~~~~
+    // solving the interior is just interpolating
+
+//    // Get normalised distance
+//    pointScalarField distance
+//    (
+//        IOobject
+//        (
+//            "distance",
+//            mesh().time().timeName(),
+//            mesh(),
+//            IOobject::NO_READ,
+//            IOobject::NO_WRITE,
+//            false
+//        ),
+//        pointMesh::New(mesh()),
+//        dimensionedScalar("distance", dimLength, 0.0)
+//    );
+//    forAll(distance, pointI)
+//    {
+//        if (isZonePoint[pointI])
+//        {
+//            scalar d1 = patchDist[0][pointI];
+//            scalar d2 = patchDist[1][pointI];
+//            if (d1+d2 > SMALL)
+//            {
+//                scalar s = d1/(d1+d2);
+//                distance[pointI] = s;
+//            }
+//        }
+//    }
+//    Info<< "Writing distance pointScalarField to " << mesh().time().timeName()
+//        << endl;
+//    distance.write();
+
+    // Average
+    forAll(pointDisplacement_, pointI)
+    {
+        if (isZonePoint[pointI])
+        {
+            scalar d1 = patchDist[0][pointI];
+            scalar d2 = patchDist[1][pointI];
+
+            scalar s = d1/(d1+d2+VSMALL);
+
+            pointDisplacement_[pointI] =
+                (1-s)*patchDisp[0][pointI]
+              + s*patchDisp[1][pointI];
+        }
+    }
+    pointDisplacement_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::displacementLayeredMotionFvMotionSolver::
+displacementLayeredMotionFvMotionSolver
+(
+    const polyMesh& mesh,
+    Istream& is
+)
+:
+    displacementFvMotionSolver(mesh, is),
+    pointDisplacement_
+    (
+        IOobject
+        (
+            "pointDisplacement",
+            fvMesh_.time().timeName(),
+            fvMesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pointMesh::New(fvMesh_)
+    )
+{
+    pointDisplacement_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::displacementLayeredMotionFvMotionSolver::
+~displacementLayeredMotionFvMotionSolver()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::pointField>
+Foam::displacementLayeredMotionFvMotionSolver::curPoints() const
+{
+    tmp<pointField> tcurPoints
+    (
+        points0() + pointDisplacement_.internalField()
+    );
+
+    twoDCorrectPoints(tcurPoints());
+
+//    const pointField& pts = tcurPoints();
+//    forAll(pts, pointI)
+//    {
+//        Info<< "    from:" << mesh().points()[pointI]
+//            << " to:" << pts[pointI]
+//            << endl;
+//    }
+
+
+    return tcurPoints;
+}
+
+
+void Foam::displacementLayeredMotionFvMotionSolver::solve()
+{
+    const dictionary& ms = mesh().lookupObject<motionSolver>("dynamicMeshDict");
+    const dictionary& solverDict = ms.subDict(typeName + "Coeffs");
+
+    // Apply all regions (=cellZones)
+
+    const dictionary& regionDicts = solverDict.subDict("regions");
+    forAllConstIter(dictionary, regionDicts, regionIter)
+    {
+        const word& cellZoneName = regionIter().keyword();
+        const dictionary& regionDict = regionIter().dict();
+
+        label zoneI = mesh().cellZones().findZoneID(cellZoneName);
+
+        Info<< "solve : zone:" << cellZoneName << " index:" << zoneI
+            << endl;
+
+        if (zoneI == -1)
+        {
+            FatalIOErrorIn
+            (
+                "displacementLayeredMotionFvMotionSolver::solve(..)",
+                *this
+            )   << "Cannot find cellZone " << cellZoneName
+                << endl << "Valid zones are " << mesh().cellZones().names()
+                << exit(FatalIOError);
+        }
+
+        cellZoneSolve(zoneI, regionDict);
+    }
+}
+
+
+void Foam::displacementLayeredMotionFvMotionSolver::updateMesh
+(
+    const mapPolyMesh& mpm
+)
+{
+    displacementFvMotionSolver::updateMesh(mpm);
+}
+
+
+// ************************************************************************* //
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionFvMotionSolver.H b/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionFvMotionSolver.H
new file mode 100644
index 0000000000..96bae2c12f
--- /dev/null
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionFvMotionSolver.H
@@ -0,0 +1,187 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::displacementLayeredMotionFvMotionSolver
+
+Description
+    Mesh motion solver for an (multi-block) extruded fvMesh. Gets given the
+    structure of the mesh blocks and boundary conditions on these blocks.
+
+    Note: should not be an fvMotionSolver but just a motionSolver. Only here
+    so we can reuse displacementFvMotionSolver functionality (e.g. surface
+    following boundary conditions)
+
+
+    The displacementLayeredMotionCoeffs subdict of dynamicMeshDict specifies
+    per region (=cellZone) the boundary conditions on two opposing patches
+    (=faceZones). It then interpolates the boundary values using topological
+    walking so requires the cellZone to be a layered mesh.
+
+    The boundary conditions on faceZones are currently:
+
+    follow: the faceZone follows the overall mesh displacement.
+            Use this for faceZones on boundary faces (so it uses the
+            proper boundary conditions on the pointDisplacement).
+
+    fixedValue: fixed value.
+
+    timeVaryingUniformFixedValue: table-driven fixed value.
+
+    slip: the second of a pair of faceZones follows the tangential movement
+          specified by the first faceZone. (= removes the normal component).
+
+SourceFiles
+    displacementLayeredMotionFvMotionSolver.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef displacementLayeredMotionFvMotionSolver_H
+#define displacementLayeredMotionFvMotionSolver_H
+
+#include "displacementFvMotionSolver.H"
+#include "PackedBoolList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward class declarations
+
+/*---------------------------------------------------------------------------*\
+             Class displacementLayeredMotionFvMotionSolver Declaration
+\*---------------------------------------------------------------------------*/
+
+class displacementLayeredMotionFvMotionSolver
+:
+    public displacementFvMotionSolver
+{
+    // Private data
+
+        //- Point motion field
+        mutable pointVectorField pointDisplacement_;
+
+
+    // Private Member Functions
+
+        void calcZoneMask
+        (
+            const label cellZoneI, 
+            PackedBoolList& isZonePoint,
+            PackedBoolList& isZoneEdge
+        ) const;
+
+        void walkStructured
+        (
+            const label cellZoneI,
+            const PackedBoolList& isZonePoint,
+            const PackedBoolList& isZoneEdge,
+            const labelList& seedPoints,
+            const vectorField& seedData,
+            scalarField& distance,
+            vectorField& data
+        ) const;
+
+        tmp<vectorField> faceZoneEvaluate
+        (
+            const faceZone& fz,
+            const labelList& meshPoints,
+            const dictionary& dict,
+            const PtrList<pointVectorField>& patchDisp,
+            const label patchI
+        ) const;
+
+        void cellZoneSolve
+        (
+            const label cellZoneI,
+            const dictionary& zoneDict
+        );
+
+
+        //- Disallow default bitwise copy construct
+        displacementLayeredMotionFvMotionSolver
+        (
+            const displacementLayeredMotionFvMotionSolver&
+        );
+
+        //- Disallow default bitwise assignment
+        void operator=(const displacementLayeredMotionFvMotionSolver&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("displacementLayeredMotion");
+
+
+    // Constructors
+
+        //- Construct from polyMesh and data stream
+        displacementLayeredMotionFvMotionSolver
+        (
+            const polyMesh&,
+            Istream& msDataUnused
+        );
+
+
+    // Destructor
+
+        ~displacementLayeredMotionFvMotionSolver();
+
+
+    // Member Functions
+
+        //- Return reference to the point motion displacement field
+        pointVectorField& pointDisplacement()
+        {
+            return pointDisplacement_;
+        }
+
+        //- Return const reference to the point motion displacement field
+        const pointVectorField& pointDisplacement() const
+        {
+            return pointDisplacement_;
+        }
+
+        //- Return point location obtained from the current motion field
+        virtual tmp<pointField> curPoints() const;
+
+        //- Solve for motion
+        virtual void solve();
+
+        //- Update topology
+        virtual void updateMesh(const mapPolyMesh&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.C b/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.C
new file mode 100644
index 0000000000..ec56628f02
--- /dev/null
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.C
@@ -0,0 +1,54 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "pointEdgeStructuredWalk.H"
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Foam::Ostream& os,
+    const Foam::pointEdgeStructuredWalk& wDist
+)
+{
+    return os
+        << wDist.inZone_ << wDist.previousPoint_
+        << wDist.dist_ << wDist.data_;
+}
+
+Foam::Istream& Foam::operator>>
+(
+    Foam::Istream& is,
+    Foam::pointEdgeStructuredWalk& wDist
+)
+{
+    return is
+        >> wDist.inZone_ >> wDist.previousPoint_
+        >> wDist.dist_ >> wDist.data_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.H b/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.H
new file mode 100644
index 0000000000..83454236e6
--- /dev/null
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.H
@@ -0,0 +1,224 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::pointEdgeStructuredWalk
+
+Description
+    Determines length of string of edges walked to point.
+
+SourceFiles
+    pointEdgeStructuredWalkI.H
+    pointEdgeStructuredWalk.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pointEdgeStructuredWalk_H
+#define pointEdgeStructuredWalk_H
+
+#include "point.H"
+#include "label.H"
+#include "scalar.H"
+#include "tensor.H"
+#include "pTraits.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class polyPatch;
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class pointEdgeStructuredWalk Declaration
+\*---------------------------------------------------------------------------*/
+
+class pointEdgeStructuredWalk
+{
+    // Private data
+
+        //- Is this point/edge inside the walk zone
+        bool inZone_;
+
+        //- Previuos point
+        point previousPoint_;
+
+        //- Sum of distance
+        scalar dist_;
+
+        //- Passive data
+        vector data_;
+
+    // Private Member Functions
+
+        //- Evaluate distance to point.
+        inline bool update
+        (
+            const point&,
+            const pointEdgeStructuredWalk& w2,
+            const scalar tol
+        );
+
+        //- Combine current with w2. Update distSqr, origin if w2 has smaller
+        //  quantities and returns true.
+        inline bool update
+        (
+            const pointEdgeStructuredWalk& w2,
+            const scalar tol
+        );
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        inline pointEdgeStructuredWalk();
+
+        //- Construct from components
+        inline pointEdgeStructuredWalk
+        (
+            const bool,
+            const point&,
+            const scalar,
+            const vector&
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            inline bool inZone() const;
+
+            inline bool& inZone();
+
+            inline const point& previousPoint() const;
+
+            inline scalar dist() const;
+
+            inline const vector& data() const;
+
+
+        // Needed by meshWave
+
+            //- Check whether still contains original (invalid) value.
+            inline bool valid() const;
+
+            //- Check for identical geometrical data. Used for cyclics checking.
+            inline bool sameGeometry
+            (
+                const pointEdgeStructuredWalk&,
+                const scalar tol
+            ) const;
+
+            //- Convert origin to relative vector to leaving point
+            //  (= point coordinate)
+            inline void leaveDomain
+            (
+                const polyPatch& patch,
+                const label patchPointI,
+                const point& pos
+            );
+
+            //- Convert relative origin to absolute by adding entering point
+            inline void enterDomain
+            (
+                const polyPatch& patch,
+                const label patchPointI,
+                const point& pos
+            );
+
+            //- Apply rotation matrix to origin
+            inline void transform(const tensor& rotTensor);
+
+            //- Influence of edge on point
+            inline bool updatePoint
+            (
+                const polyMesh& mesh,
+                const label pointI,
+                const label edgeI,
+                const pointEdgeStructuredWalk& edgeInfo,
+                const scalar tol
+            );
+
+            //- Influence of different value on same point.
+            //  Merge new and old info.
+            inline bool updatePoint
+            (
+                const polyMesh& mesh,
+                const label pointI,
+                const pointEdgeStructuredWalk& newPointInfo,
+                const scalar tol
+            );
+
+            //- Influence of different value on same point.
+            //  No information about current position whatsoever.
+            inline bool updatePoint
+            (
+                const pointEdgeStructuredWalk& newPointInfo,
+                const scalar tol
+            );
+
+            //- Influence of point on edge.
+            inline bool updateEdge
+            (
+                const polyMesh& mesh,
+                const label edgeI,
+                const label pointI,
+                const pointEdgeStructuredWalk& pointInfo,
+                const scalar tol
+            );
+
+
+    // Member Operators
+
+        //Note: Used to determine whether to call update.
+        inline bool operator==(const pointEdgeStructuredWalk&) const;
+
+        inline bool operator!=(const pointEdgeStructuredWalk&) const;
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<<(Ostream&, const pointEdgeStructuredWalk&);
+        friend Istream& operator>>(Istream&, pointEdgeStructuredWalk&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "pointEdgeStructuredWalkI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/pointEdgeStructuredWalkI.H b/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/pointEdgeStructuredWalkI.H
new file mode 100644
index 0000000000..c6cc23f99f
--- /dev/null
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/pointEdgeStructuredWalkI.H
@@ -0,0 +1,302 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "polyMesh.H"
+#include "transform.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Update this with w2.
+inline bool Foam::pointEdgeStructuredWalk::update
+(
+    const point& pt,
+    const pointEdgeStructuredWalk& w2,
+    const scalar tol
+)
+{
+    if (!valid())
+    {
+        // current not yet set. Walked from w2 to here (=pt)
+        dist_ = w2.dist_ + mag(pt-w2.previousPoint_);
+        previousPoint_ = pt;
+        data_ = w2.data_;
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+    
+// Update this with w2 (information on same point)
+inline bool Foam::pointEdgeStructuredWalk::update
+(
+    const pointEdgeStructuredWalk& w2,
+    const scalar tol
+)
+{
+    if (!valid())
+    {
+        // current not yet set so use neighbour
+        dist_ = w2.dist_;
+        previousPoint_ = w2.previousPoint_;
+        data_ = w2.data_;
+
+        return true;
+    }
+    else
+    {
+        // already nearer to pt
+        return false;
+    }
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Null constructor
+inline Foam::pointEdgeStructuredWalk::pointEdgeStructuredWalk()
+:
+    inZone_(false),
+    previousPoint_(vector::max),
+    dist_(0),
+    data_(vector::zero)
+{}
+
+
+// Construct from origin, distance
+inline Foam::pointEdgeStructuredWalk::pointEdgeStructuredWalk
+(
+    const bool inZone,
+    const point& previousPoint, 
+    const scalar dist,
+    const vector& data
+)
+:
+    inZone_(inZone),
+    previousPoint_(previousPoint),
+    dist_(dist),
+    data_(data)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline bool Foam::pointEdgeStructuredWalk::inZone() const
+{
+    return inZone_;
+}
+
+
+inline bool& Foam::pointEdgeStructuredWalk::inZone()
+{
+    return inZone_;
+}
+
+
+inline const Foam::point& Foam::pointEdgeStructuredWalk::previousPoint() const
+{
+    return previousPoint_;
+}
+
+
+inline Foam::scalar Foam::pointEdgeStructuredWalk::dist() const
+{
+    return dist_;
+}
+
+
+inline const Foam::vector& Foam::pointEdgeStructuredWalk::data() const
+{
+    return data_;
+}
+
+
+inline bool Foam::pointEdgeStructuredWalk::valid() const
+{
+    return previousPoint_ != vector::max;
+}
+
+
+// Checks for cyclic points
+inline bool Foam::pointEdgeStructuredWalk::sameGeometry
+(
+    const pointEdgeStructuredWalk& w2,
+    const scalar tol
+) const
+{
+    scalar diff = Foam::mag(dist() - w2.dist());
+
+    if (diff < SMALL)
+    {
+        return true;
+    }
+    else
+    {
+        if ((dist() > SMALL) && ((diff/dist()) < tol))
+        {
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    }
+}
+
+
+inline void Foam::pointEdgeStructuredWalk::leaveDomain
+(
+    const polyPatch& patch,
+    const label patchPointI,
+    const point& coord
+)
+{
+    previousPoint_ -= coord;
+}
+
+
+inline void Foam::pointEdgeStructuredWalk::transform(const tensor& rotTensor)
+{
+    previousPoint_ = Foam::transform(rotTensor, previousPoint_);
+}
+
+
+// Update absolute geometric quantities. Note that distance (dist_)
+// is not affected by leaving/entering domain.
+inline void Foam::pointEdgeStructuredWalk::enterDomain
+(
+    const polyPatch& patch,
+    const label patchPointI,
+    const point& coord
+)
+{
+    // back to absolute form
+    previousPoint_ += coord;
+}
+
+
+// Update this with information from connected edge
+inline bool Foam::pointEdgeStructuredWalk::updatePoint
+(
+    const polyMesh& mesh,
+    const label pointI,
+    const label edgeI,
+    const pointEdgeStructuredWalk& edgeInfo,
+    const scalar tol
+)
+{
+    if (inZone_)
+    {
+        return update(mesh.points()[pointI], edgeInfo, tol);
+    }
+    else
+    {
+        return false;
+    }
+}    
+
+
+// Update this with new information on same point
+inline bool Foam::pointEdgeStructuredWalk::updatePoint
+(
+    const polyMesh& mesh,
+    const label pointI,
+    const pointEdgeStructuredWalk& newPointInfo,
+    const scalar tol
+)
+{
+    if (inZone_)
+    {
+        return update(mesh.points()[pointI], newPointInfo, tol);
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// Update this with new information on same point. No extra information.
+inline bool Foam::pointEdgeStructuredWalk::updatePoint
+(
+    const pointEdgeStructuredWalk& newPointInfo,
+    const scalar tol
+)
+{
+    return update(newPointInfo, tol);
+}    
+
+
+// Update this with information from connected point
+inline bool Foam::pointEdgeStructuredWalk::updateEdge
+(
+    const polyMesh& mesh,
+    const label edgeI,
+    const label pointI,
+    const pointEdgeStructuredWalk& pointInfo,
+    const scalar tol
+)
+{
+    if (inZone_)
+    {
+        const pointField& points = mesh.points();
+
+        const edge& e = mesh.edges()[edgeI];
+
+        const point edgeMid(0.5*(points[e[0]] + points[e[1]]));
+
+        return update(edgeMid, pointInfo, tol);
+    }
+    else
+    {
+        return false;
+    }
+}    
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline bool Foam::pointEdgeStructuredWalk::operator==
+(
+    const Foam::pointEdgeStructuredWalk& rhs
+) const
+{
+    return previousPoint() == rhs.previousPoint();
+}
+
+
+inline bool Foam::pointEdgeStructuredWalk::operator!=
+(
+    const Foam::pointEdgeStructuredWalk& rhs
+) const
+{
+    return !(*this == rhs);
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
index 92e8694874..894fc13d03 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
@@ -30,8 +30,6 @@ License
 #include "volPointInterpolation.H"
 #include "addToRunTimeSelectionTable.H"
 #include "fvMesh.H"
-#include "isoSurface.H"
-// #include "isoSurfaceCell.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -244,6 +242,14 @@ void Foam::distanceSurface::createGeometry()
             distance_,
             regularise_
         )
+        //new isoSurfaceCell
+        //(
+        //    fvm,
+        //    cellDistance,
+        //    pointDistance_,
+        //    distance_,
+        //    regularise_
+        //)
     );
 
     if (debug)
diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.H b/src/sampling/sampledSurface/distanceSurface/distanceSurface.H
index 9fd2ffd9f0..3e25b7c8a0 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurface.H
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.H
@@ -85,6 +85,7 @@ class distanceSurface
         scalarField pointDistance_;
 
         //- Constructed iso surface
+        //autoPtr<isoSurfaceCell> isoSurfPtr_;
         autoPtr<isoSurface> isoSurfPtr_;
 
         //- triangles converted to faceList
@@ -169,6 +170,7 @@ public:
         }
 
 
+        //const isoSurfaceCell& surface() const
         const isoSurface& surface() const
         {
             return isoSurfPtr_();
diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C b/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C
index ea54214fb6..df5e1befb9 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C
@@ -25,7 +25,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "distanceSurface.H"
-#include "isoSurface.H"
 #include "volFieldsFwd.H"
 #include "pointFields.H"
 #include "volPointInterpolation.H"
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.C b/src/sampling/sampledSurface/isoSurface/isoSurface.C
index 8bdce180e6..a3888f856e 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.C
@@ -2054,61 +2054,68 @@ Foam::isoSurface::isoSurface
     }
 
 
-    //if (false)
-    //{
-    List<FixedList<label, 3> > faceEdges;
-    labelList edgeFace0, edgeFace1;
-    Map<labelList> edgeFacesRest;
-
-
-    while (true)
+    if (false)
     {
-        // Calculate addressing
-        calcAddressing(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
+        List<FixedList<label, 3> > faceEdges;
+        labelList edgeFace0, edgeFace1;
+        Map<labelList> edgeFacesRest;
 
-        // See if any dangling triangles
-        boolList keepTriangles;
-        label nDangling = markDanglingTriangles
-        (
-            faceEdges,
-            edgeFace0,
-            edgeFace1,
-            edgeFacesRest,
-            keepTriangles
-        );
 
-        if (debug)
+        while (true)
         {
-            Pout<< "isoSurface : detected " << nDangling
-                << " dangling triangles." << endl;
-        }
-
-        if (nDangling == 0)
-        {
-            break;
-        }
-
-        // Create face map (new to old)
-        labelList subsetTriMap(findIndices(keepTriangles, true));
-
-        labelList subsetPointMap;
-        labelList reversePointMap;
-        triSurface::operator=
-        (
-            subsetMesh
+            // Calculate addressing
+            calcAddressing
             (
                 *this,
-                subsetTriMap,
-                reversePointMap,
-                subsetPointMap
-            )
-        );
-        meshCells_ = labelField(meshCells_, subsetTriMap);
-        inplaceRenumber(reversePointMap, triPointMergeMap_);
-    }
+                faceEdges,
+                edgeFace0,
+                edgeFace1,
+                edgeFacesRest
+            );
 
-    orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
-    //}
+            // See if any dangling triangles
+            boolList keepTriangles;
+            label nDangling = markDanglingTriangles
+            (
+                faceEdges,
+                edgeFace0,
+                edgeFace1,
+                edgeFacesRest,
+                keepTriangles
+            );
+
+            if (debug)
+            {
+                Pout<< "isoSurface : detected " << nDangling
+                    << " dangling triangles." << endl;
+            }
+
+            if (nDangling == 0)
+            {
+                break;
+            }
+
+            // Create face map (new to old)
+            labelList subsetTriMap(findIndices(keepTriangles, true));
+
+            labelList subsetPointMap;
+            labelList reversePointMap;
+            triSurface::operator=
+            (
+                subsetMesh
+                (
+                    *this,
+                    subsetTriMap,
+                    reversePointMap,
+                    subsetPointMap
+                )
+            );
+            meshCells_ = labelField(meshCells_, subsetTriMap);
+            inplaceRenumber(reversePointMap, triPointMergeMap_);
+        }
+
+        orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
+    }
 
 
     if (debug)
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
index 542d31a376..d534603b0f 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
@@ -1398,6 +1398,8 @@ Foam::isoSurfaceCell::isoSurfaceCell
 )
 :
     mesh_(mesh),
+    cVals_(cVals),
+    pVals_(pVals),
     iso_(iso),
     mergeDistance_(mergeTol*mesh.bounds().mag())
 {
@@ -1563,59 +1565,68 @@ Foam::isoSurfaceCell::isoSurfaceCell
     }
 
 
-    List<FixedList<label, 3> > faceEdges;
-    labelList edgeFace0, edgeFace1;
-    Map<labelList> edgeFacesRest;
-
-
-    while (true)
+    if (false)
     {
-        // Calculate addressing
-        calcAddressing(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
+        List<FixedList<label, 3> > faceEdges;
+        labelList edgeFace0, edgeFace1;
+        Map<labelList> edgeFacesRest;
 
-        // See if any dangling triangles
-        boolList keepTriangles;
-        label nDangling = markDanglingTriangles
-        (
-            faceEdges,
-            edgeFace0,
-            edgeFace1,
-            edgeFacesRest,
-            keepTriangles
-        );
 
-        if (debug)
+        while (true)
         {
-            Pout<< "isoSurfaceCell : detected " << nDangling
-                << " dangling triangles." << endl;
-        }
-
-        if (nDangling == 0)
-        {
-            break;
-        }
-
-        // Create face map (new to old)
-        labelList subsetTriMap(findIndices(keepTriangles, true));
-
-        labelList subsetPointMap;
-        labelList reversePointMap;
-        triSurface::operator=
-        (
-            subsetMesh
+            // Calculate addressing
+            calcAddressing
             (
                 *this,
-                subsetTriMap,
-                reversePointMap,
-                subsetPointMap
-            )
-        );
-        meshCells_ = labelField(meshCells_, subsetTriMap);
-        inplaceRenumber(reversePointMap, triPointMergeMap_);
+                faceEdges,
+                edgeFace0,
+                edgeFace1,
+                edgeFacesRest
+            );
+
+            // See if any dangling triangles
+            boolList keepTriangles;
+            label nDangling = markDanglingTriangles
+            (
+                faceEdges,
+                edgeFace0,
+                edgeFace1,
+                edgeFacesRest,
+                keepTriangles
+            );
+
+            if (debug)
+            {
+                Pout<< "isoSurfaceCell : detected " << nDangling
+                    << " dangling triangles." << endl;
+            }
+
+            if (nDangling == 0)
+            {
+                break;
+            }
+
+            // Create face map (new to old)
+            labelList subsetTriMap(findIndices(keepTriangles, true));
+
+            labelList subsetPointMap;
+            labelList reversePointMap;
+            triSurface::operator=
+            (
+                subsetMesh
+                (
+                    *this,
+                    subsetTriMap,
+                    reversePointMap,
+                    subsetPointMap
+                )
+            );
+            meshCells_ = labelField(meshCells_, subsetTriMap);
+            inplaceRenumber(reversePointMap, triPointMergeMap_);
+        }
+
+        orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
     }
-
-    orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
-
     //combineCellTriangles();
 }
 
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
index ea96b42b3a..ecb41ef7bc 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
@@ -85,6 +85,10 @@ class isoSurfaceCell
         //- Reference to mesh
         const polyMesh& mesh_;
 
+        const scalarField& cVals_;
+
+        const scalarField& pVals_;
+
         //- isoSurfaceCell value
         const scalar iso_;
 
@@ -352,6 +356,14 @@ public:
             const Field<Type>& cCoords,
             const Field<Type>& pCoords
         ) const;
+
+        //- Interpolates cCoords,pCoords.
+        template <class Type>
+        tmp<Field<Type> > interpolate
+        (
+            const Field<Type>& cCoords,
+            const Field<Type>& pCoords
+        ) const;
 };
 
 
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
index ca6ef109fa..16e48399ec 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
@@ -381,4 +381,56 @@ Foam::isoSurfaceCell::interpolate
 }
 
 
+template <class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::isoSurfaceCell::interpolate
+(
+    const Field<Type>& cCoords,
+    const Field<Type>& pCoords
+) const
+{
+    DynamicList<Type> triPoints(nCutCells_);
+    DynamicList<label> triMeshCells(nCutCells_);
+
+    // Dummy snap data
+    DynamicList<Type> snappedPoints;
+    labelList snappedCc(mesh_.nCells(), -1);
+    labelList snappedPoint(mesh_.nPoints(), -1);
+
+
+    generateTriPoints
+    (
+        cVals_,
+        pVals_,
+
+        cCoords,
+        pCoords,
+
+        snappedPoints,
+        snappedCc,
+        snappedPoint,
+
+        triPoints,
+        triMeshCells
+    );
+
+
+    // One value per point
+    tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
+    Field<Type>& values = tvalues();
+
+    forAll(triPoints, i)
+    {
+        label mergedPointI = triPointMergeMap_[i];
+
+        if (mergedPointI >= 0)
+        {
+            values[mergedPointI] = triPoints[i];
+        }
+    }
+
+    return tvalues;
+}
+
+
 // ************************************************************************* //
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
index 5e369ed533..0a8724bbd1 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
@@ -30,7 +30,6 @@ License
 #include "volPointInterpolation.H"
 #include "addToRunTimeSelectionTable.H"
 #include "fvMesh.H"
-#include "isoSurface.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -173,7 +172,7 @@ void Foam::sampledCuttingPlane::createGeometry()
 
                 forAll(fld, i)
                 {
-                    fld[i] =  (cc[i] - plane_.refPoint()) & plane_.normal();
+                    fld[i] = (cc[i] - plane_.refPoint()) & plane_.normal();
                 }
             }
         }
@@ -231,6 +230,14 @@ void Foam::sampledCuttingPlane::createGeometry()
             0.0,
             regularise_
         )
+        //new isoSurfaceCell
+        //(
+        //    fvm,
+        //    cellDistance,
+        //    pointDistance_,
+        //    0.0,
+        //    regularise_
+        //)
     );
 
     if (debug)
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
index 54e0f6627c..03eb4cfdab 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
@@ -38,6 +38,7 @@ SourceFiles
 
 #include "sampledSurface.H"
 #include "isoSurface.H"
+//#include "isoSurfaceCell.H"
 #include "plane.H"
 #include "ZoneIDs.H"
 #include "fvMeshSubset.H"
@@ -89,6 +90,7 @@ class sampledCuttingPlane
         scalarField pointDistance_;
 
         //- Constructed iso surface
+        //autoPtr<isoSurfaceCell> isoSurfPtr_;
         autoPtr<isoSurface> isoSurfPtr_;
 
         //- triangles converted to faceList
@@ -173,6 +175,7 @@ public:
         }
 
 
+        //const isoSurfaceCell& surface() const
         const isoSurface& surface() const
         {
             return isoSurfPtr_();
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C
index ea0ed6bf47..f594e1a78c 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C
@@ -26,7 +26,6 @@ License
 
 #include "volPointInterpolation.H"
 #include "sampledCuttingPlane.H"
-#include "isoSurface.H"
 #include "volFieldsFwd.H"
 #include "pointFields.H"
 #include "volPointInterpolation.H"
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
index 7e9274616e..9e35ee1edc 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
@@ -25,11 +25,10 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "sampledTriSurfaceMesh.H"
-#include "dictionary.H"
-#include "polyMesh.H"
-#include "polyPatch.H"
-#include "volFields.H"
+#include "treeDataPoint.H"
 #include "meshSearch.H"
+#include "Tuple2.H"
+#include "globalIndex.H"
 
 #include "addToRunTimeSelectionTable.H"
 
@@ -44,6 +43,25 @@ namespace Foam
         sampledTriSurfaceMesh,
         word
     );
+
+    //- Private class for finding nearest
+    //  - global index
+    //  - sqr(distance)
+    typedef Tuple2<scalar, label> nearInfo;
+
+    class nearestEqOp
+    {
+
+    public:
+
+        void operator()(nearInfo& x, const nearInfo& y) const
+        {
+            if (y.first() < x.first())
+            {
+                x = y;
+            }
+        }
+    };
 }
 
 
@@ -63,16 +81,16 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
         (
             name,
             mesh.time().constant(), // instance
-            "triSurface",       // local
-            mesh,               // registry
+            "triSurface",           // local
+            mesh,                   // registry
             IOobject::MUST_READ,
-            IOobject::NO_WRITE
+            IOobject::NO_WRITE,
+            false
         )
     ),
     needsUpdate_(true),
     cellLabels_(0),
-    pointMap_(0),
-    faceMap_(0)
+    pointToFace_(0)
 {}
 
 
@@ -90,16 +108,16 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
         (
             dict.lookup("surface"),
             mesh.time().constant(), // instance
-            "triSurface",       // local
-            mesh,               // registry
+            "triSurface",           // local
+            mesh,                   // registry
             IOobject::MUST_READ,
-            IOobject::NO_WRITE
+            IOobject::NO_WRITE,
+            false
         )
     ),
     needsUpdate_(true),
     cellLabels_(0),
-    pointMap_(0),
-    faceMap_(0)
+    pointToFace_(0)
 {}
 
 
@@ -128,8 +146,7 @@ bool Foam::sampledTriSurfaceMesh::expire()
     sampledSurface::clearGeom();
     MeshStorage::clear();
     cellLabels_.clear();
-    pointMap_.clear();
-    faceMap_.clear();
+    pointToFace_.clear();
 
     needsUpdate_ = true;
     return true;
@@ -148,93 +165,137 @@ bool Foam::sampledTriSurfaceMesh::update()
     // Does approximation by looking at the face centres only
     const pointField& fc = surface_.faceCentres();
 
-    cellLabels_.setSize(fc.size());
-
     meshSearch meshSearcher(mesh(), false);
 
+    const indexedOctree<treeDataPoint>& cellCentreTree =
+        meshSearcher.cellCentreTree();
+
+
+    // Global numbering for cells - only used to uniquely identify local cells.
+    globalIndex globalCells(mesh().nCells());
+    List<nearInfo> nearest(fc.size());
+    forAll(nearest, i)
+    {
+        nearest[i].first() = GREAT;
+        nearest[i].second() = labelMax;
+    }
+
+    // Search triangles using nearest on local mesh
     forAll(fc, triI)
     {
-        cellLabels_[triI] = meshSearcher.findCell
+        pointIndexHit nearInfo = cellCentreTree.findNearest
         (
             fc[triI],
-            -1,                 // seed
-            true                // use tree
+            sqr(GREAT)
         );
-    }
-
-
-    boolList include(surface_.size(), true);
-
-    label nNotFound = 0;
-
-    forAll(cellLabels_, triI)
-    {
-        if (cellLabels_[triI] == -1)
+        if (nearInfo.hit())
         {
-            include[triI] = false;
-            nNotFound++;
+            nearest[triI].first() = magSqr(nearInfo.hitPoint()-fc[triI]);
+            nearest[triI].second() = globalCells.toGlobal(nearInfo.index());
         }
     }
-    label nTotalNotFound = returnReduce(nNotFound, sumOp<label>());
+
+    // See which processor has the nearest.
+    Pstream::listCombineGather(nearest, nearestEqOp());
+    Pstream::listCombineScatter(nearest);
+
+    boolList include(surface_.size(), false);
+
+    cellLabels_.setSize(fc.size());
+    cellLabels_ = -1;
+
+    label nFound = 0;
+    forAll(nearest, triI)
+    {
+        if (nearest[triI].second() == labelMax)
+        {
+            // Not found on any processor. How to map?
+        }
+        else if (globalCells.isLocal(nearest[triI].second()))
+        {
+            cellLabels_[triI] = globalCells.toLocal(nearest[triI].second());
+
+            include[triI] = true;
+            nFound++;
+        }
+    }
+
 
     if (debug)
     {
-        Pout<< "surface:" << surface_.size()
-            << "  included:" << surface_.size()-nNotFound
-            << "  total not found" << nTotalNotFound << endl;
+        Pout<< "Local out of faces:" << cellLabels_.size()
+            << " keeping:" << nFound << endl;
     }
 
+    // Now subset the surface. Do not use triSurface::subsetMesh since requires
+    // original surface to be in compact numbering.
+
+    const triSurface& s = surface_;
+
+    // Compact to original triangle
+    labelList faceMap(s.size());
+    // Compact to original points
+    labelList pointMap(s.points().size());
+    // From original point to compact points
+    labelList reversePointMap(s.points().size(), -1);
 
-    // Add to master all triangles are outside all meshes.
     {
-        boolList onAnyProc(surface_.size(), false);
-        Pstream::listCombineGather(onAnyProc, orEqOp<bool>());
+        label newPointI = 0;
+        label newTriI = 0;
 
-        if (Pstream::master())
+        forAll(s, triI)
         {
-            label nAdded = 0;
-            forAll(onAnyProc, triI)
+            if (include[triI])
             {
-                if (!onAnyProc[triI])
+                faceMap[newTriI++] = triI;
+
+                const labelledTri& f = s[triI];
+                forAll(f, fp)
                 {
-                    include[triI] = true;
-                    nAdded++;
+                    if (reversePointMap[f[fp]] == -1)
+                    {
+                        pointMap[newPointI] = f[fp];
+                        reversePointMap[f[fp]] = newPointI++;
+                    }
                 }
             }
+        }
+        faceMap.setSize(newTriI);
+        pointMap.setSize(newPointI);
+    }
 
-            if (debug)
-            {
-                Pout<< "nAdded to master:" << nAdded << endl;
-            }
+    // Subset cellLabels
+    cellLabels_ = UIndirectList<label>(cellLabels_, faceMap)();
+
+    // Store any face per point
+    pointToFace_.setSize(pointMap.size());
+
+    // Create faces and points for subsetted surface
+    faceList& faces = this->storedFaces();
+    faces.setSize(faceMap.size());
+    forAll(faceMap, i)
+    {
+        const triFace& f = s[faceMap[i]];
+        triFace newF
+        (
+            reversePointMap[f[0]],
+            reversePointMap[f[1]],
+            reversePointMap[f[2]]
+        );
+        faces[i] = newF.triFaceFace();
+
+        forAll(newF, fp)
+        {
+            pointToFace_[newF[fp]] = i;
         }
     }
 
-
-    // Now subset the surface
-    triSurface localSurface
-    (
-        surface_.subsetMesh
-        (
-            include,
-            pointMap_,
-            faceMap_
-        )
-    );
-
-    // And convert into faces.
-    faceList& faces = this->storedFaces();
-    faces.setSize(localSurface.size());
-    forAll(localSurface, i)
-    {
-        faces[i] = localSurface[i].triFaceFace();
-    }
-
-    this->storedPoints() = localSurface.points();
+    this->storedPoints() = pointField(s.points(), pointMap);
 
     if (debug)
     {
         print(Pout);
-        Pout << endl;
+        Pout<< endl;
     }
 
     needsUpdate_ = false;
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
index eef8a16224..10ab6edb9d 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
@@ -29,15 +29,13 @@ Description
     A sampledSurface from a triSurfaceMesh. It samples on the points/triangles
     of the triSurface.
 
+    It samples using the cell nearest to the triangle centre so does
+    not check the cell the centre is actually in ...
+
     In parallel every processor just operates on the part of the surface
     where the face centres are inside the mesh. It is then up to the
     caller to stitch the partial surfaces together.
 
-    No check is done for triangle centres being on multiple processors
-    (should never happen but truncation errors ...)
-
-    Any value on a triangle outside the mesh will be set to 0.
-
 SourceFiles
     sampledTriSurfaceMesh.C
 
@@ -75,14 +73,11 @@ class sampledTriSurfaceMesh
         //- Track if the surface needs an update
         mutable bool needsUpdate_;
 
-        //- Local cell labels
+        //- From local surface triangle to mesh cell.
         labelList cellLabels_;
 
         //- From local surface back to surface_
-        labelList pointMap_;
-
-        //- From local surface back to surface_
-        labelList faceMap_;
+        labelList pointToFace_;
 
 
     // Private Member Functions
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
index 0dd08f8f26..b2a8ee956c 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
@@ -36,20 +36,12 @@ Foam::sampledTriSurfaceMesh::sampleField
 ) const
 {
     // One value per face
-    tmp<Field<Type> > tvalues(new Field<Type>(faceMap_.size()));
+    tmp<Field<Type> > tvalues(new Field<Type>(cellLabels_.size()));
     Field<Type>& values = tvalues();
 
-    forAll(faceMap_, i)
+    forAll(cellLabels_, triI)
     {
-        label cellI = cellLabels_[faceMap_[i]];
-        if (cellI != -1)
-        {
-            values[i] = vField[cellI];
-        }
-        else
-        {
-            values[i] = pTraits<Type>::zero;
-        }
+        values[triI] = vField[cellLabels_[triI]];
     }
 
     return tvalues;
@@ -64,24 +56,15 @@ Foam::sampledTriSurfaceMesh::interpolateField
 ) const
 {
     // One value per vertex
-    tmp<Field<Type> > tvalues(new Field<Type>(pointMap_.size()));
+    tmp<Field<Type> > tvalues(new Field<Type>(pointToFace_.size()));
     Field<Type>& values = tvalues();
 
-    forAll(pointMap_, i)
+    forAll(pointToFace_, pointI)
     {
-        label origPointI = pointMap_[i];
-        label origTriI = surface_.pointFaces()[origPointI][0];
+        label triI = pointToFace_[pointI];
+        label cellI = cellLabels_[triI];
 
-        label cellI = cellLabels_[origTriI];
-
-        if (cellI != -1)
-        {
-            values[i] = interpolator.interpolate(points()[i], cellI);
-        }
-        else
-        {
-            values[i] = pTraits<Type>::zero;
-        }
+        values[pointI] = interpolator.interpolate(points()[pointI], cellI);
     }
 
     return tvalues;