diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
index 0103e844fa..16020e060b 100644
--- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
+++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
@@ -564,13 +564,12 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::writeFileHeader
 }
 
 
-template<>
-Foam::scalar Foam::functionObjects::fieldValues::surfaceFieldValue::
-processValues
+bool Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
 (
     const Field<scalar>& values,
     const vectorField& Sf,
-    const scalarField& weightField
+    const scalarField& weightField,
+    scalar& result
 ) const
 {
     switch (operation_)
@@ -578,31 +577,59 @@ processValues
         case operationType::sumDirection:
         {
             vector n(dict_.lookup("direction"));
-            return sum(pos0(values*(Sf & n))*mag(values));
+            result = sum(pos0(values*(Sf & n))*mag(values));
+            return true;
         }
         case operationType::sumDirectionBalance:
         {
             vector n(dict_.lookup("direction"));
             const scalarField nv(values*(Sf & n));
-
-            return sum(pos0(nv)*mag(values) - neg(nv)*mag(values));
+            result = sum(pos0(nv)*mag(values) - neg(nv)*mag(values));
+            return true;
         }
         default:
         {
-            // Fall through to other operations
-            return processSameTypeValues(values, Sf, weightField);
+            // Fall through to same-type operations
+            return processValuesTypeType(values, Sf, weightField, result);
         }
     }
 }
 
 
-template<>
-Foam::vector Foam::functionObjects::fieldValues::surfaceFieldValue::
-processValues
+bool Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
 (
     const Field<vector>& values,
     const vectorField& Sf,
-    const scalarField& weightField
+    const scalarField& weightField,
+    scalar& result
+) const
+{
+    switch (operation_)
+    {
+        case operationType::areaNormalAverage:
+        {
+            result = sum(values & Sf)/sum(mag(Sf));
+            return true;
+        }
+        case operationType::areaNormalIntegrate:
+        {
+            result = sum(values & Sf);
+            return true;
+        }
+        default:
+        {
+            return false;
+        }
+    }
+}
+
+
+bool Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
+(
+    const Field<vector>& values,
+    const vectorField& Sf,
+    const scalarField& weightField,
+    vector& result
 ) const
 {
     switch (operation_)
@@ -612,31 +639,21 @@ processValues
             vector n(dict_.lookup("direction"));
             n /= mag(n) + rootVSmall;
             const scalarField nv(n & values);
-
-            return sum(pos0(nv)*n*(nv));
+            result = sum(pos0(nv)*n*(nv));
+            return true;
         }
         case operationType::sumDirectionBalance:
         {
             vector n(dict_.lookup("direction"));
             n /= mag(n) + rootVSmall;
             const scalarField nv(n & values);
-
-            return sum(pos0(nv)*n*(nv));
-        }
-        case operationType::areaNormalAverage:
-        {
-            scalar result = sum(values & Sf)/sum(mag(Sf));
-            return vector(result, 0.0, 0.0);
-        }
-        case operationType::areaNormalIntegrate:
-        {
-            scalar result = sum(values & Sf);
-            return vector(result, 0.0, 0.0);
+            result = sum(pos0(nv)*n*(nv));
+            return true;
         }
         default:
         {
-            // Fall through to other operations
-            return processSameTypeValues(values, Sf, weightField);
+            // Fall through to same-type operations
+            return processValuesTypeType(values, Sf, weightField, result);
         }
     }
 }
@@ -792,8 +809,7 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::write()
         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<sphericalTensor>(fieldName, weightField, orient);
         ok = ok || writeValues<symmTensor>(fieldName, weightField, orient);
         ok = ok || writeValues<tensor>(fieldName, weightField, orient);
 
diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H
index 7578c50960..8f7f9c64c7 100644
--- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H
+++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -332,24 +332,67 @@ protected:
             const bool applyOrientation = false
         ) const;
 
-        //- Apply the 'operation' to the values. Operation has to
-        //  preserve Type.
-        template<class Type>
-        Type processSameTypeValues
+        //- Apply the operation to the values, and return true if successful.
+        //  Does nothing unless overloaded below.
+        template<class Type, class ResultType>
+        bool processValues
         (
             const Field<Type>& values,
             const vectorField& Sf,
-            const scalarField& weightField
+            const scalarField& weightField,
+            ResultType& result
         ) const;
 
-        //- Apply the 'operation' to the values. Wrapper around
-        //  processSameTypeValues. See also template specialisation below.
+        //- Apply Type -> Type operation to the values. Calls
+        //  processValuesTypeType.
         template<class Type>
-        Type processValues
+        bool processValues
         (
             const Field<Type>& values,
             const vectorField& Sf,
-            const scalarField& weightField
+            const scalarField& weightField,
+            Type& 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 vectorField& Sf,
+            const scalarField& weightField,
+            scalar& result
+        ) const;
+
+        //- Apply vector -> vector operation to the values.
+        bool processValues
+        (
+            const Field<vector>& values,
+            const vectorField& Sf,
+            const scalarField& weightField,
+            scalar& result
+        ) const;
+
+        //- Apply vector -> vector operation to the values. Tries to apply
+        //  vector -> vector specific operations, otherwise calls
+        //  processValuesTypeType.
+        bool processValues
+        (
+            const Field<vector>& values,
+            const vectorField& Sf,
+            const scalarField& weightField,
+            vector& result
+        ) const;
+
+        //- Apply a Type -> Type operation to the values
+        template<class Type>
+        bool processValuesTypeType
+        (
+            const Field<Type>& values,
+            const vectorField& Sf,
+            const scalarField& weightField,
+            Type& result
         ) const;
 
         //- Output file header information
@@ -411,6 +454,16 @@ public:
             const bool orient
         );
 
+        //- Templated helper function to output field values
+        template<class Type, class ResultType>
+        bool writeValues
+        (
+            const word& fieldName,
+            const scalarField& weightField,
+            const Field<Type>& values,
+            const vectorField& Sf
+        );
+
         //- Filter a surface field according to faceIds
         template<class Type>
         tmp<Field<Type>> filterField
@@ -435,26 +488,6 @@ public:
 };
 
 
-//- Specialisation for scalar
-template<>
-scalar surfaceFieldValue::processValues
-(
-    const Field<scalar>& values,
-    const vectorField& Sf,
-    const scalarField& weightField
-) const;
-
-
-//- Specialisation for vector
-template<>
-vector surfaceFieldValue::processValues
-(
-    const Field<vector>& values,
-    const vectorField& Sf,
-    const scalarField& weightField
-) const;
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace fieldValues
diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C
index 32f716baee..365428da15 100644
--- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C
+++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -132,22 +132,48 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::getFieldValues
 }
 
 
-template<class Type>
-Type Foam::functionObjects::fieldValues::surfaceFieldValue::
-processSameTypeValues
+template<class Type, class ResultType>
+bool Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
 (
     const Field<Type>& values,
     const vectorField& Sf,
-    const scalarField& weightField
+    const scalarField& weightField,
+    ResultType& result
+) const
+{
+    return false;
+}
+
+
+template<class Type>
+bool Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
+(
+    const Field<Type>& values,
+    const vectorField& Sf,
+    const scalarField& weightField,
+    Type& result
+) const
+{
+    return processValuesTypeType(values, Sf, weightField, result);
+}
+
+
+template<class Type>
+bool Foam::functionObjects::fieldValues::surfaceFieldValue::
+processValuesTypeType
+(
+    const Field<Type>& values,
+    const vectorField& Sf,
+    const scalarField& weightField,
+    Type& result
 ) const
 {
-    Type result = Zero;
     switch (operation_)
     {
         case operationType::sum:
         {
             result = sum(values);
-            break;
+            return true;
         }
         case operationType::weightedSum:
         {
@@ -159,39 +185,17 @@ processSameTypeValues
             {
                 result = sum(values);
             }
-            break;
+            return true;
         }
         case operationType::sumMag:
         {
             result = sum(cmptMag(values));
-            break;
-        }
-        case operationType::sumDirection:
-        {
-            FatalErrorInFunction
-                << "Operation " << operationTypeNames_[operation_]
-                << " not available for values of type "
-                << pTraits<Type>::typeName
-                << exit(FatalError);
-
-            result = Zero;
-            break;
-        }
-        case operationType::sumDirectionBalance:
-        {
-            FatalErrorInFunction
-                << "Operation " << operationTypeNames_[operation_]
-                << " not available for values of type "
-                << pTraits<Type>::typeName
-                << exit(FatalError);
-
-            result = Zero;
-            break;
+            return true;
         }
         case operationType::average:
         {
             result = sum(values)/values.size();
-            break;
+            return true;
         }
         case operationType::weightedAverage:
         {
@@ -205,14 +209,14 @@ processSameTypeValues
             {
                 result = sum(values)/values.size();
             }
-            break;
+            return true;
         }
         case operationType::areaAverage:
         {
             const scalarField magSf(mag(Sf));
 
             result = sum(magSf*values)/sum(magSf);
-            break;
+            return true;
         }
         case operationType::weightedAreaAverage:
         {
@@ -228,14 +232,14 @@ processSameTypeValues
             {
                 result = sum(magSf*values)/sum(magSf);
             }
-            break;
+            return true;
         }
         case operationType::areaIntegrate:
         {
             const scalarField magSf(mag(Sf));
 
             result = sum(magSf*values);
-            break;
+            return true;
         }
         case operationType::weightedAreaIntegrate:
         {
@@ -249,17 +253,17 @@ processSameTypeValues
             {
                 result = sum(magSf*values);
             }
-            break;
+            return true;
         }
         case operationType::min:
         {
             result = min(values);
-            break;
+            return true;
         }
         case operationType::max:
         {
             result = max(values);
-            break;
+            return true;
         }
         case operationType::CoV:
         {
@@ -278,33 +282,20 @@ processSameTypeValues
                 res = sqrt(sum(magSf*sqr(vals - mean))/sum(magSf))/mean;
             }
 
-            break;
+            return true;
         }
-        case operationType::areaNormalAverage:
-        {}
-        case operationType::areaNormalIntegrate:
-        {}
         case operationType::none:
-        {}
+        {
+            return true;
+        }
+        default:
+        {
+            return false;
+        }
     }
-
-    return result;
 }
 
 
-template<class Type>
-Type Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
-(
-    const Field<Type>& values,
-    const vectorField& Sf,
-    const scalarField& weightField
-) const
-{
-    return processSameTypeValues(values, Sf, weightField);
-}
-
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
@@ -374,16 +365,26 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::writeValues
 
             if (Pstream::master())
             {
-                Type result = processValues(values, Sf, weightField);
-
-                // Add to result dictionary, over-writing any previous entry
-                resultDict_.add(fieldName, result, true);
-
-                file() << tab << result;
-
-                Log << "    " << operationTypeNames_[operation_]
-                    << "(" << regionName_ << ") of " << fieldName
-                    <<  " = " << result << endl;
+                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);
+                }
             }
         }
     }
@@ -392,6 +393,35 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::writeValues
 }
 
 
+template<class Type, class ResultType>
+bool Foam::functionObjects::fieldValues::surfaceFieldValue::writeValues
+(
+    const word& fieldName,
+    const scalarField& weightField,
+    const Field<Type>& values,
+    const vectorField& Sf
+)
+{
+    ResultType result;
+
+    if (processValues(values, Sf, weightField, result))
+    {
+        // Add to result dictionary, over-writing any previous entry
+        resultDict_.add(fieldName, result, true);
+
+        file() << tab << result;
+
+        Log << "    " << operationTypeNames_[operation_]
+            << "(" << regionName_ << ") of " << fieldName
+            <<  " = " << result << endl;
+
+        return true;
+    }
+
+    return false;
+}
+
+
 template<class Type>
 Foam::tmp<Foam::Field<Type>>
 Foam::functionObjects::fieldValues::surfaceFieldValue::filterField