diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index cc975eb970..205c767600 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@@ -274,6 +274,14 @@ addLayersControls
 
     // Create buffer region for new layer terminations
     nBufferCellsNoExtrude 0;
+
+
+    // Overall max number of layer addition iterations
+    nLayerIter 50;
+
+    // Max number of iterations after which relaxed meshQuality controls
+    // get used.
+    nRelaxedIter 20;
 }
 
 
@@ -335,6 +343,16 @@ meshQualityControls
     nSmoothScale 4;
     //- amount to scale back displacement at error points
     errorReduction 0.75;
+
+
+
+    // Optional : some meshing phases allow usage of relaxed rules.
+    // See e.g. addLayersControls::nRelaxedIter.
+    relaxed
+    {
+        //- Maximum non-orthogonality allowed. Set to 180 to disable.
+        maxNonOrtho 75;
+    }
 }
 
 
diff --git a/applications/utilities/surface/surfaceCheck/surfaceCheck.C b/applications/utilities/surface/surfaceCheck/surfaceCheck.C
index e7b7890074..4731ed3ec2 100644
--- a/applications/utilities/surface/surfaceCheck/surfaceCheck.C
+++ b/applications/utilities/surface/surfaceCheck/surfaceCheck.C
@@ -171,11 +171,11 @@ int main(int argc, char *argv[])
 
     argList::validArgs.clear();
     argList::validArgs.append("surface file");
-    argList::validOptions.insert("noSelfIntersection", "");
+    argList::validOptions.insert("checkSelfIntersection", "");
     argList::validOptions.insert("verbose", "");
     argList args(argc, argv);
 
-    bool checkSelfIntersection = !args.options().found("noSelfIntersection");
+    bool checkSelfIntersection = args.options().found("checkSelfIntersection");
     bool verbose = args.options().found("verbose");
 
     fileName surfFileName(args.additionalArgs()[0]);
diff --git a/src/OpenFOAM/meshes/boundBox/boundBox.H b/src/OpenFOAM/meshes/boundBox/boundBox.H
index e41b6efa46..38b852ef66 100644
--- a/src/OpenFOAM/meshes/boundBox/boundBox.H
+++ b/src/OpenFOAM/meshes/boundBox/boundBox.H
@@ -174,7 +174,7 @@ public:
 
         // Query
 
-            //- Completely contains other boundingBox? (inside or on edge)
+            //- Overlaps/touches boundingBox?
             bool overlaps(const boundBox& bb) const
             {
                 return
diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.H b/src/OpenFOAM/meshes/meshShapes/face/face.H
index 07dd39f93f..25a4f8fbce 100644
--- a/src/OpenFOAM/meshes/meshShapes/face/face.H
+++ b/src/OpenFOAM/meshes/meshShapes/face/face.H
@@ -222,16 +222,16 @@ public:
         ) const;
 
         //- Fast intersection with a ray.
-        //  For a hit, the pointHit.distance() is the line parameter t :
-        //  intersection=p+t*q. Only defined for FULL_RAY or
-        //  HALF_RAY.
+        //  Does face-center decomposition and returns triangle intersection
+        //  point closest to p. See triangle::intersection for details.
         pointHit intersection
         (
             const point& p,
             const vector& q,
             const point& ctr,
             const pointField& meshPoints,
-            const intersection::algorithm alg
+            const intersection::algorithm alg,
+            const scalar tol = 0.0
         ) const;
 
         //- Return nearest point to face
diff --git a/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C b/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C
index b514af7749..c5f9a80250 100644
--- a/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C
+++ b/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C
@@ -136,7 +136,8 @@ Foam::pointHit Foam::face::intersection
     const vector& q,
     const point& ctr,
     const pointField& meshPoints,
-    const intersection::algorithm alg
+    const intersection::algorithm alg,
+    const scalar tol
 ) const
 {
     scalar nearestHitDist = VGREAT;
@@ -154,7 +155,7 @@ Foam::pointHit Foam::face::intersection
             meshPoints[f[pI]],
             meshPoints[f[fcIndex(pI)]],
             ctr
-        ).intersection(p, q, alg);
+        ).intersection(p, q, alg, tol);
 
         if (curHit.hit())
         {
diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
index b24149300f..1165046636 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
@@ -153,7 +153,7 @@ public:
             inline scalar sweptVol(const triangle& t) const;
 
             //- Return point intersection with a ray.
-            //  For a hit, the distance is signed.  Positive number
+            //  For a hit, the distance is signed. Positive number
             //  represents the point in front of triangle.
             //  In case of miss pointHit is set to nearest point
             //  on triangle and its distance to the distance between
@@ -166,15 +166,20 @@ public:
                 const intersection::direction dir = intersection::VECTOR
             ) const;
 
-            //- Fast intersection with a ray.
-            //  For a hit, the pointHit.distance() is the line parameter t :
-            //  intersection=p+t*q. Only defined for VISIBLE, FULL_RAY or
-            //  HALF_RAY.
+            //- Fast intersection with a ray. Distance is normalised distance
+            //  to the intersection (normalised with the vector magnitude).
+            //  For a hit, the distance is signed. Positive number
+            //  represents the point in front of triangle. In case of miss
+            //  pointHit position is undefined.
+            //  Only defined for VISIBLE, FULL_RAY or
+            //  HALF_RAY. tol increases the virtual size of the triangle
+            //  by a relative factor.
             inline pointHit intersection
             (
                 const point& p,
                 const vector& q,
-                const intersection::algorithm alg
+                const intersection::algorithm alg,
+                const scalar tol = 0.0
             ) const;
 
             //- Return nearest point to p on triangle
diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
index ad7ed2ea8b..33b3cab15e 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
@@ -439,7 +439,8 @@ inline pointHit triangle<Point, PointRef>::intersection
 (
     const point& orig,
     const vector& dir,
-    const intersection::algorithm alg
+    const intersection::algorithm alg,
+    const scalar tol
 ) const
 {
     const vector edge1 = b_ - a_;
@@ -457,99 +458,61 @@ inline pointHit triangle<Point, PointRef>::intersection
     if (alg == intersection::VISIBLE)
     {
         // Culling branch
-        if (det < SMALL)
+        if (det < ROOTVSMALL)
         {
-            // return miss
+            // Ray on wrong side of triangle. Return miss
             return intersection;
         }
-        /* calculate distance from a_ to ray origin */
-        const vector tVec = orig-a_;
-
-        /* calculate U parameter and test bounds */
-        scalar u = tVec & pVec;
-
-        if (u < 0.0 || u > det)
-        {
-            // return miss
-            return intersection;
-        }
-
-        /* prepare to test V parameter */
-        const vector qVec = tVec ^ edge1;
-
-        /* calculate V parameter and test bounds */
-        scalar v = dir & qVec;
-
-        if (v < 0.0 || u + v > det)
-        {
-            // return miss
-            return intersection;
-        }
-
-        /* calculate t, scale parameters, ray intersects triangle */
-        scalar t = edge2 & qVec;
-        scalar inv_det = 1.0 / det;
-        t *= inv_det;
-        u *= inv_det;
-        v *= inv_det;
-
-        intersection.setHit();
-        intersection.setPoint(a_ + u*edge1 + v*edge2);
-        intersection.setDistance(t);
     }
     else if (alg == intersection::HALF_RAY || alg == intersection::FULL_RAY)
     {
         // Non-culling branch
-        if (det > -SMALL && det < SMALL)
+        if (det > -ROOTVSMALL && det < ROOTVSMALL)
         {
-            // return miss
+            // Ray parallel to triangle. Return miss
             return intersection;
         }
-        const scalar inv_det = 1.0 / det;
-
-        /* calculate distance from a_ to ray origin */
-        const vector tVec = orig - a_;
-        /* calculate U parameter and test bounds */
-        const scalar u = (tVec & pVec)*inv_det;
-
-        if (u < 0.0 || u > 1.0)
-        {
-            // return miss
-            return intersection;
-        }
-        /* prepare to test V parameter */
-        const vector qVec = tVec ^ edge1;
-        /* calculate V parameter and test bounds */
-        const scalar v = (dir & qVec) * inv_det;
-
-        if (v < 0.0 || u + v > 1.0)
-        {
-            // return miss
-            return intersection;
-        }
-        /* calculate t, ray intersects triangle */
-        const scalar t = (edge2 & qVec) * inv_det;
-
-        if (alg == intersection::HALF_RAY && t < 0)
-        {
-            // return miss
-            return intersection;
-        }
-
-        intersection.setHit();
-        intersection.setPoint(a_ + u*edge1 + v*edge2);
-        intersection.setDistance(t);
     }
-    else
+
+    const scalar inv_det = 1.0 / det;
+
+    /* calculate distance from a_ to ray origin */
+    const vector tVec = orig-a_;
+
+    /* calculate U parameter and test bounds */
+    const scalar u = (tVec & pVec)*inv_det;
+
+    if (u < -tol || u > 1.0+tol)
     {
-        FatalErrorIn
-        (
-            "triangle<Point, PointRef>::intersection(const point&"
-            ", const vector&, const intersection::algorithm)"
-        )   << "intersection only defined for VISIBLE, FULL_RAY or HALF_RAY"
-            << abort(FatalError);
+        // return miss
+        return intersection;
     }
 
+    /* prepare to test V parameter */
+    const vector qVec = tVec ^ edge1;
+
+    /* calculate V parameter and test bounds */
+    const scalar v = (dir & qVec) * inv_det;
+
+    if (v < -tol || u + v > 1.0+tol)
+    {
+        // return miss
+        return intersection;
+    }
+
+    /* calculate t, scale parameters, ray intersects triangle */
+    const scalar t = (edge2 & qVec) * inv_det;
+
+    if (alg == intersection::HALF_RAY && t < -tol)
+    {
+        // Wrong side of orig. Return miss
+        return intersection;
+    }
+
+    intersection.setHit();
+    intersection.setPoint(a_ + u*edge1 + v*edge2);
+    intersection.setDistance(t);
+
     return intersection;
 }
 
@@ -622,7 +585,7 @@ inline bool triangle<Point, PointRef>::classify
 
     bool hit = false;
 
-    if (Foam::mag(u1) < SMALL)
+    if (Foam::mag(u1) < ROOTVSMALL)
     {
         beta = u0/u2;
 
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.H b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.H
index 9ae910b5a7..c0b7668d0a 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.H
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.H
@@ -38,7 +38,6 @@ SourceFiles
 
 #include "autoPtr.H"
 #include "dictionary.H"
-#include "boolList.H"
 #include "wallPoint.H"
 #include "searchableSurfaces.H"
 #include "refinementSurfaces.H"
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index 8d3c0e1e76..0859e37d33 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -2292,7 +2292,7 @@ bool Foam::autoLayerDriver::cellsUseFace
 Foam::label Foam::autoLayerDriver::checkAndUnmark
 (
     const addPatchCellLayer& addLayer,
-    const dictionary& motionDict,
+    const dictionary& meshQualityDict,
     const indirectPrimitivePatch& pp,
     const fvMesh& newMesh,
 
@@ -2304,7 +2304,7 @@ Foam::label Foam::autoLayerDriver::checkAndUnmark
     // Check the resulting mesh for errors
     Info<< nl << "Checking mesh with layer ..." << endl;
     faceSet wrongFaces(newMesh, "wrongFaces", newMesh.nFaces()/1000);
-    motionSmoother::checkMesh(false, newMesh, motionDict, wrongFaces);
+    motionSmoother::checkMesh(false, newMesh, meshQualityDict, wrongFaces);
     Info<< "Detected " << returnReduce(wrongFaces.size(), sumOp<label>())
         << " illegal faces"
         << " (concave, zero area or negative cell pyramid volume)"
@@ -2474,8 +2474,8 @@ void Foam::autoLayerDriver::mergePatchFacesUndo
         << "      (cos:" << minCos << ')' << nl
         << "    - as long as the resulting face doesn't become concave"
         << " by more than "
-        << layerParams.concaveAngle()
-        << " degrees (0=straight, 180=fully concave)" << nl
+        << layerParams.concaveAngle() << " degrees" << nl
+        << "      (0=straight, 180=fully concave)" << nl
         << endl;
 
     label nChanged = mergePatchFacesUndo(minCos, concaveCos, motionDict);
@@ -2709,8 +2709,12 @@ void Foam::autoLayerDriver::addLayers
     boolList flaggedCells;
     boolList flaggedFaces;
 
-    while (true)
+    for (label iteration = 0; iteration < layerParams.nLayerIter(); iteration++)
     {
+        Info<< nl
+            << "Layer addition iteration " << iteration << nl
+            << "--------------------------" << endl;
+
         // Make sure displacement is equal on both sides of coupled patches.
         syncPatchDisplacement
         (
@@ -2951,10 +2955,23 @@ void Foam::autoLayerDriver::addLayers
         }
 
         // Unset the extrusion at the pp.
+        const dictionary& meshQualityDict =
+        (
+            iteration < layerParams.nRelaxedIter()
+          ? motionDict
+          : motionDict.subDict("relaxed")
+        );
+
+        if (iteration >= layerParams.nRelaxedIter())
+        {
+            Info<< "Switched to relaxed meshQuality constraints." << endl;
+        }
+
+
         label nTotChanged = checkAndUnmark
         (
             addLayer,
-            motionDict,
+            meshQualityDict,
             pp,
             newMesh,
 
@@ -2976,6 +2993,8 @@ void Foam::autoLayerDriver::addLayers
         // Reset mesh points and start again
         meshMover.movePoints(oldPoints);
         meshMover.correct();
+
+        Info<< endl;
     }
 
 
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
index e891cf761c..9c4e0a17c5 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
@@ -95,7 +95,7 @@ Foam::label Foam::autoRefineDriver::readFeatureEdges
         Info<< "Refinement level " << featureLevels[i]
             << " for all cells crossed by feature " << featFileName
             << " (" << featureMeshes[i].points().size() << " points, "
-            << featureMeshes[i].edges().size() << ")." << endl;
+            << featureMeshes[i].edges().size() << " edges)." << endl;
     }
 
     Info<< "Read feature lines in = "
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
index ddf162dd0f..dfdc16e919 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
@@ -211,8 +211,24 @@ Foam::layerParameters::layerParameters
     (
         readLabel(dict.lookup("nBufferCellsNoExtrude"))
     ),
-    nSnap_(readLabel(dict.lookup("nSnap")))
-{}
+    nSnap_(readLabel(dict.lookup("nSnap"))),
+    nLayerIter_(readLabel(dict.lookup("nLayerIter"))),
+    nRelaxedIter_(labelMax)
+{
+    if (dict.found("nRelaxedIter"))
+    {
+        dict.lookup("nRelaxedIter") >> nRelaxedIter_;
+    }
+
+    if (nLayerIter_ < 0 || nRelaxedIter_ < 0)
+    {
+        FatalErrorIn("layerParameters::layerParameters(..)")
+            << "Layer iterations should be >= 0." << endl
+            << "nLayerIter:" << nLayerIter_
+            << " nRelaxedIter:" << nRelaxedIter_
+            << exit(FatalError);
+    }
+}
 
 
 // Construct from dictionary
@@ -276,8 +292,24 @@ Foam::layerParameters::layerParameters
     (
         readLabel(dict.lookup("nBufferCellsNoExtrude"))
     ),
-    nSnap_(readLabel(dict.lookup("nRelaxIter")))
+    nSnap_(readLabel(dict.lookup("nRelaxIter"))),
+    nLayerIter_(readLabel(dict.lookup("nLayerIter"))),
+    nRelaxedIter_(labelMax)
 {
+    if (dict.found("nRelaxedIter"))
+    {
+        dict.lookup("nRelaxedIter") >> nRelaxedIter_;
+    }
+    if (nLayerIter_ < 0 || nRelaxedIter_ < 0)
+    {
+        FatalErrorIn("layerParameters::layerParameters(..)")
+            << "Layer iterations should be >= 0." << endl
+            << "nLayerIter:" << nLayerIter_
+            << " nRelaxedIter:" << nRelaxedIter_
+            << exit(FatalError);
+    }
+
+
     const dictionary& layersDict = dict.subDict("layers");
 
     forAll(boundaryMesh, patchI)
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H b/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H
index d96ce97cbd..e3d3702bfd 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H
@@ -105,6 +105,10 @@ class layerParameters
 
         label nSnap_;
 
+        label nLayerIter_;
+
+        label nRelaxedIter_;
+
 
     // Private Member Functions
 
@@ -196,26 +200,20 @@ public:
                 return nGrow_;
             }
 
-            // Number of smoothing iterations of surface normals 
+            //- Number of smoothing iterations of surface normals 
             label nSmoothSurfaceNormals() const
             {
                 return nSmoothSurfaceNormals_;
             }
 
-            // Number of smoothing iterations of interior mesh movement
-            // direction  
+            //- Number of smoothing iterations of interior mesh movement
+            //  direction  
             label nSmoothNormals() const
             {
                 return nSmoothNormals_;
             }
 
-            // Smooth layer thickness over surface patches
-            label nSmoothThickness() const
-            {
-                return nSmoothThickness_;
-            }
-
-            // Stop layer growth on highly warped cells 
+            //- Stop layer growth on highly warped cells 
             scalar maxFaceThicknessRatio() const
             {
                 return maxFaceThicknessRatio_;
@@ -226,20 +224,26 @@ public:
                 return layerTerminationCos_;
             }
 
-            // Reduce layer growth where ratio thickness to medial 
-            // distance is large 
+            //- Smooth layer thickness over surface patches
+            label nSmoothThickness() const
+            {
+                return nSmoothThickness_;
+            }
+
+            //- Reduce layer growth where ratio thickness to medial 
+            //  distance is large 
             scalar maxThicknessToMedialRatio() const
             {
                 return maxThicknessToMedialRatio_;
             }
 
-            // Angle used to pick up medial axis points
+            //- Angle used to pick up medial axis points
             scalar minMedianAxisAngleCos() const
             {
                 return minMedianAxisAngleCos_;
             }
 
-            // Create buffer region for new layer terminations
+            //- Create buffer region for new layer terminations
             label nBufferCellsNoExtrude() const
             {
                 return nBufferCellsNoExtrude_;
@@ -249,6 +253,24 @@ public:
             {
                 return nSnap_;
             }
+
+            // Overall
+
+                //- Number of overall layer addition iterations
+                label nLayerIter() const
+                {
+                    return nLayerIter_;
+                }
+
+                //- Number of iterations after which relaxed motion rules
+                //  are to be used.
+                label nRelaxedIter() const
+                {
+                    return nRelaxedIter_;
+                }
+
+
+
 };
 
 
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.H b/src/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.H
index ffb8a15d05..d27c62929d 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.H
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.H
@@ -42,9 +42,7 @@ SourceFiles
 
 #include "point.H"
 #include "label.H"
-#include "scalar.H"
 #include "tensor.H"
-#include "pTraits.H"
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointDataI.H b/src/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointDataI.H
index 234d569f85..9832bbf877 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointDataI.H
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointDataI.H
@@ -22,8 +22,6 @@ 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"
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 0025f33332..06004c8a6e 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -47,7 +47,6 @@ SourceFiles
 #include "hexRef8.H"
 #include "mapPolyMesh.H"
 #include "autoPtr.H"
-#include "pointIOField.H"
 #include "labelPair.H"
 #include "indirectPrimitivePatch.H"
 #include "pointFieldsFwd.H"
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index d42ac05b22..7aabbc6321 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -328,7 +328,7 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
     // Database to pass into trackedParticle::move
     trackedParticle::trackData td(cloud, maxFeatureLevel);
 
-    // Track all particles to their end position.
+    // Track all particles to their end position (= starting feature point)
     cloud.move(td);
 
     // Reset level
diff --git a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
index e8fd361401..cee48bbe1f 100644
--- a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
+++ b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
@@ -38,7 +38,6 @@ SourceFiles
 #ifndef refinementSurfaces_H
 #define refinementSurfaces_H
 
-#include "PackedBoolList.H"
 #include "triSurfaceGeoMesh.H"
 #include "triSurfaceFields.H"
 #include "vectorList.H"
diff --git a/src/autoMesh/autoHexMesh/trackedParticle/ExactParticle.C b/src/autoMesh/autoHexMesh/trackedParticle/ExactParticle.C
index bcf92e9747..237e35276c 100644
--- a/src/autoMesh/autoHexMesh/trackedParticle/ExactParticle.C
+++ b/src/autoMesh/autoHexMesh/trackedParticle/ExactParticle.C
@@ -102,23 +102,7 @@ Foam::scalar Foam::ExactParticle<ParticleType>::trackToFace
         }
     }
 
-
-    if (hitFacei != -1)
-    {
-        if (trackFraction > 1.0)
-        {
-            // Nearest intersection beyond endPosition so we hit endPosition.
-            this->position_ = endPosition;
-            this->facei_ = -1;
-            return 1.0;
-        }
-        else
-        {
-            this->position_ = intersection;
-            this->facei_ = hitFacei;
-        }
-    }
-    else
+    if (hitFacei == -1)
     {
         // Did not find any intersection. Fall back to original approximate
         // algorithm
@@ -129,8 +113,23 @@ Foam::scalar Foam::ExactParticle<ParticleType>::trackToFace
         );
     }
 
+    if (trackFraction >= (1.0-SMALL))
+    {
+        // Nearest intersection beyond endPosition so we hit endPosition.
+        trackFraction = 1.0;
+        this->position_ = endPosition;
+        this->facei_ = -1;
+        return 1.0;
+    }
+    else
+    {
+        this->position_ = intersection;
+        this->facei_ = hitFacei;
+    }
 
-    // Normal situation (trackFraction 0..1)
+
+    // Normal situation (trackFraction 0..1). Straight copy
+    // of Particle::trackToFace.
 
     bool internalFace = this->cloud().internalFace(this->facei_);
 
diff --git a/src/dynamicMesh/motionSmoother/motionSmoother.C b/src/dynamicMesh/motionSmoother/motionSmoother.C
index 18f5c85478..4f40827f4b 100644
--- a/src/dynamicMesh/motionSmoother/motionSmoother.C
+++ b/src/dynamicMesh/motionSmoother/motionSmoother.C
@@ -825,6 +825,28 @@ bool Foam::motionSmoother::scaleMesh
     const bool smoothMesh,
     const label nAllowableErrors
 )
+{
+    return scaleMesh
+    (
+        checkFaces,
+        baffles,
+        paramDict_,
+        paramDict_,
+        smoothMesh,
+        nAllowableErrors
+    );
+}
+
+
+bool Foam::motionSmoother::scaleMesh
+(
+    labelList& checkFaces,
+    const List<labelPair>& baffles,
+    const dictionary& paramDict,
+    const dictionary& meshQualityDict,
+    const bool smoothMesh,
+    const label nAllowableErrors
+)
 {
     if (!smoothMesh && adaptPatchIDs_.empty())
     {
@@ -859,9 +881,9 @@ bool Foam::motionSmoother::scaleMesh
     }
 
     const scalar errorReduction =
-        readScalar(paramDict_.lookup("errorReduction"));
+        readScalar(paramDict.lookup("errorReduction"));
     const label nSmoothScale =
-        readLabel(paramDict_.lookup("nSmoothScale"));
+        readLabel(paramDict.lookup("nSmoothScale"));
 
 
     // Note: displacement_ should already be synced already from setDisplacement
@@ -921,7 +943,7 @@ bool Foam::motionSmoother::scaleMesh
 
     // Check. Returns parallel number of incorrect faces.
     faceSet wrongFaces(mesh_, "wrongFaces", mesh_.nFaces()/100+100);
-    checkMesh(false, mesh_, paramDict_, checkFaces, baffles, wrongFaces);
+    checkMesh(false, mesh_, meshQualityDict, checkFaces, baffles, wrongFaces);
 
     if (returnReduce(wrongFaces.size(), sumOp<label>()) <= nAllowableErrors)
     {
diff --git a/src/dynamicMesh/motionSmoother/motionSmoother.H b/src/dynamicMesh/motionSmoother/motionSmoother.H
index e5777a3664..e4aeac6f50 100644
--- a/src/dynamicMesh/motionSmoother/motionSmoother.H
+++ b/src/dynamicMesh/motionSmoother/motionSmoother.H
@@ -52,7 +52,7 @@ Description
     @endverbatim
 
 Note
-    Shared points (parallel): a processor can have points which are part of
+    - Shared points (parallel): a processor can have points which are part of
     pp on another processor but have no pp itself (i.e. it has points
     and/or edges but no faces of pp). Hence we have to be careful when e.g.
     synchronising displacements that the value from the processor which has
@@ -61,10 +61,13 @@ Note
     else. The combine operator used will give preference to non-zero
     values.
 
-    Various routines take baffles. These are sets of boundary faces that
+    - Various routines take baffles. These are sets of boundary faces that
     are treated as a single internal face. This is a hack used to apply
     movement to internal faces.
 
+    - Mesh constraints are looked up from the supplied dictionary. (uses
+    recursive lookup)
+
 SourceFiles
     motionSmoother.C
     motionSmootherTemplates.C
@@ -395,6 +398,17 @@ public:
                 const label nAllow = 0
             );
 
+            //- Move mesh with externally provided mesh constraints
+            bool scaleMesh
+            (
+                labelList& checkFaces,
+                const List<labelPair>& baffles,
+                const dictionary& paramDict,
+                const dictionary& meshQualityDict,
+                const bool smoothMesh = true,
+                const label nAllow = 0
+            );
+
             //- Update topology
             void updateMesh();
 
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherCheck.C b/src/dynamicMesh/motionSmoother/motionSmootherCheck.C
index 06996c541f..7fcb9f1e7a 100644
--- a/src/dynamicMesh/motionSmoother/motionSmootherCheck.C
+++ b/src/dynamicMesh/motionSmoother/motionSmootherCheck.C
@@ -28,19 +28,6 @@ License
 #include "polyMeshGeometry.H"
 #include "IOmanip.H"
 
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool Foam::motionSmoother::checkMesh
@@ -74,17 +61,50 @@ bool Foam::motionSmoother::checkMesh
     labelHashSet& wrongFaces
 )
 {
-    const scalar maxNonOrtho(readScalar(dict.lookup("maxNonOrtho")));
-    const scalar minVol(readScalar(dict.lookup("minVol")));
-    const scalar maxConcave(readScalar(dict.lookup("maxConcave")));
-    const scalar minArea(readScalar(dict.lookup("minArea")));
-    const scalar maxIntSkew(readScalar(dict.lookup("maxInternalSkewness")));
-    const scalar maxBounSkew(readScalar(dict.lookup("maxBoundarySkewness")));
-    const scalar minWeight(readScalar(dict.lookup("minFaceWeight")));
-    const scalar minVolRatio(readScalar(dict.lookup("minVolRatio")));
-    const scalar minTwist(readScalar(dict.lookup("minTwist")));
-    const scalar minTriangleTwist(readScalar(dict.lookup("minTriangleTwist")));
-    const scalar minDet(readScalar(dict.lookup("minDeterminant")));
+    const scalar maxNonOrtho
+    (
+        readScalar(dict.lookup("maxNonOrtho", true))
+    );
+    const scalar minVol
+    (
+        readScalar(dict.lookup("minVol", true))
+    );
+    const scalar maxConcave
+    (
+        readScalar(dict.lookup("maxConcave", true))
+    );
+    const scalar minArea
+    (
+        readScalar(dict.lookup("minArea", true))
+    );
+    const scalar maxIntSkew
+    (
+        readScalar(dict.lookup("maxInternalSkewness", true))
+    );
+    const scalar maxBounSkew
+    (
+        readScalar(dict.lookup("maxBoundarySkewness", true))
+    );
+    const scalar minWeight
+    (
+        readScalar(dict.lookup("minFaceWeight", true))
+    );
+    const scalar minVolRatio
+    (
+        readScalar(dict.lookup("minVolRatio", true))
+    );
+    const scalar minTwist
+    (
+        readScalar(dict.lookup("minTwist", true))
+    );
+    const scalar minTriangleTwist
+    (
+        readScalar(dict.lookup("minTriangleTwist", true))
+    );
+    const scalar minDet
+    (
+        readScalar(dict.lookup("minDeterminant", true))
+    );
 
     label nWrongFaces = 0;
 
@@ -389,17 +409,50 @@ bool Foam::motionSmoother::checkMesh
     labelHashSet& wrongFaces
 )
 {
-    const scalar maxNonOrtho(readScalar(dict.lookup("maxNonOrtho")));
-    const scalar minVol(readScalar(dict.lookup("minVol")));
-    const scalar maxConcave(readScalar(dict.lookup("maxConcave")));
-    const scalar minArea(readScalar(dict.lookup("minArea")));
-    const scalar maxIntSkew(readScalar(dict.lookup("maxInternalSkewness")));
-    const scalar maxBounSkew(readScalar(dict.lookup("maxBoundarySkewness")));
-    const scalar minWeight(readScalar(dict.lookup("minFaceWeight")));
-    const scalar minVolRatio(readScalar(dict.lookup("minVolRatio")));
-    const scalar minTwist(readScalar(dict.lookup("minTwist")));
-    const scalar minTriangleTwist(readScalar(dict.lookup("minTriangleTwist")));
-    const scalar minDet(readScalar(dict.lookup("minDeterminant")));
+    const scalar maxNonOrtho
+    (
+        readScalar(dict.lookup("maxNonOrtho", true))
+    );
+    const scalar minVol
+    (
+        readScalar(dict.lookup("minVol", true))
+    );
+    const scalar maxConcave
+    (
+        readScalar(dict.lookup("maxConcave", true))
+    );
+    const scalar minArea
+    (
+        readScalar(dict.lookup("minArea", true))
+    );
+    const scalar maxIntSkew
+    (
+        readScalar(dict.lookup("maxInternalSkewness", true))
+    );
+    const scalar maxBounSkew
+    (
+        readScalar(dict.lookup("maxBoundarySkewness", true))
+    );
+    const scalar minWeight
+    (
+        readScalar(dict.lookup("minFaceWeight", true))
+    );
+    const scalar minVolRatio
+    (
+        readScalar(dict.lookup("minVolRatio", true))
+    );
+    const scalar minTwist
+    (
+        readScalar(dict.lookup("minTwist", true))
+    );
+    const scalar minTriangleTwist
+    (
+        readScalar(dict.lookup("minTriangleTwist", true))
+    );
+    const scalar minDet
+    (
+        readScalar(dict.lookup("minDeterminant", true))
+    );
 
     label nWrongFaces = 0;
 
diff --git a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C
index 9dbe7e8c13..61697be9b3 100644
--- a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C
+++ b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C
@@ -31,6 +31,9 @@ License
 #include "mapDistribute.H"
 #include "meshTools.H"
 #include "OFstream.H"
+#include "Random.H"
+#include "treeDataFace.H"
+#include "indexedOctree.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/meshTools/indexedOctree/indexedOctree.C b/src/meshTools/indexedOctree/indexedOctree.C
index 2131c3aaf6..4f55564999 100644
--- a/src/meshTools/indexedOctree/indexedOctree.C
+++ b/src/meshTools/indexedOctree/indexedOctree.C
@@ -29,6 +29,12 @@ License
 #include "meshTools.H"
 #include "OFstream.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template <class Type>
+Foam::scalar Foam::indexedOctree<Type>::perturbTol_ = 10*SMALL;
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 // Does bb intersect a sphere around sample? Or is any corner point of bb
@@ -203,7 +209,7 @@ Foam::indexedOctree<Type>::divide
      || bb.min()[2] >= bb.max()[2]
     )
     {
-        FatalErrorIn("indexedOctree<Type>::divide")
+        FatalErrorIn("indexedOctree<Type>::divide(..)")
             << "Badly formed bounding box:" << bb
             << abort(FatalError);
     }
@@ -674,61 +680,768 @@ void Foam::indexedOctree<Type>::findNearest
 }
 
 
+template <class Type>
+Foam::treeBoundBox Foam::indexedOctree<Type>::subBbox
+(
+    const label parentNodeI,
+    const direction octant
+) const
+{
+    // Get type of node at octant
+    const node& nod = nodes_[parentNodeI];
+    labelBits index = nod.subNodes_[octant];
+
+    if (isNode(index))
+    {
+        // Use stored bb
+        return nodes_[getNode(index)].bb_;
+    }
+    else
+    {
+        // Calculate subBb
+        return nod.bb_.subBbox(octant);
+    }
+}
+
+
+// Takes a bb and a point on/close to the edge of the bb and pushes the point
+// inside by a small fraction. 
+template <class Type>
+Foam::point Foam::indexedOctree<Type>::pushPoint
+(
+    const treeBoundBox& bb,
+    const point& pt,
+    const bool pushInside
+)
+{
+    // Get local length scale.
+    const vector perturbVec = perturbTol_*(bb.span());
+
+    point perturbedPt(pt);
+
+    // Modify all components which are close to any face of the bb to be
+    // well inside/outside them.
+
+    if (pushInside)
+    {
+        for (direction dir = 0; dir < vector::nComponents; dir++)
+        {
+            if (mag(pt[dir]-bb.min()[dir]) < mag(perturbVec[dir]))
+            {
+                // Close to 'left' side. Push well beyond left side.
+                scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
+                perturbedPt[dir] = bb.min()[dir] + perturbDist;
+            }
+            else if (mag(pt[dir]-bb.max()[dir]) < mag(perturbVec[dir]))
+            {
+                // Close to 'right' side. Push well beyond right side.
+                scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
+                perturbedPt[dir] = bb.max()[dir] - perturbDist;
+            }
+        }
+    }
+    else
+    {
+        for (direction dir = 0; dir < vector::nComponents; dir++)
+        {
+            if (mag(pt[dir]-bb.min()[dir]) < mag(perturbVec[dir]))
+            {
+                scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
+                perturbedPt[dir] = bb.min()[dir] - perturbDist;
+            }
+            else if (mag(pt[dir]-bb.max()[dir]) < mag(perturbVec[dir]))
+            {
+                scalar perturbDist = perturbVec[dir] + ROOTVSMALL;
+                perturbedPt[dir] = bb.max()[dir] + perturbDist;
+            }
+        }
+    }
+
+    if (debug)
+    {
+        if (pushInside != bb.contains(perturbedPt))
+        {
+            FatalErrorIn("indexedOctree<Type>::pushPoint(..)")
+                << "pushed point:" << pt
+                << " to:" << perturbedPt
+                << " wanted side:" << pushInside
+                << " obtained side:" << bb.contains(perturbedPt)
+                << " of bb:" << bb
+                << abort(FatalError);
+        }
+    }
+
+    return perturbedPt;
+}
+
+
+// Takes a bb and a point on the edge of the bb and pushes the point
+// outside by a small fraction. 
+template <class Type>
+Foam::point Foam::indexedOctree<Type>::pushPoint
+(
+    const treeBoundBox& bb,
+    const direction faceID,
+    const point& pt,
+    const bool pushInside
+)
+{
+    // Get local length scale.
+    const vector perturbVec = perturbTol_*bb.span();
+
+    point perturbedPt(pt);
+
+    // Modify all components which are close to any face of the bb to be
+    // well outside them.
+
+    if (faceID == 0)
+    {
+        FatalErrorIn("indexedOctree<Type>::pushPoint(..)")
+            << abort(FatalError);
+    }
+
+    if (faceID & treeBoundBox::LEFTBIT)
+    {
+        if (pushInside)
+        {
+            perturbedPt[0] = bb.min()[0] + (perturbVec[0] + ROOTVSMALL);
+        }
+        else
+        {
+            perturbedPt[0] = bb.min()[0] - (perturbVec[0] + ROOTVSMALL);
+        }
+    }
+    else if (faceID & treeBoundBox::RIGHTBIT)
+    {
+        if (pushInside)
+        {
+            perturbedPt[0] = bb.max()[0] - (perturbVec[0] + ROOTVSMALL);
+        }
+        else
+        {
+            perturbedPt[0] = bb.max()[0] + (perturbVec[0] + ROOTVSMALL);
+        }
+    }
+
+    if (faceID & treeBoundBox::BOTTOMBIT)
+    {
+        if (pushInside)
+        {
+            perturbedPt[1] = bb.min()[1] + (perturbVec[1] + ROOTVSMALL);
+        }
+        else
+        {
+            perturbedPt[1] = bb.min()[1] - (perturbVec[1] + ROOTVSMALL);
+        }
+    }
+    else if (faceID & treeBoundBox::TOPBIT)
+    {
+        if (pushInside)
+        {
+            perturbedPt[1] = bb.max()[1] - (perturbVec[1] + ROOTVSMALL);
+        }
+        else
+        {
+            perturbedPt[1] = bb.max()[1] + (perturbVec[1] + ROOTVSMALL);
+        }
+    }
+
+    if (faceID & treeBoundBox::BACKBIT)
+    {
+        if (pushInside)
+        {
+            perturbedPt[2] = bb.min()[2] + (perturbVec[2] + ROOTVSMALL);
+        }
+        else
+        {
+            perturbedPt[2] = bb.min()[2] - (perturbVec[2] + ROOTVSMALL);
+        }
+    }
+    else if (faceID & treeBoundBox::FRONTBIT)
+    {
+        if (pushInside)
+        {
+            perturbedPt[2] = bb.max()[2] - (perturbVec[2] + ROOTVSMALL);
+        }
+        else
+        {
+            perturbedPt[2] = bb.max()[2] + (perturbVec[2] + ROOTVSMALL);
+        }
+    }
+
+    if (debug)
+    {
+        if (pushInside != bb.contains(perturbedPt))
+        {
+            FatalErrorIn("indexedOctree<Type>::pushPoint(..)")
+                << "pushed point:" << pt << " on face:" << faceString(faceID)
+                << " to:" << perturbedPt
+                << " wanted side:" << pushInside
+                << " obtained side:" << bb.contains(perturbedPt)
+                << " of bb:" << bb
+                << abort(FatalError);
+        }
+    }
+
+    return perturbedPt;
+}
+
+
+// Guarantees that if pt is on a face it gets perturbed so it is away
+// from the face edges.
+// If pt is not on a face does nothing.
+template <class Type>
+Foam::point Foam::indexedOctree<Type>::pushPointIntoFace
+(
+    const treeBoundBox& bb,
+    const vector& dir,          // end-start
+    const point& pt
+)
+{
+    if (debug)
+    {
+        if (bb.posBits(pt) != 0)
+        {
+            FatalErrorIn("indexedOctree<Type>::pushPointIntoFace(..)")
+                << " bb:" << bb << endl
+                << "does not contain point " << pt << abort(FatalError);
+        }
+    }
+
+
+    // Handle two cases:
+    // - point exactly on multiple faces. Push away from all but one.
+    // - point on a single face. Push away from edges of face.
+
+    direction ptFaceID = bb.faceBits(pt);
+
+    direction nFaces = 0;
+    FixedList<direction, 3> faceIndices;
+
+    if (ptFaceID & treeBoundBox::LEFTBIT)
+    {
+        faceIndices[nFaces++] = treeBoundBox::LEFT;
+    }
+    else if (ptFaceID & treeBoundBox::RIGHTBIT)
+    {
+        faceIndices[nFaces++] = treeBoundBox::RIGHT;
+    }
+
+    if (ptFaceID & treeBoundBox::BOTTOMBIT)
+    {
+        faceIndices[nFaces++] = treeBoundBox::BOTTOM;
+    }
+    else if (ptFaceID & treeBoundBox::TOPBIT)
+    {
+        faceIndices[nFaces++] = treeBoundBox::TOP;
+    }
+
+    if (ptFaceID & treeBoundBox::BACKBIT)
+    {
+        faceIndices[nFaces++] = treeBoundBox::BACK;
+    }
+    else if (ptFaceID & treeBoundBox::FRONTBIT)
+    {
+        faceIndices[nFaces++] = treeBoundBox::FRONT;
+    }
+
+
+    // Determine the face we want to keep the point on
+
+    direction keepFaceID;
+
+    if (nFaces == 0)
+    {
+        // Return original point
+        return pt;
+    }
+    else if (nFaces == 1)
+    {
+        keepFaceID = faceIndices[0];
+    }
+    else
+    {
+        // Determine best face out of faceIndices[0 .. nFaces-1].
+        // The best face is the one most perpendicular to the ray direction.
+
+        keepFaceID = faceIndices[0];
+        scalar maxInproduct = mag(treeBoundBox::faceNormals[keepFaceID] & dir);
+
+        for (direction i = 1; i < nFaces; i++)
+        {
+            direction face = faceIndices[i];
+            scalar s = mag(treeBoundBox::faceNormals[face] & dir);
+            if (s > maxInproduct)
+            {
+                maxInproduct = s;
+                keepFaceID = face;
+            }
+        }
+    }
+
+
+    // 1. Push point into bb, away from all corners
+
+    point facePoint(pushPoint(bb, pt, true));
+    direction faceID = 0;
+
+    // 2. Snap it back onto the preferred face
+
+    if (keepFaceID == treeBoundBox::LEFT)
+    {
+        facePoint.x() = bb.min().x();
+        faceID = treeBoundBox::LEFTBIT;
+    }
+    else if (keepFaceID == treeBoundBox::RIGHT)
+    {
+        facePoint.x() = bb.max().x();
+        faceID = treeBoundBox::RIGHTBIT;
+    }
+    else if (keepFaceID == treeBoundBox::BOTTOM)
+    {
+        facePoint.y() = bb.min().y();
+        faceID = treeBoundBox::BOTTOMBIT;
+    }
+    else if (keepFaceID == treeBoundBox::TOP)
+    {
+        facePoint.y() = bb.max().y();
+        faceID = treeBoundBox::TOPBIT;
+    }
+    else if (keepFaceID == treeBoundBox::BACK)
+    {
+        facePoint.z() = bb.min().z();
+        faceID = treeBoundBox::BACKBIT;
+    }
+    else if (keepFaceID == treeBoundBox::FRONT)
+    {
+        facePoint.z() = bb.max().z();
+        faceID = treeBoundBox::FRONTBIT;
+    }
+
+
+    if (debug)
+    {
+        if (faceID != bb.faceBits(facePoint))
+        {
+            FatalErrorIn("indexedOctree<Type>::pushPointIntoFace(..)")
+                << "Pushed point from " << pt
+                << " on face:" << ptFaceID << " of bb:" << bb << endl
+                << "onto " << facePoint
+                << " on face:" << faceID
+                << " which is not consistent with geometric face "
+                << bb.faceBits(facePoint)
+                << abort(FatalError);
+        }
+        if (bb.posBits(facePoint) != 0)
+        {
+            FatalErrorIn("indexedOctree<Type>::pushPointIntoFace(..)")
+                << " bb:" << bb << endl
+                << "does not contain perturbed point "
+                << facePoint << abort(FatalError);
+        }
+    }
+
+
+    return facePoint;   
+}
+
+
+//// Takes a bb and a point on the outside of the bb. Checks if on multiple faces
+//// and if so perturbs point so it is only on one face.
+//template <class Type>
+//void Foam::indexedOctree<Type>::checkMultipleFaces
+//(
+//    const treeBoundBox& bb,
+//    const vector& dir,          // end-start
+//    pointIndexHit& faceHitInfo,
+//    direction& faceID
+//)
+//{
+//    // Do the quick elimination of no or one face.
+//    if
+//    (
+//        (faceID == 0)
+//     || (faceID == treeBoundBox::LEFTBIT)
+//     || (faceID == treeBoundBox::RIGHTBIT)
+//     || (faceID == treeBoundBox::BOTTOMBIT)
+//     || (faceID == treeBoundBox::TOPBIT)
+//     || (faceID == treeBoundBox::BACKBIT)
+//     || (faceID == treeBoundBox::FRONTBIT)
+//    )
+//    {
+//        return;
+//    }
+//
+//
+//    // Check the direction of vector w.r.t. faces being intersected.
+//    FixedList<scalar, 6> inproducts(-GREAT);
+//
+//    direction nFaces = 0;
+//
+//    if (faceID & treeBoundBox::LEFTBIT)
+//    {
+//        inproducts[treeBoundBox::LEFT] = mag
+//        (
+//            treeBoundBox::faceNormals[treeBoundBox::LEFT]
+//          & dir
+//        );
+//        nFaces++;
+//    }
+//    if (faceID & treeBoundBox::RIGHTBIT)
+//    {
+//        inproducts[treeBoundBox::RIGHT] = mag
+//        (
+//            treeBoundBox::faceNormals[treeBoundBox::RIGHT]
+//          & dir
+//        );
+//        nFaces++;
+//    }
+//
+//    if (faceID & treeBoundBox::BOTTOMBIT)
+//    {
+//        inproducts[treeBoundBox::BOTTOM] = mag
+//        (
+//            treeBoundBox::faceNormals[treeBoundBox::BOTTOM]
+//          & dir
+//        );
+//        nFaces++;
+//    }
+//    if (faceID & treeBoundBox::TOPBIT)
+//    {
+//        inproducts[treeBoundBox::TOP] = mag
+//        (
+//            treeBoundBox::faceNormals[treeBoundBox::TOP]
+//          & dir
+//        );
+//        nFaces++;
+//    }
+//
+//    if (faceID & treeBoundBox::BACKBIT)
+//    {
+//        inproducts[treeBoundBox::BACK] = mag
+//        (
+//            treeBoundBox::faceNormals[treeBoundBox::BACK]
+//          & dir
+//        );
+//        nFaces++;
+//    }
+//    if (faceID & treeBoundBox::FRONTBIT)
+//    {
+//        inproducts[treeBoundBox::FRONT] = mag
+//        (
+//            treeBoundBox::faceNormals[treeBoundBox::FRONT]
+//          & dir
+//        );
+//        nFaces++;
+//    }
+//
+//    if (nFaces == 0 || nFaces == 1 || nFaces > 3)
+//    {
+//        FatalErrorIn("indexedOctree<Type>::checkMultipleFaces(..)")
+//            << "Problem : nFaces:" << nFaces << abort(FatalError);
+//    }
+//
+//    // Keep point on most perpendicular face; shift it away from the aligned
+//    // ones.
+//    // E.g. line hits top and left face:
+//    //     a
+//    // ----+----+
+//    //     |    |
+//    //     |    |
+//    //     +----+
+//    // Shift point down (away from top):
+//    //     
+//    //    a+----+
+//    // ----|    |
+//    //     |    |
+//    //     +----+
+//
+//    label maxIndex = -1;
+//    scalar maxInproduct = -GREAT;
+//
+//    for (direction i = 0; i < 6; i++)
+//    {
+//        if (inproducts[i] > maxInproduct)
+//        {
+//            maxInproduct = inproducts[i];
+//            maxIndex = i;
+//        }
+//    }
+//
+//    if (maxIndex == -1)
+//    {
+//        FatalErrorIn("indexedOctree<Type>::checkMultipleFaces(..)")
+//            << "Problem maxIndex:" << maxIndex << " inproducts:" << inproducts
+//            << abort(FatalError);
+//    }
+//
+//    const point oldPoint(faceHitInfo.rawPoint());
+//    const direction oldFaceID = faceID;
+//
+//    // 1. Push point into bb, away from all corners
+//
+//    faceHitInfo.rawPoint() = pushPoint(bb, oldFaceID, oldPoint, true);
+//
+//    // 2. Snap it back onto the preferred face
+//
+//    if (maxIndex == treeBoundBox::LEFT)
+//    {
+//        faceHitInfo.rawPoint().x() = bb.min().x();
+//        faceID = treeBoundBox::LEFTBIT;
+//    }
+//    else if (maxIndex == treeBoundBox::RIGHT)
+//    {
+//        faceHitInfo.rawPoint().x() = bb.max().x();
+//        faceID = treeBoundBox::RIGHTBIT;
+//    }
+//    else if (maxIndex == treeBoundBox::BOTTOM)
+//    {
+//        faceHitInfo.rawPoint().y() = bb.min().y();
+//        faceID = treeBoundBox::BOTTOMBIT;
+//    }
+//    else if (maxIndex == treeBoundBox::TOP)
+//    {
+//        faceHitInfo.rawPoint().y() = bb.max().y();
+//        faceID = treeBoundBox::TOPBIT;
+//    }
+//    else if (maxIndex == treeBoundBox::BACK)
+//    {
+//        faceHitInfo.rawPoint().z() = bb.min().z();
+//        faceID = treeBoundBox::BACKBIT;
+//    }
+//    else if (maxIndex == treeBoundBox::FRONT)
+//    {
+//        faceHitInfo.rawPoint().z() = bb.max().z();
+//        faceID = treeBoundBox::FRONTBIT;
+//    }
+//
+//    Pout<< "From ray:" << dir
+//        << " from point:" << oldPoint
+//        << " on faces:" << faceString(oldFaceID)
+//        << " of bb:" << bb
+//        << " with inprods:" << inproducts
+//        << " maxIndex:" << maxIndex << endl
+//        << "perturbed to point:" << faceHitInfo.rawPoint()
+//        << " on face:" << faceString(faceID)
+//        << endl;
+//
+//
+//    if (debug)
+//    {
+//        if (faceID != bb.faceBits(faceHitInfo.rawPoint()))
+//        {
+//            FatalErrorIn("indexedOctree<Type>::checkMultipleFaces(..)")
+//                << "Pushed point from " << oldPoint
+//                << " on face:" << oldFaceID << " of bb:" << bb << endl
+//                << "onto " << faceHitInfo.rawPoint()
+//                << " on face:" << faceID
+//                << " which is not consistent with geometric face "
+//                <<  bb.faceBits(faceHitInfo.rawPoint())
+//                << abort(FatalError);
+//        }
+//    }
+//}
+
+
+// Get parent node and octant. Return false if top of tree reached.
+template <class Type>
+bool Foam::indexedOctree<Type>::walkToParent
+(
+    const label nodeI,
+    const direction octant,
+
+    label& parentNodeI,
+    label& parentOctant
+) const
+{
+    parentNodeI = nodes_[nodeI].parent_;
+
+    if (parentNodeI == -1)
+    {
+        // Reached edge of tree
+        return false;
+    }
+
+    const node& parentNode = nodes_[parentNodeI];
+
+    // Find octant nodeI is in.
+    parentOctant = 255;
+
+    for (direction i = 0; i < parentNode.subNodes_.size(); i++)
+    {
+        labelBits index = parentNode.subNodes_[i];
+
+        if (isNode(index) && getNode(index) == nodeI)
+        {
+            parentOctant = i;
+            break;
+        }
+    }
+
+    if (parentOctant == 255)
+    {
+        FatalErrorIn("walkToParent(..)")
+            << "Problem: no parent found for octant:" << octant
+            << " node:" << nodeI
+            << abort(FatalError);
+    }
+
+    return true;
+}
+
+
 // Walk tree to neighbouring node. Gets current position as
 // node and octant in this node and walks in the direction given by
-// the faceID (one of treeBoundBox::LEFTBIT, RIGHTBIT etc.)
+// the facePointBits (combination of treeBoundBox::LEFTBIT, TOPBIT etc.)
 // Returns false if edge of tree hit.
 template <class Type>
 bool Foam::indexedOctree<Type>::walkToNeighbour
 (
     const point& facePoint,
-    const direction faceID,         // direction to walk in
+    const direction faceID,  // face(s) that facePoint is on
     label& nodeI,
     direction& octant
 ) const
 {
+    label oldNodeI = nodeI;
+    direction oldOctant = octant;
+
     // Find out how to test for octant. Say if we want to go left we need
     // to traverse up the tree until we hit a node where our octant is
     // on the right.
 
+    // Coordinate direction to test
+    const direction X = treeBoundBox::RIGHTHALF;
+    const direction Y = treeBoundBox::TOPHALF;
+    const direction Z = treeBoundBox::FRONTHALF;
+
     direction octantMask = 0;
-    direction valueMask = 0;
+    direction wantedValue = 0;
 
     if ((faceID & treeBoundBox::LEFTBIT) != 0)
     {
-        // We want to go left so check if in right octant.
-        octantMask |= treeBoundBox::RIGHTHALF;
-        valueMask |= treeBoundBox::RIGHTHALF;
+        // We want to go left so check if in right octant (i.e. x-bit is set)
+        octantMask |= X;
+        wantedValue |= X;
     }
     else if ((faceID & treeBoundBox::RIGHTBIT) != 0)
     {
-        octantMask |= treeBoundBox::RIGHTHALF;  // valueMask already 0
+        octantMask |= X;  // wantedValue already 0
     }
 
     if ((faceID & treeBoundBox::BOTTOMBIT) != 0)
     {
-        octantMask |= treeBoundBox::TOPHALF;
-        valueMask |= treeBoundBox::TOPHALF;
+        // Want to go down so check for y-bit set.
+        octantMask |= Y;
+        wantedValue |= Y;
     }
     else if ((faceID & treeBoundBox::TOPBIT) != 0)
     {
-        octantMask |= treeBoundBox::TOPHALF;
+        // Want to go up so check for y-bit not set.
+        octantMask |= Y;
     }
 
     if ((faceID & treeBoundBox::BACKBIT) != 0)
     {
-        octantMask |= treeBoundBox::FRONTHALF;
-        valueMask |= treeBoundBox::FRONTHALF;
+        octantMask |= Z;
+        wantedValue |= Z;
     }
     else if ((faceID & treeBoundBox::FRONTBIT) != 0)
     {
-        octantMask |= treeBoundBox::FRONTHALF;
+        octantMask |= Z;
     }
 
+    // So now we have the coordinate directions in the octant we need to check
+    // and the resulting value.
+
+    /*
+    // +---+---+
+    // |   |   |
+    // |   |   |
+    // |   |   |
+    // +---+-+-+
+    // |   | | |
+    // |  a+-+-+
+    // |   |\| |
+    // +---+-+-+
+    //        \ 
+    //
+    // e.g. ray is at (a) in octant 0(or 4) with faceIDs : LEFTBIT+TOPBIT.
+    // If we would be in octant 1(or 5) we could go to the correct octant
+    // in the same node by just flipping the x and y bits (exoring).
+    // But if we are not in octant 1/5 we have to go up until we are.
+    // In general for leftbit+topbit:
+    // - we have to check for x and y : octantMask  = 011
+    // - the checked bits have to be  : wantedValue = ?01
+    */
+
+    //Pout<< "For point " << facePoint << endl;
+
     // Go up until we have chance to cross to the wanted direction
-    while (valueMask != (octant & octantMask))
+    while (wantedValue != (octant & octantMask))
     {
-        label parentNodeI = nodes_[nodeI].parent_;
+        // Go up to the parent.
+
+        // Remove the directions that are not on the boundary of the parent.
+        // See diagram above
+        if (wantedValue & X)    // && octantMask&X
+        {
+            // Looking for right octant.
+            if (octant & X)
+            {
+                // My octant is one of the left ones so punch out the x check
+                octantMask &= ~X;
+                wantedValue &= ~X;
+            }
+        }
+        else
+        {
+            if (!(octant & X))
+            {
+                // My octant is right but I want to go left.
+                octantMask &= ~X;
+                wantedValue &= ~X;
+            }
+        }
+
+        if (wantedValue & Y)
+        {
+            if (octant & Y)
+            {
+                octantMask &= ~Y;
+                wantedValue &= ~Y;
+            }
+        }
+        else
+        {
+            if (!(octant & Y))
+            {
+                octantMask &= ~Y;
+                wantedValue &= ~Y;
+            }
+        }
+
+        if (wantedValue & Z)
+        {
+            if (octant & Z)
+            {
+                octantMask &= ~Z;
+                wantedValue &= ~Z;
+            }
+        }
+        else
+        {
+            if (!(octant & Z))
+            {
+                octantMask &= ~Z;
+                wantedValue &= ~Z;
+            }
+        }
+
+
+        label parentNodeI;
+        label parentOctant;
+        walkToParent(nodeI, octant, parentNodeI, parentOctant);
 
         if (parentNodeI == -1)
         {
@@ -736,38 +1449,46 @@ bool Foam::indexedOctree<Type>::walkToNeighbour
             return false;
         }
 
-        const node& parentNode = nodes_[parentNodeI];
+        //Pout<< "    walked from node:" << nodeI << " octant:" << octant
+        //    << " bb:" << nodes_[nodeI].bb_.subBbox(octant) << endl
+        //    << "    to:" << parentNodeI << " octant:" << parentOctant
+        //    << " bb:" << nodes_[parentNodeI].bb_.subBbox(parentOctant)
+        //    << endl;
+        //
+        //Pout<< "    octantMask:" << octantMask
+        //    << " wantedValue:" << wantedValue << endl;
 
-        // Find octant nodeI is in.
-        direction parentOctant = 255;
-
-        for (direction i = 0; i < parentNode.subNodes_.size(); i++)
-        {
-            labelBits index = parentNode.subNodes_[i];
-
-            if (isNode(index) && getNode(index) == nodeI)
-            {
-                parentOctant = i;
-                break;
-            }
-        }
-
-        if (parentOctant == 255)
-        {
-            FatalErrorIn("walkToNeighbour")
-                << abort(FatalError);
-        }
         nodeI = parentNodeI;
         octant = parentOctant;
     }
 
-    // So now we hit the correct parent node. Switch to the correct octant
+    // So now we hit the correct parent node. Switch to the correct octant.
+    // Is done by jumping to the 'other half' so if e.g. in x direction in
+    // right half we now jump to the left half.
     octant ^= octantMask;
 
+    //Pout<< "    to node:" << nodeI << " octant:" << octant
+    //    << " subBb:" <<subBbox(nodeI, octant) << endl;
+
+
+    if (debug)
+    {
+        const treeBoundBox subBb(subBbox(nodeI, octant));
+
+        if (!subBb.contains(facePoint))
+        {
+            FatalErrorIn("indexedOctree<Type>::walkToNeighbour(..)")
+                << "When searching for " << facePoint
+                << " ended up in node:" << nodeI
+                << " octant:" << octant
+                << " with bb:" << subBb
+                << abort(FatalError);
+        }
+    }
+
+
     // See if we need to travel down. Note that we already go into the
-    // the first level ourselves (instead of having findNode decide) since
-    // the facePoint is now exactly on the mid of the node so there could
-    // be truncation problems.
+    // the first level ourselves (instead of having findNode decide)
     labelBits index = nodes_[nodeI].subNodes_[octant];
 
     if (isNode(index))
@@ -778,49 +1499,82 @@ bool Foam::indexedOctree<Type>::walkToNeighbour
         octant = getOctant(node);
     }
 
+
+    if (debug)
+    {
+        const treeBoundBox subBb(subBbox(nodeI, octant));
+
+        if (nodeI == oldNodeI && octant == oldOctant)
+        {
+            FatalErrorIn("indexedOctree<Type>::walkToNeighbour(..)")
+                << "Did not go to neighbour when searching for " << facePoint
+                << endl
+                << "    starting from face:" << faceString(faceID)
+                << " node:" << nodeI
+                << " octant:" << octant
+                << " bb:" << subBb
+                << abort(FatalError);
+        }
+
+        if (!subBb.contains(facePoint))
+        {
+            FatalErrorIn("indexedOctree<Type>::walkToNeighbour(..)")
+                << "When searching for " << facePoint
+                << " ended up in node:" << nodeI
+                << " octant:" << octant
+                << " bb:" << subBb
+                << abort(FatalError);
+        }
+    }
+
+
     return true;
 }
 
 
-// Return unique number for the face of a bounding box a point is on.
-// (number is single bit but not really nessecary)
-// Return 0 if point not on any face of bb.
 template <class Type>
-Foam::direction Foam::indexedOctree<Type>::getFace
+Foam::word Foam::indexedOctree<Type>::faceString
 (
-    const treeBoundBox& bb,
-    const point& pt
+    const direction faceID
 )
 {
-    direction faceID = 0;
+    word desc;
 
-    if (pt.x() <= bb.min().x())
+    if (faceID == 0)
     {
-        faceID |= treeBoundBox::LEFTBIT;
+        desc = "noFace";
     }
-    if (pt.x() >= bb.max().x())
+    if (faceID & treeBoundBox::LEFTBIT)
     {
-        faceID |= treeBoundBox::RIGHTBIT;
+        if (!desc.empty()) desc += "+";
+        desc += "left";
     }
-
-    if (pt.y() <= bb.min().y())
+    if (faceID & treeBoundBox::RIGHTBIT)
     {
-        faceID |= treeBoundBox::BOTTOMBIT;
+        if (!desc.empty()) desc += "+";
+        desc += "right";
     }
-    if (pt.y() >= bb.max().y())
+    if (faceID & treeBoundBox::BOTTOMBIT)
     {
-        faceID |= treeBoundBox::TOPBIT;
+        if (!desc.empty()) desc += "+";
+        desc += "bottom";
     }
-
-    if (pt.z() <= bb.min().z())
+    if (faceID & treeBoundBox::TOPBIT)
     {
-        faceID |= treeBoundBox::BACKBIT;
+        if (!desc.empty()) desc += "+";
+        desc += "top";
     }
-    if (pt.z() >= bb.max().z())
+    if (faceID & treeBoundBox::BACKBIT)
     {
-        faceID |= treeBoundBox::FRONTBIT;
+        if (!desc.empty()) desc += "+";
+        desc += "back";
     }
-    return faceID;
+    if (faceID & treeBoundBox::FRONTBIT)
+    {
+        if (!desc.empty()) desc += "+";
+        desc += "front";
+    }
+    return desc;
 }
 
 
@@ -829,21 +1583,37 @@ Foam::direction Foam::indexedOctree<Type>::getFace
 //  hitInfo.point = point on shape
 // Else return a miss and the bounding box face hit:
 //  hitInfo.point = coordinate of intersection of ray with bounding box
-//  faceID  = index of bounding box face
+//  hitBits  = posbits of point on bounding box
 template <class Type>
 void Foam::indexedOctree<Type>::traverseNode
 (
     const bool findAny,
+    const point& treeStart,
+    const vector& treeVec,
+
     const point& start,
     const point& end,
-    const vector& dir,
     const label nodeI,
     const direction octant,
 
     pointIndexHit& hitInfo,
-    direction& faceID
+    direction& hitBits
 ) const
 {
+    if (debug)
+    {
+        const treeBoundBox octantBb(subBbox(nodeI, octant));
+
+        if (octantBb.posBits(start) != 0)
+        {
+            FatalErrorIn("indexedOctree<Type>::traverseNode(..)")
+                << "Node:" << nodeI << " octant:" << octant
+                << " bb:" << octantBb << endl
+                << "does not contain point " << start << abort(FatalError);
+        }
+    }
+
+
     const node& nod = nodes_[nodeI];
 
     labelBits index = nod.subNodes_[octant];
@@ -909,49 +1679,72 @@ void Foam::indexedOctree<Type>::traverseNode
     // Nothing intersected in this node
     // Traverse to end of node. Do by ray tracing back from end and finding
     // intersection with bounding box of node.
+    // start point is now considered to be inside bounding box.
+
+    const treeBoundBox octantBb(subBbox(nodeI, octant));
 
     point pt;
+    bool intersected = octantBb.intersects
+    (
+        end,            //treeStart,
+        (start-end),    //treeVec,
 
-    if (isNode(index))
+        end,
+        start,
+
+        pt,
+        hitBits
+    );
+
+
+    if (intersected)
     {
-        const treeBoundBox& subBb = nodes_[getNode(index)].bb_;
-
-        if (!subBb.intersects(end, start, pt))
-        {
-            faceID = 0;
-
-            FatalErrorIn("indexedOctree<Type>::traverseNode(..)")
-                << "Did not hit side of box " << subBb
-                << " with ray from " << start << " to " << end
-                << abort(FatalError);
-        }
-        else
-        {
-            faceID = getFace(subBb, pt);
-        }
+        // Return miss. Set misspoint to face.
+        hitInfo.setPoint(pt);
     }
     else
     {
-        const treeBoundBox subBb(nod.bb_.subBbox(octant));
+        // Rare case: did not find intersection of ray with octantBb. Can
+        // happen if end is on face/edge of octantBb. Do like
+        // lagrangian and re-track with perturbed vector (slightly pushed out
+        // of bounding box)
 
-        if (!subBb.intersects(end, start, pt))
-        {
-            faceID = 0;
+        point perturbedEnd(pushPoint(octantBb, end, false));
 
-            FatalErrorIn("indexedOctree<Type>::traverseNode(..)")
-                << "Did not hit side of box " << subBb
-                << " with ray from " << start << " to " << end
-                << abort(FatalError);
-        }
-        else
+
+        //if (debug)
         {
-            faceID = getFace(subBb, pt);
+            // Dump octantBb to obj
+            writeOBJ(nodeI, octant);
+            // Dump ray to obj as well
+            {
+                OFstream str("ray.obj");
+                meshTools::writeOBJ(str, start);
+                meshTools::writeOBJ(str, end);
+                str << "l 1 2" << nl;
+            }
+            WarningIn("indexedOctree<Type>::traverseNode(..)")
+                << "Did not intersect ray from endpoint:" << end
+                << " to startpoint:" << start
+                << " with bounding box:" << octantBb << nl
+                << "Re-intersecting with perturbed endpoint:" << perturbedEnd
+                << endl;
         }
+
+        traverseNode
+        (
+            findAny,
+            treeStart,
+            treeVec,
+            start,
+            perturbedEnd,
+            nodeI,
+            octant,
+
+            hitInfo,
+            hitBits
+        );
     }
-
-
-    // Return miss. Set misspoint to face.
-    hitInfo.setPoint(pt);
 }
 
 
@@ -963,40 +1756,75 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
     const point& treeStart,
     const point& treeEnd,
     const label startNodeI,
-    const direction startOctant
+    const direction startOctant,
+    const bool verbose
 ) const
 {
-    const vector dir(treeEnd - treeStart);
+    const vector treeVec(treeEnd - treeStart);
 
     // Current node as parent+octant
     label nodeI = startNodeI;
     direction octant = startOctant;
 
+    if (verbose)
+    {
+        Pout<< "findLine : treeStart:" << treeStart
+            << " treeEnd:" << treeEnd << endl
+            << "node:" << nodeI
+            << " octant:" << octant
+            << " bb:" << subBbox(nodeI, octant) << endl;
+    }
+
     // Current position. Initialize to miss
     pointIndexHit hitInfo(false, treeStart, -1);
 
-    // Side of current bounding box current point is on
-    direction faceID = 0;
-
-    while (true)
+    //while (true)
+    label i = 0;
+    for (; i < 100000; i++)
     {
+        if (verbose)
+        {
+            Pout<< "iter:" << i
+                << " at startPoint:" << hitInfo.rawPoint() << endl
+                << "    node:" << nodeI
+                << " octant:" << octant
+                << " bb:" << subBbox(nodeI, octant) << endl;
+        }
+
+
         // Ray-trace to end of current node. Updates point (either on triangle
         // in case of hit or on node bounding box in case of miss)
 
-        point startPoint(hitInfo.rawPoint());
+        const treeBoundBox octantBb(subBbox(nodeI, octant));
+
+        // Make sure point is away from any edges/corners
+        point startPoint
+        (
+            pushPointIntoFace
+            (
+                octantBb,
+                treeVec,
+                hitInfo.rawPoint()
+            )
+        );
+
+        // Faces of current bounding box current point is on
+        direction hitFaceID = 0;
 
         traverseNode
         (
             findAny,
+            treeStart,
+            treeVec,
+
             startPoint,     // Note: pass in copy since hitInfo
                             // also used in return value.
-            treeEnd,
-            dir,
+            treeEnd,        // pass in overall end so is nicely outside bb
             nodeI,
             octant,
 
             hitInfo,
-            faceID
+            hitFaceID
         );
 
         // Did we hit a triangle?
@@ -1005,29 +1833,48 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
             break;
         }
 
-        if (faceID == 0)
+        if (hitFaceID == 0)
         {
-            // Was never inside the tree. Return miss.
+            // endpoint inside the tree. Return miss.
             break;
         }
 
-        //Pout<< "findLine : treeStart:" << treeStart
-        //    << " treeEnd:" << treeEnd
-        //    << " tracked from " << startPoint << " to edge of nodeBb:"
-        //    << hitInfo.missPoint()
-        //    << " node:" << nodeI << " octant:" << octant
-        //    << " subBb:"
-        //    << nodes_[nodeI].bb_.subBbox(octant)
-        //    << endl;
+
+        if (verbose)
+        {
+            Pout<< "    iter:" << i
+                << " hit face:" << faceString(hitFaceID)
+                << " at:" << hitInfo.rawPoint() << nl
+                << "    node:" << nodeI
+                << " octant:" << octant
+                << " bb:" << subBbox(nodeI, octant) << nl
+                << "    walking to neighbour containing:"
+                <<  pushPoint
+                    (
+                        octantBb,
+                        hitFaceID,
+                        hitInfo.rawPoint(),
+                        false                   // push outside of octantBb
+                    )
+                << endl;
+        }
 
 
         // Nothing hit so we are on face of bounding box (given as node+octant+
-        // faceID). Traverse to neighbouring node.
+        // position bits). Traverse to neighbouring node. Use slightly perturbed
+        // point.
 
         bool ok = walkToNeighbour
         (
-            hitInfo.rawPoint(), // point on face
-            faceID,             // direction to walk in
+            pushPoint
+            (
+                octantBb,
+                hitFaceID,
+                hitInfo.rawPoint(),
+                false                   // push outside of octantBb
+            ),
+            hitFaceID,  // face(s) that hitInfo is on
+
             nodeI,
             octant
         );
@@ -1037,7 +1884,53 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
             // Hit the edge of the tree. Return miss.
             break;
         }
+
+        if (verbose)
+        {
+            const treeBoundBox octantBb(subBbox(nodeI, octant));
+            Pout<< "    walked for point:" << hitInfo.rawPoint() << endl
+                << "    to neighbour node:" << nodeI
+                << " octant:" << octant
+                << " face:" << faceString(octantBb.faceBits(hitInfo.rawPoint()))
+                << " of octantBb:" << octantBb << endl
+                << endl;
+        }
     }
+
+    if (i == 100000)
+    {
+        // Probably in loop.
+        if (!verbose)
+        {
+            // Redo intersection but now with verbose flag switched on.
+            return findLine
+            (
+                findAny,
+                treeStart,
+                treeEnd,
+                startNodeI,
+                startOctant,
+                true            //verbose
+            );
+        }
+        if (debug)
+        {
+            FatalErrorIn("indexedOctree<Type>::findLine(..)")
+                << "Got stuck in loop raytracing from:" << treeStart
+                << " to:" << treeEnd << endl
+                << "inside top box:" << subBbox(startNodeI, startOctant)
+                << abort(FatalError);
+        }
+        else
+        {
+            WarningIn("indexedOctree<Type>::findLine(..)")
+                << "Got stuck in loop raytracing from:" << treeStart
+                << " to:" << treeEnd << endl
+                << "inside top box:" << subBbox(startNodeI, startOctant)
+                << endl;
+        }
+    }
+
     return hitInfo;
 }
 
@@ -1057,6 +1950,9 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
     {
         const treeBoundBox& treeBb = nodes_[0].bb_;
 
+        // No effort is made to deal with points which are on edge of tree
+        // bounding box for now.
+
         direction startBit = treeBb.posBits(start);
         direction endBit = treeBb.posBits(end);
 
@@ -1066,6 +1962,9 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
             return pointIndexHit(false, vector::zero, -1);
         }
 
+
+        // Trim segment to treeBb
+
         point trackStart(start);
         point trackEnd(end);
 
@@ -1087,6 +1986,7 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
             }
         }
 
+
         // Find lowest level tree node that start is in.
         labelBits index = findNode(0, trackStart);
 
@@ -1234,32 +2134,6 @@ void Foam::indexedOctree<Type>::writeOBJ
 
         str<< "l " << e[0]+pointVertI+1 << ' ' << e[1]+pointVertI+1 << nl;
     }
-
-
-    //// Dump triangles
-    //if (isContent(index))
-    //{
-    //    const labelList& indices = contents_[getContent(index)];
-    //    const triSurface& surf = shapes_.surface();
-    //    const pointField& points = surf.points();
-    //
-    //    forAll(indices, i)
-    //    {
-    //        label shapeI = indices[i];
-    //
-    //        const labelledTri& f = surf[shapeI];
-    //
-    //        meshTools::writeOBJ(str, points[f[0]]);
-    //        vertI++;
-    //        meshTools::writeOBJ(str, points[f[1]]);
-    //        vertI++;
-    //        meshTools::writeOBJ(str, points[f[2]]);
-    //        vertI++;
-    //
-    //        str<< "l " << vertI-2 << ' ' << vertI-1 << ' ' << vertI << ' '
-    //            << vertI-2 << nl;
-    //    }
-    //}
 }
 
 
@@ -1305,6 +2179,14 @@ Foam::indexedOctree<Type>::indexedOctree
     contents_(0),
     nodeTypes_(0)
 {
+    if (debug)
+    {
+        Pout<< "indexedOctree<Type>::indexedOctree:" << nl
+            << "    shapes:" << shapes.size() << nl
+            << "    bb:" << bb << nl
+            << endl;
+    }
+
     if (shapes.size() == 0)
     {
         return;
@@ -1409,7 +2291,6 @@ Foam::indexedOctree<Type>::indexedOctree
     nodes_.transfer(nodes);
     nodes.clear();
 
-
     if (debug)
     {
         label nEntries = 0;
@@ -1418,14 +2299,18 @@ Foam::indexedOctree<Type>::indexedOctree
             nEntries += contents_[i].size();
         }
 
-        Pout<< "indexedOctree<Type>::indexedOctree : finished construction:"
+        Pout<< "indexedOctree<Type>::indexedOctree"
+            << " : finished construction of tree of:" << shapes.typeName
             << nl
+            << "    bb:" << this->bb() << nl
             << "    shapes:" << shapes.size() << nl
             << "    nLevels:" << nLevels << nl
             << "    treeNodes:" << nodes_.size() << nl
             << "    nEntries:" << nEntries << nl
-            << "        per treeLeaf:" << nEntries/contents.size() << nl
-            << "        per shape (duplicity):" << nEntries/shapes.size() << nl
+            << "        per treeLeaf:"
+            << scalar(nEntries)/contents.size() << nl
+            << "        per shape (duplicity):"
+            << scalar(nEntries)/shapes.size() << nl
             << endl;
     }
 }
@@ -1447,6 +2332,13 @@ Foam::indexedOctree<Type>::indexedOctree
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template <class Type>
+Foam::scalar& Foam::indexedOctree<Type>::perturbTol()
+{
+    return perturbTol_;
+}
+
+
 template <class Type>
 Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
 (
@@ -1570,6 +2462,16 @@ Foam::labelBits Foam::indexedOctree<Type>::findNode
 
     const node& nod = nodes_[nodeI];
 
+    if (debug)
+    {
+        if (!nod.bb_.contains(sample))
+        {
+            FatalErrorIn("findNode(..)")
+                << "Cannot find " << sample << " in node " << nodeI
+                << abort(FatalError);
+        }
+    }
+
     direction octant = nod.bb_.subOctant(sample);
 
     labelBits index = nod.subNodes_[octant];
diff --git a/src/meshTools/indexedOctree/indexedOctree.H b/src/meshTools/indexedOctree/indexedOctree.H
index b4aad9a305..6c2bfa3c6e 100644
--- a/src/meshTools/indexedOctree/indexedOctree.H
+++ b/src/meshTools/indexedOctree/indexedOctree.H
@@ -128,6 +128,15 @@ public:
 
 private:
 
+    // Static data
+
+        //- Relative peturbation tolerance. Determines when point is
+        //  considered to be close to face/edge of bb of node.
+        //  The tolerance is relative to the bounding box of the smallest
+        //  node.
+        static scalar perturbTol_;
+
+
     // Private data
 
         //- Underlying shapes for geometric queries.
@@ -144,8 +153,9 @@ private:
 
     // Private Member Functions
 
-        //- Like above but now bb is implicitly provided as parent bb + mid
-        //  + octant
+        //- Helper: does bb intersect a sphere around sample? Or is any
+        //  corner point of bb closer than nearestDistSqr to sample.
+        //  (bb is implicitly provided as parent bb + octant)
         static bool overlaps
         (
             const treeBoundBox& parentBb,
@@ -215,6 +225,60 @@ private:
                 point& nearestPoint
             ) const;
 
+            //- Return bbox of octant
+            treeBoundBox subBbox
+            (
+                const label parentNodeI,
+                const direction octant
+            ) const;
+
+            //- Helper: take a point on/close to face of bb and push it
+            //  inside or outside of bb.
+            static point pushPoint
+            (
+                const treeBoundBox&,
+                const point&,
+                const bool pushInside
+            );
+
+            //- Helper: take a point on face of bb and push it
+            //  inside or outside of bb.
+            static point pushPoint
+            (
+                const treeBoundBox&,
+                const direction,
+                const point&,
+                const bool pushInside
+            );
+
+            //- Helper: take point on face(s) of bb and push it away from
+            //  edges of face.
+            static point pushPointIntoFace
+            (
+                const treeBoundBox& bb,
+                const vector& dir,          // end-start
+                const point& pt
+            );
+
+            ////- Push point on multiple faces away from any corner/edge.
+            //static void checkMultipleFaces
+            //(
+            //    const treeBoundBox& bb,
+            //    const vector& dir,          // end-start
+            //    pointIndexHit& faceHitInfo,
+            //    direction& faceID
+            //);
+
+            //- Walk to parent of node+octant.
+            bool walkToParent
+            (
+                const label nodeI,
+                const direction octant,
+
+                label& parentNodeI,
+                label& parentOctant
+            ) const;
+
             //- Walk tree to neighbouring node. Return false if edge of tree
             //  hit.
             bool walkToNeighbour
@@ -225,8 +289,8 @@ private:
                 direction& octant
             ) const;
 
-            //- Categorize point on bounding box.
-            static direction getFace(const treeBoundBox&, const point&);
+            //- Debug: return verbose the bounding box faces
+            static word faceString(const direction faceID);
 
             //- Traverse a node. If intersects a triangle return first
             // intersection point.
@@ -235,9 +299,11 @@ private:
             void traverseNode
             (
                 const bool findAny,
+                const point& treeStart,
+                const vector& treeVec,
+
                 const point& start,
                 const point& end,
-                const vector& dir,
                 const label nodeI,
                 const direction octantI,
 
@@ -252,7 +318,8 @@ private:
                 const point& treeStart,
                 const point& treeEnd,
                 const label startNodeI,
-                const direction startOctantI
+                const direction startOctantI,
+                const bool verbose = false
             ) const;
 
             //- Find any or nearest intersection of line between start and end.
@@ -310,6 +377,10 @@ private:
 
 public:
 
+        //- Get the perturbation tolerance
+        static scalar& perturbTol();
+
+
     // Constructors
 
         //- Construct null
@@ -505,6 +576,7 @@ public:
                 const point& sample
             );
 
+
         // Write
 
             //- Print tree. Either print all indices (printContent = true) or
diff --git a/src/meshTools/indexedOctree/treeDataCell.H b/src/meshTools/indexedOctree/treeDataCell.H
index 4f09188942..100dcb62cf 100644
--- a/src/meshTools/indexedOctree/treeDataCell.H
+++ b/src/meshTools/indexedOctree/treeDataCell.H
@@ -37,9 +37,7 @@ SourceFiles
 #ifndef treeDataCell_H
 #define treeDataCell_H
 
-#include "treeBoundBox.H"
 #include "treeBoundBoxList.H"
-#include "labelList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/meshTools/indexedOctree/treeDataEdge.C b/src/meshTools/indexedOctree/treeDataEdge.C
index c4bb08f374..12f347c164 100644
--- a/src/meshTools/indexedOctree/treeDataEdge.C
+++ b/src/meshTools/indexedOctree/treeDataEdge.C
@@ -22,16 +22,15 @@ 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 "treeDataEdge.H"
 #include "indexedOctree.H"
-#include "polyMesh.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
+defineTypeNameAndDebug(Foam::treeDataEdge, 0);
+
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
diff --git a/src/meshTools/indexedOctree/treeDataEdge.H b/src/meshTools/indexedOctree/treeDataEdge.H
index 916bad870e..5dad366b3b 100644
--- a/src/meshTools/indexedOctree/treeDataEdge.H
+++ b/src/meshTools/indexedOctree/treeDataEdge.H
@@ -36,8 +36,6 @@ SourceFiles
 #ifndef treeDataEdge_H
 #define treeDataEdge_H
 
-#include "treeBoundBox.H"
-#include "pointField.H"
 #include "treeBoundBoxList.H"
 #include "linePointRef.H"
 
diff --git a/src/meshTools/indexedOctree/treeDataPoint.C b/src/meshTools/indexedOctree/treeDataPoint.C
index c64fc04ee5..5303156518 100644
--- a/src/meshTools/indexedOctree/treeDataPoint.C
+++ b/src/meshTools/indexedOctree/treeDataPoint.C
@@ -32,6 +32,11 @@ Description
 #include "polyMesh.H"
 #include "triangleFuncs.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(Foam::treeDataPoint, 0);
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
diff --git a/src/meshTools/indexedOctree/treeDataTriSurface.C b/src/meshTools/indexedOctree/treeDataTriSurface.C
index 06a496ee88..ea519291e3 100644
--- a/src/meshTools/indexedOctree/treeDataTriSurface.C
+++ b/src/meshTools/indexedOctree/treeDataTriSurface.C
@@ -22,8 +22,6 @@ 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 "treeDataTriSurface.H"
@@ -226,7 +224,7 @@ Foam::label Foam::treeDataTriSurface::getVolumeType
 
     if (!pHit.hit())
     {
-        FatalErrorIn("treeDataTriSurface::getVolumeType")
+        FatalErrorIn("treeDataTriSurface::getVolumeType(..)")
             << "treeBb:" << treeBb
             << " sample:" << sample
             << " pHit:" << pHit
@@ -238,7 +236,8 @@ Foam::label Foam::treeDataTriSurface::getVolumeType
         surface_,
         sample,
         pHit.index(),
-        pHit.hitPoint()
+        pHit.hitPoint(),
+        indexedOctree<treeDataTriSurface>::perturbTol()
     );
 
     if (t == triSurfaceTools::UNKNOWN)
@@ -255,12 +254,13 @@ Foam::label Foam::treeDataTriSurface::getVolumeType
     }
     else
     {
-        FatalErrorIn("treeDataTriSurface::getVolumeType")
+        FatalErrorIn("treeDataTriSurface::getVolumeType(..)")
             << "problem" << abort(FatalError);
         return indexedOctree<treeDataTriSurface>::UNKNOWN;
     }
 }
 
+
 // Check if any point on triangle is inside cubeBb.
 bool Foam::treeDataTriSurface::overlaps
 (
@@ -305,6 +305,7 @@ bool Foam::treeDataTriSurface::overlaps
     // Now we have the difficult case: all points are outside but connecting
     // edges might go through cube. Use fast intersection of bounding box.
 
+    //return triangleFuncs::intersectBbExact(p0, p1, p2, cubeBb);
     return triangleFuncs::intersectBb(p0, p1, p2, cubeBb);
 }
 
@@ -445,7 +446,15 @@ bool Foam::treeDataTriSurface::intersects
 
     const vector dir(end - start);
 
-    pointHit inter = tri.intersection(start, dir, intersection::HALF_RAY);
+    // Use relative tolerance (from octree) to determine intersection.
+
+    pointHit inter = tri.intersection
+    (
+        start,
+        dir,
+        intersection::HALF_RAY,
+        indexedOctree<treeDataTriSurface>::perturbTol()
+    );
 
     if (inter.hit() && inter.distance() <= 1)
     {
@@ -459,6 +468,16 @@ bool Foam::treeDataTriSurface::intersects
     {
         return false;
     }
+
+
+    //- Using exact intersections
+    //pointHit info = f.tri(points).intersectionExact(start, end);
+    //
+    //if (info.hit())
+    //{
+    //    intersectionPoint = info.hitPoint();
+    //}
+    //return info.hit();
 }
 
 
diff --git a/src/meshTools/meshSearch/meshSearch.C b/src/meshTools/meshSearch/meshSearch.C
index ef1aa46f67..2c803a2772 100644
--- a/src/meshTools/meshSearch/meshSearch.C
+++ b/src/meshTools/meshSearch/meshSearch.C
@@ -27,9 +27,11 @@ License
 #include "meshSearch.H"
 #include "polyMesh.H"
 #include "indexedOctree.H"
-#include "pointIndexHit.H"
 #include "DynamicList.H"
 #include "demandDrivenData.H"
+#include "treeDataCell.H"
+#include "treeDataFace.H"
+#include "treeDataPoint.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
diff --git a/src/meshTools/meshSearch/meshSearch.H b/src/meshTools/meshSearch/meshSearch.H
index 7b99a7159c..c8d44bdd48 100644
--- a/src/meshTools/meshSearch/meshSearch.H
+++ b/src/meshTools/meshSearch/meshSearch.H
@@ -37,9 +37,6 @@ SourceFiles
 #ifndef meshSearch_H
 #define meshSearch_H
 
-#include "treeDataCell.H"
-#include "treeDataFace.H"
-#include "treeDataPoint.H"
 #include "pointIndexHit.H"
 #include "Cloud.H"
 #include "passiveParticle.H"
@@ -51,6 +48,10 @@ namespace Foam
 
 // Forward declaration of classes
 class polyMesh;
+class treeDataCell;
+class treeDataFace;
+class treeDataPoint;
+template<class Type> class indexedOctree;
 
 /*---------------------------------------------------------------------------*\
                            Class meshSearch Declaration
diff --git a/src/meshTools/octree/treeBoundBox.C b/src/meshTools/octree/treeBoundBox.C
index ce62c47192..6133bc3106 100644
--- a/src/meshTools/octree/treeBoundBox.C
+++ b/src/meshTools/octree/treeBoundBox.C
@@ -309,12 +309,14 @@ bool Foam::treeBoundBox::overlaps
 //     This makes sure that posBits tests 'inside'
 bool Foam::treeBoundBox::intersects
 (
+    const point& overallStart,
+    const vector& overallVec,
     const point& start,
     const point& end,
-    point& pt
+    point& pt,
+    direction& ptOnFaces
 ) const
 {
-    const vector vec(end - start);
     const direction endBits = posBits(end);
     pt = start;
 
@@ -325,80 +327,106 @@ bool Foam::treeBoundBox::intersects
         if (ptBits == 0)
         {
             // pt inside bb
+            ptOnFaces = faceBits(pt);
             return true;
         }
 
         if ((ptBits & endBits) != 0)
         {
             // pt and end in same block outside of bb
+            ptOnFaces = faceBits(pt);
             return false;
         }
 
         if (ptBits & LEFTBIT)
         {
             // Intersect with plane V=min, n=-1,0,0
-            if (Foam::mag(vec.x()) > VSMALL)
+            if (Foam::mag(overallVec.x()) > VSMALL)
             {
-                scalar s = (min().x() - pt.x())/vec.x();
+                scalar s = (min().x() - overallStart.x())/overallVec.x();
+                pt.x() = min().x();
+                pt.y() = overallStart.y() + overallVec.y()*s;
+                pt.z() = overallStart.z() + overallVec.z()*s;
+            }
+            else
+            {
+                // Vector not in x-direction. But still intersecting bb planes.
+                // So must be close - just snap to bb.
                 pt.x() = min().x();
-                pt.y() = pt.y() + vec.y()*s;
-                pt.z() = pt.z() + vec.z()*s;
             }
         }
-        if (ptBits & RIGHTBIT)
+        else if (ptBits & RIGHTBIT)
         {
             // Intersect with plane V=max, n=1,0,0
-            if (Foam::mag(vec.x()) > VSMALL)
+            if (Foam::mag(overallVec.x()) > VSMALL)
+            {
+                scalar s = (max().x() - overallStart.x())/overallVec.x();
+                pt.x() = max().x();
+                pt.y() = overallStart.y() + overallVec.y()*s;
+                pt.z() = overallStart.z() + overallVec.z()*s;
+            }
+            else
             {
-                scalar s = (max().x() - pt.x())/vec.x();
                 pt.x() = max().x();
-                pt.y() = pt.y() + vec.y()*s;
-                pt.z() = pt.z() + vec.z()*s;
             }
         }
-
-        if (ptBits & BOTTOMBIT)
+        else if (ptBits & BOTTOMBIT)
         {
             // Intersect with plane V=min, n=0,-1,0
-            if (Foam::mag(vec.y()) > VSMALL)
+            if (Foam::mag(overallVec.y()) > VSMALL)
             {
-                scalar s = (min().y() - pt.y())/vec.y();
-                pt.x() = pt.x() + vec.x()*s;
+                scalar s = (min().y() - overallStart.y())/overallVec.y();
+                pt.x() = overallStart.x() + overallVec.x()*s;
                 pt.y() = min().y();
-                pt.z() = pt.z() + vec.z()*s;
+                pt.z() = overallStart.z() + overallVec.z()*s;
+            }
+            else
+            {
+                pt.x() = min().y();
             }
         }
-        if (ptBits & TOPBIT)
+        else if (ptBits & TOPBIT)
         {
             // Intersect with plane V=max, n=0,1,0
-            if (Foam::mag(vec.y()) > VSMALL)
+            if (Foam::mag(overallVec.y()) > VSMALL)
+            {
+                scalar s = (max().y() - overallStart.y())/overallVec.y();
+                pt.x() = overallStart.x() + overallVec.x()*s;
+                pt.y() = max().y();
+                pt.z() = overallStart.z() + overallVec.z()*s;
+            }
+            else
             {
-                scalar s = (max().y() - pt.y())/vec.y();
-                pt.x() = pt.x() + vec.x()*s;
                 pt.y() = max().y();
-                pt.z() = pt.z() + vec.z()*s;
             }
         }
-
-        if (ptBits & BACKBIT)
+        else if (ptBits & BACKBIT)
         {
             // Intersect with plane V=min, n=0,0,-1
-            if (Foam::mag(vec.z()) > VSMALL)
+            if (Foam::mag(overallVec.z()) > VSMALL)
+            {
+                scalar s = (min().z() - overallStart.z())/overallVec.z();
+                pt.x() = overallStart.x() + overallVec.x()*s;
+                pt.y() = overallStart.y() + overallVec.y()*s;
+                pt.z() = min().z();
+            }
+            else
             {
-                scalar s = (min().z() - pt.z())/vec.z();
-                pt.x() = pt.x() + vec.x()*s;
-                pt.y() = pt.y() + vec.y()*s;
                 pt.z() = min().z();
             }
         }
-        if (ptBits & FRONTBIT)
+        else if (ptBits & FRONTBIT)
         {
             // Intersect with plane V=max, n=0,0,1
-            if (Foam::mag(vec.z()) > VSMALL)
+            if (Foam::mag(overallVec.z()) > VSMALL)
+            {
+                scalar s = (max().z() - overallStart.z())/overallVec.z();
+                pt.x() = overallStart.x() + overallVec.x()*s;
+                pt.y() = overallStart.y() + overallVec.y()*s;
+                pt.z() = max().z();
+            }
+            else
             {
-                scalar s = (max().z() - pt.z())/vec.z();
-                pt.x() = pt.x() + vec.x()*s;
-                pt.y() = pt.y() + vec.y()*s;
                 pt.z() = max().z();
             }
         }
@@ -406,6 +434,18 @@ bool Foam::treeBoundBox::intersects
 }
 
 
+bool Foam::treeBoundBox::intersects
+(
+    const point& start,
+    const point& end,
+    point& pt
+) const
+{
+    direction ptBits;
+    return intersects(start, end-start, start, end, pt, ptBits);
+}
+
+
 // this.bb fully contains bb
 bool Foam::treeBoundBox::contains(const treeBoundBox& bb) const
 {
@@ -452,6 +492,40 @@ bool Foam::treeBoundBox::contains(const vector& dir, const point& pt) const
 }
 
 
+// Code position of pt on bounding box faces
+Foam::direction Foam::treeBoundBox::faceBits(const point& pt) const
+{
+    direction faceBits = 0;
+    if (pt.x() == min().x())
+    {
+        faceBits |= LEFTBIT;
+    }
+    else if (pt.x() == max().x())
+    {
+        faceBits |= RIGHTBIT;
+    }
+
+    if (pt.y() == min().y())
+    {
+        faceBits |= BOTTOMBIT;
+    }
+    else if (pt.y() == max().y())
+    {
+        faceBits |= TOPBIT;
+    }
+
+    if (pt.z() == min().z())
+    {
+        faceBits |= BACKBIT;
+    }
+    else if (pt.z() == max().z())
+    {
+        faceBits |= FRONTBIT;
+    }
+    return faceBits;
+}
+
+
 // Code position of point relative to box
 Foam::direction Foam::treeBoundBox::posBits(const point& pt) const
 {
@@ -461,7 +535,7 @@ Foam::direction Foam::treeBoundBox::posBits(const point& pt) const
     {
         posBits |= LEFTBIT;
     }
-    if (pt.x() > max().x())
+    else if (pt.x() > max().x())
     {
         posBits |= RIGHTBIT;
     }
@@ -470,7 +544,7 @@ Foam::direction Foam::treeBoundBox::posBits(const point& pt) const
     {
         posBits |= BOTTOMBIT;
     }
-    if (pt.y() > max().y())
+    else if (pt.y() > max().y())
     {
         posBits |= TOPBIT;
     }
@@ -479,7 +553,7 @@ Foam::direction Foam::treeBoundBox::posBits(const point& pt) const
     {
         posBits |= BACKBIT;
     }
-    if (pt.z() > max().z())
+    else if (pt.z() > max().z())
     {
         posBits |= FRONTBIT;
     }
diff --git a/src/meshTools/octree/treeBoundBox.H b/src/meshTools/octree/treeBoundBox.H
index 77498b2935..69791e04a8 100644
--- a/src/meshTools/octree/treeBoundBox.H
+++ b/src/meshTools/octree/treeBoundBox.H
@@ -120,12 +120,12 @@ public:
         enum faceBit
         {
             NOFACE    = 0,
-            LEFTBIT   = 0x1 << LEFT,
-            RIGHTBIT  = 0x1 << RIGHT,
-            BOTTOMBIT = 0x1 << BOTTOM,
-            TOPBIT    = 0x1 << TOP,
-            BACKBIT   = 0x1 << BACK,
-            FRONTBIT  = 0x1 << FRONT,
+            LEFTBIT   = 0x1 << LEFT,    //1
+            RIGHTBIT  = 0x1 << RIGHT,   //2
+            BOTTOMBIT = 0x1 << BOTTOM,  //4
+            TOPBIT    = 0x1 << TOP,     //8
+            BACKBIT   = 0x1 << BACK,    //16
+            FRONTBIT  = 0x1 << FRONT,   //32
         };
 
         //- Edges codes.
@@ -265,9 +265,25 @@ public:
             //- Overlaps boundingSphere (centre + sqr(radius))?
             bool overlaps(const point&, const scalar radiusSqr) const;
 
-            //- Intersects segment; set point to intersection position,
+            //- Intersects segment; set point to intersection position and face,
             //  return true if intersection found.
-            // (intPt argument used during calculation even if not intersecting)
+            //  (pt argument used during calculation even if not intersecting).
+            //  Calculates intersections from outside supplied vector
+            //  (overallStart, overallVec). This is so when
+            //  e.g. tracking through lots of consecutive boxes
+            // (typical octree) we're not accumulating truncation errors. Set
+            // to start, (end-start) if not used.
+            bool intersects
+            (
+                const point& overallStart,
+                const vector& overallVec,
+                const point& start,
+                const point& end,
+                point& pt,
+                direction& ptBits
+            ) const;
+
+            //- Like above but does not return faces point is on
             bool intersects
             (
                 const point& start,
@@ -285,6 +301,9 @@ public:
             //  dir would cause it to go inside.
             bool contains(const vector& dir, const point&) const;
 
+            //- Code position of pt on bounding box faces
+            direction faceBits(const point& pt) const;
+
             //- Position of point relative to bounding box
             direction posBits(const point&) const;
 
diff --git a/src/meshTools/triSurface/orientedSurface/orientedSurface.C b/src/meshTools/triSurface/orientedSurface/orientedSurface.C
index 10e1fa8f2f..113bdf58f4 100644
--- a/src/meshTools/triSurface/orientedSurface/orientedSurface.C
+++ b/src/meshTools/triSurface/orientedSurface/orientedSurface.C
@@ -238,7 +238,8 @@ void Foam::orientedSurface::propagateOrientation
         s,
         samplePoint,
         nearestFaceI,
-        nearestPt
+        nearestPt,
+        10*SMALL
     );
 
     if (side == triSurfaceTools::UNKNOWN)
diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
index 08138f450b..b89afdab63 100644
--- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
+++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
@@ -2203,7 +2203,8 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
     const triSurface& surf,
     const point& sample,
     const label nearestFaceI,   // nearest face
-    const point& nearestPoint   // nearest point on nearest face
+    const point& nearestPoint,  // nearest point on nearest face
+    const scalar tol
 )
 {
     const labelledTri& f = surf[nearestFaceI];
@@ -2217,7 +2218,7 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
         points[f[0]],
         points[f[1]],
         points[f[2]]
-    ).classify(nearestPoint, 1E-6, nearType, nearLabel);
+    ).classify(nearestPoint, tol, nearType, nearLabel);
 
     if (nearType == triPointRef::NONE)
     {
diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H
index 8ddac22a29..b88d460b77 100644
--- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H
+++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H
@@ -454,13 +454,14 @@ public:
 
         //- Given nearest point (to sample) on surface determines which side
         //  sample is. Uses either face normal, edge normal or point normal
-        //  (non-trivial). Feed in output of e.g. triangle::classify.
+        //  (non-trivial). Uses triangle::classify.
         static sideType surfaceSide
         (
             const triSurface& surf,
             const point& sample,
             const label nearestFaceI,   // nearest face
-            const point& nearestPt      // nearest point on nearest face
+            const point& nearestPt,     // nearest point on nearest face
+            const scalar tol            // tolerance for nearness test.
         );
 
     // Triangulation of faces
diff --git a/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C b/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C
index be83b6326c..1bef529b03 100644
--- a/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C
+++ b/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C
@@ -22,8 +22,6 @@ 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 "triangleFuncs.H"
@@ -31,8 +29,6 @@ Description
 #include "treeBoundBox.H"
 #include "SortableList.H"
 #include "boolList.H"
-#include "scalar.H"
-
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -170,7 +166,6 @@ bool Foam::triangleFuncs::intersectAxesBundle
 }
 
 
-
 // Intersect triangle with bounding box. Return true if
 // any of the faces of bb intersect triangle.
 // Note: so returns false if triangle inside bb.
@@ -262,6 +257,272 @@ bool Foam::triangleFuncs::intersectBb
 }
 
 
+//// Intersect triangle with bounding box. Return true if
+//// any of the faces of bb intersect triangle.
+//// Note: so returns false if triangle inside bb.
+//bool Foam::triangleFuncs::intersectBbExact
+//(
+//    const point& p0,
+//    const point& p1,
+//    const point& p2,
+//    const treeBoundBox& cubeBb
+//)
+//{
+//    const point& min = cubeBb.min();
+//    const point& max = cubeBb.max();
+//
+//    const point& cube0 = min;
+//    const point  cube1(min.x(), min.y(), max.z());
+//    const point  cube2(max.x(), min.y(), max.z());
+//    const point  cube3(max.x(), min.y(), min.z());
+//
+//    const point  cube4(min.x(), max.y(), min.z());
+//    const point  cube5(min.x(), max.y(), max.z());
+//    const point& cube6 = max;
+//    const point  cube7(max.x(), max.y(), min.z());
+//
+//    // Test intersection of triangle with twelve edges of box.
+//    {
+//        triPointRef tri(p0, p1, p2);
+//        if (tri.intersectionExact(cube0, cube1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(cube1, cube2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(cube2, cube3).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(cube3, cube0).hit())
+//        {   
+//            return true;
+//        }
+//
+//        if (tri.intersectionExact(cube4, cube5).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(cube5, cube6).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(cube6, cube7).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(cube7, cube4).hit())
+//        {
+//            return true;
+//        }
+//
+//        if (tri.intersectionExact(cube0, cube4).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(cube1, cube5).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(cube2, cube6).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(cube3, cube7).hit())
+//        {
+//            return true;
+//        }
+//    }
+//    // Test intersection of triangle edges with bounding box
+//    {
+//        triPointRef tri(cube0, cube1, cube2);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//    {
+//        triPointRef tri(cube2, cube3, cube0);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//    {
+//        triPointRef tri(cube4, cube5, cube6);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//    {
+//        triPointRef tri(cube6, cube7, cube4);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//
+//
+//    {
+//        triPointRef tri(cube4, cube5, cube1);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//    {
+//        triPointRef tri(cube1, cube0, cube4);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//    {
+//        triPointRef tri(cube7, cube6, cube2);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//    {
+//        triPointRef tri(cube2, cube3, cube7);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//
+//    {
+//        triPointRef tri(cube0, cube4, cube7);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//    {
+//        triPointRef tri(cube7, cube3, cube0);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//    {
+//        triPointRef tri(cube1, cube5, cube6);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//    {
+//        triPointRef tri(cube6, cube2, cube1);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//    return false;
+//}
+
+
 bool Foam::triangleFuncs::intersect
 (
     const point& va0,
@@ -470,145 +731,4 @@ bool Foam::triangleFuncs::intersect
 }
 
 
-bool Foam::triangleFuncs::classify
-(
-    const point& baseVertex,
-    const vector& E0,
-    const vector& E1,
-    const vector& n,
-    const point& pInter,
-    const scalar tol,
-    label& nearType,
-    label& nearLabel
-)
-{
-    // Get largest component of normal
-    scalar magX = Foam::mag(n.x());
-    scalar magY = Foam::mag(n.y());
-    scalar magZ = Foam::mag(n.z());
-
-    label i0 = -1;
-    if ((magX >= magY) && (magX >= magZ))
-    {
-        i0 = 0;
-    }
-    else if ((magY >= magX) && (magY >= magZ))
-    {
-        i0 = 1;
-    }
-    else
-    {
-        i0 = 2;
-    }
-
-    // Get other components
-    label i1 = (i0 + 1) % 3;
-    label i2 = (i1 + 1) % 3;
-
-    scalar u1 = E0[i1];
-    scalar v1 = E0[i2];
-
-    scalar u2 = E1[i1];
-    scalar v2 = E1[i2];
-
-    scalar det = v2*u1 - u2*v1;
-
-    scalar u0 = pInter[i1] - baseVertex[i1];
-    scalar v0 = pInter[i2] - baseVertex[i2];
-
-    scalar alpha = 0;
-    scalar beta = 0;
-
-    bool hit = false;
-    
-    if (Foam::mag(u1) < ROOTVSMALL)
-    {
-        beta = u0/u2;
-
-        alpha = (v0 - beta*v2)/v1;
-
-        hit = ((alpha >= 0) && ((alpha + beta) <= 1));
-    }
-    else
-    {
-        beta = (v0*u1 - u0*v1)/det;
-
-        alpha = (u0 - beta*u2)/u1;
-
-        hit = ((alpha >= 0) && ((alpha + beta) <= 1));
-    }
-
-    //
-    // Now alpha, beta are the coordinates in the triangle local coordinate
-    // system E0, E1
-    //
-
-    nearType = NONE;
-    nearLabel = -1;
-
-    if (mag(alpha+beta-1) < tol)
-    {
-        // On edge between vert 1-2 (=edge 1)
-
-        if (mag(alpha) < tol)
-        {
-            nearType = POINT;
-            nearLabel = 2;
-        }
-        else if (mag(beta) < tol)
-        {
-            nearType = POINT;
-            nearLabel = 1;
-        }
-        else if ((alpha >= 0) && (alpha <= 1) && (beta >= 0) && (beta <= 1))
-        {
-            nearType = EDGE;
-            nearLabel = 1;
-        }
-    }
-    else if (mag(alpha) < tol)
-    {
-        // On edge between vert 2-0 (=edge 2)
-
-        if (mag(beta) < tol)
-        {
-            nearType = POINT;
-            nearLabel = 0;
-        }
-        else if (mag(beta-1) < tol)
-        {
-            nearType = POINT;
-            nearLabel = 2;
-        }
-        else if ((beta >= 0) && (beta <= 1))
-        {
-            nearType = EDGE;
-            nearLabel = 2;
-        }
-    }
-    else if (mag(beta) < tol)
-    {
-        // On edge between vert 0-1 (= edge 0)
-
-        if (mag(alpha) < tol)
-        {
-            nearType = POINT;
-            nearLabel = 0;
-        }
-        else if (mag(alpha-1) < tol)
-        {
-            nearType = POINT;
-            nearLabel = 1;
-        }
-        else if ((alpha >= 0) && (alpha <= 1))
-        {
-            nearType = EDGE;
-            nearLabel = 0;
-        }
-    }
-
-    return hit;
-}
-
-
 // ************************************************************************* //
diff --git a/src/meshTools/triSurface/triangleFuncs/triangleFuncs.H b/src/meshTools/triSurface/triangleFuncs/triangleFuncs.H
index 1180e10195..0aa2a4c1c4 100644
--- a/src/meshTools/triSurface/triangleFuncs/triangleFuncs.H
+++ b/src/meshTools/triSurface/triangleFuncs/triangleFuncs.H
@@ -147,25 +147,6 @@ public:
         point& pInter0,
         point& pInter1
     );
-
-
-    //- Classify point on triangle plane w.r.t. triangle edges.
-    //  - inside/outside (returned)
-    //  - near point (0, 1, 2)
-    //  - near edge (0, 1, 2). Note: edges are counted from starting vertex
-    //    so e.g. edge 2 is from f[2] to f[0]
-    static bool classify
-    (
-        const point& baseVertex,
-        const vector& E0,
-        const vector& E1,
-        const vector& n,
-        const point& pInter,
-        const scalar tol,
-        label& nearType,
-        label& nearLabel
-    );
-
 };
 
 
diff --git a/tutorials/mesh/snappyHexMesh/iglooWithFridges/constant/triSurface/fridgeA.eMesh b/tutorials/mesh/snappyHexMesh/iglooWithFridges/constant/triSurface/fridgeA.eMesh
new file mode 100644
index 0000000000..ff7aab4b8f
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/iglooWithFridges/constant/triSurface/fridgeA.eMesh
@@ -0,0 +1,58 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.3.1                                 |
+|   \\  /    A nd           | Web:      http://www.openfoam.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+FoamFile
+{
+    version 2.0;
+    format ascii;
+
+    root "..";
+    case "cavity";
+    instance ""constant"";
+    local "polyMesh";
+
+    class featureEdgeMesh;
+    object points;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Points
+(
+(1.99 1.99 0.01)     //0
+(3.01 1.99 0.01)     //1
+(3.01 3.01 0.01)     //2
+(1.99 3.01 0.01)     //3
+
+(1.99 1.99 2.01)     //4
+(3.01 1.99 2.01)     //5
+(3.01 3.01 2.01)     //6
+(1.99 3.01 2.01)     //7
+)
+
+// Edges
+(
+(0 1)
+(1 2)
+(2 3)
+(3 0)
+
+(4 5)
+(5 6)
+(6 7)
+(7 4)
+
+(0 4)
+(1 5)
+(2 6)
+(3 7)
+)
+
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/snappyHexMesh/iglooWithFridges/system/snappyHexMeshDict b/tutorials/mesh/snappyHexMesh/iglooWithFridges/system/snappyHexMeshDict
index 19ef2afec1..be30045c6d 100644
--- a/tutorials/mesh/snappyHexMesh/iglooWithFridges/system/snappyHexMeshDict
+++ b/tutorials/mesh/snappyHexMesh/iglooWithFridges/system/snappyHexMeshDict
@@ -143,10 +143,10 @@ castellatedMeshControls
     // This is a featureEdgeMesh, read from constant/triSurface for now.
     features
     (
-//        {
-//            file "fridgeA.eMesh";
-//            level 3;
-//        }
+        {
+            file "fridgeA.eMesh";
+            level 3;
+        }
     );
 
 
@@ -308,6 +308,10 @@ addLayersControls
 
     // Create buffer region for new layer terminations
     nBufferCellsNoExtrude 0;
+
+
+    // Overall max number of layer addition iterations
+    nLayerIter 50;
 }
 
 
diff --git a/tutorials/mesh/snappyHexMesh/motorBike/system/snappyHexMeshDict b/tutorials/mesh/snappyHexMesh/motorBike/system/snappyHexMeshDict
index 3f0ff0417e..0d84177d8b 100644
--- a/tutorials/mesh/snappyHexMesh/motorBike/system/snappyHexMeshDict
+++ b/tutorials/mesh/snappyHexMesh/motorBike/system/snappyHexMeshDict
@@ -160,7 +160,7 @@ snapControls
     //- Relative distance for points to be attracted by surface feature point
     //  or edge. True distance is this factor times local
     //  maximum edge length.
-    tolerance       4;
+    tolerance 4.0;
 
     //- Number of mesh displacement relaxation iterations.
     nSolveIter      30;
@@ -182,344 +182,14 @@ addLayersControls
             nSurfaceLayers  1;
         }
 
-        motorBike_frt-fairing:001%1
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_windshield:002%2
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rr-wh-rim:005%5
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rr-wh-rim:010%10
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_fr-wh-rim:011%11
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_fr-wh-brake-disk:012%12
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_frame:016-shadow%13
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rear-susp:014%14
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rear-susp:014-shadow%15
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_frame:016%16
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rr-wh-rim:005-shadow%17
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rr-wh-chain-hub:022%22
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rearseat%24
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_frt-fairing%25
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_windshield%26
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_headlights%27
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_driversseat%28
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rear-body%29
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_fuel-tank%30
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_exhaust%31
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rr-wh-rim%32
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_fr-mud-guard%33
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_fr-wh-rim%34
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_fr-wh-brake-disk%35
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_fr-brake-caliper%36
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_fr-wh-tyre%37
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_hbars%38
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_fr-forks%39
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_chain%40
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rr-wh-tyre%41
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_square-dial%42
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_round-dial%43
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_dial-holder%44
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rear-susp%45
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rear-brake-lights%46
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rear-light-bracket%47
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_frame%48
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rear-mud-guard%49
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rear-susp-spring-damp%50
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_fairing-inner-plate%51
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_clutch-housing%52
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_radiator%53
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_water-pipe%54
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_water-pump%55
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_engine%56
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rear-shock-link%57
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rear-brake-fluid-pot-bracket%58
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rear-brake-fluid-pot%59
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_footpeg%60
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rr-wh-chain-hub%61
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rear-brake-caliper%62
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rider-helmet%65
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rider-visor%66
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rider-boots%67
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rider-gloves%68
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rider-body%69
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_frame:0%70
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_frt-fairing:001-shadow%74
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_windshield-shadow%75
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_fr-mud-guard-shadow%81
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_fr-wh-brake-disk-shadow%83
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rear-mud-guard-shadow%84
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rear-susp-spring-damp-shadow%85
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_radiator-shadow%86
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rear-shock-link-shadow%87
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rear-brake-fluid-pot-bracket-shadow%88
-        {
-            nSurfaceLayers  1;
-        }
-
-        motorBike_rr-wh-chain-hub-shadow%89
+        "motorBike_.*"
         {
             nSurfaceLayers  1;
         }
     }
 
     // Expansion factor for layer mesh
-    expansionRatio  1;
+    expansionRatio 1.0;
 
     //- Wanted thickness of final added cell layer. If multiple layers
     //  is the
@@ -569,6 +239,10 @@ addLayersControls
 
     // Create buffer region for new layer terminations
     nBufferCellsNoExtrude 0;
+
+
+    // Overall max number of layer addition iterations
+    nLayerIter 50;
 }