diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C
index 5f04825a94..a3d68f55ca 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C
@@ -38,7 +38,8 @@ Foam::scalar Foam::compressibleCourantNo
     scalar CoNum = 0.0;
     scalar meanCoNum = 0.0;
 
-    if (mesh.nInternalFaces())
+    //- Can have fluid domains with 0 cells so do not test.
+    //if (mesh.nInternalFaces())
     {
         surfaceScalarField SfUfbyDelta =
             mesh.surfaceInterpolation::deltaCoeffs()
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index 58f324160e..cc975eb970 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@@ -131,6 +131,12 @@ castellatedMeshControls
                     level (3 3);
                 }
             }
+
+            // Optional angle to detect small-large cell situation perpendicular
+            // to the surface. Is the angle of face w.r.t the local surface
+            // normal. Use on flat(ish) surfaces only. Otherwise
+            // leave out or set to negative number.
+            //perpendicularAngle 10;
         }
     }
 
@@ -317,9 +323,10 @@ meshQualityControls
     minTriangleTwist -1;
 
     //- if >0 : preserve single cells with all points on the surface if the
-    //  resulting volume after snapping is larger than minVolFraction times old
-    //  volume. If <0 : delete always.
-    minVolFraction 0.1;
+    //  resulting volume after snapping (by approximation) is larger than
+    //  minVolCollapseRatio times old volume (i.e. not collapsed to flat cell).
+    //  If <0 : delete always.
+    //minVolCollapseRatio 0.5;
 
 
     // Advanced
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
index de69375595..14f5704836 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
@@ -102,6 +102,67 @@ Foam::label Foam::checkTopology
         }
     }
 
+    if (allTopology)
+    {
+        labelList nInternalFaces(mesh.nCells(), 0);
+
+        for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+        {
+            nInternalFaces[mesh.faceOwner()[faceI]]++;
+            nInternalFaces[mesh.faceNeighbour()[faceI]]++;
+        }
+        const polyBoundaryMesh& patches = mesh.boundaryMesh();
+        forAll(patches, patchI)
+        {
+            if (patches[patchI].coupled())
+            {
+                const unallocLabelList& owners = patches[patchI].faceCells();
+
+                forAll(owners, i)
+                {
+                    nInternalFaces[owners[i]]++;
+                }
+            }
+        }
+
+        faceSet oneCells(mesh, "oneInternalFaceCells", mesh.nCells()/100);
+        faceSet twoCells(mesh, "twoInternalFacesCells", mesh.nCells()/100);
+
+        forAll(nInternalFaces, cellI)
+        {
+            if (nInternalFaces[cellI] <= 1)
+            {
+                oneCells.insert(cellI);
+            }
+            else if (nInternalFaces[cellI] == 2)
+            {
+                twoCells.insert(cellI);
+            }
+        }
+
+        label nOneCells = returnReduce(oneCells.size(), sumOp<label>());
+
+        if (nOneCells > 0)
+        {
+            Info<< "  <<Writing " << nOneCells
+                << " cells with with single non-boundary face to set "
+                << oneCells.name()
+                << endl;
+            oneCells.write();
+        }
+
+        label nTwoCells = returnReduce(twoCells.size(), sumOp<label>());
+
+        if (nTwoCells > 0)
+        {
+            Info<< "  <<Writing " << nTwoCells
+                << " cells with with single non-boundary face to set "
+                << twoCells.name()
+                << endl;
+            twoCells.write();
+        }
+    }
+
     {
         regionSplit rs(mesh);
 
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.H b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.H
index 985f6a3d17..927415c005 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.H
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.H
@@ -60,10 +60,11 @@ SourceFiles
 
 #include "className.H"
 #include "fileName.H"
-#include "volPointInterpolation.H"
 #include "stringList.H"
 #include "wordList.H"
 #include "primitivePatch.H"
+#include "PrimitivePatchInterpolation.H"
+#include "volPointInterpolation.H"
 
 // * * * * * * * * * * * * * Forward Declarations  * * * * * * * * * * * * * //
 
@@ -491,7 +492,6 @@ class vtkPV3Foam
             void convertVolFields
             (
                 const fvMesh&,
-                const volPointInterpolation&,
                 const PtrList<PrimitivePatchInterpolation<primitivePatch> >&,
                 const IOobjectList&,
                 vtkMultiBlockDataSet* output
@@ -502,7 +502,7 @@ class vtkPV3Foam
             void convertVolFieldBlock
             (
                 const GeometricField<Type, fvPatchField, volMesh>&,
-                const GeometricField<Type, pointPatchField, pointMesh>&,
+                autoPtr<GeometricField<Type, pointPatchField, pointMesh> >&,
                 vtkMultiBlockDataSet* output,
                 const partInfo& selector,
                 const List<polyDecomp>& decompLst
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamFields.C b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamFields.C
index 4da5f7de7e..9708df9ac1 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamFields.C
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamFields.C
@@ -106,8 +106,6 @@ void Foam::vtkPV3Foam::convertVolFields
         printMemory();
     }
 
-    // Construct interpolation on the raw mesh
-    volPointInterpolation pInterp(mesh);
 
     PtrList<PrimitivePatchInterpolation<primitivePatch> >
         ppInterpList(mesh.boundaryMesh().size());
@@ -127,23 +125,23 @@ void Foam::vtkPV3Foam::convertVolFields
 
     convertVolFields<scalar>
     (
-        mesh, pInterp, ppInterpList, objects, output
+        mesh, ppInterpList, objects, output
     );
     convertVolFields<vector>
     (
-        mesh, pInterp, ppInterpList, objects, output
+        mesh, ppInterpList, objects, output
     );
     convertVolFields<sphericalTensor>
     (
-        mesh, pInterp, ppInterpList, objects, output
+        mesh, ppInterpList, objects, output
     );
     convertVolFields<symmTensor>
     (
-        mesh, pInterp, ppInterpList, objects, output
+        mesh, ppInterpList, objects, output
     );
     convertVolFields<tensor>
     (
-        mesh, pInterp, ppInterpList, objects, output
+        mesh, ppInterpList, objects, output
     );
 
     if (debug)
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamVolFields.H b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamVolFields.H
index 6c8d32f470..a7cfc94e2c 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamVolFields.H
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamVolFields.H
@@ -34,6 +34,8 @@ InClass
 #include "emptyFvPatchField.H"
 #include "wallPolyPatch.H"
 #include "faceSet.H"
+#include "volPointInterpolation.H"
+
 #include "vtkPV3FoamFaceField.H"
 #include "vtkPV3FoamPatchField.H"
 
@@ -43,7 +45,6 @@ template<class Type>
 void Foam::vtkPV3Foam::convertVolFields
 (
     const fvMesh& mesh,
-    const volPointInterpolation& pInterp,
     const PtrList<PrimitivePatchInterpolation<primitivePatch> >& ppInterpList,
     const IOobjectList& objects,
     vtkMultiBlockDataSet* output
@@ -63,23 +64,22 @@ void Foam::vtkPV3Foam::convertVolFields
             continue;
         }
 
+        // Load field
         GeometricField<Type, fvPatchField, volMesh> tf
         (
             *iter(),
             mesh
         );
 
-        tmp<GeometricField<Type, pointPatchField, pointMesh> > tptf
-        (
-            pInterp.interpolate(tf)
-        );
+        // Interpolated field (demand driven)
+        autoPtr<GeometricField<Type, pointPatchField, pointMesh> > ptfPtr;
 
 
         // Convert activated internalMesh regions
         convertVolFieldBlock
         (
             tf,
-            tptf(),
+            ptfPtr,
             output,
             partInfoVolume_,
             regionPolyDecomp_
@@ -89,7 +89,7 @@ void Foam::vtkPV3Foam::convertVolFields
         convertVolFieldBlock
         (
             tf,
-            tptf(),
+            ptfPtr,
             output,
             partInfoCellZones_,
             zonePolyDecomp_
@@ -99,7 +99,7 @@ void Foam::vtkPV3Foam::convertVolFields
         convertVolFieldBlock
         (
             tf,
-            tptf(),
+            ptfPtr,
             output,
             partInfoCellSets_,
             csetPolyDecomp_
@@ -258,46 +258,60 @@ void Foam::vtkPV3Foam::convertVolFields
     }
 }
 
+
 template<class Type>
 void Foam::vtkPV3Foam::convertVolFieldBlock
 (
     const GeometricField<Type, fvPatchField, volMesh>& tf,
-    const GeometricField<Type, pointPatchField, pointMesh>& ptf,
+    autoPtr<GeometricField<Type, pointPatchField, pointMesh> >& ptfPtr,
     vtkMultiBlockDataSet* output,
     const partInfo& selector,
     const List<polyDecomp>& decompLst
 )
 {
-   for (int partId = selector.start(); partId < selector.end(); ++partId)
-   {
-       const label datasetNo = partDataset_[partId];
+    for (int partId = selector.start(); partId < selector.end(); ++partId)
+    {
+        const label datasetNo = partDataset_[partId];
 
-       if (datasetNo >= 0 && partStatus_[partId])
-       {
-           convertVolField
-           (
-               tf,
-               output,
-               selector,
-               datasetNo,
-               decompLst[datasetNo]
-           );
+        if (datasetNo >= 0 && partStatus_[partId])
+        {
+            convertVolField
+            (
+                tf,
+                output,
+                selector,
+                datasetNo,
+                decompLst[datasetNo]
+            );
 
-           convertPointField
-           (
-               ptf,
-               tf,
-               output,
-               selector,
-               datasetNo,
-               decompLst[datasetNo]
-           );
-       }
-   }
+            if (!ptfPtr.valid())
+            {
+                if (debug)
+                {
+                    Info<< "convertVolFieldBlock interpolating:" << tf.name()
+                        << endl;
+                }
+
+                ptfPtr.reset
+                (
+                    volPointInterpolation::New(tf.mesh()).interpolate(tf).ptr()
+                );
+            }
+
+            convertPointField
+            (
+                ptfPtr(),
+                tf,
+                output,
+                selector,
+                datasetNo,
+                decompLst[datasetNo]
+            );
+        }
+    }
 }
 
 
-
 template<class Type>
 void Foam::vtkPV3Foam::convertVolField
 (
diff --git a/applications/utilities/surface/surfaceConvert/surfaceConvert.C b/applications/utilities/surface/surfaceConvert/surfaceConvert.C
index 884cd360ea..311dab28f6 100644
--- a/applications/utilities/surface/surfaceConvert/surfaceConvert.C
+++ b/applications/utilities/surface/surfaceConvert/surfaceConvert.C
@@ -115,7 +115,7 @@ int main(int argc, char *argv[])
     }
 
     Info<< "writing " << exportName;
-    if (scaleFactor > 0)
+    if (scaleFactor <= 0)
     {
         Info<< " without scaling" << endl;
     }
diff --git a/src/OpenFOAM/containers/Lists/FixedList/FixedListI.H b/src/OpenFOAM/containers/Lists/FixedList/FixedListI.H
index 564742d491..c3b09850f9 100644
--- a/src/OpenFOAM/containers/Lists/FixedList/FixedListI.H
+++ b/src/OpenFOAM/containers/Lists/FixedList/FixedListI.H
@@ -162,7 +162,7 @@ inline void Foam::FixedList<T, Size>::checkIndex(const label i) const
             << "attempt to access element from zero-sized list"
             << abort(FatalError);
     }
-    else if (i < 0 || i >= Size)
+    else if (i < 0 || i >= label(Size))
     {
         FatalErrorIn("FixedList<T, Size>::checkIndex(const label)")
             << "index " << i << " out of range 0 ... " << (Size-1)
diff --git a/src/OpenFOAM/db/objectRegistry/objectRegistry.C b/src/OpenFOAM/db/objectRegistry/objectRegistry.C
index 1038ecffa1..55cf75fb0f 100644
--- a/src/OpenFOAM/db/objectRegistry/objectRegistry.C
+++ b/src/OpenFOAM/db/objectRegistry/objectRegistry.C
@@ -55,7 +55,8 @@ Foam::objectRegistry::objectRegistry
     HashTable<regIOobject*>(nIoObjects),
     time_(t),
     parent_(t),
-    dbDir_(name())
+    dbDir_(name()),
+    event_(1)
 {}
 
 
@@ -69,7 +70,8 @@ Foam::objectRegistry::objectRegistry
     HashTable<regIOobject*>(nIoObjects),
     time_(io.time()),
     parent_(io.db()),
-    dbDir_(parent_.dbDir()/local()/name())
+    dbDir_(parent_.dbDir()/local()/name()),
+    event_(1)
 {
     writeOpt() = IOobject::AUTO_WRITE;
 }
@@ -135,6 +137,42 @@ const Foam::objectRegistry& Foam::objectRegistry::subRegistry
 }
 
 
+Foam::label Foam::objectRegistry::getEvent() const
+{
+    label curEvent = event_++;
+
+    if (event_ == labelMax)
+    {
+        WarningIn("objectRegistry::getEvent() const")
+            << "Event counter has overflowed. Resetting counter on all"
+            << " dependent objects." << endl
+            << "This might cause extra evaluations." << endl;
+
+        // Reset event counter
+        curEvent = 1;
+        event_ = 2;
+
+        for (const_iterator iter = begin(); iter != end(); ++iter)
+        {
+            const regIOobject& io = *iter();
+
+            if (objectRegistry::debug)
+            {
+                Pout<< "objectRegistry::getEvent() : "
+                    << "resetting count on " << iter.key() << endl;
+            }
+
+            if (io.eventNo() != 0)
+            {
+                const_cast<regIOobject&>(io).eventNo() = curEvent;
+            }
+        }
+    }
+
+    return curEvent;
+}
+
+
 bool Foam::objectRegistry::checkIn(regIOobject& io) const
 {
     if (objectRegistry::debug)
diff --git a/src/OpenFOAM/db/objectRegistry/objectRegistry.H b/src/OpenFOAM/db/objectRegistry/objectRegistry.H
index 6c701e1d85..061f452c0b 100644
--- a/src/OpenFOAM/db/objectRegistry/objectRegistry.H
+++ b/src/OpenFOAM/db/objectRegistry/objectRegistry.H
@@ -64,6 +64,9 @@ class objectRegistry
         //- Local directory path of this objectRegistry relative to time
         fileName dbDir_;
 
+        //- Current event
+        mutable label event_;
+
 
     // Private Member Functions
 
@@ -151,6 +154,9 @@ public:
             template<class Type>
             const Type& lookupObject(const word& name) const;
 
+            //- Return new event number.
+            label getEvent() const;
+
 
         // Edit
 
diff --git a/src/OpenFOAM/db/regIOobject/regIOobject.C b/src/OpenFOAM/db/regIOobject/regIOobject.C
index 1d4d407955..fde596e89f 100644
--- a/src/OpenFOAM/db/regIOobject/regIOobject.C
+++ b/src/OpenFOAM/db/regIOobject/regIOobject.C
@@ -47,6 +47,7 @@ Foam::regIOobject::regIOobject(const IOobject& io)
     registered_(false),
     ownedByRegistry_(false),
     lastModified_(0),
+    eventNo_(db().getEvent()),
     isPtr_(NULL)
 {
     // Register with objectRegistry if requested
@@ -64,6 +65,7 @@ Foam::regIOobject::regIOobject(const regIOobject& rio)
     registered_(false),
     ownedByRegistry_(false),
     lastModified_(rio.lastModified_),
+    eventNo_(db().getEvent()),
     isPtr_(NULL)
 {
     // Do not register copy with objectRegistry
@@ -78,6 +80,7 @@ Foam::regIOobject::regIOobject(const regIOobject& rio, bool registerCopy)
     registered_(false),
     ownedByRegistry_(false),
     lastModified_(rio.lastModified_),
+    eventNo_(db().getEvent()),
     isPtr_(NULL)
 {
     if (registerCopy && rio.registered_)
@@ -153,6 +156,91 @@ bool Foam::regIOobject::checkOut()
 }
 
 
+bool Foam::regIOobject::uptodate(const word& a) const
+{
+    if (db().lookupObject<regIOobject>(a).eventNo() >= eventNo_)
+    {
+        return false;
+    }
+    else
+    {
+        return true;
+    }
+}
+
+
+bool Foam::regIOobject::uptodate(const word& a, const word& b) const
+{
+    if
+    (
+        db().lookupObject<regIOobject>(a).eventNo() >= eventNo_
+     || db().lookupObject<regIOobject>(b).eventNo() >= eventNo_
+    )
+    {
+        return false;
+    }
+    else
+    {
+        return true;
+    }
+}
+
+
+bool Foam::regIOobject::uptodate
+(
+    const word& a,
+    const word& b,
+    const word& c
+) const
+{
+    if
+    (
+        db().lookupObject<regIOobject>(a).eventNo() >= eventNo_
+     || db().lookupObject<regIOobject>(b).eventNo() >= eventNo_
+     || db().lookupObject<regIOobject>(c).eventNo() >= eventNo_
+    )
+    {
+        return false;
+    }
+    else
+    {
+        return true;
+    }
+}
+
+
+bool Foam::regIOobject::uptodate
+(
+    const word& a,
+    const word& b,
+    const word& c,
+    const word& d
+) const
+{
+    if
+    (
+        db().lookupObject<regIOobject>(a).eventNo() >= eventNo_
+     || db().lookupObject<regIOobject>(b).eventNo() >= eventNo_
+     || db().lookupObject<regIOobject>(c).eventNo() >= eventNo_
+     || db().lookupObject<regIOobject>(d).eventNo() >= eventNo_
+    )
+    {
+        return false;
+    }
+    else
+    {
+        return true;
+    }
+}
+
+
+//- Flag me as uptodate
+void Foam::regIOobject::setUptodate()
+{
+    eventNo_ = db().getEvent();
+}
+
+
 // Rename object and re-register with objectRegistry under new name
 void Foam::regIOobject::rename(const word& newName)
 {
diff --git a/src/OpenFOAM/db/regIOobject/regIOobject.H b/src/OpenFOAM/db/regIOobject/regIOobject.H
index 6b50fd7842..9c6ae89065 100644
--- a/src/OpenFOAM/db/regIOobject/regIOobject.H
+++ b/src/OpenFOAM/db/regIOobject/regIOobject.H
@@ -45,6 +45,7 @@ SourceFiles
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+
 namespace Foam
 {
 
@@ -70,6 +71,9 @@ private:
         //- Time of last modification
         mutable time_t lastModified_;
 
+        //- eventNo of last update
+        label eventNo_;
+
         //- Istream for reading
         Istream* isPtr_;
 
@@ -141,6 +145,24 @@ public:
             inline void release();
 
 
+        // Dependency checking
+
+            //- Event number at last update.
+            inline label eventNo() const;
+
+            //- Event number at last update.
+            inline label& eventNo();
+
+            //- Am I uptodate with respect to other regIOobjects
+            bool uptodate(const word&) const;
+            bool uptodate(const word&, const word&) const;
+            bool uptodate(const word&, const word&, const word&) const;
+            bool uptodate(const word&, const word&, const word&, const word&)
+             const;
+
+            //- Flag me as uptodate
+            void setUptodate();
+
         // Edit
 
             //- Rename
diff --git a/src/OpenFOAM/db/regIOobject/regIOobjectI.H b/src/OpenFOAM/db/regIOobject/regIOobjectI.H
index 467127e0aa..033d3d50a0 100644
--- a/src/OpenFOAM/db/regIOobject/regIOobjectI.H
+++ b/src/OpenFOAM/db/regIOobject/regIOobjectI.H
@@ -80,4 +80,15 @@ void Foam::regIOobject::release()
 }
 
 
+Foam::label Foam::regIOobject::eventNo() const
+{
+    return eventNo_;
+}
+
+Foam::label& Foam::regIOobject::eventNo()
+{
+    return eventNo_;
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.C b/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.C
index 1bea3b3a53..3043b103aa 100644
--- a/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.C
+++ b/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.C
@@ -97,6 +97,68 @@ slicedBoundaryField
 }
 
 
+template
+<
+    class Type,
+    template<class> class PatchField,
+    template<class> class SlicedPatchField,
+    class GeoMesh
+>
+Foam::tmp<Foam::FieldField<PatchField, Type> >
+Foam::SlicedGeometricField<Type, PatchField, SlicedPatchField, GeoMesh>::
+slicedBoundaryField
+(
+    const Mesh& mesh,
+    const FieldField<PatchField, Type>& bField,
+    const bool preserveCouples
+)
+{
+    tmp<FieldField<PatchField, Type> > tbf
+    (
+        new FieldField<PatchField, Type>(mesh.boundary().size())
+    );
+
+    FieldField<PatchField, Type>& bf = tbf();
+
+    forAll (mesh.boundary(), patchi)
+    {
+        if (preserveCouples && mesh.boundary()[patchi].coupled())
+        {
+            // For coupled patched construct the correct patch field type
+            bf.set
+            (
+                patchi,
+                PatchField<Type>::New
+                (
+                    mesh.boundary()[patchi].type(),
+                    mesh.boundary()[patchi],
+                    *this
+                )
+            );
+
+            // Assign field
+            bf[patchi] == bField[patchi];
+        }
+        else
+        {
+            // Create unallocated copy of patch field
+            bf.set
+            (
+                patchi,
+                new SlicedPatchField<Type>
+                (
+                    mesh.boundary()[patchi],
+                    DimensionedField<Type, GeoMesh>::null()
+                )
+            );
+            bf[patchi].UList<Type>::operator=(bField[patchi]);
+        }
+    }
+
+    return tbf;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template
@@ -204,6 +266,64 @@ SlicedGeometricField
 }
 
 
+template
+<
+    class Type,
+    template<class> class PatchField,
+    template<class> class SlicedPatchField,
+    class GeoMesh
+>
+Foam::SlicedGeometricField<Type, PatchField, SlicedPatchField, GeoMesh>::
+SlicedGeometricField
+(
+    const IOobject& io,
+    const GeometricField<Type, PatchField, GeoMesh>& gf,
+    const bool preserveCouples
+)
+:
+    GeometricField<Type, PatchField, GeoMesh>
+    (
+        io,
+        gf.mesh(),
+        gf.dimensions(),
+        Field<Type>(),
+        slicedBoundaryField(gf.mesh(), gf.boundaryField(), preserveCouples)
+    )
+{
+    // Set the internalField to the supplied internal field
+    UList<Type>::operator=(gf.internalField());
+
+    correctBoundaryConditions();
+}
+
+
+template
+<
+    class Type,
+    template<class> class PatchField,
+    template<class> class SlicedPatchField,
+    class GeoMesh
+>
+Foam::SlicedGeometricField<Type, PatchField, SlicedPatchField, GeoMesh>::
+SlicedGeometricField
+(
+    const SlicedGeometricField<Type, PatchField, SlicedPatchField, GeoMesh>& gf
+)
+:
+    GeometricField<Type, PatchField, GeoMesh>
+    (
+        gf,
+        gf.mesh(),
+        gf.dimensions(),
+        Field<Type>(),
+        slicedBoundaryField(gf.mesh(), gf.boundaryField(), true)
+    )
+{
+    // Set the internalField to the supplied internal field
+    UList<Type>::operator=(gf.internalField());
+}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 template
diff --git a/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.H b/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.H
index 10f75917b5..81c41fd505 100644
--- a/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.H
+++ b/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.H
@@ -87,8 +87,17 @@ private:
             const bool preserveCouples
         );
 
-        //- Disallow default bitwise copy construct
-        SlicedGeometricField(const SlicedGeometricField&);
+        //- Slice the given field and a create a PtrList of SlicedPatchField
+        //  from which the boundary field is built
+        tmp<FieldField<PatchField, Type> >  slicedBoundaryField
+        (
+            const Mesh& mesh,
+            const FieldField<PatchField, Type>& bField,
+            const bool preserveCouples
+        );
+
+        ////- Disallow default bitwise copy construct
+        //SlicedGeometricField(const SlicedGeometricField&);
 
         //- Disallow default bitwise assignment
         void operator=(const SlicedGeometricField&);
@@ -128,6 +137,27 @@ public:
             const bool preserveCouples=true
         );
 
+        //- Construct from GeometricField. Reuses full internal and
+        //  patch fields except on couples (preserveCouples=true).
+        SlicedGeometricField
+        (
+            const IOobject&,
+            const GeometricField<Type, PatchField, GeoMesh>&,
+            const bool preserveCouples=true
+        );
+
+        //- Construct as copy
+        SlicedGeometricField
+        (
+            const SlicedGeometricField
+            <
+                Type,
+                PatchField,
+                SlicedPatchField,
+                GeoMesh
+            >&
+        );
+
 
     // Destructor
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshIO.C b/src/OpenFOAM/meshes/polyMesh/polyMeshIO.C
index 4032a7b10b..b8b133ed01 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshIO.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshIO.C
@@ -268,7 +268,12 @@ Foam::polyMesh::readUpdateState Foam::polyMesh::readUpdate()
         // Calculate the geometry for the patches (transformation tensors etc.)
         boundary_.calcGeometry();
 
+        // Derived info
+        bounds_ = boundBox(points_);
+        geometricD_ = Vector<label>::zero;
+        solutionD_ = Vector<label>::zero;
 
+        // Zones
         pointZoneMesh newPointZones
         (
             IOobject
@@ -418,6 +423,13 @@ Foam::polyMesh::readUpdateState Foam::polyMesh::readUpdate()
                 false
             )
         );
+
+        // Derived info
+        bounds_ = boundBox(points_);
+
+        // Rotation can cause direction vector to change
+        geometricD_ = Vector<label>::zero;
+        solutionD_ = Vector<label>::zero;
         
         return polyMesh::POINTS_MOVED;
     }
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
index 0c6655862c..30c5839cc4 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
@@ -58,35 +58,6 @@ defineTypeNameAndDebug(autoSnapDriver, 0);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::autoSnapDriver::getZonedSurfaces
-(
-    labelList& zonedSurfaces,
-    labelList& unzonedSurfaces
-) const
-{    // Surfaces with zone information
-    const wordList& faceZoneNames = meshRefiner_.surfaces().faceZoneNames();
-
-    zonedSurfaces.setSize(faceZoneNames.size());
-    label zonedI = 0;
-    unzonedSurfaces.setSize(faceZoneNames.size());
-    label unzonedI = 0;
-
-    forAll(faceZoneNames, surfI)
-    {
-        if (faceZoneNames[surfI].size())
-        {
-            zonedSurfaces[zonedI++] = surfI;
-        }
-        else
-        {
-            unzonedSurfaces[unzonedI++] = surfI;
-        }
-    }
-    zonedSurfaces.setSize(zonedI);
-    unzonedSurfaces.setSize(unzonedI);
-}
-
-
 // Get faces to repatch. Returns map from face to patch.
 Foam::Map<Foam::label> Foam::autoSnapDriver::getZoneBafflePatches
 (
@@ -711,9 +682,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::createZoneBaffles
     List<labelPair>& baffles
 )
 {
-    labelList zonedSurfaces;
-    labelList unzonedSurfaces;
-    getZonedSurfaces(zonedSurfaces, unzonedSurfaces);
+    labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
 
     autoPtr<mapPolyMesh> map;
 
@@ -798,9 +767,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::mergeZoneBaffles
     const List<labelPair>& baffles
 )
 {
-    labelList zonedSurfaces;
-    labelList unzonedSurfaces;
-    getZonedSurfaces(zonedSurfaces, unzonedSurfaces);
+    labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
 
     autoPtr<mapPolyMesh> map;
 
@@ -1048,9 +1015,10 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
         labelList snapSurf(localPoints.size(), -1);
 
         // Divide surfaces into zoned and unzoned
-        labelList zonedSurfaces;
-        labelList unzonedSurfaces;
-        getZonedSurfaces(zonedSurfaces, unzonedSurfaces);
+        labelList zonedSurfaces =
+            meshRefiner_.surfaces().getNamedSurfaces();
+        labelList unzonedSurfaces =
+            meshRefiner_.surfaces().getUnnamedSurfaces();
 
 
         // 1. All points to non-interface surfaces
@@ -1368,9 +1336,8 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::repatchToSurface
     indirectPrimitivePatch& pp = ppPtr();
 
     // Divide surfaces into zoned and unzoned
-    labelList zonedSurfaces;
-    labelList unzonedSurfaces;
-    getZonedSurfaces(zonedSurfaces, unzonedSurfaces);
+    labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
+    labelList unzonedSurfaces = meshRefiner_.surfaces().getUnnamedSurfaces();
 
 
     // Faces that do not move
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
index b2d0f28911..5d9b2bc357 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
@@ -84,9 +84,6 @@ class autoSnapDriver
 
         // Snapping
 
-             //- Split surfaces into non-zoned and zones ones
-            void getZonedSurfaces(labelList&, labelList&) const;
-
             //- Get faces to repatch. Returns map from face to patch.
             Map<label> getZoneBafflePatches(const bool allowBoundary) const;
 
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 7720386219..0025f33332 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -374,11 +374,28 @@ private:
                 const labelList&
             ) const;
 
+            bool isCollapsedFace
+            (
+                const pointField&,
+                const pointField& neiCc,
+                const scalar minFaceArea,
+                const scalar maxNonOrtho,
+                const label faceI
+            ) const;
+
+            bool isCollapsedCell
+            (
+                const pointField&,
+                const scalar volFraction,
+                const label cellI
+            ) const;
+
             //- Returns list with for every internal face -1 or the patch
             //  they should be baffled into. If removeEdgeConnectedCells is set
             //  removes cells based on perpendicularAngle.
             labelList markFacesOnProblemCells
             (
+                const dictionary& motionDict,
                 const bool removeEdgeConnectedCells,
                 const scalarField& perpendicularAngle,
                 const labelList& globalToPatch
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index 798d6d5ae4..c5941aeadc 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -1503,6 +1503,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
         (
             markFacesOnProblemCells
             (
+                motionDict,
                 removeEdgeConnectedCells,
                 perpendicularAngle,
                 globalToPatch
@@ -1548,6 +1549,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
             (
                 markFacesOnProblemCells
                 (
+                    motionDict,
                     removeEdgeConnectedCells,
                     perpendicularAngle,
                     globalToPatch
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
index 37d38cac48..ea1e8a083c 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
@@ -278,13 +278,105 @@ Foam::Map<Foam::label> Foam::meshRefinement::findEdgeConnectedProblemCells
 }
 
 
+// Check if moving face to new points causes a 'collapsed' face.
+// Uses new point position only for the face, not the neighbouring
+// cell centres
+bool Foam::meshRefinement::isCollapsedFace
+(
+    const pointField& points,
+    const pointField& neiCc,
+    const scalar minFaceArea,
+    const scalar maxNonOrtho,
+    const label faceI
+) const
+{
+    vector s = mesh_.faces()[faceI].normal(points);
+    scalar magS = mag(s);
+
+    // Check face area
+    if (magS < minFaceArea)
+    {
+        return true;
+    }
+
+    // Check orthogonality
+    const point& ownCc = mesh_.cellCentres()[mesh_.faceOwner()[faceI]];
+
+    if (mesh_.isInternalFace(faceI))
+    {
+        label nei = mesh_.faceNeighbour()[faceI];
+        vector d = ownCc - mesh_.cellCentres()[nei];
+
+        scalar dDotS = (d & s)/(mag(d)*magS + VSMALL);
+        
+        if (dDotS < maxNonOrtho)
+        {
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    }
+    else
+    {
+        label patchI = mesh_.boundaryMesh().whichPatch(faceI);
+
+        if (mesh_.boundaryMesh()[patchI].coupled())
+        {
+            vector d = ownCc - neiCc[faceI-mesh_.nInternalFaces()];
+
+            scalar dDotS = (d & s)/(mag(d)*magS + VSMALL);
+
+            if (dDotS < maxNonOrtho)
+            {
+                return true;
+            }
+            else
+            {
+                return false;
+            }
+        }
+        else
+        {
+            // Collapsing normal boundary face does not cause problems with
+            // non-orthogonality
+            return true;
+        }
+    }
+}
+
+
+// Check if moving cell to new points causes it to collapse.
+bool Foam::meshRefinement::isCollapsedCell
+(
+    const pointField& points,
+    const scalar volFraction,
+    const label cellI
+) const
+{
+    scalar vol = mesh_.cells()[cellI].mag(points, mesh_.faces());
+
+    if (vol/mesh_.cellVolumes()[cellI] < volFraction)
+    {
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
 // Returns list with for every internal face -1 or the patch they should
-// be baffled into. Gets run after createBaffles so all the surface
-// intersections have already been turned into baffles. Used to remove cells
-// by baffling all their faces and have the splitMeshRegions chuck away non
-// used regions.
+// be baffled into. Gets run after createBaffles so all the unzoned surface
+// intersections have already been turned into baffles. (Note: zoned surfaces
+// are not baffled at this stage)
+// Used to remove cells by baffling all their faces and have the
+// splitMeshRegions chuck away non used regions.
 Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
 (
+    const dictionary& motionDict,
     const bool removeEdgeConnectedCells,
     const scalarField& perpendicularAngle,
     const labelList& globalToPatch
@@ -294,6 +386,10 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
     const labelList& pointLevel = meshCutter_.pointLevel();
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
+    // Swap neighbouring cell centres and cell level
+    labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
+    pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
+    calcNeighbourData(neiLevel, neiCc);
 
     // Per internal face (boundary faces not used) the patch that the
     // baffle should get (or -1)
@@ -334,6 +430,9 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
     // Count of faces marked for baffling
     label nBaffleFaces = 0;
 
+    // Count of faces not baffled since would not cause a collapse
+    label nPrevented = 0;
+
     if (removeEdgeConnectedCells && max(perpendicularAngle) >= 0)
     {
         Info<< "markFacesOnProblemCells :"
@@ -420,6 +519,74 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
     );
 
 
+    // See if checking for collapse
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Collapse checking parameters
+    scalar volFraction = -1;
+    if (motionDict.found("minVolCollapseRatio"))
+    {
+        volFraction = readScalar(motionDict.lookup("minVolCollapseRatio"));
+    }
+    const bool checkCollapse = (volFraction > 0);
+    scalar minArea = -1;
+    scalar maxNonOrtho = -1;
+
+
+    // Find nearest (non-baffle) surface
+    pointField newPoints;
+
+    if (checkCollapse)
+    {
+        minArea = readScalar(motionDict.lookup("minArea"));
+        maxNonOrtho = readScalar(motionDict.lookup("maxNonOrtho"));
+
+        Info<< "markFacesOnProblemCells :"
+            << " Deleting all-anchor surface cells only if"
+            << "snapping them violates mesh quality constraints:" << nl
+            << "    snapped/original cell volume < " << volFraction << nl
+            << "    face area                    < " << minArea << nl
+            << "    non-orthogonality            > " << maxNonOrtho << nl
+            << endl;
+
+        // Construct addressing engine.
+        autoPtr<indirectPrimitivePatch> ppPtr
+        (
+            meshRefinement::makePatch
+            (
+                mesh_,
+                adaptPatchIDs
+            )
+        );
+        const indirectPrimitivePatch& pp = ppPtr();
+        const pointField& localPoints = pp.localPoints();
+        const labelList& meshPoints = pp.meshPoints();
+
+        List<pointIndexHit> hitInfo;
+        labelList hitSurface;
+        surfaces_.findNearest
+        (
+            surfaces_.getUnnamedSurfaces(),
+            localPoints,
+            scalarField(localPoints.size(), sqr(GREAT)),    // sqr of attraction
+            hitSurface,
+            hitInfo
+        );
+
+        // Start of from current points
+        newPoints = mesh_.points();
+
+        forAll(hitInfo, i)
+        {
+            if (hitInfo[i].hit())
+            {
+                newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
+            }
+        }
+    }
+
+
+
     // For each cell count the number of anchor points that are on
     // the boundary:
     // 8 : check the number of (baffle) boundary faces. If 3 or more block
@@ -430,6 +597,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
     //     3 or more of these cells. (note that on a flat surface a boundary
     //     point will only have 4 cells connected to it)
 
+
     // Does cell have exactly 7 of its 8 anchor points on the boundary?
     PackedBoolList hasSevenBoundaryAnchorPoints(mesh_.nCells());
     // If so what is the remaining non-boundary anchor point?
@@ -494,23 +662,45 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
             // Eugene: delete cell no matter what.
             //if (nBfaces > 1)
             {
-                forAll(cFaces, cf)
+                if
+                (
+                    checkCollapse
+                && !isCollapsedCell(newPoints, volFraction, cellI)
+                )
                 {
-                    label faceI = cFaces[cf];
-
-                    if (facePatch[faceI] == -1 && mesh_.isInternalFace(faceI))
+                    nPrevented++;
+                    //Pout<< "Preventing collapse of 8 anchor point cell "
+                    //    << cellI << " at " << mesh_.cellCentres()[cellI]
+                    //    << " since new volume "
+                    //    << mesh_.cells()[cellI].mag(newPoints, mesh_.faces())
+                    //    << " old volume " << mesh_.cellVolumes()[cellI]
+                    //    << endl;
+                }
+                else
+                {
+                    // Block all faces of cell
+                    forAll(cFaces, cf)
                     {
-                        facePatch[faceI] = getBafflePatch(facePatch, faceI);
-                        nBaffleFaces++;
+                        label faceI = cFaces[cf];
 
-                        // Mark face as a 'boundary'
-                        markBoundaryFace
+                        if
                         (
-                            faceI,
-                            isBoundaryFace,
-                            isBoundaryEdge,
-                            isBoundaryPoint
-                        );
+                            facePatch[faceI] == -1
+                         && mesh_.isInternalFace(faceI)
+                        )
+                        {
+                            facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                            nBaffleFaces++;
+
+                            // Mark face as a 'boundary'
+                            markBoundaryFace
+                            (
+                                faceI,
+                                isBoundaryFace,
+                                isBoundaryEdge,
+                                isBoundaryPoint
+                            );
+                        }
                     }
                 }
             }
@@ -556,29 +746,52 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
 
                 if (hasSevenBoundaryAnchorPoints.get(cellI) == 1u)
                 {
-                    const cell& cFaces = mesh_.cells()[cellI];
-
-                    forAll(cFaces, cf)
+                    if
+                    (
+                        checkCollapse
+                    && !isCollapsedCell(newPoints, volFraction, cellI)
+                    )
                     {
-                        label faceI = cFaces[cf];
+                        nPrevented++;
+                        //Pout<< "Preventing collapse of 7 anchor cell "
+                        //    << cellI
+                        //    << " at " << mesh_.cellCentres()[cellI]
+                        //    << " since new volume "
+                        //    << mesh_.cells()[cellI].mag
+                        //        (newPoints, mesh_.faces())
+                        //    << " old volume " << mesh_.cellVolumes()[cellI]
+                        //    << endl;
+                    }
+                    else
+                    {
+                        const cell& cFaces = mesh_.cells()[cellI];
 
-                        if
-                        (
-                            facePatch[faceI] == -1
-                         && mesh_.isInternalFace(faceI)
-                        )
+                        forAll(cFaces, cf)
                         {
-                            facePatch[faceI] = getBafflePatch(facePatch, faceI);
-                            nBaffleFaces++;
+                            label faceI = cFaces[cf];
 
-                            // Mark face as a 'boundary'
-                            markBoundaryFace
+                            if
                             (
-                                faceI,
-                                isBoundaryFace,
-                                isBoundaryEdge,
-                                isBoundaryPoint
-                            );
+                                facePatch[faceI] == -1
+                             && mesh_.isInternalFace(faceI)
+                            )
+                            {
+                                facePatch[faceI] = getBafflePatch
+                                (
+                                    facePatch,
+                                    faceI
+                                );
+                                nBaffleFaces++;
+
+                                // Mark face as a 'boundary'
+                                markBoundaryFace
+                                (
+                                    faceI,
+                                    isBoundaryFace,
+                                    isBoundaryEdge,
+                                    isBoundaryPoint
+                                );
+                            }
                         }
                     }
                 }
@@ -635,11 +848,33 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
 
             if (nFaceBoundaryEdges == fEdges.size())
             {
-                facePatch[faceI] = getBafflePatch(facePatch, faceI);
-                nBaffleFaces++;
+                if
+                (
+                    checkCollapse
+                && !isCollapsedFace
+                    (
+                        newPoints,
+                        neiCc,
+                        minArea,
+                        maxNonOrtho,
+                        faceI
+                    )
+                )
+                {
+                    nPrevented++;
+                    //Pout<< "Preventing collapse of face " << faceI
+                    //    << " with all boundary edges "
+                    //    << " at " << mesh_.faceCentres()[faceI]
+                    //    << endl;
+                }
+                else
+                {
+                    facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                    nBaffleFaces++;
 
-                // Do NOT update boundary data since this would grow blocked
-                // faces across gaps.
+                    // Do NOT update boundary data since this would grow blocked
+                    // faces across gaps.
+                }
             }
         }
     }
@@ -669,11 +904,34 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
 
                     if (nFaceBoundaryEdges == fEdges.size())
                     {
-                        facePatch[faceI] = getBafflePatch(facePatch, faceI);
-                        nBaffleFaces++;
+                        if
+                        (
+                            checkCollapse
+                        && !isCollapsedFace
+                            (
+                                newPoints,
+                                neiCc,
+                                minArea,
+                                maxNonOrtho,
+                                faceI
+                            )
+                        )
+                        {
+                            nPrevented++;
+                            //Pout<< "Preventing collapse of coupled face "
+                            //    << faceI
+                            //    << " with all boundary edges "
+                            //    << " at " << mesh_.faceCentres()[faceI]
+                            //    << endl;
+                        }
+                        else
+                        {
+                            facePatch[faceI] = getBafflePatch(facePatch, faceI);
+                            nBaffleFaces++;
 
-                        // Do NOT update boundary data since this would grow
-                        // blocked faces across gaps.
+                            // Do NOT update boundary data since this would grow
+                            // blocked faces across gaps.
+                        }
                     }
                 }
 
@@ -687,266 +945,181 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
         << " additional internal faces to be converted into baffles."
         << endl;
 
-    return facePatch;
-}
-//XXXXXXXXXXXXXX
-// Mark faces to be baffled to prevent snapping problems. Does
-// test to find nearest surface and checks which faces would get squashed.
-Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
-(
-    const dictionary& motionDict,
-    const labelList& globalToPatch
-) const
-{
-    // Get the labels of added patches.
-    labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
-
-    // Construct addressing engine.
-    autoPtr<indirectPrimitivePatch> ppPtr
-    (
-        meshRefinement::makePatch
-        (
-            mesh_,
-            adaptPatchIDs
-        )
-    );
-    const indirectPrimitivePatch& pp = ppPtr();
-    const pointField& localPoints = pp.localPoints();
-    const labelList& meshPoints = pp.meshPoints();
-
-    // Find nearest (non-baffle) surface
-    pointField newPoints(mesh_.points());
+    if (checkCollapse)
     {
-        List<pointIndexHit> hitInfo;
-        labelList hitSurface;
-        surfaces_.findNearest
-        (
-            surfaces_.getUnnamedSurfaces(),
-            localPoints,
-            scalarField(localPoints.size(), sqr(GREAT)),    // sqr of attraction
-            hitSurface,
-            hitInfo
-        );
-
-        forAll(hitInfo, i)
-        {
-            if (hitInfo[i].hit())
-            {
-                //label pointI = meshPoints[i];
-                //Pout<< "   " << pointI << " moved from "
-                //    << mesh_.points()[pointI] << " by "
-                //    << mag(hitInfo[i].hitPoint()-mesh_.points()[pointI])
-                //    << endl;
-                newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
-            }
-        }
+        Info<< "markFacesOnProblemCells : prevented "
+            << returnReduce(nPrevented, sumOp<label>())
+            << " internal faces fom getting converted into baffles."
+            << endl;
     }
 
-    // Per face (internal or coupled!) the patch that the
-    // baffle should get (or -1).
-    labelList facePatch(mesh_.nFaces(), -1);
-    // Count of baffled faces
-    label nBaffleFaces = 0;
+    return facePatch;
+}
 
 
-//    // Sync position? Or not since same face on both side so just sync
-//    // result of baffle.
+//// Mark faces to be baffled to prevent snapping problems. Does
+//// test to find nearest surface and checks which faces would get squashed.
+//Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
+//(
+//    const dictionary& motionDict,
+//    const labelList& globalToPatch
+//) const
+//{
+//    // Get the labels of added patches.
+//    labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
 //
-//    const scalar minArea(readScalar(motionDict.lookup("minArea")));
+//    // Construct addressing engine.
+//    autoPtr<indirectPrimitivePatch> ppPtr
+//    (
+//        meshRefinement::makePatch
+//        (
+//            mesh_,
+//            adaptPatchIDs
+//        )
+//    );
+//    const indirectPrimitivePatch& pp = ppPtr();
+//    const pointField& localPoints = pp.localPoints();
+//    const labelList& meshPoints = pp.meshPoints();
 //
-//    Pout<< "markFacesOnProblemCellsGeometric : Comparing to minArea:"
-//        << minArea << endl;
-//
-//    pointField facePoints;
-//    for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
+//    // Find nearest (non-baffle) surface
+//    pointField newPoints(mesh_.points());
 //    {
-//        const face& f = mesh_.faces()[faceI];
+//        List<pointIndexHit> hitInfo;
+//        labelList hitSurface;
+//        surfaces_.findNearest
+//        (
+//            surfaces_.getUnnamedSurfaces(),
+//            localPoints,
+//            scalarField(localPoints.size(), sqr(GREAT)),// sqr of attraction
+//            hitSurface,
+//            hitInfo
+//        );
 //
-//        bool usesPatchPoint = false;
-//
-//        facePoints.setSize(f.size());
-//        forAll(f, fp)
+//        forAll(hitInfo, i)
 //        {
-//            Map<label>::const_iterator iter = pp.meshPointMap().find(f[fp]);
-//
-//            if (iter != pp.meshPointMap().end())
+//            if (hitInfo[i].hit())
 //            {
-//                facePoints[fp] = newPosition[iter()];
-//                usesPatchPoint = true;
+//                //label pointI = meshPoints[i];
+//                //Pout<< "   " << pointI << " moved from "
+//                //    << mesh_.points()[pointI] << " by "
+//                //    << mag(hitInfo[i].hitPoint()-mesh_.points()[pointI])
+//                //    << endl;
+//                newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
 //            }
-//            else
+//        }
+//    }
+//
+//    // Per face (internal or coupled!) the patch that the
+//    // baffle should get (or -1).
+//    labelList facePatch(mesh_.nFaces(), -1);
+//    // Count of baffled faces
+//    label nBaffleFaces = 0;
+//
+//    {
+//        pointField oldPoints(mesh_.points());
+//        mesh_.movePoints(newPoints);
+//        faceSet wrongFaces(mesh_, "wrongFaces", 100);
+//        {
+//            //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
+//
+//            // Just check the errors from squashing
+//            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+//
+//            const labelList allFaces(identity(mesh_.nFaces()));
+//            label nWrongFaces = 0;
+//
+//            scalar minArea(readScalar(motionDict.lookup("minArea")));
+//            if (minArea > -SMALL)
 //            {
-//                facePoints[fp] = mesh_.points()[f[fp]];
+//                polyMeshGeometry::checkFaceArea
+//                (
+//                    false,
+//                    minArea,
+//                    mesh_,
+//                    mesh_.faceAreas(),
+//                    allFaces,
+//                    &wrongFaces
+//                );
+//
+//                label nNewWrongFaces = returnReduce
+//                (
+//                    wrongFaces.size(),
+//                    sumOp<label>()
+//                );
+//
+//                Info<< "    faces with area < "
+//                    << setw(5) << minArea
+//                    << " m^2                            : "
+//                    << nNewWrongFaces-nWrongFaces << endl;
+//
+//                nWrongFaces = nNewWrongFaces;
+//            }
+//
+////            scalar minDet(readScalar(motionDict.lookup("minDeterminant")));
+//            scalar minDet = 0.01;
+//            if (minDet > -1)
+//            {
+//                polyMeshGeometry::checkCellDeterminant
+//                (
+//                    false,
+//                    minDet,
+//                    mesh_,
+//                    mesh_.faceAreas(),
+//                    allFaces,
+//                    polyMeshGeometry::affectedCells(mesh_, allFaces),
+//                    &wrongFaces
+//                );
+//
+//                label nNewWrongFaces = returnReduce
+//                (
+//                    wrongFaces.size(),
+//                    sumOp<label>()
+//                );
+//
+//                Info<< "    faces on cells with determinant < "
+//                    << setw(5) << minDet << "                : "
+//                    << nNewWrongFaces-nWrongFaces << endl;
+//
+//                nWrongFaces = nNewWrongFaces;
 //            }
 //        }
 //
-//        if (usesPatchPoint)
-//        {
-//            // Check area of face wrt original area
-//            face identFace(identity(f.size()));
 //
-//            if (identFace.mag(facePoints) < minArea)
+//        forAllConstIter(faceSet, wrongFaces, iter)
+//        {
+//            label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
+//
+//            if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
 //            {
-//                facePatch[faceI] = getBafflePatch(facePatch, faceI);
+//                facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
 //                nBaffleFaces++;
+//
+//                //Pout<< "    " << iter.key()
+//                //    //<< " on patch " << mesh_.boundaryMesh()[patchI].name()
+//                //    << " is destined for patch " << facePatch[iter.key()]
+//                //    << endl;
 //            }
 //        }
+//        // Restore points.
+//        mesh_.movePoints(oldPoints);
 //    }
 //
 //
-//    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-//    forAll(patches, patchI)
-//    {
-//        const polyPatch& pp = patches[patchI];
+//    Info<< "markFacesOnProblemCellsGeometric : marked "
+//        << returnReduce(nBaffleFaces, sumOp<label>())
+//        << " additional internal and coupled faces"
+//        << " to be converted into baffles." << endl;
 //
-//        if (pp.coupled())
-//        {
-//            forAll(pp, i)
-//            {
-//                label faceI = pp.start()+i;
+//    syncTools::syncFaceList
+//    (
+//        mesh_,
+//        facePatch,
+//        maxEqOp<label>(),
+//        false               // no separation
+//    );
 //
-//                const face& f = mesh_.faces()[faceI];
-//
-//                bool usesPatchPoint = false;
-//
-//                facePoints.setSize(f.size());
-//                forAll(f, fp)
-//                {
-//                    Map<label>::const_iterator iter =
-//                        pp.meshPointMap().find(f[fp]);
-//
-//                    if (iter != pp.meshPointMap().end())
-//                    {
-//                        facePoints[fp] = newPosition[iter()];
-//                        usesPatchPoint = true;
-//                    }
-//                    else
-//                    {
-//                        facePoints[fp] = mesh_.points()[f[fp]];
-//                    }
-//                }
-//
-//                if (usesPatchPoint)
-//                {
-//                    // Check area of face wrt original area
-//                    face identFace(identity(f.size()));
-//
-//                    if (identFace.mag(facePoints) < minArea)
-//                    {
-//                        facePatch[faceI] = getBafflePatch(facePatch, faceI);
-//                        nBaffleFaces++;
-//                    }
-//                }
-//            }
-//        }
-//    }
-
-    {
-        pointField oldPoints(mesh_.points());
-        mesh_.movePoints(newPoints);
-        faceSet wrongFaces(mesh_, "wrongFaces", 100);
-        {
-            //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
-
-            // Just check the errors from squashing
-            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-            const labelList allFaces(identity(mesh_.nFaces()));
-            label nWrongFaces = 0;
-
-            scalar minArea(readScalar(motionDict.lookup("minArea")));
-            if (minArea > -SMALL)
-            {
-                polyMeshGeometry::checkFaceArea
-                (
-                    false,
-                    minArea,
-                    mesh_,
-                    mesh_.faceAreas(),
-                    allFaces,
-                    &wrongFaces
-                );
-
-                label nNewWrongFaces = returnReduce
-                (
-                    wrongFaces.size(),
-                    sumOp<label>()
-                );
-
-                Info<< "    faces with area < "
-                    << setw(5) << minArea
-                    << " m^2                            : "
-                    << nNewWrongFaces-nWrongFaces << endl;
-
-                nWrongFaces = nNewWrongFaces;
-            }
-
-//            scalar minDet(readScalar(motionDict.lookup("minDeterminant")));
-            scalar minDet = 0.01;
-            if (minDet > -1)
-            {
-                polyMeshGeometry::checkCellDeterminant
-                (
-                    false,
-                    minDet,
-                    mesh_,
-                    mesh_.faceAreas(),
-                    allFaces,
-                    polyMeshGeometry::affectedCells(mesh_, allFaces),
-                    &wrongFaces
-                );
-
-                label nNewWrongFaces = returnReduce
-                (
-                    wrongFaces.size(),
-                    sumOp<label>()
-                );
-
-                Info<< "    faces on cells with determinant < "
-                    << setw(5) << minDet << "                : "
-                    << nNewWrongFaces-nWrongFaces << endl;
-
-                nWrongFaces = nNewWrongFaces;
-            }
-        }
-
-
-        forAllConstIter(faceSet, wrongFaces, iter)
-        {
-            label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
-
-            if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
-            {
-                facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
-                nBaffleFaces++;
-
-                //Pout<< "    " << iter.key()
-                //    //<< " on patch " << mesh_.boundaryMesh()[patchI].name()
-                //    << " is destined for patch " << facePatch[iter.key()]
-                //    << endl;
-            }
-        }
-        // Restore points.
-        mesh_.movePoints(oldPoints);
-    }
-
-
-    Info<< "markFacesOnProblemCellsGeometric : marked "
-        << returnReduce(nBaffleFaces, sumOp<label>())
-        << " additional internal and coupled faces"
-        << " to be converted into baffles." << endl;
-
-    syncTools::syncFaceList
-    (
-        mesh_,
-        facePatch,
-        maxEqOp<label>(),
-        false               // no separation
-    );
-
-    return facePatch;
-}
-//XXXXXXXX
+//    return facePatch;
+//}
 
 
 // ************************************************************************* //
diff --git a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
index 2c6fbbe8f2..e8fd361401 100644
--- a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
+++ b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
@@ -154,7 +154,7 @@ public:
                 return cellZoneNames_;
             }
 
-            //- Get indices of named surfaces (surfaces with faceZoneName)
+            //- Get indices of unnamed surfaces (surfaces without faceZoneName)
             labelList getUnnamedSurfaces() const;
 
             //- Get indices of named surfaces (surfaces with faceZoneName)
diff --git a/src/decompositionAgglomeration/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C b/src/decompositionAgglomeration/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
index 882123c526..72d39d215c 100644
--- a/src/decompositionAgglomeration/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
+++ b/src/decompositionAgglomeration/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
@@ -165,6 +165,11 @@ void Foam::hierarchGeomDecomp::findBinary
     // (one beyond) index of highValue
     label high = values.size();
 
+    // Safeguards to avoid infinite loop.
+    label lowPrev = -1;
+    label midPrev = -1;
+    label highPrev = -1;
+
     //while (low <= high)
     while (true)
     {
@@ -197,6 +202,20 @@ void Foam::hierarchGeomDecomp::findBinary
         // Update mid, midValue
         midValue = 0.5*(lowValue+highValue);
         mid = findLower(values, midValue, low, high);
+
+        // Safeguard if same as previous.
+        bool hasNotChanged = (mid == midPrev) && (low == lowPrev) && (high == highPrev);
+        if (returnReduce(hasNotChanged, andOp<bool>()))
+        {
+            WarningIn("hierarchGeomDecomp::findBinary(..)")
+                << "unable to find desired deomposition split, making do!" 
+                << endl;
+            break;
+        }
+
+        midPrev = mid;
+        lowPrev = low;
+        highPrev = high;
     }
 }
 
diff --git a/src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.C
index 47670dea8c..8643dba429 100644
--- a/src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.C
@@ -56,13 +56,7 @@ slicedFvPatchField<Type>::slicedFvPatchField
 )
 :
     fvPatchField<Type>(p, iF)
-{
-    notImplemented
-    (
-        "slicedFvPatchField<Type>::"
-        "slicedFvPatchField(const fvPatch&, const Field<Type>&)"
-    );
-}
+{}
 
 
 template<class Type>
diff --git a/src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.H
index c712e5bb5f..3d73cfc08d 100644
--- a/src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.H
@@ -75,7 +75,7 @@ public:
             const Field<Type>&
         );
 
-        //- Construct from patch and internal field
+        //- Construct from patch and internal field. Assign value later.
         slicedFvPatchField
         (
             const fvPatch&,
diff --git a/src/finiteVolume/fields/fvsPatchFields/basic/sliced/slicedFvsPatchField.C b/src/finiteVolume/fields/fvsPatchFields/basic/sliced/slicedFvsPatchField.C
index 9dffe160d9..7e4174119b 100644
--- a/src/finiteVolume/fields/fvsPatchFields/basic/sliced/slicedFvsPatchField.C
+++ b/src/finiteVolume/fields/fvsPatchFields/basic/sliced/slicedFvsPatchField.C
@@ -56,13 +56,7 @@ slicedFvsPatchField<Type>::slicedFvsPatchField
 )
 :
     fvsPatchField<Type>(p, iF)
-{
-    notImplemented
-    (
-        "slicedFvsPatchField<Type>::"
-        "slicedFvsPatchField(const fvPatch&, const Field<Type>&)"
-    );
-}
+{}
 
 
 template<class Type>
diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C b/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C
index a2491f1ab4..1a66028f95 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C
@@ -62,13 +62,7 @@ Foam::distanceSurface::interpolateField
     );
 
     // Sample.
-    return surface().interpolate
-    (
-        cellDistancePtr_(),
-        pointDistance_,
-        volFld,
-        pointFld()
-    );
+    return surface().interpolate(volFld, pointFld());
 }
 
 
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.C b/src/sampling/sampledSurface/isoSurface/isoSurface.C
index 78081050ab..35c1bec77d 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.C
@@ -30,6 +30,7 @@ License
 #include "syncTools.H"
 #include "addToRunTimeSelectionTable.H"
 #include "slicedVolFields.H"
+#include "volFields.H"
 #include "surfaceFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -41,6 +42,243 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+bool Foam::isoSurface::noTransform(const tensor& tt) const
+{
+    return
+        (mag(tt.xx()-1) < mergeDistance_)
+     && (mag(tt.yy()-1) < mergeDistance_)
+     && (mag(tt.zz()-1) < mergeDistance_)
+     && (mag(tt.xy()) < mergeDistance_)
+     && (mag(tt.xz()) < mergeDistance_)
+     && (mag(tt.yx()) < mergeDistance_)
+     && (mag(tt.yz()) < mergeDistance_)
+     && (mag(tt.zx()) < mergeDistance_)
+     && (mag(tt.zy()) < mergeDistance_);
+}
+
+
+bool Foam::isoSurface::collocatedPatch(const polyPatch& pp)
+{
+    const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>(pp);
+
+    return cpp.parallel() && !cpp.separated();
+}
+
+
+// Calculates per face whether couple is collocated.
+Foam::PackedBoolList Foam::isoSurface::collocatedFaces
+(
+    const coupledPolyPatch& pp
+) const
+{
+    // Initialise to false
+    PackedBoolList collocated(pp.size());
+
+    const vectorField& separation = pp.separation();
+    const tensorField& forwardT = pp.forwardT();
+
+    if (forwardT.size() == 0)
+    {
+        // Parallel.
+        if (separation.size() == 0)
+        {
+            collocated = 1u;
+        }
+        else if (separation.size() == 1)
+        {
+            // Fully separate. Do not synchronise.
+        }
+        else
+        {
+            // Per face separation.
+            forAll(pp, faceI)
+            {
+                if (mag(separation[faceI]) < mergeDistance_)
+                {
+                    collocated[faceI] = 1u;
+                }
+            }
+        }
+    }
+    else if (forwardT.size() == 1)
+    {
+        // Fully transformed.
+    }
+    else
+    {
+        // Per face transformation.
+        forAll(pp, faceI)
+        {
+            if (noTransform(forwardT[faceI]))
+            {
+                collocated[faceI] = 1u;
+            }
+        }
+    }
+    return collocated;
+}
+
+
+// Insert the data for local point patchPointI into patch local values
+// and/or into the shared values field.
+void Foam::isoSurface::insertPointData
+(
+    const processorPolyPatch& pp,
+    const Map<label>& meshToShared,
+    const pointField& pointValues,
+    const label patchPointI,
+    pointField& patchValues,
+    pointField& sharedValues
+) const
+{
+    label meshPointI = pp.meshPoints()[patchPointI];
+
+    // Store in local field
+    label nbrPointI = pp.neighbPoints()[patchPointI];
+    if (nbrPointI >= 0 && nbrPointI < patchValues.size())
+    {
+        minEqOp<point>()(patchValues[nbrPointI], pointValues[meshPointI]);
+    }
+
+    // Store in shared field
+    Map<label>::const_iterator iter = meshToShared.find(meshPointI);
+    if (iter != meshToShared.end())
+    {
+        minEqOp<point>()(sharedValues[iter()], pointValues[meshPointI]);
+    }
+}
+
+
+void Foam::isoSurface::syncUnseparatedPoints
+(
+    pointField& pointValues,
+    const point& nullValue
+) const
+{
+    // Until syncPointList handles separated coupled patches with multiple
+    // transforms do our own synchronisation of non-separated patches only
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    const globalMeshData& pd = mesh_.globalData();
+    const labelList& sharedPtAddr = pd.sharedPointAddr();
+    const labelList& sharedPtLabels = pd.sharedPointLabels();
+
+    // Create map from meshPoint to globalShared index.
+    Map<label> meshToShared(2*sharedPtLabels.size());
+    forAll(sharedPtLabels, i)
+    {
+        meshToShared.insert(sharedPtLabels[i], sharedPtAddr[i]);
+    }
+
+    // Values on shared points.
+    pointField sharedInfo(pd.nGlobalPoints(), nullValue);
+
+
+    if (Pstream::parRun())
+    {
+        // Send
+        forAll(patches, patchI)
+        {
+            if
+            (
+                isA<processorPolyPatch>(patches[patchI])
+             && patches[patchI].nPoints() > 0
+            )
+            {
+                const processorPolyPatch& pp =
+                    refCast<const processorPolyPatch>(patches[patchI]);
+                const labelList& meshPts = pp.meshPoints();
+
+                pointField patchInfo(meshPts.size(), nullValue);
+
+                PackedBoolList isCollocated(collocatedFaces(pp));
+
+                forAll(isCollocated, faceI)
+                {
+                    if (isCollocated[faceI])
+                    {
+                        const face& f = pp.localFaces()[faceI];
+
+                        forAll(f, fp)
+                        {
+                            label pointI = f[fp];
+
+                            insertPointData
+                            (
+                                pp,
+                                meshToShared,
+                                pointValues,
+                                pointI,
+                                patchInfo,
+                                sharedInfo
+                            );
+                        }
+                    }
+                }
+
+                OPstream toNbr(Pstream::blocking, pp.neighbProcNo());
+                toNbr << patchInfo;
+            }
+        }
+
+        // Receive and combine.
+
+        forAll(patches, patchI)
+        {
+            if
+            (
+                isA<processorPolyPatch>(patches[patchI])
+             && patches[patchI].nPoints() > 0
+            )
+            {
+                const processorPolyPatch& pp =
+                    refCast<const processorPolyPatch>(patches[patchI]);
+
+                pointField nbrPatchInfo(pp.nPoints());
+                {
+                    // We do not know the number of points on the other side
+                    // so cannot use Pstream::read.
+                    IPstream fromNbr(Pstream::blocking, pp.neighbProcNo());
+                    fromNbr >> nbrPatchInfo;
+                }
+
+                // Null any value which is not on neighbouring processor
+                nbrPatchInfo.setSize(pp.nPoints(), nullValue);
+
+                const labelList& meshPts = pp.meshPoints();
+
+                forAll(meshPts, pointI)
+                {
+                    label meshPointI = meshPts[pointI];
+                    minEqOp<point>()
+                    (
+                        pointValues[meshPointI],
+                        nbrPatchInfo[pointI]
+                    );
+                }
+            }
+        }
+    }
+
+
+    // Don't do cyclics for now. Are almost always separated anyway.
+
+
+    // Shared points
+
+    // Combine on master.
+    Pstream::listCombineGather(sharedInfo, minEqOp<point>());
+    Pstream::listCombineScatter(sharedInfo);
+
+    // Now we will all have the same information. Merge it back with
+    // my local information. (Note assignment and not combine operator)
+    forAll(sharedPtLabels, i)
+    {
+        label meshPointI = sharedPtLabels[i];
+        pointValues[meshPointI] = sharedInfo[sharedPtAddr[i]];
+    }
+}
+
+
 Foam::scalar Foam::isoSurface::isoFraction
 (
     const scalar s0,
@@ -89,6 +327,7 @@ bool Foam::isoSurface::isEdgeOfFaceCut
 void Foam::isoSurface::getNeighbour
 (
     const labelList& boundaryRegion,
+    const volVectorField& meshC,
     const volScalarField& cVals,
     const label cellI,
     const label faceI,
@@ -98,13 +337,12 @@ void Foam::isoSurface::getNeighbour
 {
     const labelList& own = mesh_.faceOwner();
     const labelList& nei = mesh_.faceNeighbour();
-    const surfaceScalarField& weights = mesh_.weights();
 
     if (mesh_.isInternalFace(faceI))
     {
         label nbr = (own[faceI] == cellI ? nei[faceI] : own[faceI]);
         nbrValue = cVals[nbr];
-        nbrPoint = mesh_.cellCentres()[nbr];
+        nbrPoint = meshC[nbr];
     }
     else
     {
@@ -113,38 +351,8 @@ void Foam::isoSurface::getNeighbour
         const polyPatch& pp = mesh_.boundaryMesh()[patchI];
         label patchFaceI = faceI-pp.start();
 
-        if (isA<emptyPolyPatch>(pp))
-        {
-            // Assume zero gradient
-            nbrValue = cVals[own[faceI]];
-            nbrPoint = mesh_.faceCentres()[faceI];
-        }
-        else if (pp.coupled())
-        {
-            if (!refCast<const coupledPolyPatch>(pp).separated())
-            {
-                // Behave as internal face:
-                // other side value
-                nbrValue = cVals.boundaryField()[patchI][patchFaceI];
-                // other side cell centre
-                nbrPoint = mesh_.C().boundaryField()[patchI][patchFaceI];
-            }
-            else
-            {
-                // Do some interpolation for now
-                const scalarField& w = weights.boundaryField()[patchI];
-
-                nbrPoint = mesh_.faceCentres()[faceI];
-                nbrValue =
-                    (1-w[patchFaceI])*cVals[own[faceI]]
-                  + w[patchFaceI]*cVals.boundaryField()[patchI][patchFaceI];
-            }
-        }
-        else
-        {
-            nbrValue = cVals.boundaryField()[patchI][patchFaceI];
-            nbrPoint = mesh_.faceCentres()[faceI];
-        }
+        nbrValue = cVals.boundaryField()[patchI][patchFaceI];
+        nbrPoint = meshC.boundaryField()[patchI][patchFaceI];
     }
 }
 
@@ -153,6 +361,7 @@ void Foam::isoSurface::getNeighbour
 void Foam::isoSurface::calcCutTypes
 (
     const labelList& boundaryRegion,
+    const volVectorField& meshC,
     const volScalarField& cVals,
     const scalarField& pVals
 )
@@ -174,6 +383,7 @@ void Foam::isoSurface::calcCutTypes
         getNeighbour
         (
             boundaryRegion,
+            meshC,
             cVals,
             own[faceI],
             faceI,
@@ -215,6 +425,7 @@ void Foam::isoSurface::calcCutTypes
             getNeighbour
             (
                 boundaryRegion,
+                meshC,
                 cVals,
                 own[faceI],
                 faceI,
@@ -408,6 +619,7 @@ Foam::pointIndexHit Foam::isoSurface::collapseSurface
 void Foam::isoSurface::calcSnappedCc
 (
     const labelList& boundaryRegion,
+    const volVectorField& meshC,
     const volScalarField& cVals,
     const scalarField& pVals,
 
@@ -448,6 +660,7 @@ void Foam::isoSurface::calcSnappedCc
                 getNeighbour
                 (
                     boundaryRegion,
+                    meshC,
                     cVals,
                     cellI,
                     faceI,
@@ -574,6 +787,7 @@ void Foam::isoSurface::calcSnappedPoint
 (
     const PackedBoolList& isBoundaryPoint,
     const labelList& boundaryRegion,
+    const volVectorField& meshC,
     const volScalarField& cVals,
     const scalarField& pVals,
 
@@ -639,6 +853,7 @@ void Foam::isoSurface::calcSnappedPoint
             getNeighbour
             (
                 boundaryRegion,
+                meshC,
                 cVals,
                 own,
                 faceI,
@@ -747,14 +962,22 @@ void Foam::isoSurface::calcSnappedPoint
         }
     }
 
-    syncTools::syncPointList
-    (
-        mesh_,
-        collapsedPoint,
-        minEqOp<point>(),
-        greatPoint,
-        true                // are coordinates so separate
-    );
+
+//Pout<< "**hack" << endl;
+//pointField oldCollaped(collapsedPoint);
+    syncUnseparatedPoints(collapsedPoint, greatPoint);
+//forAll(collapsedPoint, pointI)
+//{
+//    if (collapsedPoint[pointI] != oldCollaped[pointI])
+//    {
+//        Pout<< "**Synced point " << pointI
+//            << " coord:" << mesh_.points()[pointI]
+//            << " from " << oldCollaped[pointI]
+//            << " to " << collapsedPoint[pointI]
+//            << endl;
+//    }
+//}
+
 
     snappedPoint.setSize(mesh_.nPoints());
     snappedPoint = -1;
@@ -1534,6 +1757,7 @@ Foam::isoSurface::isoSurface
 )
 :
     mesh_(cVals.mesh()),
+    pVals_(pVals),
     iso_(iso),
     mergeDistance_(mergeTol*mesh_.bounds().mag())
 {
@@ -1545,8 +1769,8 @@ Foam::isoSurface::isoSurface
             << min(cVals.internalField()) << " / "
             << max(cVals.internalField()) << nl
             << "    point min/max : "
-            << min(pVals) << " / "
-            << max(pVals) << nl
+            << min(pVals_) << " / "
+            << max(pVals_) << nl
             << "    isoValue      : " << iso << nl
             << "    regularise    : " << regularise << nl
             << "    mergeTol      : " << mergeTol << nl
@@ -1555,6 +1779,89 @@ Foam::isoSurface::isoSurface
 
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
+
+    // Rewrite input field
+    // ~~~~~~~~~~~~~~~~~~~
+    // Rewrite input volScalarField to have interpolated values
+    // on separated patches.
+
+    cValsPtr_.reset(adaptPatchFields(cVals).ptr());
+
+
+    // Construct cell centres field consistent with cVals
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Generate field to interpolate. This is identical to the mesh.C()
+    // except on separated coupled patches and on empty patches.
+
+    slicedVolVectorField meshC
+    (
+        IOobject
+        (
+            "C",
+            mesh_.pointsInstance(),
+            mesh_.meshSubDir,
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        mesh_,
+        dimLength,
+        mesh_.cellCentres(),
+        mesh_.faceCentres()
+    );
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        // Adapt separated coupled (proc and cyclic) patches
+        if (isA<coupledPolyPatch>(pp) && !collocatedPatch(pp))
+        {
+            fvPatchVectorField& pfld = const_cast<fvPatchVectorField&>
+            (
+                meshC.boundaryField()[patchI]
+            );
+
+            PackedBoolList isCollocated
+            (
+                collocatedFaces(refCast<const coupledPolyPatch>(pp))
+            );
+
+            forAll(isCollocated, i)
+            {
+                if (!isCollocated[i])
+                {
+                    pfld[i] = mesh_.faceCentres()[pp.start()+i];
+                }
+            }
+        }
+        else if (isA<emptyPolyPatch>(pp))
+        {
+            typedef slicedVolVectorField::GeometricBoundaryField bType;
+
+            bType& bfld = const_cast<bType&>(meshC.boundaryField());
+
+            // Clear old value. Cannot resize it since is a slice.
+            bfld.set(patchI, NULL);
+
+            // Set new value we can change
+            bfld.set
+            (
+                patchI,
+                new calculatedFvPatchField<vector>
+                (
+                    mesh_.boundary()[patchI],
+                    meshC
+                )
+            );
+
+            // Change to face centres
+            bfld[patchI] = pp.patchSlice(mesh_.faceCentres());
+        }
+    }
+
+
+
     // Pre-calculate patch-per-face to avoid whichPatch call.
     labelList boundaryRegion(mesh_.nFaces()-mesh_.nInternalFaces());
 
@@ -1572,8 +1879,9 @@ Foam::isoSurface::isoSurface
     }
 
 
+
     // Determine if any cut through face/cell
-    calcCutTypes(boundaryRegion, cVals, pVals);
+    calcCutTypes(boundaryRegion, meshC, cValsPtr_(), pVals_);
 
 
     DynamicList<point> snappedPoints(nCutCells_);
@@ -1585,8 +1893,9 @@ Foam::isoSurface::isoSurface
         calcSnappedCc
         (
             boundaryRegion,
-            cVals,
-            pVals,
+            meshC,
+            cValsPtr_(),
+            pVals_,
 
             snappedPoints,
             snappedCc
@@ -1609,48 +1918,67 @@ Foam::isoSurface::isoSurface
     label nCellSnaps = snappedPoints.size();
 
 
-    // Determine if point is on boundary. Points on boundaries are never
-    // snapped. Coupled boundaries are handled explicitly so not marked here.
-    PackedBoolList isBoundaryPoint(mesh_.nPoints());
-
-    forAll(patches, patchI)
-    {
-        // Mark all boundary points that are not physically coupled (so anything
-        // but collocated coupled patches)
-        const polyPatch& pp = patches[patchI];
-
-        if
-        (
-           !pp.coupled()
-        || refCast<const coupledPolyPatch>(pp).separated()
-        )
-        {
-            label faceI = pp.start();
-
-            forAll(pp, i)
-            {
-                const face& f = mesh_.faces()[faceI];
-
-                forAll(f, fp)
-                {
-                    isBoundaryPoint.set(f[fp], 1);
-                }
-                faceI++;
-            }
-        }
-    }
-
-
     // Per point -1 or a point inside snappedPoints.
     labelList snappedPoint;
     if (regularise)
     {
+        // Determine if point is on boundary.
+        PackedBoolList isBoundaryPoint(mesh_.nPoints());
+
+        forAll(patches, patchI)
+        {
+            // Mark all boundary points that are not physically coupled
+            // (so anything but collocated coupled patches)
+
+            if (patches[patchI].coupled())
+            {
+                if (!collocatedPatch(patches[patchI]))
+                {
+                    const coupledPolyPatch& cpp =
+                        refCast<const coupledPolyPatch>
+                        (
+                            patches[patchI]
+                        );
+
+                    PackedBoolList isCollocated(collocatedFaces(cpp));
+
+                    forAll(isCollocated, i)
+                    {
+                        if (!isCollocated[i])
+                        {
+                            const face& f = mesh_.faces()[cpp.start()+i];
+
+                            forAll(f, fp)
+                            {
+                                isBoundaryPoint.set(f[fp], 1);
+                            }
+                        }
+                    }
+                }
+            }
+            else
+            {
+                const polyPatch& pp = patches[patchI];
+
+                forAll(pp, i)
+                {
+                    const face& f = mesh_.faces()[pp.start()+i];
+
+                    forAll(f, fp)
+                    {
+                        isBoundaryPoint.set(f[fp], 1);
+                    }
+                }
+            }
+        }
+
         calcSnappedPoint
         (
             isBoundaryPoint,
             boundaryRegion,
-            cVals,
-            pVals,
+            meshC,
+            cValsPtr_(),
+            pVals_,
 
             snappedPoints,
             snappedPoint
@@ -1669,80 +1997,14 @@ Foam::isoSurface::isoSurface
     }
 
 
-    // Generate field to interpolate. This is identical to the mesh.C()
-    // except on separated coupled patches and on empty patches.
-    slicedVolVectorField meshC
-    (
-        IOobject
-        (
-            "C",
-            mesh_.pointsInstance(),
-            mesh_.meshSubDir,
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        mesh_,
-        dimLength,
-        mesh_.cellCentres(),
-        mesh_.faceCentres()
-    );
-    {
-        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-        forAll(patches, patchI)
-        {
-            const polyPatch& pp = patches[patchI];
-
-            if
-            (
-                pp.coupled()
-             && refCast<const coupledPolyPatch>(pp).separated()
-            )
-            {
-                fvPatchVectorField& pfld = const_cast<fvPatchVectorField&>
-                (
-                    meshC.boundaryField()[patchI]
-                );
-                pfld.operator==
-                (
-                    pp.patchSlice(mesh_.faceCentres())
-                );
-            }
-            else if (isA<emptyPolyPatch>(pp))
-            {
-                typedef slicedVolVectorField::GeometricBoundaryField bType;
-
-                bType& bfld = const_cast<bType&>(meshC.boundaryField());
-
-                // Clear old value. Cannot resize it since slice.
-                bfld.set(patchI, NULL);
-
-                // Set new value we can change
-                bfld.set
-                (
-                    patchI,
-                    new calculatedFvPatchField<vector>
-                    (
-                        mesh_.boundary()[patchI],
-                        meshC
-                    )
-                );
-
-                // Change to face centres
-                bfld[patchI] = pp.patchSlice(mesh_.faceCentres());
-            }
-        }
-    }
-
 
     DynamicList<point> triPoints(nCutCells_);
     DynamicList<label> triMeshCells(nCutCells_);
 
     generateTriPoints
     (
-        cVals,
-        pVals,
+        cValsPtr_(),
+        pVals_,
 
         meshC,
         mesh_.points(),
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.H b/src/sampling/sampledSurface/isoSurface/isoSurface.H
index b859524b2c..1c323a7f60 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.H
@@ -56,6 +56,7 @@ SourceFiles
 #include "pointIndexHit.H"
 #include "PackedBoolList.H"
 #include "volFields.H"
+#include "slicedVolFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -92,6 +93,11 @@ class isoSurface
         //- Reference to mesh
         const fvMesh& mesh_;
 
+        const scalarField& pVals_;
+
+        //- Input volScalarField with separated coupled patches rewritten
+        autoPtr<slicedVolScalarField> cValsPtr_;
+
         //- Isosurface value
         const scalar iso_;
 
@@ -117,6 +123,38 @@ class isoSurface
 
     // Private Member Functions
 
+        // Point synchronisation
+
+            //- Does tensor differ (to within mergeTolerance) from identity
+            bool noTransform(const tensor& tt) const;
+
+            //- Is patch a collocated (i.e. non-separated) coupled patch?
+            static bool collocatedPatch(const polyPatch&);
+
+            //- Per face whether is collocated
+            PackedBoolList collocatedFaces(const coupledPolyPatch&) const;
+
+            //- Take value at local point patchPointI and assign it to its
+            //  correct place in patchValues (for transferral) and sharedValues
+            //  (for reduction)
+            void insertPointData
+            (
+                const processorPolyPatch& pp,
+                const Map<label>& meshToShared,
+                const pointField& pointValues,
+                const label patchPointI,
+                pointField& patchValues,
+                pointField& sharedValues
+            ) const;
+
+            //- Synchonise points on all non-separated coupled patches
+            void syncUnseparatedPoints
+            (
+                pointField& collapsedPoint,
+                const point& nullValue
+            ) const;
+
+
         //- Get location of iso value as fraction inbetween s0,s1
         scalar isoFraction
         (
@@ -136,6 +174,7 @@ class isoSurface
         void getNeighbour
         (
             const labelList& boundaryRegion,
+            const volVectorField& meshC,
             const volScalarField& cVals,
             const label cellI,
             const label faceI,
@@ -147,6 +186,7 @@ class isoSurface
         void calcCutTypes
         (
             const labelList& boundaryRegion,
+            const volVectorField& meshC,
             const volScalarField& cVals,
             const scalarField& pVals
         );
@@ -170,6 +210,7 @@ class isoSurface
         void calcSnappedCc
         (
             const labelList& boundaryRegion,
+            const volVectorField& meshC,
             const volScalarField& cVals,
             const scalarField& pVals,
             DynamicList<point>& snappedPoints,
@@ -182,12 +223,23 @@ class isoSurface
         (
             const PackedBoolList& isBoundaryPoint,
             const labelList& boundaryRegion,
+            const volVectorField& meshC,
             const volScalarField& cVals,
             const scalarField& pVals,
             DynamicList<point>& snappedPoints,
             labelList& snappedPoint
         ) const;
 
+
+        //- Return input field with coupled (and empty) patch values rewritten
+        template<class Type>
+        tmp<SlicedGeometricField
+        <Type, fvPatchField, slicedFvPatchField, volMesh> >
+        adaptPatchFields
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& fld
+        ) const;
+
         //- Generate single point by interpolation or snapping
         template<class Type>
         Type generatePoint
@@ -345,8 +397,8 @@ public:
     // Constructors
 
         //- Construct from cell values and point values. Uses boundaryField
-        //  for boundary values. Requires on coupled patchfields to be set
-        //  to the opposite cell value.
+        //  for boundary values. Holds reference to cellIsoVals and
+        //  pointIsoVals.
         isoSurface
         (
             const volScalarField& cellIsoVals,
@@ -376,8 +428,6 @@ public:
         template <class Type>
         tmp<Field<Type> > interpolate
         (
-            const volScalarField& cellIsoVals,
-            const scalarField& pointIsoVals,
             const GeometricField<Type, fvPatchField, volMesh>& cCoords,
             const Field<Type>& pCoords
         ) const;
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
index 0760b06a52..d7c72facc5 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
@@ -27,9 +27,114 @@ License
 #include "isoSurface.H"
 #include "polyMesh.H"
 #include "syncTools.H"
+#include "surfaceFields.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+template<class Type>
+Foam::tmp<Foam::SlicedGeometricField
+<
+    Type,
+    Foam::fvPatchField,
+    Foam::slicedFvPatchField,
+    Foam::volMesh
+> >
+Foam::isoSurface::adaptPatchFields
+(
+    const GeometricField<Type, fvPatchField, volMesh>& fld
+) const
+{
+    typedef SlicedGeometricField
+    <
+        Type,
+        fvPatchField,
+        slicedFvPatchField,
+        volMesh
+    > FieldType;
+
+    tmp<FieldType> tsliceFld
+    (
+        new FieldType
+        (
+            IOobject
+            (
+                fld.name(),
+                fld.instance(),
+                fld.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            fld,        // internal field
+            true        // preserveCouples
+        )
+    );
+    FieldType& sliceFld = tsliceFld();
+
+    const fvMesh& mesh = fld.mesh();
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (isA<emptyPolyPatch>(pp))
+        {
+            // Clear old value. Cannot resize it since is a slice.
+            sliceFld.boundaryField().set(patchI, NULL);
+
+            // Set new value we can change
+            sliceFld.boundaryField().set
+            (
+                patchI,
+                new calculatedFvPatchField<Type>
+                (
+                    mesh.boundary()[patchI],
+                    sliceFld
+                )
+            );
+            sliceFld.boundaryField()[patchI] ==
+                mesh.boundary()[patchI].patchInternalField
+                (
+                    sliceFld
+                );
+        }
+        else if (isA<cyclicPolyPatch>(pp))
+        {
+            // Already has interpolate as value
+        }
+        else if (isA<processorPolyPatch>(pp) && !collocatedPatch(pp))
+        {
+            fvPatchField<Type>& pfld = const_cast<fvPatchField<Type>&>
+            (
+                fld.boundaryField()[patchI]
+            );
+
+            const scalarField& w = mesh.weights().boundaryField()[patchI];
+
+            tmp<Field<Type> > f =
+                w*pfld.patchInternalField()
+              + (1.0-w)*pfld.patchNeighbourField();
+
+            PackedBoolList isCollocated
+            (
+                collocatedFaces(refCast<const processorPolyPatch>(pp))
+            );
+
+            forAll(isCollocated, i)
+            {
+                if (!isCollocated[i])
+                {
+                    pfld[i] = f()[i];
+                }
+            }
+        }
+    }
+    return tsliceFld;
+}
+
+
 template<class Type>
 Type Foam::isoSurface::generatePoint
 (
@@ -389,7 +494,6 @@ void Foam::isoSurface::generateTriPoints
     }
 
 
-
     // Generate triangle points
 
     triPoints.clear();
@@ -456,16 +560,16 @@ void Foam::isoSurface::generateTriPoints
     syncTools::swapBoundaryFaceList(mesh_, neiSnappedPoint, false);
 
 
+
     forAll(patches, patchI)
     {
         const polyPatch& pp = patches[patchI];
 
-        if
-        (
-            isA<processorPolyPatch>(pp)
-         && !refCast<const processorPolyPatch>(pp).separated()
-        )
+        if (isA<processorPolyPatch>(pp) && collocatedPatch(pp))
         {
+            // Coincident processor patches. Boundary field of 
+            // cVals and cCoords is opposite cell.
+
             //if (refCast<const processorPolyPatch>(pp).owner())
             {
                 label faceI = pp.start();
@@ -474,18 +578,6 @@ void Foam::isoSurface::generateTriPoints
                 {
                     if (faceCutType_[faceI] != NOTCUT)
                     {
-                        label bFaceI = faceI-mesh_.nInternalFaces();
-                        if
-                        (
-                            neiSnapped[bFaceI]
-                         && (neiSnappedPoint[bFaceI]==pTraits<Type>::zero)
-                        )
-                        {
-                            FatalErrorIn("isoSurface::generateTriPoints(..)")
-                                << "problem:" << abort(FatalError);
-                        }
-
-
                         generateFaceTriPoints
                         (
                             cVals,
@@ -512,84 +604,6 @@ void Foam::isoSurface::generateTriPoints
                 }
             }
         }
-        else if (isA<emptyPolyPatch>(pp))
-        {
-            // Check if field is there (when generating geometry the
-            // empty patches have been rewritten to be the face centres),
-            // otherwise use zero-gradient.
-
-            label faceI = pp.start();
-
-
-            const fvPatchScalarField& fvp = cVals.boundaryField()[patchI];
-
-            // Owner value of cVals
-            scalarField internalVals;
-            if (fvp.size() == 0)
-            {
-                internalVals.setSize(pp.size());
-                forAll(pp, i)
-                {
-                    internalVals[i] = cVals[own[pp.start()+i]];
-                }
-            }
-            const scalarField& bVals =
-            (
-                fvp.size() > 0
-              ? static_cast<const scalarField&>(fvp)
-              : internalVals
-            );
-
-
-            const fvPatchField<Type>& pc = cCoords.boundaryField()[patchI];
-
-            // Owner value of cCoords
-            Field<Type> internalCoords;
-            if (pc.size() == 0)
-            {
-                internalCoords.setSize(pp.size());
-                forAll(pp, i)
-                {
-                    internalCoords[i] = cCoords[own[pp.start()+i]];
-                }
-            }
-            const Field<Type>& bCoords =
-            (
-                pc.size() > 0
-              ? static_cast<const Field<Type>&>(pc)
-              : internalCoords
-            );
-
-
-            forAll(pp, i)
-            {
-                if (faceCutType_[faceI] != NOTCUT)
-                {
-                    generateFaceTriPoints
-                    (
-                        cVals,
-                        pVals,
-
-                        cCoords,
-                        pCoords,
-
-                        snappedPoints,
-                        snappedCc,
-                        snappedPoint,
-                        faceI,
-
-                        bVals[i],
-                        bCoords[i],
-                        false,  // fc not snapped
-                        pTraits<Type>::zero,
-
-                        triPoints,
-                        triMeshCells
-                    );
-                }
-                faceI++;
-            }
-        }
         else
         {
             label faceI = pp.start();
@@ -642,12 +656,20 @@ template <class Type>
 Foam::tmp<Foam::Field<Type> >
 Foam::isoSurface::interpolate
 (
-    const volScalarField& cVals,
-    const scalarField& pVals,
     const GeometricField<Type, fvPatchField, volMesh>& cCoords,
     const Field<Type>& pCoords
 ) const
 {
+    // Recalculate boundary values
+    tmp<SlicedGeometricField
+    <
+        Type,
+        fvPatchField,
+        slicedFvPatchField,
+        volMesh
+    > > c2(adaptPatchFields(cCoords));
+
+
     DynamicList<Type> triPoints(nCutCells_);
     DynamicList<label> triMeshCells(nCutCells_);
 
@@ -658,10 +680,10 @@ Foam::isoSurface::interpolate
 
     generateTriPoints
     (
-        cVals,
-        pVals,
+        cValsPtr_(),
+        pVals_,
 
-        cCoords,
+        c2(),
         pCoords,
 
         snappedPoints,
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
index 57267ac435..a8dfa78506 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
@@ -71,13 +71,7 @@ Foam::sampledIsoSurface::interpolateField
             volPointInterpolation::New(volSubFld.mesh()).interpolate(volSubFld);
 
         // Sample.
-        return surface().interpolate
-        (
-            *volSubFieldPtr_,
-            *pointSubFieldPtr_,
-            volSubFld,
-            tpointSubFld()
-        );
+        return surface().interpolate(volSubFld, tpointSubFld());
     }
     else
     {
@@ -85,13 +79,7 @@ Foam::sampledIsoSurface::interpolateField
             volPointInterpolation::New(volFld.mesh()).interpolate(volFld);
 
         // Sample.
-        return surface().interpolate
-        (
-            *volFieldPtr_,
-            *pointFieldPtr_,
-            volFld,
-            tpointFld()
-        );
+        return surface().interpolate(volFld, tpointFld());
     }
 }
 
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C
index f83aeca4e5..ad1d97ffe9 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C
@@ -67,13 +67,7 @@ Foam::sampledCuttingPlane::interpolateField
             volPointInterpolation::New(volSubFld.mesh()).interpolate(volSubFld);
 
         // Sample.
-        return surface().interpolate
-        (
-            cellDistancePtr_(),
-            pointDistance_,
-            volSubFld,
-            tpointSubFld()
-        );
+        return surface().interpolate(volSubFld, tpointSubFld());
     }
     else
     {
@@ -83,13 +77,7 @@ Foam::sampledCuttingPlane::interpolateField
         );
 
         // Sample.
-        return surface().interpolate
-        (
-            cellDistancePtr_(),
-            pointDistance_,
-            volFld,
-            tpointFld()
-        );
+        return surface().interpolate(volFld, tpointFld());
     }
 }