diff --git a/applications/utilities/surface/surfaceLambdaMuSmooth/surfaceLambdaMuSmooth.C b/applications/utilities/surface/surfaceLambdaMuSmooth/surfaceLambdaMuSmooth.C
index c76e466f83..cc1981e1d5 100644
--- a/applications/utilities/surface/surfaceLambdaMuSmooth/surfaceLambdaMuSmooth.C
+++ b/applications/utilities/surface/surfaceLambdaMuSmooth/surfaceLambdaMuSmooth.C
@@ -33,6 +33,9 @@ Description
     Provide an edgeMesh file containing points that are not to be moved during
     smoothing in order to preserve features.
 
+    lambda/mu smoothing: G. Taubin, IBM Research report Rc-19923 (02/01/95)
+    "A signal processing approach to fair surface design"
+
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
diff --git a/src/OpenFOAM/db/IOstreams/Fstreams/OFstream.C b/src/OpenFOAM/db/IOstreams/Fstreams/OFstream.C
index 3feefdaf40..1a09b7f473 100644
--- a/src/OpenFOAM/db/IOstreams/Fstreams/OFstream.C
+++ b/src/OpenFOAM/db/IOstreams/Fstreams/OFstream.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -104,10 +104,10 @@ Foam::OFstream::OFstream
     {
         if (debug)
         {
-            Info<< "IFstream::IFstream(const fileName&,"
-                   "streamFormat format=ASCII,"
-                   "versionNumber version=currentVersion) : "
-                   "could not open file for input\n"
+            Info<< "OFstream::OFstream(const fileName&,"
+                   "streamFormat, versionNumber, compressionType) : "
+                   "could not open file " << pathname
+                << "for input\n"
                    "in stream " << info() << Foam::endl;
         }
 
diff --git a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamBuffers.C b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamBuffers.C
index a1db6529f1..dfb60891d6 100644
--- a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamBuffers.C
+++ b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamBuffers.C
@@ -145,4 +145,19 @@ void Foam::PstreamBuffers::finishedSends(labelListList& sizes, const bool block)
 }
 
 
+void Foam::PstreamBuffers::clear()
+{
+    forAll(sendBuf_, i)
+    {
+        sendBuf_[i].clear();
+    }
+    forAll(recvBuf_, i)
+    {
+        recvBuf_[i].clear();
+    }
+    recvBufPos_ = 0;
+    finishedSendsCalled_ = false;
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamBuffers.H b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamBuffers.H
index 52b9960bde..8d5998fe82 100644
--- a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamBuffers.H
+++ b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamBuffers.H
@@ -112,8 +112,6 @@ class PstreamBuffers
 
         bool finishedSendsCalled_;
 
-    // Private Member Functions
-
 public:
 
     // Static data
@@ -155,6 +153,9 @@ public:
         //  non-blocking.
         void finishedSends(labelListList& sizes, const bool block = true);
 
+        //- Clear storage and reset
+        void clear();
+
 };
 
 
diff --git a/src/edgeMesh/Make/files b/src/edgeMesh/Make/files
index ed4cec5d95..0438dd204e 100644
--- a/src/edgeMesh/Make/files
+++ b/src/edgeMesh/Make/files
@@ -35,8 +35,6 @@ eem = $(em)/extendedEdgeMesh
 $(eem)/extendedEdgeMesh.C
 $(eem)/extendedEdgeMeshNew.C
 
-$(eem)/extendedEdgeMeshFormats/obj/OBJextendedEdgeFormat.C
-$(eem)/extendedEdgeMeshFormats/obj/OBJextendedEdgeFormatRunTime.C
 $(eem)/extendedEdgeMeshFormats/extendedEdgeMeshFormat/extendedEdgeMeshFormat.C
 $(eem)/extendedEdgeMeshFormats/extendedEdgeMeshFormat/extendedEdgeMeshFormatRunTime.C
 
diff --git a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C
index 999d2f0651..5ac41224e3 100644
--- a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C
+++ b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C
@@ -30,6 +30,7 @@ License
 #include "Time.H"
 #include "OBJstream.H"
 #include "DynamicField.H"
+#include "edgeMeshFormatsCore.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -104,6 +105,68 @@ Foam::label Foam::extendedEdgeMesh::nPointTypes = 4;
 Foam::label Foam::extendedEdgeMesh::nEdgeTypes = 5;
 
 
+Foam::wordHashSet Foam::extendedEdgeMesh::readTypes()
+{
+    return wordHashSet(*fileExtensionConstructorTablePtr_);
+}
+
+
+Foam::wordHashSet Foam::extendedEdgeMesh::writeTypes()
+{
+    return wordHashSet(*writefileExtensionMemberFunctionTablePtr_);
+}
+
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+bool Foam::extendedEdgeMesh::canReadType
+(
+    const word& ext,
+    const bool verbose
+)
+{
+    return edgeMeshFormatsCore::checkSupport
+    (
+        readTypes(),
+        ext,
+        verbose,
+        "reading"
+   );
+}
+
+
+bool Foam::extendedEdgeMesh::canWriteType
+(
+    const word& ext,
+    const bool verbose
+)
+{
+    return edgeMeshFormatsCore::checkSupport
+    (
+        writeTypes(),
+        ext,
+        verbose,
+        "writing"
+    );
+}
+
+
+bool Foam::extendedEdgeMesh::canRead
+(
+    const fileName& name,
+    const bool verbose
+)
+{
+    word ext = name.ext();
+    if (ext == "gz")
+    {
+        ext = name.lessExt().ext();
+    }
+    return canReadType(ext, verbose);
+}
+
+
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
 Foam::extendedEdgeMesh::pointStatus
diff --git a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.H b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.H
index 2648321888..58e2ba8cab 100644
--- a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.H
+++ b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.H
@@ -225,6 +225,18 @@ public:
         //- Number of possible feature edge types (i.e. number of slices)
         static label nEdgeTypes;
 
+        //- Can we read this file format?
+        static bool canRead(const fileName&, const bool verbose=false);
+
+        //- Can we read this file format?
+        static bool canReadType(const word& ext, const bool verbose=false);
+
+        //- Can we write this file format type?
+        static bool canWriteType(const word& ext, const bool verbose=false);
+
+        static wordHashSet readTypes();
+        static wordHashSet writeTypes();
+
 
     // Constructors
 
diff --git a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/obj/OBJextendedEdgeFormat.C b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/obj/OBJextendedEdgeFormat.C
deleted file mode 100644
index d63ad27ed9..0000000000
--- a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/obj/OBJextendedEdgeFormat.C
+++ /dev/null
@@ -1,66 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "OBJextendedEdgeFormat.H"
-#include "OBJedgeFormat.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::fileFormats::OBJextendedEdgeFormat::OBJextendedEdgeFormat
-(
-    const fileName& filename
-)
-{
-    read(filename);
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::fileFormats::OBJextendedEdgeFormat::read(const fileName& filename)
-{
-    edgeMesh em(filename);
-
-    clear();
-
-    // Note: should transfer here instead
-    storedPoints() = em.points();
-    storedEdges() = em.edges();
-
-    return true;
-}
-
-
-void Foam::fileFormats::OBJextendedEdgeFormat::write
-(
-    const fileName& filename,
-    const extendedEdgeMesh& mesh
-)
-{
-    OBJedgeFormat::write(filename, mesh);
-}
-
-
-// ************************************************************************* //
diff --git a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/obj/OBJextendedEdgeFormat.H b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/obj/OBJextendedEdgeFormat.H
deleted file mode 100644
index a65409f446..0000000000
--- a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/obj/OBJextendedEdgeFormat.H
+++ /dev/null
@@ -1,119 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::fileFormats::OBJextendedEdgeFormat
-
-Description
-    Provide a means of reading/writing Alias/Wavefront OBJ format.
-
-    Does not handle negative vertex indices.
-
-SourceFiles
-    OBJextendedEdgeFormat.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef OBJextendedEdgeFormat_H
-#define OBJextendedEdgeFormat_H
-
-#include "extendedEdgeMesh.H"
-//#include "IFstream.H"
-//#include "Ostream.H"
-#include "OFstream.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace fileFormats
-{
-
-/*---------------------------------------------------------------------------*\
-                        Class OBJextendedEdgeFormat Declaration
-\*---------------------------------------------------------------------------*/
-
-class OBJextendedEdgeFormat
-:
-    public extendedEdgeMesh
-{
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        OBJextendedEdgeFormat(const OBJextendedEdgeFormat&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const OBJextendedEdgeFormat&);
-
-
-public:
-
-    // Constructors
-
-        //- Construct from file name
-        OBJextendedEdgeFormat(const fileName&);
-
-
-//    // Selectors
-//
-//        //- Read file and return surface
-//        static autoPtr<extendedEdgeMesh> New(const fileName& name)
-//        {
-//            return autoPtr<extendedEdgeMesh>
-//            (
-//                new OBJextendedEdgeFormat(name)
-//            );
-//        }
-//
-
-    //- Destructor
-    virtual ~OBJextendedEdgeFormat()
-    {}
-
-
-    // Member Functions
-
-        //- Write surface mesh components by proxy
-        static void write(const fileName&, const extendedEdgeMesh&);
-
-        //- Read from file
-        virtual bool read(const fileName&);
-
-        //- Write object file
-        virtual void write(const fileName& name) const
-        {
-            write(name, *this);
-        }
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace fileFormats
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/obj/OBJextendedEdgeFormatRunTime.C b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/obj/OBJextendedEdgeFormatRunTime.C
deleted file mode 100644
index c88513a104..0000000000
--- a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/obj/OBJextendedEdgeFormatRunTime.C
+++ /dev/null
@@ -1,60 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "OBJextendedEdgeFormat.H"
-
-#include "addToRunTimeSelectionTable.H"
-#include "addToMemberFunctionSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace fileFormats
-{
-
-// read extendedEdgeMesh
-addNamedToRunTimeSelectionTable
-(
-    extendedEdgeMesh,
-    OBJextendedEdgeFormat,
-    fileExtension,
-    obj
-);
-
-//// write extendedEdgeMesh
-//addNamedToMemberFunctionSelectionTable
-//(
-//    extendedEdgeMesh,
-//    OBJextendedEdgeFormat,
-//    write,
-//    fileExtension,
-//    obj
-//);
-
-}
-}
-
-// ************************************************************************* //
diff --git a/src/lagrangian/basic/Cloud/Cloud.C b/src/lagrangian/basic/Cloud/Cloud.C
index c0e43c1686..3773de6606 100644
--- a/src/lagrangian/basic/Cloud/Cloud.C
+++ b/src/lagrangian/basic/Cloud/Cloud.C
@@ -219,22 +219,33 @@ void Foam::Cloud<ParticleType>::move(TrackData& td, const scalar trackTime)
     // Reset nTrackingRescues
     nTrackingRescues_ = 0;
 
+
+    // List of lists of particles to be transfered for all of the
+    // neighbour processors
+    List<IDLList<ParticleType> > particleTransferLists
+    (
+        neighbourProcs.size()
+    );
+
+    // List of destination processorPatches indices for all of the
+    // neighbour processors
+    List<DynamicList<label> > patchIndexTransferLists
+    (
+        neighbourProcs.size()
+    );
+
+    // Allocate transfer buffers
+    PstreamBuffers pBufs(Pstream::nonBlocking);
+
+
     // While there are particles to transfer
     while (true)
     {
-        // List of lists of particles to be transfered for all of the
-        // neighbour processors
-        List<IDLList<ParticleType> > particleTransferLists
-        (
-            neighbourProcs.size()
-        );
-
-        // List of destination processorPatches indices for all of the
-        // neighbour processors
-        List<DynamicList<label> > patchIndexTransferLists
-        (
-            neighbourProcs.size()
-        );
+        particleTransferLists = IDLList<ParticleType>();
+        forAll(patchIndexTransferLists, i)
+        {
+            patchIndexTransferLists[i].clear();
+        }
 
         // Loop over all particles
         forAllIter(typename Cloud<ParticleType>, *this, pIter)
@@ -288,8 +299,9 @@ void Foam::Cloud<ParticleType>::move(TrackData& td, const scalar trackTime)
             break;
         }
 
-        // Allocate transfer buffers
-        PstreamBuffers pBufs(Pstream::nonBlocking);
+
+        // Clear transfer buffers
+        pBufs.clear();
 
         // Stream into send buffers
         forAll(particleTransferLists, i)
@@ -308,12 +320,12 @@ void Foam::Cloud<ParticleType>::move(TrackData& td, const scalar trackTime)
             }
         }
 
-        // Set up transfers when in non-blocking mode. Returns sizes (in bytes)
-        // to be sent/received.
-        labelListList allNTrans(Pstream::nProcs());
 
+        // Start sending. Sets number of bytes transferred
+        labelListList allNTrans(Pstream::nProcs());
         pBufs.finishedSends(allNTrans);
 
+
         bool transfered = false;
 
         forAll(allNTrans, i)
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
index 496f980a11..4e0f1d160b 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
@@ -310,14 +310,13 @@ class autoSnapDriver
                 //- Find point on nearest feature edge (within searchDist).
                 //  Return point and feature
                 //  and store feature-edge to mesh-point and vice versa
-                pointIndexHit findNearFeatureEdge
+                Tuple2<label, pointIndexHit> findNearFeatureEdge
                 (
                     const indirectPrimitivePatch& pp,
                     const scalarField& snapDist,
                     const label pointI,
                     const point& estimatedPt,
 
-                    label& featI,
                     List<List<DynamicList<point> > >&,
                     List<List<DynamicList<pointConstraint> > >&,
                     vectorField&,
@@ -330,7 +329,7 @@ class autoSnapDriver
                 //  If another mesh point already referring to this feature
                 //  point and further away, reset that one to a near feature
                 //  edge (using findNearFeatureEdge above)
-                labelPair findNearFeaturePoint
+                Tuple2<label, pointIndexHit> findNearFeaturePoint
                 (
                     const indirectPrimitivePatch& pp,
                     const scalarField& snapDist,
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
index 67fee3ef96..03615a76bd 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
@@ -170,7 +170,7 @@ void Foam::autoSnapDriver::smoothAndConstrain
 
                     if (isMasterEdge[meshEdges[edgeI]])
                     {
-                        label nbrPointI = edges[pEdges[i]].otherVertex(pointI);
+                        label nbrPointI = edges[edgeI].otherVertex(pointI);
                         if (constraints[nbrPointI].first() >= nConstraints)
                         {
                             dispSum[pointI] += disp[nbrPointI];
@@ -213,80 +213,6 @@ void Foam::autoSnapDriver::smoothAndConstrain
         }
     }
 }
-//XXXXXX
-//TODO: make proper parallel so coupled edges don't have double influence
-//void Foam::autoSnapDriver::smoothAndConstrain2
-//(
-//    const bool applyConstraints,
-//    const indirectPrimitivePatch& pp,
-//    const List<pointConstraint>& constraints,
-//    vectorField& disp
-//) const
-//{
-//    const fvMesh& mesh = meshRefiner_.mesh();
-//
-//    for (label avgIter = 0; avgIter < 20; avgIter++)
-//    {
-//        vectorField dispSum(pp.nPoints(), vector::zero);
-//        labelList dispCount(pp.nPoints(), 0);
-//
-//        const labelListList& pointEdges = pp.pointEdges();
-//        const edgeList& edges = pp.edges();
-//
-//        forAll(pointEdges, pointI)
-//        {
-//            const labelList& pEdges = pointEdges[pointI];
-//
-//            forAll(pEdges, i)
-//            {
-//                label nbrPointI = edges[pEdges[i]].otherVertex(pointI);
-//                dispSum[pointI] += disp[nbrPointI];
-//                dispCount[pointI]++;
-//            }
-//        }
-//
-//        syncTools::syncPointList
-//        (
-//            mesh,
-//            pp.meshPoints(),
-//            dispSum,
-//            plusEqOp<point>(),
-//            vector::zero,
-//            mapDistribute::transform()
-//        );
-//        syncTools::syncPointList
-//        (
-//            mesh,
-//            pp.meshPoints(),
-//            dispCount,
-//            plusEqOp<label>(),
-//            0,
-//            mapDistribute::transform()
-//        );
-//
-//        // Constraints
-//        forAll(constraints, pointI)
-//        {
-//            if (dispCount[pointI] > 0)// && constraints[pointI].first() <= 1)
-//            {
-//                // Mix my displacement with neighbours' displacement
-//                disp[pointI] =
-//                    0.5
-//                   *(disp[pointI] + dispSum[pointI]/dispCount[pointI]);
-//
-//                if (applyConstraints)
-//                {
-//                    disp[pointI] = transform
-//                    (
-//                        constraints[pointI].constraintTransformation(),
-//                        disp[pointI]
-//                    );
-//                }
-//            }
-//        }
-//    }
-//}
-//XXXXXX
 
 
 void Foam::autoSnapDriver::calcNearestFace
@@ -716,7 +642,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
         pNormals = List<point>(pNormals, visitOrder);
         pDisp = List<point>(pDisp, visitOrder);
         pFc = List<point>(pFc, visitOrder);
-        pFid = UIndirectList<label>(pFid, visitOrder);
+        pFid = UIndirectList<label>(pFid, visitOrder)();
     }
 }
 
@@ -1412,14 +1338,14 @@ Foam::labelPair Foam::autoSnapDriver::findDiagonalAttraction
 }
 
 
-Foam::pointIndexHit Foam::autoSnapDriver::findNearFeatureEdge
+Foam::Tuple2<Foam::label, Foam::pointIndexHit>
+Foam::autoSnapDriver::findNearFeatureEdge
 (
     const indirectPrimitivePatch& pp,
     const scalarField& snapDist,
     const label pointI,
     const point& estimatedPt,
 
-    label& featI,
     List<List<DynamicList<point> > >& edgeAttractors,
     List<List<DynamicList<pointConstraint> > >& edgeConstraints,
     vectorField& patchAttraction,
@@ -1430,16 +1356,18 @@ Foam::pointIndexHit Foam::autoSnapDriver::findNearFeatureEdge
 
     labelList nearEdgeFeat;
     List<pointIndexHit> nearEdgeInfo;
+    vectorField nearNormal;
     features.findNearestEdge
     (
         pointField(1, estimatedPt),
         scalarField(1, sqr(snapDist[pointI])),
         nearEdgeFeat,
-        nearEdgeInfo
+        nearEdgeInfo,
+        nearNormal
     );
 
     const pointIndexHit& nearInfo = nearEdgeInfo[0];
-    featI = nearEdgeFeat[0];
+    label featI = nearEdgeFeat[0];
 
     if (nearInfo.hit())
     {
@@ -1449,12 +1377,7 @@ Foam::pointIndexHit Foam::autoSnapDriver::findNearFeatureEdge
         (
             nearInfo.hitPoint()
         );
-        pointConstraint c;
-        const edge e = features[featI].edges()[nearInfo.index()];
-        vector eVec = e.vec(features[featI].points());
-        eVec /= mag(eVec)+VSMALL;
-        c.first() = 2;
-        c.second() = eVec;
+        pointConstraint c(Tuple2<label, vector>(2, nearNormal[0]));
         edgeConstraints[featI][nearInfo.index()].append(c);
 
         // Store for later use
@@ -1462,11 +1385,12 @@ Foam::pointIndexHit Foam::autoSnapDriver::findNearFeatureEdge
             nearInfo.hitPoint()-pp.localPoints()[pointI];
         patchConstraints[pointI] = c;
     }
-    return nearInfo;
+    return Tuple2<label, pointIndexHit>(featI, nearInfo);
 }
 
 
-Foam::labelPair Foam::autoSnapDriver::findNearFeaturePoint
+Foam::Tuple2<Foam::label, Foam::pointIndexHit>
+Foam::autoSnapDriver::findNearFeaturePoint
 (
     const indirectPrimitivePatch& pp,
     const scalarField& snapDist,
@@ -1487,26 +1411,23 @@ Foam::labelPair Foam::autoSnapDriver::findNearFeaturePoint
     const refinementFeatures& features = meshRefiner_.features();
 
     labelList nearFeat;
-    labelList nearIndex;
+    List<pointIndexHit> nearInfo;
     features.findNearestPoint
     (
         pointField(1, estimatedPt),
         scalarField(1, sqr(snapDist[pointI])),
         nearFeat,
-        nearIndex
+        nearInfo
     );
 
     label featI = nearFeat[0];
-    label featPointI = -1;
 
     if (featI != -1)
     {
         const point& pt = pp.localPoints()[pointI];
 
-        const treeDataPoint& shapes =
-            features.pointTrees()[featI].shapes();
-        featPointI = shapes.pointLabels()[nearIndex[0]];
-        const point& featPt = shapes.points()[featPointI];
+        label featPointI = nearInfo[0].index();
+        const point& featPt = nearInfo[0].hitPoint();
         scalar distSqr = magSqr(featPt-pt);
 
         // Check if already attracted
@@ -1537,7 +1458,6 @@ Foam::labelPair Foam::autoSnapDriver::findNearFeaturePoint
                 patchAttraction[oldPointI] = vector::zero;
                 patchConstraints[oldPointI] = pointConstraint();
 
-                label edgeFeatI;
                 findNearFeatureEdge
                 (
                     pp,
@@ -1545,7 +1465,6 @@ Foam::labelPair Foam::autoSnapDriver::findNearFeaturePoint
                     oldPointI,
                     pp.localPoints()[oldPointI],
 
-                    edgeFeatI,
                     edgeAttractors,
                     edgeConstraints,
                     patchAttraction,
@@ -1566,7 +1485,7 @@ Foam::labelPair Foam::autoSnapDriver::findNearFeaturePoint
         }
     }
 
-    return labelPair(featI, featPointI);
+    return Tuple2<label, pointIndexHit>(featI, nearInfo[0]);
 }
 
 
@@ -1637,7 +1556,6 @@ void Foam::autoSnapDriver::determineFeatures
             << featurePointStr().name() << endl;
     }
 
-    const refinementFeatures& features = meshRefiner_.features();
 
     forAll(pp.localPoints(), pointI)
     {
@@ -1664,6 +1582,7 @@ void Foam::autoSnapDriver::determineFeatures
             constraint
         );
 
+
         if
         (
             (constraint.first() > patchConstraints[pointI].first())
@@ -1696,15 +1615,14 @@ void Foam::autoSnapDriver::determineFeatures
                         // Behave like when having two surface normals so
                         // attract to nearest feature edge (with a guess for
                         // the multipatch point as starting point)
-                        label featI = -1;
-                        pointIndexHit nearInfo = findNearFeatureEdge
+                        Tuple2<label, pointIndexHit> nearInfo =
+                        findNearFeatureEdge
                         (
                             pp,
                             snapDist,
                             pointI,
                             multiPatchPt.hitPoint(),        //estimatedPt
 
-                            featI,
                             edgeAttractors,
                             edgeConstraints,
 
@@ -1712,14 +1630,15 @@ void Foam::autoSnapDriver::determineFeatures
                             patchConstraints
                         );
 
-                        if (nearInfo.hit())
+                        const pointIndexHit& info = nearInfo.second();
+                        if (info.hit())
                         {
                             // Dump
                             if (featureEdgeStr.valid())
                             {
                                 featureEdgeStr().write
                                 (
-                                    linePointRef(pt, nearInfo.hitPoint())
+                                    linePointRef(pt, info.hitPoint())
                                 );
                             }
                         }
@@ -1729,7 +1648,7 @@ void Foam::autoSnapDriver::determineFeatures
                             {
                                 missedEdgeStr().write
                                 (
-                                    linePointRef(pt, nearInfo.missPoint())
+                                    linePointRef(pt, info.missPoint())
                                 );
                             }
                         }
@@ -1746,15 +1665,13 @@ void Foam::autoSnapDriver::determineFeatures
                 // Determine nearest point on feature edge. Store constraint
                 // (calculated from feature edge, alternative would be to
                 //  use constraint calculated from both surfaceNormals)
-                label featI = -1;
-                pointIndexHit nearInfo = findNearFeatureEdge
+                Tuple2<label, pointIndexHit> nearInfo = findNearFeatureEdge
                 (
                     pp,
                     snapDist,
                     pointI,
                     estimatedPt,
 
-                    featI,
                     edgeAttractors,
                     edgeConstraints,
 
@@ -1762,14 +1679,15 @@ void Foam::autoSnapDriver::determineFeatures
                     patchConstraints
                 );
 
-                if (nearInfo.hit())
+                // Dump to obj
+                const pointIndexHit& info = nearInfo.second();
+                if (info.hit())
                 {
-                    // Dump
                     if (featureEdgeStr.valid())
                     {
                         featureEdgeStr().write
                         (
-                            linePointRef(pt, nearInfo.hitPoint())
+                            linePointRef(pt, info.hitPoint())
                         );
                     }
                 }
@@ -1779,7 +1697,7 @@ void Foam::autoSnapDriver::determineFeatures
                     {
                         missedEdgeStr().write
                         (
-                            linePointRef(pt, nearInfo.missPoint())
+                            linePointRef(pt, info.missPoint())
                         );
                     }
                 }
@@ -1789,7 +1707,7 @@ void Foam::autoSnapDriver::determineFeatures
                 // Mark point on the nearest feature point.
                 const point estimatedPt(pt + patchAttraction[pointI]);
 
-                labelPair nearInfo = findNearFeaturePoint
+                Tuple2<label, pointIndexHit> nearInfo = findNearFeaturePoint
                 (
                     pp,
                     snapDist,
@@ -1807,18 +1725,10 @@ void Foam::autoSnapDriver::determineFeatures
                     patchConstraints
                 );
 
-                if (nearInfo.first() != -1)
+                const pointIndexHit& info = nearInfo.second();
+                if (info.hit() && featurePointStr.valid())
                 {
-                    // Dump
-                    if (featurePointStr.valid())
-                    {
-                        const treeDataPoint& shapes =
-                            features.pointTrees()[nearInfo.first()].shapes();
-                        const point& featPt =
-                            shapes.points()[nearInfo.second()];
-
-                        featurePointStr().write(linePointRef(pt, featPt));
-                    }
+                    featurePointStr().write(linePointRef(pt, info.hitPoint()));
                 }
             }
         }
@@ -2029,22 +1939,20 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
 
             if (pointStatus[pointI] == 0)   // baffle edge
             {
-                label featI;
-                const pointIndexHit nearInfo = findNearFeatureEdge
+                Tuple2<label, pointIndexHit> nearInfo = findNearFeatureEdge
                 (
                     pp,
                     snapDist,
                     pointI,
                     pt,
 
-                    featI,
                     edgeAttractors,
                     edgeConstraints,
                     rawPatchAttraction,
                     rawPatchConstraints
                 );
 
-                if (!nearInfo.hit())
+                if (!nearInfo.second().hit())
                 {
                     //Pout<< "*** Failed to find close edge to point " << pt
                     //    << endl;
@@ -2053,23 +1961,21 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
             else if (pointStatus[pointI] == 1)   // baffle point
             {
                 labelList nearFeat;
-                labelList nearIndex;
+                List<pointIndexHit> nearInfo;
                 features.findNearestPoint
                 (
                     pointField(1, pt),
                     scalarField(1, sqr(snapDist[pointI])),
                     nearFeat,
-                    nearIndex
+                    nearInfo
                 );
 
                 label featI = nearFeat[0];
 
                 if (featI != -1)
                 {
-                    const treeDataPoint& shapes =
-                        features.pointTrees()[featI].shapes();
-                    label featPointI = shapes.pointLabels()[nearIndex[0]];
-                    const point& featPt = shapes.points()[featPointI];
+                    label featPointI = nearInfo[0].index();
+                    const point& featPt = nearInfo[0].hitPoint();
                     scalar distSqr = magSqr(featPt-pt);
 
                     // Check if already attracted
@@ -2099,7 +2005,6 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                             // The current point is closer so wins. Reset
                             // the old point to attract to nearest edge
                             // instead.
-                            label edgeFeatI;
                             findNearFeatureEdge
                             (
                                 pp,
@@ -2107,7 +2012,6 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                                 oldPointI,
                                 pp.localPoints()[oldPointI],
 
-                                edgeFeatI,
                                 edgeAttractors,
                                 edgeConstraints,
                                 rawPatchAttraction,
@@ -2130,7 +2034,6 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                     //    << " for baffle-feature-point " << pt
                     //    << endl;
 
-                    label featI;
                     findNearFeatureEdge
                     (
                         pp,
@@ -2138,7 +2041,6 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                         pointI,
                         pt,                     // starting point
 
-                        featI,
                         edgeAttractors,
                         edgeConstraints,
                         rawPatchAttraction,
diff --git a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C
index a6db1ff0b9..0aa74be6c9 100644
--- a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C
+++ b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C
@@ -51,7 +51,6 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-
 Foam::labelList Foam::medialAxisMeshMover::getFixedValueBCs
 (
     const pointVectorField& fld
@@ -118,63 +117,6 @@ Foam::medialAxisMeshMover::getPatch
 }
 
 
-void Foam::medialAxisMeshMover::calculateEdgeWeights
-(
-    const PackedBoolList& isMasterEdge,
-    const labelList& meshEdges,
-    const labelList& meshPoints,
-    const edgeList& edges,
-    scalarField& edgeWeights,
-    scalarField& invSumWeight
-) const
-{
-    const pointField& pts = mesh().points();
-
-    // Calculate edgeWeights and inverse sum of edge weights
-    edgeWeights.setSize(meshEdges.size());
-    invSumWeight.setSize(meshPoints.size());
-
-    forAll(edges, edgeI)
-    {
-        const edge& e = edges[edgeI];
-        scalar eMag = max
-        (
-            VSMALL,
-            mag
-            (
-                pts[meshPoints[e[1]]]
-              - pts[meshPoints[e[0]]]
-            )
-        );
-        edgeWeights[edgeI] = 1.0/eMag;
-    }
-
-    // Sum per point all edge weights
-    weightedSum
-    (
-        mesh(),
-        isMasterEdge,
-        meshEdges,
-        meshPoints,
-        edges,
-        edgeWeights,
-        scalarField(meshPoints.size(), 1.0),  // data
-        invSumWeight
-    );
-
-    // Inplace invert
-    forAll(invSumWeight, pointI)
-    {
-        scalar w = invSumWeight[pointI];
-
-        if (w > 0.0)
-        {
-            invSumWeight[pointI] = 1.0/w;
-        }
-    }
-}
-
-
 void Foam::medialAxisMeshMover::smoothPatchNormals
 (
     const label nSmoothDisp,
@@ -193,8 +135,9 @@ void Foam::medialAxisMeshMover::smoothPatchNormals
 
     scalarField edgeWeights(meshEdges.size());
     scalarField invSumWeight(meshPoints.size());
-    calculateEdgeWeights
+    meshRefinement::calculateEdgeWeights
     (
+        mesh(),
         isMasterEdge,
         meshEdges,
         meshPoints,
@@ -207,7 +150,7 @@ void Foam::medialAxisMeshMover::smoothPatchNormals
     vectorField average;
     for (label iter = 0; iter < nSmoothDisp; iter++)
     {
-        weightedSum
+        meshRefinement::weightedSum
         (
             mesh(),
             isMasterEdge,
@@ -284,8 +227,9 @@ void Foam::medialAxisMeshMover::smoothNormals
 
     scalarField edgeWeights(meshEdges.size());
     scalarField invSumWeight(meshPoints.size());
-    calculateEdgeWeights
+    meshRefinement::calculateEdgeWeights
     (
+        mesh(),
         isMasterEdge,
         meshEdges,
         meshPoints,
@@ -297,7 +241,7 @@ void Foam::medialAxisMeshMover::smoothNormals
     vectorField average;
     for (label iter = 0; iter < nSmoothDisp; iter++)
     {
-        weightedSum
+        meshRefinement::weightedSum
         (
             mesh(),
             isMasterEdge,
@@ -1059,8 +1003,9 @@ void Foam::medialAxisMeshMover::minSmoothField
 
     scalarField edgeWeights(meshEdges.size());
     scalarField invSumWeight(meshPoints.size());
-    calculateEdgeWeights
+    meshRefinement::calculateEdgeWeights
     (
+        mesh(),
         isMasterEdge,
         meshEdges,
         meshPoints,
@@ -1075,7 +1020,7 @@ void Foam::medialAxisMeshMover::minSmoothField
     for (label iter = 0; iter < nSmoothDisp; iter++)
     {
         scalarField average(pp.nPoints());
-        weightedSum
+        meshRefinement::weightedSum
         (
             mesh(),
             isMasterEdge,
@@ -1534,8 +1479,9 @@ void Foam::medialAxisMeshMover::smoothLambdaMuDisplacement
     // Calculate inverse sum of weights
     scalarField edgeWeights(meshEdges.size());
     scalarField invSumWeight(meshPoints.size());
-    calculateEdgeWeights
+    meshRefinement::calculateEdgeWeights
     (
+        mesh(),
         isMasterEdge,
         meshEdges,
         meshPoints,
@@ -1554,7 +1500,7 @@ void Foam::medialAxisMeshMover::smoothLambdaMuDisplacement
 
     for (label iter = 0; iter < nSmoothDisp; iter++)
     {
-        weightedSum
+        meshRefinement::weightedSum
         (
             mesh(),
             isMasterEdge,
@@ -1575,7 +1521,7 @@ void Foam::medialAxisMeshMover::smoothLambdaMuDisplacement
             }
         }
 
-        weightedSum
+        meshRefinement::weightedSum
         (
             mesh(),
             isMasterEdge,
diff --git a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H
index b402ecba25..1c270de274 100644
--- a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H
+++ b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H
@@ -105,31 +105,6 @@ class medialAxisMeshMover
             const labelList&
         );
 
-        //- Weighted sum (over all subset of mesh points) by
-        //  summing contribution from (master) edges
-        template<class Type>
-        static void weightedSum
-        (
-            const polyMesh& mesh,
-            const PackedBoolList& isMasterEdge,
-            const labelList& meshEdges,
-            const labelList& meshPoints,
-            const edgeList& edges,
-            const scalarField& edgeWeights,
-            const Field<Type>& data,
-            Field<Type>& sum
-        );
-
-        void calculateEdgeWeights
-        (
-            const PackedBoolList& isMasterEdge,
-            const labelList& meshEdges,
-            const labelList& meshPoints,
-            const edgeList& edges,
-            scalarField& edgeWeights,
-            scalarField& invSumWeight
-        ) const;
-
 
         // Calculation of medial axis information
 
@@ -314,12 +289,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-#   include "medialAxisMeshMoverTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMoverTemplates.C b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMoverTemplates.C
deleted file mode 100644
index edd053f63f..0000000000
--- a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMoverTemplates.C
+++ /dev/null
@@ -1,93 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "medialAxisMeshMover.H"
-#include "syncTools.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-template<class Type>
-void Foam::medialAxisMeshMover::weightedSum
-(
-    const polyMesh& mesh,
-    const PackedBoolList& isMasterEdge,
-    const labelList& meshEdges,
-    const labelList& meshPoints,
-    const edgeList& edges,
-    const scalarField& edgeWeights,
-    const Field<Type>& pointData,
-    Field<Type>& sum
-)
-{
-    if
-    (
-        mesh.nEdges() != isMasterEdge.size()
-     || edges.size() != meshEdges.size()
-     || edges.size() != edgeWeights.size()
-     || meshPoints.size() != pointData.size()
-    )
-    {
-        FatalErrorIn("medialAxisMeshMover::weightedSum(..)")
-            << "Inconsistent sizes for edge or point data:"
-            << " meshEdges:" << meshEdges.size()
-            << " isMasterEdge:" << isMasterEdge.size()
-            << " edgeWeights:" << edgeWeights.size()
-            << " edges:" << edges.size()
-            << " pointData:" << pointData.size()
-            << " meshPoints:" << meshPoints.size()
-            << abort(FatalError);
-    }
-
-    sum.setSize(meshPoints.size());
-    sum = pTraits<Type>::zero;
-
-    forAll(edges, edgeI)
-    {
-        if (isMasterEdge.get(meshEdges[edgeI]) == 1)
-        {
-            const edge& e = edges[edgeI];
-
-            scalar eWeight = edgeWeights[edgeI];
-
-            label v0 = e[0];
-            label v1 = e[1];
-
-            sum[v0] += eWeight*pointData[v1];
-            sum[v1] += eWeight*pointData[v0];
-        }
-    }
-
-    syncTools::syncPointList
-    (
-        mesh,
-        meshPoints,
-        sum,
-        plusEqOp<Type>(),
-        pTraits<Type>::zero     // null value
-    );
-}
-
-
-// ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
index d93a5ca6d3..8d41eebd25 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
@@ -1949,6 +1949,64 @@ void Foam::meshRefinement::checkCoupledFaceZones(const polyMesh& mesh)
 }
 
 
+void Foam::meshRefinement::calculateEdgeWeights
+(
+    const polyMesh& mesh,
+    const PackedBoolList& isMasterEdge,
+    const labelList& meshEdges,
+    const labelList& meshPoints,
+    const edgeList& edges,
+    scalarField& edgeWeights,
+    scalarField& invSumWeight
+)
+{
+    const pointField& pts = mesh.points();
+
+    // Calculate edgeWeights and inverse sum of edge weights
+    edgeWeights.setSize(meshEdges.size());
+    invSumWeight.setSize(meshPoints.size());
+
+    forAll(edges, edgeI)
+    {
+        const edge& e = edges[edgeI];
+        scalar eMag = max
+        (
+            VSMALL,
+            mag
+            (
+                pts[meshPoints[e[1]]]
+              - pts[meshPoints[e[0]]]
+            )
+        );
+        edgeWeights[edgeI] = 1.0/eMag;
+    }
+
+    // Sum per point all edge weights
+    weightedSum
+    (
+        mesh,
+        isMasterEdge,
+        meshEdges,
+        meshPoints,
+        edges,
+        edgeWeights,
+        scalarField(meshPoints.size(), 1.0),  // data
+        invSumWeight
+    );
+
+    // Inplace invert
+    forAll(invSumWeight, pointI)
+    {
+        scalar w = invSumWeight[pointI];
+
+        if (w > 0.0)
+        {
+            invSumWeight[pointI] = 1.0/w;
+        }
+    }
+}
+
+
 Foam::label Foam::meshRefinement::appendPatch
 (
     fvMesh& mesh,
@@ -2167,9 +2225,7 @@ Foam::label Foam::meshRefinement::addMeshedPatch
 //        );
 
         // Store
-        label sz = meshedPatches_.size();
-        meshedPatches_.setSize(sz+1);
-        meshedPatches_[sz] = name;
+        meshedPatches_.append(name);
 
         return patchI;
     }
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index f8d6258ba7..57114bd740 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -729,6 +729,33 @@ public:
             //- Helper function: check that face zones are synced
             static void checkCoupledFaceZones(const polyMesh&);
 
+            //- Helper: calculate edge weights (1/length)
+            static void calculateEdgeWeights
+            (
+                const polyMesh& mesh,
+                const PackedBoolList& isMasterEdge,
+                const labelList& meshEdges,
+                const labelList& meshPoints,
+                const edgeList& edges,
+                scalarField& edgeWeights,
+                scalarField& invSumWeight
+            );
+
+            //- Helper: weighted sum (over all subset of mesh points) by
+            //  summing contribution from (master) edges
+            template<class Type>
+            static void weightedSum
+            (
+                const polyMesh& mesh,
+                const PackedBoolList& isMasterEdge,
+                const labelList& meshEdges,
+                const labelList& meshPoints,
+                const edgeList& edges,
+                const scalarField& edgeWeights,
+                const Field<Type>& data,
+                Field<Type>& sum
+            );
+
 
         // Refinement
 
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index a7084e9fc6..701c0dbbc0 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -458,10 +458,18 @@ void Foam::meshRefinement::markFeatureCellLevel
     // Database to pass into trackedParticle::move
     trackedParticle::trackingData td(startPointCloud, maxFeatureLevel);
 
+
     // Track all particles to their end position (= starting feature point)
     // Note that the particle might have started on a different processor
     // so this will transfer across nicely until we can start tracking proper.
     scalar maxTrackLen = 2.0*mesh_.bounds().mag();
+
+    if (debug&meshRefinement::FEATURESEEDS)
+    {
+        Pout<< "Tracking " << startPointCloud.size()
+            << " particles over distance " << maxTrackLen
+            << " to find the starting cell" << endl;
+    }
     startPointCloud.move(td, maxTrackLen);
 
 
@@ -485,6 +493,11 @@ void Foam::meshRefinement::markFeatureCellLevel
         IDLList<trackedParticle>()
     );
 
+    if (debug&meshRefinement::FEATURESEEDS)
+    {
+        Pout<< "Constructing cloud for cell marking" << endl;
+    }
+
     forAllIter(Cloud<trackedParticle>, startPointCloud, iter)
     {
         const trackedParticle& startTp = iter();
@@ -532,11 +545,14 @@ void Foam::meshRefinement::markFeatureCellLevel
     while (true)
     {
         // Track all particles to their end position.
+        if (debug&meshRefinement::FEATURESEEDS)
+        {
+            Pout<< "Tracking " << cloud.size()
+                << " particles over distance " << maxTrackLen
+                << " to mark cells" << endl;
+        }
         cloud.move(td, maxTrackLen);
 
-
-        label nParticles = 0;
-
         // Make particle follow edge.
         forAllIter(Cloud<trackedParticle>, cloud, iter)
         {
@@ -578,15 +594,15 @@ void Foam::meshRefinement::markFeatureCellLevel
                 // seeded. Delete particle.
                 cloud.deleteParticle(tp);
             }
-            else
-            {
-                // Keep particle
-                nParticles++;
-            }
         }
 
-        reduce(nParticles, sumOp<label>());
-        if (nParticles == 0)
+
+        if (debug&meshRefinement::FEATURESEEDS)
+        {
+            Pout<< "Remaining particles " << cloud.size() << endl;
+        }
+
+        if (returnReduce(cloud.size(), sumOp<label>()) == 0)
         {
             break;
         }
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C
index 45fcee4d5d..528e812197 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C
@@ -26,16 +26,12 @@ License
 #include "meshRefinement.H"
 #include "fvMesh.H"
 #include "globalIndex.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
+#include "syncTools.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 // Add a T entry
-template<class T> void meshRefinement::updateList
+template<class T> void Foam::meshRefinement::updateList
 (
     const labelList& newToOld,
     const T& nullValue,
@@ -59,7 +55,7 @@ template<class T> void meshRefinement::updateList
 
 
 template<class T>
-T meshRefinement::gAverage
+T Foam::meshRefinement::gAverage
 (
     const polyMesh& mesh,
     const PackedBoolList& isMasterElem,
@@ -109,7 +105,7 @@ T meshRefinement::gAverage
 
 
 template<class T>
-T meshRefinement::gAverage
+T Foam::meshRefinement::gAverage
 (
     const polyMesh& mesh,
     const PackedBoolList& isMasterElem,
@@ -162,7 +158,7 @@ T meshRefinement::gAverage
 
 // Compare two lists over all boundary faces
 template<class T>
-void meshRefinement::testSyncBoundaryFaceList
+void Foam::meshRefinement::testSyncBoundaryFaceList
 (
     const scalar tol,
     const string& msg,
@@ -222,7 +218,7 @@ void meshRefinement::testSyncBoundaryFaceList
 // Print list sorted by coordinates. Used for comparing non-parallel v.s.
 // parallel operation
 template<class T>
-void meshRefinement::collectAndPrint
+void Foam::meshRefinement::collectAndPrint
 (
     const UList<point>& points,
     const UList<T>& data
@@ -268,7 +264,11 @@ void meshRefinement::collectAndPrint
 
 //template<class T, class Mesh>
 template<class GeoField>
-void meshRefinement::addPatchFields(fvMesh& mesh, const word& patchFieldType)
+void Foam::meshRefinement::addPatchFields
+(
+    fvMesh& mesh,
+    const word& patchFieldType
+)
 {
     HashTable<GeoField*> flds
     (
@@ -298,7 +298,11 @@ void meshRefinement::addPatchFields(fvMesh& mesh, const word& patchFieldType)
 
 // Reorder patch field
 template<class GeoField>
-void meshRefinement::reorderPatchFields(fvMesh& mesh, const labelList& oldToNew)
+void Foam::meshRefinement::reorderPatchFields
+(
+    fvMesh& mesh,
+    const labelList& oldToNew
+)
 {
     HashTable<GeoField*> flds
     (
@@ -316,7 +320,11 @@ void meshRefinement::reorderPatchFields(fvMesh& mesh, const labelList& oldToNew)
 
 
 template<class Enum>
-int meshRefinement::readFlags(const Enum& namedEnum, const wordList& words)
+int Foam::meshRefinement::readFlags
+(
+    const Enum& namedEnum,
+    const wordList& words
+)
 {
     int flags = 0;
 
@@ -330,8 +338,66 @@ int meshRefinement::readFlags(const Enum& namedEnum, const wordList& words)
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+template<class Type>
+void Foam::meshRefinement::weightedSum
+(
+    const polyMesh& mesh,
+    const PackedBoolList& isMasterEdge,
+    const labelList& meshEdges,
+    const labelList& meshPoints,
+    const edgeList& edges,
+    const scalarField& edgeWeights,
+    const Field<Type>& pointData,
+    Field<Type>& sum
+)
+{
+    if
+    (
+        mesh.nEdges() != isMasterEdge.size()
+     || edges.size() != meshEdges.size()
+     || edges.size() != edgeWeights.size()
+     || meshPoints.size() != pointData.size()
+    )
+    {
+        FatalErrorIn("medialAxisMeshMover::weightedSum(..)")
+            << "Inconsistent sizes for edge or point data:"
+            << " meshEdges:" << meshEdges.size()
+            << " isMasterEdge:" << isMasterEdge.size()
+            << " edgeWeights:" << edgeWeights.size()
+            << " edges:" << edges.size()
+            << " pointData:" << pointData.size()
+            << " meshPoints:" << meshPoints.size()
+            << abort(FatalError);
+    }
+
+    sum.setSize(meshPoints.size());
+    sum = pTraits<Type>::zero;
+
+    forAll(edges, edgeI)
+    {
+        if (isMasterEdge.get(meshEdges[edgeI]) == 1)
+        {
+            const edge& e = edges[edgeI];
+
+            scalar eWeight = edgeWeights[edgeI];
+
+            label v0 = e[0];
+            label v1 = e[1];
+
+            sum[v0] += eWeight*pointData[v1];
+            sum[v1] += eWeight*pointData[v0];
+        }
+    }
+
+    syncTools::syncPointList
+    (
+        mesh,
+        meshPoints,
+        sum,
+        plusEqOp<Type>(),
+        pTraits<Type>::zero     // null value
+    );
+}
 
-} // End namespace Foam
 
 // ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
index a5109bc3c8..f2b730bc9d 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
+++ b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
@@ -26,6 +26,7 @@ License
 #include "refinementFeatures.H"
 #include "Time.H"
 #include "Tuple2.H"
+#include "DynamicField.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -41,7 +42,40 @@ void Foam::refinementFeatures::read
 
         fileName featFileName(dict.lookup("file"));
 
+
+        // Try reading extendedEdgeMesh first
+
+        IOobject extFeatObj
+        (
+            featFileName,                       // name
+            io.time().constant(),               // instance
+            "extendedFeatureEdgeMesh",          // local
+            io.time(),                          // registry
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        );
+
+        const fileName fName(extFeatObj.filePath());
+
+        if (!fName.empty() && extendedEdgeMesh::canRead(fName))
         {
+            autoPtr<extendedEdgeMesh> eMeshPtr = extendedEdgeMesh::New
+            (
+                fName
+            );
+
+            Info<< "Read extendedFeatureEdgeMesh " << extFeatObj.name()
+                << nl << incrIndent;
+            eMeshPtr().writeStats(Info);
+            Info<< decrIndent << endl;
+
+            set(featI, new extendedFeatureEdgeMesh(extFeatObj, eMeshPtr()));
+        }
+        else
+        {
+            // Try reading edgeMesh
+
             IOobject featObj
             (
                 featFileName,                       // name
@@ -53,21 +87,106 @@ void Foam::refinementFeatures::read
                 false
             );
 
-            autoPtr<edgeMesh> eMeshPtr = edgeMesh::New(featObj.filePath());
+            const fileName fName(featObj.filePath());
 
-            set
-            (
-                featI,
-                new featureEdgeMesh
+            if (fName.empty())
+            {
+                FatalIOErrorIn
                 (
-                    featObj,
-                    eMeshPtr->points(),
-                    eMeshPtr->edges()
-                )
+                    "refinementFeatures::read"
+                    "(const objectRegistry&"
+                    ", const PtrList<dictionary>&)",
+                    dict
+                )   << "Could not open " << featObj.objectPath()
+                    << exit(FatalIOError);
+            }
+
+
+            // Read as edgeMesh
+            autoPtr<edgeMesh> eMeshPtr = edgeMesh::New(fName);
+            const edgeMesh& eMesh = eMeshPtr();
+
+            Info<< "Read edgeMesh " << featObj.name() << nl
+                << incrIndent;
+            eMesh.writeStats(Info);
+            Info<< decrIndent << endl;
+
+
+            // Analyse for feature points. These are all classified as mixed
+            // points for lack of anything better
+            const labelListList& pointEdges = eMesh.pointEdges();
+
+            labelList oldToNew(eMesh.points().size(), -1);
+            DynamicField<point> newPoints(eMesh.points().size());
+            forAll(pointEdges, pointI)
+            {
+                if (pointEdges[pointI].size() > 2)
+                {
+                    oldToNew[pointI] = newPoints.size();
+                    newPoints.append(eMesh.points()[pointI]);
+                }
+                //else if (pointEdges[pointI].size() == 2)
+                //MEJ: do something based on a feature angle?
+            }
+            label nFeatures = newPoints.size();
+            forAll(oldToNew, pointI)
+            {
+                if (oldToNew[pointI] == -1)
+                {
+                    oldToNew[pointI] = newPoints.size();
+                    newPoints.append(eMesh.points()[pointI]);
+                }
+            }
+
+
+            const edgeList& edges = eMesh.edges();
+            edgeList newEdges(edges.size());
+            forAll(edges, edgeI)
+            {
+                const edge& e = edges[edgeI];
+                newEdges[edgeI] = edge
+                (
+                    oldToNew[e[0]],
+                    oldToNew[e[1]]
+                );
+            }
+
+            // Construct an extendedEdgeMesh with
+            // - all points on more than 2 edges : mixed feature points
+            // - all edges : external edges
+
+            extendedEdgeMesh eeMesh
+            (
+                newPoints,          // pts
+                newEdges,           // eds
+                0,                  // (point) concaveStart
+                0,                  // (point) mixedStart
+                nFeatures,          // (point) nonFeatureStart
+                edges.size(),       // (edge) internalStart
+                edges.size(),       // (edge) flatStart
+                edges.size(),       // (edge) openStart
+                edges.size(),       // (edge) multipleStart
+                vectorField(0),     // normals
+                List<extendedEdgeMesh::sideVolumeType>(0),// normalVolumeTypes
+                vectorField(0),     // edgeDirections
+                labelListList(0),   // normalDirections
+                labelListList(0),   // edgeNormals
+                labelListList(0),   // featurePointNormals
+                labelListList(0),   // featurePointEdges
+                labelList(0)        // regionEdges
             );
+
+
+            Info<< "Constructed extendedFeatureEdgeMesh " << featObj.name()
+                << nl << incrIndent;
+            eeMesh.writeStats(Info);
+            Info<< decrIndent << endl;
+
+
+            set(featI, new extendedFeatureEdgeMesh(featObj, eeMesh));
         }
 
-        const featureEdgeMesh& eMesh = operator[](featI);
+        const edgeMesh& eMesh = operator[](featI);
 
         //eMesh.mergePoints(meshRefiner_.mergeDistance());
 
@@ -145,7 +264,7 @@ void Foam::refinementFeatures::buildTrees
     const labelList& featurePoints
 )
 {
-    const featureEdgeMesh& eMesh = operator[](featI);
+    const edgeMesh& eMesh = operator[](featI);
     const pointField& points = eMesh.points();
     const edgeList& edges = eMesh.edges();
 
@@ -277,7 +396,7 @@ Foam::refinementFeatures::refinementFeatures
     const PtrList<dictionary>& featDicts
 )
 :
-    PtrList<featureEdgeMesh>(featDicts.size()),
+    PtrList<extendedFeatureEdgeMesh>(featDicts.size()),
     distances_(featDicts.size()),
     levels_(featDicts.size()),
     edgeTrees_(featDicts.size()),
@@ -289,7 +408,7 @@ Foam::refinementFeatures::refinementFeatures
     // Search engines
     forAll(*this, i)
     {
-        const featureEdgeMesh& eMesh = operator[](i);
+        const extendedEdgeMesh& eMesh = operator[](i);
         const labelListList& pointEdges = eMesh.pointEdges();
 
         DynamicList<label> featurePoints;
@@ -303,81 +422,82 @@ Foam::refinementFeatures::refinementFeatures
 
         Info<< "Detected " << featurePoints.size()
             << " featurePoints out of " << pointEdges.size()
-            << " points on feature " << eMesh.name() << endl;
+            << " points on feature " << operator[](i).name()
+            << endl;
 
         buildTrees(i, featurePoints);
     }
 }
 
 
-Foam::refinementFeatures::refinementFeatures
-(
-    const objectRegistry& io,
-    const PtrList<dictionary>& featDicts,
-    const scalar minCos
-)
-:
-    PtrList<featureEdgeMesh>(featDicts.size()),
-    distances_(featDicts.size()),
-    levels_(featDicts.size()),
-    edgeTrees_(featDicts.size()),
-    pointTrees_(featDicts.size())
-{
-    // Read features
-    read(io, featDicts);
-
-    // Search engines
-    forAll(*this, i)
-    {
-        const featureEdgeMesh& eMesh = operator[](i);
-        const pointField& points = eMesh.points();
-        const edgeList& edges = eMesh.edges();
-        const labelListList& pointEdges = eMesh.pointEdges();
-
-        DynamicList<label> featurePoints;
-        forAll(pointEdges, pointI)
-        {
-            const labelList& pEdges = pointEdges[pointI];
-            if (pEdges.size() > 2)
-            {
-                featurePoints.append(pointI);
-            }
-            else if (pEdges.size() == 2)
-            {
-                // Check the angle
-                const edge& e0 = edges[pEdges[0]];
-                const edge& e1 = edges[pEdges[1]];
-
-                const point& p = points[pointI];
-                const point& p0 = points[e0.otherVertex(pointI)];
-                const point& p1 = points[e1.otherVertex(pointI)];
-
-                vector v0 = p-p0;
-                scalar v0Mag = mag(v0);
-
-                vector v1 = p1-p;
-                scalar v1Mag = mag(v1);
-
-                if
-                (
-                    v0Mag > SMALL
-                 && v1Mag > SMALL
-                 && ((v0/v0Mag & v1/v1Mag) < minCos)
-                )
-                {
-                    featurePoints.append(pointI);
-                }
-            }
-        }
-
-        Info<< "Detected " << featurePoints.size()
-            << " featurePoints out of " << points.size()
-            << " points on feature " << eMesh.name()
-            << " when using feature cos " << minCos << endl;
-
-        buildTrees(i, featurePoints);
-    }
-}
+//Foam::refinementFeatures::refinementFeatures
+//(
+//    const objectRegistry& io,
+//    const PtrList<dictionary>& featDicts,
+//    const scalar minCos
+//)
+//:
+//    PtrList<extendedFeatureEdgeMesh>(featDicts.size()),
+//    distances_(featDicts.size()),
+//    levels_(featDicts.size()),
+//    edgeTrees_(featDicts.size()),
+//    pointTrees_(featDicts.size())
+//{
+//    // Read features
+//    read(io, featDicts);
+//
+//    // Search engines
+//    forAll(*this, i)
+//    {
+//        const edgeMesh& eMesh = operator[](i);
+//        const pointField& points = eMesh.points();
+//        const edgeList& edges = eMesh.edges();
+//        const labelListList& pointEdges = eMesh.pointEdges();
+//
+//        DynamicList<label> featurePoints;
+//        forAll(pointEdges, pointI)
+//        {
+//            const labelList& pEdges = pointEdges[pointI];
+//            if (pEdges.size() > 2)
+//            {
+//                featurePoints.append(pointI);
+//            }
+//            else if (pEdges.size() == 2)
+//            {
+//                // Check the angle
+//                const edge& e0 = edges[pEdges[0]];
+//                const edge& e1 = edges[pEdges[1]];
+//
+//                const point& p = points[pointI];
+//                const point& p0 = points[e0.otherVertex(pointI)];
+//                const point& p1 = points[e1.otherVertex(pointI)];
+//
+//                vector v0 = p-p0;
+//                scalar v0Mag = mag(v0);
+//
+//                vector v1 = p1-p;
+//                scalar v1Mag = mag(v1);
+//
+//                if
+//                (
+//                    v0Mag > SMALL
+//                 && v1Mag > SMALL
+//                 && ((v0/v0Mag & v1/v1Mag) < minCos)
+//                )
+//                {
+//                    featurePoints.append(pointI);
+//                }
+//            }
+//        }
+//
+//        Info<< "Detected " << featurePoints.size()
+//            << " featurePoints out of " << points.size()
+//            << " points on feature " << i   //eMesh.name()
+//            << " when using feature cos " << minCos << endl;
+//
+//        buildTrees(i, featurePoints);
+//    }
+//}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -387,12 +507,16 @@ void Foam::refinementFeatures::findNearestEdge
     const pointField& samples,
     const scalarField& nearestDistSqr,
     labelList& nearFeature,
-    List<pointIndexHit>& nearInfo
+    List<pointIndexHit>& nearInfo,
+    vectorField& nearNormal
 ) const
 {
     nearFeature.setSize(samples.size());
     nearFeature = -1;
     nearInfo.setSize(samples.size());
+    nearInfo = pointIndexHit();
+    nearNormal.setSize(samples.size());
+    nearNormal = vector::zero;
 
     forAll(edgeTrees_, featI)
     {
@@ -418,8 +542,18 @@ void Foam::refinementFeatures::findNearestEdge
 
                 if (info.hit())
                 {
-                    nearInfo[sampleI] = info;
                     nearFeature[sampleI] = featI;
+                    nearInfo[sampleI] = pointIndexHit
+                    (
+                        info.hit(),
+                        info.hitPoint(),
+                        tree.shapes().edgeLabels()[info.index()]
+                    );
+
+                    const treeDataEdge& td = tree.shapes();
+                    const edge& e = td.edges()[nearInfo[sampleI].index()];
+                    nearNormal[sampleI] =  e.vec(td.points());
+                    nearNormal[sampleI] /= mag(nearNormal[sampleI])+VSMALL;
                 }
             }
         }
@@ -427,18 +561,71 @@ void Foam::refinementFeatures::findNearestEdge
 }
 
 
+//void Foam::refinementFeatures::findNearestPoint
+//(
+//    const pointField& samples,
+//    const scalarField& nearestDistSqr,
+//    labelList& nearFeature,
+//    labelList& nearIndex
+//) const
+//{
+//    nearFeature.setSize(samples.size());
+//    nearFeature = -1;
+//    nearIndex.setSize(samples.size());
+//    nearIndex = -1;
+//
+//    forAll(pointTrees_, featI)
+//    {
+//        const indexedOctree<treeDataPoint>& tree = pointTrees_[featI];
+//
+//        if (tree.shapes().pointLabels().size() > 0)
+//        {
+//            forAll(samples, sampleI)
+//            {
+//                const point& sample = samples[sampleI];
+//
+//                scalar distSqr;
+//                if (nearFeature[sampleI] != -1)
+//                {
+//                    label nearFeatI = nearFeature[sampleI];
+//                    const indexedOctree<treeDataPoint>& nearTree =
+//                        pointTrees_[nearFeatI];
+//                    label featPointI =
+//                        nearTree.shapes().pointLabels()[nearIndex[sampleI]];
+//                    const point& featPt =
+//                        operator[](nearFeatI).points()[featPointI];
+//                    distSqr = magSqr(featPt-sample);
+//                }
+//                else
+//                {
+//                    distSqr = nearestDistSqr[sampleI];
+//                }
+//
+//                pointIndexHit info = tree.findNearest(sample, distSqr);
+//
+//                if (info.hit())
+//                {
+//                    nearFeature[sampleI] = featI;
+//                    nearIndex[sampleI] = info.index();
+//                }
+//            }
+//        }
+//    }
+//}
+
+
 void Foam::refinementFeatures::findNearestPoint
 (
     const pointField& samples,
     const scalarField& nearestDistSqr,
     labelList& nearFeature,
-    labelList& nearIndex
+    List<pointIndexHit>& nearInfo
 ) const
 {
     nearFeature.setSize(samples.size());
     nearFeature = -1;
-    nearIndex.setSize(samples.size());
-    nearIndex = -1;
+    nearInfo.setSize(samples.size());
+    nearInfo = pointIndexHit();
 
     forAll(pointTrees_, featI)
     {
@@ -453,14 +640,7 @@ void Foam::refinementFeatures::findNearestPoint
                 scalar distSqr;
                 if (nearFeature[sampleI] != -1)
                 {
-                    label nearFeatI = nearFeature[sampleI];
-                    const indexedOctree<treeDataPoint>& nearTree =
-                        pointTrees_[nearFeatI];
-                    label featPointI =
-                        nearTree.shapes().pointLabels()[nearIndex[sampleI]];
-                    const point& featPt =
-                        operator[](nearFeatI).points()[featPointI];
-                    distSqr = magSqr(featPt-sample);
+                    distSqr = magSqr(nearInfo[sampleI].hitPoint()-sample);
                 }
                 else
                 {
@@ -472,7 +652,12 @@ void Foam::refinementFeatures::findNearestPoint
                 if (info.hit())
                 {
                     nearFeature[sampleI] = featI;
-                    nearIndex[sampleI] = info.index();
+                    nearInfo[sampleI] = pointIndexHit
+                    (
+                        info.hit(),
+                        info.hitPoint(),
+                        tree.shapes().pointLabels()[info.index()]
+                    );
                 }
             }
         }
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H
index aa2807bac1..4cfaf20ffa 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H
+++ b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,7 +35,7 @@ SourceFiles
 #ifndef refinementFeatures_H
 #define refinementFeatures_H
 
-#include "featureEdgeMesh.H"
+#include "extendedFeatureEdgeMesh.H"
 #include "indexedOctree.H"
 #include "treeDataEdge.H"
 #include "treeDataPoint.H"
@@ -51,7 +51,7 @@ namespace Foam
 
 class refinementFeatures
 :
-    public PtrList<featureEdgeMesh>
+    public PtrList<extendedFeatureEdgeMesh>
 {
 private:
 
@@ -86,6 +86,20 @@ private:
             labelList& maxLevel
         ) const;
 
+
+protected:
+
+        const PtrList<indexedOctree<treeDataEdge> >& edgeTrees() const
+        {
+            return edgeTrees_;
+        }
+
+        const PtrList<indexedOctree<treeDataPoint> >& pointTrees() const
+        {
+            return pointTrees_;
+        }
+
+
 public:
 
     // Constructors
@@ -97,15 +111,6 @@ public:
             const PtrList<dictionary>& featDicts
         );
 
-        //- Construct from description and do geometric analysis to determine
-        //  feature points
-        refinementFeatures
-        (
-            const objectRegistry& io,
-            const PtrList<dictionary>& featDicts,
-            const scalar minCos
-        );
-
 
     // Member Functions
 
@@ -123,42 +128,38 @@ public:
                 return distances_;
             }
 
-            const PtrList<indexedOctree<treeDataEdge> >& edgeTrees() const
-            {
-                return edgeTrees_;
-            }
-
-            const PtrList<indexedOctree<treeDataPoint> >& pointTrees() const
-            {
-                return pointTrees_;
-            }
-
 
         // Query
 
             //- Highest distance of all features
             scalar maxDistance() const;
 
-            //- Find nearest point on nearest feature edge
+            //- Find nearest point on nearest feature edge. Sets
+            //  - nearFeature: index of feature mesh
+            //  - nearInfo   : location on feature edge and edge index
+            //                 (note: not feature edge index but index into
+            //                  edges() directly)
+            //  - nearNormal : local feature edge normal
             void findNearestEdge
             (
                 const pointField& samples,
                 const scalarField& nearestDistSqr,
                 labelList& nearFeature,
-                List<pointIndexHit>& nearInfo
+                List<pointIndexHit>& nearInfo,
+                vectorField& nearNormal
             ) const;
 
-            //- Find nearest feature point. Is an index into feature points
-            //  which itself is an index into the edgeMesh points.
-            //  So the point index is
-            //      pointTrees()[nearFeature].shapes().pointLabels()[nearIndex]
-            //  Wip.
+            //- Find nearest feature point. Sets
+            //  - nearFeature: index of feature mesh
+            //  - nearInfo   : location on feature point and point index.
+            //                 (note: not index into shapes().pointLabels() but
+            //                  index into points() directly)
             void findNearestPoint
             (
                 const pointField& samples,
                 const scalarField& nearestDistSqr,
                 labelList& nearFeature,
-                labelList& nearIndex
+                List<pointIndexHit>& nearInfo
             ) const;
 
             //- Find shell level higher than ptLevel
diff --git a/tutorials/incompressible/icoFoam/Allclean b/tutorials/incompressible/icoFoam/Allclean
index 79a05d40a4..587ffb1d53 100755
--- a/tutorials/incompressible/icoFoam/Allclean
+++ b/tutorials/incompressible/icoFoam/Allclean
@@ -30,4 +30,7 @@ do
     removeCase $caseName
 done
 
+
+(cd elbow && ./Allclean)
+
 # ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/incompressible/icoFoam/Allrun b/tutorials/incompressible/icoFoam/Allrun
index 44d04afe47..acdcaabaee 100755
--- a/tutorials/incompressible/icoFoam/Allrun
+++ b/tutorials/incompressible/icoFoam/Allrun
@@ -18,12 +18,6 @@ runMapFieldsConsistent()
     mapFieldsNew $1 -case $2 -sourceTime latestTime -consistent > $2/log.mapFields 2>&1
 }
 
-runFluentMeshToFoam()
-{
-    echo "fluentMeshToFoam: converting mesh $2"
-    fluentMeshToFoam $2 -case $1 > $1/log.fluentMeshToFoam 2>&1
-}
-
 copySolutionDirs()
 {
     echo "Copying $2/0* directory to $1"
@@ -95,16 +89,8 @@ do
     ( cd $caseName && runApplication `getApplication` )
 done
 
+
 # elbow case for testing Fluent-FOAM conversion tools
-
-runFluentMeshToFoam elbow elbow/elbow.msh
-
-(
-    cd elbow || exit
-
-    runApplication `getApplication`
-    runApplication foamMeshToFluent
-    runApplication foamDataToFluent
-)
+(cd elbow && ./Allrun)
 
 # ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/incompressible/icoFoam/elbow/Allclean b/tutorials/incompressible/icoFoam/elbow/Allclean
new file mode 100755
index 0000000000..02a68a6b17
--- /dev/null
+++ b/tutorials/incompressible/icoFoam/elbow/Allclean
@@ -0,0 +1,11 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+rm -f constant/polyMesh/boundary
+rm -rf fluentInterface
+cleanCase
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/incompressible/icoFoam/elbow/Allrun b/tutorials/incompressible/icoFoam/elbow/Allrun
new file mode 100755
index 0000000000..07a1ffac39
--- /dev/null
+++ b/tutorials/incompressible/icoFoam/elbow/Allrun
@@ -0,0 +1,15 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+# Get application directory
+application=`getApplication`
+
+runApplication fluentMeshToFoam elbow.msh
+runApplication "$application"
+runApplication foamMeshToFluent
+runApplication foamDataToFluent
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/incompressible/icoFoam/elbow/constant/polyMesh/boundary b/tutorials/incompressible/icoFoam/elbow/constant/polyMesh/boundary
deleted file mode 100644
index fbdbc176b3..0000000000
--- a/tutorials/incompressible/icoFoam/elbow/constant/polyMesh/boundary
+++ /dev/null
@@ -1,58 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       polyBoundaryMesh;
-    location    "constant/polyMesh";
-    object      boundary;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-6
-(
-    wall-4
-    {
-        type            wall;
-        nFaces          100;
-        startFace       1300;
-    }
-    velocity-inlet-5
-    {
-        type            patch;
-        nFaces          8;
-        startFace       1400;
-    }
-    velocity-inlet-6
-    {
-        type            patch;
-        nFaces          4;
-        startFace       1408;
-    }
-    pressure-outlet-7
-    {
-        type            patch;
-        nFaces          8;
-        startFace       1412;
-    }
-    wall-8
-    {
-        type            wall;
-        nFaces          34;
-        startFace       1420;
-    }
-    frontAndBackPlanes
-    {
-        type            empty;
-        nFaces          1836;
-        startFace       1454;
-    }
-)
-
-// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/pitzDailyExptInlet/0/epsilon b/tutorials/incompressible/simpleFoam/pitzDailyExptInlet/0/epsilon
index b70f15f119..01dccefae7 100644
--- a/tutorials/incompressible/simpleFoam/pitzDailyExptInlet/0/epsilon
+++ b/tutorials/incompressible/simpleFoam/pitzDailyExptInlet/0/epsilon
@@ -26,41 +26,6 @@ boundaryField
         type            timeVaryingMappedFixedValue;
         setAverage      0;
         offset          0;
-        value           nonuniform List<scalar>
-30
-(
-9813.84
-8665.24
-1866.31
-755.118
-205.654
-76.6694
-28.4518
-16.0868
-15.9867
-11.0187
-7.95753
-5.26064
-3.44136
-2.55317
-2.27183
-2.33608
-2.9115
-3.59492
-3.0497
-2.716
-2.9325
-3.88456
-6.91821
-14.9754
-37.5461
-217.022
-2043.58
-4864.22
-6244
-6334.7
-)
-;
     }
     outlet
     {
diff --git a/tutorials/incompressible/simpleFoam/pitzDailyExptInlet/0/k b/tutorials/incompressible/simpleFoam/pitzDailyExptInlet/0/k
index f5f11a1aed..c0fa4d80d4 100644
--- a/tutorials/incompressible/simpleFoam/pitzDailyExptInlet/0/k
+++ b/tutorials/incompressible/simpleFoam/pitzDailyExptInlet/0/k
@@ -26,41 +26,6 @@ boundaryField
         type            timeVaryingMappedFixedValue;
         setAverage      0;
         offset          0;
-        value           nonuniform List<scalar>
-30
-(
-2.95219
-2.54219
-0.725449
-0.486465
-0.353566
-0.240375
-0.172984
-0.147052
-0.146827
-0.135658
-0.12147
-0.0942189
-0.0833465
-0.0828453
-0.0955983
-0.0920838
-0.0967682
-0.0990811
-0.100866
-0.101556
-0.0967155
-0.0841739
-0.0904567
-0.130411
-0.194046
-0.219327
-0.975528
-2.22578
-3.12421
-2.28104
-)
-;
     }
     outlet
     {