ENH: newly structured surface 'sources' for sampling.

- extend the sampling concept to include surfMeshes and surfFields
  for storage.

- Note the createOnRead switch in surfMeshSamplers can be desirable in
  some situations to force creation of the surface faces within the
  constructor.

src/sampling/Make/files

@@ -21,6 +21,12 @@ surface/cuttingPlane/cuttingPlane.C
 surface/isoSurface/isoSurface.C
 surface/isoSurface/isoSurfaceCell.C
 surface/thresholdCellFaces/thresholdCellFaces.C
+surface/triSurfaceMesh/discreteSurface.C
+
+surfMeshSampler/surfMeshSampler/surfMeshSampler.C
+surfMeshSampler/surfMeshSamplers/surfMeshSamplers.C
+surfMeshSampler/plane/surfMeshPlaneSampler.C
+surfMeshSampler/triSurfaceMesh/surfMeshDiscreteSampler.C
 
 sampledSurface/sampledPatch/sampledPatch.C
 sampledSurface/sampledPatchInternalField/sampledPatchInternalField.C
@@ -35,6 +41,10 @@ sampledSurface/sampledSurfaces/sampledSurfacesGrouping.C
 sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
 sampledSurface/thresholdCellFaces/sampledThresholdCellFaces.C
 
+/* Proof-of-concept: */
+/* sampledSurface/triSurfaceMesh/sampledDiscreteSurface.C */
+
+
 surfWriters = sampledSurface/writers
 
 $(surfWriters)/surfaceWriter.C

src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H

@@ -64,6 +64,9 @@ Description
     }
     \endverbatim
 
+See also
+    Foam::surfMeshSamplers
+
 SourceFiles
     sampledSurfaces.C
 

src/sampling/sampledSurface/triSurfaceMesh/sampledDiscreteSurface.C

@@ -0,0 +1,234 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     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"
+#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 volScalarField& vField
+) const
+{
+    return SurfaceSource::sampleField(vField);
+}
+
+
+Foam::tmp<Foam::vectorField> Foam::sampledDiscreteSurface::sample
+(
+    const volVectorField& vField
+) const
+{
+    return SurfaceSource::sampleField(vField);
+}
+
+Foam::tmp<Foam::sphericalTensorField> Foam::sampledDiscreteSurface::sample
+(
+    const volSphericalTensorField& vField
+) const
+{
+    return SurfaceSource::sampleField(vField);
+}
+
+
+Foam::tmp<Foam::symmTensorField> Foam::sampledDiscreteSurface::sample
+(
+    const volSymmTensorField& vField
+) const
+{
+    return SurfaceSource::sampleField(vField);
+}
+
+
+Foam::tmp<Foam::tensorField> Foam::sampledDiscreteSurface::sample
+(
+    const volTensorField& vField
+) const
+{
+    return SurfaceSource::sampleField(vField);
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::sampledDiscreteSurface::interpolate
+(
+    const interpolation<scalar>& interpolator
+) const
+{
+    return SurfaceSource::interpolateField(interpolator);
+}
+
+
+Foam::tmp<Foam::vectorField> Foam::sampledDiscreteSurface::interpolate
+(
+    const interpolation<vector>& interpolator
+) const
+{
+    return SurfaceSource::interpolateField(interpolator);
+}
+
+Foam::tmp<Foam::sphericalTensorField> Foam::sampledDiscreteSurface::interpolate
+(
+    const interpolation<sphericalTensor>& interpolator
+) const
+{
+    return SurfaceSource::interpolateField(interpolator);
+}
+
+
+Foam::tmp<Foam::symmTensorField> Foam::sampledDiscreteSurface::interpolate
+(
+    const interpolation<symmTensor>& interpolator
+) const
+{
+    return SurfaceSource::interpolateField(interpolator);
+}
+
+
+Foam::tmp<Foam::tensorField> Foam::sampledDiscreteSurface::interpolate
+(
+    const interpolation<tensor>& interpolator
+) const
+{
+    return SurfaceSource::interpolateField(interpolator);
+}
+
+
+void Foam::sampledDiscreteSurface::print(Ostream& os) const
+{
+    os  << "sampledDiscreteSurface: " << name() << " :";
+    SurfaceSource::print(os);
+}
+
+
+// ************************************************************************* //

src/sampling/sampledSurface/triSurfaceMesh/sampledDiscreteSurface.H

@@ -0,0 +1,283 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     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/>.
+
+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 field on faces
+        template<class Type>
+        tmp<Field<Type>> sampleField
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& vField
+        ) const;
+
+
+        template<class Type>
+        tmp<Field<Type>>
+        interpolateField(const interpolation<Type>&) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("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();
+        }
+
+        //- Const access to 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
+        bool keepIds() const
+        {
+            return MeshStorage::keepIds();
+        }
+
+        //- List of element ids/order of the original surface,
+        //  when keepIds is active.
+        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 field on surface
+        virtual tmp<scalarField> sample
+        (
+            const volScalarField&
+        ) const;
+
+        //- Sample field on surface
+        virtual tmp<vectorField> sample
+        (
+            const volVectorField&
+        ) const;
+
+        //- Sample field on surface
+        virtual tmp<sphericalTensorField> sample
+        (
+            const volSphericalTensorField&
+        ) const;
+
+        //- Sample field on surface
+        virtual tmp<symmTensorField> sample
+        (
+            const volSymmTensorField&
+        ) const;
+
+        //- Sample field on surface
+        virtual tmp<tensorField> sample
+        (
+            const volTensorField&
+        ) const;
+
+
+        //- Interpolate field on surface
+        virtual tmp<scalarField> interpolate
+        (
+            const interpolation<scalar>&
+        ) const;
+
+
+        //- Interpolate field on surface
+        virtual tmp<vectorField> interpolate
+        (
+            const interpolation<vector>&
+        ) const;
+
+        //- Interpolate field on surface
+        virtual tmp<sphericalTensorField> interpolate
+        (
+            const interpolation<sphericalTensor>&
+        ) const;
+
+        //- Interpolate field on surface
+        virtual tmp<symmTensorField> interpolate
+        (
+            const interpolation<symmTensor>&
+        ) const;
+
+        //- Interpolate field on surface
+        virtual tmp<tensorField> interpolate
+        (
+            const interpolation<tensor>&
+        ) const;
+
+        //- Write information
+        virtual void print(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// #ifdef NoRepository
+//     #include "sampledDiscreteSurfaceTemplates.C"
+// #endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/sampledSurface/triSurfaceMesh/sampledDiscreteSurfaceTemplates.C

@@ -0,0 +1,32 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     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 "sampledDiscreteSurface.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+
+// ************************************************************************* //

src/sampling/surfMeshSampler/plane/surfMeshPlaneSampler.C

@@ -0,0 +1,204 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     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 "surfMeshPlaneSampler.H"
+#include "dictionary.H"
+#include "polyMesh.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(surfMeshPlaneSampler, 0);
+    addNamedToRunTimeSelectionTable
+    (
+        surfMeshSampler,
+        surfMeshPlaneSampler,
+        word,
+        plane
+    );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::surfMeshPlaneSampler::transferContent()
+{
+    SurfaceSource& src = static_cast<SurfaceSource&>(*this);
+    surfMesh& dst = getOrCreateSurfMesh();
+
+    dst.transfer(src);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::surfMeshPlaneSampler::surfMeshPlaneSampler
+(
+    const word& name,
+    const polyMesh& mesh,
+    const plane& planeDesc,
+    const keyType& zoneKey,
+    const bool triangulate
+)
+:
+    surfMeshSampler(name, mesh),
+    SurfaceSource(planeDesc),
+    zoneKey_(zoneKey),
+    triangulate_(triangulate),
+    needsUpdate_(true)
+{
+    if (debug && zoneKey_.size() && mesh.cellZones().findIndex(zoneKey_) < 0)
+    {
+        Info<< "cellZone " << zoneKey_
+            << " not found - using entire mesh" << endl;
+    }
+}
+
+
+Foam::surfMeshPlaneSampler::surfMeshPlaneSampler
+(
+    const word& name,
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    surfMeshSampler(name, mesh, dict),
+    SurfaceSource(plane(dict)),
+    zoneKey_(keyType::null),
+    triangulate_(dict.lookupOrDefault("triangulate", true)),
+    needsUpdate_(true)
+{
+    // Make plane relative to the coordinateSystem (Cartesian)
+    // allow lookup from global coordinate systems
+    if (dict.found("coordinateSystem"))
+    {
+        coordinateSystem cs(mesh, dict.subDict("coordinateSystem"));
+
+        point  base = cs.globalPosition(planeDesc().refPoint());
+        vector norm = cs.globalVector(planeDesc().normal());
+
+        // Assign the plane description
+        static_cast<plane&>(*this) = plane(base, norm);
+    }
+
+    dict.readIfPresent("zone", zoneKey_);
+
+    if (debug && zoneKey_.size() && mesh.cellZones().findIndex(zoneKey_) < 0)
+    {
+        Info<< "cellZone " << zoneKey_
+            << " not found - using entire mesh" << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::surfMeshPlaneSampler::~surfMeshPlaneSampler()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::surfMeshPlaneSampler::needsUpdate() const
+{
+    return needsUpdate_;
+}
+
+
+bool Foam::surfMeshPlaneSampler::expire()
+{
+    // Already marked as expired
+    if (needsUpdate_)
+    {
+        return false;
+    }
+
+    needsUpdate_ = true;
+    return true;
+}
+
+
+bool Foam::surfMeshPlaneSampler::update()
+{
+    if (!needsUpdate_)
+    {
+        return false;
+    }
+
+    labelList selectedCells = mesh().cellZones().findMatching(zoneKey_).used();
+    if (selectedCells.empty())
+    {
+        reCut(mesh(), triangulate_);
+    }
+    else
+    {
+        Foam::sort(selectedCells);
+        reCut(mesh(), triangulate_, selectedCells);
+    }
+
+    if (debug)
+    {
+        print(Pout);
+        Pout<< endl;
+    }
+
+    transferContent();
+
+    needsUpdate_ = false;
+    return true;
+}
+
+
+bool Foam::surfMeshPlaneSampler::sample
+(
+    const word& fieldName
+) const
+{
+    return
+    (
+        sampleType<scalar>(fieldName)
+     || sampleType<vector>(fieldName)
+     || sampleType<sphericalTensor>(fieldName)
+     || sampleType<symmTensor>(fieldName)
+     || sampleType<tensor>(fieldName)
+    );
+}
+
+
+void Foam::surfMeshPlaneSampler::print(Ostream& os) const
+{
+    os  << "surfMeshPlaneSampler: " << name() << " :"
+        << "  base:" << cuttingPlane::refPoint()
+        << "  normal:" << cuttingPlane::normal()
+        << "  triangulate:" << triangulate_
+        << "  faces:"  << SurfaceSource::surfFaces().size()
+        << "  points:" << SurfaceSource::points().size();
+}
+
+
+// ************************************************************************* //

src/sampling/surfMeshSampler/plane/surfMeshPlaneSampler.H

@@ -0,0 +1,166 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     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::surfMeshPlaneSampler
+
+Description
+    Sampling surfFields onto a surfMesh based on a plane.
+    The cuttingPlane algorithm 'cuts' the mesh.
+    The plane is triangulated by default.
+
+Note
+    Does not actually cut until update() called.
+
+SourceFiles
+    surfMeshPlaneSampler.C
+    surfMeshPlaneSamplerTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef surfMeshPlaneSampler_H
+#define surfMeshPlaneSampler_H
+
+#include "surfMeshSampler.H"
+#include "cuttingPlane.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class surfMeshPlaneSampler Declaration
+\*---------------------------------------------------------------------------*/
+
+class surfMeshPlaneSampler
+:
+    public surfMeshSampler,
+    private cuttingPlane
+{
+    // Private typedefs for convenience
+    typedef cuttingPlane SurfaceSource;
+
+    // Private data
+
+        //- If restricted to zones, name of this zone or a regular expression
+        keyType zoneKey_;
+
+        //- Triangulated faces or keep faces as is
+        const bool triangulate_;
+
+        //- Track if the surface needs an update
+        mutable bool needsUpdate_;
+
+
+    // Private Member Functions
+
+        //- Transfer mesh content from SurfaceSource to surfMesh
+        void transferContent();
+
+
+        //- Sample field on surface.
+        template<class Type>
+        bool sampleType
+        (
+            const word& fieldName
+        ) const;
+
+
+        //- Sample field on surface
+        template<class Type>
+        tmp<Field<Type>> sampleField
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& vField
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("surfMeshPlaneSampler");
+
+
+    // Constructors
+
+        //- Construct from components
+        surfMeshPlaneSampler
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const plane& planeDesc,
+            const keyType& zoneKey = word::null,
+            const bool triangulate = true
+        );
+
+        //- Construct from dictionary
+        surfMeshPlaneSampler
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~surfMeshPlaneSampler();
+
+
+    // 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();
+
+        //- Sample the volume field onto surface
+        virtual bool sample(const word& fieldName) const;
+
+
+        //- Write
+        virtual void print(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "surfMeshPlaneSamplerTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/surfMeshSampler/plane/surfMeshPlaneSamplerTemplates.C

@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     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 "surfMeshPlaneSampler.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+bool Foam::surfMeshPlaneSampler::sampleType
+(
+    const word& fieldName
+) const
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
+
+    if (!mesh().foundObject<VolFieldType>(fieldName))
+    {
+        return false;
+    }
+
+    const VolFieldType& fld = mesh().lookupObject<VolFieldType>(fieldName);
+
+    getOrCreateSurfField<Type>(fld).field()
+        = SurfaceSource::sample<Type>(fld);
+
+    return true;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::surfMeshPlaneSampler::sampleField
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vField
+) const
+{
+    return SurfaceSource::sample<Type>(vField);
+}
+
+
+// ************************************************************************* //

src/sampling/surfMeshSampler/surfMeshSampler/surfMeshSampler.C

@@ -0,0 +1,224 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     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 "surfMeshSampler.H"
+#include "MeshedSurfaces.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(surfMeshSampler, 0);
+    defineRunTimeSelectionTable(surfMeshSampler, word);
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+Foam::surfMesh&
+Foam::surfMeshSampler::getOrCreateSurfMesh() const
+{
+    if (!mesh().foundObject<surfMesh>(name()))
+    {
+        surfMesh* ptr = new surfMesh
+        (
+            IOobject
+            (
+                name(),
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            xferCopy(pointField()), // initially no points
+            xferCopy(faceList()),   // initially no faces
+            name()
+        );
+        ptr->setWriteOption(IOobject::NO_WRITE);
+
+        mesh().objectRegistry::store(ptr);
+    }
+
+    return const_cast<surfMesh&>(surface());
+}
+
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::surfMeshSampler>
+Foam::surfMeshSampler::New
+(
+    const word& name,
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+{
+    const word sampleType(dict.lookup("type"));
+
+    wordConstructorTable::iterator cstrIter =
+        wordConstructorTablePtr_->find(sampleType);
+
+    if (cstrIter == wordConstructorTablePtr_->end())
+    {
+        FatalErrorInFunction
+            << "Unknown sample type "
+            << sampleType << nl << nl
+            << "Valid sample types : " << endl
+            << wordConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<surfMeshSampler>(cstrIter()(name, mesh, dict));
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::surfMeshSampler::surfMeshSampler
+(
+    const word& name,
+    const polyMesh& mesh
+)
+:
+    name_(name),
+    mesh_(mesh)
+{}
+
+
+Foam::surfMeshSampler::surfMeshSampler
+(
+    const word& name,
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    name_(name),
+    mesh_(mesh)
+{
+    dict.readIfPresent("name", name_);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::surfMeshSampler::~surfMeshSampler()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::surfMeshSampler::create() const
+{
+    getOrCreateSurfMesh();
+}
+
+
+const Foam::surfMesh& Foam::surfMeshSampler::surface() const
+{
+    return mesh().lookupObject<surfMesh>(name());
+}
+
+
+// Demonstration of using separate tmp registry
+// Foam::label Foam::surfMeshSampler::sample
+// (
+//     const objectRegistry& store,
+//     const UList<word>& fields
+// ) const
+// {
+//     label count = 0;
+//     forAll(fields, fieldi)
+//     {
+//         const word& fieldName = fields[fieldi];
+//         bool ok =
+//         (
+//             sampleAndStore(store, fieldName)
+//          &&
+//             (
+//                 transferField<scalar>(store, fieldName)
+//              || transferField<vector>(store, fieldName)
+//              || transferField<sphericalTensor>(store, fieldName)
+//              || transferField<symmTensor>(store, fieldName)
+//              || transferField<tensor>(store, fieldName)
+//             )
+//         );
+//
+//         if (ok)
+//         {
+//             ++count;
+//         }
+//     }
+//
+//     return count;
+// }
+
+
+Foam::label Foam::surfMeshSampler::sample
+(
+    const UList<word>& fields
+) const
+{
+    label count = 0;
+    forAll(fields, fieldi)
+    {
+        if (sample(fields[fieldi]))
+        {
+            ++count;
+        }
+    }
+
+    return count;
+}
+
+
+Foam::label Foam::surfMeshSampler::write(const wordReList& select) const
+{
+    label count =
+    (
+        writeFields<scalar>(select)
+      + writeFields<vector>(select)
+      + writeFields<sphericalTensor>(select)
+      + writeFields<symmTensor>(select)
+      + writeFields<tensor>(select)
+    );
+
+    if (count)
+    {
+        surface().write();
+    }
+
+    return count;
+}
+
+
+
+void Foam::surfMeshSampler::print(Ostream& os) const
+{
+    os << type();
+}
+
+
+// ************************************************************************* //

src/sampling/surfMeshSampler/surfMeshSampler/surfMeshSampler.H

@@ -0,0 +1,285 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     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::surfMeshSampler
+
+Group
+    grpUtilitiesFunctionObjects
+
+Description
+    An abstract class for surfMesh with sampling.
+
+SourceFiles
+    surfMeshSampler.C
+    surfMeshSamplerTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef surfMeshSampler_H
+#define surfMeshSampler_H
+
+#include "surfMesh.H"
+#include "surfFields.H"
+
+#include "pointField.H"
+#include "word.H"
+#include "labelList.H"
+#include "faceList.H"
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "autoPtr.H"
+#include "volFieldsFwd.H"
+#include "surfaceFieldsFwd.H"
+#include "polyMesh.H"
+#include "coordinateSystems.H"
+#include "interpolation.H"
+#include "error.H"
+#include "IOobjectList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class surfMeshSampler;
+
+/*---------------------------------------------------------------------------*\
+                       Class surfMeshSampler Declaration
+\*---------------------------------------------------------------------------*/
+
+class surfMeshSampler
+{
+    // Private data
+
+        //- Name for surfMesh lookup
+        word name_;
+
+        //- Mesh reference
+        const polyMesh& mesh_;
+
+
+    // Private Member Functions
+
+    // Service methods
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Get existing or create new surfMesh
+        surfMesh& getOrCreateSurfMesh() const;
+
+        //- Get existing or create new surfField.
+        //  Create with same name and dimensions as the 'parent' volField.
+        template<class Type>
+        DimensionedField<Type, surfGeoMesh>&
+        getOrCreateSurfField
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& vField
+        ) const;
+
+
+        // //- Transfer field from object registry to surfField
+        // template<class Type>
+        // bool transferField
+        // (
+        //     const objectRegistry& store,
+        //     const word& fieldName
+        // ) const;
+
+
+        //- Write the given fields
+        template<class Type>
+        label writeFields(const wordReList& select) const;
+
+public:
+
+    //- Runtime type information
+    TypeName("surfMeshSampler");
+
+
+    //- Declare run-time constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        surfMeshSampler,
+        word,
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const dictionary& dict
+        ),
+        (name, mesh, dict)
+    );
+
+
+   // iNew helper class
+
+        //- Class used for the PtrLists read-construction
+        class iNew
+        {
+            //- Reference to the volume mesh
+            const polyMesh& mesh_;
+
+        public:
+
+            iNew(const polyMesh& mesh)
+            :
+                mesh_(mesh)
+            {}
+
+            autoPtr<surfMeshSampler> operator()(Istream& is) const
+            {
+                word name(is);
+                dictionary dict(is);
+
+                return surfMeshSampler::New(name, mesh_, dict);
+            }
+        };
+
+
+    // Constructors
+
+        //- Construct from name, mesh
+        surfMeshSampler
+        (
+            const word& name,
+            const polyMesh& mesh
+        );
+
+        //- Construct from dictionary
+        surfMeshSampler
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+
+        //- Clone
+        autoPtr<surfMeshSampler> clone() const
+        {
+            NotImplemented;
+            return autoPtr<surfMeshSampler>(nullptr);
+        }
+
+
+
+    // Selectors
+
+        //- Return a reference to the selected surface
+        static autoPtr<surfMeshSampler> New
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~surfMeshSampler();
+
+
+    // Member Functions
+
+      // Access
+
+        //- Access to the underlying mesh
+        const polyMesh& mesh() const
+        {
+            return mesh_;
+        }
+
+        //- Name of surface
+        const word& name() const
+        {
+            return name_;
+        }
+
+
+        //- Create surfMesh if required.
+        void create() const;
+
+        //- Return existing surfMesh.
+        virtual const surfMesh& surface() const;
+
+
+        //- Does the surface need an update?
+        virtual bool needsUpdate() const = 0;
+
+        //- 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() = 0;
+
+        //- Update the surface as required.
+        //  Do nothing (and return false) if no update was needed
+        virtual bool update() = 0;
+
+
+      // Edit
+
+        //- Rename
+        virtual void rename(const word& newName)
+        {
+            name_ = newName;
+        }
+
+
+      // Execute
+
+        //- Sample from specified volume field to surface
+        virtual bool sample(const word& fieldName) const = 0;
+
+        //- Sample from volume fields to specified surface fields.
+        virtual label sample(const UList<word>& fields) const;
+
+
+      // Write
+
+        //- Write specified fields
+        virtual label write(const wordReList& fieldSelection) const;
+
+        //- Write
+        virtual void print(Ostream&) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "surfMeshSamplerTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/surfMeshSampler/surfMeshSampler/surfMeshSamplerTemplates.C

@@ -0,0 +1,125 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     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 "surfMeshSampler.H"
+#include "dimensionedType.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+Foam::DimensionedField<Type, Foam::surfGeoMesh>&
+Foam::surfMeshSampler::getOrCreateSurfField
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vField
+) const
+{
+    typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
+
+    const surfMesh& surf  = surface();
+    const word& fieldName = vField.name();
+
+    SurfFieldType* ptr = surf.lookupObjectRefPtr<SurfFieldType>(fieldName);
+    if (!ptr)
+    {
+        ptr = new SurfFieldType
+        (
+            IOobject
+            (
+                fieldName,
+                surf.time().timeName(),
+                surf,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            surf,
+            dimensioned<Type>("0", vField.dimensions(), Zero)
+        );
+        ptr->writeOpt() = IOobject::NO_WRITE;
+
+        surf.store(ptr);
+    }
+
+    return *ptr;
+}
+
+
+// // Older code for transferring an IOField to a surfField between
+// // different registries
+// template<class Type>
+// bool Foam::surfMeshSampler::transferField
+// (
+//     const objectRegistry& store,
+//     const word& fieldName
+// ) const
+// {
+//     typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
+//     typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
+//     typedef IOField<Type> TmpFieldType;
+//
+//     // foundObject includes a type check
+//     bool ok =
+//     (
+//         mesh_.foundObject<VolFieldType>(fieldName)
+//      && store.foundObject<TmpFieldType>(fieldName)
+//     );
+//
+//     if (ok)
+//     {
+//         SurfFieldType& sfield = getOrCreateSurfField
+//         (
+//             mesh_.lookupObject<VolFieldType>(fieldName)
+//         );
+//
+//         TmpFieldType& iofield =
+//             store.lookupObjectRef<TmpFieldType>(fieldName);
+//
+//         sfield.transfer(iofield);
+//         store.checkOut(iofield);
+//     }
+//
+//     return ok;
+// }
+
+
+template<class Type>
+Foam::label Foam::surfMeshSampler::writeFields
+(
+    const wordReList& select
+) const
+{
+    typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
+    const surfMesh& s = surface();
+
+    wordList names = s.sortedNames<SurfFieldType>(select);
+    forAll(names, namei)
+    {
+        s.lookupObject<SurfFieldType>(names[namei]).write();
+    }
+
+    return names.size();
+}
+
+
+// ************************************************************************* //

src/sampling/surfMeshSampler/surfMeshSamplers/surfMeshSamplers.C

@@ -0,0 +1,421 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     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 "surfMeshSamplers.H"
+#include "volFields.H"
+#include "dictionary.H"
+#include "Time.H"
+#include "IOmanip.H"
+#include "volPointInterpolation.H"
+#include "PatchTools.H"
+#include "mapPolyMesh.H"
+#include "wordReListMatcher.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(surfMeshSamplers, 0);
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        surfMeshSamplers,
+        dictionary
+    );
+}
+
+
+bool Foam::surfMeshSamplers::verbose_ = false;
+
+void Foam::surfMeshSamplers::checkOutNames
+(
+    const objectRegistry& registry,
+    const UList<word>& names
+)
+{
+    objectRegistry& reg = const_cast<objectRegistry&>(registry);
+
+    forAll(names, namei)
+    {
+        objectRegistry::iterator iter = reg.find(names[namei]);
+        if (iter != reg.end())
+        {
+            registry.checkOut(*iter());
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// //- Temporary object registry for passing around values
+// const objectRegistry& tmpRegistry() const;
+//
+// //- Temporary object registry for passing around values
+// const objectRegistry& tmpRegistry(const word& subName) const;
+
+// const Foam::objectRegistry&
+// Foam::surfMeshSamplers::tmpRegistry() const
+// {
+//     // Sub-registry for sampling, choose name for fewer collisions
+//     return mesh_.thisDb().subRegistry
+//     (
+//         "$tmp$" + type() + "$" + name(),
+//         true,
+//         false
+//     );
+// }
+//
+//
+// const Foam::objectRegistry&
+// Foam::surfMeshSamplers::tmpRegistry(const word& subName) const
+// {
+//     return tmpRegistry().subRegistry(subName, true, false);
+// }
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::surfMeshSamplers::surfMeshSamplers
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    functionObjects::regionFunctionObject(name, runTime, dict),
+    PtrList<surfMeshSampler>(),
+    mesh_(refCast<const fvMesh>(obr_)),
+    fieldSelection_()
+{
+    read(dict);
+}
+
+
+Foam::surfMeshSamplers::surfMeshSamplers
+(
+    const word& name,
+    const objectRegistry& obr,
+    const dictionary& dict
+)
+:
+    functionObjects::regionFunctionObject(name, obr, dict),
+    PtrList<surfMeshSampler>(),
+    mesh_(refCast<const fvMesh>(obr)),
+    fieldSelection_(),
+    derivedNames_()
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::surfMeshSamplers::~surfMeshSamplers()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::surfMeshSamplers::verbose(const bool verbosity)
+{
+    verbose_ = verbosity;
+}
+
+
+bool Foam::surfMeshSamplers::execute()
+{
+    if (size())
+    {
+        const objectRegistry& db = mesh_.thisDb();
+
+        // Manage derived names
+        DynamicList<word> added(derivedNames_.size());
+        DynamicList<word> cleanup(derivedNames_.size());
+
+        forAll(derivedNames_, namei)
+        {
+            const word& derivedName = derivedNames_[namei];
+
+            if (derivedName == "rhoU")
+            {
+                added.append(derivedName);
+
+                if (!db.foundObject<volVectorField>(derivedName))
+                {
+                    cleanup.append(derivedName);
+
+                    db.store
+                    (
+                        new volVectorField
+                        (
+                            derivedName,
+                            // rhoU = rho * U
+                            (
+                                mesh_.lookupObject<volScalarField>("rho")
+                              * mesh_.lookupObject<volVectorField>("U")
+                            )
+                        )
+                    );
+                }
+            }
+            else if (derivedName == "pTotal")
+            {
+                added.append(derivedName);
+
+                if (!db.foundObject<volScalarField>(derivedName))
+                {
+                    cleanup.append(derivedName);
+
+                    db.store
+                    (
+                        new volScalarField
+                        (
+                            derivedName,
+                            // pTotal = p + U^2 / 2
+                            (
+                                mesh_.lookupObject<volScalarField>("p")
+                              + 0.5
+                              * mesh_.lookupObject<volScalarField>("rho")
+                              * magSqr(mesh_.lookupObject<volVectorField>("U"))
+                            )
+                        )
+                    );
+                }
+            }
+            else
+            {
+                WarningInFunction
+                    << "unknown derived name: " << derivedName << nl
+                    << "Use one of 'rhoU', 'pTotal'" << nl
+                    << endl;
+            }
+        }
+
+        // The acceptable fields
+        wordHashSet acceptable(added);
+        acceptable.insert(acceptType<scalar>());
+        acceptable.insert(acceptType<vector>());
+        acceptable.insert(acceptType<sphericalTensor>());
+        acceptable.insert(acceptType<symmTensor>());
+        acceptable.insert(acceptType<tensor>());
+
+        const wordList fields = acceptable.sortedToc();
+        if (!fields.empty())
+        {
+            forAll(*this, surfI)
+            {
+                surfMeshSampler& s = operator[](surfI);
+
+                // Potentially monitor the update for writing geometry?
+                if (s.needsUpdate())
+                {
+                    s.update();
+                }
+
+                s.sample(fields);
+            }
+        }
+
+        checkOutNames(db, cleanup);
+    }
+
+    return true;
+}
+
+
+bool Foam::surfMeshSamplers::write()
+{
+    // Write sampled fields (on surface)
+    //
+    // Doesn't bother checking which fields have been generated here
+    // or elsewhere
+
+    // This could be more efficient
+    wordReList select(fieldSelection_.size() + derivedNames_.size());
+
+    label nElem = 0;
+    forAll(fieldSelection_, i)
+    {
+        select[nElem++] = fieldSelection_[i];
+    }
+    forAll(derivedNames_, i)
+    {
+        select[nElem++] = derivedNames_[i];
+    }
+
+    // avoid duplicate entries
+    select = wordReListMatcher::uniq(select);
+
+    forAll(*this, surfI)
+    {
+        const surfMeshSampler& s = operator[](surfI);
+        s.write(select);
+    }
+
+    return true;
+}
+
+
+bool Foam::surfMeshSamplers::read(const dictionary& dict)
+{
+    derivedNames_.clear();
+
+    const bool createOnRead =
+        dict.lookupOrDefault<Switch>("createOnRead", false);
+
+    if (dict.found("surfaces"))
+    {
+        fieldSelection_ = wordReListMatcher::uniq
+        (
+            wordReList(dict.lookup("fields"))
+        );
+        Info<< type() << " fields: " << fieldSelection_ << nl;
+
+        if (dict.readIfPresent("derived", derivedNames_))
+        {
+            Info<< type() << " derived: " << derivedNames_ << nl;
+        }
+
+        PtrList<surfMeshSampler> newList
+        (
+            dict.lookup("surfaces"),
+            surfMeshSampler::iNew(mesh_)
+        );
+        transfer(newList);
+
+        // Ensure all surfaces and merge information are expired
+        expire();
+
+        // Need to initialize corresponding surfMesh for others in the chain.
+        // This can simply be a zero-sized placeholder, or the real surface with
+        // faces.
+        if (this->size())
+        {
+            Info<< "Reading surface description:" << nl;
+            forAll(*this, surfI)
+            {
+                surfMeshSampler& s = operator[](surfI);
+
+                Info<< "    " << s.name() << nl;
+                if (createOnRead)
+                {
+                    s.update();
+                }
+                else
+                {
+                    s.create();
+                }
+            }
+            Info<< endl;
+        }
+    }
+
+    return true;
+}
+
+
+void Foam::surfMeshSamplers::updateMesh(const mapPolyMesh& mpm)
+{
+    if (&mpm.mesh() == &mesh_)
+    {
+        expire();
+    }
+
+    // pointMesh and interpolation will have been reset in mesh.update
+}
+
+
+void Foam::surfMeshSamplers::movePoints(const polyMesh& m)
+{
+    if (&m == &mesh_)
+    {
+        expire();
+    }
+}
+
+
+void Foam::surfMeshSamplers::readUpdate(const polyMesh::readUpdateState state)
+{
+    if (state != polyMesh::UNCHANGED)
+    {
+        expire();
+    }
+}
+
+
+bool Foam::surfMeshSamplers::needsUpdate() const
+{
+    forAll(*this, surfI)
+    {
+        if (operator[](surfI).needsUpdate())
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
+
+bool Foam::surfMeshSamplers::expire()
+{
+    bool justExpired = false;
+
+    forAll(*this, surfI)
+    {
+        if (operator[](surfI).expire())
+        {
+            justExpired = true;
+        }
+    }
+
+    // True if any surfaces just expired
+    return justExpired;
+}
+
+
+bool Foam::surfMeshSamplers::update()
+{
+    if (!needsUpdate())
+    {
+        return false;
+    }
+
+    bool updated = false;
+    forAll(*this, surfI)
+    {
+        if (operator[](surfI).update())
+        {
+            updated = true;
+        }
+    }
+
+    return updated;
+}
+
+
+// ************************************************************************* //

src/sampling/surfMeshSampler/surfMeshSamplers/surfMeshSamplers.H

@@ -0,0 +1,228 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 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::surfMeshSamplers
+
+Description
+    Set of surfaces to sample into a surfMesh/surfField.
+
+    The execute() method is used to sample, and the write() method to write.
+    It is fairly common to use for sampling only and have the write disabled.
+
+    \verbatim
+    surfaces
+    {
+        type    surfMeshes;
+        libs    ("libsampling.so");
+
+        // Sample at every time-step
+        executeControl  timeStep;
+        executeInterval   1;
+
+        // Disable writing (or write at same frequency as fields)
+        writeControl    none;
+        writeInterval   1;
+
+        // Fields to be sampled
+        fields          (p U);
+
+        // Optional: pre-defined derived fields to be sampled
+        derived         (rhoU pTotal);
+
+        // Optional: create surface immediately on read
+        // The default is to create a placeholder without any faces.
+        createOnRead    false;
+
+        surfaces
+        (
+            f0surf
+            {
+                type        sampledTriSurfaceMesh;
+                surface     f0surf.obj;
+                source      cells;
+            }
+        );
+    }
+    \endverbatim
+
+
+See also
+    Foam::sampledSurfaces
+
+SourceFiles
+    surfMeshSamplers.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef surfMeshSamplers_H
+#define surfMeshSamplers_H
+
+#include "regionFunctionObject.H"
+#include "surfMeshSampler.H"
+#include "volFieldsFwd.H"
+#include "surfaceFieldsFwd.H"
+#include "wordReList.H"
+#include "IOobjectList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class Time;
+class fvMesh;
+class dictionary;
+
+/*---------------------------------------------------------------------------*\
+                       Class surfMeshSamplers Declaration
+\*---------------------------------------------------------------------------*/
+
+class surfMeshSamplers
+:
+    public functionObjects::regionFunctionObject,
+    public PtrList<surfMeshSampler>
+{
+    // Static data members
+
+        //- output verbosity
+        static bool verbose_;
+
+
+    // Private data
+
+        //- Const reference to fvMesh
+        const fvMesh& mesh_;
+
+
+      // Read from dictonary
+
+        //- Names of fields to sample
+        wordReList fieldSelection_;
+
+        //- Names of derived fields to create and sample
+        wordList derivedNames_;
+
+
+    // Private Member Functions
+
+        //- Remove items by name from objectRegistry
+        static void checkOutNames
+        (
+            const objectRegistry& registry,
+            const UList<word>& names
+        );
+
+
+        //- Filter acceptable fields types
+        template<class Type>
+        wordList acceptType() const;
+
+
+        //- Disallow default bitwise copy construct and assignment
+        surfMeshSamplers(const surfMeshSamplers&) = delete;
+        void operator=(const surfMeshSamplers&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("surfMeshes");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        surfMeshSamplers
+        (
+            const word& name,
+            const Time& time,
+            const dictionary& dict
+        );
+
+        //- Construct for given objectRegistry and dictionary
+        surfMeshSamplers
+        (
+            const word& name,
+            const objectRegistry&,
+            const dictionary&
+        );
+
+
+    //- Destructor
+    virtual ~surfMeshSamplers();
+
+
+    // Member Functions
+
+        //- Do any of the surfaces need an update?
+        virtual bool needsUpdate() const;
+
+        //- Mark the surfaces as needing an update.
+        //  May also free up unneeded data.
+        //  Return false if all surfaces were already marked as expired.
+        virtual bool expire();
+
+        //- Update the surfaces as required and merge surface points (parallel).
+        //  Return false if no surfaces required an update.
+        virtual bool update();
+
+        //- set verbosity level
+        void verbose(const bool verbosity = true);
+
+        //- Read the surfMeshSamplers dictionary
+        virtual bool read(const dictionary&);
+
+        //- Execute, does sampling
+        virtual bool execute();
+
+        //- Write sampled values
+        virtual bool write();
+
+        //- Update for changes of mesh - expires the surfaces
+        virtual void updateMesh(const mapPolyMesh&);
+
+        //- Update for mesh point-motion - expires the surfaces
+        virtual void movePoints(const polyMesh&);
+
+        //- Update for changes of mesh due to readUpdate - expires the surfaces
+        virtual void readUpdate(const polyMesh::readUpdateState state);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "surfMeshSamplersTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/surfMeshSampler/surfMeshSamplers/surfMeshSamplersTemplates.C

@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     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 "surfMeshSamplers.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+Foam::wordList
+Foam::surfMeshSamplers::acceptType() const
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
+
+    return mesh_.names<VolFieldType>(fieldSelection_);
+}
+
+#if 0
+template<class Type>
+Foam::wordList
+Foam::surfMeshSamplers::acceptType
+(
+    const IOobjectList& objects,
+    bool fromFiles
+) const
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
+
+    if (fromFiles_)
+    {
+        // This should actually be in the caller:
+        // IOobjectList objects1 = objects.lookup(fieldSelection_);
+
+        return objects.names(VolFieldType::typeName, fieldSelection_);
+    }
+    else
+    {
+        return mesh_.names<VolFieldType>(fieldSelection_);
+    }
+}
+#endif
+
+
+// ************************************************************************* //

src/sampling/surfMeshSampler/triSurfaceMesh/surfMeshDiscreteSampler.C

@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     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 "surfMeshDiscreteSampler.H"
+#include "MeshedSurfaces.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(surfMeshDiscreteSampler, 0);
+
+    // Add under name "sampledTriSurfaceMesh"
+    // for symmetry with normal sampledSurface
+    addNamedToRunTimeSelectionTable
+    (
+        surfMeshSampler,
+        surfMeshDiscreteSampler,
+        word,
+        sampledTriSurfaceMesh
+    );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::surfMeshDiscreteSampler::transferContent()
+{
+    SurfaceSource& src = static_cast<SurfaceSource&>(*this);
+    surfMesh& dst = getOrCreateSurfMesh();
+
+    dst.transfer(src);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::surfMeshDiscreteSampler::surfMeshDiscreteSampler
+(
+    const word& name,
+    const polyMesh& mesh,
+    const word& surfaceName,
+    const discreteSurface::samplingSource sampleSource
+)
+:
+    surfMeshSampler(name, mesh),
+    SurfaceSource(mesh, surfaceName, sampleSource, false) // no interpolate
+{}
+
+
+Foam::surfMeshDiscreteSampler::surfMeshDiscreteSampler
+(
+    const word& name,
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    surfMeshSampler(name, mesh),
+    SurfaceSource(mesh, dict, false) // no interpolate
+{}
+
+
+Foam::surfMeshDiscreteSampler::surfMeshDiscreteSampler
+(
+    const word& name,
+    const polyMesh& mesh,
+    const triSurface& surface,
+    const word& sampleSourceName
+)
+:
+    surfMeshSampler(name, mesh),
+    SurfaceSource(name, mesh, surface, sampleSourceName, false)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::surfMeshDiscreteSampler::~surfMeshDiscreteSampler()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::surfMeshDiscreteSampler::needsUpdate() const
+{
+    return SurfaceSource::needsUpdate();
+}
+
+
+bool Foam::surfMeshDiscreteSampler::expire()
+{
+    return SurfaceSource::expire();
+}
+
+
+bool Foam::surfMeshDiscreteSampler::update()
+{
+    if (SurfaceSource::update())
+    {
+        transferContent();
+        return true;
+    }
+
+    return false;
+}
+
+
+bool Foam::surfMeshDiscreteSampler::update(const treeBoundBox& bb)
+{
+    if (SurfaceSource::update(bb))
+    {
+        transferContent();
+        return true;
+    }
+
+    return false;
+}
+
+
+bool Foam::surfMeshDiscreteSampler::sample
+(
+    const word& fieldName
+) const
+{
+    return
+    (
+        sampleType<scalar>(fieldName)
+     || sampleType<vector>(fieldName)
+     || sampleType<sphericalTensor>(fieldName)
+     || sampleType<symmTensor>(fieldName)
+     || sampleType<tensor>(fieldName)
+    );
+}
+
+
+// ************************************************************************* //

src/sampling/surfMeshSampler/triSurfaceMesh/surfMeshDiscreteSampler.H

@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     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::surfMeshDiscreteSampler
+
+Description
+    Sampling surfFields onto a surfMesh based on a triSurfaceMesh.
+
+See Also
+    discreteSurface, surfMeshSampler
+
+SourceFiles
+    surfMeshDiscreteSampler.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef surfMeshDiscreteSampler_H
+#define surfMeshDiscreteSampler_H
+
+#include "surfMeshSampler.H"
+#include "discreteSurface.H"
+#include "triSurfaceMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class surfMeshDiscreteSampler;
+
+/*---------------------------------------------------------------------------*\
+                       Class surfMeshDiscreteSampler Declaration
+\*---------------------------------------------------------------------------*/
+
+class surfMeshDiscreteSampler
+:
+    public surfMeshSampler,
+    private discreteSurface
+{
+    // Private typedefs for convenience
+    typedef discreteSurface SurfaceSource;
+
+
+    // Private Member Functions
+
+        //- Transfer mesh content from SurfaceSource to surfMesh
+        void transferContent();
+
+public:
+
+    //- Runtime type information
+    TypeName("surfMeshDiscreteSampler");
+
+
+    // Constructors
+
+        //- Construct from components
+        surfMeshDiscreteSampler
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const word& surfaceName,
+            const discreteSurface::samplingSource sampleSource
+        );
+
+        //- Construct from dictionary
+        surfMeshDiscreteSampler
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Construct from triSurface
+        surfMeshDiscreteSampler
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const triSurface& surface,
+            const word& sampleSourceName
+        );
+
+
+    //- Destructor
+    virtual ~surfMeshDiscreteSampler();
+
+
+    // 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&);
+
+
+        //- Sample the volume field onto surface
+        virtual bool sample(const word& fieldName) const;
+
+        //- Sample field on surface.
+        template<class Type>
+        bool sampleType(const word& fieldName) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "surfMeshDiscreteSamplerTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/surfMeshSampler/triSurfaceMesh/surfMeshDiscreteSamplerTemplates.C

@@ -0,0 +1,97 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     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 "surfMeshDiscreteSampler.H"
+#include "dimensionedType.H"
+#include "error.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+bool Foam::surfMeshDiscreteSampler::sampleType
+(
+    const word& fieldName
+) const
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
+
+    const polyMesh& mesh = SurfaceSource::mesh();
+
+    if (!mesh.foundObject<VolFieldType>(fieldName))
+    {
+        return false;
+    }
+
+    const VolFieldType& fld = mesh.lookupObject<VolFieldType>(fieldName);
+
+    getOrCreateSurfField<Type>(fld).field()
+        = SurfaceSource::sampleField(fld);
+
+    return true;
+}
+
+
+// template<class Type>
+// Foam::tmp<Foam::DimensionedField<Type, Foam::surfGeoMesh>>
+// Foam::surfMeshDiscreteSampler::sampleField
+// (
+//     const GeometricField<Type, fvPatchField, volMesh>& fld
+// ) const
+// {
+//     typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
+//
+//     tmp<Field<Type>> tfield = SurfaceSource::sampleField(fld);
+//     SurfFieldType& result = getOrCreateSurfField<Type>(fld);
+//
+//     // need to verify if this will be needed (in the future)
+//     const surfMesh& s = surface();
+//     if (result.size() != s.size())
+//     {
+//         // maybe resampling changed the surfMesh,
+//         // but the existing surfField wasn't updated
+//
+//         result.setSize(s.size(), Zero);
+//     }
+//
+//     if (result.size() != sampleElements().size())
+//     {
+//         FatalErrorInFunction
+//             << "mismatch in field/mesh sizes "
+//             << result.name() << nl
+//             << " field has " << result.size()
+//             << " but sampleElements has " << sampleElements().size() << nl
+//             << endl
+//             << exit(FatalError);
+//
+//         Info<< "WARNING: "
+//             << endl;
+//     }
+//
+//     result = tfield;
+//     return result;
+// }
+
+
+// ************************************************************************* //

src/sampling/surface/triSurfaceMesh/discreteSurface.C

@@ -0,0 +1,868 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     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 "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 * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(discreteSurface, 0);
+
+    template<>
+    const char* NamedEnum<discreteSurface::samplingSource, 3>::names[] =
+    {
+        "cells",
+        "insideCells",
+        "boundaryFaces"
+    };
+
+    const NamedEnum<discreteSurface::samplingSource, 3>
+    discreteSurface::samplingSourceNames_;
+
+
+    //- 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;
+    forAll(zoneLst, zonei)
+    {
+        const surfZone& zn = zoneLst[zonei];
+        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(mesh().nFaces()-mesh().nInternalFaces());
+        label bndI = 0;
+        forAll(patches, patchi)
+        {
+            const polyPatch& pp = patches[patchi];
+            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()
+                  ? Foam::name("patch%d", 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();
+
+                Map<label>::iterator fnd = zoneSizes.find(regionid);
+                if (fnd != zoneSizes.end())
+                {
+                    fnd()++;
+                }
+                else
+                {
+                    // This shouldn't happen
+                    zoneSizes.insert(regionid, 1);
+                    zoneNames.set
+                    (
+                        regionid,
+                        Foam::name("patch%d", 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;
+    forAllIter(Map<label>, zoneSizes, iter)
+    {
+        // No negative regionids, so Map<label> sorts properly
+        const label regionid = iter.key();
+
+        word name;
+        Map<word>::const_iterator fnd = zoneNames.find(regionid);
+        if (fnd != zoneNames.end())
+        {
+            name = fnd();
+        }
+        if (name.empty())
+        {
+            name = ::Foam::name("patch%d", 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();
+        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 = UIndirectList<label>(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,                   // registry
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        )
+    ),
+    sampleSource_(sampleSource),
+    needsUpdate_(true),
+    keepIds_(false),
+    originalIds_(),
+    zoneIds_(),
+    sampleElements_(0),
+    samplePoints_(0)
+{}
+
+
+Foam::discreteSurface::discreteSurface
+(
+    const polyMesh& mesh,
+    const dictionary& dict,
+    const bool allowInterpolate
+)
+:
+    MeshStorage(),
+    mesh_(mesh),
+    allowInterpolate_(allowInterpolate),
+    interpolate_
+    (
+        allowInterpolate
+     && dict.lookupOrDefault<Switch>("interpolate", false)
+    ),
+    surface_
+    (
+        IOobject
+        (
+            dict.lookup("surface"),
+            mesh.time().constant(), // instance
+            "triSurface",           // local
+            mesh,                   // registry
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        )
+    ),
+    sampleSource_(samplingSourceNames_[dict.lookup("source")]),
+    needsUpdate_(true),
+    keepIds_(dict.lookupOrDefault<Switch>("keepIds", false)),
+    originalIds_(),
+    zoneIds_(),
+    sampleElements_(0),
+    samplePoints_(0)
+{}
+
+
+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,                   // registry
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        surface
+    ),
+    sampleSource_(samplingSourceNames_[sampleSourceName]),
+    needsUpdate_(true),
+    keepIds_(false),
+    originalIds_(),
+    zoneIds_(),
+    sampleElements_(0),
+    samplePoints_(0)
+{}
+
+
+// * * * * * * * * * * * * * * * * 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
+{
+    return
+    (
+        sampleType<scalar>(store, fieldName)
+     || sampleType<vector>(store, fieldName)
+     || sampleType<sphericalTensor>(store, fieldName)
+     || sampleType<symmTensor>(store, fieldName)
+     || sampleType<tensor>(store, fieldName)
+    );
+}
+
+
+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();
+}
+
+
+// ************************************************************************* //

src/sampling/surface/triSurfaceMesh/discreteSurface.H

@@ -0,0 +1,341 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     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/>.
+
+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 NamedEnum<samplingSource, 3> 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);
+
+
+        //- 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;
+
+
+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&);
+
+
+        //- Const access to per-face zone/region information
+        virtual const labelList& zoneIds() const
+        {
+            return zoneIds_;
+        }
+
+        //- If element ids/order of the original surface are kept
+        bool keepIds() const
+        {
+            return keepIds_;
+        }
+
+        //- List of element ids/order of the original surface,
+        //  when keepIds is active.
+        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;
+
+
+        //- Sample field on surface. The tmp field is 'valid' on success.
+        template<class Type>
+        tmp<Field<Type>> sampleType
+        (
+            const word& fieldName
+        ) const;
+
+
+        //- Sample field on surface
+        template<class Type>
+        tmp<Field<Type>> sampleField
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+
+
+        //- Interpolate field on surface
+        template<class Type>
+        tmp<Field<Type>> interpolateField
+        (
+            const interpolation<Type>&
+        ) const;
+
+
+        //- Write information
+        virtual void print(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "discreteSurfaceTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/surface/triSurfaceMesh/discreteSurfaceTemplates.C

@@ -0,0 +1,231 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     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 "discreteSurface.H"
+#include "surfFields.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+bool Foam::discreteSurface::sampleType
+(
+    const objectRegistry& obr,
+    const word& fieldName
+) const
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
+    typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
+    typedef IOField<Type> TmpFieldType;
+
+    if (!mesh().foundObject<VolFieldType>(fieldName))
+    {
+        return false;
+    }
+
+    const VolFieldType& fld = mesh().lookupObject<VolFieldType>(fieldName);
+    tmp<Field<Type>> tfield = sampleField(fld);
+
+    // The rest could be moved into a separate helper
+    if (isA<surfMesh>(obr))
+    {
+        const surfMesh& surf = dynamicCast<const surfMesh>(obr);
+
+        SurfFieldType* ptr = surf.lookupObjectRefPtr<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>("0", fld.dimensions(), Zero)
+            );
+            ptr->writeOpt() = IOobject::NO_WRITE;
+
+            surf.store(ptr);
+        }
+
+        ptr->field() = tfield;
+    }
+    else
+    {
+        TmpFieldType* ptr = obr.lookupObjectRefPtr<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
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
+
+    if (mesh().foundObject<VolFieldType>(fieldName))
+    {
+        const VolFieldType& fld = mesh().lookupObject<VolFieldType>(fieldName);
+        return sampleField(fld);
+    }
+    else
+    {
+        return tmp<Field<Type>>();
+    }
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::discreteSurface::sampleField
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vField
+) const
+{
+    // One value per face
+    tmp<Field<Type>> tvalues(new Field<Type>(sampleElements_.size()));
+    Field<Type>& values = tvalues.ref();
+
+    if (!onBoundary())
+    {
+        // Sample cells
+
+        forAll(sampleElements_, triI)
+        {
+            values[triI] = vField[sampleElements_[triI]];
+        }
+    }
+    else
+    {
+        // Sample boundary faces
+
+        const polyBoundaryMesh& pbm = mesh().boundaryMesh();
+        const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
+
+        // Create flat boundary field
+
+        Field<Type> bVals(nBnd, Zero);
+
+        forAll(vField.boundaryField(), patchi)
+        {
+            label bFacei = pbm[patchi].start() - mesh().nInternalFaces();
+
+            SubList<Type>
+            (
+                bVals,
+                vField.boundaryField()[patchi].size(),
+                bFacei
+            ) = vField.boundaryField()[patchi];
+        }
+
+        // Sample in flat boundary field
+
+        forAll(sampleElements_, triI)
+        {
+            label facei = sampleElements_[triI];
+            values[triI] = bVals[facei-mesh().nInternalFaces()];
+        }
+    }
+
+    return tvalues;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::discreteSurface::interpolateField
+(
+    const interpolation<Type>& interpolator
+) const
+{
+    // One value per vertex
+    tmp<Field<Type>> tvalues(new Field<Type>(sampleElements_.size()));
+    Field<Type>& 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)
+        {
+            label facei = sampleElements_[pointi];
+
+            values[pointi] = interpolator.interpolate
+            (
+                samplePoints_[pointi],
+                mesh().faceOwner()[facei],
+                facei
+            );
+        }
+    }
+
+    return tvalues;
+}
+
+
+// ************************************************************************* //