diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index 2cde386cea..7c2727ce53 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -33,7 +33,6 @@ surface/isoSurface/isoSurfaceBase.C
 surface/isoSurface/isoSurfaceCell.C
 surface/isoSurface/isoSurfaceTopo.C
 surface/thresholdCellFaces/thresholdCellFaces.C
-surface/triSurfaceMesh/discreteSurface.C
 
 sampledSurface/sampledNone/sampledNone.C
 sampledSurface/sampledPatch/sampledPatch.C
@@ -52,9 +51,6 @@ sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
 sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshNormal.C
 sampledSurface/thresholdCellFaces/sampledThresholdCellFaces.C
 
-/* Proof-of-concept: */
-/* sampledSurface/triSurfaceMesh/sampledDiscreteSurface.C */
-
 readers = sampledSurface/readers
 
 $(readers)/surfaceReader.C
diff --git a/src/sampling/sampledSurface/triSurfaceMesh/sampledDiscreteSurface.C b/src/sampling/sampledSurface/triSurfaceMesh/sampledDiscreteSurface.C
deleted file mode 100644
index ca6d87d15a..0000000000
--- a/src/sampling/sampledSurface/triSurfaceMesh/sampledDiscreteSurface.C
+++ /dev/null
@@ -1,236 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-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 "sampledDiscreteSurface.H"
-#include "meshSearch.H"
-#include "Tuple2.H"
-#include "globalIndex.H"
-#include "treeDataCell.H"
-#include "treeDataFace.H"
-#include "meshTools.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(sampledDiscreteSurface, 0);
-    addToRunTimeSelectionTable
-    (
-        sampledSurface,
-        sampledDiscreteSurface,
-        word
-    );
-}
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::sampledDiscreteSurface::sampledDiscreteSurface
-(
-    const word& name,
-    const polyMesh& mesh,
-    const word& surfaceName,
-    const discreteSurface::samplingSource sampleSource
-)
-:
-    sampledSurface(name, mesh),
-    SurfaceSource(mesh, surfaceName, sampleSource)
-{}
-
-
-Foam::sampledDiscreteSurface::sampledDiscreteSurface
-(
-    const word& name,
-    const polyMesh& mesh,
-    const dictionary& dict
-)
-:
-    sampledSurface(name, mesh, dict),
-    SurfaceSource(mesh, dict)
-{}
-
-
-Foam::sampledDiscreteSurface::sampledDiscreteSurface
-(
-    const word& name,
-    const polyMesh& mesh,
-    const triSurface& surface,
-    const word& sampleSourceName
-)
-:
-    sampledSurface(name, mesh),
-    SurfaceSource(name, mesh, surface, sampleSourceName)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::sampledDiscreteSurface::~sampledDiscreteSurface()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::sampledDiscreteSurface::needsUpdate() const
-{
-    return SurfaceSource::needsUpdate();
-}
-
-
-bool Foam::sampledDiscreteSurface::expire()
-{
-    if (SurfaceSource::expire())
-    {
-        // merged information etc
-        sampledSurface::clearGeom();
-
-        return true;
-    }
-
-    return false;
-}
-
-
-bool Foam::sampledDiscreteSurface::update()
-{
-    return SurfaceSource::update();
-}
-
-
-bool Foam::sampledDiscreteSurface::update(const treeBoundBox& bb)
-{
-    return SurfaceSource::update(bb);
-}
-
-
-bool Foam::sampledDiscreteSurface::sampleAndStore
-(
-    const objectRegistry& store,
-    const word& fieldName
-) const
-{
-    return SurfaceSource::sampleAndStore(store, fieldName);
-}
-
-
-Foam::tmp<Foam::scalarField> Foam::sampledDiscreteSurface::sample
-(
-    const interpolation<scalar>& sampler
-) const
-{
-    return SurfaceSource::sampleOnFaces(sampler);
-}
-
-
-Foam::tmp<Foam::vectorField> Foam::sampledDiscreteSurface::sample
-(
-    const interpolation<vector>& sampler
-) const
-{
-    return SurfaceSource::sampleOnFaces(sampler);
-}
-
-Foam::tmp<Foam::sphericalTensorField> Foam::sampledDiscreteSurface::sample
-(
-    const interpolation<sphericalTensor>& sampler
-) const
-{
-    return SurfaceSource::sampleOnFaces(sampler);
-}
-
-
-Foam::tmp<Foam::symmTensorField> Foam::sampledDiscreteSurface::sample
-(
-    const interpolation<symmTensor>& sampler
-) const
-{
-    return SurfaceSource::sampleOnFaces(sampler);
-}
-
-
-Foam::tmp<Foam::tensorField> Foam::sampledDiscreteSurface::sample
-(
-    const interpolation<tensor>& sampler
-) const
-{
-    return SurfaceSource::sampleOnFaces(sampler);
-}
-
-
-Foam::tmp<Foam::scalarField> Foam::sampledDiscreteSurface::interpolate
-(
-    const interpolation<scalar>& interpolator
-) const
-{
-    return SurfaceSource::sampleOnPoints(interpolator);
-}
-
-
-Foam::tmp<Foam::vectorField> Foam::sampledDiscreteSurface::interpolate
-(
-    const interpolation<vector>& interpolator
-) const
-{
-    return SurfaceSource::sampleOnPoints(interpolator);
-}
-
-Foam::tmp<Foam::sphericalTensorField> Foam::sampledDiscreteSurface::interpolate
-(
-    const interpolation<sphericalTensor>& interpolator
-) const
-{
-    return SurfaceSource::sampleOnPoints(interpolator);
-}
-
-
-Foam::tmp<Foam::symmTensorField> Foam::sampledDiscreteSurface::interpolate
-(
-    const interpolation<symmTensor>& interpolator
-) const
-{
-    return SurfaceSource::sampleOnPoints(interpolator);
-}
-
-
-Foam::tmp<Foam::tensorField> Foam::sampledDiscreteSurface::interpolate
-(
-    const interpolation<tensor>& interpolator
-) const
-{
-    return SurfaceSource::sampleOnPoints(interpolator);
-}
-
-
-void Foam::sampledDiscreteSurface::print(Ostream& os) const
-{
-    os  << "sampledDiscreteSurface: " << name() << " :";
-    SurfaceSource::print(os);
-}
-
-
-// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/triSurfaceMesh/sampledDiscreteSurface.H b/src/sampling/sampledSurface/triSurfaceMesh/sampledDiscreteSurface.H
deleted file mode 100644
index adb28da795..0000000000
--- a/src/sampling/sampledSurface/triSurfaceMesh/sampledDiscreteSurface.H
+++ /dev/null
@@ -1,294 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-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::sampledDiscreteSurface
-
-Description
-    A sampledSurface from a triSurfaceMesh.
-    It samples on the points/triangles of a triSurfaceMesh.
-
-See Also
-    discreteSurface, sampledSurface
-
-SourceFiles
-    sampledDiscreteSurface.C
-    sampledDiscreteSurfaceTemplates.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef sampledDiscreteSurface_H
-#define sampledDiscreteSurface_H
-
-#include "sampledSurface.H"
-#include "discreteSurface.H"
-#include "triSurfaceMesh.H"
-#include "MeshedSurface.H"
-#include "MeshedSurfacesFwd.H"
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-class sampledDiscreteSurface;
-
-/*---------------------------------------------------------------------------*\
-                   Class sampledDiscreteSurface Declaration
-\*---------------------------------------------------------------------------*/
-
-class sampledDiscreteSurface
-:
-    public sampledSurface,
-    public discreteSurface
-{
-    //- Private typedefs for convenience
-    typedef discreteSurface MeshStorage;
-    typedef discreteSurface SurfaceSource;
-
-
-    // Private Member Functions
-
-        //- 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("sampledDiscreteSurface");
-
-
-    // Constructors
-
-        //- Construct from components
-        sampledDiscreteSurface
-        (
-            const word& name,
-            const polyMesh& mesh,
-            const word& surfaceName,
-            const discreteSurface::samplingSource sampleSource
-        );
-
-        //- Construct from dictionary
-        sampledDiscreteSurface
-        (
-            const word& name,
-            const polyMesh& mesh,
-            const dictionary& dict
-        );
-
-        //- Construct from triSurface
-        sampledDiscreteSurface
-        (
-            const word& name,
-            const polyMesh& mesh,
-            const triSurface& surface,
-            const word& sampleSourceName
-        );
-
-
-    //- Destructor
-    virtual ~sampledDiscreteSurface();
-
-
-    // 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();
-
-        //- Update the surface using a bound box to limit the searching.
-        //  For direct use, i.e. not through sample.
-        //  Do nothing (and return false) if no update was needed
-        bool update(const treeBoundBox&);
-
-        //- Points of surface
-        virtual const pointField& points() const
-        {
-            return MeshStorage::points();
-        }
-
-        //- Faces of surface
-        virtual const faceList& faces() const
-        {
-            return MeshStorage::surfFaces();
-        }
-
-        //- Per-face zone/region information
-        virtual const labelList& zoneIds() const
-        {
-            return MeshStorage::zoneIds();
-        }
-
-        //- Face area vectors
-        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();
-        }
-
-        //- If element ids/order of the original surface are kept
-        virtual bool hasFaceIds() const
-        {
-            return MeshStorage::hasFaceIds();
-        }
-
-        //- List of element ids/order of the original surface,
-        //  when keepIds is active.
-        virtual const labelList& originalIds() const
-        {
-            return MeshStorage::originalIds();
-        }
-
-        //- Sample the volume field onto surface,
-        //  store it (temporarily) onto the given registry
-        virtual bool sampleAndStore
-        (
-            const objectRegistry& store,
-            const word& fieldName
-        ) const;
-
-
-    // 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;
-
-
-    // Output
-
-        //- Write
-        virtual void print(Ostream& os) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-// #ifdef NoRepository
-//     #include "sampledDiscreteSurfaceTemplates.C"
-// #endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/triSurfaceMesh/sampledDiscreteSurfaceTemplates.C b/src/sampling/sampledSurface/triSurfaceMesh/sampledDiscreteSurfaceTemplates.C
deleted file mode 100644
index 2d4ff8f9dd..0000000000
--- a/src/sampling/sampledSurface/triSurfaceMesh/sampledDiscreteSurfaceTemplates.C
+++ /dev/null
@@ -1,33 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2016 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 "sampledDiscreteSurface.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-
-// ************************************************************************* //
diff --git a/src/sampling/surface/triSurfaceMesh/discreteSurface.C b/src/sampling/surface/triSurfaceMesh/discreteSurface.C
deleted file mode 100644
index 700899d419..0000000000
--- a/src/sampling/surface/triSurfaceMesh/discreteSurface.C
+++ /dev/null
@@ -1,866 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2020 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 "discreteSurface.H"
-#include "meshSearch.H"
-#include "Tuple2.H"
-#include "globalIndex.H"
-#include "treeDataCell.H"
-#include "treeDataFace.H"
-#include "meshTools.H"
-
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-const Foam::Enum
-<
-    Foam::discreteSurface::samplingSource
->
-Foam::discreteSurface::samplingSourceNames_
-({
-    { samplingSource::cells, "cells" },
-    { samplingSource::insideCells, "insideCells" },
-    { samplingSource::boundaryFaces, "boundaryFaces" },
-});
-
-
-namespace Foam
-{
-    defineTypeNameAndDebug(discreteSurface, 0);
-
-    //- Private class for finding nearest
-    //  Comprising:
-    //  - global index
-    //  - sqr(distance)
-    typedef Tuple2<scalar, label> nearInfo;
-
-    class nearestEqOp
-    {
-    public:
-
-        void operator()(nearInfo& x, const nearInfo& y) const
-        {
-            if (y.first() < x.first())
-            {
-                x = y;
-            }
-        }
-    };
-}
-
-
-// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
-
-void Foam::discreteSurface::setZoneMap
-(
-    const surfZoneList& zoneLst,
-    labelList& zoneIds
-)
-{
-    label sz = 0;
-    for (const surfZone& zn : zoneLst)
-    {
-        sz += zn.size();
-    }
-
-    zoneIds.setSize(sz);
-    forAll(zoneLst, zonei)
-    {
-        const surfZone& zn = zoneLst[zonei];
-
-        // Assign sub-zone Ids
-        SubList<label>(zoneIds, zn.size(), zn.start()) = zonei;
-    }
-}
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-const Foam::indexedOctree<Foam::treeDataFace>&
-Foam::discreteSurface::nonCoupledboundaryTree() const
-{
-    // Variant of meshSearch::boundaryTree() that only does non-coupled
-    // boundary faces.
-
-    if (!boundaryTreePtr_.valid())
-    {
-        // all non-coupled boundary faces (not just walls)
-        const polyBoundaryMesh& patches = mesh().boundaryMesh();
-
-        labelList bndFaces(patches.nFaces());
-        label bndI = 0;
-        for (const polyPatch& pp : patches)
-        {
-            if (!pp.coupled())
-            {
-                forAll(pp, i)
-                {
-                    bndFaces[bndI++] = pp.start()+i;
-                }
-            }
-        }
-        bndFaces.setSize(bndI);
-
-
-        treeBoundBox overallBb(mesh().points());
-        Random rndGen(123456);
-        // Extend slightly and make 3D
-        overallBb = overallBb.extend(rndGen, 1e-4);
-        overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-        overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-
-        boundaryTreePtr_.reset
-        (
-            new indexedOctree<treeDataFace>
-            (
-                treeDataFace    // all information needed to search faces
-                (
-                    false,                      // do not cache bb
-                    mesh(),
-                    bndFaces                    // boundary faces only
-                ),
-                overallBb,                      // overall search domain
-                8,                              // maxLevel
-                10,                             // leafsize
-                3.0                             // duplicity
-            )
-        );
-    }
-
-    return *boundaryTreePtr_;
-}
-
-
-bool Foam::discreteSurface::update(const meshSearch& meshSearcher)
-{
-    // Find the cells the triangles of the surface are in.
-    // Does approximation by looking at the face centres only
-    const pointField& fc = surface_.faceCentres();
-
-    List<nearInfo> nearest(fc.size());
-
-    // Global numbering for cells/faces - only used to uniquely identify local
-    // elements
-    globalIndex globalCells(onBoundary() ? mesh().nFaces() : mesh().nCells());
-
-    forAll(nearest, i)
-    {
-        nearest[i].first()  = GREAT;
-        nearest[i].second() = labelMax;
-    }
-
-    if (sampleSource_ == cells)
-    {
-        // Search for nearest cell
-
-        const indexedOctree<treeDataCell>& cellTree = meshSearcher.cellTree();
-
-        forAll(fc, triI)
-        {
-            pointIndexHit nearInfo = cellTree.findNearest
-            (
-                fc[triI],
-                sqr(GREAT)
-            );
-            if (nearInfo.hit())
-            {
-                nearest[triI].first()  = magSqr(nearInfo.hitPoint()-fc[triI]);
-                nearest[triI].second() = globalCells.toGlobal(nearInfo.index());
-            }
-        }
-    }
-    else if (sampleSource_ == insideCells)
-    {
-        // Search for cell containing point
-
-        const indexedOctree<treeDataCell>& cellTree = meshSearcher.cellTree();
-
-        forAll(fc, triI)
-        {
-            if (cellTree.bb().contains(fc[triI]))
-            {
-                const label index = cellTree.findInside(fc[triI]);
-                if (index != -1)
-                {
-                    nearest[triI].first()  = 0.0;
-                    nearest[triI].second() = globalCells.toGlobal(index);
-                }
-            }
-        }
-    }
-    else
-    {
-        // Search for nearest boundaryFace
-
-        ////- Search on all (including coupled) boundary faces
-        //const indexedOctree<treeDataFace>& bTree = meshSearcher.boundaryTree()
-        //- Search on all non-coupled boundary faces
-        const indexedOctree<treeDataFace>& bTree = nonCoupledboundaryTree();
-
-        forAll(fc, triI)
-        {
-            pointIndexHit nearInfo = bTree.findNearest
-            (
-                fc[triI],
-                sqr(GREAT)
-            );
-            if (nearInfo.hit())
-            {
-                nearest[triI].first()  = magSqr(nearInfo.hitPoint()-fc[triI]);
-                nearest[triI].second() = globalCells.toGlobal
-                (
-                    bTree.shapes().faceLabels()[nearInfo.index()]
-                );
-            }
-        }
-    }
-
-
-    // See which processor has the nearest. Mark and subset
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    Pstream::listCombineGather(nearest, nearestEqOp());
-    Pstream::listCombineScatter(nearest);
-
-    labelList cellOrFaceLabels(fc.size(), -1);
-
-    label nFound = 0;
-    forAll(nearest, triI)
-    {
-        if (nearest[triI].second() == labelMax)
-        {
-            // Not found on any processor. How to map?
-        }
-        else if (globalCells.isLocal(nearest[triI].second()))
-        {
-            cellOrFaceLabels[triI] = globalCells.toLocal
-            (
-                nearest[triI].second()
-            );
-            nFound++;
-        }
-    }
-
-
-    if (debug)
-    {
-        Pout<< "Local out of faces:" << cellOrFaceLabels.size()
-            << " keeping:" << nFound << endl;
-    }
-
-    // Now subset the surface. Do not use triSurface::subsetMesh since requires
-    // original surface to be in compact numbering.
-
-    const triSurface& s = surface_;
-
-    // Compact to original triangle
-    labelList faceMap(s.size());
-    // Compact to original points
-    labelList pointMap(s.points().size());
-    // From original point to compact points
-    labelList reversePointMap(s.points().size(), -1);
-
-    // Handle region-wise sorting (makes things slightly more complicated)
-    zoneIds_.setSize(s.size(), -1);
-
-    // Better not to use triSurface::sortedZones here,
-    // since we'll sort ourselves
-
-    // Get zone/region sizes used, store under the original region Id
-    Map<label> zoneSizes;
-
-    // Recover region names from the input surface
-    Map<word> zoneNames;
-    {
-        const geometricSurfacePatchList& patches = s.patches();
-
-        forAll(patches, patchi)
-        {
-            zoneNames.set
-            (
-                patchi,
-                (
-                    patches[patchi].name().empty()
-                  ? geometricSurfacePatch::defaultName(patchi)
-                  : patches[patchi].name()
-                )
-            );
-
-            zoneSizes.set(patchi, 0);
-        }
-    }
-
-
-    {
-        label newPointi = 0;
-        label newFacei = 0;
-
-        forAll(s, facei)
-        {
-            if (cellOrFaceLabels[facei] != -1)
-            {
-                const triSurface::FaceType& f = s[facei];
-                const label regionid = f.region();
-
-                auto fnd = zoneSizes.find(regionid);
-                if (fnd.found())
-                {
-                    ++(*fnd);
-                }
-                else
-                {
-                    // This shouldn't happen
-                    zoneSizes.insert(regionid, 1);
-                    zoneNames.set
-                    (
-                        regionid,
-                        geometricSurfacePatch::defaultName(regionid)
-                    );
-                }
-
-                // Store new faces compact
-                faceMap[newFacei] = facei;
-                zoneIds_[newFacei] = regionid;
-                ++newFacei;
-
-                // Renumber face points
-                forAll(f, fp)
-                {
-                    const label labI = f[fp];
-
-                    if (reversePointMap[labI] == -1)
-                    {
-                        pointMap[newPointi] = labI;
-                        reversePointMap[labI] = newPointi++;
-                    }
-                }
-            }
-        }
-
-        // Trim
-        faceMap.setSize(newFacei);
-        zoneIds_.setSize(newFacei);
-        pointMap.setSize(newPointi);
-    }
-
-
-    // Assign start/size (and name) to the newZones
-    // re-use the lookup to map (zoneId => zoneI)
-    surfZoneList zoneLst(zoneSizes.size());
-    label start = 0;
-    label zoneI = 0;
-    forAllIters(zoneSizes, iter)
-    {
-        // No negative regionids, so Map<label> sorts properly
-        const label regionid = iter.key();
-
-        word name(zoneNames.lookup(regionid, word::null));
-        if (name.empty())
-        {
-            name = geometricSurfacePatch::defaultName(regionid);
-        }
-
-        zoneLst[zoneI] = surfZone
-        (
-            name,
-            0,           // initialize with zero size
-            start,
-            zoneI
-        );
-
-        // Adjust start for the next zone and save (zoneId => zoneI) mapping
-        start += iter.val();
-        iter() = zoneI++;
-    }
-
-
-    // At this stage:
-    // - faceMap to map the (unsorted) compact to original triangle
-    // - zoneIds for the next sorting
-    // - zoneSizes contains region -> count information
-
-    // Rebuild the faceMap for the sorted order
-    labelList sortedFaceMap(faceMap.size());
-
-    forAll(zoneIds_, facei)
-    {
-        label zonei = zoneIds_[facei];
-        label sortedFacei = zoneLst[zonei].start() + zoneLst[zonei].size()++;
-        sortedFaceMap[sortedFacei] = faceMap[facei];
-    }
-
-    // zoneIds are now simply flat values
-    setZoneMap(zoneLst, zoneIds_);
-
-    // Replace the faceMap with the properly sorted face map
-    faceMap.transfer(sortedFaceMap);
-
-    if (keepIds_)
-    {
-        originalIds_ = faceMap;
-    }
-
-    // Subset cellOrFaceLabels (for compact faces)
-    cellOrFaceLabels = labelUIndList(cellOrFaceLabels, faceMap)();
-
-    // Store any face per point (without using pointFaces())
-    labelList pointToFace(pointMap.size());
-
-    // Create faces and points for subsetted surface
-    faceList& surfFaces = this->storedFaces();
-    surfFaces.setSize(faceMap.size());
-
-    this->storedZones().transfer(zoneLst);
-
-    forAll(faceMap, facei)
-    {
-        const labelledTri& origF = s[faceMap[facei]];
-        face& f = surfFaces[facei];
-
-        f = triFace
-        (
-            reversePointMap[origF[0]],
-            reversePointMap[origF[1]],
-            reversePointMap[origF[2]]
-        );
-
-        forAll(f, fp)
-        {
-            pointToFace[f[fp]] = facei;
-        }
-    }
-
-    this->storedPoints() = pointField(s.points(), pointMap);
-
-    if (debug)
-    {
-        print(Pout);
-        Pout<< endl;
-    }
-
-
-    // Collect the samplePoints and sampleElements
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    if (interpolate())
-    {
-        samplePoints_.setSize(pointMap.size());
-        sampleElements_.setSize(pointMap.size(), -1);
-
-        if (sampleSource_ == cells)
-        {
-            // samplePoints_   : per surface point a location inside the cell
-            // sampleElements_ : per surface point the cell
-
-            forAll(points(), pointi)
-            {
-                const point& pt = points()[pointi];
-                label celli = cellOrFaceLabels[pointToFace[pointi]];
-                sampleElements_[pointi] = celli;
-
-                // Check if point inside cell
-                if
-                (
-                    mesh().pointInCell
-                    (
-                        pt,
-                        sampleElements_[pointi],
-                        meshSearcher.decompMode()
-                    )
-                )
-                {
-                    samplePoints_[pointi] = pt;
-                }
-                else
-                {
-                    // Find nearest point on faces of cell
-                    const cell& cFaces = mesh().cells()[celli];
-
-                    scalar minDistSqr = VGREAT;
-
-                    forAll(cFaces, i)
-                    {
-                        const face& f = mesh().faces()[cFaces[i]];
-                        pointHit info = f.nearestPoint(pt, mesh().points());
-                        if (info.distance() < minDistSqr)
-                        {
-                            minDistSqr = info.distance();
-                            samplePoints_[pointi] = info.rawPoint();
-                        }
-                    }
-                }
-            }
-        }
-        else if (sampleSource_ == insideCells)
-        {
-            // samplePoints_   : per surface point a location inside the cell
-            // sampleElements_ : per surface point the cell
-
-            forAll(points(), pointi)
-            {
-                const point& pt = points()[pointi];
-                label celli = cellOrFaceLabels[pointToFace[pointi]];
-                sampleElements_[pointi] = celli;
-                samplePoints_[pointi] = pt;
-            }
-        }
-        else
-        {
-            // samplePoints_   : per surface point a location on the boundary
-            // sampleElements_ : per surface point the boundary face containing
-            //                   the location
-
-            forAll(points(), pointi)
-            {
-                const point& pt = points()[pointi];
-                label facei = cellOrFaceLabels[pointToFace[pointi]];
-                sampleElements_[pointi] = facei;
-                samplePoints_[pointi] = mesh().faces()[facei].nearestPoint
-                (
-                    pt,
-                    mesh().points()
-                ).rawPoint();
-            }
-        }
-    }
-    else
-    {
-        // if sampleSource_ == cells:
-        //      sampleElements_ : per surface triangle the cell
-        //      samplePoints_   : n/a
-        // if sampleSource_ == insideCells:
-        //      sampleElements_ : -1 or per surface triangle the cell
-        //      samplePoints_   : n/a
-        // else:
-        //      sampleElements_ : per surface triangle the boundary face
-        //      samplePoints_   : n/a
-        sampleElements_.transfer(cellOrFaceLabels);
-        samplePoints_.clear();
-    }
-
-
-    if (debug)
-    {
-        this->clearOut();
-        OFstream str(mesh().time().path()/"surfaceToMesh.obj");
-        Info<< "Dumping correspondence from local surface (points or faces)"
-            << " to mesh (cells or faces) to " << str.name() << endl;
-
-        const vectorField& centres =
-        (
-            onBoundary()
-          ? mesh().faceCentres()
-          : mesh().cellCentres()
-        );
-
-        if (interpolate())
-        {
-            label vertI = 0;
-            forAll(samplePoints_, pointi)
-            {
-                meshTools::writeOBJ(str, points()[pointi]);
-                vertI++;
-
-                meshTools::writeOBJ(str, samplePoints_[pointi]);
-                vertI++;
-
-                label elemi = sampleElements_[pointi];
-                meshTools::writeOBJ(str, centres[elemi]);
-                vertI++;
-                str << "l " << vertI-2 << ' ' << vertI-1 << ' ' << vertI << nl;
-            }
-        }
-        else
-        {
-            label vertI = 0;
-            forAll(sampleElements_, triI)
-            {
-                meshTools::writeOBJ(str, faceCentres()[triI]);
-                vertI++;
-
-                label elemi = sampleElements_[triI];
-                meshTools::writeOBJ(str, centres[elemi]);
-                vertI++;
-                str << "l " << vertI-1 << ' ' << vertI << nl;
-            }
-        }
-    }
-
-    needsUpdate_ = false;
-    return true;
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::discreteSurface::discreteSurface
-(
-    const polyMesh& mesh,
-    const word& surfaceName,
-    const samplingSource sampleSource,
-    const bool allowInterpolate
-)
-:
-    MeshStorage(),
-    mesh_(mesh),
-    allowInterpolate_(allowInterpolate),
-    interpolate_(false),
-    surface_
-    (
-        IOobject
-        (
-            surfaceName,
-            mesh.time().constant(), // instance
-            "triSurface",           // local
-            mesh.time(),            // registry
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false
-        )
-    ),
-    sampleSource_(sampleSource),
-    needsUpdate_(true),
-    keepIds_(false),
-    originalIds_(),
-    zoneIds_(),
-    sampleElements_(),
-    samplePoints_()
-{}
-
-
-Foam::discreteSurface::discreteSurface
-(
-    const polyMesh& mesh,
-    const dictionary& dict,
-    const bool allowInterpolate
-)
-:
-    MeshStorage(),
-    mesh_(mesh),
-    allowInterpolate_(allowInterpolate),
-    interpolate_
-    (
-        allowInterpolate
-     && dict.lookupOrDefault("interpolate", false)
-    ),
-    surface_
-    (
-        IOobject
-        (
-            dict.get<word>("surface"),
-            mesh.time().constant(), // instance
-            "triSurface",           // local
-            mesh.time(),            // registry
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false
-        )
-    ),
-    sampleSource_(samplingSourceNames_.get("source", dict)),
-    needsUpdate_(true),
-    keepIds_(dict.lookupOrDefault("keepIds", false)),
-    originalIds_(),
-    zoneIds_(),
-    sampleElements_(),
-    samplePoints_()
-{}
-
-
-Foam::discreteSurface::discreteSurface
-(
-    const word& name,
-    const polyMesh& mesh,
-    const triSurface& surface,
-    const word& sampleSourceName,
-    const bool allowInterpolate
-)
-:
-    MeshStorage(),
-    mesh_(mesh),
-    allowInterpolate_(allowInterpolate),
-    interpolate_(false),
-    surface_
-    (
-        IOobject
-        (
-            name,
-            mesh.time().constant(), // instance
-            "triSurface",           // local
-            mesh.time(),            // registry
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        surface
-    ),
-    sampleSource_(samplingSourceNames_[sampleSourceName]),
-    needsUpdate_(true),
-    keepIds_(false),
-    originalIds_(),
-    zoneIds_(),
-    sampleElements_(),
-    samplePoints_()
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::discreteSurface::~discreteSurface()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::discreteSurface::interpolate() const
-{
-    return interpolate_;
-}
-
-
-bool Foam::discreteSurface::interpolate(bool b)
-{
-    if (interpolate_ == b)
-    {
-        return false;
-    }
-
-    if (b && !allowInterpolate_)
-    {
-        return false;
-    }
-
-    // Value changed, old sampling information is invalid
-    interpolate_ = b;
-    expire();
-
-    return true;
-}
-
-
-bool Foam::discreteSurface::needsUpdate() const
-{
-    return needsUpdate_;
-}
-
-
-bool Foam::discreteSurface::expire()
-{
-    // already marked as expired
-    if (needsUpdate_)
-    {
-        return false;
-    }
-
-    MeshStorage::clear();
-    zoneIds_.clear();
-
-    originalIds_.clear();
-    boundaryTreePtr_.clear();
-    sampleElements_.clear();
-    samplePoints_.clear();
-
-    needsUpdate_ = true;
-    return true;
-}
-
-
-bool Foam::discreteSurface::update()
-{
-    if (!needsUpdate_)
-    {
-        return false;
-    }
-
-    // Calculate surface and mesh overlap bounding box
-    treeBoundBox bb
-    (
-        surface_.triSurface::points(),
-        surface_.triSurface::meshPoints()
-    );
-    bb.min() = max(bb.min(), mesh().bounds().min());
-    bb.max() = min(bb.max(), mesh().bounds().max());
-
-    // Extend a bit
-    const vector span(bb.span());
-
-    bb.min() -= 0.5*span;
-    bb.max() += 0.5*span;
-
-    bb.inflate(1e-6);
-
-    // Mesh search engine, no triangulation of faces.
-    meshSearch meshSearcher(mesh(), bb, polyMesh::FACE_PLANES);
-
-    return update(meshSearcher);
-}
-
-
-bool Foam::discreteSurface::update(const treeBoundBox& bb)
-{
-    if (!needsUpdate_)
-    {
-        return false;
-    }
-
-    // Mesh search engine on subset, no triangulation of faces.
-    meshSearch meshSearcher(mesh(), bb, polyMesh::FACE_PLANES);
-
-    return update(meshSearcher);
-}
-
-
-bool Foam::discreteSurface::sampleAndStore
-(
-    const objectRegistry& store,
-    const word& fieldName,
-    const word& sampleScheme
-) const
-{
-    return
-    (
-        sampleType<scalar>(store, fieldName, sampleScheme)
-     || sampleType<vector>(store, fieldName, sampleScheme)
-     || sampleType<sphericalTensor>(store, fieldName, sampleScheme)
-     || sampleType<symmTensor>(store, fieldName, sampleScheme)
-     || sampleType<tensor>(store, fieldName, sampleScheme)
-    );
-}
-
-
-void Foam::discreteSurface::print(Ostream& os) const
-{
-    os  << "discreteSurface:"
-        << " surface:" << surface_.objectRegistry::name()
-        << " faces:"   << this->surfFaces().size()
-        << " points:"  << this->points().size()
-        << " zoneids:" << this->zoneIds().size();
-}
-
-
-// ************************************************************************* //
diff --git a/src/sampling/surface/triSurfaceMesh/discreteSurface.H b/src/sampling/surface/triSurfaceMesh/discreteSurface.H
deleted file mode 100644
index cf6f4b2bad..0000000000
--- a/src/sampling/surface/triSurfaceMesh/discreteSurface.H
+++ /dev/null
@@ -1,350 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-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::discreteSurface
-
-Description
-    The basis for sampling from triSurfaceMesh.
-    It samples on the points/triangles of the triSurface.
-
-    - it either samples cells or (non-coupled) boundary faces
-
-    - 6 different modes:
-        - source=cells, interpolate=false:
-            finds per triangle centre the nearest cell centre and uses its value
-        - source=cells, interpolate=true
-            finds per triangle centre the nearest cell centre.
-            Per surface point checks if this nearest cell is the one containing
-            point; otherwise projects the point onto the nearest point on
-            the boundary of the cell (to make sure interpolateCellPoint
-            gets a valid location)
-
-        - source=insideCells, interpolate=false:
-            finds per triangle centre the cell containing it and uses its value.
-            Trims triangles outside mesh.
-        - source=insideCells, interpolate=true
-            Per surface point interpolate cell containing it.
-
-        - source=boundaryFaces, interpolate=false:
-            finds per triangle centre the nearest point on the boundary
-            (uncoupled faces only) and uses the value (or 0 if the nearest
-            is on an empty boundary)
-        - source=boundaryFaces, interpolate=true:
-            finds per triangle centre the nearest point on the boundary
-            (uncoupled faces only).
-            Per surface point projects the point onto this boundary face
-            (to make sure interpolateCellPoint gets a valid location)
-
-    - since it finds a nearest per triangle each triangle is guaranteed
-    to be on one processor only. So after stitching (by sampledSurfaces)
-    the original surface should be complete.
-
-SourceFiles
-    discreteSurface.C
-    discreteSurfaceTemplates.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef discreteSurface_H
-#define discreteSurface_H
-
-#include "sampledSurface.H"
-#include "triSurfaceMesh.H"
-#include "MeshedSurface.H"
-#include "MeshedSurfacesFwd.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-class treeDataFace;
-class meshSearch;
-
-/*---------------------------------------------------------------------------*\
-                       Class discreteSurface Declaration
-\*---------------------------------------------------------------------------*/
-
-class discreteSurface
-:
-    public meshedSurface
-{
-public:
-        //- Types of communications
-        enum samplingSource
-        {
-            cells,
-            insideCells,
-            boundaryFaces
-        };
-
-private:
-
-    //- Private typedefs for convenience
-        typedef meshedSurface MeshStorage;
-
-
-    // Private data
-
-        static const Enum<samplingSource> samplingSourceNames_;
-
-        //- Reference to mesh
-        const polyMesh& mesh_;
-
-        //- Do we allow interpolation?
-        const bool allowInterpolate_;
-
-        //- Do we intend to interpolate the information?
-        bool interpolate_;
-
-        //- Surface to sample on
-        const triSurfaceMesh surface_;
-
-        //- Whether to sample internal cell values or boundary values
-        const samplingSource sampleSource_;
-
-        //- Track if the surface needs an update
-        mutable bool needsUpdate_;
-
-        //- Retain element ids/order of original surface
-        bool keepIds_;
-
-        //- List of element ids/order of the original surface,
-        //  when keepIds is active.
-        labelList originalIds_;
-
-        //- Search tree for all non-coupled boundary faces
-        mutable autoPtr<indexedOctree<treeDataFace>> boundaryTreePtr_;
-
-        //- For compatibility with the meshSurf interface
-        labelList zoneIds_;
-
-        //- From local surface triangle to mesh cell/face.
-        labelList sampleElements_;
-
-        //- Local points to sample per point
-        pointField samplePoints_;
-
-
-    // Private Member Functions
-
-        //- Get tree of all non-coupled boundary faces
-        const indexedOctree<treeDataFace>& nonCoupledboundaryTree() const;
-
-        bool update(const meshSearch& meshSearcher);
-
-
-
-protected:
-
-    // Protected Member Functions
-
-        //- Sample the volume field onto surface,
-        //  store it (temporarily) onto the given registry
-        template<class Type>
-        bool sampleType
-        (
-            const objectRegistry& store,
-            const word& fieldName,
-            const word& sampleScheme
-        ) const;
-
-
-        //- Sample field on surface faces
-        template<class Type>
-        tmp<Field<Type>> sampleOnFaces
-        (
-            const interpolation<Type>& sampler
-        ) const;
-
-
-        //- Sample field on surface points
-        template<class Type>
-        tmp<Field<Type>> sampleOnPoints
-        (
-            const interpolation<Type>& interpolator
-        ) const;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("discreteSurface");
-
-
-    // Constructors
-
-        //- Construct from components
-        discreteSurface
-        (
-            const polyMesh& mesh,
-            const word& surfaceName,
-            const samplingSource sampleSource,
-            const bool allowInterpolate = true
-        );
-
-        //- Construct from dictionary
-        discreteSurface
-        (
-            const polyMesh& mesh,
-            const dictionary& dict,
-            const bool allowInterpolate = true
-        );
-
-        //- Construct from triSurface
-        discreteSurface
-        (
-            const word& name,
-            const polyMesh& mesh,
-            const triSurface& surface,
-            const word& sampleSourceName,
-            const bool allowInterpolate = true
-        );
-
-
-    //- Destructor
-    virtual ~discreteSurface();
-
-
-    // Member Functions
-
-        //- Set new zoneIds list based on the surfZoneList information
-        static void setZoneMap(const surfZoneList&, labelList& zoneIds);
-
-
-    // Access
-
-        //- Access to the underlying mesh
-        const polyMesh& mesh() const
-        {
-            return mesh_;
-        }
-
-        //- Interpolation requested for surface
-        bool interpolate() const;
-
-        //- Change interpolation request for surface.
-        //  Return false if the value did not change.
-        //  Use with caution.
-        bool interpolate(bool b);
-
-
-        //- 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();
-
-        //- Update the surface using a bound box to limit the searching.
-        //  For direct use, i.e. not through sample.
-        //  Do nothing (and return false) if no update was needed
-        bool update(const treeBoundBox&);
-
-
-        //- Per-face zone/region information
-        virtual const labelList& zoneIds() const
-        {
-            return zoneIds_;
-        }
-
-        //- If element ids/order of the original surface are kept
-        virtual bool hasFaceIds() const
-        {
-            return keepIds_;
-        }
-
-        //- List of element ids/order of the original surface,
-        //  when keepIds is active.
-        virtual const labelList& originalIds() const
-        {
-            return originalIds_;
-        }
-
-        //- Sampling boundary values instead of cell values
-        bool onBoundary() const
-        {
-            return sampleSource_ == boundaryFaces;
-        }
-
-        //- From local surface face to mesh cell/face.
-        const labelList& sampleElements() const
-        {
-            return sampleElements_;
-        }
-
-        //- Local points to sample per point
-        const pointField& samplePoints() const
-        {
-            return samplePoints_;
-        }
-
-
-        //- Sample the volume field onto surface,
-        //  store it (temporarily) onto the given registry
-        virtual bool sampleAndStore
-        (
-            const objectRegistry& store,
-            const word& fieldName,
-            const word& sampleScheme
-        ) const;
-
-
-        //- Sample field on surface. The tmp field is 'valid' on success.
-        template<class Type>
-        tmp<Field<Type>> sampleType
-        (
-            const word& fieldName,
-            const word& sampleScheme
-        ) const;
-
-
-        //- Write information
-        virtual void print(Ostream& os) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-    #include "discreteSurfaceTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/sampling/surface/triSurfaceMesh/discreteSurfaceTemplates.C b/src/sampling/surface/triSurfaceMesh/discreteSurfaceTemplates.C
deleted file mode 100644
index f724182c01..0000000000
--- a/src/sampling/surface/triSurfaceMesh/discreteSurfaceTemplates.C
+++ /dev/null
@@ -1,252 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2016-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 "discreteSurface.H"
-#include "surfFields.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-template<class Type>
-bool Foam::discreteSurface::sampleType
-(
-    const objectRegistry& obr,
-    const word& fieldName,
-    const word& sampleScheme
-) const
-{
-    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
-    typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
-    typedef IOField<Type> TmpFieldType;
-
-    const auto* volFldPtr = mesh().findObject<VolFieldType>(fieldName);
-
-    if (!volFldPtr)
-    {
-        return false;
-    }
-
-    const auto& volFld = *volFldPtr;
-    auto samplerPtr = interpolation<Type>::New(sampleScheme, *volFldPtr);
-
-    tmp<Field<Type>> tfield = sampleOnFaces(*samplerPtr);
-
-    // The rest could be moved into a separate helper
-    if (isA<surfMesh>(obr))
-    {
-        const surfMesh& surf = dynamicCast<const surfMesh>(obr);
-
-        SurfFieldType* ptr = surf.getObjectPtr<SurfFieldType>(fieldName);
-        if (!ptr)
-        {
-            // Doesn't exist or the wrong type
-            ptr = new SurfFieldType
-            (
-                IOobject
-                (
-                    fieldName,
-                    surf.time().timeName(),
-                    surf,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                surf,
-                dimensioned<Type>(volFld.dimensions(), Zero)
-            );
-            ptr->writeOpt() = IOobject::NO_WRITE;
-
-            surf.store(ptr);
-        }
-
-        ptr->field() = tfield;
-    }
-    else
-    {
-        TmpFieldType* ptr = obr.getObjectPtr<TmpFieldType>(fieldName);
-        if (!ptr)
-        {
-            // Doesn't exist or the wrong type
-            ptr = new TmpFieldType
-            (
-                IOobject
-                (
-                    fieldName,
-                    obr.time().timeName(),
-                    obr,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                tfield().size()
-            );
-            ptr->writeOpt() = IOobject::NO_WRITE;
-
-            obr.store(ptr);
-        }
-
-        *ptr = tfield;
-    }
-
-    return true;
-}
-
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>>
-Foam::discreteSurface::sampleType
-(
-    const word& fieldName,
-    const word& sampleScheme
-) const
-{
-    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
-
-    const auto* volFldPtr = mesh().findObject<VolFieldType>(fieldName);
-
-    if (volFldPtr)
-    {
-        auto samplerPtr = interpolation<Type>::New(sampleScheme, *volFldPtr);
-        return sampleOnFaces(*samplerPtr);
-    }
-
-    return nullptr;
-}
-
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>>
-Foam::discreteSurface::sampleOnFaces
-(
-    const interpolation<Type>& sampler
-) const
-{
-    const labelList& elements = sampleElements_;
-
-    const auto& vField = sampler.psi();
-    const label len = elements.size();
-
-    // One value per face
-    auto tvalues = tmp<Field<Type>>::New(len);
-    auto& values = tvalues.ref();
-
-    if (!onBoundary())
-    {
-        // Sample cells
-
-        const faceList& fcs = this->surfFaces();
-        const pointField& pts = points();
-
-        for (label i=0; i < len; ++i)
-        {
-            const label celli = elements[i];
-            const point pt = fcs[i].centre(pts);
-
-            values[i] = sampler.interpolate(pt, celli);
-        }
-    }
-    else
-    {
-        // Sample boundary faces
-
-        const polyBoundaryMesh& pbm = mesh().boundaryMesh();
-        const label nBnd = mesh().nBoundaryFaces();
-
-        // Create flat boundary field
-
-        Field<Type> bVals(nBnd, Zero);
-
-        const auto& bField = vField.boundaryField();
-
-        forAll(bField, patchi)
-        {
-            const label bFacei = pbm[patchi].start() - mesh().nInternalFaces();
-
-            SubList<Type>
-            (
-                bVals,
-                bField[patchi].size(),
-                bFacei
-            ) = bField[patchi];
-        }
-
-        // Sample in flat boundary field
-
-        for (label i=0; i < len; ++i)
-        {
-            const label bFacei = (elements[i] - mesh().nInternalFaces());
-            values[i] = bVals[bFacei];
-        }
-    }
-
-    return tvalues;
-}
-
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>>
-Foam::discreteSurface::sampleOnPoints
-(
-    const interpolation<Type>& interpolator
-) const
-{
-    // One value per vertex
-    auto tvalues = tmp<Field<Type>>::New(sampleElements_.size());
-    auto& values = tvalues.ref();
-
-    if (!onBoundary())
-    {
-        // Sample cells.
-
-        forAll(sampleElements_, pointi)
-        {
-            values[pointi] = interpolator.interpolate
-            (
-                samplePoints_[pointi],
-                sampleElements_[pointi]
-            );
-        }
-    }
-    else
-    {
-        // Sample boundary faces.
-
-        forAll(samplePoints_, pointi)
-        {
-            const label facei = sampleElements_[pointi];
-
-            values[pointi] = interpolator.interpolate
-            (
-                samplePoints_[pointi],
-                mesh().faceOwner()[facei],
-                facei
-            );
-        }
-    }
-
-    return tvalues;
-}
-
-
-// ************************************************************************* //