diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index 3e42081cc0..dbef3b3987 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -29,6 +29,7 @@ surface/cutting/cuttingSurfaceBaseSelection.C
 surface/distanceSurface/distanceSurface.C
 surface/isoSurface/isoSurface.C
 surface/isoSurface/isoSurfaceCell.C
+surface/isoSurface/isoSurfaceTopo.C
 surface/thresholdCellFaces/thresholdCellFaces.C
 surface/triSurfaceMesh/discreteSurface.C
 
@@ -44,6 +45,7 @@ sampledSurface/sampledPatchInternalField/sampledPatchInternalField.C
 sampledSurface/sampledPlane/sampledPlane.C
 sampledSurface/isoSurface/sampledIsoSurface.C
 sampledSurface/isoSurface/sampledIsoSurfaceCell.C
+sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
 sampledSurface/distanceSurface/sampledDistanceSurface.C
 sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
 sampledSurface/sampledCuttingSurface/sampledCuttingSurface.C
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
new file mode 100644
index 0000000000..6f3ff764f0
--- /dev/null
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
@@ -0,0 +1,328 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "sampledIsoSurfaceTopo.H"
+#include "dictionary.H"
+#include "volFields.H"
+#include "volPointInterpolation.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvMesh.H"
+#include "isoSurfaceTopo.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(sampledIsoSurfaceTopo, 0);
+    addNamedToRunTimeSelectionTable
+    (
+        sampledSurface,
+        sampledIsoSurfaceTopo,
+        word,
+        isoSurfaceTopo
+    );
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::sampledIsoSurfaceTopo::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 derived data
+    sampledSurface::clearGeom();
+
+    // Use field from database, or try to read it in
+
+    const auto* cellFldPtr = fvm.findObject<volScalarField>(isoField_);
+
+    if (debug)
+    {
+        if (cellFldPtr)
+        {
+            InfoInFunction << "Lookup " << isoField_ << endl;
+        }
+        else
+        {
+            InfoInFunction
+                << "Reading " << isoField_
+                << " from time " << fvm.time().timeName()
+                << endl;
+        }
+    }
+
+    // For holding the volScalarField read in.
+    autoPtr<volScalarField> fieldReadPtr;
+
+    if (!cellFldPtr)
+    {
+        // Bit of a hack. Read field and store.
+
+        fieldReadPtr = autoPtr<volScalarField>::New
+        (
+            IOobject
+            (
+                isoField_,
+                fvm.time().timeName(),
+                fvm,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            fvm
+        );
+    }
+
+    const volScalarField& cellFld =
+        (fieldReadPtr.valid() ? *fieldReadPtr : *cellFldPtr);
+
+    auto tpointFld = volPointInterpolation::New(fvm).interpolate(cellFld);
+
+    //- Direct from cell field and point field. Gives bad continuity.
+    isoSurfaceTopo surf
+    (
+        fvm,
+        cellFld.primitiveField(),
+        tpointFld().primitiveField(),
+        isoVal_,
+        (regularise_ ? isoSurfaceTopo::DIAGCELL : isoSurfaceTopo::NONE)
+    );
+
+    MeshedSurface<face>& mySurface = const_cast<sampledIsoSurfaceTopo&>(*this);
+
+    mySurface.transfer(static_cast<meshedSurface&>(surf));
+    meshCells_ = std::move(surf.meshCells());
+
+    // triangulate uses remapFaces()
+    // - this is somewhat less efficient since it recopies the faces
+    // that we just created, but we probably don't want to do this
+    // too often anyhow.
+    if (triangulate_)
+    {
+        labelList faceMap;
+        mySurface.triangulate(faceMap);
+        meshCells_ = UIndirectList<label>(meshCells_, faceMap)();
+    }
+
+    if (debug)
+    {
+        Pout<< "sampledIsoSurfaceTopo::updateGeometry() : constructed iso:"
+            << nl
+            << "    regularise     : " << regularise_ << nl
+            << "    triangulate    : " << triangulate_ << nl
+            << "    isoField       : " << isoField_ << nl
+            << "    isoValue       : " << isoVal_ << nl
+            << "    points         : " << points().size() << nl
+            << "    faces          : " << MeshStorage::size() << nl
+            << "    cut cells      : " << meshCells_.size() << endl;
+    }
+
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sampledIsoSurfaceTopo::sampledIsoSurfaceTopo
+(
+    const word& name,
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    sampledSurface(name, mesh, dict),
+    MeshStorage(),
+    isoField_(dict.get<word>("isoField")),
+    isoVal_(dict.get<scalar>("isoValue")),
+    regularise_(dict.lookupOrDefault("regularise", true)),
+    triangulate_(dict.lookupOrDefault("triangulate", false)),
+    prevTimeIndex_(-1),
+    meshCells_()
+{
+    if (triangulate_ && !regularise_)
+    {
+        FatalIOErrorInFunction(dict) << "Cannot both use regularise"
+            << " and triangulate" << exit(FatalIOError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::sampledIsoSurfaceTopo::~sampledIsoSurfaceTopo()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::sampledIsoSurfaceTopo::needsUpdate() const
+{
+    const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
+
+    return fvm.time().timeIndex() != prevTimeIndex_;
+}
+
+
+bool Foam::sampledIsoSurfaceTopo::expire()
+{
+    // Clear derived data
+    sampledSurface::clearGeom();
+
+    // Already marked as expired
+    if (prevTimeIndex_ == -1)
+    {
+        return false;
+    }
+
+    // Force update
+    prevTimeIndex_ = -1;
+    return true;
+}
+
+
+bool Foam::sampledIsoSurfaceTopo::update()
+{
+    return updateGeometry();
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::sampledIsoSurfaceTopo::sample
+(
+    const interpolation<scalar>& sampler
+) const
+{
+    return sampleOnFaces(sampler);
+}
+
+
+Foam::tmp<Foam::vectorField>
+Foam::sampledIsoSurfaceTopo::sample
+(
+    const interpolation<vector>& sampler
+) const
+{
+    return sampleOnFaces(sampler);
+}
+
+
+Foam::tmp<Foam::sphericalTensorField>
+Foam::sampledIsoSurfaceTopo::sample
+(
+    const interpolation<sphericalTensor>& sampler
+) const
+{
+    return sampleOnFaces(sampler);
+}
+
+
+Foam::tmp<Foam::symmTensorField>
+Foam::sampledIsoSurfaceTopo::sample
+(
+    const interpolation<symmTensor>& sampler
+) const
+{
+    return sampleOnFaces(sampler);
+}
+
+
+Foam::tmp<Foam::tensorField>
+Foam::sampledIsoSurfaceTopo::sample
+(
+    const interpolation<tensor>& sampler
+) const
+{
+    return sampleOnFaces(sampler);
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::sampledIsoSurfaceTopo::interpolate
+(
+    const interpolation<scalar>& interpolator
+) const
+{
+    return sampleOnPoints(interpolator);
+}
+
+
+Foam::tmp<Foam::vectorField>
+Foam::sampledIsoSurfaceTopo::interpolate
+(
+    const interpolation<vector>& interpolator
+) const
+{
+    return sampleOnPoints(interpolator);
+}
+
+Foam::tmp<Foam::sphericalTensorField>
+Foam::sampledIsoSurfaceTopo::interpolate
+(
+    const interpolation<sphericalTensor>& interpolator
+) const
+{
+    return sampleOnPoints(interpolator);
+}
+
+
+Foam::tmp<Foam::symmTensorField>
+Foam::sampledIsoSurfaceTopo::interpolate
+(
+    const interpolation<symmTensor>& interpolator
+) const
+{
+    return sampleOnPoints(interpolator);
+}
+
+
+Foam::tmp<Foam::tensorField>
+Foam::sampledIsoSurfaceTopo::interpolate
+(
+    const interpolation<tensor>& interpolator
+) const
+{
+    return sampleOnPoints(interpolator);
+}
+
+
+void Foam::sampledIsoSurfaceTopo::print(Ostream& os) const
+{
+    os  << "sampledIsoSurfaceTopo: " << name() << " :"
+        << "  field:" << isoField_
+        << "  value:" << isoVal_;
+        //<< "  faces:" << faces().size()   // possibly no geom yet
+        //<< "  points:" << points().size();
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.H
new file mode 100644
index 0000000000..31612b9231
--- /dev/null
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.H
@@ -0,0 +1,291 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::sampledIsoSurfaceTopo
+
+Description
+    A sampledSurface defined by a surface of iso value.
+    To be used in sampleSurfaces / functionObjects. Recalculates iso surface
+    only if time changes.
+
+    This is often embedded as part of a sampled surfaces function object.
+
+Usage
+    Example of function object partial specification:
+    \verbatim
+    surfaces
+    (
+        surface1
+        {
+            type    isoSurfaceTopo;
+            isoField        p;
+            isoValue        0.0;
+        }
+    );
+    \endverbatim
+
+    Where the sub-entries comprise:
+    \table
+        Property | Description                             | Required | Default
+        type     | isoSurfaceTopo                          | yes      |
+        isoField | field name for obtaining iso-surface    | yes      |
+        isoValue | value of iso-surface                    | yes      |
+        regularise | filter faces                          | no       | true
+        triangulate | triangulate faces (if regularise)    | no       | false
+    \endtable
+
+Note
+    Does not currently support cell zones.
+
+SourceFiles
+    sampledIsoSurfaceTopo.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sampledIsoSurfaceTopo_H
+#define sampledIsoSurfaceTopo_H
+
+#include "sampledSurface.H"
+#include "MeshedSurface.H"
+#include "MeshedSurfacesFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class sampledIsoSurfaceTopo Declaration
+\*---------------------------------------------------------------------------*/
+
+class sampledIsoSurfaceTopo
+:
+    public sampledSurface,
+    public MeshedSurface<face>
+{
+    // Private typedefs for convenience
+    typedef MeshedSurface<face> MeshStorage;
+
+    // Private data
+
+        //- Field to get isoSurface of
+        const word isoField_;
+
+        //- Iso value
+        const scalar isoVal_;
+
+        //- Whether to coarse
+        const bool regularise_;
+
+        //- Whether to triangulate
+        const bool triangulate_;
+
+    // Recreated for every isoSurface
+
+        //- Time at last call, also track it surface needs an update
+        mutable label prevTimeIndex_;
+
+        //- For every triangle/face 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 volume field onto surface faces
+        template<class Type>
+        tmp<Field<Type>> sampleOnFaces
+        (
+            const interpolation<Type>& sampler
+        ) const;
+
+        //- Interpolate volume field onto surface points
+        template<class Type>
+        tmp<Field<Type>> sampleOnPoints
+        (
+            const interpolation<Type>& interpolator
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("sampledIsoSurfaceTopo");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        sampledIsoSurfaceTopo
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~sampledIsoSurfaceTopo();
+
+
+    // 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 MeshStorage::points();
+        }
+
+        //- Faces of surface
+        virtual const faceList& faces() const
+        {
+            return *this;
+        }
+
+        //- Const access to per-face zone/region information
+        virtual const labelList& zoneIds() const
+        {
+            return Foam::emptyLabelList;
+        }
+
+        //- Face area magnitudes
+        virtual const vectorField& Sf() const
+        {
+            return MeshStorage::Sf();
+        }
+
+        //- Face area magnitudes
+        virtual const scalarField& magSf() const
+        {
+            return MeshStorage::magSf();
+        }
+
+        //- Face centres
+        virtual const vectorField& Cf() const
+        {
+            return MeshStorage::Cf();
+        }
+
+
+    // Sample
+
+        //- Sample volume field onto surface faces
+        virtual tmp<scalarField> sample
+        (
+            const interpolation<scalar>& sampler
+        ) const;
+
+        //- Sample volume field onto surface faces
+        virtual tmp<vectorField> sample
+        (
+            const interpolation<vector>& sampler
+        ) const;
+
+        //- Sample volume field onto surface faces
+        virtual tmp<sphericalTensorField> sample
+        (
+            const interpolation<sphericalTensor>& sampler
+        ) const;
+
+        //- Sample volume field onto surface faces
+        virtual tmp<symmTensorField> sample
+        (
+            const interpolation<symmTensor>& sampler
+        ) const;
+
+        //- Sample volume field onto surface faces
+        virtual tmp<tensorField> sample
+        (
+            const interpolation<tensor>& sampler
+        ) const;
+
+
+    // Interpolate
+
+        //- Interpolate volume field onto surface points
+        virtual tmp<scalarField> interpolate
+        (
+            const interpolation<scalar>& interpolator
+        ) const;
+
+        //- Interpolate volume field onto surface points
+        virtual tmp<vectorField> interpolate
+        (
+            const interpolation<vector>& interpolator
+        ) const;
+
+        //- Interpolate volume field onto surface points
+        virtual tmp<sphericalTensorField> interpolate
+        (
+            const interpolation<sphericalTensor>& interpolator
+        ) const;
+
+        //- Interpolate volume field onto surface points
+        virtual tmp<symmTensorField> interpolate
+        (
+            const interpolation<symmTensor>& interpolator
+        ) const;
+
+        //- Interpolate volume field onto surface points
+        virtual tmp<tensorField> interpolate
+        (
+            const interpolation<tensor>& interpolator
+        ) const;
+
+        //- Write
+        virtual void print(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "sampledIsoSurfaceTopoTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopoTemplates.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopoTemplates.C
new file mode 100644
index 0000000000..5517417f47
--- /dev/null
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopoTemplates.C
@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "sampledIsoSurfaceTopo.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::sampledIsoSurfaceTopo::sampleOnFaces
+(
+    const interpolation<Type>& sampler
+) const
+{
+    updateGeometry();  // Recreate geometry if time has changed
+
+    return sampledSurface::sampleOnFaces
+    (
+        sampler,
+        meshCells_,
+        faces(),
+        points()
+    );
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::sampledIsoSurfaceTopo::sampleOnPoints
+(
+    const interpolation<Type>& interpolator
+) const
+{
+    updateGeometry();  // Recreate geometry if time has changed
+
+    const labelList& elements = meshCells_;
+
+    // One value per point
+    auto tvalues = tmp<Field<Type>>::New(points().size());
+    auto& values = tvalues.ref();
+
+    const faceList& fcs = faces();
+    const pointField& pts = points();
+
+    bitSet pointDone(points().size());
+
+    forAll(faces(), cutFacei)
+    {
+        const face& f = fcs[cutFacei];
+        const label celli = elements[cutFacei];
+
+        for (const label pointi : f)
+        {
+            if (pointDone.set(pointi))
+            {
+                values[pointi] = interpolator.interpolate
+                (
+                    pts[pointi],
+                    celli
+                );
+            }
+        }
+    }
+
+    return tvalues;
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTopo.C b/src/sampling/surface/isoSurface/isoSurfaceTopo.C
new file mode 100644
index 0000000000..fb0531e423
--- /dev/null
+++ b/src/sampling/surface/isoSurface/isoSurfaceTopo.C
@@ -0,0 +1,1234 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "isoSurfaceTopo.H"
+#include "polyMesh.H"
+#include "tetMatcher.H"
+#include "tetPointRef.H"
+#include "DynamicField.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(isoSurfaceTopo, 0);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::isoSurfaceTopo::isTriCut
+(
+    const triFace& tri,
+    const scalarField& pointValues
+) const
+{
+    const bool aLower = (pointValues[tri[0]] < iso_);
+    const bool bLower = (pointValues[tri[1]] < iso_);
+    const bool cLower = (pointValues[tri[2]] < iso_);
+
+    return !(aLower == bLower && aLower == cLower);
+}
+
+
+Foam::isoSurfaceTopo::cellCutType Foam::isoSurfaceTopo::calcCutType
+(
+    const bool isTet,
+    const label celli
+) const
+{
+    const cell& cFaces = mesh_.cells()[celli];
+
+    if (isTet)
+    {
+        for (const label facei : cFaces)
+        {
+            const face& f = mesh_.faces()[facei];
+
+            for (label fp = 1; fp < f.size() - 1; ++fp)
+            {
+                const triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
+
+                if (isTriCut(tri, pVals_))
+                {
+                    return CUT;
+                }
+            }
+        }
+        return NOTCUT;
+    }
+    else
+    {
+        const bool cellLower = (cVals_[celli] < iso_);
+
+        // First check if there is any cut in cell
+        bool edgeCut = false;
+
+        for (const label facei : cFaces)
+        {
+            const face& f = mesh_.faces()[facei];
+
+            // Check pyramids cut
+            for (const label pointi : f)
+            {
+                if ((pVals_[pointi] < 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]), pVals_))
+                {
+                    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 nPyrEdgeCuts = 0;
+
+            for (const label pointi : cPoints)
+            {
+                if ((pVals_[pointi] < iso_) != cellLower)
+                {
+                    ++nPyrEdgeCuts;
+                }
+            }
+
+            if (nPyrEdgeCuts == cPoints.size())
+            {
+                return SPHERE;
+            }
+            else if (nPyrEdgeCuts)
+            {
+                return CUT;
+            }
+        }
+        return NOTCUT;
+    }
+}
+
+
+Foam::label Foam::isoSurfaceTopo::calcCutTypes
+(
+    tetMatcher& tet,
+    List<cellCutType>& cellCutTypes
+)
+{
+    const cellList& cells = mesh_.cells();
+
+    cellCutTypes.setSize(cells.size());
+    label nCutCells = 0;
+    forAll(cells, celli)
+    {
+        cellCutTypes[celli] = calcCutType(tet.isA(mesh_, celli), celli);
+
+        if (cellCutTypes[celli] == CUT)
+        {
+            ++nCutCells;
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "isoSurfaceCell : candidate cut cells "
+            << nCutCells << " / " << mesh_.nCells() << endl;
+    }
+    return nCutCells;
+}
+
+
+Foam::label Foam::isoSurfaceTopo::generatePoint
+(
+    const label facei,
+    const bool edgeIsDiag,
+    const edge& vertices,
+
+    DynamicList<edge>& pointToVerts,
+    DynamicList<label>& pointToFace,
+    DynamicList<bool>& pointFromDiag,
+    EdgeMap<label>& vertsToPoint
+) const
+{
+    EdgeMap<label>::const_iterator edgeFnd = vertsToPoint.find(vertices);
+    if (edgeFnd != vertsToPoint.end())
+    {
+        return edgeFnd();
+    }
+    else
+    {
+        // Generate new point
+        label pointi = pointToVerts.size();
+
+        pointToVerts.append(vertices);
+        pointToFace.append(facei);
+        pointFromDiag.append(edgeIsDiag);
+        vertsToPoint.insert(vertices, pointi);
+        return pointi;
+    }
+}
+
+
+void Foam::isoSurfaceTopo::generateTriPoints
+(
+    const label facei,
+    const FixedList<scalar, 4>& s,
+    const FixedList<point, 4>& p,
+    const FixedList<label, 4>& pIndex,
+    const FixedList<bool, 6>& edgeIsDiag,// Per tet edge whether is face diag
+
+    DynamicList<edge>& pointToVerts,
+    DynamicList<label>& pointToFace,
+    DynamicList<bool>& pointFromDiag,
+
+    EdgeMap<label>& vertsToPoint,
+    DynamicList<label>& verts       // Every three verts is new triangle
+) const
+{
+    int triIndex = 0;
+    if (s[0] < iso_)
+    {
+        triIndex |= 1;
+    }
+    if (s[1] < iso_)
+    {
+        triIndex |= 2;
+    }
+    if (s[2] < iso_)
+    {
+        triIndex |= 4;
+    }
+    if (s[3] < iso_)
+    {
+        triIndex |= 8;
+    }
+
+    // Form the vertices of the triangles for each case
+    switch (triIndex)
+    {
+        case 0x00:
+        case 0x0F:
+        break;
+
+        case 0x01:
+        case 0x0E:
+        {
+            verts.append
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[0],
+                    edge(pIndex[0], pIndex[1]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+            verts.append
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[1],
+                    edge(pIndex[0], pIndex[2]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+            verts.append
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[2],
+                    edge(pIndex[0], pIndex[3]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+
+            if (triIndex == 0x0E)
+            {
+                // Flip normals
+                label sz = verts.size();
+                Swap(verts[sz-2], verts[sz-1]);
+            }
+        }
+        break;
+
+        case 0x02:
+        case 0x0D:
+        {
+            verts.append
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[0],
+                    edge(pIndex[1], pIndex[0]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+            verts.append
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[3],
+                    edge(pIndex[1], pIndex[3]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+            verts.append
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[4],
+                    edge(pIndex[1], pIndex[2]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+
+            if (triIndex == 0x0D)
+            {
+                // Flip normals
+                label sz = verts.size();
+                Swap(verts[sz-2], verts[sz-1]);
+            }
+        }
+        break;
+
+        case 0x03:
+        case 0x0C:
+        {
+            label p0p2
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[1],
+                    edge(pIndex[0], pIndex[2]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+            label p1p3
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[3],
+                    edge(pIndex[1], pIndex[3]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+
+            verts.append
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[2],
+                    edge(pIndex[0], pIndex[3]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+            verts.append(p1p3);
+            verts.append(p0p2);
+            verts.append(p1p3);
+            verts.append
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[4],
+                    edge(pIndex[1], pIndex[2]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+            verts.append(p0p2);
+
+            if (triIndex == 0x0C)
+            {
+                // Flip normals
+                label sz = verts.size();
+                Swap(verts[sz-5], verts[sz-4]);
+                Swap(verts[sz-2], verts[sz-1]);
+            }
+        }
+        break;
+
+        case 0x04:
+        case 0x0B:
+        {
+            verts.append
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[1],
+                    edge(pIndex[2], pIndex[0]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+            verts.append
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[4],
+                    edge(pIndex[2], pIndex[1]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+            verts.append
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[5],
+                    edge(pIndex[2], pIndex[3]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+
+            if (triIndex == 0x0B)
+            {
+                // Flip normals
+                label sz = verts.size();
+                Swap(verts[sz-2], verts[sz-1]);
+            }
+        }
+        break;
+
+        case 0x05:
+        case 0x0A:
+        {
+            label p0p1
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[0],
+                    edge(pIndex[0], pIndex[1]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+            label p2p3
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[5],
+                    edge(pIndex[2], pIndex[3]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+
+            verts.append(p0p1);
+            verts.append(p2p3);
+            verts.append
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[2],
+                    edge(pIndex[0], pIndex[3]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+            verts.append(p0p1);
+            verts.append
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[4],
+                    edge(pIndex[1], pIndex[2]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+            verts.append(p2p3);
+
+            if (triIndex == 0x0A)
+            {
+                // Flip normals
+                label sz = verts.size();
+                Swap(verts[sz-5], verts[sz-4]);
+                Swap(verts[sz-2], verts[sz-1]);
+            }
+        }
+        break;
+
+        case 0x06:
+        case 0x09:
+        {
+            label p0p1
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[0],
+                    edge(pIndex[0], pIndex[1]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+            label p2p3
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[5],
+                    edge(pIndex[2], pIndex[3]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+
+            verts.append(p0p1);
+            verts.append
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[3],
+                    edge(pIndex[1], pIndex[3]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+            verts.append(p2p3);
+            verts.append(p0p1);
+            verts.append(p2p3);
+            verts.append
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[1],
+                    edge(pIndex[0], pIndex[2]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+
+            if (triIndex == 0x09)
+            {
+                // Flip normals
+                label sz = verts.size();
+                Swap(verts[sz-5], verts[sz-4]);
+                Swap(verts[sz-2], verts[sz-1]);
+            }
+        }
+        break;
+
+        case 0x08:
+        case 0x07:
+        {
+            verts.append
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[2],
+                    edge(pIndex[3], pIndex[0]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+            verts.append
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[5],
+                    edge(pIndex[3], pIndex[2]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+            verts.append
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiag[3],
+                    edge(pIndex[3], pIndex[1]),
+                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                )
+            );
+            if (triIndex == 0x07)
+            {
+                // Flip normals
+                label sz = verts.size();
+                Swap(verts[sz-2], verts[sz-1]);
+            }
+        }
+        break;
+    }
+}
+
+
+void Foam::isoSurfaceTopo::generateTriPoints
+(
+    const polyMesh& mesh,
+    const label celli,
+    const bool isTet,
+
+    DynamicList<edge>& pointToVerts,
+    DynamicList<label>& pointToFace,
+    DynamicList<bool>& pointFromDiag,
+
+    EdgeMap<label>& vertsToPoint,
+    DynamicList<label>& verts,
+    DynamicList<label>& faceLabels
+) const
+{
+    const cell& cFaces = mesh.cells()[celli];
+
+    if (isTet)
+    {
+        // For tets don't do cell-centre decomposition, just use the
+        // tet points and values
+
+        label facei = cFaces[0];
+        const face& f0 = mesh_.faces()[facei];
+
+        // 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;
+            }
+        }
+
+
+        label p0 = f0[0];
+        label p1 = f0[1];
+        label p2 = f0[2];
+        FixedList<bool, 6> edgeIsDiag(false);
+
+        if (mesh.faceOwner()[facei] == celli)
+        {
+            Swap(p1, p2);
+        }
+
+        tetPointRef tet
+        (
+            mesh.points()[p0],
+            mesh.points()[p1],
+            mesh.points()[p2],
+            mesh.points()[oppositeI]
+        );
+
+        label startTrii = verts.size();
+        generateTriPoints
+        (
+            facei,
+            FixedList<scalar, 4>
+            ({
+                pVals_[p0],
+                pVals_[p1],
+                pVals_[p2],
+                pVals_[oppositeI]
+            }),
+            FixedList<point, 4>
+            ({
+                mesh.points()[p0],
+                mesh.points()[p1],
+                mesh.points()[p2],
+                mesh.points()[oppositeI]
+            }),
+            FixedList<label, 4>
+            ({
+                p0,
+                p1,
+                p2,
+                oppositeI
+            }),
+            edgeIsDiag,
+
+            pointToVerts,
+            pointToFace,
+            pointFromDiag,
+            vertsToPoint,
+            verts       // Every three verts is new triangle
+        );
+
+        label nTris = (verts.size()-startTrii)/3;
+        for (label i = 0; i < nTris; ++i)
+        {
+            faceLabels.append(facei);
+        }
+    }
+    else
+    {
+        for (const label facei : cFaces)
+        {
+            const face& f = mesh.faces()[facei];
+
+            label fp0 = mesh.tetBasePtIs()[facei];
+
+            label startTrii = verts.size();
+
+            // Skip undefined tets
+            if (fp0 < 0)
+            {
+                fp0 = 0;
+            }
+
+            label fp = f.fcIndex(fp0);
+            for (label i = 2; i < f.size(); ++i)
+            {
+                label nextFp = f.fcIndex(fp);
+
+                FixedList<bool, 6> edgeIsDiag(false);
+
+                label p0 = f[fp0];
+                label p1 = f[fp];
+                label p2 = f[nextFp];
+                if (mesh.faceOwner()[facei] == celli)
+                {
+                    Swap(p1, p2);
+                    if (i != 2) edgeIsDiag[1] = true;
+                    if (i != f.size()-1) edgeIsDiag[0] = true;
+                }
+                else
+                {
+                    if (i != 2) edgeIsDiag[0] = true;
+                    if (i != f.size()-1) edgeIsDiag[1] = true;
+                }
+
+                tetPointRef tet
+                (
+                    mesh.points()[p0],
+                    mesh.points()[p1],
+                    mesh.points()[p2],
+                    mesh.cellCentres()[celli]
+                );
+
+                generateTriPoints
+                (
+                    facei,
+                    FixedList<scalar, 4>
+                    ({
+                        pVals_[p0],
+                        pVals_[p1],
+                        pVals_[p2],
+                        cVals_[celli]
+                    }),
+                    FixedList<point, 4>
+                    ({
+                        mesh.points()[p0],
+                        mesh.points()[p1],
+                        mesh.points()[p2],
+                        mesh.cellCentres()[celli]
+                    }),
+                    FixedList<label, 4>
+                    ({
+                        p0,
+                        p1,
+                        p2,
+                        mesh.nPoints()+celli
+                    }),
+                    edgeIsDiag,
+
+                    pointToVerts,
+                    pointToFace,
+                    pointFromDiag,
+                    vertsToPoint,
+                    verts       // Every three verts is new triangle
+                );
+
+                fp = nextFp;
+            }
+
+            label nTris = (verts.size()-startTrii)/3;
+            for (label i = 0; i < nTris; ++i)
+            {
+                faceLabels.append(facei);
+            }
+        }
+    }
+}
+
+
+void Foam::isoSurfaceTopo::triangulateOutside
+(
+    const bool filterDiag,
+    const primitivePatch& pp,
+    const boolList& pointFromDiag,
+    const labelList& pointToFace,
+    const label cellID,
+
+    DynamicList<face>& compactFaces,
+    DynamicList<label>& compactCellIDs
+) const
+{
+    // We can form pockets:
+    // - 1. triangle on face
+    // - 2. multiple triangles on interior (from diag edges)
+    // - the edge loop will be pocket since it is only the diag
+    //   edges that give it volume?
+
+    // Retriangulate the exterior loops
+    const labelListList& edgeLoops = pp.edgeLoops();
+    const labelList& mp = pp.meshPoints();
+
+    for (const labelList& loop : edgeLoops)
+    {
+        if (loop.size() > 2)
+        {
+            compactFaces.append(face(0));
+            face& f = compactFaces.last();
+
+            f.setSize(loop.size());
+            label fpi = 0;
+            forAll(f, i)
+            {
+                label pointi = mp[loop[i]];
+                if (filterDiag && pointFromDiag[pointi])
+                {
+                    label prevPointi = mp[loop[loop.fcIndex(i)]];
+                    if
+                    (
+                        pointFromDiag[prevPointi]
+                     && (pointToFace[pointi] != pointToFace[prevPointi])
+                    )
+                    {
+                        f[fpi++] = pointi;
+                    }
+                    else
+                    {
+                        // Filter out diagonal point
+                    }
+                }
+                else
+                {
+                    f[fpi++] = pointi;
+                }
+            }
+
+            if (fpi > 2)
+            {
+                f.setSize(fpi);
+            }
+            else
+            {
+                // Keep original face
+                forAll(f, i)
+                {
+                    label pointi = mp[loop[i]];
+                    f[i] = pointi;
+                }
+            }
+            compactCellIDs.append(cellID);
+        }
+    }
+}
+
+
+Foam::MeshedSurface<Foam::face> Foam::isoSurfaceTopo::removeInsidePoints
+(
+    const bool filterDiag,
+    const MeshStorage& s,
+    const boolList& pointFromDiag,
+    const labelList& pointToFace,
+    const labelList& start,                 // Per cell the starting triangle
+    DynamicList<label>& pointCompactMap,    // Per returned point the original
+    DynamicList<label>& compactCellIDs      // Per returned tri the cellID
+) const
+{
+    const pointField& points = s.points();
+
+    pointCompactMap.clear();
+    compactCellIDs.clear();
+
+    DynamicList<face> compactFaces(s.size()/8);
+
+    for (label celli = 0; celli < start.size()-1; ++celli)
+    {
+        // All triangles of the current cell
+
+        label nTris = start[celli+1]-start[celli];
+
+        if (nTris)
+        {
+            const SubList<face> cellFaces(s, nTris, start[celli]);
+            const primitivePatch pp(cellFaces, points);
+
+            triangulateOutside
+            (
+                filterDiag,
+                pp,
+                pointFromDiag,
+                pointToFace,
+                //protectedFace,
+                celli,
+
+                compactFaces,
+                compactCellIDs
+            );
+        }
+    }
+
+
+    // Compact out unused points
+    // Pick up the used vertices
+    labelList oldToCompact(points.size(), -1);
+    DynamicField<point> compactPoints(points.size());
+    pointCompactMap.clear();
+
+    for (face& f : compactFaces)
+    {
+        forAll(f, fp)
+        {
+            label pointi = f[fp];
+            label compacti = oldToCompact[pointi];
+            if (compacti == -1)
+            {
+                compacti = compactPoints.size();
+                oldToCompact[pointi] = compacti;
+                compactPoints.append(points[pointi]);
+                pointCompactMap.append(pointi);
+            }
+            f[fp] = compacti;
+        }
+    }
+
+
+    MeshStorage cpSurface
+    (
+        std::move(compactPoints),
+        std::move(compactFaces),
+        s.surfZones()
+    );
+
+    return cpSurface;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::isoSurfaceTopo::isoSurfaceTopo
+(
+    const polyMesh& mesh,
+    const scalarField& cVals,
+    const scalarField& pVals,
+    const scalar iso,
+    const filterType filter
+)
+:
+    mesh_(mesh),
+    cVals_(cVals),
+    pVals_(pVals),
+    iso_(iso)
+{
+    if (debug)
+    {
+        Pout<< "isoSurfaceTopo : iso:" << iso_ << " filter:" << filter << endl;
+    }
+
+    tetMatcher tet;
+
+    // Determine if any cut through cell
+    List<cellCutType> cellCutTypes;
+    const label nCutCells = calcCutTypes(tet, cellCutTypes);
+
+    // Per cell: 5 pyramids cut, each generating 2 triangles
+    //  - pointToVerts : from generated iso point to originating mesh verts
+    DynamicList<edge> pointToVerts(10*nCutCells);
+    //  - pointToFace : from generated iso point to originating mesh face
+    DynamicList<label> pointToFace(10*nCutCells);
+    //  - pointToFace : from generated iso point whether is on face diagonal
+    DynamicList<bool> pointFromDiag(10*nCutCells);
+
+    // Per cell: number of intersected edges:
+    //          - four faces cut so 4 mesh edges + 4 face-diagonal edges
+    //          - 4 of the pyramid edges
+    EdgeMap<label> vertsToPoint(12*nCutCells);
+    DynamicList<label> verts(12*nCutCells);
+    // Per cell: 5 pyramids cut (since only one pyramid not cut)
+    DynamicList<label> faceLabels(5*nCutCells);
+    DynamicList<label> cellLabels(5*nCutCells);
+
+
+    labelList startTri(mesh_.nCells()+1, 0);
+
+    for (label celli = 0; celli < mesh_.nCells(); ++celli)
+    {
+        startTri[celli] = faceLabels.size();
+        if (cellCutTypes[celli] != NOTCUT)
+        {
+            generateTriPoints
+            (
+                mesh,
+                celli,
+                tet.isA(mesh_, celli),
+
+                pointToVerts,
+                pointToFace,
+                pointFromDiag,
+
+                vertsToPoint,
+                verts,
+                faceLabels
+            );
+
+            for (label i = startTri[celli]; i < faceLabels.size(); ++i)
+            {
+                cellLabels.append(celli);
+            }
+        }
+    }
+    startTri[mesh_.nCells()] = faceLabels.size();
+
+
+    pointToVerts_.transfer(pointToVerts);
+    meshCells_.transfer(cellLabels);
+    pointToFace_.transfer(pointToFace);
+
+    tmp<pointField> allPoints
+    (
+        interpolate
+        (
+            mesh_.cellCentres(),
+            mesh_.points()
+        )
+    );
+
+
+    // Assign to MeshedSurface
+    faceList allTris(faceLabels.size());
+    label verti = 0;
+    for (face& allTri : allTris)
+    {
+        allTri.setSize(3);
+        allTri[0] = verts[verti++];
+        allTri[1] = verts[verti++];
+        allTri[2] = verts[verti++];
+    }
+
+
+    surfZoneList allZones(1);
+    allZones[0] = surfZone
+    (
+        "allFaces",
+        allTris.size(),     // Size
+        0,                  // Start
+        0                   // Index
+    );
+
+    MeshStorage updated
+    (
+        std::move(allPoints),
+        std::move(allTris),
+        std::move(allZones)
+    );
+    MeshStorage::transfer(updated);
+
+    // Now:
+    // - generated faces and points are assigned to *this
+    // - per point we know:
+    //  - pointOnDiag: whether it is on a face-diagonal edge
+    //  - pointToFace_: from what pyramid (cell+face) it was produced
+    //    (note that the pyramid faces are shared between multiple mesh faces)
+    //  - pointToVerts_ : originating mesh vertex or cell centre
+
+
+    if (debug)
+    {
+        Pout<< "isoSurfaceTopo : generated " << size() << " faces." << endl;
+    }
+
+
+    if (filter != NONE)
+    {
+        // Triangulate outside (filter edges to cell centres and optionally
+        // face diagonals)
+        DynamicList<label> pointCompactMap(size()); // Back to original point
+        DynamicList<label> compactCellIDs(size());  // Per tri the cell
+        MeshStorage::operator=
+        (
+            removeInsidePoints
+            (
+                (filter == DIAGCELL ? true : false),
+                *this,
+                pointFromDiag,
+                pointToFace_,
+                startTri,
+                pointCompactMap,
+                compactCellIDs
+            )
+        );
+
+        pointToVerts_ = UIndirectList<edge>(pointToVerts_, pointCompactMap)();
+        pointToFace_ = UIndirectList<label>(pointToFace_, pointCompactMap)();
+        pointFromDiag = UIndirectList<bool>(pointFromDiag, pointCompactMap)();
+        meshCells_.transfer(compactCellIDs);
+
+        if (debug)
+        {
+            Pout<< "isoSurfaceTopo :"
+                << " after removing cell centre and face-diag triangles : "
+                << size() << endl;
+        }
+
+
+        if (filter == DIAGCELL)
+        {
+            // We remove verts on face diagonals. This is in fact just
+            // straightening the edges of the face through the cell. This can
+            // close off 'pockets' of triangles and create open or
+            // multiply-connected triangles
+
+            // Solved by eroding open-edges
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+            // Mark points on mesh outside. Note that we extend with nCells
+            // so we can easily index with pointToVerts_.
+            PackedBoolList isBoundaryPoint(mesh.nPoints() + mesh.nCells());
+            for
+            (
+                label facei = mesh.nInternalFaces();
+                facei < mesh.nFaces();
+                ++facei
+            )
+            {
+                isBoundaryPoint.set(mesh.faces()[facei]);
+            }
+
+
+            while (true)
+            {
+                const labelList& mp = meshPoints();
+
+                PackedBoolList removeFace(this->size());
+                label nFaces = 0;
+                {
+                    const labelListList& edgeFaces =
+                        MeshStorage::edgeFaces();
+                    forAll(edgeFaces, edgei)
+                    {
+                        const labelList& eFaces = edgeFaces[edgei];
+                        if (eFaces.size() == 1)
+                        {
+                            // Open edge. Check that vertices do not originate
+                            // from a boundary face
+                            const edge& e = edges()[edgei];
+                            const edge& verts0 = pointToVerts_[mp[e[0]]];
+                            const edge& verts1 = pointToVerts_[mp[e[1]]];
+                            if
+                            (
+                                isBoundaryPoint[verts0[0]]
+                             && isBoundaryPoint[verts0[1]]
+                             && isBoundaryPoint[verts1[0]]
+                             && isBoundaryPoint[verts1[1]]
+                            )
+                            {
+                                // Open edge on boundary face. Keep
+                            }
+                            else
+                            {
+                                // Open edge. Mark for erosion
+                                if (removeFace.set(eFaces[0]))
+                                {
+                                    ++nFaces;
+                                }
+                            }
+                        }
+                    }
+                }
+
+                if (debug)
+                {
+                    Pout<< "isoSurfaceTopo :"
+                        << " removing " << nFaces
+                        << " faces since on open edges" << endl;
+                }
+
+                if (returnReduce(nFaces, sumOp<label>()) == 0)
+                {
+                    break;
+                }
+
+                // Remove the faces
+                labelHashSet keepFaces(2*size());
+                forAll(removeFace, facei)
+                {
+                    if (!removeFace[facei])
+                    {
+                        keepFaces.insert(facei);
+                    }
+                }
+
+                labelList pointMap;
+                labelList faceMap;
+                MeshStorage filteredSurf
+                (
+                    MeshStorage::subsetMesh
+                    (
+                        keepFaces,
+                        pointMap,
+                        faceMap
+                    )
+                );
+                MeshStorage::transfer(filteredSurf);
+
+                pointToVerts_ = UIndirectList<edge>(pointToVerts_, pointMap)();
+                pointToFace_ = UIndirectList<label>(pointToFace_, pointMap)();
+                pointFromDiag = UIndirectList<bool>(pointFromDiag, pointMap)();
+                meshCells_ = UIndirectList<label>(meshCells_, faceMap)();
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTopo.H b/src/sampling/surface/isoSurface/isoSurfaceTopo.H
new file mode 100644
index 0000000000..322eab8d7c
--- /dev/null
+++ b/src/sampling/surface/isoSurface/isoSurfaceTopo.H
@@ -0,0 +1,266 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::isoSurfaceTopo
+
+Description
+    Marching tet iso surface algorithm with optional filtering to keep only
+    points originating from mesh edges.
+
+SourceFiles
+    isoSurfaceTopo.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef isoSurfaceTopo_H
+#define isoSurfaceTopo_H
+
+#include "labelPair.H"
+#include "pointIndexHit.H"
+#include "PackedBoolList.H"
+#include "MeshedSurface.H"
+#include "edgeList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyMesh;
+class tetMatcher;
+
+/*---------------------------------------------------------------------------*\
+                       Class isoSurfaceTopo Declaration
+\*---------------------------------------------------------------------------*/
+
+class isoSurfaceTopo
+:
+    public MeshedSurface<face>
+{
+    // Private typedefs for convenience
+    typedef MeshedSurface<face> MeshStorage;
+
+public:
+
+        enum filterType
+        {
+            NONE,       // No filtering
+            DIAGCELL,   // Remove points from face-diagonal and pyramid
+                        // (vertex to cell-centre) edges
+            CELL        // Only remove points from pyramid edges
+        };
+
+
+private:
+
+    // Private data
+
+        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_;
+
+        //- Iso value
+        const scalar iso_;
+
+        //- Per point: originating mesh vertex/cc. See encoding above
+        edgeList pointToVerts_;
+
+        //- For every face the original cell in mesh
+        labelList meshCells_;
+
+        //- For every point the originating face in mesh
+        labelList pointToFace_;
+
+
+    // Private Member Functions
+
+        //- 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 bool isTet,
+            const label
+        ) const;
+
+        //- Determine for all mesh whether cell is cut
+        label calcCutTypes
+        (
+            tetMatcher& tet,
+            List<cellCutType>& cellCutTypes
+        );
+
+        //- Generate single point on edge
+        label generatePoint
+        (
+            const label facei,
+            const bool edgeIsDiag,
+            const edge& vertices,
+
+            DynamicList<edge>& pointToVerts,
+            DynamicList<label>& pointToFace,
+            DynamicList<bool>& pointFromDiag,
+            EdgeMap<label>& vertsToPoint
+        ) const;
+
+        //- Generate triangles from tet
+        void generateTriPoints
+        (
+            const label facei,
+            const FixedList<scalar, 4>& s,
+            const FixedList<point, 4>& p,
+            const FixedList<label, 4>& pIndex,
+            const FixedList<bool, 6>& edgeIsDiag,
+
+            DynamicList<edge>& pointToVerts,
+            DynamicList<label>& pointToFace,
+            DynamicList<bool>& pointFromDiag,
+
+            EdgeMap<label>& vertsToPoint,
+            DynamicList<label>& verts
+        ) const;
+
+        //- Generate triangles from cell
+        void generateTriPoints
+        (
+            const polyMesh& mesh,
+            const label celli,
+            const bool isTet,
+
+            DynamicList<edge>& pointToVerts,
+            DynamicList<label>& pointToFace,
+            DynamicList<bool>& pointFromDiag,
+
+            EdgeMap<label>& vertsToPoint,
+            DynamicList<label>& verts,
+            DynamicList<label>& faceLabels
+        ) const;
+
+
+        // Simplification
+
+            void triangulateOutside
+            (
+                const bool filterDiag,
+                const PrimitivePatch<face, SubList, const pointField&>& pp,
+                const boolList& pointFromDiag,
+                const labelList& pointToFace,
+                const label cellID,
+
+                DynamicList<face>& compactFaces,
+                DynamicList<label>& compactCellIDs
+            ) const;
+
+            MeshStorage removeInsidePoints
+            (
+                const bool filterDiag,
+                const MeshStorage& s,
+                const boolList& pointFromDiag,
+                const labelList& pointToFace,
+                const labelList& start,              // Per cell:starting tri
+                DynamicList<label>& pointCompactMap, // Per point the original
+                DynamicList<label>& compactCellIDs   // Per face the cellID
+            ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("isoSurfaceTopo");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        isoSurfaceTopo
+        (
+            const polyMesh& mesh,
+            const scalarField& cellValues,
+            const scalarField& pointValues,
+            const scalar iso,
+            const filterType filter = DIAGCELL
+        );
+
+
+    // Member Functions
+
+        //- For every face original cell in mesh
+        const labelList& meshCells() const
+        {
+            return meshCells_;
+        }
+
+        //- For every point originating face (pyramid) in mesh
+        const labelList& pointToFace() const
+        {
+            return pointToFace_;
+        }
+
+        //- Per point: originating mesh vertex/cc. See encoding above�
+        const edgeList& pointToVerts() const
+        {
+            return pointToVerts_;
+        }
+
+        //- 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 "isoSurfaceTopoTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTopoTemplates.C b/src/sampling/surface/isoSurface/isoSurfaceTopoTemplates.C
new file mode 100644
index 0000000000..631d948e4c
--- /dev/null
+++ b/src/sampling/surface/isoSurface/isoSurfaceTopoTemplates.C
@@ -0,0 +1,91 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::isoSurfaceTopo::interpolate
+(
+    const Field<Type>& cellCoords,
+    const Field<Type>& pointCoords
+) const
+{
+    tmp<Field<Type>> tfld(new Field<Type>(pointToVerts_.size()));
+    Field<Type>& fld = tfld.ref();
+
+    forAll(pointToVerts_, i)
+    {
+        scalar s0;
+        Type p0;
+        {
+            label v0 = pointToVerts_[i][0];
+            if (v0 < mesh_.nPoints())
+            {
+                s0 = pVals_[v0];
+                p0 = pointCoords[v0];
+            }
+            else
+            {
+                label celli = v0-mesh_.nPoints();
+                s0 = cVals_[celli];
+                p0 = cellCoords[celli];
+            }
+        }
+
+        scalar s1;
+        Type p1;
+        {
+            label v1 = pointToVerts_[i][1];
+            if (v1 < mesh_.nPoints())
+            {
+                s1 = pVals_[v1];
+                p1 = pointCoords[v1];
+            }
+            else
+            {
+                label celli = v1-mesh_.nPoints();
+                s1 = cVals_[celli];
+                p1 = cellCoords[celli];
+            }
+        }
+
+        scalar d = s1-s0;
+        if (mag(d) > VSMALL)
+        {
+            scalar s = (iso_-s0)/d;
+            fld[i] = s*p1+(1.0-s)*p0;
+        }
+        else
+        {
+            fld[i] = 0.5*(p0+p1);
+        }
+    }
+
+    return tfld;
+}
+
+
+// ************************************************************************* //