diff --git a/src/dynamicMesh/motionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C b/src/dynamicMesh/motionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C
index a79511ba08..da93c3d292 100644
--- a/src/dynamicMesh/motionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C
+++ b/src/dynamicMesh/motionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2015-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -59,6 +59,26 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+const Foam::coordSystem::cylindrical&
+Foam::displacementLayeredMotionMotionSolver::getCylindrical
+(
+    const label cellZoneI,
+    const dictionary& zoneDict
+)
+{
+    auto iter = cylSystems_.cfind(cellZoneI);
+
+    if (iter.found())
+    {
+        return *(iter.val());
+    }
+
+    cylSystems_.set(cellZoneI, new coordSystem::cylindrical(zoneDict));
+
+    return *cylSystems_[cellZoneI];
+}
+
+
 void Foam::displacementLayeredMotionMotionSolver::calcZoneMask
 (
     const label cellZoneI,
@@ -66,68 +86,49 @@ void Foam::displacementLayeredMotionMotionSolver::calcZoneMask
     bitSet& isZoneEdge
 ) const
 {
+    isZonePoint.resize(mesh().nPoints());
+    isZoneEdge.resize(mesh().nEdges());
+
     if (cellZoneI == -1)
     {
-        isZonePoint.setSize(mesh().nPoints());
         isZonePoint = true;
-
-        isZoneEdge.setSize(mesh().nEdges());
         isZoneEdge = true;
+        return;
     }
-    else
+
+
+    isZonePoint.reset();
+    isZoneEdge.reset();
+
+    const cellZone& cz = mesh().cellZones()[cellZoneI];
+
+    // Mark points, edges inside cellZone
+    for (const label celli : cz)
     {
-        const cellZone& cz = mesh().cellZones()[cellZoneI];
-
-        label nPoints = 0;
-        forAll(cz, i)
-        {
-            const labelList& cPoints = mesh().cellPoints(cz[i]);
-            forAll(cPoints, cPointi)
-            {
-                if (isZonePoint.set(cPoints[cPointi]))
-                {
-                    ++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.set(cEdges[cEdgeI]))
-                {
-                    ++nEdges;
-                }
-            }
-        }
-        syncTools::syncEdgeList
-        (
-            mesh(),
-            isZoneEdge,
-            orEqOp<unsigned int>(),
-            0
-        );
-
-        if (debug)
-        {
-            Info<< "On cellZone " << cz.name()
-                << " marked " << returnReduce(nPoints, sumOp<label>())
-                << " points and " << returnReduce(nEdges, sumOp<label>())
-                << " edges." << endl;
-        }
+        isZonePoint.set(mesh().cellPoints(celli));
+        isZoneEdge.set(mesh().cellEdges(celli));
     }
+
+    syncTools::syncPointList
+    (
+        mesh(),
+        isZonePoint,
+        orEqOp<unsigned int>(),
+        0
+    );
+
+    syncTools::syncEdgeList
+    (
+        mesh(),
+        isZoneEdge,
+        orEqOp<unsigned int>(),
+        0
+    );
+
+    DebugInfo
+        << "On cellZone " << cz.name() << " marked "
+        << returnReduce(isZonePoint.count(), sumOp<label>()) << " points and "
+        << returnReduce(isZoneEdge.count(), sumOp<label>()) << " edges" << nl;
 }
 
 
@@ -140,19 +141,23 @@ void Foam::displacementLayeredMotionMotionSolver::walkStructured
     const labelList& seedPoints,
     const vectorField& seedData,
     scalarField& distance,
-    vectorField& data
+    vectorField& data,
+    labelField& patchPoints
 ) const
 {
     List<pointEdgeStructuredWalk> seedInfo(seedPoints.size());
 
     forAll(seedPoints, i)
     {
+        const label pointi = seedPoints[i];
+
         seedInfo[i] = pointEdgeStructuredWalk
         (
-            points0()[seedPoints[i]],  // location of data
-            points0()[seedPoints[i]],  // previous location
+            points0()[pointi],  // location of data
+            points0()[pointi],  // previous location
             0.0,
-            seedData[i]
+            seedData[i],
+            pointi              // pass thru seed point id
         );
     }
 
@@ -162,36 +167,31 @@ void Foam::displacementLayeredMotionMotionSolver::walkStructured
     // Mark points inside cellZone.
     // Note that we use points0, not mesh.points()
     // so as not to accumulate errors.
-    forAll(isZonePoint, pointi)
+    for (const label pointi : isZonePoint)
     {
-        if (isZonePoint[pointi])
-        {
-            allPointInfo[pointi] = pointEdgeStructuredWalk
-            (
-                points0()[pointi],  // location of data
-                vector::max,        // not valid
-                0.0,
-                Zero        // passive data
-            );
-        }
+        allPointInfo[pointi] = pointEdgeStructuredWalk
+        (
+            points0()[pointi],  // location of data
+            point::max,         // not valid
+            0.0,
+            Zero                // passive data
+        );
     }
 
+
     // Current info on edges
     List<pointEdgeStructuredWalk> allEdgeInfo(mesh().nEdges());
 
     // Mark edges inside cellZone
-    forAll(isZoneEdge, edgeI)
+    for (const label edgei : isZoneEdge)
     {
-        if (isZoneEdge[edgeI])
-        {
-            allEdgeInfo[edgeI] = pointEdgeStructuredWalk
-            (
-                mesh().edges()[edgeI].centre(points0()),    // location of data
-                vector::max,                                // not valid
-                0.0,
-                Zero
-            );
-        }
+        allEdgeInfo[edgei] = pointEdgeStructuredWalk
+        (
+            mesh().edges()[edgei].centre(points0()),    // location of data
+            point::max,                                 // not valid
+            0.0,
+            Zero
+        );
     }
 
     // Walk
@@ -207,12 +207,17 @@ void Foam::displacementLayeredMotionMotionSolver::walkStructured
     );
 
     // Extract distance and passive data
-    forAll(allPointInfo, pointi)
+    for (const label pointi : isZonePoint)
     {
-        if (isZonePoint[pointi])
+        const auto& pointInfo = allPointInfo[pointi];
+
+        distance[pointi] = pointInfo.dist();
+        data[pointi] = pointInfo.data();
+
+        // Optional information
+        if (patchPoints.size())
         {
-            distance[pointi] = allPointInfo[pointi].dist();
-            data[pointi] = allPointInfo[pointi].data();
+            patchPoints[pointi] = pointInfo.index();
         }
     }
 }
@@ -229,8 +234,8 @@ Foam::displacementLayeredMotionMotionSolver::faceZoneEvaluate
     const label patchi
 ) const
 {
-    tmp<vectorField> tfld(new vectorField(meshPoints.size()));
-    vectorField& fld = tfld.ref();
+    auto tfld = tmp<vectorField>::New(meshPoints.size());
+    auto& fld = tfld.ref();
 
     const word type(dict.get<word>("type"));
 
@@ -265,7 +270,7 @@ Foam::displacementLayeredMotionMotionSolver::faceZoneEvaluate
         // Reads name of name of patch. Then get average point displacement on
         // patch. That becomes the value of fld.
         const word patchName(dict.get<word>("patch"));
-        label patchID = mesh().boundaryMesh().findPatchID(patchName);
+        const label patchID = mesh().boundaryMesh().findPatchID(patchName);
         pointField pdf
         (
             pointDisplacement_.boundaryField()[patchID].patchInternalField()
@@ -275,9 +280,11 @@ Foam::displacementLayeredMotionMotionSolver::faceZoneEvaluate
     else
     {
         FatalIOErrorInFunction(*this)
-            << "Unknown faceZonePatch type " << type << " for faceZone "
-            << fz.name() << exit(FatalIOError);
+            << "Unknown faceZonePatch type " << type
+            << " for faceZone " << fz.name() << nl
+            << exit(FatalIOError);
     }
+
     return tfld;
 }
 
@@ -288,8 +295,7 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
     const dictionary& zoneDict
 )
 {
-    bitSet isZonePoint(mesh().nPoints());
-    bitSet isZoneEdge(mesh().nEdges());
+    bitSet isZonePoint, isZoneEdge;
     calcZoneMask(cellZoneI, isZonePoint, isZoneEdge);
 
     const dictionary& patchesDict = zoneDict.subDict("boundaryField");
@@ -303,6 +309,7 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
             << exit(FatalIOError);
     }
 
+    PtrList<labelField> patchPoints(patchesDict.size());
     PtrList<scalarField> patchDist(patchesDict.size());
     PtrList<pointVectorField> patchDisp(patchesDict.size());
 
@@ -311,7 +318,7 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
     for (const entry& dEntry : patchesDict)
     {
         const word& faceZoneName = dEntry.keyword();
-        label zoneI = mesh().faceZones().findZoneID(faceZoneName);
+        const label zoneI = mesh().faceZones().findZoneID(faceZoneName);
         if (zoneI == -1)
         {
             FatalIOErrorInFunction(*this)
@@ -323,6 +330,7 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
         // Determine the points of the faceZone within the cellZone
         const faceZone& fz = mesh().faceZones()[zoneI];
 
+        patchPoints.set(patchi, new labelField(mesh().nPoints(), label(-1)));
         patchDist.set(patchi, new scalarField(mesh().nPoints()));
         patchDisp.set
         (
@@ -342,7 +350,7 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
             )
         );
 
-        patchi++;
+        ++patchi;
     }
 
 
@@ -382,17 +390,16 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
             patchi
         );
 
-        if (debug)
-        {
-            Info<< "For cellZone:" << cellZoneI
-                << " for faceZone:" << fz.name()
-                << " nPoints:" << tseed().size()
-                << " have patchField:"
-                << " max:" << gMax(tseed())
-                << " min:" << gMin(tseed())
-                << " avg:" << gAverage(tseed())
-                << endl;
-        }
+        DebugInfo
+            << "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
@@ -404,13 +411,14 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
             meshPoints,
             tseed,
             patchDist[patchi],
-            patchDisp[patchi]
+            patchDisp[patchi],
+            patchPoints[patchi]
         );
 
         // Implement real bc.
         patchDisp[patchi].correctBoundaryConditions();
 
-        patchi++;
+        ++patchi;
     }
 
 
@@ -435,14 +443,13 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
             dimensionedScalar(dimLength, Zero)
         );
 
-        forAll(distance, pointi)
+        for (const label pointi : isZonePoint)
         {
-            scalar d1 = patchDist[0][pointi];
-            scalar d2 = patchDist[1][pointi];
+            const scalar d1 = patchDist[0][pointi];
+            const scalar d2 = patchDist[1][pointi];
             if (d1 + d2 > SMALL)
             {
-                scalar s = d1/(d1 + d2);
-                distance[pointi] = s;
+                distance[pointi] = d1/(d1 + d2);
             }
         }
 
@@ -457,36 +464,43 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
 
     if (interpolationScheme == "oneSided")
     {
-        forAll(pointDisplacement_, pointi)
+        for (const label pointi : isZonePoint)
         {
-            if (isZonePoint[pointi])
-            {
-                pointDisplacement_[pointi] = patchDisp[0][pointi];
-            }
+            pointDisplacement_[pointi] = patchDisp[0][pointi];
         }
     }
     else if (interpolationScheme == "linear")
     {
-        forAll(pointDisplacement_, pointi)
+        for (const label pointi : isZonePoint)
         {
-            if (isZonePoint[pointi])
-            {
-                scalar d1 = patchDist[0][pointi];
-                scalar d2 = patchDist[1][pointi];
-                scalar s = d1/(d1 + d2 + VSMALL);
+            const scalar d1 = patchDist[0][pointi];
+            const scalar d2 = patchDist[1][pointi];
+            const scalar s = d1/(d1 + d2 + VSMALL);
 
-                const vector& pd1 = patchDisp[0][pointi];
-                const vector& pd2 = patchDisp[1][pointi];
+            const vector& pd1 = patchDisp[0][pointi];
+            const vector& pd2 = patchDisp[1][pointi];
 
-                pointDisplacement_[pointi] = (1 - s)*pd1 + s*pd2;
-            }
+            pointDisplacement_[pointi] = (1 - s)*pd1 + s*pd2;
         }
     }
+    else if (interpolationScheme == "cylindrical")
+    {
+        const coordSystem::cylindrical& cs =
+            this->getCylindrical(cellZoneI, zoneDict);
+
+        // May wish to implement alternative distance calculation here
+
+
+        FatalErrorInFunction
+            << "cylindrical interpolation not yet available" << nl
+            << "coordinate system " << cs << nl
+            << exit(FatalError);
+    }
     else
     {
         FatalErrorInFunction
             << "Invalid interpolationScheme: " << interpolationScheme
-            << ". Valid schemes are 'oneSided' and 'linear'"
+            << ". Valid schemes: (oneSided linear cylindrical)" << nl
             << exit(FatalError);
     }
 }
@@ -518,13 +532,6 @@ displacementLayeredMotionMotionSolver
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::displacementLayeredMotionMotionSolver::
-~displacementLayeredMotionMotionSolver()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::pointField>
@@ -553,7 +560,7 @@ void Foam::displacementLayeredMotionMotionSolver::solve()
         const word& cellZoneName = dEntry.keyword();
         const dictionary& regionDict = dEntry.dict();
 
-        label zoneI = mesh().cellZones().findZoneID(cellZoneName);
+        const label zoneI = mesh().cellZones().findZoneID(cellZoneName);
 
         Info<< "solving for zone: " << cellZoneName << endl;
 
@@ -586,11 +593,11 @@ void Foam::displacementLayeredMotionMotionSolver::updateMesh
 
     forAll(points0_, pointi)
     {
-        label oldPointi = mpm.pointMap()[pointi];
+        const label oldPointi = mpm.pointMap()[pointi];
 
         if (oldPointi >= 0)
         {
-            label masterPointi = mpm.reversePointMap()[oldPointi];
+            const label masterPointi = mpm.reversePointMap()[oldPointi];
 
             if ((masterPointi != pointi))
             {
diff --git a/src/dynamicMesh/motionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.H b/src/dynamicMesh/motionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.H
index 42213aa2e5..b4339eceaa 100644
--- a/src/dynamicMesh/motionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.H
+++ b/src/dynamicMesh/motionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.H
@@ -63,14 +63,17 @@ SourceFiles
 #define displacementLayeredMotionMotionSolver_H
 
 #include "displacementMotionSolver.H"
-#include "bitSet.H"
+#include "cylindricalCS.H"
+#include "PtrMap.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-// Forward class declarations
+// Forward Declarations
+class bitSet;
+
 
 /*---------------------------------------------------------------------------*\
              Class displacementLayeredMotionMotionSolver Declaration
@@ -80,8 +83,21 @@ class displacementLayeredMotionMotionSolver
 :
     public displacementMotionSolver
 {
+    // Private Data
+
+       //- Cylindrical coordinate systems (per cellZone)
+       PtrMap<coordSystem::cylindrical> cylSystems_;
+
+
     // Private Member Functions
 
+        //- Get existing or create new cylindrical system
+        const coordSystem::cylindrical& getCylindrical
+        (
+            const label cellZoneI,
+            const dictionary& zoneDict
+        );
+
         void calcZoneMask
         (
             const label cellZoneI,
@@ -97,7 +113,8 @@ class displacementLayeredMotionMotionSolver
             const labelList& seedPoints,
             const vectorField& seedData,
             scalarField& distance,
-            vectorField& data
+            vectorField& data,
+            labelField& patchPoints
         ) const;
 
         tmp<vectorField> faceZoneEvaluate
@@ -137,23 +154,23 @@ public:
         //- Construct from polyMesh and IOdictionary
         displacementLayeredMotionMotionSolver
         (
-            const polyMesh&,
-            const IOdictionary&
+            const polyMesh& mesh,
+            const IOdictionary& dict
         );
 
-        //- Construct from polyMesh, IOdictionary and displacement and
-        //  zero-time fields
+        //- Construct from polyMesh, IOdictionary, displacement and
+        //- zero-time fields
         displacementLayeredMotionMotionSolver
         (
-            const polyMesh&,
-            const IOdictionary&,
-            const pointVectorField&,
-            const pointIOField&
+            const polyMesh& mesh,
+            const IOdictionary& dict,
+            const pointVectorField& pointDisplacement,
+            const pointIOField& points0
         );
 
 
     //- Destructor
-    ~displacementLayeredMotionMotionSolver();
+    ~displacementLayeredMotionMotionSolver() = default;
 
 
     // Member Functions
@@ -165,7 +182,7 @@ public:
         virtual void solve();
 
         //- Update topology
-        virtual void updateMesh(const mapPolyMesh&);
+        virtual void updateMesh(const mapPolyMesh& mpm);
 };
 
 
diff --git a/src/dynamicMesh/motionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.C b/src/dynamicMesh/motionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.C
index 1066f0656e..c1843b8314 100644
--- a/src/dynamicMesh/motionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.C
+++ b/src/dynamicMesh/motionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -31,24 +31,25 @@ License
 
 Foam::Ostream& Foam::operator<<
 (
-    Foam::Ostream& os,
-    const Foam::pointEdgeStructuredWalk& wDist
+    Ostream& os,
+    const pointEdgeStructuredWalk& wDist
 )
 {
     return os
         << wDist.point0_ << wDist.previousPoint_
-        << wDist.dist_ << wDist.data_;
+        << wDist.dist_ << wDist.data_ << wDist.index_;
 }
 
+
 Foam::Istream& Foam::operator>>
 (
-    Foam::Istream& is,
-    Foam::pointEdgeStructuredWalk& wDist
+    Istream& is,
+    pointEdgeStructuredWalk& wDist
 )
 {
     return is
         >> wDist.point0_ >> wDist.previousPoint_
-        >> wDist.dist_ >> wDist.data_;
+        >> wDist.dist_ >> wDist.data_ >> wDist.index_;
 }
 
 
diff --git a/src/dynamicMesh/motionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.H b/src/dynamicMesh/motionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.H
index 204934b313..8fc9b6a06b 100644
--- a/src/dynamicMesh/motionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.H
+++ b/src/dynamicMesh/motionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -46,13 +46,9 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of classes
+// Forward declarations
 class polyPatch;
 class polyMesh;
-
-
-// Forward declaration of friend functions and operators
-
 class pointEdgeStructuredWalk;
 
 Istream& operator>>(Istream&, pointEdgeStructuredWalk&);
@@ -60,12 +56,12 @@ Ostream& operator<<(Ostream&, const pointEdgeStructuredWalk&);
 
 
 /*---------------------------------------------------------------------------*\
-                           Class pointEdgeStructuredWalk Declaration
+                   Class pointEdgeStructuredWalk Declaration
 \*---------------------------------------------------------------------------*/
 
 class pointEdgeStructuredWalk
 {
-    // Private data
+    // Private Data
 
         //- Starting location
         point point0_;
@@ -79,6 +75,10 @@ class pointEdgeStructuredWalk
         //- Passive data
         vector data_;
 
+        //- Index of passive data (optional)
+        label index_;
+
+
     // Private Member Functions
 
         //- Evaluate distance to point.
@@ -100,22 +100,29 @@ public:
         //- Construct from components
         inline pointEdgeStructuredWalk
         (
-            const point&,
-            const point&,
-            const scalar,
-            const vector&
+            const point& point0,
+            const point& previousPoint,
+            const scalar dist,
+            const vector& data,
+            const label index = -1
         );
 
 
     // Member Functions
 
-        // Access
+    // Access
 
-            inline bool inZone() const;
+        //- True if starting point is valid (ie, not point::max)
+        inline bool inZone() const;
 
-            inline scalar dist() const;
+        //- The distance information
+        inline scalar dist() const;
 
-            inline const vector& data() const;
+        //- Tracking data
+        inline const vector& data() const;
+
+        //- Index (if any) associated with data
+        inline label index() const;
 
 
         // Needed by meshWave
diff --git a/src/dynamicMesh/motionSolvers/displacement/layeredSolver/pointEdgeStructuredWalkI.H b/src/dynamicMesh/motionSolvers/displacement/layeredSolver/pointEdgeStructuredWalkI.H
index 811856b48c..e0c1cd7beb 100644
--- a/src/dynamicMesh/motionSolvers/displacement/layeredSolver/pointEdgeStructuredWalkI.H
+++ b/src/dynamicMesh/motionSolvers/displacement/layeredSolver/pointEdgeStructuredWalkI.H
@@ -30,7 +30,6 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-// Update this with w2.
 template<class TrackingData>
 inline bool Foam::pointEdgeStructuredWalk::update
 (
@@ -45,6 +44,7 @@ inline bool Foam::pointEdgeStructuredWalk::update
         dist_ = w2.dist_ + mag(point0_-w2.previousPoint_);
         previousPoint_ = point0_;
         data_ = w2.data_;
+        index_ = w2.index_;
 
         return true;
     }
@@ -55,29 +55,30 @@ inline bool Foam::pointEdgeStructuredWalk::update
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Null constructor
 inline Foam::pointEdgeStructuredWalk::pointEdgeStructuredWalk()
 :
-    point0_(vector::max),
-    previousPoint_(vector::max),
+    point0_(point::max),
+    previousPoint_(point::max),
     dist_(0),
-    data_(Zero)
+    data_(Zero),
+    index_(-1)
 {}
 
 
-// Construct from origin, distance
 inline Foam::pointEdgeStructuredWalk::pointEdgeStructuredWalk
 (
     const point& point0,
     const point& previousPoint,
     const scalar dist,
-    const vector& data
+    const vector& data,
+    const label index
 )
 :
     point0_(point0),
     previousPoint_(previousPoint),
     dist_(dist),
-    data_(data)
+    data_(data),
+    index_(index)
 {}
 
 
@@ -85,16 +86,10 @@ inline Foam::pointEdgeStructuredWalk::pointEdgeStructuredWalk
 
 inline bool Foam::pointEdgeStructuredWalk::inZone() const
 {
-    return point0_ != vector::max;
+    return point0_ != point::max;
 }
 
 
-//inline const Foam::point& Foam::pointEdgeStructuredWalk::previousPoint() const
-//{
-//    return previousPoint_;
-//}
-
-
 inline Foam::scalar Foam::pointEdgeStructuredWalk::dist() const
 {
     return dist_;
@@ -107,10 +102,16 @@ inline const Foam::vector& Foam::pointEdgeStructuredWalk::data() const
 }
 
 
+inline Foam::label Foam::pointEdgeStructuredWalk::index() const
+{
+    return index_;
+}
+
+
 template<class TrackingData>
 inline bool Foam::pointEdgeStructuredWalk::valid(TrackingData& td) const
 {
-    return previousPoint_ != vector::max;
+    return previousPoint_ != point::max;
 }