From 1d398d8cae95e65c5e8273a9db3a8a7a209c584d Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Tue, 8 Sep 2020 14:49:34 +0200
Subject: [PATCH 1/8] ENH: handle abaqus SURFACE input (#1600)

- read surfaces which are defined in terms of solid element sides. Eg,

```
  *ELEMENT, TYPE=C3D4, ELSET=...
  1, ...
  2, ...

  *SURFACE, NAME=Things, TYPE=ELEMENT
  1, S1
  2, S1
```

  The element and side number are encoded as a synthetic face id
  according to

      -(10 * elemId + sideNum)

  but the underlying solid geometry is discarded, since there is no
  reasonable way to pass it through the surface sampling mechanism.
---
 src/fileFormats/abaqus/ABAQUSCore.C           | 191 ++++++++++++++++--
 src/fileFormats/abaqus/ABAQUSCore.H           |  66 +++++-
 .../abaqus/ABAQUSsurfaceFormat.C              |   1 +
 3 files changed, 242 insertions(+), 16 deletions(-)

diff --git a/src/fileFormats/abaqus/ABAQUSCore.C b/src/fileFormats/abaqus/ABAQUSCore.C
index 356bffa3c1..c99f217b03 100644
--- a/src/fileFormats/abaqus/ABAQUSCore.C
+++ b/src/fileFormats/abaqus/ABAQUSCore.C
@@ -30,8 +30,9 @@ License
 #include "ListOps.H"
 #include "stringOps.H"
 #include "UIListStream.H"
-
-#undef Foam_readAbaqusSurface
+#include "cellModel.H"
+#include <algorithm>
+#include <cctype>
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -319,11 +320,10 @@ Foam::fileFormats::ABAQUSCore::readHelper::readPoints
     const label initialCount = points_.size();
 
     char sep; // Comma separator (dummy)
+    string line;
     label id;
     point p;
 
-    string line;
-
     // Read nodes (points) until next "*Section"
     while (is.peek() != '*' && is.peek() != EOF)
     {
@@ -368,10 +368,10 @@ Foam::fileFormats::ABAQUSCore::readHelper::readElements
     const label initialCount = elemTypes_.size();
 
     char sep; // Comma separator (dummy)
-    label id;
-    labelList elemNodes(nNodes, Zero);
-
     string line;
+    label id;
+
+    labelList elemNodes(nNodes, Zero);
 
     // Read element connectivity until next "*Section"
 
@@ -403,6 +403,130 @@ Foam::fileFormats::ABAQUSCore::readHelper::readElements
 }
 
 
+Foam::label
+Foam::fileFormats::ABAQUSCore::readHelper::readSurfaceElements
+(
+    ISstream& is,
+    const label setId
+)
+{
+    // Info<< "*Surface" << nl;
+
+    // Models for supported solids (need to face mapping)
+    const cellModel& tet   = cellModel::ref(cellModel::TET);
+    const cellModel& prism = cellModel::ref(cellModel::PRISM);
+    const cellModel& hex   = cellModel::ref(cellModel::HEX);
+
+    // Face mapping from Abaqus cellModel to OpenFOAM cellModel
+    const auto& abqToFoamFaceMap = abaqusToFoamFaceAddr();
+
+    const label initialCount = elemTypes_.size();
+
+    char sep; // Comma separator (dummy)
+    string line;
+    label id;
+
+    // Read until next "*Section"
+
+    // Parse for elemId, sideId.
+    // Eg, "1235, S1"
+    while (is.peek() != '*' && is.peek() != EOF)
+    {
+        is >> id >> sep;
+        is.getLine(line);
+
+        const word sideName(word::validate(stringOps::upper(line)));
+
+        if
+        (
+            sideName.size() != 2
+         || sideName[0] != 'S'
+         || !std::isdigit(sideName[1])
+        )
+        {
+            Info<< "Abaqus reader: unsupported surface element side "
+                << id << ", " << sideName << nl;
+            continue;
+        }
+
+        const label index = elemIds_.find(id);
+        if (id <= 0 || index < 0)
+        {
+            Info<< "Abaqus reader: unsupported surface element "
+                << id << nl;
+            continue;
+        }
+
+        const auto faceIdIter = abqToFoamFaceMap.cfind(elemTypes_[index]);
+        if (!faceIdIter.found())
+        {
+            Info<< "Abaqus reader: reject non-solid shape: " << nl;
+        }
+
+        // The abaqus element side number (1-based)
+        const label sideNum = (sideName[1] - '0');
+
+        const label foamFaceNum = (*faceIdIter)[sideNum - 1];
+
+        const labelList& connect = connectivity_[index];
+
+        // Nodes for the derived shell element
+        labelList elemNodes;
+
+        switch (elemTypes_[index])
+        {
+            case shapeType::abaqusTet:
+            {
+                elemNodes = labelList(connect, tet.modelFaces()[foamFaceNum]);
+                break;
+            }
+            case shapeType::abaqusPrism:
+            {
+                elemNodes = labelList(connect, prism.modelFaces()[foamFaceNum]);
+                break;
+            }
+            case shapeType::abaqusHex:
+            {
+                elemNodes = labelList(connect, hex.modelFaces()[foamFaceNum]);
+                break;
+            }
+            default:
+                break;
+        }
+
+        enum shapeType shape = shapeType::abaqusUnknownShape;
+
+        if (elemNodes.size() == 3)
+        {
+            shape = shapeType::abaqusTria;
+        }
+        else if (elemNodes.size() == 4)
+        {
+            shape = shapeType::abaqusQuad;
+        }
+        else
+        {
+            // Cannot happen
+            FatalErrorInFunction
+                << "Could not map face side for "
+                << id << ", " << sideName << nl
+                << exit(FatalError);
+        }
+
+        // Synthesize face Id from solid element Id and side Id
+        const label newElemId = ABAQUSCore::encodeSolidId(id, sideNum);
+
+        // Further checks?
+        connectivity_.append(std::move(elemNodes));
+        elemTypes_.append(shape);
+        elemIds_.append(newElemId);
+        elsetIds_.append(setId);
+    }
+
+    return (elemTypes_.size() - initialCount);
+}
+
+
 void Foam::fileFormats::ABAQUSCore::readHelper::read
 (
     ISstream& is
@@ -427,9 +551,11 @@ void Foam::fileFormats::ABAQUSCore::readHelper::read
         // Some abaqus files use upper-case or mixed-case for section names,
         // convert all to upper-case for ease.
 
-        string upperLine(stringOps::upper(line));
+        const string upperLine(stringOps::upper(line));
 
+        //
         // "*Nodes" section
+        //
         if (upperLine.starts_with("*NODE"))
         {
             // Ignore "NSET=...", we cannot do anything useful with it
@@ -446,14 +572,16 @@ void Foam::fileFormats::ABAQUSCore::readHelper::read
             continue;
         }
 
+        //
         // "*Element" section
+        //
         if (upperLine.starts_with("*ELEMENT,"))
         {
             // Must have "TYPE=..."
-            auto elemTypeName = getIdentifier("TYPE", line);
+            const string elemTypeName(getIdentifier("TYPE", line));
 
             // May have "ELSET=..." on the same line
-            string elsetName(getIdentifier("ELSET", line));
+            const string elsetName(getIdentifier("ELSET", line));
 
             const shapeType shape(getElementType(elemTypeName));
 
@@ -485,18 +613,42 @@ void Foam::fileFormats::ABAQUSCore::readHelper::read
             continue;
         }
 
-
+        //
         // "*Surface" section
+        //
         if (upperLine.starts_with("*SURFACE,"))
         {
-            #ifdef Foam_readAbaqusSurface
+            // Require "NAME=..." on the same line
+            const string elsetName(getIdentifier("NAME", line));
+
+            // May have "TYPE=..." on the same line.
+            // If missing, default is ELEMENT.
+            const string surfTypeName(getIdentifier("TYPE", line));
+
+            if
+            (
+                !surfTypeName.empty()
+             && stringOps::upper(surfTypeName) != "ELEMENT"
+            )
+            {
+                Info<< "Reading abaqus surface type "
+                    << surfTypeName << " is not implemented" << nl;
+                continue;
+            }
+
+            // Treat like an element set
+            const label elsetId = addNewElset(elsetName);
 
             skipComments(is);
 
-            #else
-            Info<< "Reading of abaqus surfaces not implemented" << nl;
-            #endif
+            nread = readSurfaceElements(is, elsetId);
 
+            if (verbose_)
+            {
+                InfoErr
+                    << "Read " << nread << " *SURFACE entries for "
+                    << elsetName << nl;
+            }
             continue;
         }
     }
@@ -623,6 +775,15 @@ void Foam::fileFormats::ABAQUSCore::readHelper::compact_nodes()
 }
 
 
+void Foam::fileFormats::ABAQUSCore::readHelper::renumber_elements_1to0()
+{
+    for (label& elemId : elemIds_)
+    {
+        renumber0_elemId(elemId);
+    }
+}
+
+
 void Foam::fileFormats::ABAQUSCore::writePoints
 (
     Ostream& os,
diff --git a/src/fileFormats/abaqus/ABAQUSCore.H b/src/fileFormats/abaqus/ABAQUSCore.H
index 3f99516c85..4149d59621 100644
--- a/src/fileFormats/abaqus/ABAQUSCore.H
+++ b/src/fileFormats/abaqus/ABAQUSCore.H
@@ -90,7 +90,6 @@ SourceFiles
 
 namespace Foam
 {
-
 namespace fileFormats
 {
 
@@ -148,11 +147,62 @@ public:
             return (tag & 0x80);
         }
 
+        //- Is a combined (synthetic) face Id?
+        inline static bool isEncodedSolidId(const label combinedId)
+        {
+            return (combinedId < 0);
+        }
+
+        //- Combine solid element Id and side Id into synthetic face Id
+        inline static label encodeSolidId(const label id, const label side)
+        {
+            return -(10 * id + side);
+        }
+
+        //- Entangle solid element id from synthetic face Id
+        inline static label decodeSolidElementId(const label combinedId)
+        {
+            return
+            (
+                isEncodedSolidId(combinedId)
+              ? (-combinedId / 10)
+              : combinedId
+            );
+        }
+
+        //- Entangle solid side id from synthetic face Id
+        //- Synthesize faceId from solid element Id and sideId
+        inline static label decodeSolidSideNum(const label combinedId)
+        {
+            return
+            (
+                isEncodedSolidId(combinedId)
+              ? (-combinedId % 10)
+              : 0
+            );
+        }
+
 
 protected:
 
     // Protected Member Functions
 
+        //- From 1-based to 0-based
+        inline static void renumber0_elemId(label& elemId)
+        {
+            if (isEncodedSolidId(elemId))
+            {
+                // Eg,
+                // 1-based (solid 4, side 2) is -42
+                // 0-based (solid 3, side 2) is -32
+                elemId += 10;
+            }
+            else
+            {
+                --elemId;
+            }
+        }
+
         //- Face addressing from ABAQUS faces to OpenFOAM faces.
         //  For hex, prism, tet primitive shapes.
         static const Map<labelList>& abaqusToFoamFaceAddr();
@@ -249,12 +299,26 @@ protected:
                     const label setId = 0
                 );
 
+                //- Read elements within an "*Surface" section.
+                //  If the shape is known/supported, appends to
+                //  connectivity, elemType, elemIds lists.
+                //
+                //  \return the number of elements read
+                label readSurfaceElements
+                (
+                    ISstream& is,
+                    const label setId = 0
+                );
+
 
                 //- Remove non-shell elements and compact the points
                 void purge_solids();
 
                 //- Compact unused points and relabel connectivity
                 void compact_nodes();
+
+                //- Renumber elements from 1-based to 0-based
+                void renumber_elements_1to0();
         };
 
 
diff --git a/src/surfMesh/surfaceFormats/abaqus/ABAQUSsurfaceFormat.C b/src/surfMesh/surfaceFormats/abaqus/ABAQUSsurfaceFormat.C
index 72ebcfc65e..cfcd5357dc 100644
--- a/src/surfMesh/surfaceFormats/abaqus/ABAQUSsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/abaqus/ABAQUSsurfaceFormat.C
@@ -154,6 +154,7 @@ bool Foam::fileFormats::ABAQUSsurfaceFormat<Face>::read
     // No solids
     reader.purge_solids();
     reader.compact_nodes();
+    reader.renumber_elements_1to0();  // Convert 1-based -> 0-based
 
 
     // Convert connectivity to faces

From f82af7cb2abdba8175ab5e270f76525bf761e6ad Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Mon, 14 Sep 2020 16:39:36 +0200
Subject: [PATCH 2/8] ENH: initial support for abaqus surface sampled output
 (#1600)

- supports geometry and field-specific scaling, separate geometry and
  fields. Beta-feature for suppressing geometry output entirely.
---
 src/surfMesh/Make/files                       |   1 +
 .../writers/abaqus/abaqusSurfaceWriter.C      | 355 ++++++++++++++++++
 .../writers/abaqus/abaqusSurfaceWriter.H      | 239 ++++++++++++
 .../writers/abaqus/abaqusSurfaceWriterImpl.C  | 308 +++++++++++++++
 4 files changed, 903 insertions(+)
 create mode 100644 src/surfMesh/writers/abaqus/abaqusSurfaceWriter.C
 create mode 100644 src/surfMesh/writers/abaqus/abaqusSurfaceWriter.H
 create mode 100644 src/surfMesh/writers/abaqus/abaqusSurfaceWriterImpl.C

diff --git a/src/surfMesh/Make/files b/src/surfMesh/Make/files
index 7c68022eb8..eac20cc94a 100644
--- a/src/surfMesh/Make/files
+++ b/src/surfMesh/Make/files
@@ -63,6 +63,7 @@ writers = writers
 
 $(writers)/surfaceWriter.C
 $(writers)/caching/surfaceWriterCaching.C
+$(writers)/abaqus/abaqusSurfaceWriter.C
 $(writers)/boundaryData/boundaryDataSurfaceWriter.C
 $(writers)/ensight/ensightSurfaceWriter.C
 $(writers)/foam/foamSurfaceWriter.C
diff --git a/src/surfMesh/writers/abaqus/abaqusSurfaceWriter.C b/src/surfMesh/writers/abaqus/abaqusSurfaceWriter.C
new file mode 100644
index 0000000000..8932d90684
--- /dev/null
+++ b/src/surfMesh/writers/abaqus/abaqusSurfaceWriter.C
@@ -0,0 +1,355 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "abaqusSurfaceWriter.H"
+#include "ABAQUSCore.H"
+#include "IOmanip.H"
+#include "ListOps.H"
+#include "OSspecific.H"
+#include "surfaceWriterMethods.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace surfaceWriters
+{
+    defineTypeName(abaqusWriter);
+    addToRunTimeSelectionTable(surfaceWriter, abaqusWriter, word);
+    addToRunTimeSelectionTable(surfaceWriter, abaqusWriter, wordDict);
+}
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Field writing implementation
+#include "abaqusSurfaceWriterImpl.C"
+
+// Field writing methods
+defineSurfaceWriterWriteFields(Foam::surfaceWriters::abaqusWriter);
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Write connectivity as CSV list
+inline static void writeConnectivity
+(
+    Ostream& os,
+    const label elemId,
+    const labelUList& elem
+)
+{
+    os  << "  " << elemId;
+
+    for (const label vert : elem)
+    {
+        os << ", " << (vert + 1);
+    }
+
+    os << nl;
+}
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::surfaceWriters::abaqusWriter::writeFace
+(
+    Ostream& os,
+    const labelUList& f,
+    const label elemId,
+    const label propId,
+    bool header
+) const
+{
+    // Only called with 3 or 4 points!
+
+    if (header)
+    {
+        os  << "*ELEMENT, TYPE=S" << f.size();
+
+        if (propId >= 0)
+        {
+            os  << ", ELSET=_" << propId;
+        }
+
+        os  << nl;
+    }
+
+    writeConnectivity(os, elemId, f);
+}
+
+
+void Foam::surfaceWriters::abaqusWriter::writeGeometry
+(
+    Ostream& os,
+    const meshedSurf& surf,
+    labelList& decompOffsets,
+    DynamicList<face>& decompFaces
+) const
+{
+    const pointField& points = surf.points();
+    const faceList&    faces = surf.faces();
+    const labelList&   zones = surf.zoneIds();
+    const labelList& elemIds = surf.faceIds();
+
+    // Possible to use faceIds?
+    bool useOrigFaceIds =
+    (
+        elemIds.size() == faces.size()
+     && !ListOps::found(elemIds, lessOp1<label>(0))
+    );
+
+    if (useOrigFaceIds)
+    {
+        // Not possible with on-the-fly face decomposition
+        for (const auto& f : faces)
+        {
+            if (f.size() > 4)
+            {
+                useOrigFaceIds = false;
+                break;
+            }
+        }
+    }
+
+
+    os  << "** Geometry" << nl;
+
+    os  << nl
+        << "**" << nl
+        << "** Points" << nl
+        << "**" << nl;
+
+    fileFormats::ABAQUSCore::writePoints(os, points, geometryScale_);
+
+
+    // Write faces, with on-the-fly decomposition (triangulation)
+    decompOffsets.resize(faces.size()+1);
+    decompFaces.clear();
+
+    decompOffsets[0] = 0; // The first offset is always zero
+
+    os  << "**" << nl
+        << "** Faces" << nl
+        << "**" << nl;
+
+    // Simple tracking for change of element type/set
+    labelPair prevOutput(-1, -1);
+
+    label elemId = 0;  // The element-id
+    forAll(faces, facei)
+    {
+        const face& f = faces[facei];
+
+        if (useOrigFaceIds)
+        {
+            // When available and not decomposed
+            elemId = elemIds[facei];
+        }
+
+        // 1-offset for PID
+        const label propId = 1 + (facei < zones.size() ? zones[facei] : 0);
+
+        const label n = f.size();
+
+        bool header =
+            (prevOutput.first() != n || prevOutput.second() != propId);
+
+        if (header)
+        {
+            // Update values
+            prevOutput.first() = n;
+            prevOutput.second() = propId;
+        }
+
+        if (n == 3 || n == 4)
+        {
+            writeFace(os, f, ++elemId, propId, header);
+        }
+        else
+        {
+            // Decompose into tris
+            prevOutput.first() = 3;
+
+            f.triangles(points, decompFaces);
+
+            for
+            (
+                label decompi = decompOffsets[facei];
+                decompi < decompFaces.size();
+                ++decompi
+            )
+            {
+                writeFace
+                (
+                    os,
+                    decompFaces[decompi],
+                    ++elemId,
+                    propId,
+                    header
+                );
+
+                header = false;
+            }
+        }
+
+        // The end offset, which is the next begin offset
+        decompOffsets[facei+1] = decompFaces.size();
+    }
+
+    os  << "**" << nl
+        << "**" << nl;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::surfaceWriters::abaqusWriter::abaqusWriter()
+:
+    surfaceWriter(),
+    geometryScale_(1),
+    fieldScale_(),
+    noGeometry_(false),
+    outputLayout_(outputLayoutType::BY_FIELD)
+{}
+
+
+Foam::surfaceWriters::abaqusWriter::abaqusWriter
+(
+    const dictionary& options
+)
+:
+    surfaceWriter(options),
+    geometryScale_(options.getOrDefault<scalar>("scale", 1)),
+    fieldScale_(options.subOrEmptyDict("fieldScale")),
+    noGeometry_(options.getOrDefault("noGeometry", false)),
+    outputLayout_(outputLayoutType::BY_FIELD)
+{}
+
+
+Foam::surfaceWriters::abaqusWriter::abaqusWriter
+(
+    const meshedSurf& surf,
+    const fileName& outputPath,
+    bool parallel,
+    const dictionary& options
+)
+:
+    abaqusWriter(options)
+{
+    open(surf, outputPath, parallel);
+}
+
+
+Foam::surfaceWriters::abaqusWriter::abaqusWriter
+(
+    const pointField& points,
+    const faceList& faces,
+    const fileName& outputPath,
+    bool parallel,
+    const dictionary& options
+)
+:
+    abaqusWriter(options)
+{
+    open(points, faces, outputPath, parallel);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::fileName Foam::surfaceWriters::abaqusWriter::write()
+{
+    checkOpen();
+
+    // Geometry:
+    // 1) rootdir/<TIME>/surfaceName.abq
+    // 2) rootdir/geometry/surfaceName_<TIME>.abq
+
+    fileName outputFile;
+
+    switch (outputLayout_)
+    {
+        case outputLayoutType::BY_TIME:
+        {
+            outputFile = outputPath_;
+            if (useTimeDir() && !timeName().empty())
+            {
+                // Splice in time-directory
+                outputFile =
+                    outputPath_.path() / timeName() / outputPath_.name();
+            }
+            break;
+        }
+        case outputLayoutType::BY_FIELD:
+        {
+            outputFile = outputPath_ / "geometry" / outputPath_.name();
+            if (!timeName().empty())
+            {
+                // Append time information to file name
+                outputFile += '_' + timeName();
+            }
+            break;
+        }
+    }
+    outputFile.ext("abq");
+
+    if (verbose_)
+    {
+        Info<< "Writing abaqus geometry to " << outputFile << endl;
+    }
+
+
+    const meshedSurf& surf = surface();
+
+    if (Pstream::master() || !parallel_)
+    {
+        if (!isDir(outputFile.path()))
+        {
+            mkDir(outputFile.path());
+        }
+
+        OFstream os(outputFile);
+
+        labelList decompOffsets;
+        DynamicList<face> decompFaces;
+
+        writeGeometry(os, surf, decompOffsets, decompFaces);
+    }
+
+    wroteGeom_ = true;
+    return outputFile;
+}
+
+
+// ************************************************************************* //
diff --git a/src/surfMesh/writers/abaqus/abaqusSurfaceWriter.H b/src/surfMesh/writers/abaqus/abaqusSurfaceWriter.H
new file mode 100644
index 0000000000..d740986298
--- /dev/null
+++ b/src/surfMesh/writers/abaqus/abaqusSurfaceWriter.H
@@ -0,0 +1,239 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::surfaceWriters::abaqusWriter
+
+Description
+    A surface writer for the ABAQUS file format - both surface mesh and fields
+
+    The formatOptions for abaqus:
+    \table
+        Property   | Description                            | Required | Default
+        scale      | output geometry scaling                | no  | 1
+        fieldScale | output field scaling (dictionary)      | no  | empty
+        noGeometry | Suppress geometry output (beta feature) | no  | false
+    \endtable
+
+    For example,
+    \verbatim
+    formatOptions
+    {
+        abaqus
+        {
+            scale   1000;     // [m] -> [mm]
+            fieldScale
+            {
+               "p.*"   0.01;  // [Pa] -> [mbar]
+            }
+        }
+    }
+    \endverbatim
+
+    \heading Output file locations
+
+    The \c rootdir normally corresponds to something like
+    \c postProcessing/\<name\>
+
+    \subheading Geometry
+    \verbatim
+    rootdir
+    `-- surfaceName0
+        `-- geometry
+            `-- surfaceName0_<time>.abq
+    \endverbatim
+
+    \subheading Fields
+    \verbatim
+    rootdir
+    `-- surfaceName0
+        `-- field0
+        |   `-- surfaceName0_<time>.{inp}
+        `-- field1
+            `-- surfaceName0_<time>.{inp}
+    \endverbatim
+
+Note
+    Output variable scaling does not apply to integer types such as Ids.
+
+SourceFiles
+    abaqusSurfaceWriter.C
+    abaqusSurfaceWriterImpl.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef abaqusSurfaceWriter_H
+#define abaqusSurfaceWriter_H
+
+#include "surfaceWriter.H"
+#include "ABAQUSCore.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace surfaceWriters
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class abaqusWriter Declaration
+\*---------------------------------------------------------------------------*/
+
+class abaqusWriter
+:
+    public surfaceWriter
+{
+    // Private Data
+
+        //- Output geometry scaling
+        const scalar geometryScale_;
+
+        //- Output field scaling
+        const dictionary fieldScale_;
+
+        //- BETA feature
+        bool noGeometry_;
+
+        //- Output directory layouts
+        enum class outputLayoutType
+        {
+            BY_TIME = 0,
+            BY_FIELD,
+        };
+
+        //- Output directory layout
+        outputLayoutType outputLayout_;
+
+
+    // Private Member Functions
+
+        //- Write face element (tri/quad)
+        void writeFace
+        (
+            Ostream& os,
+            const labelUList& f,
+            const label elemId,     //!< 1-based Element Id
+            const label propId,     //!< 1-based Property Id
+            const bool header = true
+        ) const;
+
+        //- Write the surface mesh geometry, tracking face decomposition
+        //
+        //  \param decompOffsets  begin/end offsets (size+1) into decompFaces
+        //  \param decompFaces    Non tri/quad decomposed into triangles
+        void writeGeometry
+        (
+            Ostream& os,
+            const meshedSurf& surf,
+            labelList& decompOffsets,
+            DynamicList<face>& decompFaces
+        ) const;
+
+        //- Write a face-based value
+        template<class Type>
+        Ostream& writeFaceValue
+        (
+            Ostream& os,
+            const Type& value,
+            const label elemId      //!< 0-based Element Id
+        ) const;
+
+        //- Templated write operation
+        template<class Type>
+        fileName writeTemplate
+        (
+            const word& fieldName,          //!< Name of field
+            const Field<Type>& localValues  //!< Local field values to write
+        );
+
+
+public:
+
+    //- Declare type-name, virtual type (with debug switch)
+    TypeNameNoDebug("abaqus");
+
+
+    // Constructors
+
+        //- Default construct
+        abaqusWriter();
+
+        //- Construct with some output options
+        explicit abaqusWriter(const dictionary& options);
+
+        //- Construct from components
+        abaqusWriter
+        (
+            const meshedSurf& surf,
+            const fileName& outputPath,
+            bool parallel = Pstream::parRun(),
+            const dictionary& options = dictionary()
+        );
+
+        //- Construct from components
+        abaqusWriter
+        (
+            const pointField& points,
+            const faceList& faces,
+            const fileName& outputPath,
+            bool parallel = Pstream::parRun(),
+            const dictionary& options = dictionary()
+        );
+
+
+    //- Destructor
+    virtual ~abaqusWriter() = default;
+
+
+    // Member Functions
+
+        //- Format uses faceIds as part of its output
+        virtual bool usesFaceIds() const // override
+        {
+            return true;
+        }
+
+        //- Write surface geometry to file.
+        virtual fileName write(); // override
+
+        declareSurfaceWriterWriteMethod(label);
+        declareSurfaceWriterWriteMethod(scalar);
+        declareSurfaceWriterWriteMethod(vector);
+        declareSurfaceWriterWriteMethod(sphericalTensor);
+        declareSurfaceWriterWriteMethod(symmTensor);
+        declareSurfaceWriterWriteMethod(tensor);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceWriters
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/surfMesh/writers/abaqus/abaqusSurfaceWriterImpl.C b/src/surfMesh/writers/abaqus/abaqusSurfaceWriterImpl.C
new file mode 100644
index 0000000000..b8d17ed306
--- /dev/null
+++ b/src/surfMesh/writers/abaqus/abaqusSurfaceWriterImpl.C
@@ -0,0 +1,308 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "OFstream.H"
+#include "IOmanip.H"
+#include "OSspecific.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+Foam::Ostream& Foam::surfaceWriters::abaqusWriter::writeFaceValue
+(
+    Ostream& os,
+    const Type& value,
+    const label elemId
+) const
+{
+    if (fileFormats::ABAQUSCore::isEncodedSolidId(elemId))
+    {
+        const label solidId =
+            fileFormats::ABAQUSCore::decodeSolidElementId(elemId);
+
+        const label sideNum =
+            fileFormats::ABAQUSCore::decodeSolidSideNum(elemId);
+
+        // Element-id: 0-based to 1-based
+        // Side-id: always 1-based
+        os  << (solidId + 1) << ", P" << sideNum;
+    }
+    else
+    {
+        // Element-id: 0-based to 1-based
+        os  << (elemId + 1) << ", P";
+    }
+
+    os << ", ";
+    if (pTraits<Type>::nComponents == 1)
+    {
+        os << value;
+    }
+    else
+    {
+        os << mag(value);
+    }
+    os  << nl;
+
+    return os;
+}
+
+
+template<class Type>
+Foam::fileName Foam::surfaceWriters::abaqusWriter::writeTemplate
+(
+    const word& fieldName,
+    const Field<Type>& localValues
+)
+{
+    checkOpen();
+
+    // Field:
+    // 1) rootdir/<TIME>/<field>_surfaceName.inp
+    // 2) rootdir/<field>/surfaceName_<TIME>.inp
+
+    fileName outputFile;
+
+    switch (outputLayout_)
+    {
+        case outputLayoutType::BY_TIME:
+        {
+            outputFile = outputPath_;
+            if (useTimeDir() && !timeName().empty())
+            {
+                // Splice in time-directory
+                outputFile =
+                    outputPath_.path() / timeName() / outputPath_.name();
+            }
+
+            // Append <field>_surfaceName
+            outputFile /= fieldName + '_' + outputPath_.name();
+            break;
+        }
+        case outputLayoutType::BY_FIELD:
+        {
+            outputFile = outputPath_ / fieldName / outputPath_.name();
+            if (!timeName().empty())
+            {
+                // Append time information to file name
+                outputFile += '_' + timeName();
+            }
+            break;
+        }
+    }
+    outputFile.ext("inp");
+
+
+    // Output scaling for the variable, but not for integer types.
+    // could also solve with clever templating
+
+    const scalar varScale =
+    (
+        std::is_integral<Type>::value
+      ? scalar(1)
+      : fieldScale_.getOrDefault<scalar>(fieldName, 1)
+    );
+
+    if (verbose_)
+    {
+        Info<< "Writing field " << fieldName;
+        if (!equal(varScale, 1))
+        {
+            Info<< " (scaling " << varScale << ')';
+        }
+        Info<< " to " << outputFile << endl;
+    }
+
+
+    // Implicit geometry merge()
+    tmp<Field<Type>> tfield = mergeField(localValues) * varScale;
+
+    const meshedSurf& surf = surface();
+
+    if (Pstream::master() || !parallel_)
+    {
+        const auto& values = tfield();
+
+        if (!isDir(outputFile.path()))
+        {
+            mkDir(outputFile.path());
+        }
+
+        // const scalar timeValue(0);
+
+
+        // Additional bookkeeping for decomposing non tri/quad
+        labelList decompOffsets;
+        DynamicList<face> decompFaces;
+
+
+        OFstream os(outputFile);
+
+        if (noGeometry_ || wroteGeom_)
+        {
+            // Geometry already written (or suppressed)
+            // - still need decomposition information
+
+            fileFormats::ABAQUSCore::faceDecomposition
+            (
+                surf.points(),
+                surf.faces(),
+                decompOffsets,
+                decompFaces
+            );
+        }
+        else
+        {
+            // Write geometry (separate file)
+
+            OFstream osGeom(outputFile.lessExt().ext("abq"));
+            writeGeometry(osGeom, surf, decompOffsets, decompFaces);
+        }
+
+        // Write field
+        // *DLOAD
+        // Element-based: // elemId, P, 1000
+        // Surface-based: // elemId, P4, 1000
+
+        os  << "**" << nl
+            << "** field = " << fieldName << nl
+            << "** type = " << pTraits<Type>::typeName << nl;
+
+        if (useTimeDir() && !timeName().empty())
+        {
+            os  << "** time = " << timeName() << nl;
+        }
+
+        os  << "**" << nl
+            << "*DLOAD" << nl;
+
+
+        // Regular (undecomposed) faces
+        const faceList&    faces = surf.faces();
+        const labelList& elemIds = surf.faceIds();
+
+        // Possible to use faceIds?
+        const bool useOrigFaceIds =
+        (
+            elemIds.size() == faces.size()
+         && decompFaces.empty()
+        );
+
+
+        label elemId = 0;
+
+        if (this->isPointData())
+        {
+            forAll(faces, facei)
+            {
+                if (useOrigFaceIds)
+                {
+                    // When available and not decomposed
+                    elemId = elemIds[facei];
+                }
+
+                const label beginElemId = elemId;
+
+                // Any face decomposition
+                for
+                (
+                    label decompi = decompOffsets[facei];
+                    decompi < decompOffsets[facei+1];
+                    ++decompi
+                )
+                {
+                    const face& f = decompFaces[decompi];
+
+                    Type v = Zero;
+                    for (const label verti : f)
+                    {
+                        v += values[verti];
+                    }
+                    v /= f.size();
+
+                    writeFaceValue(os, v, elemId);  // 0-based
+                    ++elemId;
+                }
+
+
+                // Face not decomposed
+                if (beginElemId == elemId)
+                {
+                    const face& f = faces[facei];
+
+                    Type v = Zero;
+                    for (const label verti : f)
+                    {
+                        v += values[verti];
+                    }
+                    v /= f.size();
+
+                    writeFaceValue(os, v, elemId);  // 0-based
+                    ++elemId;
+                }
+            }
+        }
+        else
+        {
+            auto valIter = values.cbegin();
+
+            forAll(faces, facei)
+            {
+                if (useOrigFaceIds)
+                {
+                    // When available and not decomposed
+                    elemId = elemIds[facei];
+                }
+
+                const Type v(*valIter);
+                ++valIter;
+
+                label nValues =
+                    max
+                    (
+                        label(1),
+                        (decompOffsets[facei+1] - decompOffsets[facei])
+                    );
+
+                while (nValues--)
+                {
+                    writeFaceValue(os, v, elemId);  // 0-based
+                    ++elemId;
+                }
+            }
+        }
+
+        os  << "**" << nl
+            << "**" << nl;
+    }
+
+    wroteGeom_ = true;
+    return outputFile;
+}
+
+
+// ************************************************************************* //

From 6cf815181733b657f3c316f69d4c4b0dc0e2eeb5 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Tue, 13 Oct 2020 15:48:15 +0200
Subject: [PATCH 3/8] ENH: support surfaceFieldValue on multiple faceZones or
 patches (#1874)

- additional "names" entry to specify a word/regex list of selections
  For example,
  {
      type    patch;
      name    inlets;
      names   ("inlet_[0-9].*" inlet);
  }

- if "names" exists AND contains a literal (non-regex) that can be used
  as a suitable value for "name", the "name" entry becomes optional.
  For example,
  {
      type    patch;
      names   ("inlet_[0-9].*" inlet);

      // inferred name = inlet
  }

- reduce some overhead in surfaceFieldValue

TUT: surfaceFieldValue on patches : reactingParcelFoam/verticalChannel
---
 .../surfaceFieldValue/surfaceFieldValue.C     | 397 ++++++++++++------
 .../surfaceFieldValue/surfaceFieldValue.H     |  67 +--
 .../reactingParcelFoam/verticalChannel/Allrun |   7 +-
 .../verticalChannel/Allrun-parallel           |  16 +
 .../verticalChannel/system/controlDict        |  22 +-
 .../verticalChannel/system/decomposeParDict   |  30 ++
 .../verticalChannel/system/controlDict        |   2 +-
 7 files changed, 374 insertions(+), 167 deletions(-)
 create mode 100755 tutorials/lagrangian/reactingParcelFoam/verticalChannel/Allrun-parallel
 create mode 100644 tutorials/lagrangian/reactingParcelFoam/verticalChannel/system/decomposeParDict

diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
index 1d7cd88a29..5f429e1214 100644
--- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
+++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
@@ -32,7 +32,7 @@ License
 #include "coupledPolyPatch.H"
 #include "sampledSurface.H"
 #include "mergePoints.H"
-#include "indirectPrimitivePatch.H"
+#include "uindirectPrimitivePatch.H"
 #include "PatchTools.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -130,115 +130,225 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::obr() const
 
 void Foam::functionObjects::fieldValues::surfaceFieldValue::setFaceZoneFaces()
 {
-    const label zoneId = mesh_.faceZones().findZoneID(regionName_);
+    // Indices for all matches, already sorted
+    const labelList zoneIds
+    (
+        mesh_.faceZones().indices(selectionNames_)
+    );
 
-    if (zoneId < 0)
+    // Total number of faces selected
+    label numFaces = 0;
+    for (const label zoneId : zoneIds)
     {
-        FatalErrorInFunction
-            << type() << " " << name() << ": "
-            << regionTypeNames_[regionType_] << '(' << regionName_ << "):" << nl
-            << "    Unknown face zone name: " << regionName_
-            << ". Valid face zones are: " << mesh_.faceZones().names()
-            << nl << exit(FatalError);
+        numFaces += mesh_.faceZones()[zoneId].size();
     }
 
-    const faceZone& fZone = mesh_.faceZones()[zoneId];
-
-    DynamicList<label> faceIds(fZone.size());
-    DynamicList<label> facePatchIds(fZone.size());
-    DynamicList<bool> faceFlip(fZone.size());
-
-    forAll(fZone, i)
+    if (zoneIds.empty())
     {
-        const label facei = fZone[i];
+        FatalErrorInFunction
+            << type() << ' ' << name() << ": "
+            << regionTypeNames_[regionType_] << '(' << regionName_ << "):" << nl
+            << "    No matching face zone(s): "
+            << flatOutput(selectionNames_)  << nl
+            << "    Known face zones: "
+            << flatOutput(mesh_.faceZones().names()) << nl
+            << exit(FatalError);
+    }
 
-        label faceId = -1;
-        label facePatchId = -1;
-        if (mesh_.isInternalFace(facei))
+    // Could also check this
+    #if 0
+    if (!returnReduce(bool(numFaces), orOp<bool>()))
+    {
+        WarningInFunction
+            << type() << ' ' << name() << ": "
+            << regionTypeNames_[regionType_] << '(' << regionName_ << "):" << nl
+            << "    The faceZone specification: "
+            << flatOutput(selectionNames_) << nl
+            << "    resulted in 0 faces" << nl
+            << exit(FatalError);
+    }
+    #endif
+
+    faceId_.resize(numFaces);
+    facePatchId_.resize(numFaces);
+    faceFlip_.resize(numFaces);
+
+    numFaces = 0;
+
+    for (const label zoneId : zoneIds)
+    {
+        const faceZone& fZone = mesh_.faceZones()[zoneId];
+
+        forAll(fZone, i)
         {
-            faceId = facei;
-            facePatchId = -1;
-        }
-        else
-        {
-            facePatchId = mesh_.boundaryMesh().whichPatch(facei);
-            const polyPatch& pp = mesh_.boundaryMesh()[facePatchId];
-            if (isA<coupledPolyPatch>(pp))
+            const label meshFacei = fZone[i];
+            const bool isFlip = fZone.flipMap()[i];
+
+            // Internal faces
+            label faceId = meshFacei;
+            label facePatchId = -1;
+
+            // Boundary faces
+            if (!mesh_.isInternalFace(meshFacei))
             {
-                if (refCast<const coupledPolyPatch>(pp).owner())
+                facePatchId = mesh_.boundaryMesh().whichPatch(meshFacei);
+                const polyPatch& pp = mesh_.boundaryMesh()[facePatchId];
+
+                if (isA<coupledPolyPatch>(pp))
                 {
-                    faceId = pp.whichFace(facei);
+                    if (refCast<const coupledPolyPatch>(pp).owner())
+                    {
+                        faceId = pp.whichFace(meshFacei);
+                    }
+                    else
+                    {
+                        faceId = -1;
+                    }
+                }
+                else if (!isA<emptyPolyPatch>(pp))
+                {
+                    faceId = meshFacei - pp.start();
                 }
                 else
                 {
                     faceId = -1;
+                    facePatchId = -1;
                 }
             }
-            else if (!isA<emptyPolyPatch>(pp))
-            {
-                faceId = facei - pp.start();
-            }
-            else
-            {
-                faceId = -1;
-                facePatchId = -1;
-            }
-        }
 
-        if (faceId >= 0)
-        {
-            faceIds.append(faceId);
-            facePatchIds.append(facePatchId);
-            faceFlip.append(fZone.flipMap()[i] ? true : false);
+            if (faceId >= 0)
+            {
+                faceId_[numFaces] = faceId;
+                facePatchId_[numFaces] = facePatchId;
+                faceFlip_[numFaces] = isFlip;
+
+                ++numFaces;
+            }
         }
     }
 
-    faceId_.transfer(faceIds);
-    facePatchId_.transfer(facePatchIds);
-    faceFlip_.transfer(faceFlip);
+    // Shrink to size used
+    faceId_.resize(numFaces);
+    facePatchId_.resize(numFaces);
+    faceFlip_.resize(numFaces);
     nFaces_ = returnReduce(faceId_.size(), sumOp<label>());
-
-    if (debug)
-    {
-        Pout<< "Original face zone size = " << fZone.size()
-            << ", new size = " << faceId_.size() << endl;
-    }
 }
 
 
 void Foam::functionObjects::fieldValues::surfaceFieldValue::setPatchFaces()
 {
-    const label patchid = mesh_.boundaryMesh().findPatchID(regionName_);
+    // Patch indices for all matches
+    labelList patchIds;
 
-    if (patchid < 0)
+    // Total number of faces selected
+    label numFaces = 0;
+
+    labelList selected
+    (
+        mesh_.boundaryMesh().patchSet
+        (
+            selectionNames_,
+            false  // warnNotFound - we do that ourselves
+        ).sortedToc()
+    );
+
+    DynamicList<label> bad;
+    for (const label patchi : selected)
+    {
+        const polyPatch& pp = mesh_.boundaryMesh()[patchi];
+
+        if (isA<emptyPolyPatch>(pp))
+        {
+            bad.append(patchi);
+        }
+        else
+        {
+            numFaces += pp.size();
+        }
+    }
+
+    if (bad.size())
+    {
+        label nGood = (selected.size() - bad.size());
+
+        auto& os = (nGood > 0 ? WarningInFunction : FatalErrorInFunction);
+
+        os  << "Cannot sample an empty patch" << nl;
+
+        for (const label patchi : bad)
+        {
+            os  << "    "
+                << mesh_.boundaryMesh()[patchi].name() << nl;
+        }
+
+        if (nGood)
+        {
+            os  << "No non-empty patches selected" << endl
+                << exit(FatalError);
+        }
+        else
+        {
+            os  << "Selected " << nGood << " non-empty patches" << nl;
+        }
+
+        patchIds.resize(nGood);
+        nGood = 0;
+        for (const label patchi : selected)
+        {
+            if (!bad.found(patchi))
+            {
+                patchIds[nGood] = patchi;
+                ++nGood;
+            }
+        }
+    }
+    else
+    {
+        patchIds = std::move(selected);
+    }
+
+    if (patchIds.empty())
     {
         FatalErrorInFunction
-            << type() << " " << name() << ": "
+            << type() << ' ' << name() << ": "
             << regionTypeNames_[regionType_] << '(' << regionName_ << "):" << nl
-            << "    Unknown patch name: " << regionName_
-            << ". Valid patch names are: "
+            << "    No matching patch name(s): "
+            << flatOutput(selectionNames_)  << nl
+            << "    Known patch names:" << nl
             << mesh_.boundaryMesh().names() << nl
             << exit(FatalError);
     }
 
-    const polyPatch& pp = mesh_.boundaryMesh()[patchid];
-
-    label nFaces = pp.size();
-    if (isA<emptyPolyPatch>(pp))
+    // Could also check this
+    #if 0
+    if (!returnReduce(bool(numFaces), orOp<bool>()))
     {
-        nFaces = 0;
+        WarningInFunction
+            << type() << ' ' << name() << ": "
+            << regionTypeNames_[regionType_] << '(' << regionName_ << "):" << nl
+            << "    The patch specification: "
+            << flatOutput(selectionNames_) << nl
+            << "    resulted in 0 faces" << nl
+            << exit(FatalError);
     }
+    #endif
 
-    faceId_.setSize(nFaces);
-    facePatchId_.setSize(nFaces);
-    faceFlip_.setSize(nFaces);
+    faceId_.resize(numFaces);
+    facePatchId_.resize(numFaces);
+    faceFlip_.resize(numFaces);
     nFaces_ = returnReduce(faceId_.size(), sumOp<label>());
 
-    forAll(faceId_, facei)
+    numFaces = 0;
+    for (const label patchi : patchIds)
     {
-        faceId_[facei] = facei;
-        facePatchId_[facei] = patchid;
-        faceFlip_[facei] = false;
+        const polyPatch& pp = mesh_.boundaryMesh()[patchi];
+        const label len = pp.size();
+
+        SubList<label>(faceId_, len, numFaces) = identity(len);
+        SubList<label>(facePatchId_, len, numFaces) = patchi;
+        SubList<bool>(faceFlip_, len, numFaces) = false;
+
+        numFaces += len;
     }
 }
 
@@ -252,24 +362,25 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::combineMeshGeometry
     List<faceList> allFaces(Pstream::nProcs());
     List<pointField> allPoints(Pstream::nProcs());
 
-    labelList globalFacesIs(faceId_);
-    forAll(globalFacesIs, i)
     {
-        if (facePatchId_[i] != -1)
+        IndirectList<face> selectedFaces(mesh_.faces(), labelList(faceId_));
+        labelList& meshFaceIds = selectedFaces.addressing();
+
+        forAll(meshFaceIds, i)
         {
             const label patchi = facePatchId_[i];
-            globalFacesIs[i] += mesh_.boundaryMesh()[patchi].start();
+            if (patchi != -1)
+            {
+                meshFaceIds[i] += mesh_.boundaryMesh()[patchi].start();
+            }
         }
-    }
 
-    // Add local faces and points to the all* lists
-    indirectPrimitivePatch pp
-    (
-        IndirectList<face>(mesh_.faces(), globalFacesIs),
-        mesh_.points()
-    );
-    allFaces[Pstream::myProcNo()] = pp.localFaces();
-    allPoints[Pstream::myProcNo()] = pp.localPoints();
+        // Add local faces and points to the all* lists
+        uindirectPrimitivePatch pp(selectedFaces, mesh_.points());
+
+        allFaces[Pstream::myProcNo()] = pp.localFaces();
+        allPoints[Pstream::myProcNo()] = pp.localPoints();
+    }
 
     Pstream::gatherList(allFaces);
     Pstream::gatherList(allPoints);
@@ -283,53 +394,41 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::combineMeshGeometry
         nPoints += allPoints[proci].size();
     }
 
-    faces.setSize(nFaces);
-    points.setSize(nPoints);
+    faces.resize(nFaces);
+    points.resize(nPoints);
 
     nFaces = 0;
     nPoints = 0;
 
-    // My own data first
+    // My data first
     {
-        const faceList& fcs = allFaces[Pstream::myProcNo()];
-        for (const face& f : fcs)
+        for (const face& f : allFaces[Pstream::myProcNo()])
         {
-            face& newF = faces[nFaces++];
-            newF.setSize(f.size());
-            forAll(f, fp)
-            {
-                newF[fp] = f[fp] + nPoints;
-            }
+            faces[nFaces++] = offsetOp<face>()(f, nPoints);
         }
 
-        const pointField& pts = allPoints[Pstream::myProcNo()];
-        for (const point& pt : pts)
+        for (const point& p : allPoints[Pstream::myProcNo()])
         {
-            points[nPoints++] = pt;
+            points[nPoints++] = p;
         }
     }
 
     // Other proc data follows
     forAll(allFaces, proci)
     {
-        if (proci != Pstream::myProcNo())
+        if (proci == Pstream::myProcNo())
         {
-            const faceList& fcs = allFaces[proci];
-            for (const face& f : fcs)
-            {
-                face& newF = faces[nFaces++];
-                newF.setSize(f.size());
-                forAll(f, fp)
-                {
-                    newF[fp] = f[fp] + nPoints;
-                }
-            }
+            continue;
+        }
 
-            const pointField& pts = allPoints[proci];
-            for (const point& pt : pts)
-            {
-                points[nPoints++] = pt;
-            }
+        for (const face& f : allFaces[proci])
+        {
+            faces[nFaces++] = offsetOp<face>()(f, nPoints);
+        }
+
+        for (const point& p : allPoints[proci])
+        {
+            points[nPoints++] = p;
         }
     }
 
@@ -383,11 +482,7 @@ combineSurfaceGeometry
             PatchTools::gatherAndMerge
             (
                 mergeDim,
-                primitivePatch
-                (
-                    SubList<face>(s.faces(), s.faces().size()),
-                    s.points()
-                ),
+                primitivePatch(SubList<face>(s.faces()), s.points()),
                 points,
                 faces,
                 pointsMap
@@ -413,11 +508,7 @@ combineSurfaceGeometry
             PatchTools::gatherAndMerge
             (
                 mergeDim,
-                primitivePatch
-                (
-                    SubList<face>(s.faces(), s.faces().size()),
-                    s.points()
-                ),
+                primitivePatch(SubList<face>(s.faces()), s.points()),
                 points,
                 faces,
                 pointsMap
@@ -522,7 +613,7 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::update()
     if (nFaces_ == 0)
     {
         FatalErrorInFunction
-            << type() << " " << name() << ": "
+            << type() << ' ' << name() << ": "
             << regionTypeNames_[regionType_] << '(' << regionName_ << "):" << nl
             << "    Region has no faces" << exit(FatalError);
     }
@@ -548,7 +639,7 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::writeFileHeader
     if (canWriteHeader() && (operation_ != opNone))
     {
         writeCommented(os, "Region type : ");
-        os << regionTypeNames_[regionType_] << " " << regionName_ << endl;
+        os << regionTypeNames_[regionType_] << ' ' << regionName_ << nl;
 
         writeHeaderValue(os, "Faces", nFaces_);
         writeHeaderValue(os, "Area", totalArea_);
@@ -570,7 +661,7 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::writeFileHeader
         for (const word& fieldName : fields_)
         {
             os  << tab << operationTypeNames_[operation_]
-                << "(" << fieldName << ")";
+                << '(' << fieldName << ')';
         }
 
         os  << endl;
@@ -821,9 +912,10 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::surfaceFieldValue
             true  // Failsafe behaviour
         )
     ),
-    weightFieldName_("none"),
     needsUpdate_(true),
     writeArea_(false),
+    selectionNames_(),
+    weightFieldName_("none"),
     totalArea_(0),
     nFaces_(0),
     faceId_(),
@@ -854,9 +946,10 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::surfaceFieldValue
             true  // Failsafe behaviour
         )
     ),
-    weightFieldName_("none"),
     needsUpdate_(true),
     writeArea_(false),
+    selectionNames_(),
+    weightFieldName_("none"),
     totalArea_(0),
     nFaces_(0),
     faceId_(),
@@ -876,18 +969,54 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::read
 {
     fieldValue::read(dict);
 
-    weightFieldName_ = "none";
     needsUpdate_ = true;
     writeArea_ = dict.getOrDefault("writeArea", false);
+    weightFieldName_ = "none";
     totalArea_ = 0;
     nFaces_ = 0;
     faceId_.clear();
     facePatchId_.clear();
     faceFlip_.clear();
-    sampledPtr_.clear();
-    surfaceWriterPtr_.clear();
+    sampledPtr_.reset(nullptr);
+    surfaceWriterPtr_.reset(nullptr);
+
+    // Can have "name" (word) and/or "names" (wordRes)
+    //
+    // If "names" exists AND contains a literal (non-regex) that can be used
+    // as a suitable value for "name", the "name" entry becomes optional.
+
+    regionName_.clear();
+    selectionNames_.clear();
+
+    {
+        dict.readIfPresent("names", selectionNames_);
+
+        for (const auto& item : selectionNames_)
+        {
+            if (item.isLiteral())
+            {
+                regionName_ = item;
+                break;
+            }
+        }
+
+        // Mandatory if we didn't pick up a value from selectionNames_
+        dict.readEntry
+        (
+            "name",
+            regionName_,
+            keyType::LITERAL,
+            regionName_.empty()
+        );
+
+        // Ensure there is always content for selectionNames_
+        if (selectionNames_.empty())
+        {
+            selectionNames_.resize(1);
+            selectionNames_.first() = regionName_;
+        }
+    }
 
-    dict.readEntry("name", regionName_);
 
     // Create sampled surface, but leave 'expired' (ie, no update) since it
     // may depend on fields or data that do not yet exist
@@ -901,7 +1030,7 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::read
         );
     }
 
-    Info<< type() << " " << name() << ":" << nl
+    Info<< type() << ' ' << name() << ':' << nl
         << "    operation     = ";
 
     if (postOperation_ != postOpNone)
diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H
index 1dde638e55..5687e9f4a4 100644
--- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H
+++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H
@@ -31,17 +31,17 @@ Group
     grpFieldFunctionObjects
 
 Description
-    Provides a 'face regionType' variant of the \c fieldValues function object.
+    A \c face regionType variant of the \c fieldValues function object.
 
     Given a list of user-specified fields and a selection of mesh (or general
     surface) faces, a number of operations can be performed, such as sums,
     averages and integrations.
 
     For example, to calculate the volumetric or mass flux across a patch,
-    apply the 'sum' operator to the flux field (typically \c phi).
+    apply the \c sum operator to the flux field (typically \c phi).
 
 Usage
-    Minimal example by using \c system/controlDict.functions:
+    A minimal example:
     \verbatim
     surfaceFieldValuePatch1
     {
@@ -56,6 +56,8 @@ Usage
         name            <patch>;
 
         // Optional entries (runtime modifiable)
+        names           (<patch-name> <patch-regex>);
+
         postOperation   none;
         weightField     alpha1;
         scaleFactor     1.0;
@@ -79,6 +81,8 @@ Usage
         name            <faceZone>;
 
         // Optional entries (runtime modifiable)
+        names           (<zone-name> <zone-regex>);
+
         postOperation   none;
         weightField     alpha1;
         scaleFactor     1.0;
@@ -92,12 +96,13 @@ Usage
 
     where the entries mean:
     \table
-      Property     | Description                        | Type | Req'd | Dflt
+      Property     | Description                        | Type | Reqd | Dflt
       type         | Type name: surfaceFieldValue       | word |  yes  | -
-      libs         | Library name: fieldFunctionObjects | word |  yes  | -
-      fields       | Names of operand fields            | wordList | yes | -
+      libs         | Libraries: fieldFunctionObjects    | wordList |  yes  | -
       regionType   | Face regionType: see below         | word |  yes  | -
-      name         | Name for regionType                | word |  yes  | -
+      fields       | Names of operand fields            | wordList | yes | -
+      name         | Name of the regionType             | word |  yes  | -
+      names        | Extended selection                 | word/regex list | no | -
       operation    | Operation type: see below          | word |  yes  | -
       postOperation | Post-operation type: see below    | word |  no   | none
       weightField  | Name of field to apply weighting   | word |  no   | none
@@ -112,9 +117,9 @@ Usage
 
     Options for the \c regionType entry:
     \plaintable
-      faceZone              | The \b name entry to specify the faceZone
-      patch                 | The \b name entry to specify the patch
-      functionObjectSurface | The \b name entry to specify a polySurface
+      faceZone      | The \b name entry specifies a faceZone. Supports \b names
+      patch         | The \b name entry specifies a patch. Supports \b names
+      functionObjectSurface | The \b name entry specifies a polySurface
       sampledSurface        | A \b sampledSurfaceDict sub-dictionary and \b name
     \endplaintable
 
@@ -156,12 +161,16 @@ Usage
     Usage by the \c postProcess utility is not available.
 
 Note
+    - Some types (eg, patch or faceZone) support the selection of multiple
+      entries, which can be specified as list of words or regular expressions
+      by \b names entry.<br>
+      If the \b names enty exists \em and contains a literal that can be used
+      as a suitable value for \b name, the \b name entry becomes optional
+      instead of being mandatory.
     - The values reported by the \c areaNormalAverage and \c areaNormalIntegrate
       operations are written as the first component of a field with the same
       rank as the input field.
-    - Faces on empty patches get ignored.
-    - If the field is a volField the \c faceZone
-      can only consist of boundary faces.
+    - Faces on empty patches are ignored.
     - Using \c functionObjectSurface:
       - The keyword %subRegion should not be used to select surfaces.
         Instead specify the regionType 'functionObjectSurface' and provide
@@ -242,8 +251,8 @@ public:
         //- Region type enumeration
         enum regionTypes
         {
-            stFaceZone = 0x01,   //!< Calculate on a faceZone
-            stPatch    = 0x02,   //!< Calculate on a patch
+            stFaceZone = 0x01,   //!< Calculate with faceZone(s)
+            stPatch    = 0x02,   //!< Calculate with patch(es)
             stObject   = 0x11,   //!< Calculate with function object surface
             stSampled  = 0x12    //!< Sample onto surface and calculate
         };
@@ -385,15 +394,18 @@ protected:
         //- Optional post-evaluation operation
         postOperationType postOperation_;
 
-        //- Weight field name - optional
-        word weightFieldName_;
-
         //- Track if the surface needs an update
         bool needsUpdate_;
 
         //- Optionally write the area of the surfaceFieldValue
         bool writeArea_;
 
+        //- Extended selections
+        wordRes selectionNames_;
+
+        //- Weight field name - optional
+        word weightFieldName_;
+
         //- Total area of the surfaceFieldValue
         scalar totalArea_;
 
@@ -401,17 +413,20 @@ protected:
         label nFaces_;
 
 
-        // If operating on mesh faces (faceZone, patch)
+    // If operating on mesh faces (faceZone, patch)
 
-            //- Local list of face IDs
-            labelList faceId_;
+        //- Local list of face IDs
+        labelList faceId_;
 
-            //- Local list of patch ID per face
-            labelList facePatchId_;
+        //- Local list of patch ID per face
+        labelList facePatchId_;
 
-            //- List representing the face flip map
-            //  (false: use as-is, true: negate)
-            boolList faceFlip_;
+        //- List representing the face flip map
+        //  (false: use as-is, true: negate)
+        boolList faceFlip_;
+
+
+    // Demand-driven
 
         //- The sampledSurface (when operating on sampledSurface)
         autoPtr<sampledSurface> sampledPtr_;
diff --git a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/Allrun b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/Allrun
index 9ec0596bab..a5d9a2cb05 100755
--- a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/Allrun
+++ b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/Allrun
@@ -3,15 +3,12 @@ cd "${0%/*}" || exit                                # Run from this directory
 . ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
 #------------------------------------------------------------------------------
 
-# create mesh
-runApplication blockMesh
-
 restore0Dir
 
-# initialise with potentialFoam solution
+runApplication blockMesh
+
 runApplication potentialFoam
 
-# run the solver
 runApplication $(getApplication)
 
 #------------------------------------------------------------------------------
diff --git a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/Allrun-parallel b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/Allrun-parallel
new file mode 100755
index 0000000000..8bb64db819
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/Allrun-parallel
@@ -0,0 +1,16 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+restore0Dir
+
+runApplication blockMesh
+
+runApplication potentialFoam
+
+runApplication decomposePar
+
+runParallel $(getApplication)
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/system/controlDict b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/system/controlDict
index 9ac52b9bae..3442d41a3a 100644
--- a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/system/controlDict
+++ b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/system/controlDict
@@ -53,7 +53,27 @@ maxDeltaT       1e-03;
 
 functions
 {
-    surfaceFieldValue1
+    avgInlets
+    {
+        type            surfaceFieldValue;
+        libs            (fieldFunctionObjects);
+        enabled         yes;
+        writeControl    writeTime;
+        log             yes;
+        writeFields     off;
+        surfaceFormat   vtk;
+        regionType      patch;
+        // name            inlet;
+        names           (inlet "inlet.*");
+        operation       weightedAverage;
+        weightField     phi;
+        fields
+        (
+            T
+        );
+    }
+
+    avgOutlets
     {
         type            surfaceFieldValue;
         libs            (fieldFunctionObjects);
diff --git a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/system/decomposeParDict b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/system/decomposeParDict
new file mode 100644
index 0000000000..e12fff0f4b
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/system/decomposeParDict
@@ -0,0 +1,30 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2006                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 4;
+
+method          scotch;
+
+coeffs
+{
+    n           (2 2 1);
+}
+
+distributed     no;
+
+roots           ( );
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/system/controlDict b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/system/controlDict
index b12d88a3d2..ae2bb7accb 100644
--- a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/system/controlDict
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/system/controlDict
@@ -48,7 +48,7 @@ runTimeModifiable yes;
 
 functions
 {
-    surfaceFieldValue1
+    avgOutlets
     {
         type            surfaceFieldValue;
         libs            (fieldFunctionObjects);

From 90fe49e8e53870ea484ed3f3952a39aca6e33437 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Tue, 13 Oct 2020 15:52:20 +0200
Subject: [PATCH 4/8] ENH: adjust sampledPatch (#1874) for similarity with
 surfaceFieldValue

---
 .../sampledPatch/sampledPatch.C               | 142 ++++++++++++------
 .../sampledPatch/sampledPatch.H               |  42 +++---
 .../sampledPatch/sampledPatchTemplates.C      |  17 ++-
 3 files changed, 128 insertions(+), 73 deletions(-)

diff --git a/src/sampling/sampledSurface/sampledPatch/sampledPatch.C b/src/sampling/sampledSurface/sampledPatch/sampledPatch.C
index eb176bf654..360000f6c4 100644
--- a/src/sampling/sampledSurface/sampledPatch/sampledPatch.C
+++ b/src/sampling/sampledSurface/sampledPatch/sampledPatch.C
@@ -32,6 +32,7 @@ License
 #include "polyPatch.H"
 #include "volFields.H"
 #include "surfaceFields.H"
+#include "uindirectPrimitivePatch.H"
 
 #include "addToRunTimeSelectionTable.H"
 
@@ -55,7 +56,7 @@ Foam::sampledPatch::sampledPatch
 )
 :
     sampledSurface(name, mesh),
-    patchNames_(patchNames),
+    selectionNames_(patchNames),
     triangulate_(triangulate),
     needsUpdate_(true)
 {}
@@ -69,7 +70,7 @@ Foam::sampledPatch::sampledPatch
 )
 :
     sampledSurface(name, mesh, dict),
-    patchNames_(dict.get<wordRes>("patches")),
+    selectionNames_(dict.get<wordRes>("patches")),
     triangulate_(dict.getOrDefault("triangulate", false)),
     needsUpdate_(true)
 {}
@@ -81,8 +82,63 @@ const Foam::labelList& Foam::sampledPatch::patchIDs() const
 {
     if (patchIDs_.empty())
     {
-        patchIDs_ = mesh().boundaryMesh().patchSet(patchNames_).sortedToc();
+        labelList selected
+        (
+            mesh().boundaryMesh().patchSet(selectionNames_).sortedToc()
+        );
+
+        DynamicList<label> bad;
+        for (const label patchi : selected)
+        {
+            const polyPatch& pp = mesh().boundaryMesh()[patchi];
+
+            if (isA<emptyPolyPatch>(pp))
+            {
+                bad.append(patchi);
+            }
+        }
+
+        if (bad.size())
+        {
+            label nGood = (selected.size() - bad.size());
+
+            auto& os = nGood > 0 ? WarningInFunction : FatalErrorInFunction;
+
+            os  << "Cannot sample an empty patch" << nl;
+
+            for (const label patchi : bad)
+            {
+                os  << "    "
+                    << mesh().boundaryMesh()[patchi].name() << nl;
+            }
+
+            if (nGood)
+            {
+                os  << "No non-empty patches selected" << endl
+                    << exit(FatalError);
+            }
+            else
+            {
+                os  << "Selected " << nGood << " non-empty patches" << nl;
+            }
+
+            patchIDs_.resize(nGood);
+            nGood = 0;
+            for (const label patchi : selected)
+            {
+                if (!bad.found(patchi))
+                {
+                    patchIDs_[nGood] = patchi;
+                    ++nGood;
+                }
+            }
+        }
+        else
+        {
+            patchIDs_ = std::move(selected);
+        }
     }
+
     return patchIDs_;
 }
 
@@ -102,11 +158,13 @@ bool Foam::sampledPatch::expire()
     }
 
     sampledSurface::clearGeom();
-    MeshStorage::clear();
+    Mesh::clear();
+
     patchIDs_.clear();
+    patchStart_.clear();
+
     patchIndex_.clear();
     patchFaceLabels_.clear();
-    patchStart_.clear();
 
     needsUpdate_ = true;
     return true;
@@ -120,52 +178,45 @@ bool Foam::sampledPatch::update()
         return false;
     }
 
-    label sz = 0;
+    // Total number of faces selected
+    label numFaces = 0;
     for (const label patchi : patchIDs())
     {
         const polyPatch& pp = mesh().boundaryMesh()[patchi];
-
-        if (isA<emptyPolyPatch>(pp))
-        {
-            FatalErrorInFunction
-                << "Cannot sample an empty patch. Patch " << pp.name()
-                << exit(FatalError);
-        }
-
-        sz += pp.size();
+        numFaces += pp.size();
     }
 
-    // For every face (or triangle) the originating patch and local face in the
-    // patch.
-    patchIndex_.setSize(sz);
-    patchFaceLabels_.setSize(sz);
-    patchStart_.setSize(patchIDs().size());
-    labelList meshFaceLabels(sz);
+    patchStart_.resize(patchIDs().size());
 
-    sz = 0;
+    // The originating patch and local face in the patch.
+    patchIndex_.resize(numFaces);
+    patchFaceLabels_.resize(numFaces);
 
-    forAll(patchIDs(), i)
+    IndirectList<face> selectedFaces(mesh().faces(), labelList());
+    labelList& meshFaceIds = selectedFaces.addressing();
+    meshFaceIds.resize(numFaces);
+
+    numFaces = 0;
+
+    forAll(patchIDs(), idx)
     {
-        const label patchi = patchIDs()[i];
-
-        patchStart_[i] = sz;
-
+        const label patchi = patchIDs()[idx];
         const polyPatch& pp = mesh().boundaryMesh()[patchi];
+        const label len = pp.size();
 
-        forAll(pp, j)
-        {
-            patchIndex_[sz] = i;
-            patchFaceLabels_[sz] = j;
-            meshFaceLabels[sz] = pp.start()+j;
-            ++sz;
-        }
+        patchStart_[idx] = numFaces;
+
+        SubList<label>(patchIndex_, len, numFaces) = idx;
+
+        SubList<label>(patchFaceLabels_, len, numFaces) = identity(len);
+
+        SubList<label>(meshFaceIds, len, numFaces) = identity(len, pp.start());
+
+        numFaces += len;
     }
 
-    indirectPrimitivePatch allPatches
-    (
-        IndirectList<face>(mesh().faces(), meshFaceLabels),
-        mesh().points()
-    );
+
+    uindirectPrimitivePatch allPatches(selectedFaces, mesh().points());
 
     this->storedPoints() = allPatches.localPoints();
     this->storedFaces()  = allPatches.localFaces();
@@ -177,7 +228,7 @@ bool Foam::sampledPatch::update()
     // too often anyhow.
     if (triangulate_)
     {
-        MeshStorage::triangulate();
+        Mesh::triangulate();
     }
 
     if (debug)
@@ -194,10 +245,9 @@ bool Foam::sampledPatch::update()
 // remap action on triangulation
 void Foam::sampledPatch::remapFaces(const labelUList& faceMap)
 {
-    // Recalculate the cells cut
     if (!faceMap.empty())
     {
-        MeshStorage::remapFaces(faceMap);
+        Mesh::remapFaces(faceMap);
         patchFaceLabels_ = labelList
         (
             labelUIndList(patchFaceLabels_, faceMap)
@@ -207,11 +257,11 @@ void Foam::sampledPatch::remapFaces(const labelUList& faceMap)
             labelUIndList(patchIndex_, faceMap)
         );
 
-        // Redo patchStart.
-        if (patchIndex_.size() > 0)
+        // Update patchStart
+        if (patchIndex_.size())
         {
             patchStart_[patchIndex_[0]] = 0;
-            for (label i = 1; i < patchIndex_.size(); i++)
+            for (label i = 1; i < patchIndex_.size(); ++i)
             {
                 if (patchIndex_[i] != patchIndex_[i-1])
                 {
@@ -367,7 +417,7 @@ Foam::tmp<Foam::tensorField> Foam::sampledPatch::interpolate
 void Foam::sampledPatch::print(Ostream& os) const
 {
     os  << "sampledPatch: " << name() << " :"
-        << "  patches:" << patchNames()
+        << "  patches: " << flatOutput(selectionNames_)
         << "  faces:" << faces().size()
         << "  points:" << points().size();
 }
diff --git a/src/sampling/sampledSurface/sampledPatch/sampledPatch.H b/src/sampling/sampledSurface/sampledPatch/sampledPatch.H
index b2e035eebe..2d47876496 100644
--- a/src/sampling/sampledSurface/sampledPatch/sampledPatch.H
+++ b/src/sampling/sampledSurface/sampledPatch/sampledPatch.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -62,8 +62,7 @@ SourceFiles
 #define sampledPatch_H
 
 #include "sampledSurface.H"
-#include "MeshedSurface.H"
-#include "wordRes.H"
+#include "MeshedSurfaces.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -76,19 +75,19 @@ namespace Foam
 
 class sampledPatch
 :
-    public MeshedSurface<face>,
+    public meshedSurface,
     public sampledSurface
 {
-    //- Private typedefs for convenience
-    typedef MeshedSurface<face> MeshStorage;
+    //- Mesh storage type
+    typedef meshedSurface Mesh;
 
 
-    // Private data
+    // Private Data
 
-        //- Name of patches
-        const wordRes patchNames_;
+        //- Selection (word/regex) of patches
+        const wordRes selectionNames_;
 
-        //- Corresponding patchIDs
+        //- The IDs for selected patches. Sorted and no emptyPolyPatch
         mutable labelList patchIDs_;
 
         //- Triangulated faces or keep faces as is
@@ -97,15 +96,15 @@ class sampledPatch
         //- Track if the surface needs an update
         mutable bool needsUpdate_;
 
+        //- Start indices (in patchFaceLabels_) of patches
+        labelList patchStart_;
+
         //- For every face (or triangle) the originating patch
         labelList patchIndex_;
 
         //- For every face (or triangle) the index in the originating patch
         labelList patchFaceLabels_;
 
-        //- Start indices (in patchFaceLabels_) of patches
-        labelList patchStart_;
-
 
     // Private Member Functions
 
@@ -137,23 +136,28 @@ class sampledPatch
 
 protected:
 
+        //- The selection (word/regex) of patches
         const wordRes& patchNames() const
         {
-            return patchNames_;
+            return selectionNames_;
         }
 
+        //- The patches selected
         const labelList& patchIDs() const;
 
+        //- The offset into patchIndex, patchFaceLabels
         const labelList& patchStart() const
         {
             return patchStart_;
         }
 
+        //- For each face, the patch ID
         const labelList& patchIndex() const
         {
             return patchIndex_;
         }
 
+        //- For each face, the patch local face ID
         const labelList& patchFaceLabels() const
         {
             return patchFaceLabels_;
@@ -208,13 +212,13 @@ public:
         //- Points of surface
         virtual const pointField& points() const
         {
-            return MeshStorage::points();
+            return Mesh::points();
         }
 
         //- Faces of surface
         virtual const faceList& faces() const
         {
-            return MeshStorage::surfFaces();
+            return Mesh::surfFaces();
         }
 
         //- Per-face zone/region information
@@ -226,19 +230,19 @@ public:
         //- Face area vectors
         virtual const vectorField& Sf() const
         {
-            return MeshStorage::Sf();
+            return Mesh::Sf();
         }
 
         //- Face area magnitudes
         virtual const scalarField& magSf() const
         {
-            return MeshStorage::magSf();
+            return Mesh::magSf();
         }
 
         //- Face centres
         virtual const vectorField& Cf() const
         {
-            return MeshStorage::Cf();
+            return Mesh::Cf();
         }
 
 
diff --git a/src/sampling/sampledSurface/sampledPatch/sampledPatchTemplates.C b/src/sampling/sampledSurface/sampledPatch/sampledPatchTemplates.C
index 8ff986a05a..b818c2d860 100644
--- a/src/sampling/sampledSurface/sampledPatch/sampledPatchTemplates.C
+++ b/src/sampling/sampledSurface/sampledPatch/sampledPatchTemplates.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -93,20 +93,21 @@ Foam::sampledPatch::sampleOnPoints
 
     bitSet pointDone(points().size());
 
-    forAll(faces(), cutFacei)
+    forAll(faces(), i)
     {
-        const label patchi = patchIDs_[patchIndex_[cutFacei]];
+        const face& f = faces()[i];
+        const label patchi = patchIDs_[patchIndex_[i]];
+        const label patchFacei = patchFaceLabels_[i];
+
         const polyPatch& pp = mesh().boundaryMesh()[patchi];
-        const label patchFacei = patchFaceLabels()[cutFacei];
-        const face& f = faces()[cutFacei];
+
+        const label facei = patchFacei + pp.start();
+        const label celli = own[facei];
 
         for (const label pointi : f)
         {
             if (pointDone.set(pointi))
             {
-                const label facei = patchFacei + pp.start();
-                const label celli = own[facei];
-
                 values[pointi] = interpolator.interpolate
                 (
                     points()[pointi],

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 5/8] 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 2d2796dcca..5a9670931b 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 0000000000..949fc14141
--- /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 0000000000..69b648385e
--- /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 0000000000..5798324490
--- /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 9d907d286a..9830373ee5 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 bc4bdfb094..bfc4ca4c69 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 a7e2611540..940382c4de 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 f4a3fab3a9..597eb33fbc 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 68777ae5cf..492ce24a65 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 75a495fdcc..5c430ef5da 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;
+    }
 );
 
 

From f8d08a805ba9ac59e4a64e84ea0657dfc8956a6c Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Wed, 21 Oct 2020 22:10:41 +0200
Subject: [PATCH 6/8] ENH: additional handling for out-of-range sampling
 (#1891)

- when sampling onto a meshed surface, the sampling surface may be
  outside of the mesh region, or simply too far away to be considered
  reasonable.

  Can now specify a max search distance and default values for samples
  that are too distant.
  If a default value is not specified, uses Type(Zero).

  Eg,

      maxDistance     0.005;
      defaultValue
      {
          "p.*"   1e5;
          T       273.15;
          U       (-100 -100 -100);
      }
---
 .../sampledMeshedSurface.C                    | 144 +++++++++++-------
 .../sampledMeshedSurface.H                    |  67 +++++---
 .../sampledMeshedSurfaceTemplates.C           |  93 +++++++----
 .../sampledSurface/sampledSurface.H           |   8 +-
 .../sampledSurface/sampledSurfaceTemplates.C  |  16 +-
 .../rhoSimpleFoam/squareBend/system/sampling  |  38 +++++
 6 files changed, 250 insertions(+), 116 deletions(-)

diff --git a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C
index 9830373ee5..2eed7d3d43 100644
--- a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C
+++ b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C
@@ -170,7 +170,7 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
     // Does approximation by looking at the face centres only
     const pointField& fc = surface_.faceCentres();
 
-    List<nearInfo> nearest(fc.size(), nearInfo(GREAT, labelMax));
+    List<nearInfo> nearest(fc.size(), nearInfo(Foam::sqr(GREAT), labelMax));
 
     if (sampleSource_ == samplingSource::cells)
     {
@@ -181,13 +181,14 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
         forAll(fc, facei)
         {
             const point& pt = fc[facei];
+            auto& near = nearest[facei];
 
-            pointIndexHit info = cellTree.findNearest(pt, sqr(GREAT));
+            pointIndexHit info = cellTree.findNearest(pt, Foam::sqr(GREAT));
 
             if (info.hit())
             {
-                nearest[facei].first()  = magSqr(info.hitPoint()-pt);
-                nearest[facei].second() = globalCells.toGlobal(info.index());
+                near.first()  = magSqr(info.hitPoint()-pt);
+                near.second() = globalCells.toGlobal(info.index());
             }
         }
     }
@@ -200,14 +201,15 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
         forAll(fc, facei)
         {
             const point& pt = fc[facei];
+            auto& near = nearest[facei];
 
             if (cellTree.bb().contains(pt))
             {
                 const label index = cellTree.findInside(pt);
                 if (index != -1)
                 {
-                    nearest[facei].first()  = 0;
-                    nearest[facei].second() = globalCells.toGlobal(index);
+                    near.first()  = 0;
+                    near.second() = globalCells.toGlobal(index);
                 }
             }
         }
@@ -222,13 +224,14 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
         forAll(fc, facei)
         {
             const point& pt = fc[facei];
+            auto& near = nearest[facei];
 
             pointIndexHit info = bndTree.findNearest(pt, sqr(GREAT));
 
             if (info.hit())
             {
-                nearest[facei].first()  = magSqr(info.hitPoint()-pt);
-                nearest[facei].second() =
+                near.first()  = magSqr(info.hitPoint()-pt);
+                near.second() =
                     globalCells.toGlobal
                     (
                         bndTree.shapes().faceLabels()[info.index()]
@@ -250,16 +253,26 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
 
     forAll(nearest, facei)
     {
-        const label index = nearest[facei].second();
+        const auto& near = nearest[facei];
+
+        const label index = near.second();
 
         if (index == labelMax)
         {
-            // Not found on any processor. How to map?
+            // Not found on any processor, or simply too far.
+            // How to map?
         }
         else if (globalCells.isLocal(index))
         {
-            cellOrFaceLabels[facei] = globalCells.toLocal(index);
             facesToSubset.set(facei);
+
+            // Mark as special if too far away
+            cellOrFaceLabels[facei] =
+            (
+                (near.first() < maxDistanceSqr_)
+              ? globalCells.toLocal(index)
+              : -1
+            );
         }
     }
 
@@ -295,7 +308,7 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
     {
         // With point interpolation
 
-        samplePoints_.resize(pointMap.size());
+        samplePoints_ = points();
         sampleElements_.resize(pointMap.size(), -1);
 
         // Store any face per point (without using pointFaces())
@@ -314,32 +327,25 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
 
         if (sampleSource_ == samplingSource::cells)
         {
-            // samplePoints_   : per surface point a location inside the cell
-            // sampleElements_ : per surface point the cell
+            // sampleElements_ : per surface point, the cell
+            // samplePoints_   : per surface point, a location inside the cell
 
-            forAll(points(), pointi)
+            forAll(samplePoints_, pointi)
             {
-                const point& pt = points()[pointi];
+                // Use _copy_ of point
+                const point pt = samplePoints_[pointi];
 
                 const label celli = cellOrFaceLabels[pointToFace[pointi]];
 
                 sampleElements_[pointi] = celli;
 
-                // Check if point inside cell
                 if
                 (
-                    mesh().pointInCell
-                    (
-                        pt,
-                        sampleElements_[pointi],
-                        meshSearcher.decompMode()
-                    )
+                    celli >= 0
+                 && !mesh().pointInCell(pt, celli, meshSearcher.decompMode())
                 )
                 {
-                    samplePoints_[pointi] = pt;
-                }
-                else
-                {
+                    // Point not inside cell
                     // Find nearest point on faces of cell
 
                     scalar minDistSqr = VGREAT;
@@ -348,7 +354,12 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
                     {
                         const face& f = mesh().faces()[facei];
 
-                        pointHit info = f.nearestPoint(pt, mesh().points());
+                        pointHit info =
+                            f.nearestPoint
+                            (
+                                pt,
+                                mesh().points()
+                            );
 
                         if (info.distance() < minDistSqr)
                         {
@@ -361,50 +372,52 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
         }
         else if (sampleSource_ == samplingSource::insideCells)
         {
-            // samplePoints_   : per surface point a location inside the cell
             // sampleElements_ : per surface point the cell
+            // samplePoints_   : per surface point a location inside the cell
 
-            forAll(points(), pointi)
+            forAll(samplePoints_, pointi)
             {
-                const point& pt = points()[pointi];
-
                 const label celli = cellOrFaceLabels[pointToFace[pointi]];
 
                 sampleElements_[pointi] = celli;
-                samplePoints_[pointi] = pt;
             }
         }
         else  // samplingSource::boundaryFaces
         {
-            // samplePoints_   : per surface point a location on the boundary
-            // sampleElements_ : per surface point the boundary face containing
+            // sampleElements_ : per surface point, the boundary face containing
             //                   the location
+            // samplePoints_   : per surface point, a location on the boundary
 
-            forAll(points(), pointi)
+            forAll(samplePoints_, pointi)
             {
-                const point& pt = points()[pointi];
+                const point& pt = samplePoints_[pointi];
 
                 const label facei = cellOrFaceLabels[pointToFace[pointi]];
 
                 sampleElements_[pointi] = facei;
-                samplePoints_[pointi] = mesh().faces()[facei].nearestPoint
-                (
-                    pt,
-                    mesh().points()
-                ).rawPoint();
+
+                if (facei >= 0)
+                {
+                    samplePoints_[pointi] =
+                        mesh().faces()[facei].nearestPoint
+                        (
+                            pt,
+                            mesh().points()
+                        ).rawPoint();
+                }
             }
         }
     }
     else
     {
         // if sampleSource_ == cells:
-        //      sampleElements_ : per surface triangle the cell
+        //      sampleElements_ : per surface face, the cell
         //      samplePoints_   : n/a
         // if sampleSource_ == insideCells:
-        //      sampleElements_ : -1 or per surface triangle the cell
+        //      sampleElements_ : -1 or per surface face, the cell
         //      samplePoints_   : n/a
-        // else:
-        //      sampleElements_ : per surface triangle the boundary face
+        // if sampleSource_ == boundaryFaces:
+        //      sampleElements_ : per surface face, the boundary face
         //      samplePoints_   : n/a
 
         sampleElements_.transfer(cellOrFaceLabels);
@@ -432,14 +445,14 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
             forAll(samplePoints_, pointi)
             {
                 meshTools::writeOBJ(str, points()[pointi]);
-                vertI++;
+                ++vertI;
 
                 meshTools::writeOBJ(str, samplePoints_[pointi]);
-                vertI++;
+                ++vertI;
 
                 label elemi = sampleElements_[pointi];
                 meshTools::writeOBJ(str, centres[elemi]);
-                vertI++;
+                ++vertI;
                 str << "l " << vertI-2 << ' ' << vertI-1 << ' ' << vertI << nl;
             }
         }
@@ -449,11 +462,11 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
             forAll(sampleElements_, triI)
             {
                 meshTools::writeOBJ(str, faceCentres()[triI]);
-                vertI++;
+                ++vertI;
 
                 label elemi = sampleElements_[triI];
                 meshTools::writeOBJ(str, centres[elemi]);
-                vertI++;
+                ++vertI;
                 str << "l " << vertI-1 << ' ' << vertI << nl;
             }
         }
@@ -475,7 +488,7 @@ Foam::sampledMeshedSurface::sampledMeshedSurface
 )
 :
     sampledSurface(name, mesh),
-    MeshStorage(),
+    meshedSurface(),
     surfaceName_(surfaceName),
     surface_
     (
@@ -487,7 +500,9 @@ Foam::sampledMeshedSurface::sampledMeshedSurface
     keepIds_(true),
     zoneIds_(),
     sampleElements_(),
-    samplePoints_()
+    samplePoints_(),
+    maxDistanceSqr_(Foam::sqr(GREAT)),
+    defaultValues_()
 {}
 
 
@@ -499,7 +514,7 @@ Foam::sampledMeshedSurface::sampledMeshedSurface
 )
 :
     sampledSurface(name, mesh, dict),
-    MeshStorage(),
+    meshedSurface(),
     surfaceName_
     (
         meshedSurface::findFile
@@ -518,8 +533,25 @@ Foam::sampledMeshedSurface::sampledMeshedSurface
     keepIds_(dict.getOrDefault("keepIds", true)),
     zoneIds_(),
     sampleElements_(),
-    samplePoints_()
+    samplePoints_(),
+    maxDistanceSqr_(Foam::sqr(GREAT)),
+    defaultValues_(dict.subOrEmptyDict("defaultValue"))
 {
+    if (dict.readIfPresent("maxDistance", maxDistanceSqr_))
+    {
+        // Info<< "Limit to maxDistance " << maxDistanceSqr_ << nl;
+        // if (defaultValues_.empty())
+        // {
+        //     Info<< "defaultValues = zero" << nl;
+        // }
+        // else
+        // {
+        //     defaultValues_.writeEntry(Info);
+        // }
+
+        maxDistanceSqr_ = Foam::sqr(maxDistanceSqr_);
+    }
+
     wordRes includePatches;
     dict.readIfPresent("patches", includePatches);
     includePatches.uniq();
@@ -597,7 +629,7 @@ bool Foam::sampledMeshedSurface::expire()
     }
 
     sampledSurface::clearGeom();
-    MeshStorage::clear();
+    meshedSurface::clear();
     zoneIds_.clear();
 
     sampleElements_.clear();
diff --git a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.H b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.H
index bfc4ca4c69..5b61d0ef18 100644
--- a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.H
+++ b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.H
@@ -35,33 +35,35 @@ Description
 
     - 6 different modes:
         - source=cells, interpolate=false:
-            finds per triangle centre the \em nearest cell centre and uses its value
+            finds per surface face centre the \em nearest cell centre
+            and uses its value
         - source=cells, interpolate=true:
-            finds per triangle centre the \em nearest cell centre.
+            finds per surface face 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
             gets a valid location)
 
         - source=insideCells, interpolate=false:
-            finds per triangle centre the cell containing it and uses its value.
-            Trims triangles outside mesh.
+            finds per surface face centre the cell containing it
+            and uses its value.
+            Trims surface faces outside of the 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
+            finds per surface face centre the \em 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
+            finds per surface face centre the \em 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.
+    - since it finds the nearest per surface face, each surface face
+      is guaranteed to be on one processor only.
+      So after stitching the original surface should be complete.
 
     This is often embedded as part of a sampled surfaces function object.
 
@@ -75,6 +77,16 @@ Usage
             type    meshedSurface;
             surface something.obj;
             source  cells;
+
+            //- Max sampling distance
+            maxDistance 0.005;
+
+            //- Fallback for missed sampling in 'cells' mode
+            defaultValue
+            {
+                "p.*"   1e5;
+                T       273.15;
+            }
         }
     }
     \endverbatim
@@ -90,8 +102,11 @@ Usage
         file     | Alternative file name                    | no  |
         fileType | The surface format                       | no  | (extension)
         scale    | Surface scaling factor                   | no  | 0
+        maxDistance | Max search distance                   | no  | GREAT
+        defaultValue | Value beyond max distance (dictionary) | no  | empty
     \endtable
 
+
 SourceFiles
     sampledMeshedSurface.C
     sampledMeshedSurfaceTemplates.C
@@ -102,8 +117,7 @@ SourceFiles
 #define sampledMeshedSurface_H
 
 #include "sampledSurface.H"
-#include "MeshedSurface.H"
-#include "MeshedSurfacesFwd.H"
+#include "MeshedSurfaces.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -123,18 +137,19 @@ class sampledMeshedSurface
     public meshedSurface
 {
 public:
+
         //- Types of sampling regions
-        enum class samplingSource
+        enum class samplingSource : unsigned char
         {
             cells,          //!< Use nearest cell value
-            insideCells,    //!< Face within a cell
-            boundaryFaces   //!< Use boundary face values
+            insideCells,    //!< Surface face within a cell, or trim
+            boundaryFaces   //!< Use nearest boundary face values
         };
 
 private:
 
-        //- Private typedefs for convenience
-        typedef meshedSurface MeshStorage;
+    //- The concrete storage type for faces and points
+    typedef meshedSurface Mesh;
 
 
     // Private Data
@@ -165,6 +180,12 @@ private:
         //- Local points to sample per point
         pointField samplePoints_;
 
+        //- Max search distance squared
+        scalar maxDistanceSqr_;
+
+        //- Out-of-range fallback values (when beyond the max search distance)
+        dictionary defaultValues_;
+
 
     // Private Member Functions
 
@@ -239,13 +260,13 @@ public:
         //- Points of surface
         virtual const pointField& points() const
         {
-            return MeshStorage::points();
+            return Mesh::points();
         }
 
         //- Faces of surface
         virtual const faceList& faces() const
         {
-            return MeshStorage::surfFaces();
+            return Mesh::surfFaces();
         }
 
         //- Per-face zone/region information
@@ -257,31 +278,31 @@ public:
         //- Per-face identifier (eg, element Id)
         virtual const labelList& faceIds() const
         {
-            return MeshStorage::faceIds();
+            return Mesh::faceIds();
         }
 
         //- Face area vectors
         virtual const vectorField& Sf() const
         {
-            return MeshStorage::Sf();
+            return Mesh::Sf();
         }
 
         //- Face area magnitudes
         virtual const scalarField& magSf() const
         {
-            return MeshStorage::magSf();
+            return Mesh::magSf();
         }
 
         //- Face centres
         virtual const vectorField& Cf() const
         {
-            return MeshStorage::Cf();
+            return Mesh::Cf();
         }
 
         //- If element ids/order of the original surface are kept
         virtual bool hasFaceIds() const
         {
-            return keepIds_ && !MeshStorage::faceIds().empty();
+            return keepIds_ && !Mesh::faceIds().empty();
         }
 
         //- Sampling boundary values instead of cell values
diff --git a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceTemplates.C b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceTemplates.C
index dd8cd7c259..6833dd2bda 100644
--- a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceTemplates.C
@@ -37,8 +37,14 @@ Foam::sampledMeshedSurface::sampleOnFaces
     const interpolation<Type>& sampler
 ) const
 {
+    const Type deflt
+    (
+        defaultValues_.getOrDefault<Type>(sampler.psi().name(), Zero)
+    );
+
     const labelList& elements = sampleElements_;
 
+
     if (!onBoundary())
     {
         // Sample cells
@@ -48,7 +54,8 @@ Foam::sampledMeshedSurface::sampleOnFaces
             sampler,
             elements,
             faces(),
-            points()
+            points(),
+            deflt
         );
     }
 
@@ -60,33 +67,33 @@ Foam::sampledMeshedSurface::sampleOnFaces
     auto tvalues = tmp<Field<Type>>::New(elements.size());
     auto& values = tvalues.ref();
 
-    const polyBoundaryMesh& pbm = mesh().boundaryMesh();
-    const label nBnd = mesh().nBoundaryFaces();
-
     // Create flat boundary field
 
-    Field<Type> bVals(nBnd, Zero);
+    const polyBoundaryMesh& pbm = mesh().boundaryMesh();
+
+    Field<Type> bVals(mesh().nBoundaryFaces(), Zero);
 
     const auto& bField = sampler.psi().boundaryField();
 
     forAll(bField, patchi)
     {
-        const label bFacei = (pbm[patchi].start() - mesh().nInternalFaces());
-
-        SubList<Type>
-        (
-            bVals,
-            bField[patchi].size(),
-            bFacei
-        ) = bField[patchi];
+        SubList<Type>(bVals, pbm[patchi].range()) = bField[patchi];
     }
 
-    // Sample in flat boundary field
+    // Sample within the flat boundary field
 
     forAll(elements, i)
     {
         const label bFacei = (elements[i] - mesh().nInternalFaces());
-        values[i] = bVals[bFacei];
+
+        if (bFacei < 0)
+        {
+            values[i] = deflt;
+        }
+        else
+        {
+            values[i] = bVals[bFacei];
+        }
     }
 
     return tvalues;
@@ -100,34 +107,60 @@ Foam::sampledMeshedSurface::sampleOnPoints
     const interpolation<Type>& interpolator
 ) const
 {
-    // One value per vertex
-    auto tvalues = tmp<Field<Type>>::New(sampleElements_.size());
+    const Type deflt
+    (
+        defaultValues_.getOrDefault<Type>(interpolator.psi().name(), Zero)
+    );
+
+    const labelList& elements = sampleElements_;
+
+
+    // One value per vertex.
+    // - sampleElements_ and samplePoints_ have the identical size
+    auto tvalues = tmp<Field<Type>>::New(elements.size());
     auto& values = tvalues.ref();
 
     if (!onBoundary())
     {
         // Sample cells
 
-        forAll(sampleElements_, pointi)
+        forAll(elements, i)
         {
-            values[pointi] = interpolator.interpolate
-            (
-                samplePoints_[pointi],
-                sampleElements_[pointi]
-            );
+            const label celli = elements[i];
+
+            if (celli < 0)
+            {
+                values[i] = deflt;
+            }
+            else
+            {
+                values[i] = interpolator.interpolate
+                (
+                    samplePoints_[i],
+                    celli
+                );
+            }
         }
     }
-    else
+
+
+    //
+    // Sample boundary faces
+    //
+
+    forAll(elements, i)
     {
-        // Sample boundary faces
+        const label facei = elements[i];
 
-        forAll(samplePoints_, pointi)
+        if (facei < 0)
         {
-            const label facei = sampleElements_[pointi];
-
-            values[pointi] = interpolator.interpolate
+            values[i] = deflt;
+        }
+        else
+        {
+            values[i] = interpolator.interpolate
             (
-                samplePoints_[pointi],
+                samplePoints_[i],
                 mesh().faceOwner()[facei],
                 facei
             );
diff --git a/src/sampling/sampledSurface/sampledSurface/sampledSurface.H b/src/sampling/sampledSurface/sampledSurface/sampledSurface.H
index b21df8349e..d135ab0693 100644
--- a/src/sampling/sampledSurface/sampledSurface/sampledSurface.H
+++ b/src/sampling/sampledSurface/sampledSurface/sampledSurface.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -132,14 +132,16 @@ protected:
 
     // Protected Member Functions
 
-        //- General loop for sampling elements to faces
+        //- General loop for sampling elements to faces.
+        //  The defaultValue is used for invalid (negative) elements
         template<class Type>
         static tmp<Field<Type>> sampleOnFaces
         (
             const interpolation<Type>& sampler,
             const labelUList& elements,
             const faceList& fcs,
-            const pointField& pts
+            const pointField& pts,
+            const Type& defaultValue = Type(Zero)
         );
 
 
diff --git a/src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C b/src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C
index 90a5e7f376..2b494d5027 100644
--- a/src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2019 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,7 +37,8 @@ Foam::sampledSurface::sampleOnFaces
     const interpolation<Type>& sampler,
     const labelUList& elements,
     const faceList& fcs,
-    const pointField& pts
+    const pointField& pts,
+    const Type& defaultValue
 )
 {
     const label len = elements.size();
@@ -57,9 +58,16 @@ Foam::sampledSurface::sampleOnFaces
     for (label i=0; i < len; ++i)
     {
         const label celli = elements[i];
-        const point pt = fcs[i].centre(pts);
+        if (celli < 0)
+        {
+            values[i] = defaultValue;
+        }
+        else
+        {
+            const point pt = fcs[i].centre(pts);
 
-        values[i] = sampler.interpolate(pt, celli);
+            values[i] = sampler.interpolate(pt, celli);
+        }
     }
 
     return tvalues;
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling b/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling
index b5f9317229..529a5a6c06 100644
--- a/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling
+++ b/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling
@@ -121,6 +121,44 @@ sampled
 }
 
 
+// Sample and write for meshed surfaces
+sampleMesh
+{
+    type    surfaces;
+    libs    (sampling);
+    log     true;
+
+    writeControl    writeTime;
+    writeInterval   1;
+
+    surfaceFormat   vtk;
+
+    fields      (p rho U T);
+
+    surfaces
+    {
+        // Oversized sampling - for general testing
+        meshed_sample
+        {
+            type        meshedSurface;
+            surface     oversized_sample.obj;
+            source      cells;
+            maxDistance 0.0025;
+            defaultValue
+            {
+                T     273;
+            }
+        }
+
+        meshed_interpolate
+        {
+            $meshed_sample;
+            interpolate true;
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * Calculations  * * * * * * * * * * * * * * * //
 
 massflow

From ab692caf7cb8afe46708c93ac2144574503a519f Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Tue, 24 Nov 2020 17:03:00 +0100
Subject: [PATCH 7/8] ENH: support multiple weights on some field function
 objects (#1930)

- weight fields are combined by multiplication

  - volFieldValue:
    * 0-N scalar fields

  - surfaceFieldValue:
    * 0-N scalar fields
    * 0-1 vector fields

  In some cases this can be used to avoid creating additional
  fields.

      weightFields  (rho U);
  vs.
      derivedFields (rhoU);
      weightField   rhoU;
---
 .../surfaceFieldValue/surfaceFieldValue.C     | 129 +++++++++++++-----
 .../surfaceFieldValue/surfaceFieldValue.H     |   9 +-
 .../surfaceFieldValueTemplates.C              |   2 +-
 .../fieldValues/volFieldValue/volFieldValue.C |  89 ++++++++++--
 .../fieldValues/volFieldValue/volFieldValue.H |   9 +-
 .../volFieldValue/volFieldValueTemplates.C    |   2 +-
 .../rhoSimpleFoam/squareBend/system/sampling  |   2 +
 7 files changed, 189 insertions(+), 53 deletions(-)

diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
index 5f429e1214..bc55e59523 100644
--- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
+++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
@@ -645,9 +645,14 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::writeFileHeader
         writeHeaderValue(os, "Area", totalArea_);
         writeHeaderValue(os, "Scale factor", scaleFactor_);
 
-        if (weightFieldName_ != "none")
+        if (weightFieldNames_.size())
         {
-            writeHeaderValue(os, "Weight field", weightFieldName_);
+            writeHeaderValue
+            (
+                os,
+                "Weight field",
+                flatOutput(weightFieldNames_, FlatOutput::BareComma{})
+            );
         }
 
         writeCommented(os, "Time");
@@ -915,7 +920,7 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::surfaceFieldValue
     needsUpdate_(true),
     writeArea_(false),
     selectionNames_(),
-    weightFieldName_("none"),
+    weightFieldNames_(),
     totalArea_(0),
     nFaces_(0),
     faceId_(),
@@ -949,7 +954,7 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::surfaceFieldValue
     needsUpdate_(true),
     writeArea_(false),
     selectionNames_(),
-    weightFieldName_("none"),
+    weightFieldNames_(),
     totalArea_(0),
     nFaces_(0),
     faceId_(),
@@ -971,7 +976,8 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::read
 
     needsUpdate_ = true;
     writeArea_ = dict.getOrDefault("writeArea", false);
-    weightFieldName_ = "none";
+    weightFieldNames_.clear();
+
     totalArea_ = 0;
     nFaces_ = 0;
     faceId_.clear();
@@ -1054,11 +1060,29 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::read
                 << exit(FatalIOError);
         }
 
-        if (dict.readIfPresent("weightField", weightFieldName_))
+        // Can have "weightFields" or "weightField"
+
+        bool missing = true;
+        if (dict.readIfPresent("weightFields", weightFieldNames_))
         {
-            Info<< "    weight field  = " << weightFieldName_ << nl;
+            missing = false;
         }
         else
+        {
+            weightFieldNames_.resize(1);
+
+            if (dict.readIfPresent("weightField", weightFieldNames_.first()))
+            {
+                missing = false;
+                if ("none" == weightFieldNames_.first())
+                {
+                    // "none" == no weighting
+                    weightFieldNames_.clear();
+                }
+            }
+        }
+
+        if (missing)
         {
             // Suggest possible alternative unweighted operation?
             FatalIOErrorInFunction(dict)
@@ -1069,6 +1093,16 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::read
                 << "or use a different operation."
                 << exit(FatalIOError);
         }
+
+        Info<< "    weight field  = ";
+        if (weightFieldNames_.empty())
+        {
+            Info<< "none" << nl;
+        }
+        else
+        {
+            Info<< flatOutput(weightFieldNames_) << nl;
+        }
     }
 
     // Backwards compatibility for v1612 and older
@@ -1177,46 +1211,79 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::write()
         }
     }
 
-    // Only a few weight types (scalar, vector)
-    if (weightFieldName_ != "none")
+
+    // Check availability and type of weight field
+    // Only support a few weight types:
+    // scalar: 0-N fields
+    // vector: 0-1 fields
+
+    // Default is a zero-size scalar weight field (ie, weight = 1)
+    scalarField scalarWeights;
+    vectorField vectorWeights;
+
+    for (const word& weightName : weightFieldNames_)
     {
-        if (validField<scalar>(weightFieldName_))
+        if (validField<scalar>(weightName))
         {
-            scalarField weightField
-            (
-                getFieldValues<scalar>(weightFieldName_, true)
-            );
+            tmp<scalarField> tfld = getFieldValues<scalar>(weightName, true);
 
-            // Process the fields
-            writeAll(Sf, weightField, points, faces);
+            if (scalarWeights.empty())
+            {
+                scalarWeights = tfld;
+            }
+            else
+            {
+                scalarWeights *= tfld;
+            }
         }
-        else if (validField<vector>(weightFieldName_))
+        else if (validField<vector>(weightName))
         {
-            vectorField weightField
-            (
-                getFieldValues<vector>(weightFieldName_, true)
-            );
+            tmp<vectorField> tfld = getFieldValues<vector>(weightName, true);
 
-            // Process the fields
-            writeAll(Sf, weightField, points, faces);
+            if (vectorWeights.empty())
+            {
+                vectorWeights = tfld;
+            }
+            else
+            {
+                FatalErrorInFunction
+                    << "weightField " << weightName
+                    << " - only one vector weight field allowed. " << nl
+                    << "weights: " << flatOutput(weightFieldNames_) << nl
+                    << abort(FatalError);
+            }
         }
-        else
+        else if (weightName != "none")
         {
+            // Silently ignore "none", flag everything else as an error
+
+            // TBD: treat missing "rho" like incompressible with rho = 1
+            // and/or provided rhoRef value
+
             FatalErrorInFunction
-                << "weightField " << weightFieldName_
-                << " not found or an unsupported type"
+                << "weightField " << weightName
+                << " not found or an unsupported type" << nl
                 << abort(FatalError);
         }
     }
+
+
+    // Process the fields
+    if (vectorWeights.size())
+    {
+        if (scalarWeights.size())
+        {
+            vectorWeights *= scalarWeights;
+        }
+
+        writeAll(Sf, vectorWeights, points, faces);
+    }
     else
     {
-        // Default is a zero-size scalar weight field (ie, weight = 1)
-        scalarField weightField;
-
-        // Process the fields
-        writeAll(Sf, weightField, points, faces);
+        writeAll(Sf, scalarWeights, points, faces);
     }
 
+
     if (operation_ != opNone)
     {
         file() << endl;
diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H
index 5687e9f4a4..1933802913 100644
--- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H
+++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H
@@ -84,7 +84,7 @@ Usage
         names           (<zone-name> <zone-regex>);
 
         postOperation   none;
-        weightField     alpha1;
+        weightFields    (rho U);
         scaleFactor     1.0;
         writeArea       false;
         surfaceFormat   none;
@@ -105,7 +105,8 @@ Usage
       names        | Extended selection                 | word/regex list | no | -
       operation    | Operation type: see below          | word |  yes  | -
       postOperation | Post-operation type: see below    | word |  no   | none
-      weightField  | Name of field to apply weighting   | word |  no   | none
+      weightField   | Name of field to apply weighting   | word |  maybe |
+      weightFields  | Names of fields to apply weighting | wordList | maybe |
       scaleFactor  | Output value scaling factor        | scalar |  no  | 1.0
       writeArea    | Write the surface area             | bool |  no   | false
       surfaceFormat | Output value format               | word <!--
@@ -403,8 +404,8 @@ protected:
         //- Extended selections
         wordRes selectionNames_;
 
-        //- Weight field name - optional
-        word weightFieldName_;
+        //- Weight field name(s) - optional
+        wordList weightFieldNames_;
 
         //- Total area of the surfaceFieldValue
         scalar totalArea_;
diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C
index 53e1dcad9a..8b08f07037 100644
--- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C
+++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C
@@ -118,7 +118,7 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::getFieldValues
     if (mandatory)
     {
         FatalErrorInFunction
-            << "Field " << fieldName << " not found in database"
+            << "Field " << fieldName << " not found in database" << nl
             << abort(FatalError);
     }
 
diff --git a/src/functionObjects/field/fieldValues/volFieldValue/volFieldValue.C b/src/functionObjects/field/fieldValues/volFieldValue/volFieldValue.C
index 9ac7caec75..f18041e68a 100644
--- a/src/functionObjects/field/fieldValues/volFieldValue/volFieldValue.C
+++ b/src/functionObjects/field/fieldValues/volFieldValue/volFieldValue.C
@@ -135,9 +135,15 @@ void Foam::functionObjects::fieldValues::volFieldValue::writeFileHeader
 ) const
 {
     volRegion::writeFileHeader(*this, os);
-    if (weightFieldName_ != "none")
+
+    if (weightFieldNames_.size())
     {
-        writeHeaderValue(os, "Weight field", weightFieldName_);
+        writeHeaderValue
+        (
+            os,
+            "Weight field",
+            flatOutput(weightFieldNames_, FlatOutput::BareComma{})
+        );
     }
 
     writeCommented(os, "Time");
@@ -210,7 +216,7 @@ Foam::functionObjects::fieldValues::volFieldValue::volFieldValue
             true  // Failsafe behaviour
         )
     ),
-    weightFieldName_("none")
+    weightFieldNames_()
 {
     read(dict);
     writeFileHeader(file());
@@ -237,7 +243,7 @@ Foam::functionObjects::fieldValues::volFieldValue::volFieldValue
             true  // Failsafe behaviour
         )
     ),
-    weightFieldName_("none")
+    weightFieldNames_()
 {
     read(dict);
 }
@@ -252,15 +258,33 @@ bool Foam::functionObjects::fieldValues::volFieldValue::read
 {
     fieldValue::read(dict);
 
-    weightFieldName_ = "none";
+    weightFieldNames_.clear();
 
     if (usesWeight())
     {
-        if (dict.readIfPresent("weightField", weightFieldName_))
+        // Can have "weightFields" or "weightField"
+
+        bool missing = true;
+        if (dict.readIfPresent("weightFields", weightFieldNames_))
         {
-            Info<< "    weight field = " << weightFieldName_;
+            missing = false;
         }
         else
+        {
+            weightFieldNames_.resize(1);
+
+            if (dict.readIfPresent("weightField", weightFieldNames_.first()))
+            {
+                missing = false;
+                if ("none" == weightFieldNames_.first())
+                {
+                    // "none" == no weighting
+                    weightFieldNames_.clear();
+                }
+            }
+        }
+
+        if (missing)
         {
             // Suggest possible alternative unweighted operation?
             FatalIOErrorInFunction(dict)
@@ -271,6 +295,16 @@ bool Foam::functionObjects::fieldValues::volFieldValue::read
                 << "or use a different operation."
                 << exit(FatalIOError);
         }
+
+        Info<< "    weight field  = ";
+        if (weightFieldNames_.empty())
+        {
+            Info<< "none" << nl;
+        }
+        else
+        {
+            Info<< flatOutput(weightFieldNames_) << nl;
+        }
     }
 
     Info<< nl << endl;
@@ -297,14 +331,45 @@ bool Foam::functionObjects::fieldValues::volFieldValue::write()
         V = filterField(fieldValue::mesh_.V());
     }
 
-    // Weight field - zero-size means weight = 1
-    scalarField weightField;
-    if (weightFieldName_ != "none")
+    // Check availability and type of weight field
+    // Only support a few weight types:
+    // scalar: 0-N fields
+
+    // Default is a zero-size scalar weight field (ie, weight = 1)
+    scalarField scalarWeights;
+
+    for (const word& weightName : weightFieldNames_)
     {
-        weightField = getFieldValues<scalar>(weightFieldName_, true);
+        if (validField<scalar>(weightName))
+        {
+            tmp<scalarField> tfld = getFieldValues<scalar>(weightName, true);
+
+            if (scalarWeights.empty())
+            {
+                scalarWeights = tfld;
+            }
+            else
+            {
+                scalarWeights *= tfld;
+            }
+        }
+        else if (weightName != "none")
+        {
+            // Silently ignore "none", flag everything else as an error
+
+            // TBD: treat missing "rho" like incompressible with rho = 1
+            // and/or provided rhoRef value
+
+            FatalErrorInFunction
+                << "weightField " << weightName
+                << " not found or an unsupported type" << nl
+                << abort(FatalError);
+        }
     }
 
-    writeAll(V, weightField);
+
+    // Process the fields
+    writeAll(V, scalarWeights);
 
     if (Pstream::master())
     {
diff --git a/src/functionObjects/field/fieldValues/volFieldValue/volFieldValue.H b/src/functionObjects/field/fieldValues/volFieldValue/volFieldValue.H
index 9b5dd4dda7..fd9ab43e01 100644
--- a/src/functionObjects/field/fieldValues/volFieldValue/volFieldValue.H
+++ b/src/functionObjects/field/fieldValues/volFieldValue/volFieldValue.H
@@ -69,7 +69,8 @@ Usage
       name         | Name for regionType                | word |  yes  | -
       operation    | Operation type: see below          | word |  yes  | -
       postOperation | Post-operation type: see below    | word |  no   | none
-      weightField  | Name of field to apply weighting   | word |  no   | none
+      weightField   | Name of field to apply weighting   | word |  maybe |
+      weightFields  | Names of fields to apply weighting | wordList | maybe |
     \endtable
 
     The inherited entries are elaborated in:
@@ -214,8 +215,8 @@ protected:
         //- Optional post-evaluation operation
         postOperationType postOperation_;
 
-        //- Weight field name - only used for weighted modes
-        word weightFieldName_;
+        //- Weight field name(s) - optional
+        wordList weightFieldNames_;
 
 
     // Protected Member Functions
@@ -233,7 +234,7 @@ protected:
         //  Checks for availability on any processor.
         inline bool canWeight(const scalarField& weightField) const;
 
-        //- True if the field name is valid (exists, and a supported type)
+        //- Return true if the field name is valid
         template<class Type>
         bool validField(const word& fieldName) const;
 
diff --git a/src/functionObjects/field/fieldValues/volFieldValue/volFieldValueTemplates.C b/src/functionObjects/field/fieldValues/volFieldValue/volFieldValueTemplates.C
index 49b93125fe..2e419e9060 100644
--- a/src/functionObjects/field/fieldValues/volFieldValue/volFieldValueTemplates.C
+++ b/src/functionObjects/field/fieldValues/volFieldValue/volFieldValueTemplates.C
@@ -71,7 +71,7 @@ Foam::functionObjects::fieldValues::volFieldValue::getFieldValues
     if (mandatory)
     {
         FatalErrorInFunction
-            << "Field " << fieldName << " not found in database"
+            << "Field " << fieldName << " not found in database" << nl
             << abort(FatalError);
     }
 
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling b/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling
index 529a5a6c06..052d6a35a8 100644
--- a/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling
+++ b/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling
@@ -182,6 +182,7 @@ areaAverage
 
     operation       weightedAreaAverage;
     weightField     rhoU;
+    weightFields    ( rho U none );  // 2012 and later
     fields          ( p pTotal );
 }
 
@@ -194,6 +195,7 @@ areaIntegrate
 
     operation       weightedAreaIntegrate;
     weightField     rhoU;
+    weightFields    ( rho U );  // 2012 and later
     fields          ( T );
 }
 

From 029b1373a37bab6e5895cb115cb2b32c6a065f61 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Mon, 30 Nov 2020 13:54:41 +0100
Subject: [PATCH 8/8] DOC: document writeFile updateHeader flag (fixes #1942)

---
 .../db/functionObjects/writeFile/writeFile.H      | 15 +++++++++++----
 1 file changed, 11 insertions(+), 4 deletions(-)

diff --git a/src/OpenFOAM/db/functionObjects/writeFile/writeFile.H b/src/OpenFOAM/db/functionObjects/writeFile/writeFile.H
index 6a5f5f3a88..11df4d403c 100644
--- a/src/OpenFOAM/db/functionObjects/writeFile/writeFile.H
+++ b/src/OpenFOAM/db/functionObjects/writeFile/writeFile.H
@@ -47,12 +47,19 @@ Usage
 
     where the entries mean:
     \table
-      Property        | Description                      | Type | Req'd | Dflt
-      writePrecision  | Number of decimal points | label | no  | \<system dflt\>
-      writeToFile     | Flag to produce text file output | bool | no    | true
-      useUserTime | Flag to use user time, e.g. degrees  | bool | no    | true
+      Property        | Description                      | Type | Reqd | Dflt
+      writePrecision  | Number of decimal points | int  | no  | \<system dflt\>
+      writeToFile     | Produce text file output?        | bool | no  | true
+      useUserTime     | Use user time (e.g. degrees)?    | bool | no  | true
+      updateHeader    | Update header on mesh changes?   | bool | no  | true
     \endtable
 
+Note
+    The file header is normally updated whenver the mesh points or
+    topology changes. In some cases, the function object is actually
+    unaffected by these changes.
+    Use the \c updateHeader flag to override the default behaviour.
+
 See also
     - Foam::functionObject
     - Foam::functionObjects::logFiles