diff --git a/applications/utilities/mesh/manipulation/createPatch/createPatch.C b/applications/utilities/mesh/manipulation/createPatch/createPatch.C
index 01df2c2f28..a993ecf92a 100644
--- a/applications/utilities/mesh/manipulation/createPatch/createPatch.C
+++ b/applications/utilities/mesh/manipulation/createPatch/createPatch.C
@@ -625,10 +625,8 @@ int main(int argc, char *argv[])
                     << " as patch " << destPatchI
                     << " from " << patchDict << endl;
 
-                patchDict.remove("nFaces");
-                patchDict.add("nFaces", 0);
-                patchDict.remove("startFace");
-                patchDict.add("startFace", startFaceI);
+                patchDict.set("nFaces", 0);
+                patchDict.set("startFace", startFaceI);
 
                 // Add an empty patch.
                 allPatches.append
diff --git a/applications/utilities/parallelProcessing/reconstructParMesh/Make/options b/applications/utilities/parallelProcessing/reconstructParMesh/Make/options
index 73f34b0f85..f9b5424edd 100644
--- a/applications/utilities/parallelProcessing/reconstructParMesh/Make/options
+++ b/applications/utilities/parallelProcessing/reconstructParMesh/Make/options
@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
diff --git a/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C b/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C
index e824c20e39..bfd6ae083b 100644
--- a/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C
+++ b/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C
@@ -46,6 +46,8 @@ Description
 #include "mapAddedPolyMesh.H"
 #include "polyMeshAdder.H"
 #include "faceCoupleInfo.H"
+#include "fvMeshAdder.H"
+#include "polyTopoChange.H"
 
 using namespace Foam;
 
@@ -203,6 +205,78 @@ autoPtr<faceCoupleInfo> determineCoupledFaces
 }
 
 
+autoPtr<mapPolyMesh> mergeSharedPoints
+(
+    const scalar mergeDist,
+    polyMesh& mesh,
+    labelListList& pointProcAddressing
+)
+{
+    // Find out which sets of points get merged and create a map from
+    // mesh point to unique point.
+    Map<label> pointToMaster
+    (
+        fvMeshAdder::findSharedPoints
+        (
+            mesh,
+            mergeDist
+        )
+    );
+
+    Info<< "mergeSharedPoints : detected " << pointToMaster.size()
+        << " points that are to be merged." << endl;
+
+    if (returnReduce(pointToMaster.size(), sumOp<label>()) == 0)
+    {
+        return autoPtr<mapPolyMesh>(NULL);
+    }
+
+    polyTopoChange meshMod(mesh);
+
+    fvMeshAdder::mergePoints(mesh, pointToMaster, meshMod);
+
+    // Change the mesh (no inflation). Note: parallel comms allowed.
+    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
+
+    // Update fields. No inflation, parallel sync.
+    mesh.updateMesh(map);
+
+    // pointProcAddressing give indices into the master mesh so adapt them
+    // for changed point numbering.
+
+    // Adapt constructMaps for merged points.
+    forAll(pointProcAddressing, procI)
+    {
+        labelList& constructMap = pointProcAddressing[procI];
+
+        forAll(constructMap, i)
+        {
+            label oldPointI = constructMap[i];
+
+            // New label of point after changeMesh.
+            label newPointI = map().reversePointMap()[oldPointI];
+
+            if (newPointI < -1)
+            {
+                constructMap[i] = -newPointI-2;
+            }
+            else if (newPointI >= 0)
+            {
+                constructMap[i] = newPointI;
+            }
+            else
+            {
+                FatalErrorIn("fvMeshDistribute::mergeSharedPoints()")
+                    << "Problem. oldPointI:" << oldPointI
+                    << " newPointI:" << newPointI << abort(FatalError);
+            }
+        }
+    }
+
+    return map;    
+}
+
+
 int main(int argc, char *argv[])
 {
     argList::noParallel();
@@ -214,7 +288,7 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 
-    Pout<< "This is an experimental tool which tries to merge"
+    Info<< "This is an experimental tool which tries to merge"
         << " individual processor" << nl
         << "meshes back into one master mesh. Use it if the original"
         << " master mesh has" << nl
@@ -246,7 +320,7 @@ int main(int argc, char *argv[])
     }
     scalar writeTol = Foam::pow(10.0, -scalar(IOstream::defaultPrecision()));
 
-    Pout<< "Merge tolerance : " << mergeTol << nl
+    Info<< "Merge tolerance : " << mergeTol << nl
         << "Write tolerance : " << writeTol << endl;
 
     if (runTime.writeFormat() == IOstream::ASCII && mergeTol < writeTol)
@@ -267,11 +341,11 @@ int main(int argc, char *argv[])
 
     if (fullMatch)
     {
-        Pout<< "Doing geometric matching on all boundary faces." << nl << endl;
+        Info<< "Doing geometric matching on all boundary faces." << nl << endl;
     }
     else
     {
-        Pout<< "Doing geometric matching on correct procBoundaries only."
+        Info<< "Doing geometric matching on correct procBoundaries only."
             << nl << "This assumes a correct decomposition." << endl;
     }
 
@@ -292,7 +366,7 @@ int main(int argc, char *argv[])
         nProcs++;
     }
 
-    Pout<< "Found " << nProcs << " processor directories" << nl << endl;
+    Info<< "Found " << nProcs << " processor directories" << nl << endl;
 
 
     // Read all databases.
@@ -300,7 +374,7 @@ int main(int argc, char *argv[])
 
     forAll (databases, procI)
     {
-        Pout<< "Reading database "
+        Info<< "Reading database "
             << args.caseName()/fileName(word("processor") + name(procI))
             << endl;
 
@@ -337,7 +411,7 @@ int main(int argc, char *argv[])
     }
 
     // Set master time
-    Pout<< "Setting master time to " << databases[0].timeName() << nl << endl;
+    Info<< "Setting master time to " << databases[0].timeName() << nl << endl;
     runTime.setTime(databases[0]);
 
 
@@ -373,7 +447,7 @@ int main(int argc, char *argv[])
                 << endl << exit(FatalError);
         }
 
-        Pout<< "Reading points from "
+        Info<< "Reading points from "
             << databases[procI].caseName()
             << " for time = " << databases[procI].timeName()
             << nl << endl;
@@ -403,7 +477,7 @@ int main(int argc, char *argv[])
     }
     const scalar mergeDist = mergeTol*mag(bb.max() - bb.min());
 
-    Pout<< "Overall mesh bounding box  : " << bb << nl
+    Info<< "Overall mesh bounding box  : " << bb << nl
         << "Relative tolerance         : " << mergeTol << nl
         << "Absolute matching distance : " << mergeDist << nl
         << endl;
@@ -422,7 +496,7 @@ int main(int argc, char *argv[])
 
     {
         // Construct empty mesh.
-        Pout<< "Constructing empty mesh to add to." << nl << endl;
+        Info<< "Constructing empty mesh to add to." << nl << endl;
         polyMesh masterMesh
         (
             IOobject
@@ -439,7 +513,7 @@ int main(int argc, char *argv[])
 
         for (label procI = 0; procI < nProcs; procI++)
         {
-            Pout<< "Reading mesh to add from "
+            Info<< "Reading mesh to add from "
                 << databases[procI].caseName()
                 << " for time = " << databases[procI].timeName()
                 << nl << endl;
@@ -475,7 +549,7 @@ int main(int argc, char *argv[])
 
 
             // Add elements to mesh
-            Pout<< "Adding to master mesh" << nl << endl;
+            Info<< "Adding to master mesh" << nl << endl;
 
             autoPtr<mapAddedPolyMesh> map = polyMeshAdder::add
             (
@@ -503,16 +577,19 @@ int main(int argc, char *argv[])
             renumber(map().addedPointMap(), pointProcAddressing[procI]);
             renumber(map().addedPatchMap(), boundaryProcAddressing[procI]);
 
-            Pout<< endl;
+            Info<< endl;
         }
 
+        // See if any points on the mastermesh have become connected
+        // because of connections through processor meshes.
+        mergeSharedPoints(mergeDist, masterMesh, pointProcAddressing);
 
         // Save some properties on the reconstructed mesh
         masterInternalFaces = masterMesh.nInternalFaces();
         masterOwner = masterMesh.faceOwner();
 
 
-        Pout<< "\nWriting merged mesh to "
+        Info<< "\nWriting merged mesh to "
             << runTime.path()/runTime.timeName()
             << nl << endl;
 
@@ -527,12 +604,12 @@ int main(int argc, char *argv[])
 
     // Write the addressing
 
-    Pout<< "Reconstructing the addressing from the processor meshes"
+    Info<< "Reconstructing the addressing from the processor meshes"
         << " to the newly reconstructed mesh" << nl << endl;
 
     forAll(databases, procI)
     {
-        Pout<< "Reading processor " << procI << " mesh from "
+        Info<< "Reading processor " << procI << " mesh from "
             << databases[procI].caseName() << endl;
 
         polyMesh procMesh
@@ -548,7 +625,7 @@ int main(int argc, char *argv[])
 
         // From processor point to reconstructed mesh point
 
-        Pout<< "Writing pointProcAddressing to "
+        Info<< "Writing pointProcAddressing to "
             << databases[procI].caseName()
               /procMesh.facesInstance()
               /polyMesh::meshSubDir
@@ -572,7 +649,7 @@ int main(int argc, char *argv[])
 
         // From processor face to reconstructed mesh face
 
-        Pout<< "Writing faceProcAddressing to "
+        Info<< "Writing faceProcAddressing to "
             << databases[procI].caseName()
               /procMesh.facesInstance()
               /polyMesh::meshSubDir
@@ -635,7 +712,7 @@ int main(int argc, char *argv[])
 
         // From processor cell to reconstructed mesh cell
 
-        Pout<< "Writing cellProcAddressing to "
+        Info<< "Writing cellProcAddressing to "
             << databases[procI].caseName()
               /procMesh.facesInstance()
               /polyMesh::meshSubDir
@@ -660,7 +737,7 @@ int main(int argc, char *argv[])
 
         // From processor patch to reconstructed mesh patch
 
-        Pout<< "Writing boundaryProcAddressing to "
+        Info<< "Writing boundaryProcAddressing to "
             << databases[procI].caseName()
               /procMesh.facesInstance()
               /polyMesh::meshSubDir
@@ -681,10 +758,10 @@ int main(int argc, char *argv[])
             boundaryProcAddressing[procI]
         ).write();
 
-        Pout<< endl;
+        Info<< endl;
     }
 
-    Pout<< "End.\n" << endl;
+    Info<< "End.\n" << endl;
 
     return 0;
 }
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndexTemplates.C b/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndexTemplates.C
index 9fc8ff0360..6f7e548c70 100644
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndexTemplates.C
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndexTemplates.C
@@ -54,9 +54,11 @@ Foam::Field<T> Foam::channelIndex::collapse
 ) const
 {
     // Average and order
+    const Field<T> summedField(regionSum(cellField));
+
     Field<T> regionField
     (
-        regionSum(cellField)
+        summedField
       / regionCount_,
         sortMap_
     );
diff --git a/src/OpenFOAM/containers/Lists/DynamicList/DynamicListI.H b/src/OpenFOAM/containers/Lists/DynamicList/DynamicListI.H
index da6bd13067..9fdfe21bb1 100644
--- a/src/OpenFOAM/containers/Lists/DynamicList/DynamicListI.H
+++ b/src/OpenFOAM/containers/Lists/DynamicList/DynamicListI.H
@@ -79,8 +79,9 @@ inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::setSize
     const label s
 )
 {
-    if (s < List<T>::size())
+    if (s <= List<T>::size())
     {
+        // shrink addressable size, leave allocated size untouched
         List<T>::size() = s;
     }
     else
@@ -100,8 +101,9 @@ inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::setSize
     const T& t
 )
 {
-    if (s < List<T>::size())
+    if (s <= List<T>::size())
     {
+        // shrink addressable size, leave allocated size untouched
         List<T>::size() = s;
     }
     else
diff --git a/src/OpenFOAM/containers/Lists/FixedList/FixedListIO.C b/src/OpenFOAM/containers/Lists/FixedList/FixedListIO.C
index 11ae6f7140..0695eeaefc 100644
--- a/src/OpenFOAM/containers/Lists/FixedList/FixedListIO.C
+++ b/src/OpenFOAM/containers/Lists/FixedList/FixedListIO.C
@@ -77,7 +77,7 @@ Foam::Istream& Foam::operator>>(Foam::Istream& is, FixedList<T, Size>& L)
         }
         else
         {
-            // Putback the openning bracket
+            // Putback the opening bracket
             is.putBack(firstToken);
         }
 
@@ -190,8 +190,9 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const FixedList<T, Size>& L)
 
         if (uniform)
         {
-            // Write size of list and start contents delimiter
-            os << token::BEGIN_BLOCK;
+            // Write size of list (so it is valid dictionary entry) and
+            // start contents delimiter
+            os << L.size() << token::BEGIN_BLOCK;
 
             // Write list contents
             os << L[0];
@@ -201,7 +202,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const FixedList<T, Size>& L)
         }
         else if (Size < 11 && contiguous<T>())
         {
-            // Write size of list and start contents delimiter
+            // Write start of contents delimiter
             os << token::BEGIN_LIST;
 
             // Write list contents
@@ -216,7 +217,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const FixedList<T, Size>& L)
         }
         else
         {
-            // Write size of list and start contents delimiter
+            // Write start of contents delimiter
             os << nl << token::BEGIN_LIST;
 
             // Write list contents
diff --git a/src/OpenFOAM/meshes/meshTools/matchPoints.C b/src/OpenFOAM/meshes/meshTools/matchPoints.C
index 943724c461..89161531ec 100644
--- a/src/OpenFOAM/meshes/meshTools/matchPoints.C
+++ b/src/OpenFOAM/meshes/meshTools/matchPoints.C
@@ -45,9 +45,19 @@ bool Foam::matchPoints
 
     bool fullMatch = true;
 
-    SortableList<scalar> pts0MagSqr(magSqr(pts0 - origin));
+    point compareOrigin = origin;
 
-    SortableList<scalar> pts1MagSqr(magSqr(pts1 - origin));
+    if (origin == point(VGREAT, VGREAT, VGREAT))
+    {
+        if (pts1.size() > 0)
+        {
+            compareOrigin = sum(pts1)/pts1.size();
+        }
+    }
+
+    SortableList<scalar> pts0MagSqr(magSqr(pts0 - compareOrigin));
+
+    SortableList<scalar> pts1MagSqr(magSqr(pts1 - compareOrigin));
 
     forAll(pts0MagSqr, i)
     {
diff --git a/src/OpenFOAM/meshes/meshTools/mergePoints.C b/src/OpenFOAM/meshes/meshTools/mergePoints.C
index cc966c8127..72099978c2 100644
--- a/src/OpenFOAM/meshes/meshTools/mergePoints.C
+++ b/src/OpenFOAM/meshes/meshTools/mergePoints.C
@@ -40,6 +40,16 @@ bool Foam::mergePoints
     const point& origin
 )
 {
+    point compareOrigin = origin;
+
+    if (origin == point(VGREAT, VGREAT, VGREAT))
+    {
+        if (points.size() > 0)
+        {
+            compareOrigin = sum(points)/points.size();
+        }
+    }
+
     // Create a old to new point mapping array
     pointMap.setSize(points.size());
     pointMap = -1;
@@ -56,7 +66,7 @@ bool Foam::mergePoints
     const scalar mergeTolSqr = sqr(mergeTol);
 
     // Sort points by magSqr
-    SortableList<scalar> sortedMagSqr(magSqr(points - origin));
+    SortableList<scalar> sortedMagSqr(magSqr(points - compareOrigin));
 
     bool hasMerged = false;
 
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
index 17ad2db042..7f7f3d2f5a 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
@@ -279,27 +279,13 @@ Foam::layerParameters::layerParameters
 {
     const dictionary& layersDict = dict.subDict("layers");
 
-    forAllConstIter(dictionary, layersDict, iter)
+    forAll(boundaryMesh, patchI)
     {
-        const word& key = iter().keyword();
+        const word& patchName = boundaryMesh[patchI].name();
 
-        if (layersDict.isDict(key))
+        if (layersDict.found(patchName))
         {
-            label patchI = boundaryMesh.findPatchID(key);
-
-            if (patchI == -1)
-            {
-                FatalErrorIn
-                (
-                    "layerParameters::layerParameters"
-                    "(const dictionary&, const polyBoundaryMesh&)"
-                )   << "Specified illegal patch " << key
-                    << " in layer dictionary." << endl
-                    << "Valid patch names are " << boundaryMesh.names()
-                    << exit(FatalError);
-            }
-
-            const dictionary& layerDict = layersDict.subDict(key);
+            const dictionary& layerDict = layersDict.subDict(patchName);
 
             numLayers_[patchI] =
                 readLabel(layerDict.lookup("nSurfaceLayers"));
diff --git a/src/dynamicMesh/fvMeshAdder/fvMeshAdder.H b/src/dynamicMesh/fvMeshAdder/fvMeshAdder.H
index d19ccf0325..95ecfd6519 100644
--- a/src/dynamicMesh/fvMeshAdder/fvMeshAdder.H
+++ b/src/dynamicMesh/fvMeshAdder/fvMeshAdder.H
@@ -40,6 +40,7 @@ SourceFiles
 
 #include "polyMeshAdder.H"
 #include "fvPatchFieldsFwd.H"
+#include "fvsPatchFieldsFwd.H"
 #include "fvPatchFieldMapper.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
index 25746f43a0..c181d2d5c1 100644
--- a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
+++ b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
@@ -586,8 +586,8 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::repatch
                 "fvMeshDistribute::repatch(const labelList&, labelListList&)"
             )   << "reverseFaceMap contains -1 at index:"
                 << index << endl
-                << "This means that the repatch operation was not just a shuffle?"
-                << abort(FatalError);
+                << "This means that the repatch operation was not just"
+                << " a shuffle?" << abort(FatalError);
         }
     }
 
@@ -622,9 +622,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::mergeSharedPoints
         )
     );
 
-    bool merged = pointToMaster.size() > 0;
-
-    if (!returnReduce(merged, orOp<bool>()))
+    if (returnReduce(pointToMaster.size(), sumOp<label>()) == 0)
     {
         return autoPtr<mapPolyMesh>(NULL);
     }
@@ -639,14 +637,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::mergeSharedPoints
     // Update fields. No inflation, parallel sync.
     mesh_.updateMesh(map);
 
-    // Move mesh (since morphing does not do this)
-    if (map().hasMotionPoints())
-    {
-        mesh_.movePoints(map().preMotionPoints());
-    }
-
     // Adapt constructMaps for merged points.
-    // 1.4.1: use reversePointMap < -1 feature.
     forAll(constructPointMap, procI)
     {
         labelList& constructMap = constructPointMap[procI];
@@ -655,16 +646,38 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::mergeSharedPoints
         {
             label oldPointI = constructMap[i];
 
-            // See if merged into other point
-            Map<label>::const_iterator iter = pointToMaster.find(oldPointI);
+            label newPointI = map().reversePointMap()[oldPointI];
 
-            if (iter != pointToMaster.end())
+            if (newPointI < -1)
             {
-                oldPointI = iter();
+                constructMap[i] = -newPointI-2;
+            }
+            else if (newPointI >= 0)
+            {
+                constructMap[i] = newPointI;
+            }
+            else
+            {
+                FatalErrorIn("fvMeshDistribute::mergeSharedPoints()")
+                    << "Problem. oldPointI:" << oldPointI
+                    << " newPointI:" << newPointI << abort(FatalError);
             }
-
-            constructMap[i] = map().reversePointMap()[oldPointI];
         }
+        //- old: use pointToMaster map.
+        //forAll(constructMap, i)
+        //{
+        //    label oldPointI = constructMap[i];
+        //
+        //    // See if merged into other point
+        //    Map<label>::const_iterator iter = pointToMaster.find(oldPointI);
+        //
+        //    if (iter != pointToMaster.end())
+        //    {
+        //        oldPointI = iter();
+        //    }
+        //
+        //    constructMap[i] = map().reversePointMap()[oldPointI];
+        //}
     }
     return map;
 }
diff --git a/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C b/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C
index 9fef59a208..b04bc9101e 100644
--- a/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C
+++ b/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C
@@ -1737,85 +1737,218 @@ Foam::autoPtr<Foam::mapAddedPolyMesh> Foam::polyMeshAdder::add
 Foam::Map<Foam::label> Foam::polyMeshAdder::findSharedPoints
 (
     const polyMesh& mesh,
-    const scalar mergeTol
+    const scalar mergeDist
 )
 {
     const labelList& sharedPointLabels = mesh.globalData().sharedPointLabels();
+    const labelList& sharedPointAddr = mesh.globalData().sharedPointAddr();
 
-    labelList sharedToMerged;
-    pointField mergedPoints;
-    bool hasMerged = Foam::mergePoints
-    (
-        pointField
-        (
-            IndirectList<point>
+    // Because of adding the missing pieces e.g. when redistributing a mesh
+    // it can be that there are multiple points on the same processor that
+    // refer to the same shared point.
+
+    // Invert point-to-shared addressing
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    Map<labelList> sharedToMesh(sharedPointLabels.size());
+
+    label nMultiple = 0;
+
+    forAll(sharedPointLabels, i)
+    {
+        label pointI = sharedPointLabels[i];
+
+        label sharedI = sharedPointAddr[i];
+
+        Map<labelList>::iterator iter = sharedToMesh.find(sharedI);
+
+        if (iter != sharedToMesh.end())
+        {
+            // sharedI already used by other point. Add this one.
+
+            nMultiple++;
+
+            labelList& connectedPointLabels = iter();
+
+            label sz = connectedPointLabels.size();
+
+            // Check just to make sure.
+            if (findIndex(connectedPointLabels, pointI) != -1)
+            {
+                FatalErrorIn("polyMeshAdder::findSharedPoints(..)")
+                    << "Duplicate point in sharedPoint addressing." << endl
+                    << "When trying to add point " << pointI << " on shared "
+                    << sharedI  << " with connected points "
+                    << connectedPointLabels
+                    << abort(FatalError);
+            }
+
+            connectedPointLabels.setSize(sz+1);
+            connectedPointLabels[sz] = pointI;
+        }
+        else
+        {
+            sharedToMesh.insert(sharedI, labelList(1, pointI));
+        }
+    }
+
+    
+    // Assign single master for every shared with multiple geometric points
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    Map<label> pointToMaster(nMultiple);
+
+    forAllConstIter(Map<labelList>, sharedToMesh, iter)
+    {
+        const labelList& connectedPointLabels = iter();
+
+        //Pout<< "For shared:" << iter.key()
+        //    << " found points:" << connectedPointLabels
+        //    << " at coords:"
+        //    <<  pointField(mesh.points(), connectedPointLabels) << endl;
+
+        if (connectedPointLabels.size() > 1)
+        {
+            const pointField connectedPoints
             (
                 mesh.points(),
-                sharedPointLabels
-            )()
-        ),
-        mergeTol,
-        false,
-        sharedToMerged,
-        mergedPoints
-    );
+                connectedPointLabels
+            );
 
-    // Find out which sets of points get merged and create a map from
-    // mesh point to unique point.
-
-    Map<label> pointToMaster(10*sharedToMerged.size());
-
-    if (hasMerged)
-    {
-        labelListList mergeSets
-        (
-            invertOneToMany
+            labelList toMergedPoints;
+            pointField mergedPoints;
+            bool hasMerged = Foam::mergePoints
             (
-                sharedToMerged.size(),
-                sharedToMerged
-            )
-        );
+                connectedPoints,
+                mergeDist,
+                false,
+                toMergedPoints,
+                mergedPoints
+            );
 
-        label nMergeSets = 0;
-
-        forAll(mergeSets, setI)
-        {
-            const labelList& mergeSet = mergeSets[setI];
-
-            if (mergeSet.size() > 1)
+            if (hasMerged)
             {
-                // Take as master the shared point with the lowest mesh
-                // point label. (rather arbitrarily - could use max or any other
-                // one of the points)
+                // Invert toMergedPoints
+                const labelListList mergeSets
+                (
+                    invertOneToMany
+                    (
+                        mergedPoints.size(),
+                        toMergedPoints
+                    )
+                );
 
-                nMergeSets++;
-
-                label masterI = labelMax;
-
-                forAll(mergeSet, i)
+                // Find master for valid merges
+                forAll(mergeSets, setI)
                 {
-                    label sharedI = mergeSet[i];
+                    const labelList& mergeSet = mergeSets[setI];
 
-                    masterI = min(masterI, sharedPointLabels[sharedI]);
-                }
+                    if (mergeSet.size() > 1)
+                    {
+                        // Pick lowest numbered point
+                        label masterPointI = labelMax;
 
-                forAll(mergeSet, i)
-                {
-                    label sharedI = mergeSet[i];
+                        forAll(mergeSet, i)
+                        {
+                            label pointI = connectedPointLabels[mergeSet[i]];
 
-                    pointToMaster.insert(sharedPointLabels[sharedI], masterI);
+                            masterPointI = min(masterPointI, pointI);
+                        }
+
+                        forAll(mergeSet, i)
+                        {
+                            label pointI = connectedPointLabels[mergeSet[i]];
+
+                            //Pout<< "Merging point " << pointI
+                            //    << " at " << mesh.points()[pointI]
+                            //    << " into master point "
+                            //    << masterPointI
+                            //    << " at " << mesh.points()[masterPointI]
+                            //    << endl;
+
+                            pointToMaster.insert(pointI, masterPointI);
+                        }
+                    }
                 }
             }
         }
-
-        //if (debug)
-        //{
-        //    Pout<< "polyMeshAdder : merging:"
-        //        << pointToMaster.size() << " into " << nMergeSets << " sets."
-        //        << endl;
-        //}
     }
 
+    //- Old: geometric merging. Causes problems for two close shared points.
+    //labelList sharedToMerged;
+    //pointField mergedPoints;
+    //bool hasMerged = Foam::mergePoints
+    //(
+    //    pointField
+    //    (
+    //        IndirectList<point>
+    //        (
+    //            mesh.points(),
+    //            sharedPointLabels
+    //        )()
+    //    ),
+    //    mergeDist,
+    //    false,
+    //    sharedToMerged,
+    //    mergedPoints
+    //);
+    //
+    //// Find out which sets of points get merged and create a map from
+    //// mesh point to unique point.
+    //
+    //Map<label> pointToMaster(10*sharedToMerged.size());
+    //
+    //if (hasMerged)
+    //{
+    //    labelListList mergeSets
+    //    (
+    //        invertOneToMany
+    //        (
+    //            sharedToMerged.size(),
+    //            sharedToMerged
+    //        )
+    //    );
+    //
+    //    label nMergeSets = 0;
+    //
+    //    forAll(mergeSets, setI)
+    //    {
+    //        const labelList& mergeSet = mergeSets[setI];
+    //
+    //        if (mergeSet.size() > 1)
+    //        {
+    //            // Take as master the shared point with the lowest mesh
+    //            // point label. (rather arbitrarily - could use max or
+    //            // any other one of the points)
+    //
+    //            nMergeSets++;
+    //
+    //            label masterI = labelMax;
+    //
+    //            forAll(mergeSet, i)
+    //            {
+    //                label sharedI = mergeSet[i];
+    //
+    //                masterI = min(masterI, sharedPointLabels[sharedI]);
+    //            }
+    //
+    //            forAll(mergeSet, i)
+    //            {
+    //                label sharedI = mergeSet[i];
+    //
+    //                pointToMaster.insert(sharedPointLabels[sharedI], masterI);
+    //            }
+    //        }
+    //    }
+    //
+    //    //if (debug)
+    //    //{
+    //    //    Pout<< "polyMeshAdder : merging:"
+    //    //        << pointToMaster.size() << " into " << nMergeSets
+    //    //        << " sets." << endl;
+    //    //}
+    //}
+
     return pointToMaster;
 }
 
@@ -1836,8 +1969,7 @@ void Foam::polyMeshAdder::mergePoints
         {
             if (iter() != pointI)
             {
-                //1.4.1: meshMod.removePoint(pointI, iter());
-                meshMod.setAction(polyRemovePoint(pointI));
+                meshMod.removePoint(pointI, iter());
             }
         }
     }
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
index 08e9566e48..4bb6e3368a 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
@@ -2617,6 +2617,13 @@ void Foam::polyTopoChange::removePoint
             << abort(FatalError);
     }
 
+    if (pointI == mergePointI)
+    {
+        FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
+            << "Cannot remove/merge point " << pointI << " onto itself."
+            << abort(FatalError);
+    }
+
     points_[pointI] = greatPoint;
     pointMap_[pointI] = -1;
     if (mergePointI >= 0)
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/refinementHistory.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/refinementHistory.C
index 15d28b5b60..45f1cb9e2a 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/refinementHistory.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/refinementHistory.C
@@ -410,8 +410,14 @@ Foam::refinementHistory::refinementHistory
     }
     else
     {
-        splitCells_.setSize(nCells, splitCell8());
-        visibleCells_ = identity(nCells);
+        visibleCells_.setSize(nCells);
+        splitCells_.setSize(nCells);
+
+        for (label cellI = 0; cellI < nCells; cellI++)
+        {
+            visibleCells_[cellI] = cellI;
+            splitCells_.append(splitCell8());
+        }
     }
 
     // Check indices.
diff --git a/src/triSurface/triSurface/interfaces/STL/readSTLASCII.L b/src/triSurface/triSurface/interfaces/STL/readSTLASCII.L
index e0652edd8c..b5f5c81358 100644
--- a/src/triSurface/triSurface/interfaces/STL/readSTLASCII.L
+++ b/src/triSurface/triSurface/interfaces/STL/readSTLASCII.L
@@ -459,6 +459,7 @@ bool triSurface::readSTLASCII(const fileName& STLfileName)
         patches_[iter()].name() = iter.key();
     }
 
+    // Fill in the missing information in the patches
     setDefaultPatches();
 
     return true;
diff --git a/src/triSurface/triSurface/triSurface.C b/src/triSurface/triSurface/triSurface.C
index a9f7b53e66..548c853890 100644
--- a/src/triSurface/triSurface/triSurface.C
+++ b/src/triSurface/triSurface/triSurface.C
@@ -631,15 +631,16 @@ surfacePatchList triSurface::calcPatches(labelList& faceMap) const
 
     faceMap = sortedRegion.indices();
 
-    // Compact regions
-
-    // Get last region
-    //label maxRegion = 0;                  // for compacted regions
+    // Extend regions
     label maxRegion = patches_.size()-1;    // for non-compacted regions
 
     if (faceMap.size() > 0)
     {
-        maxRegion = operator[](faceMap[faceMap.size() - 1]).region();
+        maxRegion = max
+        (
+            maxRegion,
+            operator[](faceMap[faceMap.size() - 1]).region()
+        );
     }
 
     // Get new region list