From 0b99b5af1a561eb178d8b258e615c3843c999548 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://openfoam.org>
Date: Tue, 13 Nov 2018 12:14:13 +0000
Subject: [PATCH] sampling: Removed deprecated isoSurfaceCell

---
 .../foamyHexMeshBackgroundMesh/Make/options   |    3 -
 .../foamyHexMeshBackgroundMesh.C              |    5 +-
 src/sampling/Make/files                       |    2 -
 .../distanceSurface/distanceSurface.C         |   51 +-
 .../distanceSurface/distanceSurface.H         |   18 +-
 .../distanceSurfaceTemplates.C                |   50 +-
 .../sampledSurface/isoSurface/isoSurface.H    |    6 -
 .../isoSurface/isoSurfaceCell.C               | 1589 -----------------
 .../isoSurface/isoSurfaceCell.H               |  352 ----
 .../isoSurface/isoSurfaceCellTemplates.C      |  566 ------
 .../isoSurface/sampledIsoSurfaceCell.C        |  368 ----
 .../isoSurface/sampledIsoSurfaceCell.H        |  246 ---
 .../sampledIsoSurfaceCellTemplates.C          |   88 -
 .../sampledCuttingPlane/sampledCuttingPlane.C |    9 -
 .../sampledCuttingPlane/sampledCuttingPlane.H |    4 -
 15 files changed, 25 insertions(+), 3332 deletions(-)
 delete mode 100644 src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
 delete mode 100644 src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
 delete mode 100644 src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
 delete mode 100644 src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
 delete mode 100644 src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
 delete mode 100644 src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C

diff --git a/applications/utilities/mesh/generation/foamyMesh/foamyHexMeshBackgroundMesh/Make/options b/applications/utilities/mesh/generation/foamyMesh/foamyHexMeshBackgroundMesh/Make/options
index 701bea14b7..61661dcca2 100644
--- a/applications/utilities/mesh/generation/foamyMesh/foamyHexMeshBackgroundMesh/Make/options
+++ b/applications/utilities/mesh/generation/foamyMesh/foamyHexMeshBackgroundMesh/Make/options
@@ -20,9 +20,6 @@ EXE_INC = \
 
 EXE_LIBS = \
     $(CGAL_LIBS) \
-    -lboost_thread \
-    -lmpfr \
-    -lgmp \
     -lconformalVoronoiMesh \
     -ldecompositionMethods /* -L$(FOAM_LIBBIN)/dummy -lscotchDecomp */ \
     -ltriSurface \
diff --git a/applications/utilities/mesh/generation/foamyMesh/foamyHexMeshBackgroundMesh/foamyHexMeshBackgroundMesh.C b/applications/utilities/mesh/generation/foamyMesh/foamyHexMeshBackgroundMesh/foamyHexMeshBackgroundMesh.C
index 7cec301dc7..a360396873 100644
--- a/applications/utilities/mesh/generation/foamyMesh/foamyHexMeshBackgroundMesh/foamyHexMeshBackgroundMesh.C
+++ b/applications/utilities/mesh/generation/foamyMesh/foamyHexMeshBackgroundMesh/foamyHexMeshBackgroundMesh.C
@@ -41,7 +41,7 @@ Description
 #include "cellShape.H"
 #include "cellModeller.H"
 #include "DynamicField.H"
-#include "isoSurfaceCell.H"
+#include "isoSurface.H"
 #include "vtkSurfaceWriter.H"
 #include "syncTools.H"
 
@@ -716,9 +716,8 @@ int main(int argc, char *argv[])
             }
         }
 
-        isoSurfaceCell iso
+        isoSurface iso
         (
-            fvm,
             cellDistance,
             pointDistance,
             0,      // distance,
diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index 31fe506ae5..5c61d41776 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -25,8 +25,6 @@ sampledSurface/sampledPatchInternalField/sampledPatchInternalField.C
 sampledSurface/sampledPlane/sampledPlane.C
 sampledSurface/isoSurface/isoSurface.C
 sampledSurface/isoSurface/sampledIsoSurface.C
-sampledSurface/isoSurface/isoSurfaceCell.C
-sampledSurface/isoSurface/sampledIsoSurfaceCell.C
 sampledSurface/distanceSurface/distanceSurface.C
 sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
 sampledSurface/sampledSurface/sampledSurface.C
diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
index b11bb8d8f1..ec0bd17127 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
@@ -50,7 +50,6 @@ void Foam::distanceSurface::createGeometry()
 
     // Clear any stored topologies
     facesPtr_.clear();
-    isoSurfCellPtr_.clear();
     isoSurfPtr_.clear();
 
     // Clear derived data
@@ -265,33 +264,16 @@ void Foam::distanceSurface::createGeometry()
 
 
     //- Direct from cell field and point field.
-    if (cell_)
-    {
-        isoSurfCellPtr_.reset
+    isoSurfPtr_.reset
+    (
+        new isoSurface
         (
-            new isoSurfaceCell
-            (
-                fvm,
-                cellDistance,
-                pointDistance_,
-                distance_,
-                regularise_
-            )
-        );
-    }
-    else
-    {
-        isoSurfPtr_.reset
-        (
-            new isoSurface
-            (
-                cellDistance,
-                pointDistance_,
-                distance_,
-                regularise_
-            )
-        );
-    }
+            cellDistance,
+            pointDistance_,
+            distance_,
+            regularise_
+        )
+    );
 
     if (debug)
     {
@@ -330,24 +312,13 @@ Foam::distanceSurface::distanceSurface
     ),
     distance_(readScalar(dict.lookup("distance"))),
     signed_(readBool(dict.lookup("signed"))),
-    cell_(dict.lookupOrDefault("cell", true)),
     regularise_(dict.lookupOrDefault("regularise", true)),
     average_(dict.lookupOrDefault("average", false)),
     zoneKey_(keyType::null),
     needsUpdate_(true),
-    isoSurfCellPtr_(nullptr),
     isoSurfPtr_(nullptr),
     facesPtr_(nullptr)
-{
-//    dict.readIfPresent("zone", zoneKey_);
-//
-//    if (debug && zoneKey_.size() && mesh.cellZones().findZoneID(zoneKey_) < 0)
-//    {
-//        Info<< "cellZone " << zoneKey_
-//            << " not found - using entire mesh" << endl;
-//    }
-}
-
+{}
 
 
 Foam::distanceSurface::distanceSurface
@@ -384,12 +355,10 @@ Foam::distanceSurface::distanceSurface
     ),
     distance_(distance),
     signed_(signedDistance),
-    cell_(cell),
     regularise_(regularise),
     average_(average),
     zoneKey_(keyType::null),
     needsUpdate_(true),
-    isoSurfCellPtr_(nullptr),
     isoSurfPtr_(nullptr),
     facesPtr_(nullptr)
 {}
diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.H b/src/sampling/sampledSurface/distanceSurface/distanceSurface.H
index 1f2ce0832d..206cb72931 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurface.H
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.H
@@ -27,8 +27,6 @@ Class
 Description
     A sampledSurface defined by a distance to a surface.
 
-    Uses either isoSurfaceCell or isoSurface.
-
 SourceFiles
     distanceSurface.C
 
@@ -39,7 +37,6 @@ SourceFiles
 
 #include "sampledSurface.H"
 #include "searchableSurface.H"
-#include "isoSurfaceCell.H"
 #include "isoSurface.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -66,9 +63,6 @@ class distanceSurface
         //- Signed distance
         const bool signed_;
 
-        //- Whether to use isoSurfaceCell or isoSurface
-        const bool cell_;
-
         //- Whether to coarsen
         const Switch regularise_;
 
@@ -88,9 +82,6 @@ class distanceSurface
         //- Distance to points
         scalarField pointDistance_;
 
-        //- Constructed iso surface
-        autoPtr<isoSurfaceCell> isoSurfCellPtr_;
-
         //- Constructed iso surface
         autoPtr<isoSurface> isoSurfPtr_;
 
@@ -191,14 +182,7 @@ public:
 
         const triSurface& surface() const
         {
-            if (cell_)
-            {
-                return isoSurfCellPtr_();
-            }
-            else
-            {
-                return isoSurfPtr_();
-            }
+            return isoSurfPtr_();
         }
 
         //- Sample field on surface
diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C b/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C
index 6dcf8fa95d..6d2bdff1d0 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C
@@ -37,20 +37,10 @@ Foam::distanceSurface::sampleField
     const GeometricField<Type, fvPatchField, volMesh>& vField
 ) const
 {
-    if (cell_)
-    {
-        return tmp<Field<Type>>
-        (
-            new Field<Type>(vField, isoSurfCellPtr_().meshCells())
-        );
-    }
-    else
-    {
-        return tmp<Field<Type>>
-        (
-            new Field<Type>(vField, isoSurfPtr_().meshCells())
-        );
-    }
+    return tmp<Field<Type>>
+    (
+        new Field<Type>(vField, isoSurfPtr_().meshCells())
+    );
 }
 
 
@@ -72,31 +62,15 @@ Foam::distanceSurface::interpolateField
         volPointInterpolation::New(fvm).interpolate(volFld)
     );
 
-    // Sample.
-    if (cell_)
-    {
-        return isoSurfCellPtr_().interpolate
+    return isoSurfPtr_().interpolate
+    (
         (
-            (
-                average_
-              ? pointAverage(pointFld())()
-              : volFld
-            ),
-            pointFld()
-        );
-    }
-    else
-    {
-        return isoSurfPtr_().interpolate
-        (
-            (
-                average_
-              ? pointAverage(pointFld())()
-              : volFld
-            ),
-            pointFld()
-        );
-    }
+            average_
+          ? pointAverage(pointFld())()
+          : volFld
+        ),
+        pointFld()
+    );
 }
 
 
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.H b/src/sampling/sampledSurface/isoSurface/isoSurface.H
index 0225493678..4fce24834b 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.H
@@ -248,8 +248,6 @@ class isoSurface
         ) const;
 
 
-        //- Note: cannot use simpler isoSurfaceCell::generateTriPoints since
-        //  the need here to sometimes pass in remote 'snappedPoints'
         template<class Type>
         void generateTriPoints
         (
@@ -345,10 +343,6 @@ class isoSurface
 
 public:
 
-    //- Declare friendship with isoSurfaceCell to share some functionality
-    friend class isoSurfaceCell;
-
-
     //- Runtime type information
     TypeName("isoSurface");
 
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
deleted file mode 100644
index 748e19b6fb..0000000000
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
+++ /dev/null
@@ -1,1589 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "isoSurfaceCell.H"
-#include "dictionary.H"
-#include "polyMesh.H"
-#include "mergePoints.H"
-#include "tetMatcher.H"
-#include "syncTools.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(isoSurfaceCell, 0);
-}
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-Foam::scalar Foam::isoSurfaceCell::isoFraction
-(
-    const scalar s0,
-    const scalar s1
-) const
-{
-    scalar d = s1-s0;
-
-    if (mag(d) > vSmall)
-    {
-        return (iso_-s0)/d;
-    }
-    else
-    {
-        return -1.0;
-    }
-}
-
-
-bool Foam::isoSurfaceCell::isTriCut
-(
-    const triFace& tri,
-    const scalarField& pointValues
-) const
-{
-    bool aLower = (pointValues[tri[0]] < iso_);
-    bool bLower = (pointValues[tri[1]] < iso_);
-    bool cLower = (pointValues[tri[2]] < iso_);
-
-    return !(aLower == bLower && aLower == cLower);
-}
-
-
-Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
-(
-    const PackedBoolList& isTet,
-    const scalarField& cellValues,
-    const scalarField& pointValues,
-    const label celli
-) const
-{
-    const cell& cFaces = mesh_.cells()[celli];
-
-    if (isTet.get(celli) == 1)
-    {
-        forAll(cFaces, cFacei)
-        {
-            const face& f = mesh_.faces()[cFaces[cFacei]];
-
-            for (label fp = 1; fp < f.size() - 1; fp++)
-            {
-                triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
-
-                if (isTriCut(tri, pointValues))
-                {
-                    return CUT;
-                }
-            }
-        }
-        return NOTCUT;
-    }
-    else
-    {
-        bool cellLower = (cellValues[celli] < iso_);
-
-        // First check if there is any cut in cell
-        bool edgeCut = false;
-
-        forAll(cFaces, cFacei)
-        {
-            label facei = cFaces[cFacei];
-            const face& f = mesh_.faces()[facei];
-
-            // Check pyramids cut
-            forAll(f, fp)
-            {
-                if ((pointValues[f[fp]] < iso_) != cellLower)
-                {
-                    edgeCut = true;
-                    break;
-                }
-            }
-
-            if (edgeCut)
-            {
-                break;
-            }
-
-            const label fp0 = mesh_.tetBasePtIs()[facei];
-            label fp = f.fcIndex(fp0);
-            for (label i = 2; i < f.size(); i++)
-            {
-                label nextFp = f.fcIndex(fp);
-
-                if (isTriCut(triFace(f[fp0], f[fp], f[nextFp]), pointValues))
-                {
-                    edgeCut = true;
-                    break;
-                }
-
-                fp = nextFp;
-            }
-
-            if (edgeCut)
-            {
-                break;
-            }
-        }
-
-        if (edgeCut)
-        {
-            // Count actual cuts (expensive since addressing needed)
-            // Note: not needed if you don't want to preserve maxima/minima
-            // centred around cellcentre. In that case just always return CUT
-
-            const labelList& cPoints = mesh_.cellPoints(celli);
-
-            label nPyrCuts = 0;
-
-            forAll(cPoints, i)
-            {
-                if ((pointValues[cPoints[i]] < iso_) != cellLower)
-                {
-                    nPyrCuts++;
-                }
-            }
-
-            if (nPyrCuts == cPoints.size())
-            {
-                return SPHERE;
-            }
-            else
-            {
-                return CUT;
-            }
-        }
-        else
-        {
-            return NOTCUT;
-        }
-    }
-}
-
-
-void Foam::isoSurfaceCell::calcCutTypes
-(
-    const PackedBoolList& isTet,
-    const scalarField& cVals,
-    const scalarField& pVals
-)
-{
-    cellCutType_.setSize(mesh_.nCells());
-    nCutCells_ = 0;
-    forAll(mesh_.cells(), celli)
-    {
-        cellCutType_[celli] = calcCutType(isTet, cVals, pVals, celli);
-
-        if (cellCutType_[celli] == CUT)
-        {
-            nCutCells_++;
-        }
-    }
-
-    if (debug)
-    {
-        Pout<< "isoSurfaceCell : detected " << nCutCells_
-            << " candidate cut cells." << endl;
-    }
-}
-
-
-Foam::labelPair Foam::isoSurfaceCell::findCommonPoints
-(
-    const labelledTri& tri0,
-    const labelledTri& tri1
-)
-{
-    labelPair common(-1, -1);
-
-    label fp0 = 0;
-    label fp1 = findIndex(tri1, tri0[fp0]);
-
-    if (fp1 == -1)
-    {
-        fp0 = 1;
-        fp1 = findIndex(tri1, tri0[fp0]);
-    }
-
-    if (fp1 != -1)
-    {
-        // So tri0[fp0] is tri1[fp1]
-
-        // Find next common point
-        label fp0p1 = tri0.fcIndex(fp0);
-        label fp1p1 = tri1.fcIndex(fp1);
-        label fp1m1 = tri1.rcIndex(fp1);
-
-        if (tri0[fp0p1] == tri1[fp1p1] || tri0[fp0p1] == tri1[fp1m1])
-        {
-            common[0] = tri0[fp0];
-            common[1] = tri0[fp0p1];
-        }
-    }
-    return common;
-}
-
-
-Foam::point Foam::isoSurfaceCell::calcCentre(const triSurface& s)
-{
-    vector sum = Zero;
-
-    forAll(s, i)
-    {
-        sum += s[i].centre(s.points());
-    }
-    return sum/s.size();
-}
-
-
-Foam::pointIndexHit Foam::isoSurfaceCell::collapseSurface
-(
-    const label celli,
-    pointField& localPoints,
-    DynamicList<labelledTri, 64>& localTris
-) const
-{
-    pointIndexHit info(false, Zero, localTris.size());
-
-    if (localTris.size() == 1)
-    {
-        const labelledTri& tri = localTris[0];
-        info.setPoint(tri.centre(localPoints));
-        info.setHit();
-    }
-    else if (localTris.size() == 2)
-    {
-        // Check if the two triangles share an edge.
-        const labelledTri& tri0 = localTris[0];
-        const labelledTri& tri1 = localTris[1];
-
-        labelPair shared = findCommonPoints(tri0, tri1);
-
-        if (shared[0] != -1)
-        {
-            const vector n0 = tri0.normal(localPoints);
-            const vector n1 = tri1.normal(localPoints);
-
-            if ((n0 & n1) < 0)
-            {
-                // Too big an angle. Do not simplify.
-            }
-            else
-            {
-                info.setPoint
-                (
-                    0.5
-                  * (
-                        tri0.centre(localPoints)
-                      + tri1.centre(localPoints)
-                    )
-                );
-                info.setHit();
-            }
-        }
-    }
-    else if (localTris.size())
-    {
-        // Check if single region. Rare situation.
-        triSurface surf
-        (
-            localTris,
-            geometricSurfacePatchList(0),
-            localPoints,
-            true
-        );
-        localTris.clearStorage();
-
-        labelList faceZone;
-        label nZones = surf.markZones
-        (
-            boolList(surf.nEdges(), false),
-            faceZone
-        );
-
-        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)
-            {
-                info.setPoint(calcCentre(surf));
-                info.setHit();
-            }
-        }
-    }
-
-    return info;
-}
-
-
-void Foam::isoSurfaceCell::calcSnappedCc
-(
-    const PackedBoolList& isTet,
-    const scalarField& cVals,
-    const scalarField& pVals,
-
-    DynamicList<point>& snappedPoints,
-    labelList& snappedCc
-) const
-{
-    const pointField& cc = mesh_.cellCentres();
-    const pointField& pts = mesh_.points();
-
-    snappedCc.setSize(mesh_.nCells());
-    snappedCc = -1;
-
-    // Work arrays
-    DynamicList<point, 64> localPoints(64);
-    DynamicList<labelledTri, 64> localTris(64);
-    Map<label> pointToLocal(64);
-
-    forAll(mesh_.cells(), celli)
-    {
-        if (cellCutType_[celli] == CUT && isTet.get(celli) == 0)
-        {
-            scalar cVal = cVals[celli];
-
-            const cell& cFaces = mesh_.cells()[celli];
-
-            localPoints.clear();
-            localTris.clear();
-            pointToLocal.clear();
-
-            // Create points for all intersections close to cell centre
-            // (i.e. from pyramid edges)
-
-            forAll(cFaces, cFacei)
-            {
-                const face& f = mesh_.faces()[cFaces[cFacei]];
-
-                forAll(f, fp)
-                {
-                    label pointi = f[fp];
-
-                    scalar s = isoFraction(cVal, pVals[pointi]);
-
-                    if (s >= 0.0 && s <= 0.5)
-                    {
-                        if (pointToLocal.insert(pointi, localPoints.size()))
-                        {
-                            localPoints.append((1.0-s)*cc[celli]+s*pts[pointi]);
-                        }
-                    }
-                }
-            }
-
-            if (localPoints.size() == 1)
-            {
-                // No need for any analysis.
-                snappedCc[celli] = snappedPoints.size();
-                snappedPoints.append(localPoints[0]);
-
-                // Pout<< "cell:" << celli
-                //    << " at " << mesh_.cellCentres()[celli]
-                //    << " collapsing " << localPoints
-                //    << " intersections down to "
-                //    << snappedPoints[snappedCc[celli]] << endl;
-            }
-            else if (localPoints.size() == 2)
-            {
-                //? No need for any analysis.???
-                snappedCc[celli] = snappedPoints.size();
-                snappedPoints.append(0.5*(localPoints[0]+localPoints[1]));
-
-                // Pout<< "cell:" << celli
-                //    << " at " << mesh_.cellCentres()[celli]
-                //    << " collapsing " << localPoints
-                //    << " intersections down to "
-                //    << snappedPoints[snappedCc[celli]] << endl;
-            }
-            else if (localPoints.size())
-            {
-                // Need to analyse
-                forAll(cFaces, cFacei)
-                {
-                    label facei = cFaces[cFacei];
-                    const face& f = mesh_.faces()[facei];
-
-                    // Do a tetrahedralisation. Each face to cc becomes pyr.
-                    // Each pyr gets split into tets by diagonalisation
-                    // of face.
-
-                    const label fp0 = mesh_.tetBasePtIs()[facei];
-                    label fp = f.fcIndex(fp0);
-                    for (label i = 2; i < f.size(); i++)
-                    {
-                        label nextFp = f.fcIndex(fp);
-                        triFace tri(f[fp0], f[fp], f[nextFp]);
-
-                        // Get fractions for the three edges to cell centre
-                        FixedList<scalar, 3> s(3);
-                        s[0] = isoFraction(cVal, pVals[tri[0]]);
-                        s[1] = isoFraction(cVal, pVals[tri[1]]);
-                        s[2] = isoFraction(cVal, pVals[tri[2]]);
-
-                        if
-                        (
-                            (s[0] >= 0.0 && s[0] <= 0.5)
-                         && (s[1] >= 0.0 && s[1] <= 0.5)
-                         && (s[2] >= 0.0 && s[2] <= 0.5)
-                        )
-                        {
-                            if
-                            (
-                                (mesh_.faceOwner()[facei] == celli)
-                             == (cVal >= pVals[tri[0]])
-                            )
-                            {
-                                localTris.append
-                                (
-                                    labelledTri
-                                    (
-                                        pointToLocal[tri[1]],
-                                        pointToLocal[tri[0]],
-                                        pointToLocal[tri[2]],
-                                        0
-                                    )
-                                );
-                            }
-                            else
-                            {
-                                localTris.append
-                                (
-                                    labelledTri
-                                    (
-                                        pointToLocal[tri[0]],
-                                        pointToLocal[tri[1]],
-                                        pointToLocal[tri[2]],
-                                        0
-                                    )
-                                );
-                            }
-                        }
-
-                        fp = nextFp;
-                    }
-                }
-
-                pointField surfPoints;
-                surfPoints.transfer(localPoints);
-                pointIndexHit info = collapseSurface
-                (
-                    celli,
-                    surfPoints,
-                    localTris
-                );
-
-                if (info.hit())
-                {
-                    snappedCc[celli] = snappedPoints.size();
-                    snappedPoints.append(info.hitPoint());
-
-                    // Pout<< "cell:" << celli
-                    //    << " at " << mesh_.cellCentres()[celli]
-                    //    << " collapsing " << surfPoints
-                    //    << " intersections down to "
-                    //    << snappedPoints[snappedCc[celli]] << endl;
-                }
-            }
-        }
-    }
-}
-
-
-void Foam::isoSurfaceCell::genPointTris
-(
-    const scalarField& cellValues,
-    const scalarField& pointValues,
-    const label pointi,
-    const label facei,
-    const label celli,
-    DynamicList<point, 64>& localTriPoints
-) const
-{
-    const pointField& cc = mesh_.cellCentres();
-    const pointField& pts = mesh_.points();
-    const face& f = mesh_.faces()[facei];
-
-    const label fp0 = mesh_.tetBasePtIs()[facei];
-    label fp = f.fcIndex(fp0);
-    for (label i = 2; i < f.size(); i++)
-    {
-        label nextFp = f.fcIndex(fp);
-        triFace tri(f[fp0], f[fp], f[nextFp]);
-
-        label index = findIndex(tri, pointi);
-
-        if (index == -1)
-        {
-            continue;
-        }
-
-        // Tet between index..index-1, index..index+1, index..cc
-        label b = tri[tri.fcIndex(index)];
-        label c = tri[tri.rcIndex(index)];
-
-        // Get fractions for the three edges emanating from point
-        FixedList<scalar, 3> s(3);
-        s[0] = isoFraction(pointValues[pointi], pointValues[b]);
-        s[1] = isoFraction(pointValues[pointi], pointValues[c]);
-        s[2] = isoFraction(pointValues[pointi], cellValues[celli]);
-
-        if
-        (
-            (s[0] >= 0.0 && s[0] <= 0.5)
-         && (s[1] >= 0.0 && s[1] <= 0.5)
-         && (s[2] >= 0.0 && s[2] <= 0.5)
-        )
-        {
-            point p0 = (1.0-s[0])*pts[pointi] + s[0]*pts[b];
-            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)
-             == (pointValues[pointi] > cellValues[celli])
-            )
-            {
-                localTriPoints.append(p0);
-                localTriPoints.append(p1);
-                localTriPoints.append(p2);
-            }
-            else
-            {
-                localTriPoints.append(p1);
-                localTriPoints.append(p0);
-                localTriPoints.append(p2);
-            }
-        }
-
-        fp = nextFp;
-    }
-}
-
-
-void Foam::isoSurfaceCell::genPointTris
-(
-    const scalarField& pointValues,
-    const label pointi,
-    const label facei,
-    const label celli,
-    DynamicList<point, 64>& localTriPoints
-) const
-{
-    const pointField& pts = mesh_.points();
-    const cell& cFaces = mesh_.cells()[celli];
-
-    FixedList<label, 4> tet;
-
-    // Make tet from this face to the 4th point (same as cellcentre in
-    // non-tet cells)
-    const face& f = mesh_.faces()[facei];
-
-    // Find 4th point
-    label ccPointi = -1;
-    forAll(cFaces, cFacei)
-    {
-        const face& f1 = mesh_.faces()[cFaces[cFacei]];
-        forAll(f1, fp)
-        {
-            label p1 = f1[fp];
-
-            if (findIndex(f, p1) == -1)
-            {
-                ccPointi = p1;
-                break;
-            }
-        }
-        if (ccPointi != -1)
-        {
-            break;
-        }
-    }
-
-
-    // Tet between index..index-1, index..index+1, index..cc
-    label index = findIndex(f, pointi);
-    label b = f[f.fcIndex(index)];
-    label c = f[f.rcIndex(index)];
-
-    // Pout<< " p0:" << pointi << " b:" << b << " c:" << c
-    //<< " d:" << ccPointi << endl;
-
-    // Get fractions for the three edges emanating from point
-    FixedList<scalar, 3> s(3);
-    s[0] = isoFraction(pointValues[pointi], pointValues[b]);
-    s[1] = isoFraction(pointValues[pointi], pointValues[c]);
-    s[2] = isoFraction(pointValues[pointi], pointValues[ccPointi]);
-
-    if
-    (
-        (s[0] >= 0.0 && s[0] <= 0.5)
-     && (s[1] >= 0.0 && s[1] <= 0.5)
-     && (s[2] >= 0.0 && s[2] <= 0.5)
-    )
-    {
-        point p0 = (1.0-s[0])*pts[pointi] + s[0]*pts[b];
-        point p1 = (1.0-s[1])*pts[pointi] + s[1]*pts[c];
-        point p2 = (1.0-s[2])*pts[pointi] + s[2]*pts[ccPointi];
-
-        if (mesh_.faceOwner()[facei] != celli)
-        {
-            localTriPoints.append(p0);
-            localTriPoints.append(p1);
-            localTriPoints.append(p2);
-        }
-        else
-        {
-            localTriPoints.append(p1);
-            localTriPoints.append(p0);
-            localTriPoints.append(p2);
-        }
-    }
-}
-
-
-void Foam::isoSurfaceCell::calcSnappedPoint
-(
-    const PackedBoolList& isTet,
-    const scalarField& cVals,
-    const scalarField& pVals,
-
-    DynamicList<point>& snappedPoints,
-    labelList& snappedPoint
-) const
-{
-    // 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());
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-    forAll(patches, patchi)
-    {
-        const polyPatch& pp = patches[patchi];
-
-        if (!pp.coupled())
-        {
-            label facei = pp.start();
-            forAll(pp, i)
-            {
-                const face& f = mesh_.faces()[facei++];
-
-                forAll(f, fp)
-                {
-                    isBoundaryPoint.set(f[fp], 1);
-                }
-            }
-        }
-    }
-
-
-    const point greatPoint(great, great, great);
-
-    pointField collapsedPoint(mesh_.nPoints(), greatPoint);
-
-
-    // Work arrays
-    DynamicList<point, 64> localTriPoints(100);
-    labelHashSet localPointCells(100);
-
-    forAll(mesh_.pointFaces(), pointi)
-    {
-        if (isBoundaryPoint.get(pointi) == 1)
-        {
-            continue;
-        }
-
-        const labelList& pFaces = mesh_.pointFaces()[pointi];
-
-        bool anyCut = false;
-
-        forAll(pFaces, i)
-        {
-            label facei = pFaces[i];
-
-            if
-            (
-                cellCutType_[mesh_.faceOwner()[facei]] == CUT
-             || (
-                    mesh_.isInternalFace(facei)
-                 && cellCutType_[mesh_.faceNeighbour()[facei]] == CUT
-                )
-            )
-            {
-                anyCut = true;
-                break;
-            }
-        }
-
-        if (!anyCut)
-        {
-            continue;
-        }
-
-
-        // Do a pointCells walk (by using pointFaces)
-
-        localPointCells.clear();
-        localTriPoints.clear();
-
-        forAll(pFaces, pFacei)
-        {
-            label facei = pFaces[pFacei];
-            label own = mesh_.faceOwner()[facei];
-
-            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);
-                }
-            }
-            else
-            {
-                genPointTris
-                (
-                    cVals,
-                    pVals,
-                    pointi,
-                    facei,
-                    own,
-                    localTriPoints
-                );
-            }
-
-            if (mesh_.isInternalFace(facei))
-            {
-                label nei = mesh_.faceNeighbour()[facei];
-
-                if (isTet.get(nei) == 1)
-                {
-                    if (localPointCells.insert(nei))
-                    {
-                        genPointTris(pVals, pointi, facei, nei, localTriPoints);
-                    }
-                }
-                else
-                {
-                    genPointTris
-                    (
-                        cVals,
-                        pVals,
-                        pointi,
-                        facei,
-                        nei,
-                        localTriPoints
-                    );
-                }
-            }
-        }
-
-        if (localTriPoints.size() == 3)
-        {
-            // Single triangle. No need for any analysis. Average points.
-            pointField points;
-            points.transfer(localTriPoints);
-            collapsedPoint[pointi] = sum(points)/points.size();
-
-            // Pout<< "    point:" << pointi
-            //    << " replacing coord:" << mesh_.points()[pointi]
-            //    << " by average:" << collapsedPoint[pointi] << endl;
-        }
-        else if (localTriPoints.size())
-        {
-            // Convert points into triSurface.
-
-            // Merge points and compact out non-valid triangles
-            labelList triMap;               // merged to unmerged triangle
-            labelList triPointReverseMap;   // unmerged to merged point
-            triSurface surf
-            (
-                stitchTriPoints
-                (
-                    false,                  // do not check for duplicate tris
-                    localTriPoints,
-                    triPointReverseMap,
-                    triMap
-                )
-            );
-
-            labelList faceZone;
-            label nZones = surf.markZones
-            (
-                boolList(surf.nEdges(), false),
-                faceZone
-            );
-
-            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)
-                {
-                    collapsedPoint[pointi] = calcCentre(surf);
-                }
-            }
-        }
-    }
-
-    syncTools::syncPointPositions
-    (
-        mesh_,
-        collapsedPoint,
-        minMagSqrEqOp<point>(),
-        greatPoint
-    );
-
-    snappedPoint.setSize(mesh_.nPoints());
-    snappedPoint = -1;
-
-    forAll(collapsedPoint, pointi)
-    {
-        // 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]);
-        }
-    }
-}
-
-
-Foam::triSurface Foam::isoSurfaceCell::stitchTriPoints
-(
-    const bool checkDuplicates,
-    const List<point>& triPoints,
-    labelList& triPointReverseMap,  // unmerged to merged point
-    labelList& triMap               // merged to unmerged triangle
-) const
-{
-    label nTris = triPoints.size()/3;
-
-    if ((triPoints.size() % 3) != 0)
-    {
-        FatalErrorInFunction
-            << "Problem: number of points " << triPoints.size()
-            << " not a multiple of 3." << abort(FatalError);
-    }
-
-    pointField newPoints;
-    mergePoints
-    (
-        triPoints,
-        mergeDistance_,
-        false,
-        triPointReverseMap,
-        newPoints
-    );
-
-    // Check that enough merged.
-    if (debug)
-    {
-        Pout<< "isoSurfaceCell : merged from " << triPoints.size()
-            << " points down to " << newPoints.size() << endl;
-
-        pointField newNewPoints;
-        labelList oldToNew;
-        bool hasMerged = mergePoints
-        (
-            newPoints,
-            mergeDistance_,
-            true,
-            oldToNew,
-            newNewPoints
-        );
-
-        if (hasMerged)
-        {
-            FatalErrorInFunction
-                << "Merged points contain duplicates"
-                << " when merging with distance " << mergeDistance_ << endl
-                << "merged:" << newPoints.size() << " re-merged:"
-                << newNewPoints.size()
-                << abort(FatalError);
-        }
-    }
-
-
-    List<labelledTri> tris;
-    {
-        DynamicList<labelledTri> dynTris(nTris);
-        label rawPointi = 0;
-        DynamicList<label> newToOldTri(nTris);
-
-        for (label oldTriI = 0; oldTriI < nTris; oldTriI++)
-        {
-            labelledTri tri
-            (
-                triPointReverseMap[rawPointi],
-                triPointReverseMap[rawPointi+1],
-                triPointReverseMap[rawPointi+2],
-                0
-            );
-            if ((tri[0] != tri[1]) && (tri[0] != tri[2]) && (tri[1] != tri[2]))
-            {
-                newToOldTri.append(oldTriI);
-                dynTris.append(tri);
-            }
-
-            rawPointi += 3;
-        }
-
-        triMap.transfer(newToOldTri);
-        tris.transfer(dynTris);
-    }
-
-
-    // Use face centres to determine 'flat hole' situation (see RMT paper).
-    // Two unconnected triangles get connected because (some of) the edges
-    // separating them get collapsed. Below only checks for duplicate triangles,
-    // not non-manifold edge connectivity.
-    if (checkDuplicates)
-    {
-        if (debug)
-        {
-            Pout<< "isoSurfaceCell : merged from " << nTris
-                << " down to " << tris.size() << " triangles." << endl;
-        }
-
-        pointField centres(tris.size());
-        forAll(tris, triI)
-        {
-            centres[triI] = tris[triI].centre(newPoints);
-        }
-
-        pointField mergedCentres;
-        labelList oldToMerged;
-        bool hasMerged = mergePoints
-        (
-            centres,
-            mergeDistance_,
-            false,
-            oldToMerged,
-            mergedCentres
-        );
-
-        if (debug)
-        {
-            Pout<< "isoSurfaceCell : detected "
-                << centres.size()-mergedCentres.size()
-                << " duplicate triangles." << endl;
-        }
-
-        if (hasMerged)
-        {
-            // Filter out duplicates.
-            label newTriI = 0;
-            DynamicList<label> newToOldTri(tris.size());
-            labelList newToMaster(mergedCentres.size(), -1);
-            forAll(tris, triI)
-            {
-                label mergedI = oldToMerged[triI];
-
-                if (newToMaster[mergedI] == -1)
-                {
-                    newToMaster[mergedI] = triI;
-                    newToOldTri.append(triMap[triI]);
-                    tris[newTriI++] = tris[triI];
-                }
-            }
-
-            triMap.transfer(newToOldTri);
-            tris.setSize(newTriI);
-        }
-    }
-
-    return triSurface(tris, geometricSurfacePatchList(0), newPoints, true);
-}
-
-
-bool Foam::isoSurfaceCell::validTri(const triSurface& surf, const label facei)
-{
-    // Simple check on indices ok.
-
-    const labelledTri& f = surf[facei];
-
-    forAll(f, fp)
-    {
-        if (f[fp] < 0 || f[fp] >= surf.points().size())
-        {
-            WarningInFunction
-                << "triangle " << facei << " vertices " << f
-                << " uses point indices outside point range 0.."
-                << surf.points().size()-1 << endl;
-
-            return false;
-        }
-    }
-
-    if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
-    {
-        WarningInFunction
-            << "triangle " << facei
-            << " uses non-unique vertices " << f
-            << endl;
-        return false;
-    }
-
-    // duplicate triangle check
-
-    const labelList& fFaces = surf.faceFaces()[facei];
-
-    // Check if faceNeighbours use same points as this face.
-    // Note: discards normal information - sides of baffle are merged.
-    forAll(fFaces, i)
-    {
-        label nbrFacei = fFaces[i];
-
-        if (nbrFacei <= facei)
-        {
-            // lower numbered faces already checked
-            continue;
-        }
-
-        const labelledTri& nbrF = surf[nbrFacei];
-
-        if
-        (
-            ((f[0] == nbrF[0]) || (f[0] == nbrF[1]) || (f[0] == nbrF[2]))
-         && ((f[1] == nbrF[0]) || (f[1] == nbrF[1]) || (f[1] == nbrF[2]))
-         && ((f[2] == nbrF[0]) || (f[2] == nbrF[1]) || (f[2] == nbrF[2]))
-        )
-        {
-            WarningInFunction
-                << "triangle " << facei << " vertices " << f
-                << " coords:" << f.points(surf.points())
-                << " has the same vertices as triangle " << nbrFacei
-                << " vertices " << nbrF
-                << endl;
-
-            return false;
-        }
-    }
-    return true;
-}
-
-
-void Foam::isoSurfaceCell::calcAddressing
-(
-    const triSurface& surf,
-    List<FixedList<label, 3>>& faceEdges,
-    labelList& edgeFace0,
-    labelList& edgeFace1,
-    Map<labelList>& edgeFacesRest
-) const
-{
-    const pointField& points = surf.points();
-
-    pointField edgeCentres(3*surf.size());
-    label edgeI = 0;
-    forAll(surf, triI)
-    {
-        const labelledTri& tri = surf[triI];
-        edgeCentres[edgeI++] = 0.5*(points[tri[0]]+points[tri[1]]);
-        edgeCentres[edgeI++] = 0.5*(points[tri[1]]+points[tri[2]]);
-        edgeCentres[edgeI++] = 0.5*(points[tri[2]]+points[tri[0]]);
-    }
-
-    pointField mergedCentres;
-    labelList oldToMerged;
-    bool hasMerged = mergePoints
-    (
-        edgeCentres,
-        mergeDistance_,
-        false,
-        oldToMerged,
-        mergedCentres
-    );
-
-    if (debug)
-    {
-        Pout<< "isoSurfaceCell : detected "
-            << mergedCentres.size()
-            << " edges on " << surf.size() << " triangles." << endl;
-    }
-
-    if (!hasMerged)
-    {
-        return;
-    }
-
-
-    // Determine faceEdges
-    faceEdges.setSize(surf.size());
-    edgeI = 0;
-    forAll(surf, triI)
-    {
-        faceEdges[triI][0] = oldToMerged[edgeI++];
-        faceEdges[triI][1] = oldToMerged[edgeI++];
-        faceEdges[triI][2] = oldToMerged[edgeI++];
-    }
-
-
-    // Determine edgeFaces
-    edgeFace0.setSize(mergedCentres.size());
-    edgeFace0 = -1;
-    edgeFace1.setSize(mergedCentres.size());
-    edgeFace1 = -1;
-    edgeFacesRest.clear();
-
-    forAll(oldToMerged, oldEdgeI)
-    {
-        label triI = oldEdgeI / 3;
-        label edgeI = oldToMerged[oldEdgeI];
-
-        if (edgeFace0[edgeI] == -1)
-        {
-            edgeFace0[edgeI] = triI;
-        }
-        else if (edgeFace1[edgeI] == -1)
-        {
-            edgeFace1[edgeI] = triI;
-        }
-        else
-        {
-            // WarningInFunction
-            //    << "Edge " << edgeI << " with centre " << mergedCentres[edgeI]
-            //    << " used by more than two triangles: " << edgeFace0[edgeI]
-            //    << ", "
-            //    << edgeFace1[edgeI] << " and " << triI << endl;
-            Map<labelList>::iterator iter = edgeFacesRest.find(edgeI);
-
-            if (iter != edgeFacesRest.end())
-            {
-                labelList& eFaces = iter();
-                label sz = eFaces.size();
-                eFaces.setSize(sz+1);
-                eFaces[sz] = triI;
-            }
-            else
-            {
-                edgeFacesRest.insert(edgeI, labelList(1, triI));
-            }
-        }
-    }
-}
-
-
-bool Foam::isoSurfaceCell::danglingTriangle
-(
-    const FixedList<label, 3>& fEdges,
-    const labelList& edgeFace1
-)
-{
-    label nOpen = 0;
-    forAll(fEdges, i)
-    {
-        if (edgeFace1[fEdges[i]] == -1)
-        {
-            nOpen++;
-        }
-    }
-
-    if (nOpen == 1 || nOpen == 2 || nOpen == 3)
-    {
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
-Foam::label Foam::isoSurfaceCell::markDanglingTriangles
-(
-    const List<FixedList<label, 3>>& faceEdges,
-    const labelList& edgeFace0,
-    const labelList& edgeFace1,
-    const Map<labelList>& edgeFacesRest,
-    boolList& keepTriangles
-)
-{
-    keepTriangles.setSize(faceEdges.size());
-    keepTriangles = true;
-
-    label nDangling = 0;
-
-    // Remove any dangling triangles
-    forAllConstIter(Map<labelList>,  edgeFacesRest, iter)
-    {
-        // These are all the non-manifold edges. Filter out all triangles
-        // with only one connected edge (= this edge)
-
-        label edgeI = iter.key();
-        const labelList& otherEdgeFaces = iter();
-
-        // Remove all dangling triangles
-        if (danglingTriangle(faceEdges[edgeFace0[edgeI]], edgeFace1))
-        {
-            keepTriangles[edgeFace0[edgeI]] = false;
-            nDangling++;
-        }
-        if (danglingTriangle(faceEdges[edgeFace1[edgeI]], edgeFace1))
-        {
-            keepTriangles[edgeFace1[edgeI]] = false;
-            nDangling++;
-        }
-        forAll(otherEdgeFaces, i)
-        {
-            label triI = otherEdgeFaces[i];
-            if (danglingTriangle(faceEdges[triI], edgeFace1))
-            {
-                keepTriangles[triI] = false;
-                nDangling++;
-            }
-        }
-    }
-    return nDangling;
-}
-
-
-Foam::triSurface Foam::isoSurfaceCell::subsetMesh
-(
-    const triSurface& s,
-    const labelList& newToOldFaces,
-    labelList& oldToNewPoints,
-    labelList& newToOldPoints
-)
-{
-    const boolList include
-    (
-        createWithValues<boolList>
-        (
-            s.size(),
-            false,
-            newToOldFaces,
-            true
-        )
-    );
-
-    newToOldPoints.setSize(s.points().size());
-    oldToNewPoints.setSize(s.points().size());
-    oldToNewPoints = -1;
-    {
-        label pointi = 0;
-
-        forAll(include, oldFacei)
-        {
-            if (include[oldFacei])
-            {
-                // Renumber labels for face
-                const triSurface::FaceType& f = s[oldFacei];
-
-                forAll(f, fp)
-                {
-                    label oldPointi = f[fp];
-
-                    if (oldToNewPoints[oldPointi] == -1)
-                    {
-                        oldToNewPoints[oldPointi] = pointi;
-                        newToOldPoints[pointi++] = oldPointi;
-                    }
-                }
-            }
-        }
-        newToOldPoints.setSize(pointi);
-    }
-
-    // Extract points
-    pointField newPoints(newToOldPoints.size());
-    forAll(newToOldPoints, i)
-    {
-        newPoints[i] = s.points()[newToOldPoints[i]];
-    }
-    // Extract faces
-    List<labelledTri> newTriangles(newToOldFaces.size());
-
-    forAll(newToOldFaces, i)
-    {
-        // Get old vertex labels
-        const labelledTri& tri = s[newToOldFaces[i]];
-
-        newTriangles[i][0] = oldToNewPoints[tri[0]];
-        newTriangles[i][1] = oldToNewPoints[tri[1]];
-        newTriangles[i][2] = oldToNewPoints[tri[2]];
-        newTriangles[i].region() = tri.region();
-    }
-
-    // Reuse storage.
-    return triSurface(newTriangles, s.patches(), newPoints, true);
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::isoSurfaceCell::isoSurfaceCell
-(
-    const polyMesh& mesh,
-    const scalarField& cVals,
-    const scalarField& pVals,
-    const scalar iso,
-    const bool regularise,
-    const scalar mergeTol
-)
-:
-    mesh_(mesh),
-    cVals_(cVals),
-    pVals_(pVals),
-    iso_(iso),
-    mergeDistance_(mergeTol*mesh.bounds().mag())
-{
-    if (debug)
-    {
-        Pout<< "isoSurfaceCell : mergeTol:" << mergeTol
-            << " mesh span:" << mesh.bounds().mag()
-            << " mergeDistance:" << mergeDistance_ << endl;
-    }
-
-    // Determine if cell is tet
-    PackedBoolList isTet(mesh_.nCells());
-    {
-        tetMatcher tet;
-
-        forAll(isTet, celli)
-        {
-            if (tet.isA(mesh_, celli))
-            {
-                isTet.set(celli, 1);
-            }
-        }
-    }
-
-
-    // Determine if any cut through cell
-    calcCutTypes(isTet, cVals, pVals);
-
-    DynamicList<point> snappedPoints(nCutCells_);
-
-    // Per cc -1 or a point inside snappedPoints.
-    labelList snappedCc;
-    if (regularise)
-    {
-        calcSnappedCc
-        (
-            isTet,
-            cVals,
-            pVals,
-            snappedPoints,
-            snappedCc
-        );
-    }
-    else
-    {
-        snappedCc.setSize(mesh_.nCells());
-        snappedCc = -1;
-    }
-
-    if (debug)
-    {
-        Pout<< "isoSurfaceCell : shifted " << snappedPoints.size()
-            << " cell centres to intersection." << endl;
-    }
-
-    snappedPoints.shrink();
-    label nCellSnaps = snappedPoints.size();
-
-    // Per point -1 or a point inside snappedPoints.
-    labelList snappedPoint;
-    if (regularise)
-    {
-        calcSnappedPoint
-        (
-            isTet,
-            cVals,
-            pVals,
-            snappedPoints,
-            snappedPoint
-        );
-    }
-    else
-    {
-        snappedPoint.setSize(mesh_.nPoints());
-        snappedPoint = -1;
-    }
-
-    if (debug)
-    {
-        Pout<< "isoSurfaceCell : shifted " << snappedPoints.size()-nCellSnaps
-            << " vertices to intersection." << endl;
-    }
-
-
-    {
-        DynamicList<point> triPoints(nCutCells_);
-        DynamicList<label> triMeshCells(nCutCells_);
-
-        generateTriPoints
-        (
-            cVals,
-            pVals,
-
-            mesh_.cellCentres(),
-            mesh_.points(),
-
-            snappedPoints,
-            snappedCc,
-            snappedPoint,
-
-            triPoints,
-            triMeshCells
-        );
-
-        if (debug)
-        {
-            Pout<< "isoSurfaceCell : generated " << triMeshCells.size()
-                << " unmerged triangles." << endl;
-        }
-
-        // Merge points and compact out non-valid triangles
-        labelList triMap;
-        triSurface::operator=
-        (
-            stitchTriPoints
-            (
-                regularise,         // check for duplicate tris
-                triPoints,
-                triPointMergeMap_,  // unmerged to merged point
-                triMap              // merged to unmerged triangle
-            )
-        );
-
-        if (debug)
-        {
-            Pout<< "isoSurfaceCell : generated " << triMap.size()
-                << " merged triangles." << endl;
-        }
-
-        meshCells_.setSize(triMap.size());
-        forAll(triMap, i)
-        {
-            meshCells_[i] = triMeshCells[triMap[i]];
-        }
-    }
-
-
-    if (debug)
-    {
-        Pout<< "isoSurfaceCell : checking " << size()
-            << " triangles for validity." << endl;
-
-        forAll(*this, triI)
-        {
-            // Copied from surfaceCheck
-            validTri(*this, triI);
-        }
-    }
-
-
-    if (regularise)
-    {
-        List<FixedList<label, 3>> faceEdges;
-        labelList edgeFace0, edgeFace1;
-        Map<labelList> edgeFacesRest;
-
-
-        while (true)
-        {
-            // Calculate addressing
-            calcAddressing
-            (
-                *this,
-                faceEdges,
-                edgeFace0,
-                edgeFace1,
-                edgeFacesRest
-            );
-
-            // See if any dangling triangles
-            boolList keepTriangles;
-            label nDangling = markDanglingTriangles
-            (
-                faceEdges,
-                edgeFace0,
-                edgeFace1,
-                edgeFacesRest,
-                keepTriangles
-            );
-
-            if (debug)
-            {
-                Pout<< "isoSurfaceCell : detected " << nDangling
-                    << " dangling triangles." << endl;
-            }
-
-            if (nDangling == 0)
-            {
-                break;
-            }
-
-            // Create face map (new to old)
-            labelList subsetTriMap(findIndices(keepTriangles, true));
-
-            labelList subsetPointMap;
-            labelList reversePointMap;
-            triSurface::operator=
-            (
-                subsetMesh
-                (
-                    *this,
-                    subsetTriMap,
-                    reversePointMap,
-                    subsetPointMap
-                )
-            );
-            meshCells_ = labelField(meshCells_, subsetTriMap);
-            inplaceRenumber(reversePointMap, triPointMergeMap_);
-        }
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
deleted file mode 100644
index 8dbda9f0ae..0000000000
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
+++ /dev/null
@@ -1,352 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::isoSurfaceCell
-
-Description
-    A surface formed by the iso value.
-    After "Polygonising A Scalar Field Using Tetrahedrons", Paul Bourke
-    (http://paulbourke.net/geometry/polygonise) and
-    "Regularised Marching Tetrahedra: improved iso-surface extraction",
-    G.M. Treece, R.W. Prager and A.H. Gee.
-
-    See isoSurface. This is a variant which does tetrahedrisation from
-    triangulation of face to cell centre instead of edge of face to two
-    neighbouring cell centres. This gives much lower quality triangles
-    but they are local to a cell.
-
-SourceFiles
-    isoSurfaceCell.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef isoSurfaceCell_H
-#define isoSurfaceCell_H
-
-#include "triSurface.H"
-#include "labelPair.H"
-#include "pointIndexHit.H"
-#include "PackedBoolList.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-class polyMesh;
-
-/*---------------------------------------------------------------------------*\
-                       Class isoSurfaceCell Declaration
-\*---------------------------------------------------------------------------*/
-
-class isoSurfaceCell
-:
-    public triSurface
-{
-    // Private data
-
-        enum segmentCutType
-        {
-            NEARFIRST,      // intersection close to e.first()
-            NEARSECOND,     //  ,,                   e.second()
-            NOTNEAR         // no intersection
-        };
-
-        enum cellCutType
-        {
-            NOTCUT,     // not cut
-            SPHERE,     // all edges to cell centre cut
-            CUT         // normal cut
-        };
-
-
-        //- Reference to mesh
-        const polyMesh& mesh_;
-
-        const scalarField& cVals_;
-
-        const scalarField& pVals_;
-
-        //- isoSurfaceCell value
-        const scalar iso_;
-
-        //- When to merge points
-        const scalar mergeDistance_;
-
-        //- Whether cell might be cut
-        List<cellCutType> cellCutType_;
-
-        //- Estimated number of cut cells
-        label nCutCells_;
-
-        //- For every triangle the original cell in mesh
-        labelList meshCells_;
-
-        //- For every unmerged triangle point the point in the triSurface
-        labelList triPointMergeMap_;
-
-
-    // Private Member Functions
-
-        //- Get location of iso value as fraction in between s0,s1
-        scalar isoFraction
-        (
-            const scalar s0,
-            const scalar s1
-        ) const;
-
-        //- Does any edge of triangle cross iso value?
-        bool isTriCut
-        (
-            const triFace& tri,
-            const scalarField& pointValues
-        ) const;
-
-        //- Determine whether cell is cut
-        cellCutType calcCutType
-        (
-            const PackedBoolList&,
-            const scalarField& cellValues,
-            const scalarField& pointValues,
-            const label
-        ) const;
-
-        void calcCutTypes
-        (
-            const PackedBoolList&,
-            const scalarField& cellValues,
-            const scalarField& pointValues
-        );
-
-        //- Return the two common points between two triangles
-        static labelPair findCommonPoints
-        (
-            const labelledTri&,
-            const labelledTri&
-        );
-
-        //- Calculate centre of surface.
-        static point calcCentre(const triSurface&);
-
-        //- Replace surface (localPoints, localTris) with single point.
-        //  Returns point.  Destroys arguments.
-        pointIndexHit collapseSurface
-        (
-            const label celli,
-            pointField& localPoints,
-            DynamicList<labelledTri, 64>& localTris
-        ) const;
-
-        //- Determine per cc whether all near cuts can be snapped to single
-        //  point.
-        void calcSnappedCc
-        (
-            const PackedBoolList& isTet,
-            const scalarField& cVals,
-            const scalarField& pVals,
-            DynamicList<point>& triPoints,
-            labelList& snappedCc
-        ) const;
-
-        //- Generate triangles for face connected to pointi
-        void genPointTris
-        (
-            const scalarField& cellValues,
-            const scalarField& pointValues,
-            const label pointi,
-            const label facei,
-            const label celli,
-            DynamicList<point, 64>& localTriPoints
-        ) const;
-
-        //- Generate triangles for tet connected to pointi
-        void genPointTris
-        (
-            const scalarField& pointValues,
-            const label pointi,
-            const label facei,
-            const label celli,
-            DynamicList<point, 64>& localTriPoints
-        ) const;
-
-        //- Determine per point whether all near cuts can be snapped to single
-        //  point.
-        void calcSnappedPoint
-        (
-            const PackedBoolList& isTet,
-            const scalarField& cVals,
-            const scalarField& pVals,
-            DynamicList<point>& triPoints,
-            labelList& snappedPoint
-        ) const;
-
-        //- Generate single point by interpolation or snapping
-        template<class Type>
-        Type generatePoint
-        (
-            const DynamicList<Type>& snappedPoints,
-            const scalar s0,
-            const Type& p0,
-            const label p0Index,
-            const scalar s1,
-            const Type& p1,
-            const label p1Index
-        ) const;
-
-        template<class Type>
-        void generateTriPoints
-        (
-            const DynamicList<Type>& snapped,
-            const scalar s0,
-            const Type& p0,
-            const label p0Index,
-            const scalar s1,
-            const Type& p1,
-            const label p1Index,
-            const scalar s2,
-            const Type& p2,
-            const label p2Index,
-            const scalar s3,
-            const Type& p3,
-            const label p3Index,
-            DynamicList<Type>& points
-        ) const;
-
-        template<class Type>
-        void generateTriPoints
-        (
-            const scalarField& cVals,
-            const scalarField& pVals,
-
-            const Field<Type>& cCoords,
-            const Field<Type>& pCoords,
-
-            const DynamicList<Type>& snappedPoints,
-            const labelList& snappedCc,
-            const labelList& snappedPoint,
-
-            DynamicList<Type>& triPoints,
-            DynamicList<label>& triMeshCells
-        ) const;
-
-        triSurface stitchTriPoints
-        (
-            const bool checkDuplicates,
-            const List<point>& triPoints,
-            labelList& triPointReverseMap,  // unmerged to merged point
-            labelList& triMap               // merged to unmerged triangle
-        ) const;
-
-        //- Check single triangle for (topological) validity
-        static bool validTri(const triSurface&, const label);
-
-        //- Determine edge-face addressing
-        void calcAddressing
-        (
-            const triSurface& surf,
-            List<FixedList<label, 3>>& faceEdges,
-            labelList& edgeFace0,
-            labelList& edgeFace1,
-            Map<labelList>& edgeFacesRest
-        ) const;
-
-        //- Is triangle (given by 3 edges) not fully connected?
-        static bool danglingTriangle
-        (
-            const FixedList<label, 3>& fEdges,
-            const labelList& edgeFace1
-        );
-
-        //- Mark all non-fully connected triangles
-        static label markDanglingTriangles
-        (
-            const List<FixedList<label, 3>>& faceEdges,
-            const labelList& edgeFace0,
-            const labelList& edgeFace1,
-            const Map<labelList>& edgeFacesRest,
-            boolList& keepTriangles
-        );
-
-        static triSurface subsetMesh
-        (
-            const triSurface& s,
-            const labelList& newToOldFaces,
-            labelList& oldToNewPoints,
-            labelList& newToOldPoints
-        );
-
-
-public:
-
-    //- Runtime type information
-    TypeName("isoSurfaceCell");
-
-
-    // Constructors
-
-        //- Construct from dictionary
-        isoSurfaceCell
-        (
-            const polyMesh& mesh,
-            const scalarField& cellValues,
-            const scalarField& pointValues,
-            const scalar iso,
-            const bool regularise,
-            const scalar mergeTol = 1e-6    // fraction of bounding box
-        );
-
-
-    // Member Functions
-
-        //- For every face original cell in mesh
-        const labelList& meshCells() const
-        {
-            return meshCells_;
-        }
-
-        //- Interpolates cCoords,pCoords.
-        template<class Type>
-        tmp<Field<Type>> interpolate
-        (
-            const Field<Type>& cCoords,
-            const Field<Type>& pCoords
-        ) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-    #include "isoSurfaceCellTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
deleted file mode 100644
index 7c486546cf..0000000000
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
+++ /dev/null
@@ -1,566 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "isoSurfaceCell.H"
-#include "polyMesh.H"
-#include "tetMatcher.H"
-#include "isoSurface.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-template<class Type>
-Type Foam::isoSurfaceCell::generatePoint
-(
-    const DynamicList<Type>& snappedPoints,
-
-    const scalar s0,
-    const Type& p0,
-    const label p0Index,
-
-    const scalar s1,
-    const Type& p1,
-    const label p1Index
-) const
-{
-    scalar d = s1-s0;
-
-    if (mag(d) > vSmall)
-    {
-        scalar s = (iso_-s0)/d;
-
-        if (s >= 0.5 && s <= 1 && p1Index != -1)
-        {
-            return snappedPoints[p1Index];
-        }
-        else if (s >= 0.0 && s <= 0.5 && p0Index != -1)
-        {
-            return snappedPoints[p0Index];
-        }
-        else
-        {
-            return s*p1 + (1.0-s)*p0;
-        }
-    }
-    else
-    {
-        scalar s = 0.4999;
-
-        return s*p1 + (1.0-s)*p0;
-    }
-}
-
-
-template<class Type>
-void Foam::isoSurfaceCell::generateTriPoints
-(
-    const DynamicList<Type>& snapped,
-
-    const scalar s0,
-    const Type& p0,
-    const label p0Index,
-
-    const scalar s1,
-    const Type& p1,
-    const label p1Index,
-
-    const scalar s2,
-    const Type& p2,
-    const label p2Index,
-
-    const scalar s3,
-    const Type& p3,
-    const label p3Index,
-
-    DynamicList<Type>& pts
-) const
-{
-    int triIndex = 0;
-    if (s0 < iso_)
-    {
-        triIndex |= 1;
-    }
-    if (s1 < iso_)
-    {
-        triIndex |= 2;
-    }
-    if (s2 < iso_)
-    {
-        triIndex |= 4;
-    }
-    if (s3 < iso_)
-    {
-        triIndex |= 8;
-    }
-
-    /* Form the vertices of the triangles for each case */
-    switch (triIndex)
-    {
-        case 0x00:
-        case 0x0F:
-        break;
-
-        case 0x01:
-        case 0x0E:
-        {
-            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)
-            );
-            if (triIndex == 0x0E)
-            {
-                // Flip normals
-                label sz = pts.size();
-                Swap(pts[sz-2], pts[sz-1]);
-            }
-        }
-        break;
-
-        case 0x02:
-        case 0x0D:
-        {
-            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)
-            );
-
-            if (triIndex == 0x0D)
-            {
-                // Flip normals
-                label sz = pts.size();
-                Swap(pts[sz-2], pts[sz-1]);
-            }
-        }
-        break;
-
-        case 0x03:
-        case 0x0C:
-        {
-            Type p0p2 =
-                generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
-            Type p1p3 =
-                generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
-
-            pts.append
-            (
-                generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)
-            );
-            pts.append(p1p3);
-            pts.append(p0p2);
-
-            pts.append(p1p3);
-            pts.append
-            (
-                generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
-            );
-            pts.append(p0p2);
-
-            if (triIndex == 0x0C)
-            {
-                // Flip normals
-                label sz = pts.size();
-                Swap(pts[sz-5], pts[sz-4]);
-                Swap(pts[sz-2], pts[sz-1]);
-            }
-        }
-        break;
-
-        case 0x04:
-        case 0x0B:
-        {
-            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)
-            );
-
-            if (triIndex == 0x0B)
-            {
-                // Flip normals
-                label sz = pts.size();
-                Swap(pts[sz-2], pts[sz-1]);
-            }
-        }
-        break;
-
-        case 0x05:
-        case 0x0A:
-        {
-            Type p0p1 =
-                generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
-            Type p2p3 =
-                generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
-
-            pts.append(p0p1);
-            pts.append(p2p3);
-            pts.append
-            (
-                generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)
-            );
-
-            pts.append(p0p1);
-            pts.append
-            (
-                generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
-            );
-            pts.append(p2p3);
-
-            if (triIndex == 0x0A)
-            {
-                // Flip normals
-                label sz = pts.size();
-                Swap(pts[sz-5], pts[sz-4]);
-                Swap(pts[sz-2], pts[sz-1]);
-            }
-        }
-        break;
-
-        case 0x06:
-        case 0x09:
-        {
-            Type p0p1 =
-                generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
-            Type p2p3 =
-                generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
-
-            pts.append(p0p1);
-            pts.append
-            (
-                generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index)
-            );
-            pts.append(p2p3);
-
-            pts.append(p0p1);
-            pts.append(p2p3);
-            pts.append
-            (
-                generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index)
-            );
-
-            if (triIndex == 0x09)
-            {
-                // Flip normals
-                label sz = pts.size();
-                Swap(pts[sz-5], pts[sz-4]);
-                Swap(pts[sz-2], pts[sz-1]);
-            }
-        }
-        break;
-
-        case 0x08:
-        case 0x07:
-        {
-            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)
-            );
-            if (triIndex == 0x07)
-            {
-                // Flip normals
-                label sz = pts.size();
-                Swap(pts[sz-2], pts[sz-1]);
-            }
-        }
-        break;
-    }
-}
-
-
-template<class Type>
-void Foam::isoSurfaceCell::generateTriPoints
-(
-    const scalarField& cVals,
-    const scalarField& pVals,
-
-    const Field<Type>& cCoords,
-    const Field<Type>& pCoords,
-
-    const DynamicList<Type>& snappedPoints,
-    const labelList& snappedCc,
-    const labelList& snappedPoint,
-
-    DynamicList<Type>& triPoints,
-    DynamicList<label>& triMeshCells
-) const
-{
-    tetMatcher tet;
-    label countNotFoundTets = 0;
-
-    forAll(mesh_.cells(), celli)
-    {
-        if (cellCutType_[celli] != NOTCUT)
-        {
-            label oldNPoints = triPoints.size();
-
-            const cell& cFaces = mesh_.cells()[celli];
-
-            if (tet.isA(mesh_, celli))
-            {
-                // For tets don't do cell-centre decomposition, just use the
-                // tet points and values
-
-                const face& f0 = mesh_.faces()[cFaces[0]];
-
-                // Get the other point
-                const face& f1 = mesh_.faces()[cFaces[1]];
-                label oppositeI = -1;
-                forAll(f1, fp)
-                {
-                    oppositeI = f1[fp];
-
-                    if (findIndex(f0, oppositeI) == -1)
-                    {
-                        break;
-                    }
-                }
-
-                // Start off from positive volume tet to make sure we
-                // generate outwards pointing tets
-                if (mesh_.faceOwner()[cFaces[0]] == celli)
-                {
-                    generateTriPoints
-                    (
-                        snappedPoints,
-
-                        pVals[f0[1]],
-                        pCoords[f0[1]],
-                        snappedPoint[f0[1]],
-
-                        pVals[f0[0]],
-                        pCoords[f0[0]],
-                        snappedPoint[f0[0]],
-
-                        pVals[f0[2]],
-                        pCoords[f0[2]],
-                        snappedPoint[f0[2]],
-
-                        pVals[oppositeI],
-                        pCoords[oppositeI],
-                        snappedPoint[oppositeI],
-
-                        triPoints
-                    );
-                }
-                else
-                {
-                    generateTriPoints
-                    (
-                        snappedPoints,
-
-                        pVals[f0[0]],
-                        pCoords[f0[0]],
-                        snappedPoint[f0[0]],
-
-                        pVals[f0[1]],
-                        pCoords[f0[1]],
-                        snappedPoint[f0[1]],
-
-                        pVals[f0[2]],
-                        pCoords[f0[2]],
-                        snappedPoint[f0[2]],
-
-                        pVals[oppositeI],
-                        pCoords[oppositeI],
-                        snappedPoint[oppositeI],
-
-                        triPoints
-                    );
-                }
-            }
-            else
-            {
-                forAll(cFaces, cFacei)
-                {
-                    label facei = cFaces[cFacei];
-                    const face& f = mesh_.faces()[facei];
-
-                    label fp0 = mesh_.tetBasePtIs()[facei];
-
-                    // Skip undefined tets
-                    if (fp0 < 0)
-                    {
-                        fp0 = 0;
-                        countNotFoundTets++;
-                    }
-
-                    label fp = f.fcIndex(fp0);
-                    for (label i = 2; i < f.size(); i++)
-                    {
-                        label nextFp = f.fcIndex(fp);
-                        triFace tri(f[fp0], f[fp], f[nextFp]);
-
-                        // Start off from positive volume tet to make sure we
-                        // generate outwards pointing tets
-                        if (mesh_.faceOwner()[facei] == celli)
-                        {
-                            generateTriPoints
-                            (
-                                snappedPoints,
-
-                                pVals[tri[1]],
-                                pCoords[tri[1]],
-                                snappedPoint[tri[1]],
-
-                                pVals[tri[0]],
-                                pCoords[tri[0]],
-                                snappedPoint[tri[0]],
-
-                                pVals[tri[2]],
-                                pCoords[tri[2]],
-                                snappedPoint[tri[2]],
-
-                                cVals[celli],
-                                cCoords[celli],
-                                snappedCc[celli],
-
-                                triPoints
-                            );
-                        }
-                        else
-                        {
-                            generateTriPoints
-                            (
-                                snappedPoints,
-
-                                pVals[tri[0]],
-                                pCoords[tri[0]],
-                                snappedPoint[tri[0]],
-
-                                pVals[tri[1]],
-                                pCoords[tri[1]],
-                                snappedPoint[tri[1]],
-
-                                pVals[tri[2]],
-                                pCoords[tri[2]],
-                                snappedPoint[tri[2]],
-
-                                cVals[celli],
-                                cCoords[celli],
-                                snappedCc[celli],
-
-                                triPoints
-                            );
-                        }
-
-                        fp = nextFp;
-                    }
-                }
-            }
-
-
-            // Every three triPoints is a cell
-            label nCells = (triPoints.size()-oldNPoints)/3;
-            for (label i = 0; i < nCells; i++)
-            {
-                triMeshCells.append(celli);
-            }
-        }
-    }
-
-    if (countNotFoundTets > 0)
-    {
-        WarningIn("Foam::isoSurfaceCell::generateTriPoints(..)")
-            << "Could not find " << countNotFoundTets
-            << " tet base points, which may lead to inverted triangles."
-            << endl;
-    }
-
-    triPoints.shrink();
-    triMeshCells.shrink();
-}
-
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>>
-Foam::isoSurfaceCell::interpolate
-(
-    const Field<Type>& cCoords,
-    const Field<Type>& pCoords
-) const
-{
-    DynamicList<Type> triPoints(3*nCutCells_);
-    DynamicList<label> triMeshCells(nCutCells_);
-
-    // Dummy snap data
-    DynamicList<Type> snappedPoints;
-    labelList snappedCc(mesh_.nCells(), -1);
-    labelList snappedPoint(mesh_.nPoints(), -1);
-
-    generateTriPoints
-    (
-        cVals_,
-        pVals_,
-
-        cCoords,
-        pCoords,
-
-        snappedPoints,
-        snappedCc,
-        snappedPoint,
-
-        triPoints,
-        triMeshCells
-    );
-
-    return isoSurface::interpolate
-    (
-        points().size(),
-        triPointMergeMap_,
-        triPoints
-    );
-}
-
-
-// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
deleted file mode 100644
index c46e7d1c79..0000000000
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
+++ /dev/null
@@ -1,368 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "sampledIsoSurfaceCell.H"
-#include "dictionary.H"
-#include "volFields.H"
-#include "volPointInterpolation.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvMesh.H"
-#include "isoSurfaceCell.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(sampledIsoSurfaceCell, 0);
-    addNamedToRunTimeSelectionTable
-    (
-        sampledSurface,
-        sampledIsoSurfaceCell,
-        word,
-        isoSurfaceCell
-    );
-}
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-bool Foam::sampledIsoSurfaceCell::updateGeometry() const
-{
-    const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
-
-    // no update needed
-    if (fvm.time().timeIndex() == prevTimeIndex_)
-    {
-        return false;
-    }
-
-    prevTimeIndex_ = fvm.time().timeIndex();
-
-    // Clear any stored topo
-    facesPtr_.clear();
-
-    // Clear derived data
-    sampledSurface::clearGeom();
-
-    // Optionally read volScalarField
-    autoPtr<volScalarField> readFieldPtr_;
-
-    // 1. see if field in database
-    // 2. see if field can be read
-    const volScalarField* cellFldPtr = nullptr;
-    if (fvm.foundObject<volScalarField>(isoField_))
-    {
-        if (debug)
-        {
-            InfoInFunction << "Lookup " << isoField_ << endl;
-        }
-
-        cellFldPtr = &fvm.lookupObject<volScalarField>(isoField_);
-    }
-    else
-    {
-        // Bit of a hack. Read field and store.
-
-        if (debug)
-        {
-            InfoInFunction
-                << "Reading " << isoField_
-                << " from time " <<fvm.time().timeName()
-                << endl;
-        }
-
-        readFieldPtr_.reset
-        (
-            new volScalarField
-            (
-                IOobject
-                (
-                    isoField_,
-                    fvm.time().timeName(),
-                    fvm,
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE,
-                    false
-                ),
-                fvm
-            )
-        );
-
-        cellFldPtr = readFieldPtr_.operator->();
-    }
-    const volScalarField& cellFld = *cellFldPtr;
-
-    tmp<pointScalarField> pointFld
-    (
-        volPointInterpolation::New(fvm).interpolate(cellFld)
-    );
-
-    if (average_)
-    {
-        //- From point field and interpolated cell.
-        scalarField cellAvg(fvm.nCells(), scalar(0));
-        labelField nPointCells(fvm.nCells(), 0);
-        {
-            for (label pointi = 0; pointi < fvm.nPoints(); pointi++)
-            {
-                const labelList& pCells = fvm.pointCells(pointi);
-
-                forAll(pCells, i)
-                {
-                    label celli = pCells[i];
-
-                    cellAvg[celli] += pointFld().primitiveField()[pointi];
-                    nPointCells[celli]++;
-                }
-            }
-        }
-        forAll(cellAvg, celli)
-        {
-            cellAvg[celli] /= nPointCells[celli];
-        }
-
-        const isoSurfaceCell iso
-        (
-            fvm,
-            cellAvg,
-            pointFld().primitiveField(),
-            isoVal_,
-            regularise_
-        );
-
-        const_cast<sampledIsoSurfaceCell&>
-        (
-            *this
-        ).triSurface::operator=(iso);
-        meshCells_ = iso.meshCells();
-    }
-    else
-    {
-        //- Direct from cell field and point field. Gives bad continuity.
-        const isoSurfaceCell iso
-        (
-            fvm,
-            cellFld.primitiveField(),
-            pointFld().primitiveField(),
-            isoVal_,
-            regularise_
-        );
-
-        const_cast<sampledIsoSurfaceCell&>
-        (
-            *this
-        ).triSurface::operator=(iso);
-        meshCells_ = iso.meshCells();
-    }
-
-
-    if (debug)
-    {
-        Pout<< "sampledIsoSurfaceCell::updateGeometry() : constructed iso:"
-            << nl
-            << "    regularise     : " << regularise_ << nl
-            << "    average        : " << average_ << nl
-            << "    isoField       : " << isoField_ << nl
-            << "    isoValue       : " << isoVal_ << nl
-            << "    points         : " << points().size() << nl
-            << "    tris           : " << triSurface::size() << nl
-            << "    cut cells      : " << meshCells_.size() << endl;
-    }
-
-    return true;
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::sampledIsoSurfaceCell::sampledIsoSurfaceCell
-(
-    const word& name,
-    const polyMesh& mesh,
-    const dictionary& dict
-)
-:
-    sampledSurface(name, mesh, dict),
-    isoField_(dict.lookup("isoField")),
-    isoVal_(readScalar(dict.lookup("isoValue"))),
-    regularise_(dict.lookupOrDefault("regularise", true)),
-    average_(dict.lookupOrDefault("average", true)),
-    zoneKey_(keyType::null),
-    facesPtr_(nullptr),
-    prevTimeIndex_(-1),
-    meshCells_(0)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::sampledIsoSurfaceCell::~sampledIsoSurfaceCell()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::sampledIsoSurfaceCell::needsUpdate() const
-{
-    const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
-
-    return fvm.time().timeIndex() != prevTimeIndex_;
-}
-
-
-bool Foam::sampledIsoSurfaceCell::expire()
-{
-    facesPtr_.clear();
-
-    // Clear derived data
-    sampledSurface::clearGeom();
-
-    // already marked as expired
-    if (prevTimeIndex_ == -1)
-    {
-        return false;
-    }
-
-    // force update
-    prevTimeIndex_ = -1;
-    return true;
-}
-
-
-bool Foam::sampledIsoSurfaceCell::update()
-{
-    return updateGeometry();
-}
-
-
-Foam::tmp<Foam::scalarField>
-Foam::sampledIsoSurfaceCell::sample
-(
-    const volScalarField& vField
-) const
-{
-    return sampleField(vField);
-}
-
-
-Foam::tmp<Foam::vectorField>
-Foam::sampledIsoSurfaceCell::sample
-(
-    const volVectorField& vField
-) const
-{
-    return sampleField(vField);
-}
-
-
-Foam::tmp<Foam::sphericalTensorField>
-Foam::sampledIsoSurfaceCell::sample
-(
-    const volSphericalTensorField& vField
-) const
-{
-    return sampleField(vField);
-}
-
-
-Foam::tmp<Foam::symmTensorField>
-Foam::sampledIsoSurfaceCell::sample
-(
-    const volSymmTensorField& vField
-) const
-{
-    return sampleField(vField);
-}
-
-
-Foam::tmp<Foam::tensorField>
-Foam::sampledIsoSurfaceCell::sample
-(
-    const volTensorField& vField
-) const
-{
-    return sampleField(vField);
-}
-
-
-Foam::tmp<Foam::scalarField>
-Foam::sampledIsoSurfaceCell::interpolate
-(
-    const interpolation<scalar>& interpolator
-) const
-{
-    return interpolateField(interpolator);
-}
-
-
-Foam::tmp<Foam::vectorField>
-Foam::sampledIsoSurfaceCell::interpolate
-(
-    const interpolation<vector>& interpolator
-) const
-{
-    return interpolateField(interpolator);
-}
-
-Foam::tmp<Foam::sphericalTensorField>
-Foam::sampledIsoSurfaceCell::interpolate
-(
-    const interpolation<sphericalTensor>& interpolator
-) const
-{
-    return interpolateField(interpolator);
-}
-
-
-Foam::tmp<Foam::symmTensorField>
-Foam::sampledIsoSurfaceCell::interpolate
-(
-    const interpolation<symmTensor>& interpolator
-) const
-{
-    return interpolateField(interpolator);
-}
-
-
-Foam::tmp<Foam::tensorField>
-Foam::sampledIsoSurfaceCell::interpolate
-(
-    const interpolation<tensor>& interpolator
-) const
-{
-    return interpolateField(interpolator);
-}
-
-
-void Foam::sampledIsoSurfaceCell::print(Ostream& os) const
-{
-    os  << "sampledIsoSurfaceCell: " << name() << " :"
-        << "  field:" << isoField_
-        << "  value:" << isoVal_;
-        //<< "  faces:" << faces().size()   // possibly no geom yet
-        //<< "  points:" << points().size();
-}
-
-
-// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
deleted file mode 100644
index a0db32a2a3..0000000000
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
+++ /dev/null
@@ -1,246 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::sampledIsoSurfaceCell
-
-Description
-    A sampledSurface defined by a surface of iso value. Always triangulated.
-    To be used in sampleSurfaces / functionObjects. Recalculates iso surface
-    only if time changes.
-
-SourceFiles
-    sampledIsoSurfaceCell.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef sampledIsoSurfaceCell_H
-#define sampledIsoSurfaceCell_H
-
-#include "sampledSurface.H"
-#include "triSurface.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                    Class sampledIsoSurfaceCell Declaration
-\*---------------------------------------------------------------------------*/
-
-class sampledIsoSurfaceCell
-:
-    public sampledSurface,
-    public triSurface
-{
-    // Private data
-
-        //- Field to get isoSurface of
-        const word isoField_;
-
-        //- Iso value
-        const scalar isoVal_;
-
-        //- Whether to coarse
-        const Switch regularise_;
-
-        //- Whether to recalculate cell values as average of point values
-        const Switch average_;
-
-        //- If restricted to zones, name of this zone or a regular expression
-        keyType zoneKey_;
-
-        //- Triangles converted to faceList
-        mutable autoPtr<faceList> facesPtr_;
-
-
-        // Recreated for every isoSurface
-
-            //- Time at last call, also track it surface needs an update
-            mutable label prevTimeIndex_;
-
-            //- For every triangle the original cell in mesh
-            mutable labelList meshCells_;
-
-
-    // Private Member Functions
-
-        //- Create iso surface (if time has changed)
-        //  Do nothing (and return false) if no update was needed
-        bool updateGeometry() const;
-
-        //- Sample field on faces
-        template<class Type>
-        tmp<Field<Type>> sampleField
-        (
-            const GeometricField<Type, fvPatchField, volMesh>& vField
-        ) const;
-
-
-        template<class Type>
-        tmp<Field<Type>>
-        interpolateField(const interpolation<Type>&) const;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("sampledIsoSurfaceCell");
-
-
-    // Constructors
-
-        //- Construct from dictionary
-        sampledIsoSurfaceCell
-        (
-            const word& name,
-            const polyMesh& mesh,
-            const dictionary& dict
-        );
-
-
-    //- Destructor
-    virtual ~sampledIsoSurfaceCell();
-
-
-    // Member Functions
-
-        //- Does the surface need an update?
-        virtual bool needsUpdate() const;
-
-        //- Mark the surface as needing an update.
-        //  May also free up unneeded data.
-        //  Return false if surface was already marked as expired.
-        virtual bool expire();
-
-        //- Update the surface as required.
-        //  Do nothing (and return false) if no update was needed
-        virtual bool update();
-
-
-        //- Points of surface
-        virtual const pointField& points() const
-        {
-            return triSurface::points();
-        }
-
-        //- Faces of surface
-        virtual const faceList& faces() const
-        {
-            if (facesPtr_.empty())
-            {
-                const triSurface& s = *this;
-
-                facesPtr_.reset(new faceList(s.size()));
-
-                forAll(s, i)
-                {
-                    facesPtr_()[i] = s[i].triFaceFace();
-                }
-            }
-            return facesPtr_;
-        }
-
-
-        //- Sample field on surface
-        virtual tmp<scalarField> sample
-        (
-            const volScalarField&
-        ) const;
-
-        //- Sample field on surface
-        virtual tmp<vectorField> sample
-        (
-            const volVectorField&
-        ) const;
-
-        //- Sample field on surface
-        virtual tmp<sphericalTensorField> sample
-        (
-            const volSphericalTensorField&
-        ) const;
-
-        //- Sample field on surface
-        virtual tmp<symmTensorField> sample
-        (
-            const volSymmTensorField&
-        ) const;
-
-        //- Sample field on surface
-        virtual tmp<tensorField> sample
-        (
-            const volTensorField&
-        ) const;
-
-
-        //- Interpolate field on surface
-        virtual tmp<scalarField> interpolate
-        (
-            const interpolation<scalar>&
-        ) const;
-
-        //- Interpolate field on surface
-        virtual tmp<vectorField> interpolate
-        (
-            const interpolation<vector>&
-        ) const;
-
-        //- Interpolate field on surface
-        virtual tmp<sphericalTensorField> interpolate
-        (
-            const interpolation<sphericalTensor>&
-        ) const;
-
-        //- Interpolate field on surface
-        virtual tmp<symmTensorField> interpolate
-        (
-            const interpolation<symmTensor>&
-        ) const;
-
-        //- Interpolate field on surface
-        virtual tmp<tensorField> interpolate
-        (
-            const interpolation<tensor>&
-        ) const;
-
-        //- Write
-        virtual void print(Ostream&) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-    #include "sampledIsoSurfaceCellTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C
deleted file mode 100644
index dc46d080e2..0000000000
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C
+++ /dev/null
@@ -1,88 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "sampledIsoSurfaceCell.H"
-#include "isoSurface.H"
-#include "volFieldsFwd.H"
-#include "pointFields.H"
-#include "volPointInterpolation.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>>
-Foam::sampledIsoSurfaceCell::sampleField
-(
-    const GeometricField<Type, fvPatchField, volMesh>& vField
-) const
-{
-    // Recreate geometry if time has changed
-    updateGeometry();
-
-    return tmp<Field<Type>>(new Field<Type>(vField, meshCells_));
-}
-
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>>
-Foam::sampledIsoSurfaceCell::interpolateField
-(
-    const interpolation<Type>& interpolator
-) const
-{
-    // Recreate geometry if time has changed
-    updateGeometry();
-
-    // One value per point
-    tmp<Field<Type>> tvalues(new Field<Type>(points().size()));
-    Field<Type>& values = tvalues.ref();
-
-    boolList pointDone(points().size(), false);
-
-    forAll(faces(), cutFacei)
-    {
-        const face& f = faces()[cutFacei];
-
-        forAll(f, faceVertI)
-        {
-            label pointi = f[faceVertI];
-
-            if (!pointDone[pointi])
-            {
-                values[pointi] = interpolator.interpolate
-                (
-                    points()[pointi],
-                    meshCells_[cutFacei]
-                );
-                pointDone[pointi] = true;
-            }
-        }
-    }
-
-    return tvalues;
-}
-
-
-// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
index 87dda85a16..f3d8f7de7e 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
@@ -233,15 +233,6 @@ void Foam::sampledCuttingPlane::createGeometry()
             regularise_,
             mergeTol_
         )
-        // new isoSurfaceCell
-        //(
-        //    mesh,
-        //    cellDistance,
-        //    pointDistance_,
-        //    0.0,
-        //    regularise_,
-        //    mergeTol_
-        //)
     );
 
     if (debug)
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
index 331479355a..9b4a4bd80b 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
@@ -38,7 +38,6 @@ SourceFiles
 
 #include "sampledSurface.H"
 #include "isoSurface.H"
-//#include "isoSurfaceCell.H"
 #include "plane.H"
 #include "ZoneIDs.H"
 #include "fvMeshSubset.H"
@@ -90,7 +89,6 @@ class sampledCuttingPlane
         scalarField pointDistance_;
 
         //- Constructed iso surface
-        // autoPtr<isoSurfaceCell> isoSurfPtr_;
         autoPtr<isoSurface> isoSurfPtr_;
 
         //- Triangles converted to faceList
@@ -173,8 +171,6 @@ public:
             return facesPtr_;
         }
 
-
-        // const isoSurfaceCell& surface() const
         const isoSurface& surface() const
         {
             return isoSurfPtr_();