functionObjects: surfaceFieldValue, volFieldValue: Various improvements

A number of changes have been made to the surfaceFieldValue and volFieldValue function objects to improve their usability and performance, and to extend them so that similar duplicate functionality elsewhere in OpenFOAM can be removed. Weighted operations have been removed. Weighting for averages and sums is now triggered simply by the existence of the "weightField" or "weightFields" entry. Multiple weight fields are now supported in both functions. The distinction between oriented and non-oriented fields has been removed from surfaceFieldValue. There is now just a single list of fields which are operated on. Instead of oriented fields, an "orientedSum" operation has been added, which should be used for flowRate calculations and other similar operations on fluxes. Operations minMag and maxMag have been added to both functions, to calculate the minimum and maximum field magnitudes respectively. The min and max operations are performed component-wise, as was the case previously. In volFieldValue, minMag and maxMag (and min and mag operations when applied to scalar fields) will report the location, cell and processor of the maximum or minimum value. There is also a "writeLocation" option which if set will write this location information into the output file. The fieldMinMax function has been made obsolete by this change, and has therefore been removed. surfaceFieldValue now operates in parallel without accumulating the entire surface on the master processor for calculation of the operation. Collecting the entire surface on the master processor is now only done if the surface itself is to be written out.
2021-07-08 10:51:34 +01:00
parent e00316c7be
commit 056cc20f34
28 changed files with 781 additions and 1302 deletions
--- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
+++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
@ -63,23 +63,22 @@ template<>
 const char* Foam::NamedEnum
 <
    Foam::functionObjects::fieldValues::surfaceFieldValue::operationType,
-    17
+    16
 >::names[] =
 {
    "none",
    "sum",
-    "weightedSum",
    "sumMag",
    "sumDirection",
    "sumDirectionBalance",
+    "orientedSum",
    "average",
-    "weightedAverage",
    "areaAverage",
-    "weightedAreaAverage",
    "areaIntegrate",
-    "weightedAreaIntegrate",
    "min",
    "max",
+    "minMag",
+    "maxMag",
    "CoV",
    "areaNormalAverage",
    "areaNormalIntegrate"
@ -94,7 +93,7 @@ const Foam::NamedEnum
 const Foam::NamedEnum
 <
    Foam::functionObjects::fieldValues::surfaceFieldValue::operationType,
-    17
+    16
 > Foam::functionObjects::fieldValues::surfaceFieldValue::operationTypeNames_;


@ -366,7 +365,7 @@ combineSurfaceGeometry
    pointField& points
 ) const
 {
-    if (surfacePtr_.valid())
+    if (regionType_ == regionTypes::sampledSurface)
    {
        const sampledSurface& s = surfacePtr_();

@ -402,18 +401,14 @@ combineSurfaceGeometry
 Foam::scalar
 Foam::functionObjects::fieldValues::surfaceFieldValue::totalArea() const
 {
-    scalar totalArea;
-
-    if (surfacePtr_.valid())
+    if (regionType_ == regionTypes::sampledSurface)
    {
-        totalArea = gSum(surfacePtr_().magSf());
+        return gSum(surfacePtr_().magSf());
    }
    else
    {
-        totalArea = gSum(filterField(mesh_.magSf(), false));
+        return gSum(filterField(mesh_.magSf()));
    }
-
-    return totalArea;
 }


@ -462,7 +457,7 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::initialise
            << "    Region has no faces" << exit(FatalError);
    }

-    if (surfacePtr_.valid())
+    if (regionType_ == regionTypes::sampledSurface)
    {
        surfacePtr_().update();
    }
@ -475,42 +470,18 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::initialise
        << "    total area   = " << totalArea_
        << nl;

-    if (dict.readIfPresent("orientedWeightField", weightFieldName_))
+    if (dict.readIfPresent("weightFields", weightFieldNames_))
    {
-        Info<< "    orientedWeightField = " << weightFieldName_ << nl;
-        orientWeightField_ = true;
-
-        if (regionType_ == regionTypes::sampledSurface)
-        {
-            FatalIOErrorInFunction(dict)
-                << "Cannot use orientedWeightField for a sampledSurface"
-                << exit(FatalIOError);
-        }
-
-        if (dict.found("weightField"))
-        {
-            FatalIOErrorInFunction(dict)
-                << "Either provide weightField or orientedWeightField"
-                << exit(FatalIOError);
-        }
+        Info<< name() << " " << operationTypeNames_[operation_]
+            << " weight fields " << weightFieldNames_;
    }
-    else if (dict.readIfPresent("weightField", weightFieldName_))
+    else if (dict.found("weightField"))
    {
-        Info<< "    weightField = " << weightFieldName_ << nl;
+        weightFieldNames_.setSize(1);
+        dict.lookup("weightField") >> weightFieldNames_[0];

-        if (regionType_ == regionTypes::sampledSurface)
-        {
-            FatalIOErrorInFunction(dict)
-                << "Cannot use weightField for a sampledSurface"
-                << exit(FatalIOError);
-        }
-    }
-
-    List<word> orientedFields;
-    if (dict.readIfPresent("orientedFields", orientedFields))
-    {
-        orientedFieldsStart_ = fields_.size();
-        fields_.append(orientedFields);
+        Info<< name() << " " << operationTypeNames_[operation_]
+            << " weight field " << weightFieldNames_[0];
    }

    if (dict.readIfPresent("scaleFactor", scaleFactor_))
@ -567,8 +538,9 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::writeFileHeader
 bool Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
 (
    const Field<scalar>& values,
+    const scalarField& signs,
+    const scalarField& weights,
    const vectorField& Sf,
-    const scalarField& weightField,
    scalar& result
 ) const
 {
@ -576,21 +548,28 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
    {
        case operationType::sumDirection:
        {
-            vector n(dict_.lookup("direction"));
-            result = sum(pos0(values*(Sf & n))*mag(values));
+            const vector n(dict_.lookup("direction"));
+            result = gSum(weights*pos0(values*(Sf & n))*mag(values));
            return true;
        }
        case operationType::sumDirectionBalance:
        {
-            vector n(dict_.lookup("direction"));
+            const vector n(dict_.lookup("direction"));
            const scalarField nv(values*(Sf & n));
-            result = sum(pos0(nv)*mag(values) - neg(nv)*mag(values));
+            result = gSum(weights*(pos0(nv)*mag(values) - neg(nv)*mag(values)));
            return true;
        }
        default:
        {
            // Fall through to same-type operations
-            return processValuesTypeType(values, Sf, weightField, result);
+            return processValuesTypeType
+            (
+                values,
+                signs,
+                weights,
+                Sf,
+                result
+            );
        }
    }
 }
@ -599,8 +578,9 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
 bool Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
 (
    const Field<vector>& values,
+    const scalarField& signs,
+    const scalarField& weights,
    const vectorField& Sf,
-    const scalarField& weightField,
    scalar& result
 ) const
 {
@ -608,12 +588,12 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
    {
        case operationType::areaNormalAverage:
        {
-            result = sum(values & Sf)/sum(mag(Sf));
+            result = gSum(weights*values & Sf)/gSum(mag(weights*Sf));
            return true;
        }
        case operationType::areaNormalIntegrate:
        {
-            result = sum(values & Sf);
+            result = gSum(weights*values & Sf);
            return true;
        }
        default:
@ -627,8 +607,9 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
 bool Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
 (
    const Field<vector>& values,
+    const scalarField& signs,
+    const scalarField& weights,
    const vectorField& Sf,
-    const scalarField& weightField,
    vector& result
 ) const
 {
@ -636,24 +617,29 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
    {
        case operationType::sumDirection:
        {
-            vector n(dict_.lookup("direction"));
-            n /= mag(n) + rootVSmall;
+            const vector n = normalised(dict_.lookup<vector>("direction"));
            const scalarField nv(n & values);
-            result = sum(pos0(nv)*n*(nv));
+            result = gSum(weights*pos0(nv)*n*(nv));
            return true;
        }
        case operationType::sumDirectionBalance:
        {
-            vector n(dict_.lookup("direction"));
-            n /= mag(n) + rootVSmall;
+            const vector n = normalised(dict_.lookup<vector>("direction"));
            const scalarField nv(n & values);
-            result = sum(pos0(nv)*n*(nv));
+            result = gSum(weights*pos0(nv)*n*(nv));
            return true;
        }
        default:
        {
            // Fall through to same-type operations
-            return processValuesTypeType(values, Sf, weightField, result);
+            return processValuesTypeType
+            (
+                values,
+                signs,
+                weights,
+                Sf,
+                result
+            );
        }
    }
 }
@ -672,10 +658,8 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::surfaceFieldValue
    surfaceWriterPtr_(nullptr),
    regionType_(regionTypeNames_.read(dict.lookup("regionType"))),
    operation_(operationTypeNames_.read(dict.lookup("operation"))),
-    weightFieldName_("none"),
-    orientWeightField_(false),
-    orientedFieldsStart_(labelMax),
-    scaleFactor_(1.0),
+    weightFieldNames_(),
+    scaleFactor_(1),
    writeArea_(dict.lookupOrDefault("writeArea", false)),
    nFaces_(0),
    faceId_(),
@ -696,10 +680,8 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::surfaceFieldValue
    surfaceWriterPtr_(nullptr),
    regionType_(regionTypeNames_.read(dict.lookup("regionType"))),
    operation_(operationTypeNames_.read(dict.lookup("operation"))),
-    weightFieldName_("none"),
-    orientWeightField_(false),
-    orientedFieldsStart_(labelMax),
-    scaleFactor_(1.0),
+    weightFieldNames_(),
+    scaleFactor_(1),
    writeArea_(dict.lookupOrDefault("writeArea", false)),
    nFaces_(0),
    faceId_(),
@ -737,7 +719,7 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::write()
        fieldValue::write();
    }

-    if (surfacePtr_.valid())
+    if (regionType_ == regionTypes::sampledSurface)
    {
        surfacePtr_().update();
    }
@ -758,12 +740,12 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::write()
    }

    // Write the surface geometry
-    if (surfaceWriterPtr_.valid())
+    if (writeFields_)
    {
        faceList faces;
        pointField points;

-        if (surfacePtr_.valid())
+        if (regionType_ == regionTypes::sampledSurface)
        {
            combineSurfaceGeometry(faces, points);
        }
@ -784,21 +766,24 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::write()
        }
    }

-    // Construct weight field. Note: zero size means weight = 1
-    scalarField weightField;
-    if (weightFieldName_ != "none")
+    // Construct the sign and weight fields and the surface normals
+    const scalarField signs
+    (
+        regionType_ == regionTypes::sampledSurface
+      ? scalarField(surfacePtr_().Sf().size(), 1)
+      : List<scalar>(faceSign_)
+    );
+    scalarField weights(signs.size(), 1);
+    forAll(weightFieldNames_, i)
    {
-        weightField =
-            getFieldValues<scalar>
-            (
-                weightFieldName_,
-                true,
-                orientWeightField_
-            );
+        weights *= getFieldValues<scalar>(weightFieldNames_[i]);
    }
-
-    // Combine onto master
-    combineFields(weightField);
+    const vectorField Sf
+    (
+        regionType_ == regionTypes::sampledSurface
+      ? surfacePtr_().Sf()
+      : (signs*filterField(mesh_.Sf()))()
+    );

    // Process the fields
    forAll(fields_, i)
@ -806,12 +791,18 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::write()
        const word& fieldName = fields_[i];
        bool ok = false;

-        bool orient = i >= orientedFieldsStart_;
-        ok = ok || writeValues<scalar>(fieldName, weightField, orient);
-        ok = ok || writeValues<vector>(fieldName, weightField, orient);
-        ok = ok || writeValues<sphericalTensor>(fieldName, weightField, orient);
-        ok = ok || writeValues<symmTensor>(fieldName, weightField, orient);
-        ok = ok || writeValues<tensor>(fieldName, weightField, orient);
+        #define writeValuesFieldType(fieldType, none)                          \
+            ok =                                                               \
+                ok                                                             \
+             || writeValues<fieldType>                                         \
+                (                                                              \
+                    fieldName,                                                 \
+                    signs,                                                     \
+                    weights,                                                   \
+                    Sf                                                         \
+                );
+        FOR_ALL_FIELD_TYPES(writeValuesFieldType);
+        #undef writeValuesFieldType

        if (!ok)
        {
--- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H
+++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H
@ -32,7 +32,7 @@ Description
    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 'orientedSum' operator to the flux field (typically \c phi)

    Examples of function object specification:
    \verbatim
@ -96,10 +96,9 @@ Usage
        name         | name of face regionType if required  | no |
        operation    | operation to perform    | yes         |
        weightField  | name of field to apply weighting | no |
-        orientedWeightField  | name of oriented field to apply weighting | no |
+        weightFields | Names of fields to apply weighting | no |
        scaleFactor  | scale factor            | no          | 1
        fields       | list of fields to operate on | yes    |
-        orientedFields | list of oriented fields to operate on | no |
    \endtable

    Where \c regionType is defined by
@ -113,41 +112,38 @@ Usage
    \plaintable
       none          | no operation
       sum           | sum
-       weightedSum   | weighted sum
       sumMag        | sum of component magnitudes
       sumDirection  | sum values which are positive in given direction
       sumDirectionBalance | sum of balance of values in given direction
+       orientedSum   | sum with face orientations
       average       | ensemble average
-       weightedAverage | weighted average
       areaAverage   | area weighted average
-       weightedAreaAverage | weighted area average
       areaIntegrate | area integral
-       weightedAreaIntegrate | weighted area integral
       min           | minimum
       max           | maximum
+       minMag        | minimum magnitude
+       maxMag        | maximum magnitude
       CoV           | coefficient of variation: standard deviation/mean
-       areaNormalAverage| area weighted average in face normal direction
-       areaNormalIntegrate | area weighted integral in face normal directon
+       areaNormalAverage | area weighted average in face normal direction
+       areaNormalIntegrate | area weighted integral in face normal direction
    \endplaintable

 Note
-    - The values reported by the areaNormalAverage and 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 on empty patches get ignored.
+    - The `oriented' operations will flip the sign of the field so that all the
+      normals point in a consistent direction. This is only of relevance when
+      summing mesh-oriented fields, such as the flux, on faceZones.
+    - If the field is a volField then a \c faceZone can only consist of
+      boundary faces, because only these faces have a value associated with
+      them. No cell-to-face interpolation is performed.
+    - If the field is a surfaceField then the region cannot be a \c
+      sampledSurface
+    - If a sampledSurface has interpolation set to false then the surface
+      face values will be taken directly from the cell that contains the
+      surface face centre
+    - If a \c sampledSurface has interpolation set to true then the field
+      will be interpolated to the vertices, then averaged onto the surface
      faces
-    - the `oriented' entries relate to mesh-oriented fields, such as the
-      flux, phi.  These fields will be oriented according to the face normals.
-    - using \c sampledSurface:
-        - not available for surface fields
-        - if interpolate=true they use \c interpolationCellPoint
-          otherwise they use cell values
-        - each triangle in \c sampledSurface is logically only in one cell
-          so interpolation will be wrong when triangles are larger than
-          cells.  This can only happen for sampling on a \c triSurfaceMesh
-        - take care when using isoSurfaces - these might have duplicate
-          triangles and so integration might be wrong

 See also
    Foam::fieldValues
@ -208,25 +204,24 @@ public:
        {
            none,
            sum,
-            weightedSum,
            sumMag,
            sumDirection,
            sumDirectionBalance,
+            orientedSum,
            average,
-            weightedAverage,
            areaAverage,
-            weightedAreaAverage,
            areaIntegrate,
-            weightedAreaIntegrate,
            min,
            max,
+            minMag,
+            maxMag,
            CoV,
            areaNormalAverage,
            areaNormalIntegrate
        };

        //- Operation type names
-        static const NamedEnum<operationType, 17> operationTypeNames_;
+        static const NamedEnum<operationType, 16> operationTypeNames_;


 private:
@ -273,14 +268,8 @@ protected:
        //- Operation to apply to values
        operationType operation_;

-        //- Weight field name - optional
-        word weightFieldName_;
-
-        //- Flag to indicate if flipMap should be applied to the weight field
-        bool orientWeightField_;
-
-        //- Start index of fields that require application of flipMap
-        label orientedFieldsStart_;
+        //- Weight field names - optional
+        wordList weightFieldNames_;

        //- Scale factor - optional
        scalar scaleFactor_;
@ -325,12 +314,7 @@ protected:

        //- Return field values by looking up field name
        template<class Type>
-        tmp<Field<Type>> getFieldValues
-        (
-            const word& fieldName,
-            const bool mustGet = false,
-            const bool applyOrientation = false
-        ) const;
+        tmp<Field<Type>> getFieldValues(const word& fieldName) const;

        //- Apply the operation to the values, and return true if successful.
        //  Does nothing unless overloaded below.
@ -338,8 +322,9 @@ protected:
        bool processValues
        (
            const Field<Type>& values,
+            const scalarField& signs,
+            const scalarField& weights,
            const vectorField& Sf,
-            const scalarField& weightField,
            ResultType& result
        ) const;

@ -349,19 +334,32 @@ protected:
        bool processValues
        (
            const Field<Type>& values,
+            const scalarField& signs,
+            const scalarField& weights,
            const vectorField& Sf,
-            const scalarField& weightField,
            Type& result
        ) const;

+        //- Apply Type -> scalar operation to the values
+        template<class Type>
+        bool processValues
+        (
+            const Field<Type>& values,
+            const scalarField& signs,
+            const scalarField& weights,
+            const vectorField& Sf,
+            scalar& result
+        ) const;
+
        //- Apply scalar -> scalar operation to the values. Tries to apply
        //  scalar -> scalar specific operations, otherwise calls
        //  processValuesTypeType.
        bool processValues
        (
            const Field<scalar>& values,
+            const scalarField& signs,
+            const scalarField& weights,
            const vectorField& Sf,
-            const scalarField& weightField,
            scalar& result
        ) const;

@ -369,8 +367,9 @@ protected:
        bool processValues
        (
            const Field<vector>& values,
+            const scalarField& signs,
+            const scalarField& weights,
            const vectorField& Sf,
-            const scalarField& weightField,
            scalar& result
        ) const;

@ -380,8 +379,9 @@ protected:
        bool processValues
        (
            const Field<vector>& values,
+            const scalarField& signs,
+            const scalarField& weights,
            const vectorField& Sf,
-            const scalarField& weightField,
            vector& result
        ) const;

@ -390,8 +390,9 @@ protected:
        bool processValuesTypeType
        (
            const Field<Type>& values,
+            const scalarField& signs,
+            const scalarField& weights,
            const vectorField& Sf,
-            const scalarField& weightField,
            Type& result
        ) const;

@ -450,8 +451,9 @@ public:
        bool writeValues
        (
            const word& fieldName,
-            const scalarField& weightField,
-            const bool orient
+            const scalarField& signs,
+            const scalarField& weights,
+            const vectorField& Sf
        );

        //- Templated helper function to output field values
@ -459,8 +461,9 @@ public:
        bool writeValues
        (
            const word& fieldName,
-            const scalarField& weightField,
            const Field<Type>& values,
+            const scalarField& signs,
+            const scalarField& weights,
            const vectorField& Sf
        );

@ -468,16 +471,14 @@ public:
        template<class Type>
        tmp<Field<Type>> filterField
        (
-            const GeometricField<Type, fvsPatchField, surfaceMesh>& field,
-            const bool applyOrientation
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& field
        ) const;

        //- Filter a volume field according to faceIds
        template<class Type>
        tmp<Field<Type>> filterField
        (
-            const GeometricField<Type, fvPatchField, volMesh>& field,
-            const bool applyOrientation
+            const GeometricField<Type, fvPatchField, volMesh>& field
        ) const;

        //- Read from dictionary
--- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C
+++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C
@ -62,42 +62,29 @@ template<class Type>
 Foam::tmp<Foam::Field<Type>>
 Foam::functionObjects::fieldValues::surfaceFieldValue::getFieldValues
 (
-    const word& fieldName,
-    const bool mustGet,
-    const bool applyOrientation
+    const word& fieldName
 ) const
 {
-    typedef GeometricField<Type, fvsPatchField, surfaceMesh> sf;
    typedef GeometricField<Type, fvPatchField, volMesh> vf;
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> sf;

-    if
-    (
-        regionType_ != regionTypes::sampledSurface
-     && obr_.foundObject<sf>(fieldName)
-    )
+    if (regionType_ == regionTypes::sampledSurface)
    {
-        return filterField(obr_.lookupObject<sf>(fieldName), applyOrientation);
-    }
-    else if (obr_.foundObject<vf>(fieldName))
-    {
-        const vf& fld = obr_.lookupObject<vf>(fieldName);
-
-        if (surfacePtr_.valid())
+        if (obr_.foundObject<vf>(fieldName))
        {
+            const vf& fld = obr_.lookupObject<vf>(fieldName);
+
            if (surfacePtr_().interpolate())
            {
+                // Interpolate the field to the surface points
                const interpolationCellPoint<Type> interp(fld);
                tmp<Field<Type>> tintFld(surfacePtr_().interpolate(interp));
                const Field<Type>& intFld = tintFld();

-                // Average
+                // Average the interpolated field onto the surface faces
                const faceList& faces = surfacePtr_().faces();
-                tmp<Field<Type>> tavg
-                (
-                    new Field<Type>(faces.size(), Zero)
-                );
+                tmp<Field<Type>> tavg(new Field<Type>(faces.size(), Zero));
                Field<Type>& avg = tavg.ref();
-
                forAll(faces, facei)
                {
                    const face& f = faces[facei];
@ -115,20 +102,34 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::getFieldValues
                return surfacePtr_().sample(fld);
            }
        }
-        else
+        else if (obr_.foundObject<sf>(fieldName))
        {
-            return filterField(fld, applyOrientation);
+            FatalErrorInFunction
+                << "Surface field " << fieldName
+                << " cannot be sampled onto surface " << surfacePtr_().name()
+                << ". Only vol fields can be sampled onto surfaces."
+                << abort(FatalError);
+        }
+    }
+    else
+    {
+        if (obr_.foundObject<vf>(fieldName))
+        {
+            const vf& fld = obr_.lookupObject<vf>(fieldName);
+            return filterField(fld);
+        }
+        else if (obr_.foundObject<sf>(fieldName))
+        {
+            const sf& fld = obr_.lookupObject<sf>(fieldName);
+            return filterField(fld);
        }
    }

-    if (mustGet)
-    {
-        FatalErrorInFunction
-            << "Field " << fieldName << " not found in database"
-            << abort(FatalError);
-    }
+    FatalErrorInFunction
+        << "Field " << fieldName << " not found in database"
+        << abort(FatalError);

-    return tmp<Field<Type>>(new Field<Type>(0));
+    return tmp<Field<Type>>(nullptr);
 }


@ -136,8 +137,9 @@ template<class Type, class ResultType>
 bool Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
 (
    const Field<Type>& values,
+    const scalarField& signs,
+    const scalarField& weights,
    const vectorField& Sf,
-    const scalarField& weightField,
    ResultType& result
 ) const
 {
@ -149,12 +151,43 @@ template<class Type>
 bool Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
 (
    const Field<Type>& values,
+    const scalarField& signs,
+    const scalarField& weights,
    const vectorField& Sf,
-    const scalarField& weightField,
    Type& result
 ) const
 {
-    return processValuesTypeType(values, Sf, weightField, result);
+    return processValuesTypeType(values, signs, weights, Sf, result);
+}
+
+
+template<class Type>
+bool Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
+(
+    const Field<Type>& values,
+    const scalarField& signs,
+    const scalarField& weights,
+    const vectorField& Sf,
+    scalar& result
+) const
+{
+    switch (operation_)
+    {
+        case operationType::minMag:
+        {
+            result = gMin(mag(values));
+            return true;
+        }
+        case operationType::maxMag:
+        {
+            result = gMax(mag(values));
+            return true;
+        }
+        default:
+        {
+            return false;
+        }
+    }
 }


@ -163,8 +196,9 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::
 processValuesTypeType
 (
    const Field<Type>& values,
+    const scalarField& signs,
+    const scalarField& weights,
    const vectorField& Sf,
-    const scalarField& weightField,
    Type& result
 ) const
 {
@ -172,104 +206,55 @@ processValuesTypeType
    {
        case operationType::sum:
        {
-            result = sum(values);
-            return true;
-        }
-        case operationType::weightedSum:
-        {
-            if (weightField.size())
-            {
-                result = sum(weightField*values);
-            }
-            else
-            {
-                result = sum(values);
-            }
+            result = gSum(weights*values);
            return true;
        }
        case operationType::sumMag:
        {
-            result = sum(cmptMag(values));
+            result = gSum(weights*cmptMag(values));
+            return true;
+        }
+        case operationType::orientedSum:
+        {
+            result = gSum(signs*weights*values);
            return true;
        }
        case operationType::average:
        {
-            result = sum(values)/values.size();
-            return true;
-        }
-        case operationType::weightedAverage:
-        {
-            if (weightField.size())
-            {
-                result =
-                    sum(weightField*values)
-                   /stabilise(sum(weightField), vSmall);
-            }
-            else
-            {
-                result = sum(values)/values.size();
-            }
+            result =
+                gSum(weights*values)
+               /stabilise(gSum(weights), vSmall);
            return true;
        }
        case operationType::areaAverage:
        {
            const scalarField magSf(mag(Sf));
-
-            result = sum(magSf*values)/sum(magSf);
-            return true;
-        }
-        case operationType::weightedAreaAverage:
-        {
-            const scalarField magSf(mag(Sf));
-
-            if (weightField.size())
-            {
-                result =
-                    sum(weightField*magSf*values)
-                   /stabilise(sum(magSf*weightField), vSmall);
-            }
-            else
-            {
-                result = sum(magSf*values)/sum(magSf);
-            }
+            result =
+                gSum(weights*magSf*values)
+               /stabilise(gSum(weights*magSf), vSmall);
            return true;
        }
        case operationType::areaIntegrate:
        {
            const scalarField magSf(mag(Sf));
-
-            result = sum(magSf*values);
-            return true;
-        }
-        case operationType::weightedAreaIntegrate:
-        {
-            const scalarField magSf(mag(Sf));
-
-            if (weightField.size())
-            {
-                result = sum(weightField*magSf*values);
-            }
-            else
-            {
-                result = sum(magSf*values);
-            }
+            result = gSum(weights*magSf*values);
            return true;
        }
        case operationType::min:
        {
-            result = min(values);
+            result = gMin(values);
            return true;
        }
        case operationType::max:
        {
-            result = max(values);
+            result = gMax(values);
            return true;
        }
        case operationType::CoV:
        {
            const scalarField magSf(mag(Sf));

-            Type meanValue = sum(values*magSf)/sum(magSf);
+            Type meanValue = gSum(values*magSf)/gSum(magSf);

            const label nComp = pTraits<Type>::nComponents;

@ -279,7 +264,7 @@ processValuesTypeType
                scalar mean = component(meanValue, d);
                scalar& res = setComponent(result, d);

-                res = sqrt(sum(magSf*sqr(vals - mean))/sum(magSf))/mean;
+                res = sqrt(gSum(magSf*sqr(vals - mean))/gSum(magSf))/mean;
            }

            return true;
@ -302,39 +287,25 @@ template<class Type>
 bool Foam::functionObjects::fieldValues::surfaceFieldValue::writeValues
 (
    const word& fieldName,
-    const scalarField& weightField,
-    const bool orient
+    const scalarField& signs,
+    const scalarField& weights,
+    const vectorField& Sf
 )
 {
    const bool ok = validField<Type>(fieldName);

    if (ok)
    {
-        Field<Type> values(getFieldValues<Type>(fieldName, true, orient));
-
-        vectorField Sf;
-        if (surfacePtr_.valid())
-        {
-            // Get oriented Sf
-            Sf = surfacePtr_().Sf();
-        }
-        else
-        {
-            // Get oriented Sf
-            Sf = filterField(mesh_.Sf(), true);
-        }
-
-        // Combine onto master
-        combineFields(values);
-        combineFields(Sf);
+        // Get the values
+        Field<Type> values(getFieldValues<Type>(fieldName));

        // Write raw values on surface if specified
-        if (surfaceWriterPtr_.valid())
+        if (writeFields_)
        {
            faceList faces;
            pointField points;

-            if (surfacePtr_.valid())
+            if (regionType_ == regionTypes::sampledSurface)
            {
                combineSurfaceGeometry(faces, points);
            }
@ -343,6 +314,9 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::writeValues
                combineMeshGeometry(faces, points);
            }

+            Field<Type> writeValues(weights*values);
+            combineFields(writeValues);
+
            if (Pstream::master())
            {
                surfaceWriterPtr_->write
@ -352,39 +326,41 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::writeValues
                    points,
                    faces,
                    fieldName,
-                    values,
+                    writeValues,
                    false
                );
            }
        }

+        // Do the operation
        if (operation_ != operationType::none)
        {
            // Apply scale factor
            values *= scaleFactor_;

-            if (Pstream::master())
+            bool ok = false;
+
+            #define writeValuesFieldType(fieldType, none)                      \
+                ok =                                                           \
+                    ok                                                         \
+                 || writeValues<Type, fieldType>                               \
+                    (                                                          \
+                        fieldName,                                             \
+                        values,                                                \
+                        signs,                                                 \
+                        weights,                                               \
+                        Sf                                                     \
+                    );
+            FOR_ALL_FIELD_TYPES(writeValuesFieldType);
+            #undef writeValuesFieldType
+
+            if (!ok)
            {
-                if
-                (
-                    !writeValues<Type, scalar>
-                    (fieldName, weightField, values, Sf)
-                 && !writeValues<Type, vector>
-                    (fieldName, weightField, values, Sf)
-                 && !writeValues<Type, sphericalTensor>
-                    (fieldName, weightField, values, Sf)
-                 && !writeValues<Type, symmTensor>
-                    (fieldName, weightField, values, Sf)
-                 && !writeValues<Type, tensor>
-                    (fieldName, weightField, values, Sf)
-                )
-                {
-                    FatalErrorInFunction
-                        << "Operation " << operationTypeNames_[operation_]
-                        << " not available for values of type "
-                        << pTraits<Type>::typeName
-                        << exit(FatalError);
-                }
+                FatalErrorInFunction
+                    << "Operation " << operationTypeNames_[operation_]
+                    << " not available for values of type "
+                    << pTraits<Type>::typeName
+                    << exit(FatalError);
            }
        }
    }
@ -397,23 +373,27 @@ template<class Type, class ResultType>
 bool Foam::functionObjects::fieldValues::surfaceFieldValue::writeValues
 (
    const word& fieldName,
-    const scalarField& weightField,
    const Field<Type>& values,
+    const scalarField& signs,
+    const scalarField& weights,
    const vectorField& Sf
 )
 {
    ResultType result;

-    if (processValues(values, Sf, weightField, result))
+    if (processValues(values, signs, weights, Sf, result))
    {
        // Add to result dictionary, over-writing any previous entry
        resultDict_.add(fieldName, result, true);

-        file() << tab << result;
+        if (Pstream::master())
+        {
+            file() << tab << result;

-        Log << "    " << operationTypeNames_[operation_]
-            << "(" << regionName_ << ") of " << fieldName
-            <<  " = " << result << endl;
+            Log << "    " << operationTypeNames_[operation_]
+                << "(" << regionName_ << ") of " << fieldName
+                <<  " = " << result << endl;
+        }

        return true;
    }
@ -426,8 +406,7 @@ template<class Type>
 Foam::tmp<Foam::Field<Type>>
 Foam::functionObjects::fieldValues::surfaceFieldValue::filterField
 (
-    const GeometricField<Type, fvPatchField, volMesh>& field,
-    const bool applyOrientation
+    const GeometricField<Type, fvPatchField, volMesh>& field
 ) const
 {
    tmp<Field<Type>> tvalues(new Field<Type>(faceId_.size()));
@ -435,8 +414,9 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::filterField

    forAll(values, i)
    {
-        label facei = faceId_[i];
-        label patchi = facePatchId_[i];
+        const label facei = faceId_[i];
+        const label patchi = facePatchId_[i];
+
        if (patchi >= 0)
        {
            values[i] = field.boundaryField()[patchi][facei];
@ -452,14 +432,6 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::filterField
        }
    }

-    if (applyOrientation)
-    {
-        forAll(values, i)
-        {
-            values[i] *= faceSign_[i];
-        }
-    }
-
    return tvalues;
 }

@ -468,8 +440,7 @@ template<class Type>
 Foam::tmp<Foam::Field<Type>>
 Foam::functionObjects::fieldValues::surfaceFieldValue::filterField
 (
-    const GeometricField<Type, fvsPatchField, surfaceMesh>& field,
-    const bool applyOrientation
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& field
 ) const
 {
    tmp<Field<Type>> tvalues(new Field<Type>(faceId_.size()));
@ -477,8 +448,9 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::filterField

    forAll(values, i)
    {
-        label facei = faceId_[i];
-        label patchi = facePatchId_[i];
+        const label facei = faceId_[i];
+        const label patchi = facePatchId_[i];
+
        if (patchi >= 0)
        {
            values[i] = field.boundaryField()[patchi][facei];
@ -489,14 +461,6 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::filterField
        }
    }

-    if (applyOrientation)
-    {
-        forAll(values, i)
-        {
-            values[i] *= faceSign_[i];
-        }
-    }
-
    return tvalues;
 }