ENH: new sampling type "surfaceCut"

- an alternative to distanceSurface (with distance zero) that uses
  a cell cutting approach instead of an iso-surface.

src/sampling/Make/files

@@ -22,6 +22,8 @@ sampledSet/shortestPath/shortestPathSet.C
 surface/cutting/cuttingPlane.C
 surface/cutting/cuttingPlaneCuts.C
 surface/cutting/cuttingPlaneSelection.C
+surface/cutting/cuttingSurface.C
+surface/cutting/cuttingSurfaceCuts.C
 surface/cutting/cuttingSurfaceBase.C
 surface/cutting/cuttingSurfaceBaseSelection.C
 surface/distanceSurface/distanceSurface.C
@@ -43,6 +45,7 @@ sampledSurface/isoSurface/sampledIsoSurface.C
 sampledSurface/isoSurface/sampledIsoSurfaceCell.C
 sampledSurface/distanceSurface/sampledDistanceSurface.C
 sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
+sampledSurface/sampledCuttingSurface/sampledCuttingSurface.C
 sampledSurface/sampledSurface/sampledSurface.C
 sampledSurface/sampledSurfaces/sampledSurfaces.C
 sampledSurface/sampledSurfaces/sampledSurfacesGrouping.C

src/sampling/sampledSurface/sampledCuttingSurface/sampledCuttingSurface.C

@@ -0,0 +1,241 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 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 "sampledCuttingSurface.H"
+#include "dictionary.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(sampledCuttingSurface, 0);
+    addToRunTimeSelectionTable
+    (
+        sampledSurface,
+        sampledCuttingSurface,
+        word
+    );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::bitSet Foam::sampledCuttingSurface::cellSelection(const bool warn) const
+{
+    boundBox meshBounds;
+
+    bitSet cellsToSelect =
+        cuttingSurface::cellSelection
+        (
+            mesh(), bounds_, zoneNames_, meshBounds
+        );
+
+    if (warn)
+    {
+        cuttingSurface::checkOverlap(name(), meshBounds, bounds_);
+    }
+
+    return cellsToSelect;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sampledCuttingSurface::sampledCuttingSurface
+(
+    const polyMesh& mesh,
+    const word& surfaceType,
+    const word& surfaceName,
+    const bool triangulate,
+    const boundBox& bounds
+)
+:
+    sampledSurface(surfaceName, mesh),
+    cuttingSurface(mesh, surfaceType, surfaceName),
+    zoneNames_(),
+    bounds_(bounds),
+    triangulate_(triangulate),
+    needsUpdate_(true)
+{}
+
+
+Foam::sampledCuttingSurface::sampledCuttingSurface
+(
+    const word& defaultSurfaceName,
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    sampledSurface(defaultSurfaceName, mesh, dict),
+    cuttingSurface(defaultSurfaceName, mesh, dict),
+    zoneNames_(),
+    bounds_(dict.lookupOrDefault("bounds", boundBox::invertedBox)),
+    triangulate_(dict.lookupOrDefault("triangulate", true)),
+    needsUpdate_(true)
+{
+    if (!dict.readIfPresent("zones", zoneNames_) && dict.found("zone"))
+    {
+        zoneNames_.resize(1);
+        dict.readEntry("zone", zoneNames_.first());
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::sampledCuttingSurface::needsUpdate() const
+{
+    return needsUpdate_;
+}
+
+
+bool Foam::sampledCuttingSurface::expire()
+{
+    // Already marked as expired
+    if (needsUpdate_)
+    {
+        return false;
+    }
+
+    sampledSurface::clearGeom();
+
+    needsUpdate_ = true;
+    return true;
+}
+
+
+bool Foam::sampledCuttingSurface::update()
+{
+    if (!needsUpdate_)
+    {
+        return false;
+    }
+
+    sampledSurface::clearGeom();
+
+    performCut(mesh(), triangulate_, cellSelection(true));
+
+    if (debug)
+    {
+        cuttingSurface::print(Pout);
+        Pout<< endl;
+    }
+
+    needsUpdate_ = false;
+    return true;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::sampledCuttingSurface::sample
+(
+    const interpolation<scalar>& sampler
+) const
+{
+    return sampleOnFaces(sampler);
+}
+
+
+Foam::tmp<Foam::vectorField> Foam::sampledCuttingSurface::sample
+(
+    const interpolation<vector>& sampler
+) const
+{
+    return sampleOnFaces(sampler);
+}
+
+
+Foam::tmp<Foam::sphericalTensorField> Foam::sampledCuttingSurface::sample
+(
+    const interpolation<sphericalTensor>& sampler
+) const
+{
+    return sampleOnFaces(sampler);
+}
+
+
+Foam::tmp<Foam::symmTensorField> Foam::sampledCuttingSurface::sample
+(
+    const interpolation<symmTensor>& sampler
+) const
+{
+    return sampleOnFaces(sampler);
+}
+
+
+Foam::tmp<Foam::tensorField> Foam::sampledCuttingSurface::sample
+(
+    const interpolation<tensor>& sampler
+) const
+{
+    return sampleOnFaces(sampler);
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::sampledCuttingSurface::interpolate
+(
+    const interpolation<scalar>& interpolator
+) const
+{
+    return sampleOnPoints(interpolator);
+}
+
+
+Foam::tmp<Foam::vectorField> Foam::sampledCuttingSurface::interpolate
+(
+    const interpolation<vector>& interpolator
+) const
+{
+    return sampleOnPoints(interpolator);
+}
+
+Foam::tmp<Foam::sphericalTensorField> Foam::sampledCuttingSurface::interpolate
+(
+    const interpolation<sphericalTensor>& interpolator
+) const
+{
+    return sampleOnPoints(interpolator);
+}
+
+
+Foam::tmp<Foam::symmTensorField> Foam::sampledCuttingSurface::interpolate
+(
+    const interpolation<symmTensor>& interpolator
+) const
+{
+    return sampleOnPoints(interpolator);
+}
+
+
+Foam::tmp<Foam::tensorField> Foam::sampledCuttingSurface::interpolate
+(
+    const interpolation<tensor>& interpolator
+) const
+{
+    return sampleOnPoints(interpolator);
+}
+
+
+// ************************************************************************* //

src/sampling/sampledSurface/sampledCuttingSurface/sampledCuttingSurface.H

@@ -0,0 +1,286 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 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::sampledCuttingSurface
+
+Description
+    A surface define by using an input surface to cut the mesh cells
+
+Usage
+    Dictionary controls:
+    \table
+        Property | Description                             | Required | Default
+        type     | surfaceCut                              | yes      |
+        surfaceType | type of surface                      | yes      |
+        surfaceName | name of surface in \c triSurface/    | no  | dict name
+        triangulate | triangulate faces                    | no       | true
+        bounds   | limit with bounding box                 | no       |
+        zone     | limit to cell zone (name or regex)      | no       |
+        zones    | limit to cell zones (names, regexs)     | no       |
+    \endtable
+
+Note
+    No attempt at resolving degenerate cases.
+    Since the cut faces can be quite ugly, they will often be triangulated.
+
+SourceFiles
+    sampledCuttingSurface.C
+    sampledCuttingSurfaceTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sampledCuttingSurface_H
+#define sampledCuttingSurface_H
+
+#include "sampledSurface.H"
+#include "cuttingSurface.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class sampledCuttingSurface Declaration
+\*---------------------------------------------------------------------------*/
+
+class sampledCuttingSurface
+:
+    public sampledSurface,
+    public cuttingSurface
+{
+    // Private Data
+
+        //- The zone or zones in which cutting is to occur
+        wordRes zoneNames_;
+
+        //- Optional bounding box to trim against
+        const boundBox bounds_;
+
+        //- Triangulate faces or not
+        const bool triangulate_;
+
+        //- Track if the surface needs an update
+        mutable bool needsUpdate_;
+
+
+    // Private Member Functions
+
+        //- Define cell selection from zones and bounding box.
+        //  Optionally check and warn if the bounding box
+        //  does not overlap with the mesh (or submesh)
+        bitSet cellSelection(const bool warn=false) const;
+
+
+        //- Sample volume field onto surface faces
+        template<class Type>
+        tmp<Field<Type>> sampleOnFaces
+        (
+            const interpolation<Type>& sampler
+        ) const;
+
+        //- Interpolate volume field onto surface points
+        template<class Type>
+        tmp<Field<Type>> sampleOnPoints
+        (
+            const interpolation<Type>& interpolator
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("surfaceCut");
+
+
+    // Constructors
+
+        //- Construct from components
+        sampledCuttingSurface
+        (
+            const polyMesh& mesh,
+            const word& surfaceType,
+            const word& surfaceName,
+            const bool triangulate = true,
+            const boundBox& bounds = boundBox::invertedBox
+        );
+
+        //- Construct from dictionary
+        sampledCuttingSurface
+        (
+            const word& defaultSurfaceName,
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~sampledCuttingSurface() = default;
+
+
+    // Member Functions
+
+        //- Does the surface need an update?
+        virtual bool needsUpdate() const;
+
+        //- Mark the surface as needing an update.
+        //  May also free up unneeded data.
+        //  Return false if surface was already marked as expired.
+        virtual bool expire();
+
+        //- Update the surface as required.
+        //  Do nothing (and return false) if no update was needed
+        virtual bool update();
+
+
+        //- Points of surface
+        virtual const pointField& points() const
+        {
+            return cuttingSurface::points();
+        }
+
+        //- Faces of surface
+        virtual const faceList& faces() const
+        {
+            return cuttingSurface::surfFaces();
+        }
+
+        //- Const access to per-face zone/region information
+        //  Could instead return meshCells or cellZoneId of the meshCells.
+        virtual const labelList& zoneIds() const
+        {
+            return Foam::emptyLabelList;
+        }
+
+        //- Face area magnitudes
+        virtual const vectorField& Sf() const
+        {
+            return cuttingSurface::Sf();
+        }
+
+        //- Face area magnitudes
+        virtual const scalarField& magSf() const
+        {
+            return cuttingSurface::magSf();
+        }
+
+        //- Face centres
+        virtual const vectorField& Cf() const
+        {
+            return cuttingSurface::Cf();
+        }
+
+        //- For each face, the original cell in mesh
+        using cuttingSurface::meshCells;
+
+
+    // Sample
+
+        //- Sample volume field onto surface faces
+        virtual tmp<scalarField> sample
+        (
+            const interpolation<scalar>& sampler
+        ) const;
+
+        //- Sample volume field onto surface faces
+        virtual tmp<vectorField> sample
+        (
+            const interpolation<vector>& sampler
+        ) const;
+
+        //- Sample volume field onto surface faces
+        virtual tmp<sphericalTensorField> sample
+        (
+            const interpolation<sphericalTensor>& sampler
+        ) const;
+
+        //- Sample volume field onto surface faces
+        virtual tmp<symmTensorField> sample
+        (
+            const interpolation<symmTensor>& sampler
+        ) const;
+
+        //- Sample volume field onto surface faces
+        virtual tmp<tensorField> sample
+        (
+            const interpolation<tensor>& sampler
+        ) const;
+
+
+    // Interpolate
+
+        //- Interpolate volume field onto surface points
+        virtual tmp<scalarField> interpolate
+        (
+            const interpolation<scalar>& interpolator
+        ) const;
+
+        //- Interpolate volume field onto surface points
+        virtual tmp<vectorField> interpolate
+        (
+            const interpolation<vector>& interpolator
+        ) const;
+
+        //- Interpolate volume field onto surface points
+        virtual tmp<sphericalTensorField> interpolate
+        (
+            const interpolation<sphericalTensor>& interpolator
+        ) const;
+
+        //- Interpolate volume field onto surface points
+        virtual tmp<symmTensorField> interpolate
+        (
+            const interpolation<symmTensor>& interpolator
+        ) const;
+
+        //- Interpolate volume field onto surface points
+        virtual tmp<tensorField> interpolate
+        (
+            const interpolation<tensor>& interpolator
+        ) const;
+
+
+    // Output
+
+//         //- Print information
+//         void print(Ostream& os) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "sampledCuttingSurfaceTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/sampledSurface/sampledCuttingSurface/sampledCuttingSurfaceTemplates.C

@@ -0,0 +1,89 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 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 "sampledCuttingSurface.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::sampledCuttingSurface::sampleOnFaces
+(
+    const interpolation<Type>& sampler
+) const
+{
+    return sampledSurface::sampleOnFaces
+    (
+        sampler,
+        meshCells(),
+        faces(),
+        points()
+    );
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::sampledCuttingSurface::sampleOnPoints
+(
+    const interpolation<Type>& interpolator
+) const
+{
+    // elements to sample
+    const labelList& elements = meshCells();
+
+    // One value per point.
+    // Initialize with Zero to handle missed/degenerate faces
+
+    auto tvalues = tmp<Field<Type>>::New(points().size(), Zero);
+    auto& values = tvalues.ref();
+
+    bitSet pointDone(points().size());
+
+    const faceList& fcs = faces();
+
+    forAll(fcs, facei)
+    {
+        const face& f = fcs[facei];
+        const label celli = elements[facei];
+
+        for (const label pointi : f)
+        {
+            if (pointDone.set(pointi))
+            {
+                values[pointi] = interpolator.interpolate
+                (
+                    points()[pointi],
+                    celli
+                );
+            }
+        }
+    }
+
+    return tvalues;
+}
+
+
+// ************************************************************************* //

src/sampling/surface/cutting/cuttingSurface.C

@@ -0,0 +1,176 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 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 "cuttingSurface.H"
+#include "dictionary.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::cuttingSurface::cuttingSurface
+(
+    const polyMesh& mesh,
+    const word& surfaceType,
+    const word& surfaceName
+)
+:
+    cuttingSurfaceBase(),
+    surfPtr_
+    (
+        searchableSurface::New
+        (
+            surfaceType,
+            IOobject
+            (
+                surfaceName,            // name
+                mesh.time().constant(), // directory
+                "triSurface",           // instance
+                mesh.time(),            // registry
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            ),
+            dictionary()
+        )
+    )
+{}
+
+
+Foam::cuttingSurface::cuttingSurface
+(
+    const word& defaultSurfaceName,
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    cuttingSurfaceBase(),
+    surfPtr_
+    (
+        searchableSurface::New
+        (
+            dict.get<word>("surfaceType"),
+            IOobject
+            (
+                dict.lookupOrDefault("surfaceName", defaultSurfaceName),
+                mesh.time().constant(), // directory
+                "triSurface",           // instance
+                mesh.time(),            // registry
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            ),
+            dict
+        )
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::cuttingSurface::performCut
+(
+    const primitiveMesh& mesh,
+    const bool triangulate,
+    bitSet&& cellIdLabels
+)
+{
+    const fvMesh& fvm = static_cast<const fvMesh&>(mesh);
+
+    MeshStorage::clear();
+    meshCells_.clear();
+
+    // Pre-populate with restriction
+    bitSet cellCuts(std::move(cellIdLabels));
+
+    if (cellCuts.size())
+    {
+        cellCuts.resize(mesh.nCells());
+    }
+
+    scalarField pointDist;
+    calcCellCuts(fvm, pointDist, cellCuts);
+
+
+    // Walk cell cuts to create faces
+
+    // Action #1:
+    // - Orient edge so it points in the positive gradient direction.
+    // - Edge intersection when it spans across point-distance == 0.
+    const auto edgeOrientIntersect =
+        [=](edge& e) -> bool
+        {
+            if (pointDist[e.last()] < pointDist[e.first()])
+            {
+                e.flip();
+            }
+
+            const scalar s0 = pointDist[e.first()];
+            const scalar s1 = pointDist[e.last()];
+
+            if
+            (
+                s0 > ROOTVSMALL         // Edge is all positive
+             || s1 < ROOTVSMALL         // Edge is all negative
+             || Foam::mag(s1-s0) < ROOTVSMALL
+            )
+            {
+                return false;
+            }
+
+            return true;
+        };
+
+    // Action #2:
+    // - The edge intersection alpha for point-distance == 0
+    // Return -1 for error.
+    // This is like the iso-fraction for an iso-surface.
+    const auto edgeAlphaIntersect =
+        [=](const edge& e) -> scalar
+        {
+            const scalar s0 = pointDist[e.first()];
+            const scalar s1 = pointDist[e.last()];
+            const scalar d = s1-s0;
+
+            return Foam::mag(d) < ROOTVSMALL ? -1 : (-s0/d);
+        };
+
+    walkCellCuts
+    (
+        mesh,
+        cellCuts,
+        edgeOrientIntersect,
+        edgeAlphaIntersect,
+        triangulate
+    );
+}
+
+
+void Foam::cuttingSurface::print(Ostream& os) const
+{
+    os  << "  surface:" << surfaceName()
+        << "  faces:" << MeshStorage::surfFaces().size()
+        << "  points:" << MeshStorage::points().size();
+}
+
+
+// ************************************************************************* //

src/sampling/surface/cutting/cuttingSurface.H

@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 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::cuttingSurface
+
+Description
+    Constructs a cutting surface through a mesh.
+
+    No attempt at resolving degenerate cases.
+    Since the cut faces can be quite ugly, they will often be triangulated.
+
+SourceFiles
+    cuttingSurface.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cuttingSurface_H
+#define cuttingSurface_H
+
+#include "searchableSurface.H"
+#include "cuttingSurfaceBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class cuttingSurface Declaration
+\*---------------------------------------------------------------------------*/
+
+class cuttingSurface
+:
+    public cuttingSurfaceBase
+{
+    // Private Data
+
+        //- The surface
+        const autoPtr<searchableSurface> surfPtr_;
+
+
+    // Private Member Functions
+
+        //- Determine the cut cells, possibly restricted to a list of cells
+        //
+        //  \param pointDist [out] For each mesh point, the (signed) distance
+        //      to the surface.
+        //  \param cellCuts [in,out] On input an empty set (ie, no restriction)
+        //      or subsetted cells. On output, the cells cut according to the
+        //      planeSides detection.
+        void calcCellCuts
+        (
+            const fvMesh& fvm,
+            scalarField& pointDist,
+            bitSet& cellCuts
+        );
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Cut mesh, restricted to a list of cells
+        using cuttingSurfaceBase::performCut;
+
+        //- Cut mesh, restricted to a list of cells
+        //  Reclaim memory for cellSelectionMask
+        virtual void performCut
+        (
+            const primitiveMesh& mesh,
+            const bool triangulate,
+            bitSet&& cellIdLabels
+        );
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        cuttingSurface
+        (
+            const polyMesh& mesh,
+            const word& surfaceType,
+            const word& surfaceName
+        );
+
+        //- Construct from dictionary
+        cuttingSurface
+        (
+            const word& defaultSurfaceName,
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~cuttingSurface() = default;
+
+
+    // Member Functions
+
+        //- The name of the underlying searchableSurface
+        const word& surfaceName() const
+        {
+            return surfPtr_->name();
+        }
+
+
+    // Output
+
+        //- Print information
+        void print(Ostream& os) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/surface/cutting/cuttingSurfaceCuts.C

@@ -0,0 +1,173 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 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 "cuttingSurface.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "pointFields.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::cuttingSurface::calcCellCuts
+(
+    const fvMesh& fvm,
+    scalarField& pointDist,
+    bitSet& cellCuts
+)
+{
+    const pointField& cc = fvm.C();
+
+    const faceList& faces = fvm.faces();
+    const cellList& cells = fvm.cells();
+    const pointField& pts = fvm.points();
+
+    const label nCells  = fvm.nCells();
+    const label nPoints = fvm.nPoints();
+
+    // The point distances
+
+    List<pointIndexHit> nearest;
+    surfPtr_().findNearest
+    (
+        pts,
+        scalarField(nPoints, GREAT),
+        nearest
+    );
+
+    vectorField norms;
+    surfPtr_().getNormal(nearest, norms);
+
+    pointDist.resize(nPoints);
+
+    for (label i=0; i < nPoints; ++i)
+    {
+        const point diff(pts[i] - nearest[i].hitPoint());
+
+        // Normal distance
+        pointDist[i] = (diff & norms[i]);
+    }
+
+
+    if (cellCuts.size())
+    {
+        cellCuts.resize(nCells);    // safety
+    }
+    else
+    {
+        cellCuts.resize(nCells, true);
+    }
+
+
+    // Check based on cell distance
+
+    surfPtr_().findNearest
+    (
+        cc,
+        scalarField(cc.size(), GREAT),
+        nearest
+    );
+
+
+    boundBox cellBb;
+    labelHashSet cpts;
+
+    for (const label celli : cellCuts)
+    {
+        cellBb.clear();
+
+        if (fvm.hasCellPoints())
+        {
+            cellBb.add(pts, fvm.cellPoints(celli));
+        }
+        else
+        {
+            // DIY version of cellPoints
+            cpts.clear();
+
+            const cell& cFaces = cells[celli];
+
+            for (const label facei : cFaces)
+            {
+                cpts.insert(faces[facei]);
+            }
+
+            cellBb.add(pts, cpts);
+        }
+
+        if (!cellBb.contains(nearest[celli].hitPoint()))
+        {
+            cellCuts.unset(celli);
+        }
+    }
+
+
+    if (debug)
+    {
+        volScalarField cCuts
+        (
+            IOobject
+            (
+                "cuttingSurface.cellCuts",
+                fvm.time().timeName(),
+                fvm.time(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            fvm,
+            dimensionedScalar(dimless, Zero)
+        );
+
+        for (const label celli : cellCuts)
+        {
+            cCuts[celli] = 1;
+        }
+
+        Pout<< "Writing cellCuts:" << cCuts.objectPath() << endl;
+        cCuts.write();
+
+        pointScalarField pDist
+        (
+            IOobject
+            (
+                "cuttingSurface.pointDistance",
+                fvm.time().timeName(),
+                fvm.time(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            pointMesh::New(fvm),
+            dimensionedScalar(dimLength, Zero)
+        );
+        pDist.primitiveFieldRef() = pointDist;
+
+        Pout<< "Writing point distance:" << pDist.objectPath() << endl;
+        pDist.write();
+    }
+}
+
+
+// ************************************************************************* //

tutorials/compressible/rhoSimpleFoam/squareBend/constant/triSurface/angledPlane1.obj

@@ -0,0 +1,25 @@
+# An angled plane with (1,1,1) normal
+o angledPlane1
+
+# Full extent
+# v -0.15 -0.025  0.025
+# v -0.10 -0.025 -0.025
+# v -0.05 -0.075 -0.025
+# v -0.10 -0.075  0.025
+
+# Embedded, ending at cell
+# v -0.135 -0.0325  0.0175
+# v -0.100 -0.0325 -0.0175
+# v -0.065 -0.0675 -0.0175
+# v -0.100 -0.0675  0.0175
+
+# Embedded, ending mid-cell
+
+v -0.1425 -0.02875  0.02125
+v -0.100  -0.02875 -0.02125
+v -0.0575 -0.07125 -0.02125
+v -0.100  -0.07125  0.02125
+
+f 1 2 4
+f 2 3 4
+# EOF

tutorials/compressible/rhoSimpleFoam/squareBend/system/samplingDebug

@@ -62,6 +62,32 @@ debug
             surfaceName angledPlane.obj;
         }
 
+        angledPlaneCut
+        {
+            type        surfaceCut;
+            triangulate false;
+            surfaceType triSurfaceMesh;
+            surfaceName angledPlane.obj;
+        }
+
+        angledPlane1
+        {
+            type        distanceSurface;
+            distance    0;
+            signed      true;
+            regularise  true;
+            surfaceType triSurfaceMesh;
+            surfaceName angledPlane1.obj;
+        }
+
+        angledPlane1Cut
+        {
+            type        surfaceCut;
+            triangulate false;
+            surfaceType triSurfaceMesh;
+            surfaceName angledPlane1.obj;
+        }
+
         disk1
         {
             ${_disk1}