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());
     }
 }