diff --git a/applications/utilities/postProcessing/sampling/sample/sampleDict b/applications/utilities/postProcessing/sampling/sample/sampleDict
index 191e006d1a..d0ce98158c 100644
--- a/applications/utilities/postProcessing/sampling/sample/sampleDict
+++ b/applications/utilities/postProcessing/sampling/sample/sampleDict
@@ -255,6 +255,7 @@ surfaces
         distance        0.0;
 
         interpolate     false;
+        regularise      false;       // Optional: do not simplify
         // mergeTol        1e-10;    // Optional: fraction of mesh bounding box
                                      // to merge points (default=1e-6)
     }
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
index baa36b8e24..53f99b67c0 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
@@ -287,16 +287,16 @@ Foam::pointIndexHit Foam::isoSurfaceCell::collapseSurface
 
         if (shared[0] != -1)
         {
-            //vector n0 = tri0.normal(localPoints);
-            //n0 /= mag(n0);
-            //vector n1 = tri1.normal(localPoints);
-            //n1 /= mag(n1);
-            //
-            //if ((n0 & n1) < 0)
-            //{
-            //    // Too big an angle. Do not simplify.
-            //}
-            //else
+            vector n0 = tri0.normal(localPoints);
+            n0 /= mag(n0);
+            vector n1 = tri1.normal(localPoints);
+            n1 /= mag(n1);
+
+            if ((n0 & n1) < 0)
+            {
+                // Too big an angle. Do not simplify.
+            }
+            else
             {
                 info.setPoint
                 (
@@ -332,18 +332,18 @@ Foam::pointIndexHit Foam::isoSurfaceCell::collapseSurface
         if (nZones == 1)
         {
             // Check that all normals make a decent angle
-            //scalar minCos = GREAT;
-            //const vector& n0 = surf.faceNormals()[0];
-            //for (label i = 1; i < surf.size(); i++)
-            //{
-            //    scalar cosAngle = (n0 & surf.faceNormals()[i]);
-            //    if (cosAngle < minCos)
-            //    {
-            //        minCos = cosAngle;
-            //    }
-            //}
-            //
-            //if (minCos > 0)
+            scalar minCos = GREAT;
+            const vector& n0 = surf.faceNormals()[0];
+            for (label i = 1; i < surf.size(); i++)
+            {
+                scalar cosAngle = (n0 & surf.faceNormals()[i]);
+                if (cosAngle < minCos)
+                {
+                    minCos = cosAngle;
+                }
+            }
+
+            if (minCos > 0)
             {
                 info.setPoint(calcCentre(surf));
                 info.setHit();
@@ -443,7 +443,7 @@ void Foam::isoSurfaceCell::calcSnappedCc
                     label faceI = cFaces[cFaceI];
                     const face& f = mesh_.faces()[faceI];
 
-                    // Do a tetrahedrisation. Each face to cc becomes pyr.
+                    // Do a tetrahedralisation. Each face to cc becomes pyr.
                     // Each pyr gets split into tets by diagonalisation
                     // of face.
 
@@ -467,14 +467,18 @@ void Foam::isoSurfaceCell::calcSnappedCc
                          && (s[2] >= 0.0 && s[2] <= 0.5)
                         )
                         {
-                            if (mesh_.faceOwner()[faceI] == cellI)
+                            if
+                            (
+                                (mesh_.faceOwner()[faceI] == cellI)
+                             == (cVal >= pVals[tri[0]])
+                            )
                             {
                                 localTris.append
                                 (
                                     labelledTri
                                     (
-                                        pointToLocal[tri[0]],
                                         pointToLocal[tri[1]],
+                                        pointToLocal[tri[0]],
                                         pointToLocal[tri[2]],
                                         0
                                     )
@@ -486,8 +490,8 @@ void Foam::isoSurfaceCell::calcSnappedCc
                                 (
                                     labelledTri
                                     (
-                                        pointToLocal[tri[1]],
                                         pointToLocal[tri[0]],
+                                        pointToLocal[tri[1]],
                                         pointToLocal[tri[2]],
                                         0
                                     )
@@ -575,7 +579,11 @@ void Foam::isoSurfaceCell::genPointTris
             point p1 = (1.0-s[1])*pts[pointI] + s[1]*pts[c];
             point p2 = (1.0-s[2])*pts[pointI] + s[2]*cc[cellI];
 
-            if (mesh_.faceOwner()[faceI] == cellI)
+            if
+            (
+                (mesh_.faceOwner()[faceI] == cellI)
+             == (pointValues[pointI] > cellValues[cellI])
+            )
             {
                 localTriPoints.append(p0);
                 localTriPoints.append(p1);
@@ -754,7 +762,7 @@ void Foam::isoSurfaceCell::calcSnappedPoint
         }
 
 
-        // Loop over all pointCells (by using pointFaces)
+        // Do a pointCells walk (by using pointFaces)
 
         localPointCells.clear();
         localTriPoints.clear();
@@ -766,6 +774,8 @@ void Foam::isoSurfaceCell::calcSnappedPoint
 
             if (isTet.get(own) == 1)
             {
+                // Since tets have no cell centre to include make sure
+                // we only generate a triangle once per point.
                 if (localPointCells.insert(own))
                 {
                     genPointTris(pVals, pointI, faceI, own, localTriPoints);
@@ -848,20 +858,19 @@ void Foam::isoSurfaceCell::calcSnappedPoint
 
             if (nZones == 1)
             {
-                //// Check that all normals make a decent angle
-                //scalar minCos = GREAT;
-                //const vector& n0 = surf.faceNormals()[0];
-                //for (label i = 1; i < surf.size(); i++)
-                //{
-                //    const vector& n = surf.faceNormals()[i];
-                //    scalar cosAngle = (n0 & n);
-                //    if (cosAngle < minCos)
-                //    {
-                //        minCos = cosAngle;
-                //    }
-                //}
-                //
-                //if (minCos > 0)
+                // Check that all normals make a decent angle
+                scalar minCos = GREAT;
+                const vector& n0 = surf.faceNormals()[0];
+                for (label i = 1; i < surf.size(); i++)
+                {
+                    const vector& n = surf.faceNormals()[i];
+                    scalar cosAngle = (n0 & n);
+                    if (cosAngle < minCos)
+                    {
+                        minCos = cosAngle;
+                    }
+                }
+                if (minCos > 0)
                 {
                     collapsedPoint[pointI] = calcCentre(surf);
                 }
@@ -882,7 +891,9 @@ void Foam::isoSurfaceCell::calcSnappedPoint
 
     forAll(collapsedPoint, pointI)
     {
-        if (collapsedPoint[pointI] != greatPoint)
+        // Cannot do == comparison since might be transformed so have
+        // truncation errors.
+        if (magSqr(collapsedPoint[pointI]) < 0.5*magSqr(greatPoint))
         {
             snappedPoint[pointI] = snappedPoints.size();
             snappedPoints.append(collapsedPoint[pointI]);
@@ -1208,142 +1219,142 @@ void Foam::isoSurfaceCell::calcAddressing
 }
 
 
-void Foam::isoSurfaceCell::walkOrientation
-(
-    const triSurface& surf,
-    const List<FixedList<label, 3> >& faceEdges,
-    const labelList& edgeFace0,
-    const labelList& edgeFace1,
-    const label seedTriI,
-    labelList& flipState
-)
-{
-    // Do walk for consistent orientation.
-    DynamicList<label> changedFaces(surf.size());
-
-    changedFaces.append(seedTriI);
-
-    while (changedFaces.size())
-    {
-        DynamicList<label> newChangedFaces(changedFaces.size());
-
-        forAll(changedFaces, i)
-        {
-            label triI = changedFaces[i];
-            const labelledTri& tri = surf[triI];
-            const FixedList<label, 3>& fEdges = faceEdges[triI];
-
-            forAll(fEdges, fp)
-            {
-                label edgeI = fEdges[fp];
-
-                // my points:
-                label p0 = tri[fp];
-                label p1 = tri[tri.fcIndex(fp)];
-
-                label nbrI =
-                (
-                    edgeFace0[edgeI] != triI
-                  ? edgeFace0[edgeI]
-                  : edgeFace1[edgeI]
-                );
-
-                if (nbrI != -1 && flipState[nbrI] == -1)
-                {
-                    const labelledTri& nbrTri = surf[nbrI];
-
-                    // nbr points
-                    label nbrFp = findIndex(nbrTri, p0);
-                    label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
-
-                    bool sameOrientation = (p1 == nbrP1);
-
-                    if (flipState[triI] == 0)
-                    {
-                        flipState[nbrI] = (sameOrientation ? 0 : 1);
-                    }
-                    else
-                    {
-                        flipState[nbrI] = (sameOrientation ? 1 : 0);
-                    }
-                    newChangedFaces.append(nbrI);
-                }
-            }
-        }
-
-        changedFaces.transfer(newChangedFaces);
-    }
-}
-
-
-void Foam::isoSurfaceCell::orientSurface
-(
-    triSurface& surf,
-    const List<FixedList<label, 3> >& faceEdges,
-    const labelList& edgeFace0,
-    const labelList& edgeFace1,
-    const Map<labelList>& edgeFacesRest
-)
-{
-    // -1 : unvisited
-    //  0 : leave as is
-    //  1 : flip
-    labelList flipState(surf.size(), -1);
-
-    label seedTriI = 0;
-
-    while (true)
-    {
-        // Find first unvisited triangle
-        for
-        (
-            ;
-            seedTriI < surf.size() && flipState[seedTriI] != -1;
-            seedTriI++
-        )
-        {}
-
-        if (seedTriI == surf.size())
-        {
-            break;
-        }
-
-        // Note: Determine orientation of seedTriI?
-        // for now assume it is ok
-        flipState[seedTriI] = 0;
-
-        walkOrientation
-        (
-            surf,
-            faceEdges,
-            edgeFace0,
-            edgeFace1,
-            seedTriI,
-            flipState
-        );
-    }
-
-    // Do actual flipping
-    surf.clearOut();
-    forAll(surf, triI)
-    {
-        if (flipState[triI] == 1)
-        {
-            labelledTri tri(surf[triI]);
-
-            surf[triI][0] = tri[0];
-            surf[triI][1] = tri[2];
-            surf[triI][2] = tri[1];
-        }
-        else if (flipState[triI] == -1)
-        {
-            FatalErrorIn
-            (
-                "isoSurfaceCell::orientSurface(triSurface&, const label)"
-            )   << "problem" << abort(FatalError);
-        }
-    }
-}
+//void Foam::isoSurfaceCell::walkOrientation
+//(
+//    const triSurface& surf,
+//    const List<FixedList<label, 3> >& faceEdges,
+//    const labelList& edgeFace0,
+//    const labelList& edgeFace1,
+//    const label seedTriI,
+//    labelList& flipState
+//)
+//{
+//    // Do walk for consistent orientation.
+//    DynamicList<label> changedFaces(surf.size());
+//
+//    changedFaces.append(seedTriI);
+//
+//    while (changedFaces.size())
+//    {
+//        DynamicList<label> newChangedFaces(changedFaces.size());
+//
+//        forAll(changedFaces, i)
+//        {
+//            label triI = changedFaces[i];
+//            const labelledTri& tri = surf[triI];
+//            const FixedList<label, 3>& fEdges = faceEdges[triI];
+//
+//            forAll(fEdges, fp)
+//            {
+//                label edgeI = fEdges[fp];
+//
+//                // my points:
+//                label p0 = tri[fp];
+//                label p1 = tri[tri.fcIndex(fp)];
+//
+//                label nbrI =
+//                (
+//                    edgeFace0[edgeI] != triI
+//                  ? edgeFace0[edgeI]
+//                  : edgeFace1[edgeI]
+//                );
+//
+//                if (nbrI != -1 && flipState[nbrI] == -1)
+//                {
+//                    const labelledTri& nbrTri = surf[nbrI];
+//
+//                    // nbr points
+//                    label nbrFp = findIndex(nbrTri, p0);
+//                    label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
+//
+//                    bool sameOrientation = (p1 == nbrP1);
+//
+//                    if (flipState[triI] == 0)
+//                    {
+//                        flipState[nbrI] = (sameOrientation ? 0 : 1);
+//                    }
+//                    else
+//                    {
+//                        flipState[nbrI] = (sameOrientation ? 1 : 0);
+//                    }
+//                    newChangedFaces.append(nbrI);
+//                }
+//            }
+//        }
+//
+//        changedFaces.transfer(newChangedFaces);
+//    }
+//}
+//
+//
+//void Foam::isoSurfaceCell::orientSurface
+//(
+//    triSurface& surf,
+//    const List<FixedList<label, 3> >& faceEdges,
+//    const labelList& edgeFace0,
+//    const labelList& edgeFace1,
+//    const Map<labelList>& edgeFacesRest
+//)
+//{
+//    // -1 : unvisited
+//    //  0 : leave as is
+//    //  1 : flip
+//    labelList flipState(surf.size(), -1);
+//
+//    label seedTriI = 0;
+//
+//    while (true)
+//    {
+//        // Find first unvisited triangle
+//        for
+//        (
+//            ;
+//            seedTriI < surf.size() && flipState[seedTriI] != -1;
+//            seedTriI++
+//        )
+//        {}
+//
+//        if (seedTriI == surf.size())
+//        {
+//            break;
+//        }
+//
+//        // Note: Determine orientation of seedTriI?
+//        // for now assume it is ok
+//        flipState[seedTriI] = 0;
+//
+//        walkOrientation
+//        (
+//            surf,
+//            faceEdges,
+//            edgeFace0,
+//            edgeFace1,
+//            seedTriI,
+//            flipState
+//        );
+//    }
+//
+//    // Do actual flipping
+//    surf.clearOut();
+//    forAll(surf, triI)
+//    {
+//        if (flipState[triI] == 1)
+//        {
+//            labelledTri tri(surf[triI]);
+//
+//            surf[triI][0] = tri[0];
+//            surf[triI][1] = tri[2];
+//            surf[triI][2] = tri[1];
+//        }
+//        else if (flipState[triI] == -1)
+//        {
+//            FatalErrorIn
+//            (
+//                "isoSurfaceCell::orientSurface(triSurface&, const label)"
+//            )   << "problem" << abort(FatalError);
+//        }
+//    }
+//}
 
 
 // Checks if triangle is connected through edgeI only.
@@ -1625,7 +1636,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
     (
         stitchTriPoints
         (
-            true,               // check for duplicate tris
+            regularise,         // check for duplicate tris
             triPoints,
             triPointMergeMap_,  // unmerged to merged point
             triMap
@@ -1657,7 +1668,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
     }
 
 
-    if (false)
+    if (regularise)
     {
         List<FixedList<label, 3> > faceEdges;
         labelList edgeFace0, edgeFace1;
@@ -1717,7 +1728,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
             inplaceRenumber(reversePointMap, triPointMergeMap_);
         }
 
-        orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
+        //orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
     }
 }
 
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
index fa157a517c..9c3c2990c8 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
@@ -214,15 +214,19 @@ class isoSurfaceCell
         void generateTriPoints
         (
             const DynamicList<Type>& snapped,
+            const scalar isoVal0,
             const scalar s0,
             const Type& p0,
             const label p0Index,
+            const scalar isoVal1,
             const scalar s1,
             const Type& p1,
             const label p1Index,
+            const scalar isoVal2,
             const scalar s2,
             const Type& p2,
             const label p2Index,
+            const scalar isoVal3,
             const scalar s3,
             const Type& p3,
             const label p3Index,
@@ -267,26 +271,26 @@ class isoSurfaceCell
             Map<labelList>& edgeFacesRest
         ) const;
 
-        //- Determine orientation
-        static void walkOrientation
-        (
-            const triSurface& surf,
-            const List<FixedList<label, 3> >& faceEdges,
-            const labelList& edgeFace0,
-            const labelList& edgeFace1,
-            const label seedTriI,
-            labelList& flipState
-        );
+        ////- Determine orientation
+        //static void walkOrientation
+        //(
+        //    const triSurface& surf,
+        //    const List<FixedList<label, 3> >& faceEdges,
+        //    const labelList& edgeFace0,
+        //    const labelList& edgeFace1,
+        //    const label seedTriI,
+        //    labelList& flipState
+        //);
 
-        //- Orient surface
-        static void orientSurface
-        (
-            triSurface&,
-            const List<FixedList<label, 3> >& faceEdges,
-            const labelList& edgeFace0,
-            const labelList& edgeFace1,
-            const Map<labelList>& edgeFacesRest
-        );
+        ////- Orient surface
+        //static void orientSurface
+        //(
+        //    triSurface&,
+        //    const List<FixedList<label, 3> >& faceEdges,
+        //    const labelList& edgeFace0,
+        //    const labelList& edgeFace1,
+        //    const Map<labelList>& edgeFacesRest
+        //);
 
         //- Is triangle (given by 3 edges) not fully connected?
         static bool danglingTriangle
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
index 9c476d62cd..8b933f5aab 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
@@ -76,23 +76,27 @@ void Foam::isoSurfaceCell::generateTriPoints
 (
     const DynamicList<Type>& snapped,
 
+    const scalar isoVal0,
     const scalar s0,
     const Type& p0,
     const label p0Index,
 
+    const scalar isoVal1,
     const scalar s1,
     const Type& p1,
     const label p1Index,
 
+    const scalar isoVal2,
     const scalar s2,
     const Type& p2,
     const label p2Index,
 
+    const scalar isoVal3,
     const scalar s3,
     const Type& p3,
     const label p3Index,
 
-    DynamicList<Type>& points
+    DynamicList<Type>& pts
 ) const
 {
     int triIndex = 0;
@@ -122,77 +126,160 @@ void Foam::isoSurfaceCell::generateTriPoints
 
         case 0x0E:
         case 0x01:
-            points.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
-            points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
-            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+        {
+            // 0 is common point. Orient such that normal points in positive
+            // gradient direction
+            if (isoVal0 >= isoVal1)
+            {
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+            }
+            else
+            {
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+            }
+        }
         break;
 
         case 0x0D:
         case 0x02:
-            points.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
-            points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
-            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+        {
+            // 1 is common point
+            if (isoVal1 >= isoVal0)
+            {
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+            }
+            else
+            {
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+            }
+        }
         break;
 
         case 0x0C:
         case 0x03:
         {
-            Type tp1 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
-            Type tp2 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
+            Type s02 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
+            Type s13 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
 
-            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
-            points.append(tp1);
-            points.append(tp2);
-            points.append(tp2);
-            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-            points.append(tp1);
+            if (isoVal0 >= isoVal3)
+            {
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+                pts.append(s02);
+                pts.append(s13);
+                pts.append(s13);
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+                pts.append(s02);
+            }
+            else
+            {
+                pts.append(s02);
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+                pts.append(s13);
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+                pts.append(s13);
+                pts.append(s02);
+            }
         }
         break;
 
         case 0x0B:
         case 0x04:
         {
-            points.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
-            points.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
-            points.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
+            // 2 is common point
+            if (isoVal2 >= isoVal0)
+            {
+                pts.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
+                pts.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
+                pts.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
+            }
+            else
+            {
+                pts.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
+                pts.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
+                pts.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
+            }
         }
         break;
 
         case 0x0A:
         case 0x05:
         {
-            Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
-            Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+            Type s01 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+            Type s23 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
 
-            points.append(tp0);
-            points.append(tp1);
-            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
-            points.append(tp0);
-            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-            points.append(tp1);
+            if (isoVal3 >= isoVal0)
+            {
+                pts.append(s01);
+                pts.append(s23);
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+                pts.append(s01);
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+                pts.append(s23);
+            }
+            else
+            {
+                pts.append(s23);
+                pts.append(s01);
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+                pts.append(s01);
+                pts.append(s23);
+            }
         }
         break;
 
         case 0x09:
         case 0x06:
         {
-            Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
-            Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+            Type s01 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+            Type s23 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
 
-            points.append(tp0);
-            points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
-            points.append(tp1);
-            points.append(tp0);
-            points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
-            points.append(tp1);
+            if (isoVal3 >= isoVal1)
+            {
+                pts.append(s01);
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
+                pts.append(s23);
+                pts.append(s01);
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
+                pts.append(s23);
+            }
+            else
+            {
+                pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
+                pts.append(s01);
+                pts.append(s23);
+                pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
+                pts.append(s01);
+                pts.append(s23);
+            }
         }
         break;
 
         case 0x07:
         case 0x08:
-            points.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
-            points.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
-            points.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
+        {
+            // 3 is common point
+            if (isoVal3 >= isoVal0)
+            {
+                pts.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
+                pts.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
+                pts.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
+            }
+            else
+            {
+                pts.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
+                pts.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
+                pts.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
+            }
+        }
         break;
     }
 }
@@ -252,18 +339,23 @@ void Foam::isoSurfaceCell::generateTriPoints
                     generateTriPoints
                     (
                         snappedPoints,
+
+                        pVals_[f0[1]],
                         pVals[f0[1]],
                         pCoords[f0[1]],
                         snappedPoint[f0[1]],
 
+                        pVals_[f0[0]],
                         pVals[f0[0]],
                         pCoords[f0[0]],
                         snappedPoint[f0[0]],
 
+                        pVals_[f0[2]],
                         pVals[f0[2]],
                         pCoords[f0[2]],
                         snappedPoint[f0[2]],
 
+                        pVals_[oppositeI],
                         pVals[oppositeI],
                         pCoords[oppositeI],
                         snappedPoint[oppositeI],
@@ -276,18 +368,23 @@ void Foam::isoSurfaceCell::generateTriPoints
                     generateTriPoints
                     (
                         snappedPoints,
+
+                        pVals_[f0[0]],
                         pVals[f0[0]],
                         pCoords[f0[0]],
                         snappedPoint[f0[0]],
 
+                        pVals_[f0[1]],
                         pVals[f0[1]],
                         pCoords[f0[1]],
                         snappedPoint[f0[1]],
 
+                        pVals_[f0[2]],
                         pVals[f0[2]],
                         pCoords[f0[2]],
                         snappedPoint[f0[2]],
 
+                        pVals_[oppositeI],
                         pVals[oppositeI],
                         pCoords[oppositeI],
                         snappedPoint[oppositeI],
@@ -321,18 +418,22 @@ void Foam::isoSurfaceCell::generateTriPoints
                             (
                                 snappedPoints,
 
+                                pVals_[tri[1]],
                                 pVals[tri[1]],
                                 pCoords[tri[1]],
                                 snappedPoint[tri[1]],
 
+                                pVals_[tri[0]],
                                 pVals[tri[0]],
                                 pCoords[tri[0]],
                                 snappedPoint[tri[0]],
 
+                                pVals_[tri[2]],
                                 pVals[tri[2]],
                                 pCoords[tri[2]],
                                 snappedPoint[tri[2]],
 
+                                cVals_[cellI],
                                 cVals[cellI],
                                 cCoords[cellI],
                                 snappedCc[cellI],
@@ -346,18 +447,22 @@ void Foam::isoSurfaceCell::generateTriPoints
                             (
                                 snappedPoints,
 
+                                pVals_[tri[0]],
                                 pVals[tri[0]],
                                 pCoords[tri[0]],
                                 snappedPoint[tri[0]],
 
+                                pVals_[tri[1]],
                                 pVals[tri[1]],
                                 pCoords[tri[1]],
                                 snappedPoint[tri[1]],
 
+                                pVals_[tri[2]],
                                 pVals[tri[2]],
                                 pCoords[tri[2]],
                                 snappedPoint[tri[2]],
 
+                                cVals_[cellI],
                                 cVals[cellI],
                                 cCoords[cellI],
                                 snappedCc[cellI],