diff --git a/applications/solvers/dsmc/dsmcFoam/Make/files b/applications/solvers/discreteMethods/dsmc/dsmcFoam/Make/files
similarity index 100%
rename from applications/solvers/dsmc/dsmcFoam/Make/files
rename to applications/solvers/discreteMethods/dsmc/dsmcFoam/Make/files
diff --git a/applications/solvers/dsmc/dsmcFoam/Make/options b/applications/solvers/discreteMethods/dsmc/dsmcFoam/Make/options
similarity index 100%
rename from applications/solvers/dsmc/dsmcFoam/Make/options
rename to applications/solvers/discreteMethods/dsmc/dsmcFoam/Make/options
diff --git a/applications/solvers/dsmc/dsmcFoam/createFields.H b/applications/solvers/discreteMethods/dsmc/dsmcFoam/createFields.H
similarity index 100%
rename from applications/solvers/dsmc/dsmcFoam/createFields.H
rename to applications/solvers/discreteMethods/dsmc/dsmcFoam/createFields.H
diff --git a/applications/solvers/dsmc/dsmcFoam/dsmcFoam.C b/applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C
similarity index 100%
rename from applications/solvers/dsmc/dsmcFoam/dsmcFoam.C
rename to applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C
diff --git a/applications/solvers/molecularDynamics/mdEquilibrationFoam/Make/files b/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/Make/files
similarity index 100%
rename from applications/solvers/molecularDynamics/mdEquilibrationFoam/Make/files
rename to applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/Make/files
diff --git a/applications/solvers/molecularDynamics/mdEquilibrationFoam/Make/options b/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/Make/options
similarity index 100%
rename from applications/solvers/molecularDynamics/mdEquilibrationFoam/Make/options
rename to applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/Make/options
diff --git a/applications/solvers/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C b/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C
similarity index 100%
rename from applications/solvers/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C
rename to applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C
diff --git a/applications/solvers/molecularDynamics/mdEquilibrationFoam/readmdEquilibrationDict.H b/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/readmdEquilibrationDict.H
similarity index 100%
rename from applications/solvers/molecularDynamics/mdEquilibrationFoam/readmdEquilibrationDict.H
rename to applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/readmdEquilibrationDict.H
diff --git a/applications/solvers/molecularDynamics/mdFoam/Make/files b/applications/solvers/discreteMethods/molecularDynamics/mdFoam/Make/files
similarity index 100%
rename from applications/solvers/molecularDynamics/mdFoam/Make/files
rename to applications/solvers/discreteMethods/molecularDynamics/mdFoam/Make/files
diff --git a/applications/solvers/molecularDynamics/mdFoam/Make/options b/applications/solvers/discreteMethods/molecularDynamics/mdFoam/Make/options
similarity index 100%
rename from applications/solvers/molecularDynamics/mdFoam/Make/options
rename to applications/solvers/discreteMethods/molecularDynamics/mdFoam/Make/options
diff --git a/applications/solvers/molecularDynamics/mdFoam/mdFoam.C b/applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C
similarity index 100%
rename from applications/solvers/molecularDynamics/mdFoam/mdFoam.C
rename to applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C
diff --git a/applications/test/volPointInterpolation/volPointInterpolationTest.C b/applications/test/volPointInterpolation/volPointInterpolationTest.C
index 07329ac349..d44dcf9fa4 100644
--- a/applications/test/volPointInterpolation/volPointInterpolationTest.C
+++ b/applications/test/volPointInterpolation/volPointInterpolationTest.C
@@ -68,8 +68,7 @@ int main(int argc, char *argv[])
         mesh
     );
 
-    pointMesh pMesh(mesh);
-    volPointInterpolation pInterp(mesh, pMesh);
+    const volPointInterpolation& pInterp = volPointInterpolation::New(mesh);
 
     pointScalarField pp(pInterp.interpolate(p));
     pp.write();
diff --git a/src/OpenFOAM/fields/pointPatchFields/constraint/processor/processorPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/constraint/processor/processorPointPatchField.C
index 02fb26d52c..c5b8c8334d 100644
--- a/src/OpenFOAM/fields/pointPatchFields/constraint/processor/processorPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/constraint/processor/processorPointPatchField.C
@@ -100,6 +100,7 @@ void processorPointPatchField<Type>::initSwapAdd(Field<Type>& pField) const
 {
     if (Pstream::parRun())
     {
+        // Get internal field into my point order
         Field<Type> pf(this->patchInternalField(pField));
 
         OPstream::write
@@ -130,11 +131,7 @@ void processorPointPatchField<Type>::swapAdd(Field<Type>& pField) const
 
         if (doTransform())
         {
-            const labelList& nonGlobalPatchPoints =
-                procPatch_.nonGlobalPatchPoints();
-
             const processorPolyPatch& ppp = procPatch_.procPolyPatch();
-            const labelListList& pointFaces = ppp.pointFaces();
             const tensorField& forwardT = ppp.forwardT();
 
             if (forwardT.size() == 1)
@@ -143,6 +140,10 @@ void processorPointPatchField<Type>::swapAdd(Field<Type>& pField) const
             }
             else
             {
+                const labelList& nonGlobalPatchPoints =
+                    procPatch_.nonGlobalPatchPoints();
+                const labelListList& pointFaces = ppp.pointFaces();
+
                 forAll(nonGlobalPatchPoints, pfi)
                 {
                     pnf[pfi] = transform
diff --git a/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/processor/processorPointPatch.C b/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/processor/processorPointPatch.C
index 2e3e001ca9..6ef8787d77 100644
--- a/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/processor/processorPointPatch.C
+++ b/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/processor/processorPointPatch.C
@@ -64,7 +64,7 @@ void Foam::processorPointPatch::initGeometry()
     }
     else
     {
-        // Slave side.  Create the reversed patch and pick up its points
+        // Slave side. Create the reversed patch and pick up its points
         // so that the order is correct
         const polyPatch& pp = patch();
 
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H
index f0bd005b47..c741b0839e 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H
@@ -330,7 +330,9 @@ public:
 
         // Addressing into mesh
 
-            //- Return labelList of mesh points in patch
+            //- Return labelList of mesh points in patch. They are constructed
+            //  walking through the faces in incremental order and not sorted
+            //  anymore.
             const labelList& meshPoints() const;
 
             //- Mesh point map.  Given the global point index find its
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C
index 9941b02459..30cffb6310 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C
@@ -67,30 +67,56 @@ calcMeshData() const
     // number of faces in the patch
     Map<label> markedPoints(4*this->size());
 
-    // if the point is used, set the mark to 1
+
+    // Important:
+    // ~~~~~~~~~~
+    // In <= 1.5 the meshPoints would be in increasing order but this gives
+    // problems in processor point synchronisation where we have to find out
+    // how the opposite side would have allocated points.
+
+    ////- 1.5 code:
+    //// if the point is used, set the mark to 1
+    //forAll (*this, faceI)
+    //{
+    //    const Face& curPoints = this->operator[](faceI);
+    //
+    //    forAll (curPoints, pointI)
+    //    {
+    //        markedPoints.insert(curPoints[pointI], -1);
+    //    }
+    //}
+    //
+    //// Create the storage and store the meshPoints.  Mesh points are
+    //// the ones marked by the usage loop above
+    //meshPointsPtr_ = new labelList(markedPoints.toc());
+    //labelList& pointPatch = *meshPointsPtr_;
+    //
+    //// Sort the list to preserve compatibility with the old ordering
+    //sort(pointPatch);
+    //
+    //// For every point in map give it its label in mesh points
+    //forAll (pointPatch, pointI)
+    //{
+    //    markedPoints.find(pointPatch[pointI])() = pointI;
+    //}
+
+    //- Unsorted version:
+    DynamicList<label> meshPoints(2*this->size());
     forAll (*this, faceI)
     {
         const Face& curPoints = this->operator[](faceI);
 
         forAll (curPoints, pointI)
         {
-            markedPoints.insert(curPoints[pointI], -1);
+            if (markedPoints.insert(curPoints[pointI], meshPoints.size()))
+            {
+                meshPoints.append(curPoints[pointI]);
+            }
         }
     }
+    // Transfer to straight list (reuses storage)
+    meshPointsPtr_ = new labelList(meshPoints, true);
 
-    // Create the storage and store the meshPoints.  Mesh points are
-    // the ones marked by the usage loop above
-    meshPointsPtr_ = new labelList(markedPoints.toc());
-    labelList& pointPatch = *meshPointsPtr_;
-
-    // Sort the list to preserve compatibility with the old ordering
-    sort(pointPatch);
-
-    // For every point in map give it its label in mesh points
-    forAll (pointPatch, pointI)
-    {
-        markedPoints.find(pointPatch[pointI])() = pointI;
-    }
 
     // Create local faces. Note that we start off from copy of original face
     // list (even though vertices are overwritten below). This is done so
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.H b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.H
index 83e7dd9c69..9ae910b5a7 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.H
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.H
@@ -38,7 +38,6 @@ SourceFiles
 
 #include "autoPtr.H"
 #include "dictionary.H"
-#include "pointField.H"
 #include "boolList.H"
 #include "wallPoint.H"
 #include "searchableSurfaces.H"
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
index aa7ecddf8b..011417c29e 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
@@ -29,7 +29,6 @@ Description
 
 #include "autoSnapDriver.H"
 #include "Time.H"
-#include "pointFields.H"
 #include "motionSmoother.H"
 #include "polyTopoChange.H"
 #include "OFstream.H"
@@ -1389,11 +1388,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::repatchToSurface
             }
         }
 
-        pointField localFaceCentres(pp.size());
-        forAll(pp, i)
-        {
-            localFaceCentres[i] = mesh.faceCentres()[pp.addressing()[i]];
-        }
+        pointField localFaceCentres(mesh.faceCentres(), pp.addressing());
 
         // Get nearest surface and region
         labelList hitSurface;
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index 24c7819d8d..d42ac05b22 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -1242,7 +1242,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::refine
 // Do refinement of consistent set of cells followed by truncation and
 // load balancing.
 Foam::autoPtr<Foam::mapDistributePolyMesh>
- Foam::meshRefinement::refineAndBalance
+Foam::meshRefinement::refineAndBalance
 (
     const string& msg,
     decompositionMethod& decomposer,
diff --git a/src/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C b/src/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C
index 9c06ccddea..856f2ba923 100644
--- a/src/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C
+++ b/src/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C
@@ -193,7 +193,7 @@ void Foam::shellSurfaces::orient()
 
     if (hasSurface)
     {
-        const point outsidePt = 2 * overallBb.span();
+        const point outsidePt = overallBb.max() + overallBb.span();
 
         //Info<< "Using point " << outsidePt << " to orient shells" << endl;
 
diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.C b/src/meshTools/searchableSurface/triSurfaceMesh.C
index 4703d69201..877f264341 100644
--- a/src/meshTools/searchableSurface/triSurfaceMesh.C
+++ b/src/meshTools/searchableSurface/triSurfaceMesh.C
@@ -120,6 +120,29 @@ const Foam::fileName& Foam::triSurfaceMesh::checkFile
 }
 
 
+bool Foam::triSurfaceMesh::addFaceToEdge
+(
+    const edge& e,
+    EdgeMap<label>& facesPerEdge
+)
+{
+    EdgeMap<label>::iterator eFnd = facesPerEdge.find(e);
+    if (eFnd != facesPerEdge.end())
+    {
+        if (eFnd() == 2)
+        {
+            return false;
+        }
+        eFnd()++;
+    }
+    else
+    {
+        facesPerEdge.insert(e, 1);
+    }
+    return true;
+}
+
+
 bool Foam::triSurfaceMesh::isSurfaceClosed() const
 {
     // Construct pointFaces. Let's hope surface has compact point
@@ -142,48 +165,41 @@ bool Foam::triSurfaceMesh::isSurfaceClosed() const
             const labelledTri& f = triSurface::operator[](pFaces[i]);
             label fp = findIndex(f, pointI);
 
-            // Forward edge
-            {
-                label p1 = f[f.fcIndex(fp)];
+            // Something weird: if I expand the code of addFaceToEdge in both
+            // below instances it gives a segmentation violation on some
+            // surfaces. Compiler (4.3.2) problem?
 
-                if (p1 > pointI)
+
+            // Forward edge
+            label nextPointI = f[f.fcIndex(fp)];
+
+            if (nextPointI > pointI)
+            {
+                bool okFace = addFaceToEdge
+                (
+                    edge(pointI, nextPointI),
+                    facesPerEdge
+                );
+
+                if (!okFace)
                 {
-                    const edge e(pointI, p1);
-                    EdgeMap<label>::iterator eFnd = facesPerEdge.find(e);
-                    if (eFnd != facesPerEdge.end())
-                    {
-                        if (eFnd() == 2)
-                        {
-                            return false;
-                        }
-                        eFnd()++;
-                    }
-                    else
-                    {
-                        facesPerEdge.insert(e, 1);
-                    }
+                    return false;
                 }
             }
             // Reverse edge
-            {
-                label p1 = f[f.rcIndex(fp)];
+            label prevPointI = f[f.rcIndex(fp)];
 
-                if (p1 > pointI)
+            if (prevPointI > pointI)
+            {
+                bool okFace = addFaceToEdge
+                (
+                    edge(pointI, prevPointI),
+                    facesPerEdge
+                );
+
+                if (!okFace)
                 {
-                    const edge e(pointI, p1);
-                    EdgeMap<label>::iterator eFnd = facesPerEdge.find(e);
-                    if (eFnd != facesPerEdge.end())
-                    {
-                        if (eFnd() == 2)
-                        {
-                            return false;
-                        }
-                        eFnd()++;
-                    }
-                    else
-                    {
-                        facesPerEdge.insert(e, 1);
-                    }
+                    return false;
                 }
             }
         }
diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.H b/src/meshTools/searchableSurface/triSurfaceMesh.H
index a4231457d5..b2ec61c14a 100644
--- a/src/meshTools/searchableSurface/triSurfaceMesh.H
+++ b/src/meshTools/searchableSurface/triSurfaceMesh.H
@@ -42,6 +42,7 @@ SourceFiles
 #include "indexedOctree.H"
 #include "treeDataTriSurface.H"
 #include "treeDataEdge.H"
+#include "EdgeMap.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -91,6 +92,13 @@ private:
             const fileName& objectName
         );
 
+        //- Helper function for isSurfaceClosed
+        static bool addFaceToEdge
+        (
+            const edge&,
+            EdgeMap<label>&
+        );
+
         //- Check whether surface is closed without calculating any permanent
         //  addressing.
         bool isSurfaceClosed() const;
diff --git a/src/meshTools/triSurface/orientedSurface/orientedSurface.C b/src/meshTools/triSurface/orientedSurface/orientedSurface.C
index 911d06747a..10e1fa8f2f 100644
--- a/src/meshTools/triSurface/orientedSurface/orientedSurface.C
+++ b/src/meshTools/triSurface/orientedSurface/orientedSurface.C
@@ -175,6 +175,51 @@ Foam::labelList Foam::orientedSurface::edgeToFace
 }
 
 
+void Foam::orientedSurface::walkSurface
+(
+    const triSurface& s,
+    const label startFaceI,
+    labelList& flipState
+)
+{
+    // List of faces that were changed in the last iteration.
+    labelList changedFaces(1, startFaceI);
+    // List of edges that were changed in the last iteration.
+    labelList changedEdges;
+
+    while(true)
+    {
+        changedEdges = faceToEdge(s, changedFaces);
+
+        if (debug)
+        {
+            Pout<< "From changedFaces:" << changedFaces.size()
+                << " to changedEdges:" << changedEdges.size()
+                << endl;
+        }
+
+        if (changedEdges.empty())
+        {
+            break;
+        }
+
+        changedFaces = edgeToFace(s, changedEdges, flipState);
+
+        if (debug)
+        {
+            Pout<< "From changedEdges:" << changedEdges.size()
+                << " to changedFaces:" << changedFaces.size()
+                << endl;
+        }
+
+        if (changedFaces.empty())
+        {
+            break;
+        }
+    }
+}
+
+
 void Foam::orientedSurface::propagateOrientation
 (
     const triSurface& s,
@@ -228,42 +273,8 @@ void Foam::orientedSurface::propagateOrientation
             << endl;
     }
 
-
-    // List of faces that were changed in the last iteration.
-    labelList changedFaces(1, nearestFaceI);
-    // List of edges that were changed in the last iteration.
-    labelList changedEdges;
-
-    while(true)
-    {
-        changedEdges = faceToEdge(s, changedFaces);
-
-        if (debug)
-        {
-            Pout<< "From changedFaces:" << changedFaces.size()
-                << " to changedEdges:" << changedEdges.size()
-                << endl;
-        }
-
-        if (changedEdges.empty())
-        {
-            break;
-        }
-
-        changedFaces = edgeToFace(s, changedEdges, flipState);
-
-        if (debug)
-        {
-            Pout<< "From changedEdges:" << changedEdges.size()
-                << " to changedFaces:" << changedFaces.size()
-                << endl;
-        }
-
-        if (changedFaces.empty())
-        {
-            break;
-        }
-    }
+    // Walk the surface from nearestFaceI, changing the flipstate.
+    walkSurface(s, nearestFaceI, flipState);
 }
 
 
@@ -352,6 +363,26 @@ bool Foam::orientedSurface::orient
     const bool orientOutside
 )
 {
+    bool anyFlipped = false;
+
+    // Do initial flipping to make triangles consistent. Otherwise if the
+    // nearest is e.g. on an edge inbetween inconsistent triangles it might
+    // make the wrong decision.
+    if (s.size() > 0)
+    {
+        // Whether face has to be flipped.
+        //      UNVISITED: unvisited
+        //      NOFLIP: no need to flip
+        //      FLIP: need to flip
+        labelList flipState(s.size(), UNVISITED);
+
+        flipState[0] = NOFLIP;
+        walkSurface(s, 0, flipState);
+
+        anyFlipped = flipSurface(s, flipState);
+    }
+
+
     // Whether face has to be flipped.
     //      UNVISITED: unvisited
     //      NOFLIP: no need to flip
@@ -410,7 +441,9 @@ bool Foam::orientedSurface::orient
     }
 
     // Now finally flip triangles according to flipState.
-    return flipSurface(s, flipState);
+    bool geomFlipped = flipSurface(s, flipState);
+
+    return anyFlipped || geomFlipped;
 }
 
 
diff --git a/src/meshTools/triSurface/orientedSurface/orientedSurface.H b/src/meshTools/triSurface/orientedSurface/orientedSurface.H
index c0b5ca8c42..812edd31bd 100644
--- a/src/meshTools/triSurface/orientedSurface/orientedSurface.H
+++ b/src/meshTools/triSurface/orientedSurface/orientedSurface.H
@@ -94,6 +94,15 @@ class orientedSurface
             labelList& flip
         );
 
+        //- Walk from face across connected faces. Change orientation to be
+        //  consistent with startFaceI.
+        static void walkSurface
+        (
+            const triSurface& s,
+            const label startFaceI,
+            labelList& flipState
+        );
+
         //- Given nearest point and face check orientation to nearest face
         //  and flip if nessecary (only marked in flipState) and propagate.
         static void propagateOrientation
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
index 5fd09bda66..017117370d 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
@@ -307,6 +307,48 @@ Foam::tmp<Foam::volScalarField> Foam::sampledIsoSurface::average
 }
 
 
+Foam::tmp<Foam::pointScalarField> Foam::sampledIsoSurface::average
+(
+    const pointMesh& pMesh,
+    const volScalarField& fld
+) const
+{
+    tmp<pointScalarField> tpointAvg
+    (
+        new pointScalarField
+        (
+            IOobject
+            (
+                "pointAvg",
+                fld.time().timeName(),
+                fld.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            pMesh,
+            dimensionedScalar("zero", dimless, scalar(0.0))
+        )
+    );
+    pointScalarField& pointAvg = tpointAvg();
+
+    for (label pointI = 0; pointI < fld.mesh().nPoints(); pointI++)
+    {
+        const labelList& pCells = fld.mesh().pointCells(pointI);
+
+        forAll(pCells, i)
+        {
+            pointAvg[pointI] += fld[pCells[i]];
+        }
+        pointAvg[pointI] /= pCells.size();
+    }
+    // Give value to calculatedFvPatchFields
+    pointAvg.correctBoundaryConditions();
+
+    return tpointAvg;
+}
+
+
 bool Foam::sampledIsoSurface::updateGeometry() const
 {
     const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
@@ -407,6 +449,7 @@ bool Foam::sampledIsoSurface::updateGeometry() const
                 (
                     *volFieldPtr_,
                     *pointFieldPtr_,
+                    //average(pointMesh::New(mesh()), *volFieldPtr_),
                     isoVal_,
                     regularise_,
                     mergeTol_
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
index 99324cdf5f..ccaa4d8750 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
@@ -124,6 +124,12 @@ class sampledIsoSurface
             const pointScalarField&
         ) const;
 
+        tmp<pointScalarField> average
+        (
+            const pointMesh&,
+            const volScalarField& fld
+        ) const;
+
         //- Create iso surface (if time has changed)
         //  Do nothing (and return false) if no update was needed
         bool updateGeometry() const;