diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C
index f68bbd476c..3ff3071a77 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,6 +35,59 @@ const Foam::scalar Foam::polyMeshTetDecomposition::minTetQuality = sqr(small);
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
+Foam::scalar Foam::polyMeshTetDecomposition::minQuality
+(
+    const polyMesh& mesh,
+    const point& cC,
+    const label fI,
+    const bool isOwner,
+    const label faceBasePtI
+)
+{
+    // Does fan decomposition of face (starting at faceBasePti) and determines
+    // min quality over all resulting tets.
+
+    const pointField& pPts = mesh.points();
+    const face& f = mesh.faces()[fI];
+    const point& tetBasePt = pPts[f[faceBasePtI]];
+
+    scalar thisBaseMinTetQuality = vGreat;
+
+    for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
+    {
+        label facePtI = (tetPtI + faceBasePtI) % f.size();
+        label otherFacePtI = f.fcIndex(facePtI);
+
+        label ptAI = -1;
+        label ptBI = -1;
+
+        if (isOwner)
+        {
+            ptAI = f[facePtI];
+            ptBI = f[otherFacePtI];
+        }
+        else
+        {
+            ptAI = f[otherFacePtI];
+            ptBI = f[facePtI];
+        }
+
+        const point& pA = pPts[ptAI];
+        const point& pB = pPts[ptBI];
+
+        tetPointRef tet(cC, tetBasePt, pA, pB);
+
+        scalar tetQuality = tet.quality();
+
+        if (tetQuality < thisBaseMinTetQuality)
+        {
+            thisBaseMinTetQuality = tetQuality;
+        }
+    }
+    return thisBaseMinTetQuality;
+}
+
+
 Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
 (
     const polyMesh& mesh,
@@ -45,7 +98,6 @@ Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
 )
 {
     const faceList& pFaces = mesh.faces();
-    const pointField& pPts = mesh.points();
     const vectorField& pC = mesh.cellCentres();
     const labelList& pOwner = mesh.faceOwner();
 
@@ -55,52 +107,12 @@ Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
 
     const point& oCc = pC[oCI];
 
-    List<scalar> tetQualities(2, 0.0);
-
     forAll(f, faceBasePtI)
     {
-        scalar thisBaseMinTetQuality = vGreat;
+        scalar minQ = minQuality(mesh, oCc, fI, true, faceBasePtI);
+        minQ = min(minQ, minQuality(mesh, nCc, fI, false, faceBasePtI));
 
-        const point& tetBasePt = pPts[f[faceBasePtI]];
-
-        for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
-        {
-            label facePtI = (tetPtI + faceBasePtI) % f.size();
-            label otherFacePtI = f.fcIndex(facePtI);
-
-            {
-                // owner cell tet
-                label ptAI = f[facePtI];
-                label ptBI = f[otherFacePtI];
-
-                const point& pA = pPts[ptAI];
-                const point& pB = pPts[ptBI];
-
-                tetPointRef tet(oCc, tetBasePt, pA, pB);
-
-                tetQualities[0] = tet.quality();
-            }
-
-            {
-                // neighbour cell tet
-                label ptAI = f[otherFacePtI];
-                label ptBI = f[facePtI];
-
-                const point& pA = pPts[ptAI];
-                const point& pB = pPts[ptBI];
-
-                tetPointRef tet(nCc, tetBasePt, pA, pB);
-
-                tetQualities[1] = tet.quality();
-            }
-
-            if (min(tetQualities) < thisBaseMinTetQuality)
-            {
-                thisBaseMinTetQuality = min(tetQualities);
-            }
-        }
-
-        if (thisBaseMinTetQuality > tol)
+        if (minQ > tol)
         {
             return faceBasePtI;
         }
@@ -140,7 +152,6 @@ Foam::label Foam::polyMeshTetDecomposition::findBasePoint
 )
 {
     const faceList& pFaces = mesh.faces();
-    const pointField& pPts = mesh.points();
     const vectorField& pC = mesh.cellCentres();
     const labelList& pOwner = mesh.faceOwner();
 
@@ -154,43 +165,9 @@ Foam::label Foam::polyMeshTetDecomposition::findBasePoint
 
     forAll(f, faceBasePtI)
     {
-        scalar thisBaseMinTetQuality = vGreat;
+        scalar quality = minQuality(mesh, cC, fI, own, faceBasePtI);
 
-        const point& tetBasePt = pPts[f[faceBasePtI]];
-
-        for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
-        {
-            label facePtI = (tetPtI + faceBasePtI) % f.size();
-            label otherFacePtI = f.fcIndex(facePtI);
-
-            label ptAI = -1;
-            label ptBI = -1;
-
-            if (own)
-            {
-                ptAI = f[facePtI];
-                ptBI = f[otherFacePtI];
-            }
-            else
-            {
-                ptAI = f[otherFacePtI];
-                ptBI = f[facePtI];
-            }
-
-            const point& pA = pPts[ptAI];
-            const point& pB = pPts[ptBI];
-
-            tetPointRef tet(cC, tetBasePt, pA, pB);
-
-            scalar tetQuality = tet.quality();
-
-            if (tetQuality < thisBaseMinTetQuality)
-            {
-                thisBaseMinTetQuality = tetQuality;
-            }
-        }
-
-        if (thisBaseMinTetQuality > tol)
+        if (quality > tol)
         {
             return faceBasePtI;
         }
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H
index 1b8b55f343..03c7b9d7a2 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -67,6 +67,17 @@ public:
 
     // Member Functions
 
+        //- Given a face and cc and starting index for triangulation determine
+        //  the worst tet quality.
+        static scalar minQuality
+        (
+            const polyMesh& mesh,
+            const point& cC,
+            const label fI,
+            const bool isOwner,
+            const label faceBasePtI
+        );
+
         //- Find the first suitable base point to use for a minimum
         //  triangle decomposition of the face, suiting owner and
         //  neighbour cells.  Finds the first base point on the face
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.C b/src/sampling/sampledSurface/isoSurface/isoSurface.C
index 94d15f1a7a..9f86c30572 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,6 +28,8 @@ License
 #include "tetMatcher.H"
 #include "tetPointRef.H"
 #include "DynamicField.H"
+#include "syncTools.H"
+#include "polyMeshTetDecomposition.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -106,7 +108,14 @@ Foam::isoSurface::cellCutType Foam::isoSurface::calcCutType
                 break;
             }
 
-            const label fp0 = mesh_.tetBasePtIs()[facei];
+            label fp0 = tetBasePtIs_[facei];
+
+            // Fall back for problem decompositions
+            if (fp0 < 0)
+            {
+                fp0 = 0;
+            }
+
             label fp = f.fcIndex(fp0);
             for (label i = 2; i < f.size(); i++)
             {
@@ -191,6 +200,184 @@ Foam::label Foam::isoSurface::calcCutTypes
 }
 
 
+Foam::scalar Foam::isoSurface::minTetQ
+(
+    const label facei,
+    const label faceBasePtI
+) const
+{
+    scalar q = polyMeshTetDecomposition::minQuality
+    (
+        mesh_,
+        mesh_.cellCentres()[mesh_.faceOwner()[facei]],
+        facei,
+        true,
+        faceBasePtI
+    );
+
+    if (mesh_.isInternalFace(facei))
+    {
+        q = min
+        (
+            q,
+            polyMeshTetDecomposition::minQuality
+            (
+                mesh_,
+                mesh_.cellCentres()[mesh_.faceNeighbour()[facei]],
+                facei,
+                false,
+                faceBasePtI
+            )
+        );
+    }
+    return q;
+}
+
+
+void Foam::isoSurface::fixTetBasePtIs()
+{
+    // Determine points used by two faces on the same cell
+    const cellList& cells = mesh_.cells();
+    const faceList& faces = mesh_.faces();
+    const labelList& faceOwner = mesh_.faceOwner();
+    const labelList& faceNeighbour = mesh_.faceNeighbour();
+
+
+    // Get face triangulation base point
+    tetBasePtIs_ = mesh_.tetBasePtIs();
+
+
+    // Pre-filter: mark all cells with illegal base points
+    labelHashSet problemCells(cells.size()/128);
+    forAll(tetBasePtIs_, facei)
+    {
+        if (tetBasePtIs_[facei] == -1)
+        {
+            problemCells.insert(faceOwner[facei]);
+            if (mesh_.isInternalFace(facei))
+            {
+                problemCells.insert(faceNeighbour[facei]);
+            }
+        }
+    }
+
+
+    label nAdapted = 0;
+
+
+    // Number of times a point occurs in a cell. Used to detect dangling
+    // vertices (count = 2)
+    Map<label> pointCount;
+
+    // Analyse problem cells for points shared by two faces only
+    forAll(cells, celli)
+    {
+        if (problemCells.found(celli))
+        {
+            const cell& cFaces = cells[celli];
+
+            pointCount.clear();
+
+            forAll(cFaces, i)
+            {
+                const label facei = cFaces[i];
+                const face& f = faces[facei];
+                forAll(f, fp)
+                {
+                    const label pointi = f[fp];
+
+                    Map<label>::iterator pointFnd = pointCount.find(pointi);
+                    if (pointFnd == pointCount.end())
+                    {
+                        pointCount.insert(pointi, 1);
+                    }
+                    else
+                    {
+                        ++pointFnd();
+                    }
+                }
+            }
+
+            // Check for any points with count 2
+            bool haveDangling = false;
+            forAllConstIter(Map<label>, pointCount, iter)
+            {
+                if (iter() == 1)
+                {
+                    FatalErrorInFunction << "point:" << iter.key()
+                        << " at:" << mesh_.points()[iter.key()]
+                        << " only used by one face" << exit(FatalError);
+                }
+                else if (iter() == 2)
+                {
+                    haveDangling = true;
+                    break;
+                }
+            }
+
+            if (haveDangling)
+            {
+                // Any point next to a dangling point should not be used
+                // as the fan base since this would cause two duplicate
+                // triangles.
+                forAll(cFaces, i)
+                {
+                    const label facei = cFaces[i];
+                    if (tetBasePtIs_[facei] == -1)
+                    {
+                        const face& f = faces[facei];
+
+                        // All the possible base points cause negative tets.
+                        // Choose the least-worst one
+                        scalar maxQ = -GREAT;
+                        label maxFp = -1;
+
+                        label prevCount = pointCount[f.last()];
+                        forAll(f, fp)
+                        {
+                            label nextCount = pointCount[f[f.fcIndex(fp)]];
+
+                            if (prevCount > 2 && nextCount > 2)
+                            {
+                                const scalar q = minTetQ(facei, fp);
+                                if (q > maxQ)
+                                {
+                                    maxQ = q;
+                                    maxFp = fp;
+                                }
+                            }
+                            prevCount = pointCount[f[fp]];
+                        }
+
+                        if (maxFp != -1)
+                        {
+                            // Least worst base point
+                            tetBasePtIs_[facei] = maxFp;
+                        }
+                        else
+                        {
+                            // No point found on face that would not result
+                            // in some duplicate triangle. Very rare. Do what?
+                            tetBasePtIs_[facei] = 0;
+                        }
+
+                        nAdapted++;
+                    }
+                }
+            }
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "isoSurface : adapted starting point of triangulation on "
+            << nAdapted << " faces." << endl;
+    }
+
+    syncTools::syncFaceList(mesh_, tetBasePtIs_, maxEqOp<label>());
+}
+
+
 Foam::label Foam::isoSurface::generatePoint
 (
     const label facei,
@@ -616,7 +803,6 @@ void Foam::isoSurface::generateTriPoints
 
 void Foam::isoSurface::generateTriPoints
 (
-    const polyMesh& mesh,
     const label celli,
     const bool isTet,
 
@@ -629,7 +815,10 @@ void Foam::isoSurface::generateTriPoints
     DynamicList<label>& faceLabels
 ) const
 {
-    const cell& cFaces = mesh.cells()[celli];
+    const faceList& faces = mesh_.faces();
+    const labelList& faceOwner = mesh_.faceOwner();
+    const pointField& points = mesh_.points();
+    const cell& cFaces = mesh_.cells()[celli];
 
     if (isTet)
     {
@@ -637,10 +826,10 @@ void Foam::isoSurface::generateTriPoints
         // tet points and values
 
         label facei = cFaces[0];
-        const face& f0 = mesh_.faces()[facei];
+        const face& f0 = faces[facei];
 
         // Get the other point
-        const face& f1 = mesh_.faces()[cFaces[1]];
+        const face& f1 = faces[cFaces[1]];
         label oppositeI = -1;
         forAll(f1, fp)
         {
@@ -657,17 +846,17 @@ void Foam::isoSurface::generateTriPoints
         label p2 = f0[2];
         FixedList<bool, 6> edgeIsDiag(false);
 
-        if (mesh.faceOwner()[facei] == celli)
+        if (faceOwner[facei] == celli)
         {
             Swap(p1, p2);
         }
 
         tetPointRef tet
         (
-            mesh.points()[p0],
-            mesh.points()[p1],
-            mesh.points()[p2],
-            mesh.points()[oppositeI]
+            points[p0],
+            points[p1],
+            points[p2],
+            points[oppositeI]
         );
 
         label startTrii = verts.size();
@@ -683,10 +872,10 @@ void Foam::isoSurface::generateTriPoints
             }),
             FixedList<point, 4>
             ({
-                mesh.points()[p0],
-                mesh.points()[p1],
-                mesh.points()[p2],
-                mesh.points()[oppositeI]
+                points[p0],
+                points[p1],
+                points[p2],
+                points[oppositeI]
             }),
             FixedList<label, 4>
             ({
@@ -712,16 +901,18 @@ void Foam::isoSurface::generateTriPoints
     }
     else
     {
+        const pointField& cellCentres = mesh_.cellCentres();
+
         forAll(cFaces, cFacei)
         {
             label facei = cFaces[cFacei];
-            const face& f = mesh.faces()[facei];
+            const face& f = faces[facei];
 
-            label fp0 = mesh.tetBasePtIs()[facei];
+            label fp0 = tetBasePtIs_[facei];
 
             label startTrii = verts.size();
 
-            // Skip undefined tets
+            // Fallback
             if (fp0 < 0)
             {
                 fp0 = 0;
@@ -737,7 +928,7 @@ void Foam::isoSurface::generateTriPoints
                 label p0 = f[fp0];
                 label p1 = f[fp];
                 label p2 = f[nextFp];
-                if (mesh.faceOwner()[facei] == celli)
+                if (faceOwner[facei] == celli)
                 {
                     Swap(p1, p2);
                     if (i != 2) edgeIsDiag[1] = true;
@@ -751,10 +942,10 @@ void Foam::isoSurface::generateTriPoints
 
                 tetPointRef tet
                 (
-                    mesh.points()[p0],
-                    mesh.points()[p1],
-                    mesh.points()[p2],
-                    mesh.cellCentres()[celli]
+                    points[p0],
+                    points[p1],
+                    points[p2],
+                    cellCentres[celli]
                 );
 
                 generateTriPoints
@@ -769,17 +960,17 @@ void Foam::isoSurface::generateTriPoints
                     }),
                     FixedList<point, 4>
                     ({
-                        mesh.points()[p0],
-                        mesh.points()[p1],
-                        mesh.points()[p2],
-                        mesh.cellCentres()[celli]
+                        points[p0],
+                        points[p1],
+                        points[p2],
+                        cellCentres[celli]
                     }),
                     FixedList<label, 4>
                     ({
                         p0,
                         p1,
                         p2,
-                        mesh.nPoints()+celli
+                        mesh_.nPoints()+celli
                     }),
                     edgeIsDiag,
 
@@ -983,6 +1174,8 @@ Foam::isoSurface::isoSurface
         Pout<< "isoSurface : iso:" << iso_ << " filter:" << filter << endl;
     }
 
+    fixTetBasePtIs();
+
     tetMatcher tet;
 
     // Determine if any cut through cell
@@ -1016,7 +1209,6 @@ Foam::isoSurface::isoSurface
         {
             generateTriPoints
             (
-                mesh,
                 celli,
                 tet.isA(mesh_, celli),
 
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.H b/src/sampling/sampledSurface/isoSurface/isoSurface.H
index 72887d4972..77437b8392 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.H
@@ -91,6 +91,9 @@ private:
         //- Iso value
         const scalar iso_;
 
+        //- Corrected version of tetBasePtIs
+        labelList tetBasePtIs_;
+
         //- Per point: originating mesh vertex/cc. See encoding above
         edgeList pointToVerts_;
 
@@ -103,6 +106,14 @@ private:
 
     // Private Member Functions
 
+        scalar minTetQ
+        (
+            const label facei,
+            const label faceBasePtI
+        ) const;
+
+        void fixTetBasePtIs();
+
         //- Does any edge of triangle cross iso value?
         bool isTriCut
         (
@@ -157,7 +168,6 @@ private:
         //- Generate triangles from cell
         void generateTriPoints
         (
-            const polyMesh& mesh,
             const label celli,
             const bool isTet,
 
@@ -230,7 +240,7 @@ public:
             return pointToFace_;
         }
 
-        //- Per point: originating mesh vertex/cc. See encoding above·
+        //- Per point: originating mesh vertex/cc. See encoding above
         const edgeList& pointToVerts() const
         {
             return pointToVerts_;