From a947e9ddcff7f45e18d40038c0098cf22fe8cac6 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Tue, 13 Oct 2020 15:52:45 +0200
Subject: [PATCH] ENH: new sampled faceZones (#1874)

---
 src/sampling/Make/files                       |   1 +
 .../sampledFaceZone/sampledFaceZone.C         | 402 ++++++++++++++++++
 .../sampledFaceZone/sampledFaceZone.H         | 345 +++++++++++++++
 .../sampledFaceZoneTemplates.C                | 153 +++++++
 .../sampledMeshedSurface.C                    |  17 +-
 .../sampledMeshedSurface.H                    |  18 +-
 .../windshieldCondensation/Allrun             |   5 +-
 .../windshieldCondensation/Allrun-parallel    |   5 +-
 .../windshieldCondensation/system/controlDict |  21 +
 .../windshieldCondensation/system/topoSetDict |  22 +
 10 files changed, 968 insertions(+), 21 deletions(-)
 create mode 100644 src/sampling/sampledSurface/sampledFaceZone/sampledFaceZone.C
 create mode 100644 src/sampling/sampledSurface/sampledFaceZone/sampledFaceZone.H
 create mode 100644 src/sampling/sampledSurface/sampledFaceZone/sampledFaceZoneTemplates.C

diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index 2d2796dcca5..5a9670931bb 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -35,6 +35,7 @@ surface/isoSurface/isoSurfaceTopo.C
 surface/thresholdCellFaces/thresholdCellFaces.C
 
 sampledSurface/sampledNone/sampledNone.C
+sampledSurface/sampledFaceZone/sampledFaceZone.C
 sampledSurface/sampledPatch/sampledPatch.C
 sampledSurface/sampledPatchInternalField/sampledPatchInternalField.C
 sampledSurface/sampledPlane/sampledPlane.C
diff --git a/src/sampling/sampledSurface/sampledFaceZone/sampledFaceZone.C b/src/sampling/sampledSurface/sampledFaceZone/sampledFaceZone.C
new file mode 100644
index 00000000000..949fc14141a
--- /dev/null
+++ b/src/sampling/sampledSurface/sampledFaceZone/sampledFaceZone.C
@@ -0,0 +1,402 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "sampledFaceZone.H"
+#include "dictionary.H"
+#include "polyMesh.H"
+#include "polyPatch.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "volPointInterpolation.H"
+#include "uindirectPrimitivePatch.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(sampledFaceZone, 0);
+    addNamedToRunTimeSelectionTable
+    (
+        sampledSurface,
+        sampledFaceZone,
+        word,
+        faceZone
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sampledFaceZone::sampledFaceZone
+(
+    const word& name,
+    const polyMesh& mesh,
+    const UList<wordRe>& zoneNames,
+    const bool triangulate
+)
+:
+    sampledSurface(name, mesh),
+    selectionNames_(zoneNames),
+    triangulate_(triangulate),
+    needsUpdate_(true)
+{}
+
+
+Foam::sampledFaceZone::sampledFaceZone
+(
+    const word& name,
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    sampledSurface(name, mesh, dict),
+    selectionNames_(dict.get<wordRes>("zones")),
+    triangulate_(dict.getOrDefault("triangulate", false)),
+    needsUpdate_(true)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::labelList& Foam::sampledFaceZone::zoneIDs() const
+{
+    if (zoneIds_.empty())
+    {
+        // Zone indices for all matches, already sorted
+        zoneIds_ = mesh().faceZones().indices(selectionNames_);
+    }
+
+    return zoneIds_;
+}
+
+
+bool Foam::sampledFaceZone::needsUpdate() const
+{
+    return needsUpdate_;
+}
+
+
+bool Foam::sampledFaceZone::expire()
+{
+    // already marked as expired
+    if (needsUpdate_)
+    {
+        return false;
+    }
+
+    sampledSurface::clearGeom();
+    Mesh::clear();
+
+    zoneIds_.clear();
+
+    faceId_.clear();
+    facePatchId_.clear();
+
+    needsUpdate_ = true;
+    return true;
+}
+
+
+bool Foam::sampledFaceZone::update()
+{
+    if (!needsUpdate_)
+    {
+        return false;
+    }
+
+    // Total number of faces selected
+    label numFaces = 0;
+    for (const label zonei : zoneIDs())
+    {
+        numFaces += mesh().faceZones()[zonei].size();
+    }
+
+    if (zoneIDs().empty())
+    {
+        WarningInFunction
+            << type() << ' ' << name() << ": "
+            << "    No matching face zone(s): "
+            << flatOutput(selectionNames_)  << nl
+            << "    Known face zones: "
+            << flatOutput(mesh().faceZones().names()) << nl;
+    }
+
+    // Could also check numFaces
+
+    // The mesh face or local patch face and the patch id
+    faceId_.resize(numFaces);
+    facePatchId_.resize(numFaces);
+
+    IndirectList<face> selectedFaces(mesh().faces(), labelList());
+    labelList& meshFaceIds = selectedFaces.addressing();
+    meshFaceIds.resize(numFaces);
+
+    numFaces = 0;
+
+    forAll(zoneIDs(), idx)
+    {
+        const label zonei = zoneIDs()[idx];
+        const faceZone& fZone = mesh().faceZones()[zonei];
+
+        for (const label meshFacei : fZone)
+        {
+            // Internal faces
+            label faceId = meshFacei;
+            label facePatchId = -1;
+
+            // Boundary faces
+            if (!mesh().isInternalFace(meshFacei))
+            {
+                facePatchId = mesh().boundaryMesh().whichPatch(meshFacei);
+                const polyPatch& pp = mesh().boundaryMesh()[facePatchId];
+
+                if (isA<emptyPolyPatch>(pp))
+                {
+                    // Do not sample an empty patch
+                    continue;
+                }
+
+                const auto* procPatch = isA<processorPolyPatch>(pp);
+                if (procPatch && !procPatch->owner())
+                {
+                    // Do not sample neighbour-side, retain owner-side only
+                    continue;
+                }
+
+                const auto* cpp = isA<coupledPolyPatch>(pp);
+                if (cpp)
+                {
+                    faceId = (cpp->owner() ? pp.whichFace(meshFacei) : -1);
+                }
+                else
+                {
+                    faceId = meshFacei - pp.start();
+                }
+            }
+
+            if (faceId >= 0)
+            {
+                faceId_[numFaces] = faceId;
+                facePatchId_[numFaces] = facePatchId;
+                meshFaceIds[numFaces] = meshFacei;
+                ++numFaces;
+            }
+        }
+    }
+
+    // Shrink to size used
+    faceId_.resize(numFaces);
+    facePatchId_.resize(numFaces);
+    meshFaceIds.resize(numFaces);
+
+    uindirectPrimitivePatch zoneFaces(selectedFaces, mesh().points());
+
+    this->storedPoints() = zoneFaces.localPoints();
+    this->storedFaces()  = zoneFaces.localFaces();
+
+    // triangulate - uses remapFaces()
+    if (triangulate_)
+    {
+        Mesh::triangulate();
+    }
+
+    needsUpdate_ = false;
+    return true;
+}
+
+
+// remap action on triangulation
+void Foam::sampledFaceZone::remapFaces(const labelUList& faceMap)
+{
+    if (!faceMap.empty())
+    {
+        Mesh::remapFaces(faceMap);
+        faceId_ = labelList
+        (
+            labelUIndList(faceId_, faceMap)
+        );
+        facePatchId_ = labelList
+        (
+            labelUIndList(facePatchId_, faceMap)
+        );
+    }
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::sampledFaceZone::sample
+(
+    const interpolation<scalar>& sampler
+) const
+{
+    return sampleOnFaces(sampler);
+}
+
+
+Foam::tmp<Foam::vectorField> Foam::sampledFaceZone::sample
+(
+    const interpolation<vector>& sampler
+) const
+{
+    return sampleOnFaces(sampler);
+}
+
+
+Foam::tmp<Foam::sphericalTensorField> Foam::sampledFaceZone::sample
+(
+    const interpolation<sphericalTensor>& sampler
+) const
+{
+    return sampleOnFaces(sampler);
+}
+
+
+Foam::tmp<Foam::symmTensorField> Foam::sampledFaceZone::sample
+(
+    const interpolation<symmTensor>& sampler
+) const
+{
+    return sampleOnFaces(sampler);
+}
+
+
+Foam::tmp<Foam::tensorField> Foam::sampledFaceZone::sample
+(
+    const interpolation<tensor>& sampler
+) const
+{
+    return sampleOnFaces(sampler);
+}
+
+
+bool Foam::sampledFaceZone::withSurfaceFields() const
+{
+    return true;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::sampledFaceZone::sample
+(
+    const surfaceScalarField& sField
+) const
+{
+    return sampleOnFaces(sField);
+}
+
+
+Foam::tmp<Foam::vectorField> Foam::sampledFaceZone::sample
+(
+    const surfaceVectorField& sField
+) const
+{
+    return sampleOnFaces(sField);
+}
+
+
+Foam::tmp<Foam::sphericalTensorField> Foam::sampledFaceZone::sample
+(
+    const surfaceSphericalTensorField& sField
+) const
+{
+    return sampleOnFaces(sField);
+}
+
+
+Foam::tmp<Foam::symmTensorField> Foam::sampledFaceZone::sample
+(
+    const surfaceSymmTensorField& sField
+) const
+{
+    return sampleOnFaces(sField);
+}
+
+
+Foam::tmp<Foam::tensorField> Foam::sampledFaceZone::sample
+(
+    const surfaceTensorField& sField
+) const
+{
+    return sampleOnFaces(sField);
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::sampledFaceZone::interpolate
+(
+    const interpolation<scalar>& interpolator
+) const
+{
+    return sampleOnPoints(interpolator);
+}
+
+
+Foam::tmp<Foam::vectorField> Foam::sampledFaceZone::interpolate
+(
+    const interpolation<vector>& interpolator
+) const
+{
+    return sampleOnPoints(interpolator);
+}
+
+Foam::tmp<Foam::sphericalTensorField>
+Foam::sampledFaceZone::interpolate
+(
+    const interpolation<sphericalTensor>& interpolator
+) const
+{
+    return sampleOnPoints(interpolator);
+}
+
+
+Foam::tmp<Foam::symmTensorField> Foam::sampledFaceZone::interpolate
+(
+    const interpolation<symmTensor>& interpolator
+) const
+{
+    return sampleOnPoints(interpolator);
+}
+
+
+Foam::tmp<Foam::tensorField> Foam::sampledFaceZone::interpolate
+(
+    const interpolation<tensor>& interpolator
+) const
+{
+    return sampleOnPoints(interpolator);
+}
+
+
+void Foam::sampledFaceZone::print(Ostream& os) const
+{
+    os  << "faceZone: " << name() << " :"
+        << "  zones: " << flatOutput(selectionNames_)
+        << "  faces:" << faces().size()
+        << "  points:" << points().size();
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/sampledFaceZone/sampledFaceZone.H b/src/sampling/sampledSurface/sampledFaceZone/sampledFaceZone.H
new file mode 100644
index 00000000000..69b648385e6
--- /dev/null
+++ b/src/sampling/sampledSurface/sampledFaceZone/sampledFaceZone.H
@@ -0,0 +1,345 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::sampledFaceZone
+
+Description
+    A sampledSurface defined by the cell faces corresponding to a threshold
+    value.
+
+    This is often embedded as part of a sampled surfaces function object.
+
+Usage
+    Example of function object partial specification:
+    \verbatim
+    surfaces
+    (
+        surface1
+        {
+            type    faceZones;
+            zones   (zone1 "sides.*");
+        }
+    );
+    \endverbatim
+
+    Where the sub-entries comprise:
+    \table
+        Property | Description                             | Required | Default
+        type     | faceZones                               | yes |
+        zones    | zone selection as word/regex list       | yes |
+        triangulate | triangulate faces                    | no  | false
+    \endtable
+
+SourceFiles
+    sampledFaceZone.C
+    sampledFaceZoneTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sampledFaceZone_H
+#define sampledFaceZone_H
+
+#include "sampledSurface.H"
+#include "MeshedSurfaces.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class sampledFaceZone Declaration
+\*---------------------------------------------------------------------------*/
+
+class sampledFaceZone
+:
+    public meshedSurface,
+    public sampledSurface
+{
+    //- Mesh storage type
+    typedef meshedSurface Mesh;
+
+
+    // Private Data
+
+        //- Selection (word/regex) of face zones
+        const wordRes selectionNames_;
+
+        //- IDs for selected face zones (sorted)
+        mutable labelList zoneIds_;
+
+        //- Triangulated faces or keep faces as is
+        bool triangulate_;
+
+        //- Track if the surface needs an update
+        mutable bool needsUpdate_;
+
+        //- Local list of face IDs
+        labelList faceId_;
+
+        //- Local list of patch ID per face. Is -1 for internal face
+        labelList facePatchId_;
+
+
+    // Private Member Functions
+
+        //- Sample volume/boundary field onto surface faces
+        template<class Type>
+        tmp<Field<Type>> sampleOnFaces
+        (
+            const interpolation<Type>& sampler
+        ) const;
+
+        //- Sample surface field onto surface faces
+        template<class Type>
+        tmp<Field<Type>> sampleOnFaces
+        (
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& sField
+        ) const;
+
+        //- Interpolate volume/boundary field onto surface points
+        template<class Type>
+        tmp<Field<Type>> sampleOnPoints
+        (
+            const interpolation<Type>& interpolator
+        ) const;
+
+
+        //- Re-map action on triangulation or cleanup
+        virtual void remapFaces(const labelUList& faceMap);
+
+
+protected:
+
+        //- The selected face zones (sorted)
+        const labelList& zoneIDs() const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("faceZone");
+
+
+    // Constructors
+
+        //- Construct from components
+        sampledFaceZone
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const UList<wordRe>& zoneNames,
+            const bool triangulate = false
+        );
+
+        //- Construct from dictionary
+        sampledFaceZone
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~sampledFaceZone() = 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 Mesh::points();
+        }
+
+        //- Faces of surface
+        virtual const faceList& faces() const
+        {
+            return Mesh::surfFaces();
+        }
+
+        //- Per-face zone/region information
+        virtual const labelList& zoneIds() const
+        {
+            return labelList::null();
+        }
+
+        //- Face area vectors (normals)
+        virtual const vectorField& Sf() const
+        {
+            return Mesh::Sf();
+        }
+
+        //- Face area magnitudes
+        virtual const scalarField& magSf() const
+        {
+            return Mesh::magSf();
+        }
+
+        //- Face centres
+        virtual const vectorField& Cf() const
+        {
+            return Mesh::Cf();
+        }
+
+
+    // Sample
+
+        //- Sample volume field onto surface faces
+        virtual tmp<scalarField> sample
+        (
+            const interpolation<scalar>& sampler
+        ) const;
+
+        //- Sample volume field onto surface faces
+        virtual tmp<vectorField> sample
+        (
+            const interpolation<vector>& sampler
+        ) const;
+
+        //- Sample volume field onto surface faces
+        virtual tmp<sphericalTensorField> sample
+        (
+            const interpolation<sphericalTensor>& sampler
+        ) const;
+
+        //- Sample volume field onto surface faces
+        virtual tmp<symmTensorField> sample
+        (
+            const interpolation<symmTensor>& sampler
+        ) const;
+
+        //- Sample volume field onto surface faces
+        virtual tmp<tensorField> sample
+        (
+            const interpolation<tensor>& sampler
+        ) const;
+
+
+        //- Can it sample surface-fields?
+        virtual bool withSurfaceFields() const;
+
+
+        //- Sample surface field on face zone
+        virtual tmp<scalarField> sample
+        (
+            const surfaceScalarField&
+        ) const;
+
+        //- Sample surface field on face zone
+        virtual tmp<vectorField> sample
+        (
+            const surfaceVectorField&
+        ) const;
+
+        //- Sample surface field on face zone
+        virtual tmp<sphericalTensorField> sample
+        (
+            const surfaceSphericalTensorField&
+        ) const;
+
+        //- Sample surface field on face zone
+        virtual tmp<symmTensorField> sample
+        (
+            const surfaceSymmTensorField&
+        ) const;
+
+        //- Sample surface field on face zone
+        virtual tmp<tensorField> sample
+        (
+            const surfaceTensorField&
+        ) 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
+        virtual void print(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "sampledFaceZoneTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/sampledFaceZone/sampledFaceZoneTemplates.C b/src/sampling/sampledSurface/sampledFaceZone/sampledFaceZoneTemplates.C
new file mode 100644
index 00000000000..5798324490f
--- /dev/null
+++ b/src/sampling/sampledSurface/sampledFaceZone/sampledFaceZoneTemplates.C
@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "sampledFaceZone.H"
+#include "volFieldsFwd.H"
+#include "pointFields.H"
+#include "volPointInterpolation.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::sampledFaceZone::sampleOnFaces
+(
+    const interpolation<Type>& sampler
+) const
+{
+    const auto& vField = sampler.psi();
+    const labelList& own = mesh().faceOwner();
+    const labelList& nei = mesh().faceNeighbour();
+
+    // One value per face
+    auto tvalues = tmp<Field<Type>>::New(faceId_.size());
+    auto& values = tvalues.ref();
+
+    forAll(faceId_, i)
+    {
+        const label facei = faceId_[i];
+        const label patchi = facePatchId_[i];
+
+        if (patchi != -1)
+        {
+            // Boundary face - face id is the patch-local face id
+            values[i] = vField.boundaryField()[patchi][facei];
+        }
+        else
+        {
+            // Internal face
+            values[i] = 0.5*(vField[own[facei]] + vField[nei[facei]]);
+        }
+    }
+
+    return tvalues;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::sampledFaceZone::sampleOnFaces
+(
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& sField
+) const
+{
+    // One value per face
+    auto tvalues = tmp<Field<Type>>::New(faceId_.size());
+    auto& values = tvalues.ref();
+
+    forAll(faceId_, i)
+    {
+        const label facei = faceId_[i];
+        const label patchi = facePatchId_[i];
+
+        if (patchi != -1)
+        {
+            // Boundary face - face id is the patch-local face id
+            values[i] = sField.boundaryField()[patchi][facei];
+        }
+        else
+        {
+            // Internal face
+            values[i] = sField[facei];
+        }
+    }
+
+    return tvalues;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::sampledFaceZone::sampleOnPoints
+(
+    const interpolation<Type>& interpolator
+) const
+{
+    // One value per vertex
+    auto tvalues = tmp<Field<Type>>::New(points().size(), Zero);
+    auto& values = tvalues.ref();
+
+    const labelList& own = mesh().faceOwner();
+
+    bitSet pointDone(points().size());
+
+    forAll(faces(), i)
+    {
+        const face& f = faces()[i];
+        label facei = faceId_[i];
+        const label patchi = facePatchId_[i];
+
+        if (patchi != -1)
+        {
+            // Boundary face. patch-local face id -> mesh face id
+            const polyPatch& pp = mesh().boundaryMesh()[patchi];
+
+            facei += pp.start();
+        }
+
+
+        const label celli = own[facei];
+
+        for (const label pointi : f)
+        {
+            if (pointDone.set(pointi))
+            {
+                values[pointi] = interpolator.interpolate
+                (
+                    points()[pointi],
+                    celli,
+                    facei
+                );
+            }
+        }
+    }
+
+    return tvalues;
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C
index 9d907d286ab..9830373ee51 100644
--- a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C
+++ b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C
@@ -172,7 +172,7 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
 
     List<nearInfo> nearest(fc.size(), nearInfo(GREAT, labelMax));
 
-    if (sampleSource_ == cells)
+    if (sampleSource_ == samplingSource::cells)
     {
         // Search for nearest cell
 
@@ -191,7 +191,7 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
             }
         }
     }
-    else if (sampleSource_ == insideCells)
+    else if (sampleSource_ == samplingSource::insideCells)
     {
         // Search for cell containing point
 
@@ -212,11 +212,11 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
             }
         }
     }
-    else
+    else  // samplingSource::boundaryFaces
     {
-        // Search for nearest boundaryFace
+        // Search for nearest boundary face
+        // on all non-coupled boundary faces
 
-        //- Search on all non-coupled boundary faces
         const auto& bndTree = meshSearcher.nonCoupledBoundaryTree();
 
         forAll(fc, facei)
@@ -312,7 +312,7 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
         }
 
 
-        if (sampleSource_ == cells)
+        if (sampleSource_ == samplingSource::cells)
         {
             // samplePoints_   : per surface point a location inside the cell
             // sampleElements_ : per surface point the cell
@@ -359,7 +359,7 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
                 }
             }
         }
-        else if (sampleSource_ == insideCells)
+        else if (sampleSource_ == samplingSource::insideCells)
         {
             // samplePoints_   : per surface point a location inside the cell
             // sampleElements_ : per surface point the cell
@@ -374,7 +374,7 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
                 samplePoints_[pointi] = pt;
             }
         }
-        else
+        else  // samplingSource::boundaryFaces
         {
             // samplePoints_   : per surface point a location on the boundary
             // sampleElements_ : per surface point the boundary face containing
@@ -406,6 +406,7 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
         // else:
         //      sampleElements_ : per surface triangle the boundary face
         //      samplePoints_   : n/a
+
         sampleElements_.transfer(cellOrFaceLabels);
         samplePoints_.clear();
     }
diff --git a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.H b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.H
index bc4bdfb094d..bfc4ca4c692 100644
--- a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.H
+++ b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.H
@@ -35,9 +35,9 @@ Description
 
     - 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.
+            finds per triangle centre the \em nearest cell centre and uses its value
+        - source=cells, interpolate=true:
+            finds per triangle centre the \em 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
@@ -46,7 +46,7 @@ Description
         - source=insideCells, interpolate=false:
             finds per triangle centre the cell containing it and uses its value.
             Trims triangles outside mesh.
-        - source=insideCells, interpolate=true
+        - source=insideCells, interpolate=true:
             Per surface point interpolate cell containing it.
 
         - source=boundaryFaces, interpolate=false:
@@ -124,11 +124,11 @@ class sampledMeshedSurface
 {
 public:
         //- Types of sampling regions
-        enum samplingSource
+        enum class samplingSource
         {
-            cells,
-            insideCells,
-            boundaryFaces
+            cells,          //!< Use nearest cell value
+            insideCells,    //!< Face within a cell
+            boundaryFaces   //!< Use boundary face values
         };
 
 private:
@@ -287,7 +287,7 @@ public:
         //- Sampling boundary values instead of cell values
         bool onBoundary() const
         {
-            return sampleSource_ == boundaryFaces;
+            return sampleSource_ == samplingSource::boundaryFaces;
         }
 
 
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/windshieldCondensation/Allrun b/tutorials/heatTransfer/chtMultiRegionFoam/windshieldCondensation/Allrun
index a7e2611540c..940382c4de3 100755
--- a/tutorials/heatTransfer/chtMultiRegionFoam/windshieldCondensation/Allrun
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/windshieldCondensation/Allrun
@@ -16,8 +16,9 @@ runApplication subsetMesh c0 -patch walls -overwrite
 runApplication splitMeshRegions -cellZones -overwrite
 
 # create register face and cell zones
-rm log.topoSet
-runApplication topoSet -region cabin -dict system/topoSetDictRegister
+rm -f log.topoSet.register
+runApplication -s register \
+    topoSet -region cabin -dict system/topoSetDictRegister
 
 # set the initial fields
 restore0Dir
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/windshieldCondensation/Allrun-parallel b/tutorials/heatTransfer/chtMultiRegionFoam/windshieldCondensation/Allrun-parallel
index f4a3fab3a95..597eb33fbc9 100755
--- a/tutorials/heatTransfer/chtMultiRegionFoam/windshieldCondensation/Allrun-parallel
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/windshieldCondensation/Allrun-parallel
@@ -16,8 +16,9 @@ runApplication subsetMesh c0 -patch walls -overwrite
 runApplication splitMeshRegions -cellZones -overwrite
 
 # create register face and cell zones
-rm log.topoSet
-runApplication topoSet -region cabin -dict system/topoSetDictRegister
+rm -f log.topoSet.register
+runApplication -o -s register \
+    topoSet -region cabin -dict system/topoSetDictRegister
 
 # set the initial fields
 restore0Dir
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/windshieldCondensation/system/controlDict b/tutorials/heatTransfer/chtMultiRegionFoam/windshieldCondensation/system/controlDict
index 68777ae5cfb..492ce24a652 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/windshieldCondensation/system/controlDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/windshieldCondensation/system/controlDict
@@ -83,6 +83,27 @@ functions
             ( 60    "<system>/solverControls.60")
         );
     }
+
+    // Verify values for sampling multiple faceZones
+    inletFaces
+    {
+        type    surfaces;
+        libs    (sampling);
+        writeControl    onEnd;
+        region          cabin;
+        surfaceFormat   vtk;
+        fields  (H2O U);
+
+        surfaces
+        {
+            inletFaces
+            {
+                type    faceZone;
+                zones   (inletFaces f1Zone);
+                interpolate false;
+            }
+        }
+    }
 }
 
 
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/windshieldCondensation/system/topoSetDict b/tutorials/heatTransfer/chtMultiRegionFoam/windshieldCondensation/system/topoSetDict
index 75a495fdcc9..5c430ef5dac 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/windshieldCondensation/system/topoSetDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/windshieldCondensation/system/topoSetDict
@@ -44,6 +44,28 @@ actions
         type    cellSet;
         action  invert;
     }
+
+    {
+        name    f0tmp;
+        type    faceSet;
+        action  new;
+        source  patchToFace;
+        patch   inlet;
+    }
+
+    {
+        name    inletFaces;
+        type    faceZoneSet;
+        action  new;
+        source  setToFaceZone;
+        faceSet f0tmp;
+    }
+
+    {
+        name    f0tmp;
+        type    faceSet;
+        action  remove;
+    }
 );
 
 
-- 
GitLab