Merge branch 'feature-field-statistics' into 'develop'

Draft: ENH: fieldStatistics: New function object for statistics of input fields

See merge request Development/openfoam!732

src/functionObjects/field/Make/files

@@ -154,4 +154,6 @@ resolutionIndex/resolutionIndexModels/CelikEtaIndex/CelikEtaIndex.C
 age/age.C
 comfort/comfort.C
 
+fieldStatistics/fieldStatistics.cxx
+
 LIB = $(FOAM_LIBBIN)/libfieldFunctionObjects

src/functionObjects/field/fieldStatistics/fieldStatistics.H

@@ -0,0 +1,356 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2025 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::functionObjects::fieldStatistics
+
+Group
+    grpFieldFunctionObjects
+
+Description
+    Calculates various statistics of the specified fields.
+
+    Operands:
+    \table
+      Operand      | Type                          | Location
+      input        | vol\<Type\>Field(s)           | Registry
+      output file  | dat | \<case\>/postProcessing/\<FO\>/\<time\>/\<file\>
+      output field | -                             | -
+    \endtable
+
+    where \c \<Type\>=Scalar/Vector/SphericalTensor/SymmTensor/Tensor.
+
+Usage
+    Minimal example by using \c system/controlDict.functions:
+    \verbatim
+    fieldStatistics1
+    {
+        // Mandatory entries
+        type             fieldStatistics;
+        libs             (fieldFunctionObjects);
+        fields           (<wordList>);
+        statistics       (<wordList>);
+
+        // Optional entries
+        mode             <word>;
+        mean             <word>;
+        extrema          <bool>;
+        internal         <bool>;
+
+        // Inherited entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                        | Type | Reqd  | Deflt
+      type         | Type name: fieldStatistics         | word | yes   | -
+      libs         | Library name: fieldFunctionObjects | word | yes   | -
+      fields       | List of operand fields             | wordList | yes | -
+      statistics   | List of operand statistics         | wordList | yes | -
+      mode   | Output format of the statistical results | word | no | magnitude
+      mean     | Type of the mean operation             | word | no | arithmetic
+      internal | Flag to use internal fields only in computing statistics <!--
+               -->                                      | bool | no | false
+      extrema | Flag to enable extrema data calculations | bool | no | false
+    \endtable
+
+    Available statistics of the operand field through the \c statistics entry:
+    \verbatim
+      min        | Minimum value
+      max        | Maximum value
+      mean       | Arithmetic mean value
+      variance   | Sample variance value (unbiased)
+    \endverbatim
+
+    Options for the \c mode entry:
+    \verbatim
+      magnitude  | Output statistics magnitude-wise
+      component  | Output statistics separately for each component
+    \endverbatim
+
+    Options for the \c meanType entry:
+    \verbatim
+      arithmetic  | Arithmetic mean (average)
+      volumetric  | Volume-weighted arithmetic mean
+    \endverbatim
+
+    The inherited entries are elaborated in:
+      - \link functionObject.H \endlink
+      - \link writeFile.H \endlink
+
+SourceFiles
+    fieldStatistics.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjects_fieldStatistics_H
+#define functionObjects_fieldStatistics_H
+
+#include "fvMeshFunctionObject.H"
+#include "writeFile.H"
+#include "volFieldSelection.H"
+#include <functional>
+#include <variant>
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class fieldStatistics Declaration
+\*---------------------------------------------------------------------------*/
+
+class fieldStatistics
+:
+    public fvMeshFunctionObject,
+    public writeFile
+{
+    // Private Enumerations
+
+        //- Options for the output format of the statistical results
+        enum modeType : char
+        {
+            mdMag = 0,  //!< "Output statistics magnitude-wise"
+            mdCmpt      //!< "Output statistics separately for each component"
+        };
+
+        //- Names for modeType
+        static const Enum<modeType> modeTypeNames_;
+
+
+        //- Options for the mean type of the specified fields
+        enum meanType : char
+        {
+            ARITHMETIC = 0,  //!< "Arithmetic mean (average)"
+            VOLUMETRIC       //!< "Volume-weighted arithmetic mean"
+        };
+
+        //- Names for meanType
+        static const Enum<meanType> meanTypeNames_;
+
+
+    // Private Classes
+
+        //- Type-safe union for input field types
+        using variantInput = std::variant
+        <
+            scalarField,
+            vectorField,
+            sphericalTensorField,
+            symmTensorField,
+            tensorField
+        >;
+
+        //- Type-safe union for output data types
+        using variantOutput = std::variant
+        <
+            scalar,
+            vector,
+            sphericalTensor,
+            symmTensor,
+            tensor
+        >;
+
+        //- Class to encapsulate information about specified statistic
+        struct statistic
+        {
+            //- Name of the statistic
+            word name_;
+
+            //- Return the value of the specified statistic
+            std::function<variantOutput(variantInput)> calc;
+        };
+
+        //- Class to encapsulate the data about minima and maxima
+        struct extremaData
+        {
+            //- Value of the extremum
+            variantOutput value_;
+
+            //- Processor index of the extremum
+            label procID_;
+
+            //- Cell index of the extremum
+            label cellID_;
+
+            //- Position vector of the extremum
+            vector position_;
+        };
+
+
+    // Private Data
+
+        //- Flag to use internal fields only in computing statistics
+        bool internal_;
+
+        //- Flag to enable extrema data calculations
+        bool extrema_;
+
+        //- Output-format mode - only applicable for tensor ranks > 0
+        modeType mode_;
+
+        //- Type of the mean of the specified fields
+        meanType mean_;
+
+        //- Operand fields on which statistics are computed
+        volFieldSelection fieldSet_;
+
+        //- List of file pointers; one file per field
+        HashPtrTable<OFstream> filePtrs_;
+
+        //- List of file pointers for extrema data; one file per field
+        HashPtrTable<OFstream> extremaFilePtrs_;
+
+        //- Hash table containing all specified statistics
+        HashTable<statistic> statistics_;
+
+        //- Hash table containing all statistical results per field
+        HashTable<HashTable<variantOutput>> results_;
+
+        //- Hash table containing the results of the extrema per field
+        HashTable<Pair<extremaData>> extremaResults_;
+
+
+    // Private Member Functions
+
+        //- Return the statistic container
+        statistic createStatistic(const word& statName, const modeType mode);
+
+        //- Compute the specified statistics of a given field
+        template<class T>
+        bool calcStat(const word& fieldName);
+
+
+        // Central tendency statistics
+
+        //- Return the arithmetic mean of the given input field
+        template<class T>
+        T calcMean(const Field<T>& field);
+
+
+        // Dispersion statistics
+
+        //- Return the minimum value of the given input field
+        //  Store the processor index, cell index and location of the minimum
+        template<class T>
+        T calcMin(const Field<T>& field);
+
+        //- Return the maximum value of the given input field
+        //  Store the processor index, cell index and location of the maximum
+        template<class T>
+        T calcMax(const Field<T>& field);
+
+        //- Return the sample variance of the given input field
+        template<class T>
+        T calcVariance(const Field<T>& field);
+
+        //-
+        template<class GeoField>
+        tmp<Field<typename GeoField::value_type>> flatten
+        (
+            const GeoField& field
+        );
+
+        //- Return the extrema data of the specified field
+        template<class GeoField>
+        Pair<extremaData> calcExtremaData(const GeoField& field);
+
+        //- Output the file header information
+        void writeFileHeader(Ostream& os, const word& fieldName);
+
+        //- Write the statistical data to the specified file
+        void writeStatData();
+
+        //- Write the statistical data to the standard stream
+        void logStatData();
+
+        //- Output the extrema-data file header information
+        void writeExtremaFileHeader(Ostream& os, const word& fieldName);
+
+        //- Write extrema data to the specified file
+        void writeExtremaData();
+
+        //- Write extrema data to the standard stream
+        void logExtremaData();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("fieldStatistics");
+
+
+    // Generated Methods
+
+        //- No copy construct
+        fieldStatistics(const fieldStatistics&) = delete;
+
+        //- No copy assignment
+        void operator=(const fieldStatistics&) = delete;
+
+
+    // Constructors
+
+        //- Construct from name, Time and dictionary
+        fieldStatistics
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~fieldStatistics() = default;
+
+
+    // Member Functions
+
+        //- Read function object settings
+        virtual bool read(const dictionary&);
+
+        //- Execute function object
+        virtual bool execute();
+
+        //- Write function object results
+        virtual bool write();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/functionObjects/field/fieldStatistics/fieldStatistics.cxx

@@ -0,0 +1,505 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2025 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 "fieldStatistics.H"
+#include "fieldTypes.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(fieldStatistics, 0);
+    addToRunTimeSelectionTable(functionObject, fieldStatistics, dictionary);
+}
+}
+
+const Foam::Enum
+<
+    Foam::functionObjects::fieldStatistics::modeType
+>
+Foam::functionObjects::fieldStatistics::modeTypeNames_
+({
+    { modeType::mdMag,  "magnitude" },
+    { modeType::mdCmpt, "component" },
+});
+
+const Foam::Enum
+<
+    Foam::functionObjects::fieldStatistics::meanType
+>
+Foam::functionObjects::fieldStatistics::meanTypeNames_
+({
+    { meanType::ARITHMETIC, "arithmetic" },
+    { meanType::VOLUMETRIC, "volumetric" },
+});
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Implementation
+#include "fieldStatisticsImpl.cxx"
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::functionObjects::fieldStatistics::statistic
+Foam::functionObjects::fieldStatistics::createStatistic
+(
+    const word& statName,
+    const modeType mode
+)
+{
+    statistic m;
+    m.name_ = statName;
+
+    m.calc = [this, statName, mode](variantInput input) -> variantOutput
+    {
+        return std::visit
+        (
+            [this, statName, mode](auto&& arg) -> variantOutput
+            {
+                using T = std::decay_t<decltype(arg)>;
+                using BaseType = typename T::value_type;
+
+                if (statName == "min")
+                {
+                    // if (mode == mdMag) return mag(calcMin<BaseType>(arg));
+                    if (mode == mdMag) return calcMin<scalar>(mag(arg));
+                    else return calcMin<BaseType>(arg);
+                }
+                else if (statName == "max")
+                {
+                    if (mode == mdMag) return calcMax<scalar>(mag(arg));
+                    else return calcMax<BaseType>(arg);
+                }
+                else if (statName == "mean")
+                {
+                    if (mode == mdMag) return calcMean<scalar>(mag(arg));
+                    else return calcMean<BaseType>(arg);
+                }
+                else if (statName == "variance")
+                {
+                    if (mode == mdMag) return calcVariance<scalar>(mag(arg));
+                    else return calcVariance<BaseType>(arg);
+                }
+                else return scalar{};  // Default case (for compiler)
+            },
+            input
+        );
+    };
+
+    return m;
+}
+
+
+void Foam::functionObjects::fieldStatistics::writeFileHeader
+(
+    Ostream& os,
+    const word& fieldName
+)
+{
+    writeHeader(os, word("Field Statistics: " + fieldName));
+    writeCommented(os, "Time");
+
+    // Number of input statistics (i.e., statistics_) should be the same with
+    // that of output statistics (i.e., results_). However, for consistency,
+    // the output file columns are based on output statistics.
+    const auto& result = results_(fieldName);
+
+    for (const auto& iter : result.csorted())
+    {
+        const word& name = iter.key();
+        writeTabbed(os, name);
+    }
+    os  << endl;
+}
+
+
+void Foam::functionObjects::fieldStatistics::writeStatData()
+{
+    for (const word& fieldName : fieldSet_.selectionNames())
+    {
+        const auto& results = results_(fieldName);
+
+        if (!results.size()) break;
+
+        OFstream& file = *filePtrs_(fieldName);
+
+        writeCurrentTime(file);
+
+        for (const auto& iter : results.csorted())
+        {
+            const variantOutput& value = iter.val();
+
+            std::visit
+            (
+                [&file](const auto& v)
+                {
+                    if constexpr
+                    (
+                        is_vectorspace_v<std::decay_t<decltype(v)>>
+                    )
+                    {
+                        for (const auto& val : v) file<< token::TAB << val;
+                    }
+                    else
+                    {
+                        file<< token::TAB << v;
+                    }
+                },
+                value
+            );
+        }
+        file<< nl;
+    }
+}
+
+
+void Foam::functionObjects::fieldStatistics::logStatData()
+{
+    for (const word& fieldName : fieldSet_.selectionNames())
+    {
+        const auto& results = results_(fieldName);
+
+        if (!results.size()) break;
+
+        const word outputName =
+            (mode_ == mdMag) ? word("mag(" + fieldName + ")") : fieldName;
+
+        Info<< nl << "    Field " << outputName << nl;
+
+        for (const auto& iter : results.csorted())
+        {
+            const word& name = iter.key();
+            const variantOutput& value = iter.val();
+
+            Info<< "    " << name;
+            std::visit
+            (
+                [](const auto& v)
+                {
+                    if constexpr
+                    (
+                        is_vectorspace_v<std::decay_t<decltype(v)>>
+                    )
+                    {
+                        for (const auto& val : v) Info<< " " << val;
+                    }
+                    else
+                    {
+                        Info<< " " << v;
+                    }
+                },
+                value
+            );
+            Info<< nl;
+        }
+    }
+    Info<< endl;
+}
+
+
+void Foam::functionObjects::fieldStatistics::writeExtremaFileHeader
+(
+    Ostream& os,
+    const word& fieldName
+)
+{
+    writeHeader(os, word("Field Extrema Data: " + fieldName));
+    writeCommented(os, "Time");
+    writeTabbed(os, "min");
+    writeTabbed(os, "min_procID");
+    writeTabbed(os, "min_cellID");
+    writeTabbed(os, "min_position");
+    writeTabbed(os, "max");
+    writeTabbed(os, "max_procID");
+    writeTabbed(os, "max_cellID");
+    writeTabbed(os, "max_position");
+    os  << endl;
+}
+
+
+void Foam::functionObjects::fieldStatistics::writeExtremaData()
+{
+    for (const word& fieldName : fieldSet_.selectionNames())
+    {
+        const auto& min = extremaResults_(fieldName).first();
+        const auto& max = extremaResults_(fieldName).second();
+
+        OFstream& file = *extremaFilePtrs_(fieldName);
+
+        writeCurrentTime(file);
+
+        file<< token::TAB;
+
+        std::visit([&file](const auto& v){ file<< v; }, min.value_);
+
+        file<< token::TAB
+            << min.procID_ << token::TAB
+            << min.cellID_ << token::TAB
+            << min.position_ << token::TAB;
+
+        std::visit([&file](const auto& v){ file<< v; }, max.value_);
+
+        file<< token::TAB
+            << max.procID_ << token::TAB
+            << max.cellID_ << token::TAB
+            << max.position_ << nl;
+    }
+}
+
+
+void Foam::functionObjects::fieldStatistics::logExtremaData()
+{
+    for (const word& fieldName : fieldSet_.selectionNames())
+    {
+        const auto& min = extremaResults_(fieldName).first();
+        const auto& max = extremaResults_(fieldName).second();
+
+        const word outputName =
+            (mode_ == mdMag) ? word("mag(" + fieldName + ")") : fieldName;
+
+        std::visit
+        (
+            [outputName](const auto& v)
+            {
+                Info<< "    min(" << outputName << ") = " << v;
+            },
+            min.value_
+        );
+
+        Info<< " in cell " << min.cellID_
+            << " at location " << min.position_
+            << " on processor " << min.procID_;
+
+        std::visit
+        (
+            [outputName](const auto& v)
+            {
+                Info<< nl << "    max(" << outputName << ") = " << v;
+            },
+            max.value_
+        );
+
+        Info<< " in cell " << max.cellID_
+            << " at location " << max.position_
+            << " on processor " << max.procID_ << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::fieldStatistics::fieldStatistics
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fvMeshFunctionObject(name, runTime, dict),
+    writeFile(mesh_, name, typeName, dict),
+    fieldSet_(mesh_)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::fieldStatistics::read(const dictionary& dict)
+{
+    if (!(fvMeshFunctionObject::read(dict) && writeFile::read(dict)))
+    {
+        return false;
+    }
+
+    internal_ = dict.getOrDefault("internal", false);
+    extrema_ = dict.getOrDefault("extrema", false);
+
+    mode_ = modeTypeNames_.getOrDefault("mode", dict, modeType::mdMag);
+    mean_ = meanTypeNames_.getOrDefault("mean", dict, meanType::ARITHMETIC);
+
+    // Reset and reprepare the input field names
+    fieldSet_.clear();
+    fieldSet_.read(dict);
+
+    // Reset and reprepare the input statistics
+    statistics_.clear();
+    const wordHashSet stats(dict.get<wordHashSet>("statistics"));
+    for (const word& m : stats)
+    {
+        statistics_.insert(m, createStatistic(m, mode_));
+    }
+
+    // Reset the output-statistics container
+    results_.clear();
+
+    // Reset the extrema-data container
+    extremaResults_.clear();
+
+    return true;
+}
+
+
+bool Foam::functionObjects::fieldStatistics::execute()
+{
+    fieldSet_.updateSelection();
+
+    // Calculate the specified statistics for the specified fields
+    for (const word& fieldName : fieldSet_.selectionNames())
+    {
+        const bool ok =
+        (
+            calcStat<scalar>(fieldName)
+         || calcStat<vector>(fieldName)
+         || calcStat<sphericalTensor>(fieldName)
+         || calcStat<symmTensor>(fieldName)
+         || calcStat<tensor>(fieldName)
+        );
+
+        if (!ok)
+        {
+            WarningInFunction
+                << "Unable to find field " << fieldName << endl;
+        }
+    }
+
+    // Store the statistical results into the state containers
+    for (const word& fieldName : fieldSet_.selectionNames())
+    {
+        const auto& results = results_(fieldName);
+
+        for (const auto& iter : results.csorted())
+        {
+            const word& name = iter.key();
+            const variantOutput& value = iter.val();
+
+            const word variableName(fieldName + "_" + name);
+
+            std::visit
+            (
+                [this, variableName](const auto& v)
+                {
+                    this->setResult(variableName, v);
+                },
+                value
+            );
+        }
+
+        if (extrema_)
+        {
+            const auto& min = extremaResults_(fieldName).first();
+            std::visit
+            (
+                [this, fieldName](const auto& v)
+                {
+                    this->setResult(word(fieldName + "_min"), v);
+                },
+                min.value_
+            );
+            this->setResult(word(fieldName + "_min_procID"), min.procID_);
+            this->setResult(word(fieldName + "_min_cellID"), min.cellID_);
+            this->setResult(word(fieldName + "_min_position"), min.position_);
+
+            const auto& max = extremaResults_(fieldName).second();
+            std::visit
+            (
+                [this, fieldName](const auto& v)
+                {
+                    this->setResult(word(fieldName + "_max"), v);
+                },
+                max.value_
+            );
+            this->setResult(word(fieldName + "_max_procID"), max.procID_);
+            this->setResult(word(fieldName + "_max_cellID"), max.cellID_);
+            this->setResult(word(fieldName + "_max_position"), max.position_);
+        }
+    }
+
+    return true;
+}
+
+
+bool Foam::functionObjects::fieldStatistics::write()
+{
+    Log << type() << " " << name() << " write:" << endl;
+
+    // Create an output file per field
+    if (writeToFile() && !writtenHeader_)
+    {
+        for (const word& fieldName : fieldSet_.selectionNames())
+        {
+            filePtrs_.set(fieldName, newFileAtStartTime(fieldName));
+
+            OFstream& file = *filePtrs_(fieldName);
+
+            writeFileHeader(file, fieldName);
+        }
+
+        if (extrema_)
+        {
+            for (const word& fieldName : fieldSet_.selectionNames())
+            {
+                extremaFilePtrs_.set
+                (
+                    fieldName,
+                    newFileAtStartTime(word("fieldMinMax_" + fieldName))
+                );
+
+                OFstream& file = *extremaFilePtrs_(fieldName);
+                writeExtremaFileHeader(file, fieldName);
+            }
+        }
+
+        writtenHeader_ = true;
+    }
+
+    // Write the statistical results to the output file if requested
+    if (writeToFile())
+    {
+        if (extrema_) writeExtremaData();
+
+        writeStatData();
+    }
+
+    // Print the statistical results to the standard stream if requested
+    if (log)
+    {
+        if (extrema_) logExtremaData();
+
+        logStatData();
+    }
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/functionObjects/field/fieldStatistics/fieldStatisticsImpl.cxx

@@ -0,0 +1,271 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2025 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 "fieldStatistics.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class GeoField>
+Foam::tmp<Foam::Field<typename GeoField::value_type>>
+Foam::functionObjects::fieldStatistics::flatten(const GeoField& fld)
+{
+    typedef typename GeoField::value_type value_type;
+    typedef Field<value_type> FieldType;
+
+
+    label n = fld.size();
+
+    if (!internal_)
+    {
+        for (const auto& pfld : fld.boundaryField())
+        {
+            if (!pfld.coupled())
+            {
+                n += pfld.size();
+            }
+        }
+    }
+
+    auto tflatFld = tmp<FieldType>::New(n);
+    auto& flatFld = tflatFld.ref();
+
+    // Insert internal values
+    flatFld.slice(0, fld.size()) = fld.primitiveField();
+
+    if (!internal_)
+    {
+        // Insert boundary values
+        n = fld.size();
+        for (const auto& pfld : fld.boundaryField())
+        {
+            if (!pfld.coupled())
+            {
+                flatFld.slice(n, pfld.size()) = pfld;
+                n += pfld.size();
+            }
+        }
+    }
+
+    return tflatFld;
+}
+
+
+template<class Type>
+bool Foam::functionObjects::fieldStatistics::calcStat(const word& fieldName)
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
+
+    const auto* fieldp = obr_.cfindObject<VolFieldType>(fieldName);
+    if (!fieldp)
+    {
+        return false;
+    }
+    const auto& field = *fieldp;
+
+    tmp<Field<Type>> tfld = flatten(field);
+    const auto& fld = tfld.cref();
+
+    HashTable<variantOutput> result;
+    for (const auto& iter : statistics_.csorted())
+    {
+        const statistic& stat = iter.val();
+
+        // Assign a new entry, overwriting existing entries
+        result.set(stat.name_, stat.calc(fld));
+    }
+
+    results_.set(fieldName, result);
+
+    if (extrema_)
+    {
+        extremaResults_.set
+        (
+            fieldName,
+            (mode_ == mdMag)
+          ? calcExtremaData(mag(field)())
+          : calcExtremaData(field)
+        );
+    }
+
+    return true;
+}
+
+
+template<class T>
+T Foam::functionObjects::fieldStatistics::calcMean(const Field<T>& field)
+{
+    if (internal_ && (mean_ == VOLUMETRIC))
+    {
+        return gWeightedAverage(mesh_.V(), field);
+    }
+
+    return gAverage(field);
+}
+
+
+template<class T>
+T Foam::functionObjects::fieldStatistics::calcMin(const Field<T>& field)
+{
+    return gMin(field);
+}
+
+
+template<class T>
+T Foam::functionObjects::fieldStatistics::calcMax(const Field<T>& field)
+{
+    return gMax(field);
+}
+
+
+template<class T>
+T Foam::functionObjects::fieldStatistics::calcVariance
+(
+    const Field<T>& field
+)
+{
+    const T average(calcMean(field));
+
+    T var = Zero;
+    for (const auto& elem : field)
+    {
+        var += (elem - average);
+    }
+
+    label n = field.size();
+
+    sumReduce(var, n);
+
+    if (n <= 1)
+    {
+        return T{};
+    }
+
+    return 1.0/(n - 1.0)*var;
+}
+
+
+template<class GeoField>
+Foam::Pair<Foam::functionObjects::fieldStatistics::extremaData>
+Foam::functionObjects::fieldStatistics::calcExtremaData
+(
+    const GeoField& field
+)
+{
+    typedef typename GeoField::value_type value_type;
+
+    const label proci = Pstream::myProcNo();
+
+    // Find the extrema data of the specified internal field
+
+    List<value_type> minVs(Pstream::nProcs(), pTraits<value_type>::max);
+    List<label> minCells(Pstream::nProcs(), Zero);
+    List<vector> minCs(Pstream::nProcs(), Zero);
+
+    List<value_type> maxVs(Pstream::nProcs(), pTraits<value_type>::min);
+    List<label> maxCells(Pstream::nProcs(), Zero);
+    List<vector> maxCs(Pstream::nProcs(), Zero);
+
+    labelPair minMaxIds = findMinMax(field);
+
+    if (label celli = minMaxIds.first(); celli >= 0)
+    {
+        minVs[proci] = field[celli];
+        minCells[proci] = celli;
+        minCs[proci] = mesh_.C()[celli];
+    }
+
+    if (label celli = minMaxIds.second(); celli >= 0)
+    {
+        maxVs[proci] = field[celli];
+        maxCells[proci] = celli;
+        maxCs[proci] = mesh_.C()[celli];
+    }
+
+    if (!internal_)
+    {
+        // Find the extrema data of the specified boundary fields
+        const auto& fieldBoundary = field.boundaryField();
+        const auto& CfBoundary = mesh_.C().boundaryField();
+
+        forAll(fieldBoundary, patchi)
+        {
+            const Field<value_type>& fp = fieldBoundary[patchi];
+            if (fp.size())
+            {
+                const vectorField& Cfp = CfBoundary[patchi];
+
+                const labelList& faceCells =
+                    fieldBoundary[patchi].patch().faceCells();
+
+                minMaxIds = findMinMax(fp);
+
+                if (label celli = minMaxIds.first(); minVs[proci] > fp[celli])
+                {
+                    minVs[proci] = fp[celli];
+                    minCells[proci] = faceCells[celli];
+                    minCs[proci] = Cfp[celli];
+                }
+
+                if (label celli = minMaxIds.second(); maxVs[proci] < fp[celli])
+                {
+                    maxVs[proci] = fp[celli];
+                    maxCells[proci] = faceCells[celli];
+                    maxCs[proci] = Cfp[celli];
+                }
+            }
+        }
+    }
+
+    // Collect info from all processors and output
+    Pstream::allGatherList(minVs);
+    Pstream::allGatherList(minCells);
+    Pstream::allGatherList(minCs);
+
+    Pstream::allGatherList(maxVs);
+    Pstream::allGatherList(maxCells);
+    Pstream::allGatherList(maxCs);
+
+    extremaData min;
+    const label minId = findMin(minVs);
+    min.value_ = minVs[minId];
+    min.procID_ = minId;
+    min.cellID_ = minCells[minId];
+    min.position_ = minCs[minId];
+
+    extremaData max;
+    const label maxId = findMax(maxVs);
+    max.value_ = maxVs[maxId];
+    max.procID_ = maxId;
+    max.cellID_ = maxCells[maxId];
+    max.position_ = maxCs[maxId];
+
+    return Pair<extremaData>(min, max);
+}
+
+
+// ************************************************************************* //