diff --git a/applications/solvers/Lagrangian/kinematicParcelFoam/Make/files b/applications/solvers/Lagrangian/kinematicParcelFoam/Make/files
index f98171fc93..a63fc64fa7 100644
--- a/applications/solvers/Lagrangian/kinematicParcelFoam/Make/files
+++ b/applications/solvers/Lagrangian/kinematicParcelFoam/Make/files
@@ -1,3 +1,3 @@
 kinematicParcelFoam.C
 
-EXE = $(FOAM_USER_APPBIN)/kinematicParcelFoam
+EXE = $(FOAM_APPBIN)/kinematicParcelFoam
diff --git a/applications/solvers/Lagrangian/kinematicParcelFoam/kinematicParcelFoam.C b/applications/solvers/Lagrangian/kinematicParcelFoam/kinematicParcelFoam.C
index 290480583c..57b3dc5d13 100644
--- a/applications/solvers/Lagrangian/kinematicParcelFoam/kinematicParcelFoam.C
+++ b/applications/solvers/Lagrangian/kinematicParcelFoam/kinematicParcelFoam.C
@@ -26,7 +26,8 @@ Application
     kinematicParcelFoam
 
 Description
-    Transient solver a single kinematicCloud.
+    Transient solver for a single kinematicCloud. Uses precalculated velocity
+    field to evolve a cloud.
 
 \*---------------------------------------------------------------------------*/
 
diff --git a/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L b/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L
index f0c8e0a6c4..a8cb5dbb98 100644
--- a/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L
+++ b/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L
@@ -697,7 +697,7 @@ endOfSection               {space}")"{space}
 
  /* ------              Ignore remaining space and \n s.             ------  */
 
-<*>{some_space}|\n {
+<*>{some_space}|\n|\r {
     }
 
 
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/writeSurfFields.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/writeSurfFields.C
index 5e40ed5871..58107bc776 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/writeSurfFields.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/writeSurfFields.C
@@ -28,7 +28,7 @@ License
 #include "OFstream.H"
 #include "floatScalar.H"
 #include "writeFuns.H"
-#include "emptyFvPatchFields.H"
+#include "emptyFvsPatchFields.H"
 #include "fvsPatchFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -100,7 +100,7 @@ void writeSurfFields
 
             const fvPatch& pp = mesh.boundary()[patchI];
 
-            if (isA<emptyFvPatchVectorField>(pf))
+            if (isA<emptyFvsPatchVectorField>(pf))
             {
                 // Note: loop over polypatch size, not fvpatch size.
                 forAll(pp.patch(), i)
diff --git a/applications/utilities/surface/surfaceSubset/surfaceSubset.C b/applications/utilities/surface/surfaceSubset/surfaceSubset.C
index cf0771d260..aa30790b5c 100644
--- a/applications/utilities/surface/surfaceSubset/surfaceSubset.C
+++ b/applications/utilities/surface/surfaceSubset/surfaceSubset.C
@@ -248,7 +248,7 @@ int main(int argc, char *argv[])
         indexedOctree<treeDataTriSurface> selectTree
         (
             treeDataTriSurface(selectSurf),
-            bb.extend(rndGen, 1E-3),    // slightly randomize bb
+            bb.extend(rndGen, 1E-4),    // slightly randomize bb
             8,      // maxLevel
             10,     // leafsize
             3.0     // duplicity
diff --git a/etc/apps/paraview3/bashrc b/etc/apps/paraview3/bashrc
index 56cdd8aacf..b8da8ba698 100644
--- a/etc/apps/paraview3/bashrc
+++ b/etc/apps/paraview3/bashrc
@@ -45,9 +45,9 @@ export ParaView_INST_DIR=$WM_THIRD_PARTY_DIR/ParaView$ParaView_VERSION
 export ParaView_DIR=$ParaView_INST_DIR/platforms/$WM_ARCH$WM_COMPILER
 
 if [ "$PYTHONPATH" ]; then
-    export PYTHONPATH=$PYTHONPATH:$ParaView_DIR/Utilities/VTKPythonWrapping
+    export PYTHONPATH=$PYTHONPATH:$ParaView_DIR/Utilities/VTKPythonWrapping:$ParaView_DIR/lib/paraview-3.3
 else
-    export PYTHONPATH=$ParaView_DIR/Utilities/VTKPythonWrapping
+    export PYTHONPATH=$ParaView_DIR/Utilities/VTKPythonWrapping:$ParaView_DIR/lib/paraview-3.3
 fi
 
 
diff --git a/etc/apps/paraview3/cshrc b/etc/apps/paraview3/cshrc
index c9fd98f4a2..636b10eb75 100644
--- a/etc/apps/paraview3/cshrc
+++ b/etc/apps/paraview3/cshrc
@@ -45,9 +45,9 @@ setenv ParaView_INST_DIR $WM_THIRD_PARTY_DIR/ParaView$ParaView_VERSION
 setenv ParaView_DIR $ParaView_INST_DIR/platforms/$WM_ARCH$WM_COMPILER
 
 if ($?PYTHONPATH) then
-    setenv PYTHONPATH ${PYTHONPATH}:$ParaView_DIR/bin:$ParaView_DIR/Utilities/VTKPythonWrapping
+    setenv PYTHONPATH ${PYTHONPATH}:$ParaView_DIR/Utilities/VTKPythonWrapping:$ParaView_DIR/lib/paraview-3.3
 else
-    setenv PYTHONPATH $ParaView_DIR/bin:$ParaView_DIR/Utilities/VTKPythonWrapping
+    setenv PYTHONPATH $ParaView_DIR/Utilities/VTKPythonWrapping:$ParaView_DIR/lib/paraview-3.3
 endif
 
 if ( -r $ParaView_INST_DIR ) then
diff --git a/src/OpenFOAM/db/Time/TimeIO.C b/src/OpenFOAM/db/Time/TimeIO.C
index c74b674389..4299d05b5e 100644
--- a/src/OpenFOAM/db/Time/TimeIO.C
+++ b/src/OpenFOAM/db/Time/TimeIO.C
@@ -86,13 +86,6 @@ void Foam::Time::readDict()
 
             purgeWrite_ = 0;
         }
-
-        if (writeControl_ != wcTimeStep && purgeWrite_ > 0)
-        {
-            FatalIOErrorIn("Time::readDict()", controlDict_)
-                << "writeControl must be set to timeStep for purgeWrite "
-                << exit(FatalIOError);
-        }
     }
 
     if (controlDict_.found("timeFormat"))
diff --git a/src/OpenFOAM/meshes/boundBox/boundBox.C b/src/OpenFOAM/meshes/boundBox/boundBox.C
index c368e1a662..c87e552c85 100644
--- a/src/OpenFOAM/meshes/boundBox/boundBox.C
+++ b/src/OpenFOAM/meshes/boundBox/boundBox.C
@@ -87,9 +87,47 @@ boundBox::boundBox(Istream& is)
 
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
-Ostream& operator<<(Ostream& os, const boundBox& b)
+Ostream& operator<<(Ostream& os, const boundBox& bb)
 {
-    return os << b.min() << token::SPACE << b.max();
+    if (os.format() == IOstream::ASCII)
+    {
+        os << bb.min_ << token::SPACE << bb.max_;
+    }
+    else
+    {
+        os.write
+        (
+            reinterpret_cast<const char*>(&bb.min_),
+            sizeof(boundBox)
+        );
+    }
+
+    // Check state of Ostream
+    os.check("Ostream& operator<<(Ostream&, const boundBox&)");
+
+    return os;
+}
+
+
+Istream& operator>>(Istream& is, boundBox& bb)
+{
+    if (is.format() == IOstream::ASCII)
+    {
+        return is >> bb.min_ >>  bb.max_;
+    }
+    else
+    {
+        is.read
+        (
+            reinterpret_cast<char*>(&bb.min_),
+            sizeof(boundBox)
+        );
+    }
+
+    // Check state of Istream
+    is.check("Istream& operator>>(Istream&, boundBox&)");
+
+    return is;
 }
 
 
diff --git a/src/OpenFOAM/meshes/boundBox/boundBox.H b/src/OpenFOAM/meshes/boundBox/boundBox.H
index 40287e53eb..1493c0fe2f 100644
--- a/src/OpenFOAM/meshes/boundBox/boundBox.H
+++ b/src/OpenFOAM/meshes/boundBox/boundBox.H
@@ -153,12 +153,31 @@ public:
             }
 
 
-    // Ostream operator
+    // Friend Operators
 
-        friend Ostream& operator<<(Ostream& os, const boundBox& b);
+        friend bool operator==(const boundBox& a, const boundBox& b)
+        {
+            return (a.min_ == b.min_) && (a.max_ == b.max_);
+        }
+
+        friend bool operator!=(const boundBox& a, const boundBox& b)
+        {
+            return !(a == b);
+        }
+
+
+    // IOstream operator
+
+        friend Istream& operator>>(Istream& is, boundBox&);
+        friend Ostream& operator<<(Ostream& os, const boundBox&);
 };
 
 
+//- Specify data associated with boundBox type is contiguous
+template<>
+inline bool contiguous<boundBox>() {return contiguous<point>();}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C b/src/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C
index fe64ea81c3..1a577acecb 100644
--- a/src/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C
+++ b/src/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C
@@ -30,6 +30,7 @@ License
 #include "triSurfaceMesh.H"
 #include "refinementSurfaces.H"
 #include "searchableSurfaces.H"
+#include "orientedSurface.H"
 #include "pointIndexHit.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -211,7 +212,27 @@ void Foam::shellSurfaces::orient()
                     refCast<const triSurfaceMesh>(s)
                 );
 
-                refinementSurfaces::orientSurface(outsidePt, shell);
+                // Flip surface so outsidePt is outside.
+                bool anyFlipped = orientedSurface::orient
+                (
+                    shell,
+                    outsidePt,
+                    true
+                );
+
+                if (anyFlipped)
+                {
+                    // orientedSurface will have done a clearOut of the surface.
+                    // we could do a clearout of the triSurfaceMeshes::trees()
+                    // but these aren't affected by orientation
+                    // (except for cached
+                    // sideness which should not be set at this point.
+                    // !!Should check!)
+
+                    Info<< "shellSurfaces : Flipped orientation of surface "
+                        << s.name()
+                        << " so point " << outsidePt << " is outside." << endl;
+                }
             }
         }
     }
diff --git a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C b/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C
index 6d97285d73..f3918578cd 100644
--- a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C
+++ b/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C
@@ -29,8 +29,8 @@ License
 #include "matchPoints.H"
 #include "indirectPrimitivePatch.H"
 #include "meshTools.H"
-#include "octreeDataFace.H"
-#include "octree.H"
+#include "treeDataFace.H"
+#include "indexedOctree.H"
 #include "OFstream.H"
 #include "IndirectList.H"
 
@@ -1011,19 +1011,29 @@ void Foam::faceCoupleInfo::findSlavesCoveringMaster
 )
 {
     // Construct octree from all mesh0 boundary faces
-    octreeDataFace shapes(mesh0);
+    labelList bndFaces(mesh0.nFaces()-mesh0.nInternalFaces());
+    forAll(bndFaces, i)
+    {
+        bndFaces[i] = mesh0.nInternalFaces() + i;
+    }
 
     treeBoundBox overallBb(mesh0.points());
 
-    octree<octreeDataFace> tree
-    (
-        overallBb,  // overall search domain
-        shapes,     // all information needed to do checks on cells
-        1,          // min levels
-        20.0,       // maximum ratio of cubes v.s. cells
-        10.0
-    );
+    Random rndGen(123456);
 
+    indexedOctree<treeDataFace> tree
+    (
+        treeDataFace    // all information needed to search faces
+        (
+            false,                      // do not cache bb
+            mesh0,
+            bndFaces                    // boundary faces only
+        ),
+        overallBb.extend(rndGen, 1E-4), // overall search domain
+        8,                              // maxLevel
+        10,                             // leafsize
+        3.0                             // duplicity
+    );
 
     if (debug)
     {
@@ -1048,17 +1058,11 @@ void Foam::faceCoupleInfo::findSlavesCoveringMaster
         // Generate face centre (prevent cellCentres() reconstruction)
         point fc(f1.centre(mesh1.points()));
 
-        // Search in bounding box of face only.
-        treeBoundBox tightest(static_cast<const pointField&>(f1.points(mesh1.points())));
+        pointIndexHit nearInfo = tree.findNearest(fc, Foam::sqr(absTol));
 
-        scalar tightestDist = GREAT;
-
-        label index = tree.findNearest(fc, tightest, tightestDist);
-
-
-        if (index != -1)
+        if (nearInfo.hit())
         {
-            label mesh0FaceI = shapes.meshFaces()[index];
+            label mesh0FaceI = tree.shapes().faceLabels()[nearInfo.index()];
 
             // Check if points of f1 actually lie on top of mesh0 face
             // This is the bit that might fail since if f0 is severely warped
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C
index 1b3bf6e19b..0c567ee312 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C
@@ -62,6 +62,7 @@ namespace Foam
     };
 }
 
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::hexRef8::reorder
@@ -80,7 +81,7 @@ void Foam::hexRef8::reorder
 
         if (newI >= len)
         {
-            FatalErrorIn("hexRef8::reorder") << abort(FatalError);
+            FatalErrorIn("hexRef8::reorder(..)") << abort(FatalError);
         }
 
         if (newI >= 0)
@@ -557,22 +558,11 @@ Foam::label Foam::hexRef8::getAnchorCell
         }
 
         // Problem.
-
-        // Pick up points of the cell
-        const labelList cPoints(cellPoints(cellI));
-
-        Perr<< "cell:" << cellI << " points:" << endl;
-        forAll(cPoints, i)
-        {
-            label pointI = cPoints[i];
-
-            Perr<< "    " << pointI << " coord:" << mesh_.points()[pointI]
-                << nl;
-        }
+        dumpCell(cellI);
         Perr<< "cell:" << cellI << " anchorPoints:" << cellAnchorPoints[cellI]
             << endl;
 
-        FatalErrorIn("hexRef8::getAnchorCell")
+        FatalErrorIn("hexRef8::getAnchorCell(..)")
             << "Could not find point " << pointI
             << " in the anchorPoints for cell " << cellI << endl
             << "Does your original mesh obey the 2:1 constraint and"
@@ -690,9 +680,50 @@ Foam::label Foam::hexRef8::countAnchors
 }
 
 
+void Foam::hexRef8::dumpCell(const label cellI) const
+{
+    OFstream str(mesh_.time().path()/"cell_" + Foam::name(cellI) + ".obj");
+    Pout<< "hexRef8 : Dumping cell as obj to " << str.name() << endl;
+
+    const cell& cFaces = mesh_.cells()[cellI];
+
+    Map<label> pointToObjVert;
+    label objVertI = 0;
+
+    forAll(cFaces, i)
+    {
+        const face& f = mesh_.faces()[cFaces[i]];
+
+        forAll(f, fp)
+        {
+            if (pointToObjVert.insert(f[fp], objVertI))
+            {
+                meshTools::writeOBJ(str, mesh_.points()[f[fp]]);
+                objVertI++;
+            }
+        }
+    }
+
+    forAll(cFaces, i)
+    {
+        const face& f = mesh_.faces()[cFaces[i]];
+
+        forAll(f, fp)
+        {
+            label pointI = f[fp];
+            label nexPointI = f[f.fcIndex(fp)];
+
+            str << "l " << pointToObjVert[pointI]+1
+                << ' ' << pointToObjVert[nexPointI]+1 << nl;
+        }
+    }
+}
+
+
 // Find point with certain pointLevel. Skip any higher levels.
 Foam::label Foam::hexRef8::findLevel
 (
+    const label faceI,
     const face& f,
     const label startFp,
     const bool searchForward,
@@ -707,7 +738,13 @@ Foam::label Foam::hexRef8::findLevel
 
         if (pointLevel_[pointI] < wantedLevel)
         {
-            FatalErrorIn("hexRef8::findLevel")
+            dumpCell(mesh_.faceOwner()[faceI]);
+            if (mesh_.isInternalFace(faceI))
+            {
+                dumpCell(mesh_.faceNeighbour()[faceI]);
+            }
+
+            FatalErrorIn("hexRef8::findLevel(..)")
                 << "face:" << f
                 << " level:" << IndirectList<label>(pointLevel_, f)()
                 << " startFp:" << startFp
@@ -729,7 +766,13 @@ Foam::label Foam::hexRef8::findLevel
         }
     }
 
-    FatalErrorIn("hexRef8::findLevel")
+    dumpCell(mesh_.faceOwner()[faceI]);
+    if (mesh_.isInternalFace(faceI))
+    {
+        dumpCell(mesh_.faceNeighbour()[faceI]);
+    }
+
+    FatalErrorIn("hexRef8::findLevel(..)")
         << "face:" << f
         << " level:" << IndirectList<label>(pointLevel_, f)()
         << " startFp:" << startFp
@@ -763,15 +806,6 @@ Foam::label Foam::hexRef8::getAnchorLevel(const label faceI) const
         }
         else
         {
-            //FatalErrorIn("hexRef8::getAnchorLevel(const label) const")
-            //    << "face:" << faceI
-            //    << " centre:" << mesh_.faceCentres()[faceI]
-            //    << " verts:" << f
-            //    << " point levels:" << IndirectList<label>(pointLevel_, f)()
-            //    << " own:" << mesh_.faceOwner()[faceI]
-            //    << " ownLevel:" << cellLevel_[mesh_.faceOwner()[faceI]]
-            //    << abort(FatalError);
-
             return -1;
         }
     }
@@ -797,7 +831,7 @@ void Foam::hexRef8::checkInternalOrientation
 
     if ((dir & n) < 0)
     {
-        FatalErrorIn("checkInternalOrientation")
+        FatalErrorIn("checkInternalOrientation(..)")
             << "cell:" << cellI << " old face:" << faceI
             << " newFace:" << newFace << endl
             << " coords:" << compactPoints
@@ -812,7 +846,7 @@ void Foam::hexRef8::checkInternalOrientation
 
     if (s < 0.1 || s > 0.9)
     {
-        FatalErrorIn("checkInternalOrientation")
+        FatalErrorIn("checkInternalOrientation(..)")
             << "cell:" << cellI << " old face:" << faceI
             << " newFace:" << newFace << endl
             << " coords:" << compactPoints
@@ -843,7 +877,7 @@ void Foam::hexRef8::checkBoundaryOrientation
 
     if ((dir & n) < 0)
     {
-        FatalErrorIn("checkBoundaryOrientation")
+        FatalErrorIn("checkBoundaryOrientation(..)")
             << "cell:" << cellI << " old face:" << faceI
             << " newFace:" << newFace
             << " coords:" << compactPoints
@@ -858,7 +892,7 @@ void Foam::hexRef8::checkBoundaryOrientation
 
     if (s < 0.7 || s > 1.3)
     {
-        WarningIn("checkBoundaryOrientation")
+        WarningIn("checkBoundaryOrientation(..)")
             << "cell:" << cellI << " old face:" << faceI
             << " newFace:" << newFace
             << " coords:" << compactPoints
@@ -1191,8 +1225,22 @@ void Foam::hexRef8::createInternalFaces
             }
 
             // Now the face mid point is the second cLevel+1 point
-            label edgeMid = findLevel(f, f.fcIndex(anchorFp), true, cLevel+1);
-            label faceMid = findLevel(f, f.fcIndex(edgeMid), true, cLevel+1);
+            label edgeMid = findLevel
+            (
+                faceI,
+                f,
+                f.fcIndex(anchorFp),
+                true,
+                cLevel+1
+            );
+            label faceMid = findLevel
+            (
+                faceI,
+                f,
+                f.fcIndex(edgeMid),
+                true,
+                cLevel+1
+            );
 
             faceMidPointI = f[faceMid];
         }
@@ -1205,7 +1253,13 @@ void Foam::hexRef8::createInternalFaces
         }
         else
         {
-            FatalErrorIn("createInternalFaces")
+            dumpCell(mesh_.faceOwner()[faceI]);
+            if (mesh_.isInternalFace(faceI))
+            {
+                dumpCell(mesh_.faceNeighbour()[faceI]);
+            }
+
+            FatalErrorIn("createInternalFaces(..)")
                 << "nAnchors:" << nAnchors
                 << " faceI:" << faceI
                 << abort(FatalError);
@@ -1243,9 +1297,11 @@ void Foam::hexRef8::createInternalFaces
 
                     if (edgeMidPointI == -1)
                     {
+                        dumpCell(cellI);
+
                         const labelList cPoints(cellPoints(cellI));
 
-                        FatalErrorIn("createInternalFaces")
+                        FatalErrorIn("createInternalFaces(..)")
                             << "cell:" << cellI << " cLevel:" << cLevel
                             << " cell points:" << cPoints
                             << " pointLevel:"
@@ -1264,7 +1320,7 @@ void Foam::hexRef8::createInternalFaces
                 else
                 {
                     // Search forward in face to clevel+1
-                    label edgeMid = findLevel(f, fp1, true, cLevel+1);
+                    label edgeMid = findLevel(faceI, f, fp1, true, cLevel+1);
 
                     edgeMidPointI = f[edgeMid];
                 }
@@ -1314,9 +1370,11 @@ void Foam::hexRef8::createInternalFaces
 
                     if (edgeMidPointI == -1)
                     {
+                        dumpCell(cellI);
+
                         const labelList cPoints(cellPoints(cellI));
 
-                        FatalErrorIn("createInternalFaces")
+                        FatalErrorIn("createInternalFaces(..)")
                             << "cell:" << cellI << " cLevel:" << cLevel
                             << " cell points:" << cPoints
                             << " pointLevel:"
@@ -1335,7 +1393,14 @@ void Foam::hexRef8::createInternalFaces
                 else
                 {
                     // Search back to clevel+1
-                    label edgeMid = findLevel(f, fpMin1, false, cLevel+1);
+                    label edgeMid = findLevel
+                    (
+                        faceI,
+                        f,
+                        fpMin1,
+                        false,
+                        cLevel+1
+                    );
 
                     edgeMidPointI = f[edgeMid];
                 }
@@ -1591,9 +1656,11 @@ void Foam::hexRef8::checkWantedRefinementLevels
 
         if (mag(ownLevel-neiLevel) > 1)
         {
+            dumpCell(own);
+            dumpCell(nei);
             FatalErrorIn
             (
-                "hexRef8::checkRefinementLevels(const labelList&)"
+                "hexRef8::checkWantedRefinementLevels(const labelList&)"
             )   << "cell:" << own
                 << " current level:" << cellLevel_[own]
                 << " level after refinement:" << ownLevel
@@ -1633,9 +1700,10 @@ void Foam::hexRef8::checkWantedRefinementLevels
         {
             label patchI = mesh_.boundaryMesh().whichPatch(faceI);
 
+            dumpCell(own);
             FatalErrorIn
             (
-                "hexRef8::checkRefinementLevels(const labelList&)"
+                "hexRef8::checkWantedRefinementLevels(const labelList&)"
             )   << "Celllevel does not satisfy 2:1 constraint."
                 << " On coupled face "
                 << faceI
@@ -1730,6 +1798,25 @@ Foam::hexRef8::hexRef8(const polyMesh& mesh)
             << abort(FatalError);
     }
 
+    if
+    (
+        cellLevel_.size() != mesh_.nCells()
+     || pointLevel_.size() != mesh_.nPoints()
+    )
+    {
+        FatalErrorIn
+        (
+            "hexRef8::hexRef8(const polyMesh&)"
+        )   << "Restarted from inconsistent cellLevel or pointLevel files."
+            << endl
+            << "Number of cells in mesh:" << mesh_.nCells()
+            << " does not equal size of cellLevel:" << cellLevel_.size() << endl
+            << "Number of points in mesh:" << mesh_.nPoints()
+            << " does not equal size of pointLevel:" << pointLevel_.size()
+            << abort(FatalError);
+    }
+
+
     // Check refinement levels for consistency
     checkRefinementLevels(-1, labelList(0));
 
@@ -1806,7 +1893,24 @@ Foam::hexRef8::hexRef8
         )   << "History enabled but number of visible cells in it "
             << history_.visibleCells().size()
             << " is not equal to the number of cells in the mesh "
-            << mesh_.nCells()
+            << mesh_.nCells() << abort(FatalError);
+    }
+
+    if
+    (
+        cellLevel_.size() != mesh_.nCells()
+     || pointLevel_.size() != mesh_.nPoints()
+    )
+    {
+        FatalErrorIn
+        (
+            "hexRef8::hexRef8(const polyMesh&, const labelList&"
+            ", const labelList&, const refinementHistory&)"
+        )   << "Incorrect cellLevel or pointLevel size." << endl
+            << "Number of cells in mesh:" << mesh_.nCells()
+            << " does not equal size of cellLevel:" << cellLevel_.size() << endl
+            << "Number of points in mesh:" << mesh_.nPoints()
+            << " does not equal size of pointLevel:" << pointLevel_.size()
             << abort(FatalError);
     }
 
@@ -1877,6 +1981,24 @@ Foam::hexRef8::hexRef8
     savedPointLevel_(0),
     savedCellLevel_(0)
 {
+    if
+    (
+        cellLevel_.size() != mesh_.nCells()
+     || pointLevel_.size() != mesh_.nPoints()
+    )
+    {
+        FatalErrorIn
+        (
+            "hexRef8::hexRef8(const polyMesh&, const labelList&"
+            ", const labelList&)"
+        )   << "Incorrect cellLevel or pointLevel size." << endl
+            << "Number of cells in mesh:" << mesh_.nCells()
+            << " does not equal size of cellLevel:" << cellLevel_.size() << endl
+            << "Number of points in mesh:" << mesh_.nPoints()
+            << " does not equal size of pointLevel:" << pointLevel_.size()
+            << abort(FatalError);
+    }
+
     // Check refinement levels for consistency
     checkRefinementLevels(-1, labelList(0));
 
@@ -2356,9 +2478,13 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement
 
             if (mag(ownLevel-neiLevel) > 1)
             {
+                dumpCell(own);
+                dumpCell(nei);
                 FatalErrorIn
                 (
                     "hexRef8::consistentSlowRefinement"
+                    "(const label, const labelList&, const labelList&"
+                    ", const label, const labelList&)"
                 )   << "cell:" << own
                     << " current level:" << cellLevel_[own]
                     << " current refData:" << allCellInfo[own]
@@ -2412,11 +2538,14 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement
 
             if (mag(ownLevel - nbrLevel) > 1)
             {
+                dumpCell(own);
                 label patchI = mesh_.boundaryMesh().whichPatch(faceI);
 
                 FatalErrorIn
                 (
-                    "hexRef8::checkRefinementLevels(const labelList&)"
+                    "hexRef8::consistentSlowRefinement"
+                    "(const label, const labelList&, const labelList&"
+                    ", const label, const labelList&)"
                 )   << "Celllevel does not satisfy 2:1 constraint."
                     << " On coupled face "
                     << faceI
@@ -2805,7 +2934,6 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
             << cellsToRefine.size() << " to " << newCellsToRefine.size()
             << " cells to refine." << endl;
 
-//XXXX
         // Check that newCellsToRefine obeys at least 2:1.
 
         {
@@ -2861,6 +2989,7 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
                  && savedRefineCell.get(cellI) == 0
                 )
                 {
+                    dumpCell(cellI);
                     FatalErrorIn
                     (
                         "hexRef8::consistentSlowRefinement2"
@@ -2873,7 +3002,6 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
                 }
             }
         }
-//XXXX
     }
 
     return newCellsToRefine;
@@ -3342,6 +3470,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
                 {
                     if (nAnchorPoints[cellI] == 8)
                     {
+                        dumpCell(cellI);
                         FatalErrorIn
                         (
                             "hexRef8::setRefinement(const labelList&"
@@ -3365,6 +3494,8 @@ Foam::labelListList Foam::hexRef8::setRefinement
             {
                 if (nAnchorPoints[cellI] != 8)
                 {
+                    dumpCell(cellI);
+
                     const labelList cPoints(cellPoints(cellI));
 
                     FatalErrorIn
@@ -3841,7 +3972,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
 
         if (minPointI != labelMax && minPointI != mesh_.nPoints())
         {
-            FatalErrorIn("hexRef8::setRefinement")
+            FatalErrorIn("hexRef8::setRefinement(..)")
                 << "Added point labels not consecutive to existing mesh points."
                 << nl
                 << "mesh_.nPoints():" << mesh_.nPoints()
@@ -3995,13 +4126,6 @@ void Foam::hexRef8::updateMesh
 
                 if (oldCellI == -1)
                 {
-                    //FatalErrorIn("hexRef8::updateMesh(const mapPolyMesh&)")
-                    //    << "Problem : cell " << newCellI
-                    //    << " at " << mesh_.cellCentres()[newCellI]
-                    //    << " does not originate from another cell"
-                    //    << " (i.e. is inflated)." << nl
-                    //    << "Hence we cannot determine the new cellLevel"
-                    //    << " for it." << abort(FatalError);
                     newCellLevel[newCellI] = -1;
                 }
                 else
@@ -4036,8 +4160,8 @@ void Foam::hexRef8::updateMesh
         //{
         //    WarningIn("hexRef8::updateMesh(const mapPolyMesh&)")
         //        << "Problem : "
-        //        << "cellLevel_ contains illegal value -1 after mapping at cell "
-        //        << findIndex(cellLevel_, -1) << endl
+        //        << "cellLevel_ contains illegal value -1 after mapping
+        //        << " at cell " << findIndex(cellLevel_, -1) << endl
         //        << "This means that another program has inflated cells"
         //        << " (created cells out-of-nothing) and hence we don't know"
         //        << " their cell level. Continuing with illegal value."
@@ -4174,7 +4298,7 @@ void Foam::hexRef8::subset
 
         if (findIndex(cellLevel_, -1) != -1)
         {
-            FatalErrorIn("hexRef8::subset")
+            FatalErrorIn("hexRef8::subset(..)")
                 << "Problem : "
                 << "cellLevel_ contains illegal value -1 after mapping:"
                 << cellLevel_
@@ -4195,7 +4319,7 @@ void Foam::hexRef8::subset
 
         if (findIndex(pointLevel_, -1) != -1)
         {
-            FatalErrorIn("hexRef8::subset")
+            FatalErrorIn("hexRef8::subset(..)")
                 << "Problem : "
                 << "pointLevel_ contains illegal value -1 after mapping:"
                 << pointLevel_
@@ -4292,6 +4416,7 @@ void Foam::hexRef8::checkMesh() const
 
                     if (!cellToFace.insert(labelPair(own, nei[bFaceI]), faceI))
                     {
+                        dumpCell(own);
                         FatalErrorIn("hexRef8::checkMesh()")
                             << "Faces do not seem to be correct across coupled"
                             << " boundaries" << endl
@@ -4335,6 +4460,8 @@ void Foam::hexRef8::checkMesh() const
                 const face& f = mesh_.faces()[faceI];
                 label patchI = mesh_.boundaryMesh().whichPatch(faceI);
 
+                dumpCell(mesh_.faceOwner()[faceI]);
+
                 FatalErrorIn("hexRef8::checkMesh()")
                     << "Faces do not seem to be correct across coupled"
                     << " boundaries" << endl
@@ -4373,6 +4500,8 @@ void Foam::hexRef8::checkMesh() const
 
             if (f.size() != nVerts[i])
             {
+                dumpCell(mesh_.faceOwner()[faceI]);
+
                 label patchI = mesh_.boundaryMesh().whichPatch(faceI);
 
                 FatalErrorIn("hexRef8::checkMesh()")
@@ -4421,6 +4550,8 @@ void Foam::hexRef8::checkMesh() const
 
             if (mag(anchorVec - anchorPoints[i]) > smallDim)
             {
+                dumpCell(mesh_.faceOwner()[faceI]);
+
                 label patchI = mesh_.boundaryMesh().whichPatch(faceI);
 
                 FatalErrorIn("hexRef8::checkMesh()")
@@ -4487,6 +4618,9 @@ void Foam::hexRef8::checkRefinementLevels
 
             if (mag(cellLevel_[own] - cellLevel_[nei]) > 1)
             {
+                dumpCell(own);
+                dumpCell(nei);
+
                 FatalErrorIn
                 (
                     "hexRef8::checkRefinementLevels(const label)"
@@ -4520,6 +4654,8 @@ void Foam::hexRef8::checkRefinementLevels
 
             if (mag(cellLevel_[own] - neiLevel[i]) > 1)
             {
+                dumpCell(own);
+
                 label patchI = mesh_.boundaryMesh().whichPatch(faceI);
 
                 FatalErrorIn
@@ -4621,6 +4757,8 @@ void Foam::hexRef8::checkRefinementLevels
                   > maxPointDiff
                 )
                 {
+                    dumpCell(cellI);
+
                     FatalErrorIn
                     (
                         "hexRef8::checkRefinementLevels(const label)"
@@ -4637,6 +4775,72 @@ void Foam::hexRef8::checkRefinementLevels
             }
         }
     }
+
+
+    // Hanging points
+    {
+        // Any patches with points having only two edges.
+
+        boolList isHangingPoint(mesh_.nPoints(), false);
+
+        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+        forAll(patches, patchI)
+        {
+            const polyPatch& pp = patches[patchI];
+
+            const labelList& meshPoints = pp.meshPoints();
+
+            forAll(meshPoints, i)
+            {
+                label pointI = meshPoints[i];
+
+                const labelList& pEdges = mesh_.pointEdges()[pointI];
+
+                if (pEdges.size() == 2)
+                {
+                    isHangingPoint[pointI] = true;
+                }
+            }
+        }
+
+        syncTools::syncPointList
+        (
+            mesh_,
+            isHangingPoint,
+            andEqOp<bool>(),        // only if all decide it is hanging point
+            true,                   // null
+            false                   // no separation
+        );
+
+        OFstream str(mesh_.time().path()/"hangingPoints.obj");
+
+        label nHanging = 0;
+
+        forAll(isHangingPoint, pointI)
+        {
+            if (isHangingPoint[pointI])
+            {
+                nHanging++;
+
+                Pout<< "Hanging boundary point " << pointI
+                    << " at " << mesh_.points()[pointI]
+                    << endl;
+
+                meshTools::writeOBJ(str, mesh_.points()[pointI]);
+            }
+        }
+
+        if (returnReduce(nHanging, sumOp<label>()) > 0)
+        {
+            FatalErrorIn
+            (
+                "hexRef8::checkRefinementLevels(const label)"
+            )   << "Detected a point used by two edges only (hanging point)"
+                << nl << "This is not allowed"
+                << abort(FatalError);
+        }
+    }
 }
 
 
@@ -4865,8 +5069,10 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
 
     if (maxSet)
     {
-        FatalErrorIn("hexRef8::consistentUnrefinement")
-            << "maxSet not implemented yet."
+        FatalErrorIn
+        (
+            "hexRef8::consistentUnrefinement(const labelList&, const bool"
+        )   << "maxSet not implemented yet."
             << abort(FatalError);
     }
 
@@ -4935,7 +5141,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
                 {
                     if (unrefineCell.get(own) == 0)
                     {
-                        FatalErrorIn("hexRef8::consistentUnrefinement")
+                        FatalErrorIn("hexRef8::consistentUnrefinement(..)")
                             << "problem" << abort(FatalError);
                     }
 
@@ -4953,7 +5159,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
                 {
                     if (unrefineCell.get(nei) == 0)
                     {
-                        FatalErrorIn("hexRef8::consistentUnrefinement")
+                        FatalErrorIn("hexRef8::consistentUnrefinement(..)")
                             << "problem" << abort(FatalError);
                     }
 
@@ -4989,7 +5195,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
                 {
                     if (unrefineCell.get(own) == 0)
                     {
-                        FatalErrorIn("hexRef8::consistentUnrefinement")
+                        FatalErrorIn("hexRef8::consistentUnrefinement(..)")
                             << "problem" << abort(FatalError);
                     }
 
@@ -5003,7 +5209,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
                 {
                     if (unrefineCell.get(own) == 1)
                     {
-                        FatalErrorIn("hexRef8::consistentUnrefinement")
+                        FatalErrorIn("hexRef8::consistentUnrefinement(..)")
                             << "problem" << abort(FatalError);
                     }
 
@@ -5086,8 +5292,10 @@ void Foam::hexRef8::setUnrefinement
 {
     if (!history_.active())
     {
-        FatalErrorIn("hexRef8::setUnrefinement()")
-            << "Only call if constructed with history capability"
+        FatalErrorIn
+        (
+            "hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
+        )   << "Only call if constructed with history capability"
             << abort(FatalError);
     }
 
@@ -5102,8 +5310,10 @@ void Foam::hexRef8::setUnrefinement
         {
             if (cellLevel_[cellI] < 0)
             {
-                FatalErrorIn("hexRef8::setUnrefinement()")
-                    << "Illegal cell level " << cellLevel_[cellI]
+                FatalErrorIn
+                (
+                    "hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
+                )   << "Illegal cell level " << cellLevel_[cellI]
                     << " for cell " << cellI
                     << abort(FatalError);
             }
@@ -5165,8 +5375,10 @@ void Foam::hexRef8::setUnrefinement
 
         if (facesToRemove.size() != splitFaces.size())
         {
-            FatalErrorIn("hexRef8::setUnrefinement")
-                << "Ininitial set of split points to unrefine does not"
+            FatalErrorIn
+            (
+                "hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
+            )   << "Ininitial set of split points to unrefine does not"
                 << " seem to be consistent or not mid points of refined cells"
                 << abort(FatalError);
         }
@@ -5187,8 +5399,10 @@ void Foam::hexRef8::setUnrefinement
         // Check
         if (pCells.size() != 8)
         {
-            FatalErrorIn("hexRef8::setUnrefinement")
-                << "splitPoint " << pointI
+            FatalErrorIn
+            (
+                "hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
+            )   << "splitPoint " << pointI
                 << " should have 8 cells using it. It has " << pCells
                 << abort(FatalError);
         }
@@ -5208,7 +5422,7 @@ void Foam::hexRef8::setUnrefinement
 
                 if (region == -1)
                 {
-                    FatalErrorIn("hexRef8::setUnrefinement")
+                    FatalErrorIn("hexRef8::setUnrefinement(..)")
                         << "Ininitial set of split points to unrefine does not"
                         << " seem to be consistent or not mid points"
                         << " of refined cells" << nl
@@ -5219,7 +5433,7 @@ void Foam::hexRef8::setUnrefinement
 
                 if (masterCellI != cellRegionMaster[region])
                 {
-                    FatalErrorIn("hexRef8::setUnrefinement")
+                    FatalErrorIn("hexRef8::setUnrefinement(..)")
                         << "cell:" << cellI << " on splitPoint:" << pointI
                         << " in region " << region
                         << " has master:" << cellRegionMaster[region]
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H
index 9d9ca90414..faccfdc271 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H
@@ -174,9 +174,12 @@ class hexRef8
         label findMaxLevel(const labelList& f) const;
         //- Count number of vertices <= anchorLevel
         label countAnchors(const labelList&, const label) const;
+        //- Debugging: dump cell as .obj file
+        void dumpCell(const label cellI) const;
         //- Find index of point with wantedLevel, starting from fp.
         label findLevel
         (
+            const label faceI,
             const face& f,
             const label startFp,
             const bool searchForward,
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.C
index 28fe2a21ee..fab2b74609 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.C
@@ -807,7 +807,7 @@ void Foam::removeFaces::setRefinement
     // Edges to remove
     labelHashSet edgesToRemove(faceLabels.size());
 
-    // Per face the region it is. -1 for removed faces, -2 for regions
+    // Per face the region it is in. -1 for removed faces, -2 for regions
     // consisting of single face only.
     labelList faceRegion(mesh_.nFaces(), -1);
 
@@ -1258,10 +1258,15 @@ void Foam::removeFaces::setRefinement
     // are only used by 2 unremoved edges.
     {
         // Usage of points by non-removed edges.
-        labelList nEdgesPerPoint(mesh_.nPoints(), labelMax);
+        labelList nEdgesPerPoint(mesh_.nPoints());
 
         const labelListList& pointEdges = mesh_.pointEdges();
 
+        forAll(pointEdges, pointI)
+        {
+            nEdgesPerPoint[pointI] = pointEdges[pointI].size();
+        }
+
         forAllConstIter(labelHashSet, edgesToRemove, iter)
         {
             // Edge will get removed.
@@ -1269,16 +1274,7 @@ void Foam::removeFaces::setRefinement
 
             forAll(e, i)
             {
-                label pointI = e[i];
-
-                if (nEdgesPerPoint[pointI] == labelMax)
-                {
-                    nEdgesPerPoint[pointI] = pointEdges[pointI].size()-1;
-                }
-                else
-                {
-                    nEdgesPerPoint[pointI]--;
-                }
+                nEdgesPerPoint[e[i]]--;
             }
         }
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C
index 3c4768855f..e9a32b05d0 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C
@@ -159,7 +159,7 @@ void directMappedFixedValueFvPatchField<Type>::updateCoeffs()
 
     // Get the scheduling information
     const List<labelPair>& schedule = mpp.schedule();
-    const labelListList& sendCellLabels = mpp.sendCellLabels();
+    const labelListList& sendLabels = mpp.sendLabels();
     const labelListList& receiveFaceLabels = mpp.receiveFaceLabels();
 
 
@@ -177,7 +177,7 @@ void directMappedFixedValueFvPatchField<Type>::updateCoeffs()
             toProc<< IndirectList<Type>
             (
                 this->internalField(),
-                sendCellLabels[recvProc]
+                sendLabels[recvProc]
             )();
         }
         else
@@ -204,7 +204,7 @@ void directMappedFixedValueFvPatchField<Type>::updateCoeffs()
         IndirectList<Type> fromFld
         (
             this->internalField(),
-            sendCellLabels[Pstream::myProcNo()]
+            sendLabels[Pstream::myProcNo()]
         );
     
         // Destination faces
diff --git a/src/lagrangian/basic/Cloud/CloudIO.C b/src/lagrangian/basic/Cloud/CloudIO.C
index ea37ea27a9..0401955801 100644
--- a/src/lagrangian/basic/Cloud/CloudIO.C
+++ b/src/lagrangian/basic/Cloud/CloudIO.C
@@ -118,7 +118,8 @@ void Foam::Cloud<ParticleType>::checkFieldIOobject
             "void Cloud<ParticleType>::checkFieldIOobject"
             "(Cloud<ParticleType>, IOField<DataType>)"
         )   << "Size of " << data.name()
-            << " field does not match the number of particles"
+            << " field " << data.size()
+            << " does not match the number of particles " << c.size()
             << abort(FatalError);
     }
 }
diff --git a/src/meshTools/Make/options b/src/meshTools/Make/options
index 3e8506b54f..ef1033b0e8 100644
--- a/src/meshTools/Make/options
+++ b/src/meshTools/Make/options
@@ -3,5 +3,5 @@ EXE_INC = \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude
 
 LIB_LIBS = \
-   -ltriSurface \
-   -llagrangian
+    -ltriSurface \
+    -llagrangian
diff --git a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C
index e519da59d9..1c3ab42ffb 100644
--- a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C
+++ b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C
@@ -40,6 +40,18 @@ namespace Foam
 
     addToRunTimeSelectionTable(polyPatch, directMappedPolyPatch, word);
     addToRunTimeSelectionTable(polyPatch, directMappedPolyPatch, dictionary);
+
+
+    template<>
+    const char* NamedEnum<directMappedPolyPatch::sampleMode, 3>::names[] =
+    {
+        "nearestCell",
+        "nearestPatchFace",
+        "nearestFace"
+    };
+
+    const NamedEnum<directMappedPolyPatch::sampleMode, 3>
+        directMappedPolyPatch::sampleModeNames_;
 }
 
 
@@ -113,107 +125,236 @@ void Foam::directMappedPolyPatch::collectSamples
 void Foam::directMappedPolyPatch::findSamples
 (
     const pointField& samples,
-    labelList& sampleCellProcs,
-    labelList& sampleCells,
-    pointField& sampleCc
+    labelList& sampleProcs,
+    labelList& sampleIndices,
+    pointField& sampleLocations
 ) const
 {
-    sampleCellProcs.setSize(samples.size());
-    sampleCells.setSize(samples.size());
-    sampleCc.setSize(samples.size());
-    sampleCc = point(-GREAT, -GREAT, -GREAT);
+    const polyMesh& mesh = boundaryMesh().mesh();
 
+    // All the info for nearest. Construct to miss
+    List<nearInfo> nearest(samples.size());
+
+    switch (mode_)
     {
-        // Octree based search engine
-        meshSearch meshSearchEngine(boundaryMesh().mesh(), false);
-
-        forAll(samples, sampleI)
+        case NEARESTCELL:
         {
-            sampleCells[sampleI] = meshSearchEngine.findCell(samples[sampleI]);
+            if (samplePatch_.size() != 0 && samplePatch_ != "none")
+            {
+                FatalErrorIn
+                (
+                    "directMappedPolyPatch::findSamples(const pointField&,"
+                    " labelList&, labelList&, pointField&) const"
+                )   << "No need to supply a patch name when in "
+                    << sampleModeNames_[mode_] << " mode." << exit(FatalError);
+            }
 
-            if (sampleCells[sampleI] == -1)
+            // Octree based search engine
+            meshSearch meshSearchEngine(mesh, false);
+
+            forAll(samples, sampleI)
             {
-                sampleCellProcs[sampleI] = -1;
+                const point& sample = samples[sampleI];
+
+                label cellI = meshSearchEngine.findCell(sample);
+
+                if (cellI == -1)
+                {
+                    nearest[sampleI].second().first() = Foam::sqr(GREAT);
+                    nearest[sampleI].second().second() = Pstream::myProcNo();
+                }
+                else
+                {
+                    const point& cc = mesh.cellCentres()[cellI];
+
+                    nearest[sampleI].first() = pointIndexHit
+                    (
+                        true,
+                        cc,
+                        cellI
+                    );
+                    nearest[sampleI].second().first() = magSqr(cc-sample);
+                    nearest[sampleI].second().second() = Pstream::myProcNo();
+                }
             }
-            else
+            break;
+        }
+
+        case NEARESTPATCHFACE:
+        {
+            label patchI = boundaryMesh().findPatchID(samplePatch_);
+
+            if (patchI == -1)
             {
-                sampleCellProcs[sampleI] = Pstream::myProcNo();
-                sampleCc[sampleI] =
-                    boundaryMesh().mesh().cellCentres()[sampleCells[sampleI]];
+                FatalErrorIn("directMappedPolyPatch::findSamples(..)")
+                    << "Cannot find patch " << samplePatch_ << endl
+                    << "Valid patches are " << boundaryMesh().names()
+                    << exit(FatalError);
             }
+
+            Random rndGen(123456);
+
+            const polyPatch& pp = boundaryMesh()[patchI];
+
+            // patch faces
+            const labelList patchFaces(identity(pp.size()) + pp.start());
+
+            const treeBoundBox patchBb
+            (
+                treeBoundBox(pp.points(), pp.meshPoints()).extend
+                (
+                    rndGen,
+                    1E-4
+                )
+            );
+
+            autoPtr<indexedOctree<treeDataFace> > boundaryTree
+            (
+                new indexedOctree<treeDataFace>
+                (
+                    treeDataFace    // all information needed to search faces
+                    (
+                        false,                      // do not cache bb
+                        mesh,
+                        patchFaces                  // boundary faces only
+                    ),
+                    patchBb,   // overall search domain
+                    8,                              // maxLevel
+                    10,                             // leafsize
+                    3.0                             // duplicity
+                )
+            );
+
+            forAll(samples, sampleI)
+            {
+                const point& sample = samples[sampleI];
+
+                pointIndexHit& nearInfo = nearest[sampleI].first();
+                nearInfo = boundaryTree().findNearest
+                (
+                    sample,
+                    magSqr(patchBb.max()-patchBb.min())
+                );
+
+                if (!nearInfo.hit())
+                {
+                    nearest[sampleI].second().first() = Foam::sqr(GREAT);
+                    nearest[sampleI].second().second() = Pstream::myProcNo();
+                }
+                else
+                {
+                    point fc(pp[nearInfo.index()].centre(pp.points()));
+                    nearInfo.setPoint(fc);
+                    nearest[sampleI].second().first() = magSqr(fc-sample);
+                    nearest[sampleI].second().second() = Pstream::myProcNo();
+                }
+            }
+            break;
+        }
+
+        case NEARESTFACE:
+        {
+            if (samplePatch_.size() != 0 && samplePatch_ != "none")
+            {
+                FatalErrorIn
+                (
+                    "directMappedPolyPatch::findSamples(const pointField&,"
+                    " labelList&, labelList&, pointField&) const"
+                )   << "No need to supply a patch name when in "
+                    << sampleModeNames_[mode_] << " mode." << exit(FatalError);
+            }
+
+            // Octree based search engine
+            meshSearch meshSearchEngine(mesh, false);
+
+            forAll(samples, sampleI)
+            {
+                const point& sample = samples[sampleI];
+
+                label faceI = meshSearchEngine.findNearestFace(sample);
+
+                if (faceI == -1)
+                {
+                    nearest[sampleI].second().first() = Foam::sqr(GREAT);
+                    nearest[sampleI].second().second() = Pstream::myProcNo();
+                }
+                else
+                {
+                    const point& fc = mesh.faceCentres()[faceI];
+
+                    nearest[sampleI].first() = pointIndexHit
+                    (
+                        true,
+                        fc,
+                        faceI
+                    );
+                    nearest[sampleI].second().first() = magSqr(fc-sample);
+                    nearest[sampleI].second().second() = Pstream::myProcNo();
+                }
+            }
+            break;
+        }
+
+        default:
+        {
+            FatalErrorIn("directMappedPolyPatch::findSamples(..)")
+                << "problem." << abort(FatalError);
         }
     }
 
-    // Use only that processor that contains the sample
-    Pstream::listCombineGather(sampleCells, maxEqOp<label>());
-    Pstream::listCombineScatter(sampleCells);
 
-    labelList minSampleCellProcs(sampleCellProcs);
-    Pstream::listCombineGather(sampleCellProcs, maxEqOp<label>());
-    Pstream::listCombineScatter(sampleCellProcs);
-
-    Pstream::listCombineGather(sampleCc, maxEqOp<point>());
-    Pstream::listCombineScatter(sampleCc);
+    // Find nearest.
+    Pstream::listCombineGather(nearest, nearestEqOp());
+    Pstream::listCombineScatter(nearest);
 
     if (debug)
     {
         Info<< "directMappedPolyPatch::findSamples : " << endl;
-        forAll(sampleCells, sampleI)
+        forAll(nearest, sampleI)
         {
-            Info<< "    " << sampleI << " coord:" << samples[sampleI]
-                << " found on processor:" << sampleCellProcs[sampleI]
-                << " in cell:" << sampleCells[sampleI]
-                << " with cc:" << sampleCc[sampleI] << endl;
+            label procI = nearest[sampleI].second().second();
+            label localI = nearest[sampleI].first().index();
+
+            Info<< "    " << sampleI << " coord:"<< samples[sampleI]
+                << " found on processor:" << procI
+                << " in local cell/face:" << localI
+                << " with cc:" << nearest[sampleI].first().rawPoint() << endl;
         }
     }
 
     // Check for samples not being found
-    forAll(sampleCells, sampleI)
+    forAll(nearest, sampleI)
     {
-        if (sampleCells[sampleI] == -1)
+        if (!nearest[sampleI].first().hit())
         {
             FatalErrorIn
             (
-                "findSamples(const pointField&, labelList&, labelList&)"
+                "directMappedPolyPatch::findSamples"
+                "(const pointField&, labelList&"
+                ", labelList&, pointField&)"
             )   << "Did not find sample " << samples[sampleI]
                 << " on any processor." << exit(FatalError);
         }
     }
 
 
-    // Check for samples being found in multiple cells
-    {
-        forAll(minSampleCellProcs, sampleI)
-        {
-            if (minSampleCellProcs[sampleI] == -1)
-            {
-                minSampleCellProcs[sampleI] = labelMax;
-            }
-        }
-        Pstream::listCombineGather(minSampleCellProcs, minEqOp<label>());
-        Pstream::listCombineScatter(minSampleCellProcs);
+    // Convert back into proc+local index
+    sampleProcs.setSize(samples.size());
+    sampleIndices.setSize(samples.size());
+    sampleLocations.setSize(samples.size());
 
-        forAll(minSampleCellProcs, sampleI)
-        {
-            if (minSampleCellProcs[sampleI] != sampleCellProcs[sampleI])
-            {
-                FatalErrorIn
-                (
-                    "directMappedPolyPatch::findSamples"
-                    "(const pointField&, labelList&, labelList&)"
-                )   << "Found sample " << samples[sampleI]
-                    << " on processor " << minSampleCellProcs[sampleI]
-                    << " and on processor " << sampleCellProcs[sampleI] << endl
-                    << "This is illegal. {lease move your sampling plane."
-                    << exit(FatalError);
-            }
-        }
+    forAll(nearest, sampleI)
+    {
+        sampleProcs[sampleI] = nearest[sampleI].second().second();
+        sampleIndices[sampleI] = nearest[sampleI].first().index();
+        sampleLocations[sampleI] = nearest[sampleI].first().hitPoint();
     }
 }
 
+
 void Foam::directMappedPolyPatch::calcMapping() const
 {
-    if (sendCellLabelsPtr_.valid())
+    if (sendLabelsPtr_.valid())
     {
         FatalErrorIn("directMappedPolyPatch::calcMapping() const")
             << "Mapping already calculated" << exit(FatalError);
@@ -236,18 +377,18 @@ void Foam::directMappedPolyPatch::calcMapping() const
     pointField patchFc;
     collectSamples(samples, patchFaceProcs, patchFaces, patchFc);
 
-    // Find processor and cell samples are in
-    labelList sampleCellProcs;
-    labelList sampleCells;
-    pointField sampleCc;
-    findSamples(samples, sampleCellProcs, sampleCells, sampleCc);
+    // Find processor and cell/face samples are in and actual location.
+    labelList sampleProcs;
+    labelList sampleIndices;
+    pointField sampleLocations;
+    findSamples(samples, sampleProcs, sampleIndices, sampleLocations);
 
 
     // Now we have all the data we need:
     // - where sample originates from (so destination when mapping):
     //   patchFaces, patchFaceProcs.
-    // - cell sample is in (so source when mapping)
-    //   sampleCells, sampleCellProcs.
+    // - cell/face sample is in (so source when mapping)
+    //   sampleIndices, sampleProcs.
 
 
     if (debug && Pstream::master())
@@ -267,18 +408,18 @@ void Foam::directMappedPolyPatch::calcMapping() const
         {
             meshTools::writeOBJ(str, patchFc[i]);
             vertI++;
-            meshTools::writeOBJ(str, sampleCc[i]);
+            meshTools::writeOBJ(str, sampleLocations[i]);
             vertI++;
             str << "l " << vertI-1 << ' ' << vertI << nl;
         }
     }
 
 
-    // Check that actual offset vector (sampleCc - patchFc) is more or less
-    // constant.
+    // Check that actual offset vector (sampleLocations - patchFc) is more or
+    // less constant.
     if (Pstream::master())
     {
-        const scalarField magOffset(mag(sampleCc - patchFc));
+        const scalarField magOffset(mag(sampleLocations - patchFc));
         const scalar avgOffset(average(magOffset));
 
         forAll(magOffset, sampleI)
@@ -291,7 +432,7 @@ void Foam::directMappedPolyPatch::calcMapping() const
                     << "    on a single plane."
                     << " This might give numerical problems." << endl
                     << "    At patchface " << patchFc[sampleI]
-                    << " the sampled cell " << sampleCc[sampleI] << endl
+                    << " the sampled cell " << sampleLocations[sampleI] << endl
                     << "    is not on a plane " << avgOffset
                     << " offset from the patch." << endl
                     << "    You might want to shift your plane offset."
@@ -304,7 +445,7 @@ void Foam::directMappedPolyPatch::calcMapping() const
 
 
     // Determine schedule.
-    mapDistribute distMap(sampleCellProcs, patchFaceProcs);
+    mapDistribute distMap(sampleProcs, patchFaceProcs);
 
     // Rework the schedule to cell data to send, face data to receive.
     schedulePtr_.reset(new List<labelPair>(distMap.schedule()));
@@ -313,17 +454,21 @@ void Foam::directMappedPolyPatch::calcMapping() const
     const labelListList& constructMap = distMap.constructMap();
 
     // Extract the particular data I need to send and receive.
-    sendCellLabelsPtr_.reset(new labelListList(subMap.size()));
-    labelListList& sendCellLabels = sendCellLabelsPtr_();
+    sendLabelsPtr_.reset(new labelListList(subMap.size()));
+    labelListList& sendLabels = sendLabelsPtr_();
 
     forAll(subMap, procI)
     {
-        sendCellLabels[procI] = IndirectList<label>(sampleCells, subMap[procI]);
+        sendLabels[procI] = IndirectList<label>
+        (
+            sampleIndices,
+            subMap[procI]
+        );
 
         if (debug)
         {
-            Pout<< "To proc:" << procI << " sending values of cells:"
-                << sendCellLabels[procI] << endl;
+            Pout<< "To proc:" << procI << " sending values of cells/faces:"
+                << sendLabels[procI] << endl;
         }
     }
 
@@ -356,9 +501,11 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
 )
 :
     polyPatch(name, size, start, index, bm),
+    mode_(NEARESTPATCHFACE),
+    samplePatch_("none"),
     offset_(vector::zero),
     schedulePtr_(NULL),
-    sendCellLabelsPtr_(NULL),
+    sendLabelsPtr_(NULL),
     receiveFaceLabelsPtr_(NULL)
 {}
 
@@ -372,9 +519,11 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
 )
 :
     polyPatch(name, dict, index, bm),
+    mode_(sampleModeNames_.read(dict.lookup("sampleMode"))),
+    samplePatch_(dict.lookup("samplePatch")),
     offset_(dict.lookup("offset")),
     schedulePtr_(NULL),
-    sendCellLabelsPtr_(NULL),
+    sendLabelsPtr_(NULL),
     receiveFaceLabelsPtr_(NULL)
 {}
 
@@ -386,9 +535,11 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
 )
 :
     polyPatch(pp, bm),
+    mode_(pp.mode_),
+    samplePatch_(pp.samplePatch_),
     offset_(pp.offset_),
     schedulePtr_(NULL),
-    sendCellLabelsPtr_(NULL),
+    sendLabelsPtr_(NULL),
     receiveFaceLabelsPtr_(NULL)
 {}
 
@@ -403,9 +554,11 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
 )
 :
     polyPatch(pp, bm, index, newSize, newStart),
+    mode_(pp.mode_),
+    samplePatch_(pp.samplePatch_),
     offset_(pp.offset_),
     schedulePtr_(NULL),
-    sendCellLabelsPtr_(NULL),
+    sendLabelsPtr_(NULL),
     receiveFaceLabelsPtr_(NULL)
 {}
 
@@ -423,7 +576,7 @@ Foam::directMappedPolyPatch::~directMappedPolyPatch()
 void Foam::directMappedPolyPatch::clearOut()
 {
     schedulePtr_.clear();
-    sendCellLabelsPtr_.clear();
+    sendLabelsPtr_.clear();
     receiveFaceLabelsPtr_.clear();
 }
 
@@ -431,6 +584,10 @@ void Foam::directMappedPolyPatch::clearOut()
 void Foam::directMappedPolyPatch::write(Ostream& os) const
 {
     polyPatch::write(os);
+    os.writeKeyword("sampleMode") << sampleModeNames_[mode_]
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("samplePatch") << samplePatch_
+        << token::END_STATEMENT << nl;
     os.writeKeyword("offset") << offset_ << token::END_STATEMENT << nl;
 }
 
diff --git a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.H b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.H
index 06837c5f9b..144994d09a 100644
--- a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.H
+++ b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.H
@@ -26,8 +26,8 @@ Class
     Foam::directMappedPolyPatch
 
 Description
-    Determines a mapping between patch face centres and mesh cell centres and
-    processors they're on.
+    Determines a mapping between patch face centres and mesh cell or face
+    centres and processors they're on.
 
 Note
     Storage is not optimal. It stores all face centres and cells on all
@@ -43,6 +43,9 @@ SourceFiles
 
 #include "polyPatch.H"
 #include "labelPair.H"
+#include "Tuple2.H"
+#include "pointIndexHit.H"
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -57,19 +60,41 @@ class directMappedPolyPatch
 :
     public polyPatch
 {
+
+public:
+
+        //- Mesh items to sample
+        enum sampleMode
+        {
+            NEARESTCELL,
+            NEARESTPATCHFACE,
+            NEARESTFACE
+        };
+
+private:
+
     // Private data
 
+        static const NamedEnum<sampleMode, 3> sampleModeNames_;
+
+        //- What to sample
+        sampleMode mode_;
+
+        //- Patch (only if NEARESTBOUNDARY)
+        const word samplePatch_;
+
         //- Offset vector
         const vector offset_;
 
 
+
         // Derived information
 
             //- Communication schedule.
             mutable autoPtr<List<labelPair> > schedulePtr_;
 
-            //- Cells to sample per processor
-            mutable autoPtr<labelListList> sendCellLabelsPtr_;
+            //- Cells/faces to sample per processor
+            mutable autoPtr<labelListList> sendLabelsPtr_;
 
             //- Patch faces to receive per processor
             mutable autoPtr<labelListList> receiveFaceLabelsPtr_;
@@ -86,19 +111,50 @@ class directMappedPolyPatch
             pointField& patchFc
         ) const;
 
-        //- Find cells containing samples
+        //- Find cells/faces containing samples
         void findSamples
         (
             const pointField&,
-            labelList& sampleCellProcs,
-            labelList& sampleCells,
-            pointField& sampleCc
+            labelList& sampleProcs,     // processor containing sample
+            labelList& sampleIndices,   // local index of cell/face
+            pointField& sampleLocations // actual representative location
         ) const;
 
         //- Calculate matching
         void calcMapping() const;
 
 
+    // Private class
+
+        //- Private class for finding nearest
+        //  - point+local index
+        //  - sqr(distance)
+        //  - processor
+        typedef Tuple2<pointIndexHit, Tuple2<scalar, label> > nearInfo;
+
+        class nearestEqOp
+        {
+
+        public:
+
+            void operator()(nearInfo& x, const nearInfo& y) const
+            {
+                if (y.first().hit())
+                {
+                    if (!x.first().hit())
+                    {
+                        x = y;
+                    }
+                    else if (y.second().first() < x.second().first())
+                    {
+                        x = y;
+                    }
+                }
+            }
+        };
+
+
+
 protected:
 
         void clearOut();
@@ -231,21 +287,25 @@ public:
             return schedulePtr_();
         }
 
-        //- Cells to sample per processor
-        const labelListList& sendCellLabels() const
+        //- Cells/faces to sample per processor
+        const labelListList& sendLabels() const
         {
-            if (!sendCellLabelsPtr_.valid())
+            Pout<< "Asking for sendLabels." << endl;
+            if (!sendLabelsPtr_.valid())
             {
+                Pout<< "Calculating mapping." << endl;
                 calcMapping();
             }
-            return sendCellLabelsPtr_();
+            return sendLabelsPtr_();
         }
 
         //- Patch faces to receive per processor
         const labelListList& receiveFaceLabels() const
         {
+            Pout<< "Asking for receiveFaceLabels." << endl;
             if (!receiveFaceLabelsPtr_.valid())
             {
+                Pout<< "Calculating mapping." << endl;
                 calcMapping();
             }
             return receiveFaceLabelsPtr_();
diff --git a/src/meshTools/indexedOctree/indexedOctree.C b/src/meshTools/indexedOctree/indexedOctree.C
index 1a55a88547..44259482ab 100644
--- a/src/meshTools/indexedOctree/indexedOctree.C
+++ b/src/meshTools/indexedOctree/indexedOctree.C
@@ -26,7 +26,7 @@ License
 
 #include "indexedOctree.H"
 #include "linePointRef.H"
-#include "triSurface.H"
+//#include "triSurface.H"
 #include "meshTools.H"
 #include "OFstream.H"
 
diff --git a/src/meshTools/indexedOctree/treeDataCell.C b/src/meshTools/indexedOctree/treeDataCell.C
index cb41899306..131f0c69be 100644
--- a/src/meshTools/indexedOctree/treeDataCell.C
+++ b/src/meshTools/indexedOctree/treeDataCell.C
@@ -26,7 +26,7 @@ License
 
 #include "treeDataCell.H"
 #include "indexedOctree.H"
-#include "polyMesh.H"
+#include "primitiveMesh.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -71,7 +71,7 @@ Foam::treeBoundBox Foam::treeDataCell::calcCellBb(const label cellI) const
 Foam::treeDataCell::treeDataCell
 (
     const bool cacheBb,
-    const polyMesh& mesh,
+    const primitiveMesh& mesh,
     const labelList& cellLabels
 )
 :
@@ -94,7 +94,7 @@ Foam::treeDataCell::treeDataCell
 Foam::treeDataCell::treeDataCell
 (
     const bool cacheBb,
-    const polyMesh& mesh
+    const primitiveMesh& mesh
 )
 :
     mesh_(mesh),
diff --git a/src/meshTools/indexedOctree/treeDataCell.H b/src/meshTools/indexedOctree/treeDataCell.H
index 1cb15db6fe..6fe5d86755 100644
--- a/src/meshTools/indexedOctree/treeDataCell.H
+++ b/src/meshTools/indexedOctree/treeDataCell.H
@@ -47,7 +47,7 @@ namespace Foam
 {
 
 // Forward declaration of classes
-class polyMesh;
+class primitiveMesh;
 template<class Type> class indexedOctree;
 
 /*---------------------------------------------------------------------------*\
@@ -58,7 +58,7 @@ class treeDataCell
 {
     // Private data
 
-        const polyMesh& mesh_;
+        const primitiveMesh& mesh_;
 
         //- Subset of cells to work on
         const labelList cellLabels_;
@@ -87,12 +87,12 @@ public:
         treeDataCell
         (
             const bool cacheBb,
-            const polyMesh&,
+            const primitiveMesh&,
             const labelList&
         );
 
         //- Construct from mesh. Uses all cells in mesh.
-        treeDataCell(const bool cacheBb, const polyMesh&);
+        treeDataCell(const bool cacheBb, const primitiveMesh&);
 
 
     // Member Functions
@@ -104,7 +104,7 @@ public:
                 return cellLabels_;
             }
 
-            const polyMesh& mesh() const
+            const primitiveMesh& mesh() const
             {
                 return mesh_;
             }
diff --git a/src/meshTools/meshSearch/meshSearch.C b/src/meshTools/meshSearch/meshSearch.C
index 58c1d55e40..8657f5c242 100644
--- a/src/meshTools/meshSearch/meshSearch.C
+++ b/src/meshTools/meshSearch/meshSearch.C
@@ -22,14 +22,11 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
-#include "polyMesh.H"
 #include "meshSearch.H"
-#include "triPointRef.H"
-#include "octree.H"
+#include "polyMesh.H"
+#include "indexedOctree.H"
 #include "pointIndexHit.H"
 #include "DynamicList.H"
 #include "demandDrivenData.H"
@@ -43,14 +40,75 @@ Foam::scalar Foam::meshSearch::tol_ = 1E-3;
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+bool Foam::meshSearch::findNearer
+(
+    const point& sample,
+    const pointField& points,
+    label& nearestI,
+    scalar& nearestDistSqr
+)
+{
+    bool nearer = false;
+
+    forAll(points, pointI)
+    {
+        scalar distSqr = magSqr(points[pointI] - sample);
+
+        if (distSqr < nearestDistSqr)
+        {
+            nearestDistSqr = distSqr;
+            nearestI = pointI;
+            nearer = true;
+        }
+    }
+
+    return nearer;
+}
+
+
+bool Foam::meshSearch::findNearer
+(
+    const point& sample,
+    const pointField& points,
+    const labelList& indices,
+    label& nearestI,
+    scalar& nearestDistSqr
+)
+{
+    bool nearer = false;
+
+    forAll(indices, i)
+    {
+        label pointI = indices[i];
+
+        scalar distSqr = magSqr(points[pointI] - sample);
+
+        if (distSqr < nearestDistSqr)
+        {
+            nearestDistSqr = distSqr;
+            nearestI = pointI;
+            nearer = true;
+        }
+    }
+
+    return nearer;
+}
+
+
 // tree based searching
 Foam::label Foam::meshSearch::findNearestCellTree(const point& location) const
 {
-    treeBoundBox tightest(treeBoundBox::greatBox);
+    const indexedOctree<treeDataPoint>& tree = cellCentreTree();
 
-    scalar tightestDist = GREAT;
+    scalar span = mag(tree.bb().max() - tree.bb().min());
 
-    return cellCentreTree().findNearest(location, tightest, tightestDist);
+    pointIndexHit info = tree.findNearest(location, Foam::sqr(span));
+
+    if (!info.hit())
+    {
+        info = tree.findNearest(location, Foam::sqr(GREAT));
+    }
+    return info.index();
 }
 
 
@@ -59,21 +117,18 @@ Foam::label Foam::meshSearch::findNearestCellLinear(const point& location) const
 {
     const vectorField& centres = mesh_.cellCentres();
 
-    label nearestCelli = 0;
-    scalar minProximity = magSqr(centres[0] - location);
+    label nearestIndex = 0;
+    scalar minProximity = magSqr(centres[nearestIndex] - location);
 
-    forAll(centres, celli)
-    {
-        scalar proximity = magSqr(centres[celli] - location);
+    findNearer
+    (
+        location,
+        centres,
+        nearestIndex,
+        minProximity
+    );
 
-        if (proximity < minProximity)
-        {
-            nearestCelli = celli;
-            minProximity = proximity;
-        }
-    }
-
-    return nearestCelli;
+    return nearestIndex;
 }
 
 
@@ -92,40 +147,145 @@ Foam::label Foam::meshSearch::findNearestCellWalk
         )   << "illegal seedCell:" << seedCellI << exit(FatalError);
     }
 
-    const vectorField& centres = mesh_.cellCentres();
-    const labelListList& cc = mesh_.cellCells();
-
-
     // Walk in direction of face that decreases distance
 
-    label curCell = seedCellI;
-    scalar distanceSqr = magSqr(centres[curCell] - location);
+    label curCellI = seedCellI;
+    scalar distanceSqr = magSqr(mesh_.cellCentres()[curCellI] - location);
 
     bool closer;
 
     do
     {
-        closer = false;
-
-        // set the current list of neighbouring cells
-        const labelList& neighbours = cc[curCell];
-
-        forAll (neighbours, nI)
-        {
-            scalar curDistSqr = magSqr(centres[neighbours[nI]] - location);
-
-            // search through all the neighbours.
-            // If the cell is closer, reset current cell and distance
-            if (curDistSqr < distanceSqr)
-            {
-                distanceSqr = curDistSqr;
-                curCell = neighbours[nI];
-                closer = true;    // a closer neighbour has been found
-            }
-        }
+        // Try neighbours of curCellI
+        closer = findNearer
+        (
+            location,
+            mesh_.cellCentres(),
+            mesh_.cellCells()[curCellI],
+            curCellI,
+            distanceSqr
+        );
     } while (closer);
 
-    return curCell;
+    return curCellI;
+}
+
+
+// tree based searching
+Foam::label Foam::meshSearch::findNearestFaceTree(const point& location) const
+{
+    // Search nearest cell centre.
+    const indexedOctree<treeDataPoint>& tree = cellCentreTree();
+
+    scalar span = mag(tree.bb().max() - tree.bb().min());
+
+    // Search with decent span
+    pointIndexHit info = tree.findNearest(location, Foam::sqr(span));
+
+    if (!info.hit())
+    {
+        // Search with desparate span
+        info = tree.findNearest(location, Foam::sqr(GREAT));
+    }
+
+
+    // Now check any of the faces of the nearest cell
+    const vectorField& centres = mesh_.faceCentres();
+    const cell& ownFaces = mesh_.cells()[info.index()];
+
+    label nearestFaceI = ownFaces[0];
+    scalar minProximity = magSqr(centres[nearestFaceI] - location);
+
+    findNearer
+    (
+        location,
+        centres,
+        ownFaces,
+        nearestFaceI,
+        minProximity
+    );
+
+    return nearestFaceI;
+}
+
+
+// linear searching
+Foam::label Foam::meshSearch::findNearestFaceLinear(const point& location) const
+{
+    const vectorField& centres = mesh_.faceCentres();
+
+    label nearestFaceI = 0;
+    scalar minProximity = magSqr(centres[nearestFaceI] - location);
+
+    findNearer
+    (
+        location,
+        centres,
+        nearestFaceI,
+        minProximity
+    );
+
+    return nearestFaceI;
+}
+
+
+// walking from seed
+Foam::label Foam::meshSearch::findNearestFaceWalk
+(
+    const point& location,
+    const label seedFaceI
+) const
+{
+    if (seedFaceI < 0)
+    {
+        FatalErrorIn
+        (
+            "meshSearch::findNearestFaceWalk(const point&, const label)"
+        )   << "illegal seedFace:" << seedFaceI << exit(FatalError);
+    }
+
+    const vectorField& centres = mesh_.faceCentres();
+
+
+    // Walk in direction of face that decreases distance
+
+    label curFaceI = seedFaceI;
+    scalar distanceSqr = magSqr(centres[curFaceI] - location);
+
+    while (true)
+    {
+        label betterFaceI = curFaceI;
+
+        findNearer
+        (
+            location,
+            centres,
+            mesh_.cells()[mesh_.faceOwner()[curFaceI]],
+            betterFaceI,
+            distanceSqr
+        );
+
+        if (mesh_.isInternalFace(curFaceI))
+        {
+            findNearer
+            (
+                location,
+                centres,
+                mesh_.cells()[mesh_.faceNeighbour()[curFaceI]],
+                betterFaceI,
+                distanceSqr
+            );
+        }
+
+        if (betterFaceI == curFaceI)
+        {
+            break;
+        }
+
+        curFaceI = betterFaceI;
+    }
+
+    return curFaceI;
 }
 
 
@@ -180,12 +340,11 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
 
     const face& f =  mesh_.faces()[curFaceI];
 
-    scalar minDist =
-        f.nearestPoint
-        (
-            location,
-            mesh_.points()
-        ).distance();
+    scalar minDist = f.nearestPoint
+    (
+        location,
+        mesh_.points()
+    ).distance();
 
     bool closer;
 
@@ -202,8 +361,7 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
 
         forAll (myEdges, myEdgeI)
         {
-            const labelList& neighbours =
-                mesh_.edgeFaces()[myEdges[myEdgeI]];
+            const labelList& neighbours = mesh_.edgeFaces()[myEdges[myEdgeI]];
 
             // Check any face which uses edge, is boundary face and
             // is not curFaceI itself.
@@ -220,12 +378,11 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
                 {
                     const face& f =  mesh_.faces()[faceI];
 
-                    pointHit curHit =
-                        f.nearestPoint
-                        (
-                            location,
-                            mesh_.points()
-                        );
+                    pointHit curHit = f.nearestPoint
+                    (
+                        location,
+                        mesh_.points()
+                    );
 
                     // If the face is closer, reset current face and distance
                     if (curHit.distance() < minDist)
@@ -283,9 +440,11 @@ Foam::meshSearch::~meshSearch()
     clearOut();
 }
 
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const Foam::octree<Foam::octreeDataFace>& Foam::meshSearch::boundaryTree() const
+const Foam::indexedOctree<Foam::treeDataFace>& Foam::meshSearch::boundaryTree()
+ const
 {
     if (!boundaryTreePtr_)
     {
@@ -294,17 +453,28 @@ const Foam::octree<Foam::octreeDataFace>& Foam::meshSearch::boundaryTree() const
         //
 
         // all boundary faces (not just walls)
-        octreeDataFace shapes(mesh_);
+        labelList bndFaces(mesh_.nFaces()-mesh_.nInternalFaces());
+        forAll(bndFaces, i)
+        {
+            bndFaces[i] = mesh_.nInternalFaces() + i;
+        }
 
         treeBoundBox overallBb(mesh_.points());
 
-        boundaryTreePtr_ = new octree<octreeDataFace>
+        Random rndGen(123456);
+
+        boundaryTreePtr_ = new indexedOctree<treeDataFace>
         (
-            overallBb,  // overall search domain
-            shapes,     // all information needed to do checks on cells
-            1,          // min levels
-            20.0,       // maximum ratio of cubes v.s. cells
-            10.0
+            treeDataFace    // all information needed to search faces
+            (
+                false,                      // do not cache bb
+                mesh_,
+                bndFaces                    // boundary faces only
+            ),
+            overallBb.extend(rndGen, 1E-3), // overall search domain
+            8,                              // maxLevel
+            10,                             // leafsize
+            3.0                             // duplicity
         );
     }
 
@@ -312,7 +482,8 @@ const Foam::octree<Foam::octreeDataFace>& Foam::meshSearch::boundaryTree() const
 }
 
 
-const Foam::octree<Foam::octreeDataCell>& Foam::meshSearch::cellTree() const
+const Foam::indexedOctree<Foam::treeDataCell>& Foam::meshSearch::cellTree()
+ const
 {
     if (!cellTreePtr_)
     {
@@ -320,17 +491,19 @@ const Foam::octree<Foam::octreeDataCell>& Foam::meshSearch::cellTree() const
         // Construct tree
         //
 
-        octreeDataCell shapes(mesh_);
-
         treeBoundBox overallBb(mesh_.points());
 
-        cellTreePtr_ = new octree<octreeDataCell>
+        cellTreePtr_ = new indexedOctree<treeDataCell>
         (
+            treeDataCell
+            (
+                false,  // not cache bb
+                mesh_
+            ),
             overallBb,  // overall search domain
-            shapes,     // all information needed to do checks on cells
-            1,          // min levels
-            20.0,       // maximum ratio of cubes v.s. cells
-            2.0
+            8,      // maxLevel
+            10,     // leafsize
+            3.0     // duplicity
         );
     }
 
@@ -339,8 +512,8 @@ const Foam::octree<Foam::octreeDataCell>& Foam::meshSearch::cellTree() const
 }
 
 
-const Foam::octree<Foam::octreeDataPoint>& Foam::meshSearch::cellCentreTree()
-    const
+const Foam::indexedOctree<Foam::treeDataPoint>&
+ Foam::meshSearch::cellCentreTree() const
 {
     if (!cellCentreTreePtr_)
     {
@@ -348,17 +521,14 @@ const Foam::octree<Foam::octreeDataPoint>& Foam::meshSearch::cellCentreTree()
         // Construct tree
         //
 
-        // shapes holds reference to cellCentres!
-        octreeDataPoint shapes(mesh_.cellCentres());
-
         treeBoundBox overallBb(mesh_.cellCentres());
 
-        cellCentreTreePtr_ = new octree<octreeDataPoint>
+        cellCentreTreePtr_ = new indexedOctree<treeDataPoint>
         (
+            treeDataPoint(mesh_.cellCentres()),
             overallBb,  // overall search domain
-            shapes,     // all information needed to do checks on cells
-            1,          // min levels
-            20.0,       // maximum ratio of cubes v.s. cells
+            10,         // max levels
+            10.0,       // maximum ratio of cubes v.s. cells
             100.0       // max. duplicity; n/a since no bounding boxes.
         );
     }
@@ -396,15 +566,14 @@ bool Foam::meshSearch::pointInCell(const point& p, label cellI) const
 
             forAll(f, fp)
             {
-                pointHit inter = 
-                    f.ray
-                    (
-                        ctr,
-                        dir,
-                        mesh_.points(),
-                        intersection::HALF_RAY,
-                        intersection::VECTOR
-                    );
+                pointHit inter = f.ray
+                (
+                    ctr,
+                    dir,
+                    mesh_.points(),
+                    intersection::HALF_RAY,
+                    intersection::VECTOR
+                );
 
                 if (inter.hit())
                 {
@@ -484,6 +653,31 @@ Foam::label Foam::meshSearch::findNearestCell
 }
 
 
+Foam::label Foam::meshSearch::findNearestFace
+(
+    const point& location,
+    const label seedFaceI,
+    const bool useTreeSearch
+) const
+{
+    if (seedFaceI == -1)
+    {
+        if (useTreeSearch)
+        {
+            return findNearestFaceTree(location);
+        }
+        else
+        {
+            return findNearestFaceLinear(location);
+        }
+    }
+    else
+    {
+        return findNearestFaceWalk(location, seedFaceI);
+    }
+}
+
+
 Foam::label Foam::meshSearch::findCell
 (
     const point& location,
@@ -607,19 +801,26 @@ Foam::label Foam::meshSearch::findNearestBoundaryFace
     {
         if (useTreeSearch)
         {
-            treeBoundBox tightest(treeBoundBox::greatBox);
+            const indexedOctree<treeDataFace>& tree =  boundaryTree();
 
-            scalar tightestDist = GREAT;
+            scalar span = mag(tree.bb().max() - tree.bb().min());
 
-            label index =
-                boundaryTree().findNearest
+            pointIndexHit info = boundaryTree().findNearest
+            (
+                location,
+                Foam::sqr(span)
+            );
+
+            if (!info.hit())
+            {
+                info = boundaryTree().findNearest
                 (
                     location,
-                    tightest,
-                    tightestDist
+                    Foam::sqr(GREAT)
                 );
+            }
 
-            return boundaryTree().shapes().meshFaces()[index];
+            return tree.shapes().faceLabels()[info.index()];
         }
         else
         {
@@ -670,7 +871,7 @@ Foam::pointIndexHit Foam::meshSearch::intersection
     if (curHit.hit())
     {
         // Change index into octreeData into face label
-        curHit.setIndex(boundaryTree().shapes().meshFaces()[curHit.index()]);
+        curHit.setIndex(boundaryTree().shapes().faceLabels()[curHit.index()]);
     }
     return curHit;
 }
@@ -733,7 +934,9 @@ Foam::List<Foam::pointIndexHit> Foam::meshSearch::intersections
 
 bool Foam::meshSearch::isInside(const point& p) const
 {
-    return boundaryTree().getSampleType(p) == octree<octreeDataFace>::INSIDE;
+    return
+        boundaryTree().getVolumeType(p)
+     == indexedOctree<treeDataFace>::INSIDE;
 }
 
 
diff --git a/src/meshTools/meshSearch/meshSearch.H b/src/meshTools/meshSearch/meshSearch.H
index c50c0b8bb0..951f41c50f 100644
--- a/src/meshTools/meshSearch/meshSearch.H
+++ b/src/meshTools/meshSearch/meshSearch.H
@@ -26,7 +26,8 @@ Class
     Foam::meshSearch
 
 Description
-    Various searches on polyMesh; uses (demand driven) octree to search.
+    Various (local, not parallel) searches on polyMesh;
+    uses (demand driven) octree to search.
 
 SourceFiles
     meshSearch.C
@@ -36,11 +37,10 @@ SourceFiles
 #ifndef meshSearch_H
 #define meshSearch_H
 
-#include "octreeDataCell.H"
-#include "octreeDataFace.H"
-#include "octreeDataPoint.H"
+#include "treeDataCell.H"
+#include "treeDataFace.H"
+#include "treeDataPoint.H"
 #include "pointIndexHit.H"
-#include "className.H"
 #include "Cloud.H"
 #include "passiveParticle.H"
 
@@ -71,46 +71,77 @@ class meshSearch
 
         //- demand driven octrees
 
-        mutable octree<octreeDataFace>* boundaryTreePtr_;
-        mutable octree<octreeDataCell>* cellTreePtr_;
-        mutable octree<octreeDataPoint>* cellCentreTreePtr_;
+        mutable indexedOctree<treeDataFace>* boundaryTreePtr_;
+        mutable indexedOctree<treeDataCell>* cellTreePtr_;
+        mutable indexedOctree<treeDataPoint>* cellCentreTreePtr_;
 
 
     // Private Member Functions
 
-        //- nearest cell centre using octree
-        label findNearestCellTree(const point& location) const;
-
-        //- nearest cell centre going through all cells
-        label findNearestCellLinear(const point& location) const;
-
-        //- walk from seed. Does not 'go around' boundary, just returns
-        //  last cell before boundary.
-        label findNearestCellWalk
+        //- Updates nearestI, nearestDistSqr from any closer ones.
+        static bool findNearer
         (
-            const point& location,
-            const label seedCellI
-        ) const;
+            const point& sample,
+            const pointField& points,
+            label& nearestI,
+            scalar& nearestDistSqr
+        );
 
-        //- cell containing location. Linear search.
-        label findCellLinear(const point& location) const;
-
-        //- walk from seed to find nearest boundary face. Gets stuck in
-        //  local minimum.
-        label findNearestBoundaryFaceWalk
+        //- Updates nearestI, nearestDistSqr from any selected closer ones.
+        static bool findNearer
         (
-            const point& location,
-            const label seedFaceI
-        ) const;
+            const point& sample,
+            const pointField& points,
+            const labelList& indices,
+            label& nearestI,
+            scalar& nearestDistSqr
+        );
 
-        //- Calculate offset vector in direction dir with as length a fraction
-        //  of the cell size (of the cell straddling boundary face)
-        vector offset
-        (
-            const point& bPoint,
-            const label bFaceI,
-            const vector& dir
-        ) const;
+
+        // Cells
+
+            //- nearest cell centre using octree
+            label findNearestCellTree(const point&) const;
+
+            //- nearest cell centre going through all cells
+            label findNearestCellLinear(const point&) const;
+
+            //- walk from seed. Does not 'go around' boundary, just returns
+            //  last cell before boundary.
+            label findNearestCellWalk(const point&, const label) const;
+
+            //- cell containing location. Linear search.
+            label findCellLinear(const point&) const;
+
+
+        // Cells
+
+            label findNearestFaceTree(const point&) const;
+
+            label findNearestFaceLinear(const point&) const;
+
+            label findNearestFaceWalk(const point&, const label) const;
+
+
+
+        // Boundary faces
+
+            //- walk from seed to find nearest boundary face. Gets stuck in
+            //  local minimum.
+            label findNearestBoundaryFaceWalk
+            (
+                const point& location,
+                const label seedFaceI
+            ) const;
+
+            //- Calculate offset vector in direction dir with as length a
+            //  fraction of the cell size (of the cell straddling boundary face)
+            vector offset
+            (
+                const point& bPoint,
+                const label bFaceI,
+                const vector& dir
+            ) const;
 
 
         //- Disallow default bitwise copy construct
@@ -154,13 +185,13 @@ public:
 
             //- Get (demand driven) reference to octree holding all
             //  boundary faces
-            const octree<octreeDataFace>& boundaryTree() const;
+            const indexedOctree<treeDataFace>& boundaryTree() const;
 
             //- Get (demand driven) reference to octree holding all cells
-            const octree<octreeDataCell>& cellTree() const;
+            const indexedOctree<treeDataCell>& cellTree() const;
 
             //- Get (demand driven) reference to octree holding all cell centres
-            const octree<octreeDataPoint>& cellCentreTree() const;
+            const indexedOctree<treeDataPoint>& cellCentreTree() const;
 
 
         // Queries
@@ -181,6 +212,13 @@ public:
                 const bool useTreeSearch = true
             ) const;
 
+            label findNearestFace
+            (
+                const point& location,
+                const label seedFaceI = -1,
+                const bool useTreeSearch = true
+            ) const;
+
             //- Find cell containing (using pointInCell) location.
             //  If seed provided walks and falls back to linear/tree search.
             //  (so handles holes correctly)s
@@ -206,7 +244,7 @@ public:
             //- Find first intersection of boundary in segment [pStart, pEnd]
             //  (so inclusive of endpoints). Always octree for now
             pointIndexHit intersection(const point& pStart, const point& pEnd)
-                const;
+            const;
 
             //- Find all intersections of boundary within segment pStart .. pEnd
             //  Always octree for now
diff --git a/src/meshTools/octree/treeBoundBox.H b/src/meshTools/octree/treeBoundBox.H
index 02dbac38af..1805d88b76 100644
--- a/src/meshTools/octree/treeBoundBox.H
+++ b/src/meshTools/octree/treeBoundBox.H
@@ -314,6 +314,7 @@ public:
 
             //- Return slightly wider bounding box
             //  Extends all dimensions with s*span*Random::scalar01()
+            //  and guarantees in any direction s*mag(span) minimum width
             inline treeBoundBox extend(Random&, const scalar s) const;
 
     // Friend Operators
diff --git a/src/meshTools/octree/treeBoundBoxI.H b/src/meshTools/octree/treeBoundBoxI.H
index c444b2808a..36143fedf1 100644
--- a/src/meshTools/octree/treeBoundBoxI.H
+++ b/src/meshTools/octree/treeBoundBoxI.H
@@ -437,7 +437,15 @@ inline treeBoundBox treeBoundBox::extend(Random& rndGen, const scalar s) const
 {
     treeBoundBox bb(*this);
 
-    const vector span(bb.max() - bb.min());
+    vector span(bb.max() - bb.min());
+
+    // Make 3D
+    scalar magSpan = Foam::mag(span);
+
+    for (direction dir = 0; dir < vector::nComponents; dir++)
+    {
+        span[dir] = Foam::max(s*magSpan, span[dir]);
+    }
 
     bb.min() -= cmptMultiply(s*rndGen.vector01(), span);
     bb.max() += cmptMultiply(s*rndGen.vector01(), span);
diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.C b/src/meshTools/searchableSurface/triSurfaceMesh.C
index 61b760699d..a426741960 100644
--- a/src/meshTools/searchableSurface/triSurfaceMesh.C
+++ b/src/meshTools/searchableSurface/triSurfaceMesh.C
@@ -230,7 +230,7 @@ const Foam::indexedOctree<Foam::treeDataTriSurface>&
             new indexedOctree<treeDataTriSurface>
             (
                 treeDataTriSurface(*this),
-                bb.extend(rndGen, 1E-3),    // slightly randomize bb
+                bb.extend(rndGen, 1E-4),    // slightly randomize bb
                 10,     // maxLevel
                 10,     // leafsize
                 3.0     // duplicity
@@ -274,7 +274,7 @@ const Foam::indexedOctree<Foam::treeDataEdge>&
                     localPoints(),  // points
                     bEdges          // selected edges
                 ),
-                bb.extend(rndGen, 1E-3),    // slightly randomize bb
+                bb.extend(rndGen, 1E-4),    // slightly randomize bb
                 8,      // maxLevel
                 10,     // leafsize
                 3.0     // duplicity
diff --git a/src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.C b/src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.C
index 0bccf0906d..30f4f3d75a 100644
--- a/src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.C
+++ b/src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.C
@@ -31,8 +31,6 @@ Description
 #include "meshSearch.H"
 #include "triSurface.H"
 #include "triSurfaceSearch.H"
-#include "octree.H"
-#include "octreeDataTriSurface.H"
 #include "cellClassification.H"
 
 #include "addToRunTimeSelectionTable.H"
diff --git a/src/meshTools/sets/pointSources/surfaceToPoint/surfaceToPoint.C b/src/meshTools/sets/pointSources/surfaceToPoint/surfaceToPoint.C
index d38de31a1d..e3c994cf31 100644
--- a/src/meshTools/sets/pointSources/surfaceToPoint/surfaceToPoint.C
+++ b/src/meshTools/sets/pointSources/surfaceToPoint/surfaceToPoint.C
@@ -30,8 +30,6 @@ Description
 #include "polyMesh.H"
 #include "meshSearch.H"
 #include "triSurfaceSearch.H"
-#include "octree.H"
-#include "octreeDataTriSurface.H"
 
 #include "addToRunTimeSelectionTable.H"
 
diff --git a/tutorials/interDyMFoam/sloshingTank2D/constant/polyMesh/blockMeshDict.m4 b/tutorials/interDyMFoam/sloshingTank2D/constant/polyMesh/blockMeshDict.m4
index ad753abc56..0473efbc15 100644
--- a/tutorials/interDyMFoam/sloshingTank2D/constant/polyMesh/blockMeshDict.m4
+++ b/tutorials/interDyMFoam/sloshingTank2D/constant/polyMesh/blockMeshDict.m4
@@ -8,7 +8,7 @@
 FoamFile
 {
     version     2.0;
-    format      ascii;
+    `format'      ascii;
     class       dictionary;
     object      blockMeshDict;
 }
@@ -16,7 +16,7 @@ FoamFile
 // General m4 macros
 
 changecom(//)changequote([,])
-define(calc, [esyscmd(perl -e 'use Math::Trig; use POSIX; print ($1)')])
+define(calc, [esyscmd(perl -e 'use Math::Trig; use POSIX; printf ($1)')])
 define(VCOUNT, 0)
 define(vlabel, [[// ]Vertex $1 = VCOUNT define($1, VCOUNT)define([VCOUNT], incr(VCOUNT))])
 
diff --git a/tutorials/oodles/pitzDailyDirectMapped/constant/polyMesh/boundary b/tutorials/oodles/pitzDailyDirectMapped/constant/polyMesh/boundary
index 89b5e56c17..c67eadf385 100644
--- a/tutorials/oodles/pitzDailyDirectMapped/constant/polyMesh/boundary
+++ b/tutorials/oodles/pitzDailyDirectMapped/constant/polyMesh/boundary
@@ -20,8 +20,9 @@ inlet
 {
     nFaces          30;
     startFace       27238;
-    type            directMappedPatch;
-    offset          ( 0.05 0 0 );
+    offset          ( 0.0495 0 0 );
+    sampleMode      nearestCell;
+    samplePatch     none;
 }
 
 outlet          
diff --git a/tutorials/oodles/pitzDailyDirectMapped/system/changeDictionaryDict b/tutorials/oodles/pitzDailyDirectMapped/system/changeDictionaryDict
index 282852ed87..7ae2253dbf 100644
--- a/tutorials/oodles/pitzDailyDirectMapped/system/changeDictionaryDict
+++ b/tutorials/oodles/pitzDailyDirectMapped/system/changeDictionaryDict
@@ -23,6 +23,8 @@ dictionaryReplacement
         {
             type            directMappedPatch;
             offset          (0.0495 0 0);
+            sampleMode      nearestCell;
+            samplePatch     none;
         }
     }
 }