diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
index 910732fa7a..bcdd163834 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
@@ -78,10 +78,8 @@ SourceFiles
 namespace Foam
 {
 
-typedef treeDataPrimitivePatch<face, List, const pointField, point>
-    treeDataBPatch;
-
 typedef PrimitivePatch<face, List, const pointField, point> bPatch;
+typedef treeDataPrimitivePatch<bPatch> treeDataBPatch;
 
 /*---------------------------------------------------------------------------*\
                   Class backgroundMeshDecomposition Declaration
diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
index 395ee77cbe..e64f6cb78f 100644
--- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
+++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
@@ -1018,6 +1018,11 @@ int main(int argc, char *argv[])
 
     // Update proc maps
     if (cellProcAddressing.headerOk())
+    if
+    (
+        cellProcAddressing.headerOk()
+     && cellProcAddressing.size() == mesh.nCells()
+    )
     {
         Info<< "Renumbering processor cell decomposition map "
             << cellProcAddressing.name() << endl;
@@ -1028,6 +1033,11 @@ int main(int argc, char *argv[])
         );
     }
     if (faceProcAddressing.headerOk())
+    if
+    (
+        faceProcAddressing.headerOk()
+     && faceProcAddressing.size() == mesh.nFaces()
+    )
     {
         Info<< "Renumbering processor face decomposition map "
             << faceProcAddressing.name() << endl;
@@ -1054,6 +1064,11 @@ int main(int argc, char *argv[])
         }
     }
     if (pointProcAddressing.headerOk())
+    if
+    (
+        pointProcAddressing.headerOk()
+     && pointProcAddressing.size() == mesh.nPoints()
+    )
     {
         Info<< "Renumbering processor point decomposition map "
             << pointProcAddressing.name() << endl;
@@ -1173,21 +1188,41 @@ int main(int argc, char *argv[])
 
     mesh.write();
     if (cellProcAddressing.headerOk())
+    if
+    (
+        cellProcAddressing.headerOk()
+     && cellProcAddressing.size() == mesh.nCells()
+    )
     {
         cellProcAddressing.instance() = mesh.facesInstance();
         cellProcAddressing.write();
     }
     if (faceProcAddressing.headerOk())
+    if
+    (
+        faceProcAddressing.headerOk()
+     && faceProcAddressing.size() == mesh.nFaces()
+    )
     {
         faceProcAddressing.instance() = mesh.facesInstance();
         faceProcAddressing.write();
     }
     if (pointProcAddressing.headerOk())
+    if
+    (
+        pointProcAddressing.headerOk()
+     && pointProcAddressing.size() == mesh.nPoints()
+    )
     {
         pointProcAddressing.instance() = mesh.facesInstance();
         pointProcAddressing.write();
     }
     if (boundaryProcAddressing.headerOk())
+    if
+    (
+        boundaryProcAddressing.headerOk()
+     && boundaryProcAddressing.size() == mesh.boundaryMesh().size()
+    )
     {
         boundaryProcAddressing.instance() = mesh.facesInstance();
         boundaryProcAddressing.write();
diff --git a/applications/utilities/preProcessing/setFields/setFields.C b/applications/utilities/preProcessing/setFields/setFields.C
index 7b690d46f2..dfde514992 100644
--- a/applications/utilities/preProcessing/setFields/setFields.C
+++ b/applications/utilities/preProcessing/setFields/setFields.C
@@ -92,7 +92,15 @@ bool setCellFieldType
                 field.boundaryField()[patchi].patchInternalField();
         }
 
-        field.write();
+        if (!field.write())
+        {
+            FatalErrorIn
+            (
+                "void setCellFieldType"
+                "(const fvMesh& mesh, const labelList& selectedCells,"
+                "Istream& fieldValueStream)"
+            ) << "Failed writing field " << fieldName << endl;
+        }
     }
     else
     {
@@ -260,7 +268,15 @@ bool setFaceFieldType
             }
         }
 
-        field.write();
+        if (!field.write())
+        {
+            FatalErrorIn
+            (
+                "void setFaceFieldType"
+                "(const fvMesh& mesh, const labelList& selectedFaces,"
+                "Istream& fieldValueStream)"
+            )   << "Failed writing field " << field.name() << exit(FatalError);
+        }
     }
     else
     {
diff --git a/src/dynamicMesh/boundaryMesh/boundaryMesh.C b/src/dynamicMesh/boundaryMesh/boundaryMesh.C
index 51be127dd7..4a705acead 100644
--- a/src/dynamicMesh/boundaryMesh/boundaryMesh.C
+++ b/src/dynamicMesh/boundaryMesh/boundaryMesh.C
@@ -924,10 +924,10 @@ Foam::labelList Foam::boundaryMesh::getNearest
     // Create the octrees
     indexedOctree
     <
-        treeDataPrimitivePatch<face, UIndirectList, const pointField&>
+        treeDataPrimitivePatch<uindirectPrimitivePatch>
     > leftTree
     (
-        treeDataPrimitivePatch<face, UIndirectList, const pointField&>
+        treeDataPrimitivePatch<uindirectPrimitivePatch>
         (
             false,          // cacheBb
             leftPatch
@@ -939,10 +939,10 @@ Foam::labelList Foam::boundaryMesh::getNearest
     );
     indexedOctree
     <
-        treeDataPrimitivePatch<face, UIndirectList, const pointField&>
+        treeDataPrimitivePatch<uindirectPrimitivePatch>
     > rightTree
     (
-        treeDataPrimitivePatch<face, UIndirectList, const pointField&>
+        treeDataPrimitivePatch<uindirectPrimitivePatch>
         (
             false,          // cacheBb
             rightPatch
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
index 2e908b7318..35d0822d2a 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -249,9 +249,10 @@ class autoSnapDriver
                     const indirectPrimitivePatch& pp,
                     const scalarField& snapDist,
 
-                    const List<List<point> >& pointFaceNormals,
+                    const List<List<point> >& pointFaceSurfNormals,
                     const List<List<point> >& pointFaceDisp,
                     const List<List<point> >& pointFaceCentres,
+                    const labelListList& pointFacePatchID,
 
                     vectorField& patchAttraction,
                     List<pointConstraint>& patchConstraints
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
index 6347d5dd16..950849e15e 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
@@ -71,14 +71,15 @@ namespace Foam
         }
     };
 
+    template<class T>
     class listPlusEqOp
     {
     public:
 
         void operator()
         (
-            List<point>& x,
-            const List<point>& y
+            List<T>& x,
+            const List<T>& y
         ) const
         {
             label sz = x.size();
@@ -486,7 +487,7 @@ void Foam::autoSnapDriver::binFeatureFaces
 
     const label pointI,
 
-    const List<List<point> >& pointFaceNormals,
+    const List<List<point> >& pointFaceSurfNormals,
     const List<List<point> >& pointFaceDisp,
     const List<List<point> >& pointFaceCentres,
 
@@ -495,7 +496,7 @@ void Foam::autoSnapDriver::binFeatureFaces
     DynamicList<label>& surfaceCount
 ) const
 {
-    const List<point>& pfNormals = pointFaceNormals[pointI];
+    const List<point>& pfNormals = pointFaceSurfNormals[pointI];
     const List<point>& pfDisp = pointFaceDisp[pointI];
     const List<point>& pfCentres = pointFaceCentres[pointI];
 
@@ -531,7 +532,7 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
     const indirectPrimitivePatch& pp,
     const scalarField& snapDist,
 
-    const List<List<point> >& pointFaceNormals,
+    const List<List<point> >& pointFaceSurfNormals,
     const List<List<point> >& pointFaceDisp,
     const List<List<point> >& pointFaceCentres,
 
@@ -590,7 +591,7 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
             snapDist,
             pointI,
 
-            pointFaceNormals,
+            pointFaceSurfNormals,
             pointFaceDisp,
             pointFaceCentres,
 
@@ -720,7 +721,7 @@ void Foam::autoSnapDriver::determineAllFeatures
     const indirectPrimitivePatch& pp,
     const scalarField& snapDist,
 
-    const List<List<point> >& pointFaceNormals,
+    const List<List<point> >& pointFaceSurfNormals,
     const List<List<point> >& pointFaceDisp,
     const List<List<point> >& pointFaceCentres,
 
@@ -873,7 +874,7 @@ void Foam::autoSnapDriver::determineFeatures
     //const vectorField& faceSurfaceNormal,
     //const vectorField& faceDisp,
     //const vectorField& faceRotation,
-    const List<List<point> >& pointFaceNormals,
+    const List<List<point> >& pointFaceSurfNormals,
     const List<List<point> >& pointFaceDisp,
     const List<List<point> >& pointFaceCentres,
 
@@ -955,7 +956,7 @@ void Foam::autoSnapDriver::determineFeatures
 
                 //faceSurfaceNormal,
                 //faceDisp,
-                pointFaceNormals,
+                pointFaceSurfNormals,
                 pointFaceDisp,
                 pointFaceCentres,
 
@@ -1132,9 +1133,10 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
     //const vectorField& faceSurfaceNormal,
     //const vectorField& faceDisp,
     //const vectorField& faceRotation,
-    const List<List<point> >& pointFaceNormals,
+    const List<List<point> >& pointFaceSurfNormals,
     const List<List<point> >& pointFaceDisp,
     const List<List<point> >& pointFaceCentres,
+    const labelListList& pointFacePatchID,
 
     vectorField& patchAttraction,
     List<pointConstraint>& patchConstraints
@@ -1181,7 +1183,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
         pp,
         snapDist,
 
-        pointFaceNormals,
+        pointFaceSurfNormals,
         pointFaceDisp,
         pointFaceCentres,
 
@@ -1340,6 +1342,131 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
     }
 
 
+
+    //MEJ: any faces that have multi-patch points only keep the multi-patch
+    //     points.
+    {
+        autoPtr<OFstream> multiPatchStr;
+        if (debug&meshRefinement::OBJINTERSECTIONS)
+        {
+            multiPatchStr.reset
+            (
+                new OFstream
+                (
+                    meshRefiner_.mesh().time().path()
+                  / "multiPatch_" + name(iter) + ".obj"
+                )
+            );
+            Pout<< nl << "Dumping removed constraints due to same-face"
+                << " multi-patch points to "
+                << multiPatchStr().name() << endl;
+        }
+
+
+        // 1. Mark points on multiple patches
+        PackedBoolList isMultiPatchPoint(pp.size());
+
+        forAll(pointFacePatchID, pointI)
+        {
+            const labelList& patches = pointFacePatchID[pointI];
+            label patch0 = patches[0];
+            for (label i = 1; i < patches.size(); i++)
+            {
+                if (patch0 != patches[i])
+                {
+                    isMultiPatchPoint[pointI] = 1u;
+                    break;
+                }
+            }
+        }
+
+        // 2. Make sure multi-patch points are also attracted
+        forAll(isMultiPatchPoint, pointI)
+        {
+            if (isMultiPatchPoint[pointI])
+            {
+                if
+                (
+                    patchConstraints[pointI].first() == 0
+                 && allPatchConstraints[pointI].first() > 0
+                )
+                {
+                    patchAttraction[pointI] = allPatchAttraction[pointI];
+                    patchConstraints[pointI] = allPatchConstraints[pointI];
+
+                    if (multiPatchStr.valid())
+                    {
+                        Pout<< "Adding constraint on multiPatchPoint:"
+                            << pp.localPoints()[pointI]
+                            << " constraint:" << patchConstraints[pointI]
+                            << " attraction:" << patchAttraction[pointI]
+                            << endl;
+                    }
+                }
+            }
+        }
+
+        // Up to here it is all parallel ok.
+
+
+        // 3. Knock out any attraction on faces with multi-patch points
+        label nChanged = 0;
+        forAll(pp.localFaces(), faceI)
+        {
+            const face& f = pp.localFaces()[faceI];
+
+            label nMultiPatchPoints = 0;
+            forAll(f, fp)
+            {
+                label pointI = f[fp];
+                if
+                (
+                    isMultiPatchPoint[pointI]
+                 && patchConstraints[pointI].first() != 0
+                )
+                {
+                    nMultiPatchPoints++;
+                }
+            }
+
+            if (nMultiPatchPoints > 0)
+            {
+                forAll(f, fp)
+                {
+                    label pointI = f[fp];
+                    if
+                    (
+                       !isMultiPatchPoint[pointI]
+                     && patchConstraints[pointI].first() != 0
+                    )
+                    {
+                        //Pout<< "Knocking out constraint"
+                        //    << " on non-multiPatchPoint:"
+                        //    << pp.localPoints()[pointI] << endl;
+                        patchAttraction[pointI] = vector::zero;
+                        patchConstraints[pointI] = pointConstraint();
+                        nChanged++;
+
+                        if (multiPatchStr.valid())
+                        {
+                            meshTools::writeOBJ
+                            (
+                                multiPatchStr(),
+                                pp.localPoints()[pointI]
+                            );
+                        }
+                    }
+                }
+            }
+        }
+
+        reduce(nChanged, sumOp<label>());
+        Info<< "Removing constraints near multi-patch points : changed "
+            << nChanged << " points" << endl;
+    }
+
+
+
     // Dump
     if (debug&meshRefinement::OBJINTERSECTIONS)
     {
@@ -1753,8 +1880,12 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     const pointField& localPoints = pp.localPoints();
     const fvMesh& mesh = meshRefiner_.mesh();
 
-    // Displacement and orientation per pp face.
+    // Displacement and orientation per pp face
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // vector from point on surface back to face centre
     vectorField faceDisp(pp.size(), vector::zero);
+    // normal of surface at point on surface
     vectorField faceSurfaceNormal(pp.size(), vector::zero);
     vectorField faceRotation(pp.size(), vector::zero);
 
@@ -1771,34 +1902,42 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     vectorField patchDisp = nearestDisp;
 
 
-    // Collect (possibly remote) face-wise data on coupled points.
+    // Collect (possibly remote) per point data of all surrounding faces
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // - faceSurfaceNormal
     // - faceDisp
     // - faceRotation
-    // - faceCentres
+    // - faceCentres&faceNormal
 
     // For now just get all surrounding face data. Expensive - should just
     // store and sync data on coupled points only (see e.g PatchToolsNormals.C)
 
-    List<List<point> > pointFaceNormals(pp.nPoints());
+    List<List<point> > pointFaceSurfNormals(pp.nPoints());
     List<List<point> > pointFaceDisp(pp.nPoints());
     List<List<point> > pointFaceCentres(pp.nPoints());
+    List<labelList>    pointFacePatchID(pp.nPoints());
 
     // Fill local data
     forAll(pp.pointFaces(), pointI)
     {
         const labelList& pFaces = pp.pointFaces()[pointI];
-        List<point>& pNormals = pointFaceNormals[pointI];
+        List<point>& pNormals = pointFaceSurfNormals[pointI];
         pNormals.setSize(pFaces.size());
         List<point>& pDisp = pointFaceDisp[pointI];
         pDisp.setSize(pFaces.size());
         List<point>& pFc = pointFaceCentres[pointI];
         pFc.setSize(pFaces.size());
+        labelList& pFid = pointFacePatchID[pointI];
+        pFid.setSize(pFaces.size());
+
         forAll(pFaces, i)
         {
-            pNormals[i] = faceSurfaceNormal[pFaces[i]];
-            pDisp[i] = faceDisp[pFaces[i]];
-            pFc[i] = pp.faceCentres()[pFaces[i]];
+            label faceI = pFaces[i];
+            pNormals[i] = faceSurfaceNormal[faceI];
+            pDisp[i] = faceDisp[faceI];
+            pFc[i] = pp.faceCentres()[faceI];
+            label meshFaceI = pp.addressing()[faceI];
+            pFid[i] = mesh.boundaryMesh().whichPatch(meshFaceI);
         }
     }
 
@@ -1806,8 +1945,8 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     (
         mesh,
         pp.meshPoints(),
-        pointFaceNormals,
-        listPlusEqOp(),
+        pointFaceSurfNormals,
+        listPlusEqOp<point>(),
         List<point>(),
         listTransform()
     );
@@ -1816,7 +1955,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         mesh,
         pp.meshPoints(),
         pointFaceDisp,
-        listPlusEqOp(),
+        listPlusEqOp<point>(),
         List<point>(),
         listTransform()
     );
@@ -1825,13 +1964,27 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         mesh,
         pp.meshPoints(),
         pointFaceCentres,
-        listPlusEqOp(),
+        listPlusEqOp<point>(),
         List<point>(),
         listTransform()
     );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        pointFacePatchID,
+        listPlusEqOp<label>(),
+        List<label>()
+    );
 
 
 
+    // Main calculation
+    // ~~~~~~~~~~~~~~~~
+    // This is the main intelligence which calculates per point the vector to
+    // attract it to the nearest surface. There are lots of possibilities
+    // here.
+
     // Nearest feature
     vectorField patchAttraction(localPoints.size(), vector::zero);
     // Constraints at feature
@@ -1845,9 +1998,10 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         pp,
         snapDist,
 
-        pointFaceNormals,
+        pointFaceSurfNormals,
         pointFaceDisp,
         pointFaceCentres,
+        pointFacePatchID,
 
         patchAttraction,
         patchConstraints
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index 2664719b26..e1b76c44e1 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -272,7 +272,11 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
     // Find all start cells of features. Is done by tracking from keepPoint.
     Cloud<trackedParticle> cloud(mesh_, IDLList<trackedParticle>());
 
-    // Create particles on whichever processor holds the keepPoint.
+
+    // Features are identical on all processors. Number them so we know
+    // what to seed. Do this on only the processor that
+    // holds the keepPoint.
+
     label cellI = -1;
     label tetFaceI = -1;
     label tetPtI = -1;
@@ -281,13 +285,25 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
 
     if (cellI != -1)
     {
+        // I am the processor that holds the keepPoint
+
         forAll(features_, featI)
         {
             const featureEdgeMesh& featureMesh = features_[featI];
             const label featureLevel = features_.levels()[featI];
-
             const labelListList& pointEdges = featureMesh.pointEdges();
 
+            // Find regions on edgeMesh
+            labelList edgeRegion;
+            label nRegions = featureMesh.regions(edgeRegion);
+
+
+            PackedBoolList regionVisited(nRegions);
+
+
+            // 1. Seed all 'knots' in edgeMesh
+
+
             forAll(pointEdges, pointI)
             {
                 if (pointEdges[pointI].size() != 2)
@@ -316,6 +332,50 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
                             pointI                          // end point
                         )
                     );
+
+                    // Mark
+                    if (pointEdges[pointI].size() > 0)
+                    {
+                        label e0 = pointEdges[pointI][0];
+                        label regionI = edgeRegion[e0];
+                        regionVisited[regionI] = 1u;
+                    }
+                }
+            }
+
+
+            // 2. Any regions that have not been visited at all? These can
+            //    only be circular regions!
+            forAll(featureMesh.edges(), edgeI)
+            {
+                if (regionVisited.set(edgeRegion[edgeI], 1u))
+                {
+                    const edge& e = featureMesh.edges()[edgeI];
+                    label pointI = e.start();
+                    if (debug)
+                    {
+                        Pout<< "Adding particle from point:" << pointI
+                            << " coord:" << featureMesh.points()[pointI]
+                            << " on circular region:" << edgeRegion[edgeI]
+                            << endl;
+                    }
+
+                    // Non-manifold point. Create particle.
+                    cloud.addParticle
+                    (
+                        new trackedParticle
+                        (
+                            mesh_,
+                            keepPoint,
+                            cellI,
+                            tetFaceI,
+                            tetPtI,
+                            featureMesh.points()[pointI],   // endpos
+                            featureLevel,                   // level
+                            featI,                          // featureMesh
+                            pointI                          // end point
+                        )
+                    );
                 }
             }
         }
@@ -329,6 +389,8 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
     trackedParticle::trackingData td(cloud, maxFeatureLevel);
 
     // Track all particles to their end position (= starting feature point)
+    // Note that the particle might have started on a different processor
+    // so this will transfer across nicely until we can start tracking proper.
     cloud.move(td, GREAT);
 
     // Reset level
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
index 3962697845..5acc8dfa4c 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
@@ -1164,50 +1164,36 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::normaliseWeights
     const word& patchName,
     const labelListList& addr,
     scalarListList& wght,
+    scalarField& wghtSum,
     const bool output
 )
 {
-    scalar minBound = VGREAT;
-    scalar maxBound = -VGREAT;
-
-    scalar tSum = 0.0;
-
     // Normalise the weights
+    wghtSum.setSize(wght.size());
     forAll(wght, faceI)
     {
         scalar s = sum(wght[faceI]);
         scalar t = s/patchAreas[faceI];
 
-        tSum += t;
-
-        if (t < minBound)
-        {
-            minBound = t;
-        }
-
-        if (t > maxBound)
-        {
-            maxBound = t;
-        }
-
         forAll(addr[faceI], i)
         {
             wght[faceI][i] /= s;
         }
+
+        wghtSum[faceI] = t;
     }
 
 
     if (output)
     {
-        const label nFace = returnReduce(wght.size(), sumOp<scalar>());
-        reduce(tSum, sumOp<scalar>());
+        const label nFace = returnReduce(wght.size(), sumOp<label>());
 
         if (nFace)
         {
             Info<< "AMI: Patch " << patchName << " weights min/max/average = "
-                << returnReduce(minBound, minOp<scalar>()) << ", "
-                << returnReduce(maxBound, maxOp<scalar>()) << ", "
-                << tSum/nFace << endl;
+                << gMin(wghtSum) << ", "
+                << gMax(wghtSum) << ", "
+                << gAverage(wghtSum) << endl;
         }
     }
 }
@@ -1227,6 +1213,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
     scalarField& srcMagSf,
     labelListList& srcAddress,
     scalarListList& srcWeights,
+    scalarField& srcWeightsSum,
     autoPtr<mapDistribute>& tgtMap
 )
 {
@@ -1468,6 +1455,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
         "source",
         srcAddress,
         srcWeights,
+        srcWeightsSum,
         false
     );
 }
@@ -1488,9 +1476,11 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
     singlePatchProc_(-999),
     srcAddress_(),
     srcWeights_(),
+    srcWeightsSum_(),
     srcNonOverlap_(),
     tgtAddress_(),
     tgtWeights_(),
+    tgtWeightsSum_(),
     treePtr_(NULL),
     triMode_(triMode),
     srcMapPtr_(NULL),
@@ -1521,9 +1511,11 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
     singlePatchProc_(-999),
     srcAddress_(),
     srcWeights_(),
+    srcWeightsSum_(),
     srcNonOverlap_(),
     tgtAddress_(),
     tgtWeights_(),
+    tgtWeightsSum_(),
     treePtr_(NULL),
     triMode_(triMode),
     srcMapPtr_(NULL),
@@ -1609,9 +1601,11 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
     singlePatchProc_(fineAMI.singlePatchProc_),
     srcAddress_(),
     srcWeights_(),
+    srcWeightsSum_(),
     srcNonOverlap_(),
     tgtAddress_(),
     tgtWeights_(),
+    tgtWeightsSum_(),
     treePtr_(NULL),
     triMode_(fineAMI.triMode_),
     srcMapPtr_(NULL),
@@ -1681,6 +1675,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
         srcMagSf_,
         srcAddress_,
         srcWeights_,
+        srcWeightsSum_,
         tgtMapPtr_
     );
 
@@ -1706,6 +1701,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
         tgtMagSf_,
         tgtAddress_,
         tgtWeights_,
+        tgtWeightsSum_,
         srcMapPtr_
     );
 
@@ -1857,8 +1853,24 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
         );
 
         // weights normalisation
-        normaliseWeights(srcMagSf_, "source", srcAddress_, srcWeights_, true);
-        normaliseWeights(tgtMagSf_, "target", tgtAddress_, tgtWeights_, true);
+        normaliseWeights
+        (
+            srcMagSf_,
+            "source",
+            srcAddress_,
+            srcWeights_,
+            srcWeightsSum_,
+            true
+        );
+        normaliseWeights
+        (
+            tgtMagSf_,
+            "target",
+            tgtAddress_,
+            tgtWeights_,
+            tgtWeightsSum_,
+            true
+        );
 
         // cache maps and reset addresses
         List<Map<label> > cMap;
@@ -1877,8 +1889,24 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
 
         calcAddressing(srcPatch, tgtPatch);
 
-        normaliseWeights(srcMagSf_, "source", srcAddress_, srcWeights_, true);
-        normaliseWeights(tgtMagSf_, "target", tgtAddress_, tgtWeights_, true);
+        normaliseWeights
+        (
+            srcMagSf_,
+            "source",
+            srcAddress_,
+            srcWeights_,
+            srcWeightsSum_,
+            true
+        );
+        normaliseWeights
+        (
+            tgtMagSf_,
+            "target",
+            tgtAddress_,
+            tgtWeights_,
+            tgtWeightsSum_,
+            true
+        );
     }
 
     if (debug)
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
index 36cd661175..b393030823 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
@@ -82,7 +82,7 @@ class AMIInterpolation
     public AMIInterpolationName
 {
     //- local typedef to octree tree-type
-    typedef treeDataPrimitivePatch<face, SubList, const pointField&> treeType;
+    typedef treeDataPrimitivePatch<TargetPatch> treeType;
 
 
     // Private data
@@ -106,6 +106,9 @@ class AMIInterpolation
             //- Weights of target faces per source face
             scalarListList srcWeights_;
 
+            //- Sum of weights of target faces per source face
+            scalarField srcWeightsSum_;
+
             //- Labels of faces that are not overlapped by any target faces
             //  (should be empty for correct functioning)
             labelList srcNonOverlap_;
@@ -122,6 +125,9 @@ class AMIInterpolation
             //- Weights of source faces per target face
             scalarListList tgtWeights_;
 
+            //- Sum of weights of source faces per target face
+            scalarField tgtWeightsSum_;
+
 
         //- Octree used to find face seeds
         autoPtr<indexedOctree<treeType> > treePtr_;
@@ -285,6 +291,7 @@ class AMIInterpolation
                 const word& patchName,
                 const labelListList& addr,
                 scalarListList& wght,
+                scalarField& wghtSum,
                 const bool output
             );
 
@@ -304,6 +311,7 @@ class AMIInterpolation
                 scalarField& srcMagSf,
                 labelListList& srcAddress,
                 scalarListList& srcWeights,
+                scalarField& srcWeightsSum,
                 autoPtr<mapDistribute>& tgtMap
             );
 
@@ -365,11 +373,14 @@ public:
                 //- Return const access to source patch weights
                 inline const scalarListList& srcWeights() const;
 
+                //- Return const access to normalisation factor of source
+                //  patch weights (i.e. the sum before normalisation)
+                inline const scalarField& srcWeightsSum() const;
+
                 //- Labels of faces that are not overlapped by any target faces
                 //  (should be empty for correct functioning)
                 inline const labelList& srcNonOverlap() const;
 
-
                 //- Source map pointer - valid only if singlePatchProc = -1
                 //  This gets source data into a form to be consumed by
                 //  tgtAddress, tgtWeights
@@ -387,7 +398,11 @@ public:
                 //- Return const access to target patch weights
                 inline const scalarListList& tgtWeights() const;
 
-                //- Target map pointer - valid only if singlePatchProc = -1
+                //- Return const access to normalisation factor of target
+                //  patch weights (i.e. the sum before normalisation)
+                inline const scalarField& tgtWeightsSum() const;
+
+                //- Target map pointer -  valid only if singlePatchProc=-1.
                 //  This gets target data into a form to be consumed by
                 //  srcAddress, srcWeights
                 inline const mapDistribute& tgtMap() const;
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H
index 3dd21b5aa3..f2782bad88 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H
@@ -63,6 +63,14 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcNonOverlap() const
 }
 
 
+template<class SourcePatch, class TargetPatch>
+inline const Foam::scalarField&
+Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcWeightsSum() const
+{
+    return srcWeightsSum_;
+}
+
+
 template<class SourcePatch, class TargetPatch>
 inline const Foam::mapDistribute&
 Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcMap() const
@@ -95,6 +103,14 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::tgtWeights() const
 }
 
 
+template<class SourcePatch, class TargetPatch>
+inline const Foam::scalarField&
+Foam::AMIInterpolation<SourcePatch, TargetPatch>::tgtWeightsSum() const
+{
+    return tgtWeightsSum_;
+}
+
+
 template<class SourcePatch, class TargetPatch>
 inline const Foam::mapDistribute&
 Foam::AMIInterpolation<SourcePatch, TargetPatch>::tgtMap() const
diff --git a/src/meshTools/indexedOctree/treeDataPrimitivePatch.C b/src/meshTools/indexedOctree/treeDataPrimitivePatch.C
index 6927cbaa50..4cd2968232 100644
--- a/src/meshTools/indexedOctree/treeDataPrimitivePatch.C
+++ b/src/meshTools/indexedOctree/treeDataPrimitivePatch.C
@@ -29,30 +29,14 @@ License
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-template
-<
-    class Face,
-    template<class> class FaceList,
-    class PointField,
-    class PointType
->
-Foam::scalar
-Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
-tolSqr = sqr(1e-6);
+template<class PatchType>
+Foam::scalar Foam::treeDataPrimitivePatch<PatchType>::tolSqr = sqr(1e-6);
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-template
-<
-    class Face,
-    template<class> class FaceList,
-    class PointField,
-    class PointType
->
-Foam::treeBoundBox
-Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
-calcBb
+template<class PatchType>
+Foam::treeBoundBox Foam::treeDataPrimitivePatch<PatchType>::calcBb
 (
     const pointField& points,
     const face& f
@@ -71,15 +55,8 @@ calcBb
 }
 
 
-template
-<
-    class Face,
-    template<class> class FaceList,
-    class PointField,
-    class PointType
->
-void Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
-update()
+template<class PatchType>
+void Foam::treeDataPrimitivePatch<PatchType>::update()
 {
     if (cacheBb_)
     {
@@ -96,18 +73,11 @@ update()
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
-template
-<
-    class Face,
-    template<class> class FaceList,
-    class PointField,
-    class PointType
->
-Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
-treeDataPrimitivePatch
+template<class PatchType>
+Foam::treeDataPrimitivePatch<PatchType>::treeDataPrimitivePatch
 (
     const bool cacheBb,
-    const PrimitivePatch<Face, FaceList, PointField, PointType>& patch
+    const PatchType& patch
 )
 :
     patch_(patch),
@@ -119,16 +89,8 @@ treeDataPrimitivePatch
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template
-<
-    class Face,
-    template<class> class FaceList,
-    class PointField,
-    class PointType
->
-Foam::pointField
-Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
-shapePoints() const
+template<class PatchType>
+Foam::pointField Foam::treeDataPrimitivePatch<PatchType>::shapePoints() const
 {
     pointField cc(patch_.size());
 
@@ -143,27 +105,10 @@ shapePoints() const
 
 //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
 //  Only makes sense for closed surfaces.
-template
-<
-    class Face,
-    template<class> class FaceList,
-    class PointField,
-    class PointType
->
-Foam::label
-Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
-getVolumeType
+template<class PatchType>
+Foam::label Foam::treeDataPrimitivePatch<PatchType>:: getVolumeType
 (
-    const indexedOctree
-    <
-        treeDataPrimitivePatch
-        <
-            Face,
-            FaceList,
-            PointField,
-            PointType
-        >
-    >& oc,
+    const indexedOctree<treeDataPrimitivePatch<PatchType> >& oc,
     const point& sample
 ) const
 {
@@ -396,16 +341,8 @@ getVolumeType
 
 
 // Check if any point on shape is inside cubeBb.
-template
-<
-    class Face,
-    template<class> class FaceList,
-    class PointField,
-    class PointType
->
-bool
-Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
-overlaps
+template<class PatchType>
+bool Foam::treeDataPrimitivePatch<PatchType>::overlaps
 (
     const label index,
     const treeBoundBox& cubeBb
@@ -462,16 +399,8 @@ overlaps
 
 
 // Check if any point on shape is inside sphere.
-template
-<
-    class Face,
-    template<class> class FaceList,
-    class PointField,
-    class PointType
->
-bool
-Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
-overlaps
+template<class PatchType>
+bool Foam::treeDataPrimitivePatch<PatchType>::overlaps
 (
     const label index,
     const point& centre,
@@ -512,16 +441,8 @@ overlaps
 
 // Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex,
 // nearestPoint.
-template
-<
-    class Face,
-    template<class> class FaceList,
-    class PointField,
-    class PointType
->
-void
-Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
-findNearest
+template<class PatchType>
+void Foam::treeDataPrimitivePatch<PatchType>::findNearest
 (
     const labelUList& indices,
     const point& sample,
@@ -552,16 +473,8 @@ findNearest
 }
 
 
-template
-<
-    class Face,
-    template<class> class FaceList,
-    class PointField,
-    class PointType
->
-bool
-Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
-intersects
+template<class PatchType>
+bool Foam::treeDataPrimitivePatch<PatchType>::intersects
 (
     const label index,
     const point& start,
diff --git a/src/meshTools/indexedOctree/treeDataPrimitivePatch.H b/src/meshTools/indexedOctree/treeDataPrimitivePatch.H
index 44f76d2c75..bbb3c74956 100644
--- a/src/meshTools/indexedOctree/treeDataPrimitivePatch.H
+++ b/src/meshTools/indexedOctree/treeDataPrimitivePatch.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -59,13 +59,7 @@ TemplateName(treeDataPrimitivePatch);
                    Class treeDataPrimitivePatch Declaration
 \*---------------------------------------------------------------------------*/
 
-template
-<
-    class Face,
-    template<class> class FaceList,
-    class PointField,
-    class PointType=point
->
+template<class PatchType>
 class treeDataPrimitivePatch
 :
     public treeDataPrimitivePatchName
@@ -78,7 +72,7 @@ class treeDataPrimitivePatch
     // Private data
 
         //- Underlying geometry
-        const PrimitivePatch<Face, FaceList, PointField, PointType>& patch_;
+        const PatchType& patch_;
 
         //- Whether to precalculate and store face bounding box
         const bool cacheBb_;
@@ -103,7 +97,7 @@ public:
         treeDataPrimitivePatch
         (
             const bool cacheBb,
-            const PrimitivePatch<Face, FaceList, PointField, PointType>&
+            const PatchType&
         );
 
 
@@ -121,8 +115,7 @@ public:
             pointField shapePoints() const;
 
             //- Return access to the underlying patch
-            const PrimitivePatch<Face, FaceList, PointField, PointType>&
-            patch() const
+            const PatchType& patch() const
             {
                 return patch_;
             }
@@ -134,16 +127,7 @@ public:
             //  Only makes sense for closed surfaces.
             label getVolumeType
             (
-                const indexedOctree
-                <
-                    treeDataPrimitivePatch
-                    <
-                        Face,
-                        FaceList,
-                        PointField,
-                        PointType
-                    >
-                >&,
+                const indexedOctree<treeDataPrimitivePatch<PatchType> >&,
                 const point&
             ) const;
 
diff --git a/src/meshTools/searchableSurface/searchableSurface.C b/src/meshTools/searchableSurface/searchableSurface.C
index b586ca191d..decd4c7259 100644
--- a/src/meshTools/searchableSurface/searchableSurface.C
+++ b/src/meshTools/searchableSurface/searchableSurface.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -95,4 +95,21 @@ Foam::searchableSurface::~searchableSurface()
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::searchableSurface::findNearest
+(
+    const pointField& sample,
+    const scalarField& nearestDistSqr,
+    List<pointIndexHit>& info,
+    vectorField& normal,
+    labelList& region
+) const
+{
+    findNearest(sample, nearestDistSqr, info);
+    getNormal(info, normal);
+    getRegion(info, region);
+}
+
+
 // ************************************************************************* //
diff --git a/src/meshTools/searchableSurface/searchableSurface.H b/src/meshTools/searchableSurface/searchableSurface.H
index 535997e900..94713fbf38 100644
--- a/src/meshTools/searchableSurface/searchableSurface.H
+++ b/src/meshTools/searchableSurface/searchableSurface.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -323,6 +323,17 @@ public:
                 List<volumeType>&
             ) const = 0;
 
+            //- Find nearest, normal and region. Can be overridden with
+            //  optimised implementation
+            virtual void findNearest
+            (
+                const pointField& sample,
+                const scalarField& nearestDistSqr,
+                List<pointIndexHit>&,
+                vectorField& normal,
+                labelList& region
+            ) const;
+
 
         // Other
 
diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.C b/src/meshTools/searchableSurface/triSurfaceMesh.C
index 6bbfd22960..86fcd5889b 100644
--- a/src/meshTools/searchableSurface/triSurfaceMesh.C
+++ b/src/meshTools/searchableSurface/triSurfaceMesh.C
@@ -827,6 +827,8 @@ void Foam::triSurfaceMesh::getNormal
     vectorField& normal
 ) const
 {
+    const triSurface& s = static_cast<const triSurface&>(*this);
+
     normal.setSize(info.size());
 
     if (minQuality_ >= 0)
@@ -834,7 +836,6 @@ void Foam::triSurfaceMesh::getNormal
         // Make sure we don't use triangles with low quality since
         // normal is not reliable.
 
-        const triSurface& s = static_cast<const triSurface&>(*this);
         const labelListList& faceFaces = s.faceFaces();
 
         forAll(info, i)
@@ -842,6 +843,7 @@ void Foam::triSurfaceMesh::getNormal
             if (info[i].hit())
             {
                 label faceI = info[i].index();
+                normal[i] = s[faceI].normal(points());
 
                 scalar qual = s[faceI].tri(points()).quality();
 
@@ -861,11 +863,8 @@ void Foam::triSurfaceMesh::getNormal
                         }
                     }
                 }
-                else
-                {
-                    normal[i] = s[faceI].normal(points());
-                }
-                normal[i] /= mag(normal[i]);
+
+                normal[i] /= mag(normal[i]) + VSMALL;
             }
             else
             {
@@ -885,7 +884,7 @@ void Foam::triSurfaceMesh::getNormal
                 //normal[i] = faceNormals()[faceI];
 
                 //- Uncached
-                normal[i] = triSurface::operator[](faceI).normal(points());
+                normal[i] = s[faceI].normal(points());
                 normal[i] /= mag(normal[i]) + VSMALL;
             }
             else