diff --git a/src/OpenFOAM/db/functionObjects/writeFile/writeFile.H b/src/OpenFOAM/db/functionObjects/writeFile/writeFile.H
index 11df4d403c7ac303980f28300364efd43d5f62c6..cb89d303ad37c49733678e552ac87487e5a82594 100644
--- a/src/OpenFOAM/db/functionObjects/writeFile/writeFile.H
+++ b/src/OpenFOAM/db/functionObjects/writeFile/writeFile.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -237,6 +237,14 @@ public:
             const string& property,
             const Type& value
         ) const;
+
+        //- Write a given value to stream with the space delimiter
+        template<class Type>
+        void writeValue
+        (
+            Ostream& os,
+            const Type& val
+        ) const;
 };
 
 
diff --git a/src/OpenFOAM/db/functionObjects/writeFile/writeFileTemplates.C b/src/OpenFOAM/db/functionObjects/writeFile/writeFileTemplates.C
index 36b66a4296dce13f0ce30248149f8ee35f149c71..2283958151c6f0c362d3f9e144872b41c5a5f0a1 100644
--- a/src/OpenFOAM/db/functionObjects/writeFile/writeFileTemplates.C
+++ b/src/OpenFOAM/db/functionObjects/writeFile/writeFileTemplates.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
+    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -41,4 +42,18 @@ void Foam::functionObjects::writeFile::writeHeaderValue
 }
 
 
+template<class Type>
+void Foam::functionObjects::writeFile::writeValue
+(
+    Ostream& os,
+    const Type& val
+) const
+{
+    for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; ++cmpt)
+    {
+        os  << ' ' << component(val, cmpt);
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/functionObjects/field/Make/files b/src/functionObjects/field/Make/files
index 59fabd85cd3c15a171028e05a2ec88f986bedaec..b7c2550efbde555b67f5c74433539b100f13c553 100644
--- a/src/functionObjects/field/Make/files
+++ b/src/functionObjects/field/Make/files
@@ -1,5 +1,11 @@
 AMIWeights/AMIWeights.C
 
+binField/binField.C
+binField/binModels/binModel/binModel.C
+binField/binModels/binModel/binModelNew.C
+binField/binModels/uniformBin/uniformBin.C
+binField/binModels/singleDirectionUniformBin/singleDirectionUniformBin.C
+
 columnAverage/columnAverage.C
 
 continuityError/continuityError.C
diff --git a/src/functionObjects/field/binField/binField.C b/src/functionObjects/field/binField/binField.C
new file mode 100644
index 0000000000000000000000000000000000000000..b805afcd9ae1289567a457c0627bcbb52dc3efce
--- /dev/null
+++ b/src/functionObjects/field/binField/binField.C
@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 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 "binField.H"
+#include "binModel.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(binField, 0);
+    addToRunTimeSelectionTable(functionObject, binField, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::binField::binField
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict,
+    bool readFields
+)
+:
+    fvMeshFunctionObject(name, runTime, dict),
+    binModelPtr_(nullptr)
+{
+    if (readFields)
+    {
+        read(dict);
+    }
+}
+
+
+Foam::functionObjects::binField::binField
+(
+    const word& name,
+    const objectRegistry& obr,
+    const dictionary& dict,
+    bool readFields
+)
+:
+    fvMeshFunctionObject(name, obr, dict),
+    binModelPtr_(nullptr)
+{
+    if (readFields)
+    {
+        read(dict);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::binField::read(const dictionary& dict)
+{
+    if (!fvMeshFunctionObject::read(dict))
+    {
+        return false;
+    }
+
+    Info<< type() << " " << name() << ":" << endl;
+
+    binModelPtr_.reset(binModel::New(dict, mesh_, name()));
+
+    return true;
+}
+
+
+bool Foam::functionObjects::binField::execute()
+{
+    Log << type() << " " << name() << ":" << nl
+        << "    Calculating bins" << nl << endl;
+
+    binModelPtr_->apply();
+
+    return true;
+}
+
+
+bool Foam::functionObjects::binField::write()
+{
+    return true;
+}
+
+
+void Foam::functionObjects::binField::updateMesh(const mapPolyMesh& mpm)
+{
+    binModelPtr_->updateMesh(mpm);
+}
+
+
+void Foam::functionObjects::binField::movePoints(const polyMesh& mesh)
+{
+    binModelPtr_->movePoints(mesh);
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/binField/binField.H b/src/functionObjects/field/binField/binField.H
new file mode 100644
index 0000000000000000000000000000000000000000..1fd57a1aedae8ed3f0142b0b28c12a2508680849
--- /dev/null
+++ b/src/functionObjects/field/binField/binField.H
@@ -0,0 +1,221 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 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::binField
+
+Group
+    grpFieldFunctionObjects
+
+Description
+    Calculates binned data, where specified patches are divided into
+    segments according to various input bin characteristics, so that
+    spatially-localised information can be output for each segment.
+
+    Operands:
+    \table
+      Operand      | Type                          | Location
+      input        | vol\<Type\>Field(s)          <!--
+               --> | \<time\>/\<inpField\>s
+      output file  | dat                           <!--
+               --> | 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
+    binField1
+    {
+        // Mandatory entries
+        type                    binField;
+        libs                    (fieldFunctionObjects);
+        binModel                <word>;
+        fields                  (<wordHashSet>);
+        patches                 (<wordRes>);
+        binData
+        {
+            // Entries of the chosen binModel
+        }
+
+        // Optional entries
+        cellZones               (<wordRes>);
+        decomposePatchValues    <bool>;
+
+        // Conditional optional entries
+
+            // Option-1, i.e. general coordinate system specification
+            coordinateSystem
+            {
+                type            cartesian;
+                origin          (0 0 0);
+                rotation
+                {
+                    type        axes;
+                    e3          (0 0 1);
+                    e1          (1 0 0); // (e1, e2) or (e2, e3) or (e3, e1)
+                }
+            }
+
+            // Option-2, i.e. the centre of rotation
+            // by inherently using e3=(0 0 1) and e1=(1 0 0)
+            CofR                (0 0 0);
+
+        // Inherited entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                           | Type   | Reqd | Deflt
+      type      | Type name: binField                   | word   | yes  | -
+      libs      | Library name: fieldFunctionObjects    | word   | yes  | -
+      binModel  | Name of the bin model                 | word   | yes  | -
+      fields    | Names of operand fields           | wordHasSet | yes  | -
+      patches   | Names of operand patches             | wordRes | yes  | -
+      binData   | Entries of the chosen bin model       | dict   | yes  | -
+      decomposePatchValues | Flag to output normal and tangential      <!--
+                --> components                          | bool   | no   | false
+      cellZones | Names of operand cell zones          | wordRes | no   | -
+      coordinateSystem | Coordinate system specifier      | dict | cndtnl | -
+      CofR    | Centre of rotation                      | vector | cndtnl | -
+    \endtable
+
+    Options for the \c binModel entry:
+    \verbatim
+      singleDirectionUniformBin    | Segments in a single direction
+      uniformBin                   | Segments in multiple directions
+    \endverbatim
+
+    The inherited entries are elaborated in:
+      - \link fvMeshFunctionObject.H \endlink
+      - \link writeFile.H \endlink
+      - \link coordinateSystem.H \endlink
+
+SourceFiles
+    binField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_functionObjects_binField_H
+#define Foam_functionObjects_binField_H
+
+#include "fvMeshFunctionObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class binModel;
+
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class binField Declaration
+\*---------------------------------------------------------------------------*/
+
+class binField
+:
+    public fvMeshFunctionObject
+{
+protected:
+
+    // Protected Data
+
+        //- Runtime-selectable bin model
+        autoPtr<binModel> binModelPtr_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("binField");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        binField
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict,
+            const bool readFields = true
+        );
+
+        //- Construct from objectRegistry and dictionary
+        binField
+        (
+            const word& name,
+            const objectRegistry& obr,
+            const dictionary& dict,
+            const bool readFields = true
+        );
+
+        //- No copy construct
+        binField(const binField&) = delete;
+
+        //- No copy assignment
+        void operator=(const binField&) = delete;
+
+
+    //- Destructor
+    virtual ~binField() = default;
+
+
+    // Member Functions
+
+        //- Read the dictionary
+        virtual bool read(const dictionary& dict);
+
+        //- Execute the function object
+        virtual bool execute();
+
+        //- Write to data files/fields and to streams
+        virtual bool write();
+
+        //- Update for changes of mesh
+        virtual void updateMesh(const mapPolyMesh& mpm);
+
+        //- Update for changes of mesh
+        virtual void movePoints(const polyMesh& mesh);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/binField/binModels/binModel/binModel.C b/src/functionObjects/field/binField/binModels/binModel/binModel.C
new file mode 100644
index 0000000000000000000000000000000000000000..5e351b078372971f04dbaa1e8f02e4ffa0a5bdd3
--- /dev/null
+++ b/src/functionObjects/field/binField/binModels/binModel/binModel.C
@@ -0,0 +1,205 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021-2022 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 "binModel.H"
+#include "fvMesh.H"
+#include "cartesianCS.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(binModel, 0);
+    defineRunTimeSelectionTable(binModel, dictionary);
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<>
+bool Foam::binModel::decomposePatchValues
+(
+    List<List<vector>>& data,
+    const label bini,
+    const vector& v,
+    const vector& n
+) const
+{
+    if (!decomposePatchValues_)
+    {
+        return false;
+    }
+
+    #ifdef FULLDEBUG
+    if (data.size() != 3)
+    {
+        FatalErrorInFunction
+            << "Inconsistent data list size - expect size 3"
+            << abort(FatalError);
+    }
+    #endif
+
+    data[1][bini] += n*(v & n);
+    data[2][bini] += v - n*(v & n);
+
+    return true;
+}
+
+
+void Foam::binModel::setCoordinateSystem
+(
+    const dictionary& dict,
+    const word& e3Name,
+    const word& e1Name
+)
+{
+    coordSysPtr_.clear();
+
+    if (dict.found(coordinateSystem::typeName_()))
+    {
+        coordSysPtr_ =
+            coordinateSystem::New
+            (
+                mesh_,
+                dict,
+                coordinateSystem::typeName_()
+            );
+
+        Info<< "Setting co-ordinate system:" << nl
+            << "    - type          : " << coordSysPtr_->name() << nl
+            << "    - origin        : " << coordSysPtr_->origin() << nl
+            << "    - e3            : " << coordSysPtr_->e3() << nl
+            << "    - e1            : " << coordSysPtr_->e1() << endl;
+    }
+    else if (dict.found("CofR"))
+    {
+        const vector origin(dict.get<point>("CofR"));
+
+        const vector e3
+        (
+            e3Name == word::null
+          ? vector(0, 0, 1)
+          : dict.get<vector>(e3Name)
+        );
+
+        const vector e1
+        (
+            e1Name == word::null
+          ? vector(1, 0, 0)
+          : dict.get<vector>(e1Name)
+        );
+
+        coordSysPtr_.reset(new coordSystem::cartesian(origin, e3, e1));
+    }
+    else
+    {
+        coordSysPtr_.reset(new coordSystem::cartesian(dict));
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::binModel::binModel
+(
+    const dictionary& dict,
+    const fvMesh& mesh,
+    const word& outputPrefix
+)
+:
+    writeFile(mesh, outputPrefix),
+    mesh_(mesh),
+    decomposePatchValues_(false),
+    cumulative_(false),
+    coordSysPtr_(),
+    nBin_(1),
+    patchSet_(),
+    fieldNames_(),
+    cellZoneIDs_(),
+    filePtrs_()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::binModel::read(const dictionary& dict)
+{
+    patchSet_ = mesh_.boundaryMesh().patchSet(dict.get<wordRes>("patches"));
+    fieldNames_ = dict.get<wordHashSet>("fields").sortedToc();
+
+    if (dict.found("cellZones"))
+    {
+        DynamicList<label> zoneIDs;
+        DynamicList<wordRe> czUnmatched;
+        for (const auto& cz : dict.get<wordRes>("cellZones"))
+        {
+            const labelList czi(mesh_.cellZones().indices(cz));
+
+            if (czi.empty())
+            {
+                czUnmatched.append(cz);
+            }
+            else
+            {
+                zoneIDs.append(czi);
+            }
+        }
+
+        if (czUnmatched.size())
+        {
+            WarningInFunction
+                << "Unable to find zone(s): " << czUnmatched << nl
+                << "Valid cellZones are : " << mesh_.cellZones().sortedNames()
+                << endl;
+        }
+
+        cellZoneIDs_.transfer(zoneIDs);
+    }
+
+    decomposePatchValues_ = dict.get<bool>("decomposePatchValues");
+
+    filePtrs_.setSize(fieldNames_.size());
+    forAll(filePtrs_, i)
+    {
+        filePtrs_.set(i, createFile(fieldNames_[i] + "Bin"));
+    }
+
+    setCoordinateSystem(dict);
+
+    return true;
+}
+
+
+void Foam::binModel::updateMesh(const mapPolyMesh& mpm)
+{}
+
+
+void Foam::binModel::movePoints(const polyMesh& mesh)
+{}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/binField/binModels/binModel/binModel.H b/src/functionObjects/field/binField/binModels/binModel/binModel.H
new file mode 100644
index 0000000000000000000000000000000000000000..3583a0a97310590b4af8c2db82b776d53f7c3892
--- /dev/null
+++ b/src/functionObjects/field/binField/binModels/binModel/binModel.H
@@ -0,0 +1,243 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021-2022 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::binModel
+
+Description
+    Base class for bin models to handle general bin characteristics.
+
+SourceFiles
+    binModel.C
+    binModelNew.C
+    binModelTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_binModel_H
+#define Foam_binModel_H
+
+#include "dictionary.H"
+#include "HashSet.H"
+#include "volFields.H"
+#include "runTimeSelectionTables.H"
+#include "OFstream.H"
+#include "coordinateSystem.H"
+#include "writeFile.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class fvMesh;
+
+/*---------------------------------------------------------------------------*\
+                            Class binModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class binModel
+:
+    public functionObjects::writeFile
+{
+protected:
+
+    // Protected Data
+
+        //- Reference to the mesh
+        const fvMesh& mesh_;
+
+        //- Decompose patch values into normal and tangential components
+        bool decomposePatchValues_;
+
+        //- Flag to accumulate bin data with
+        //- increasing distance in binning direction
+        bool cumulative_;
+
+        //- Local coordinate system of bins
+        autoPtr<coordinateSystem> coordSysPtr_;
+
+        //- Total number of bins
+        label nBin_;
+
+        //- Indices of operand patches
+        labelHashSet patchSet_;
+
+        //- Names of operand fields
+        wordList fieldNames_;
+
+        //- Indices of operand cell zones
+        labelList cellZoneIDs_;
+
+        //- List of file pointers; 1 file per field
+        PtrList<OFstream> filePtrs_;
+
+
+    // Protected Member Functions
+
+        //- Set the co-ordinate system from dictionary and axes names
+        void setCoordinateSystem
+        (
+            const dictionary& dict,
+            const word& e3Name = word::null,
+            const word& e1Name = word::null
+        );
+
+        //- Helper function to decompose patch values
+        //- into normal and tangential components
+        template<class Type>
+        bool decomposePatchValues
+        (
+            List<List<Type>>& data,
+            const label bini,
+            const Type& v,
+            const vector& n
+        ) const;
+
+        //- Helper function to construct a string description for a given type
+        template<class Type>
+        string writeComponents(const word& stem) const;
+
+        //- Write binned data to stream
+        template<class Type>
+        void writeBinnedData
+        (
+            List<List<Type>>& data,
+            Ostream& os
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("binModel");
+
+
+    // Declare runtime constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            binModel,
+            dictionary,
+            (
+                const dictionary& dict,
+                const fvMesh& mesh,
+                const word& outputPrefix
+            ),
+            (dict, mesh, outputPrefix)
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected bin model
+        static autoPtr<binModel> New
+        (
+            const dictionary& dict,
+            const fvMesh& mesh,
+            const word& outputPrefix
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        binModel
+        (
+            const dictionary& dict,
+            const fvMesh& mesh,
+            const word& outputPrefix
+        );
+
+        //- No copy construct
+        binModel(const binModel&) = delete;
+
+        //- No copy assignment
+        void operator=(const binModel&) = delete;
+
+
+    //- Destructor
+    virtual ~binModel() = default;
+
+
+    // Member Functions
+
+        //- Read the dictionary
+        virtual bool read(const dictionary& dict);
+
+
+    // Access
+
+        //- Return the total number of bins
+        label nBin() const noexcept
+        {
+            return nBin_;
+        };
+
+
+    // Evaluation
+
+        //- Initialise bin properties
+        virtual void initialise() = 0;
+
+        //- Apply bins
+        virtual void apply() = 0;
+
+        //- Update for changes of mesh
+        virtual void updateMesh(const mapPolyMesh& mpm);
+
+        //- Update for changes of mesh
+        virtual void movePoints(const polyMesh& mesh);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<>
+bool binModel::decomposePatchValues
+(
+    List<List<vector>>& data,
+    const label bini,
+    const vector& v,
+    const vector& n
+) const;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "binModelTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/binField/binModels/binModel/binModelNew.C b/src/functionObjects/field/binField/binModels/binModel/binModelNew.C
new file mode 100644
index 0000000000000000000000000000000000000000..2447b0991982aa05c01f152dd5366959cb00c346
--- /dev/null
+++ b/src/functionObjects/field/binField/binModels/binModel/binModelNew.C
@@ -0,0 +1,59 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021-2022 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 "binModel.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::binModel> Foam::binModel::New
+(
+    const dictionary& dict,
+    const fvMesh& mesh,
+    const word& outputPrefix
+)
+{
+    word modelType(dict.get<word>("binModel"));
+
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
+
+    if (!cstrIter.found())
+    {
+        FatalIOErrorInLookup
+        (
+            dict,
+            "binModel",
+            modelType,
+            *dictionaryConstructorTablePtr_
+        ) << exit(FatalIOError);
+    }
+
+    return autoPtr<binModel>(cstrIter()(dict, mesh, outputPrefix));
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/binField/binModels/binModel/binModelTemplates.C b/src/functionObjects/field/binField/binModels/binModel/binModelTemplates.C
new file mode 100644
index 0000000000000000000000000000000000000000..fdd86ffaa45fa2da087f52a256488066b8f2d220
--- /dev/null
+++ b/src/functionObjects/field/binField/binModels/binModel/binModelTemplates.C
@@ -0,0 +1,102 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021-2022 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+bool Foam::binModel::decomposePatchValues
+(
+    List<List<Type>>& data,
+    const label bini,
+    const Type& v,
+    const vector& n
+) const
+{
+    return decomposePatchValues_;
+}
+
+
+template<class Type>
+Foam::string Foam::binModel::writeComponents(const word& stem) const
+{
+    if (pTraits<Type>::nComponents == 1)
+    {
+        return stem;
+    }
+
+    string result = "";
+    for (label cmpt = 0; cmpt < pTraits<Type>::nComponents; ++cmpt)
+    {
+        if (cmpt) result += " ";
+        result += stem + "_" + word(pTraits<Type>::componentNames[cmpt]);
+    }
+    return result;
+};
+
+
+template<class Type>
+void Foam::binModel::writeBinnedData
+(
+    List<List<Type>>& data,
+    Ostream& os
+) const
+{
+    if (cumulative_)
+    {
+        for (auto& datai : data)
+        {
+            for (label bini = 1; bini < nBin_; ++bini)
+            {
+                datai[bini] += datai[bini-1];
+            }
+        }
+    }
+
+    writeCurrentTime(os);
+
+    for (label bini = 0; bini < nBin_; ++bini)
+    {
+        Type total = Zero;
+
+        for (label i = 0; i < data.size(); ++i)
+        {
+            total += data[i][bini];
+        }
+
+        writeValue(os, total);
+
+        for (label i = 0; i < data.size(); ++i)
+        {
+            writeValue(os, data[i][bini]);
+        }
+    }
+
+    os  << endl;
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/binField/binModels/singleDirectionUniformBin/singleDirectionUniformBin.C b/src/functionObjects/field/binField/binModels/singleDirectionUniformBin/singleDirectionUniformBin.C
new file mode 100644
index 0000000000000000000000000000000000000000..2140504cc04c00209a0c916f63d894e9b01f0fd6
--- /dev/null
+++ b/src/functionObjects/field/binField/binModels/singleDirectionUniformBin/singleDirectionUniformBin.C
@@ -0,0 +1,186 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021-2022 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 "singleDirectionUniformBin.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace binModels
+{
+    defineTypeNameAndDebug(singleDirectionUniformBin, 0);
+    addToRunTimeSelectionTable(binModel, singleDirectionUniformBin, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void Foam::binModels::singleDirectionUniformBin::initialise()
+{
+    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+
+    // Determine extents of patches in a given direction
+    scalar geomMin = GREAT;
+    scalar geomMax = -GREAT;
+    for (const label patchi : patchSet_)
+    {
+        const polyPatch& pp = pbm[patchi];
+        const scalarField d(pp.faceCentres() & binDir_);
+        geomMin = min(min(d), geomMin);
+        geomMax = max(max(d), geomMax);
+    }
+
+    for (const label zonei : cellZoneIDs_)
+    {
+        const cellZone& cZone = mesh_.cellZones()[zonei];
+        const vectorField cz(mesh_.C(), cZone);
+        const scalarField d(cz & binDir_);
+
+        geomMin = min(min(d), geomMin);
+        geomMax = max(max(d), geomMax);
+    }
+
+    reduce(geomMin, minOp<scalar>());
+    reduce(geomMax, maxOp<scalar>());
+
+    // Slightly boost max so that region of interest is fully within bounds
+    geomMax = 1.0001*(geomMax - geomMin) + geomMin;
+
+    // Use geometry limits if not specified by the user
+    if (binMin_ == GREAT) binMin_ = geomMin;
+    if (binMax_ == GREAT) binMax_ = geomMax;
+
+    binDx_ = (binMax_ - binMin_)/scalar(nBin_);
+
+    if (binDx_ <= 0)
+    {
+        FatalErrorInFunction
+            << "Max bound must be greater than min bound" << nl
+            << "    d           = " << binDx_ << nl
+            << "    min         = " << binMin_ << nl
+            << "    max         = " << binMax_ << nl
+            << exit(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::binModels::singleDirectionUniformBin::singleDirectionUniformBin
+(
+    const dictionary& dict,
+    const fvMesh& mesh,
+    const word& outputPrefix
+)
+:
+    binModel(dict, mesh, outputPrefix),
+    binDx_(0),
+    binMin_(GREAT),
+    binMax_(GREAT),
+    binDir_(Zero)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::binModels::singleDirectionUniformBin::read(const dictionary& dict)
+{
+    if (!binModel::read(dict))
+    {
+        return false;
+    }
+
+    Info<< "    Activating a set of single-direction bins" << endl;
+
+    const dictionary& binDict = dict.subDict("binData");
+
+    nBin_ = binDict.getCheck<label>("nBin", labelMinMax::ge(1));
+
+    Info<< "    Employing " << nBin_ << " bins" << endl;
+    if (binDict.readIfPresent("min", binMin_))
+    {
+        Info<< "    - min        : " << binMin_ << endl;
+    }
+    if (binDict.readIfPresent("max", binMax_))
+    {
+        Info<< "    - max        : " << binMax_ << endl;
+    }
+
+    cumulative_ = binDict.getOrDefault<bool>("cumulative", false);
+    Info<< "    - cumulative    : " << cumulative_ << endl;
+    Info<< "    - decomposePatchValues    : " << decomposePatchValues_ << endl;
+
+    binDir_ = binDict.get<vector>("direction");
+    binDir_.normalise();
+
+    if (mag(binDir_) == 0)
+    {
+        FatalIOErrorInFunction(dict)
+            << "Input direction should not be zero valued" << nl
+            << "    direction = " << binDir_ << nl
+            << exit(FatalIOError);
+    }
+
+    Info<< "    - direction     : " << binDir_ << nl << endl;
+
+    initialise();
+
+    return true;
+}
+
+
+void Foam::binModels::singleDirectionUniformBin::apply()
+{
+    forAll(fieldNames_, i)
+    {
+        const bool ok =
+            processField<scalar>(i)
+         || processField<vector>(i)
+         || processField<sphericalTensor>(i)
+         || processField<symmTensor>(i)
+         || processField<tensor>(i);
+
+        if (!ok)
+        {
+            WarningInFunction
+                << "Unable to find field " << fieldNames_[i]
+                << ". Avaliable objects are "
+                << mesh_.objectRegistry::sortedToc()
+                << endl;
+        }
+    }
+
+    writtenHeader_ = true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/binField/binModels/singleDirectionUniformBin/singleDirectionUniformBin.H b/src/functionObjects/field/binField/binModels/singleDirectionUniformBin/singleDirectionUniformBin.H
new file mode 100644
index 0000000000000000000000000000000000000000..1e2a631b11b0f8153f88231979db0df25324239f
--- /dev/null
+++ b/src/functionObjects/field/binField/binModels/singleDirectionUniformBin/singleDirectionUniformBin.H
@@ -0,0 +1,186 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021-2022 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::binModels::singleDirectionUniformBin
+
+Description
+    Calculates binned data in a specified direction.
+
+    For example, a 10cm-long patch extending only in the x-direction
+    can be binned into 5 bins in the same direction, so that
+    local information can be output for each 2cm-long segment.
+
+Usage
+    Minimal example by using \c system/controlDict.functions:
+    \verbatim
+    binField1
+    {
+        // Other binField entries
+        ...
+
+        // Mandatory entries
+        binModel          singleDirectionUniformBin;
+
+        binData
+        {
+            // Mandatory entries
+            nBin          <label>;
+            direction     <vector>;
+
+            // Optional entries
+            cumulative    <bool>;
+            min           <scalar>;
+            max           <scalar>;
+        }
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                           | Type   | Reqd | Deflt
+      binModel  | Type name: singleDirectionUniformBin  | word   | yes  | -
+      binData   | Entries of the chosen bin model       | dict   | yes  | -
+      nBin      | Number of bins in binning direction   | label  | yes  | -
+      direction | Binning direction                     | vector | yes  | -
+      cumulative | Flag to bin data accumulated with increasing distance <!--
+                --> in binning direction                | bool   | no   | false
+      min       | Min-bound in the binning direction with respect to   <!--
+                --> the global coordinate system's origin | scalar | no | GREAT
+      max       | Max-bound in the binning direction with respect to   <!--
+                --> the global coordinate system's origin | scalar | no | GREAT
+    \endtable
+
+SourceFiles
+    singleDirectionUniformBin.C
+    singleDirectionUniformBinTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_binModels_singleDirectionUniformBin_H
+#define Foam_binModels_singleDirectionUniformBin_H
+
+#include "binModel.H"
+#include "writeFile.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace binModels
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class singleDirectionUniformBin Declaration
+\*---------------------------------------------------------------------------*/
+
+class singleDirectionUniformBin
+:
+    public binModel
+{
+protected:
+
+    // Protected Data
+
+        //- Distance between bin divisions
+        scalar binDx_;
+
+        //- Minimum bin bound
+        scalar binMin_;
+
+        //- Maximum bin bound
+        scalar binMax_;
+
+        //- Binning direction
+        vector binDir_;
+
+
+    // Protected Member Functions
+
+        //- Write header for a binned-data file
+        template<class Type>
+        void writeFileHeader(OFstream& os) const;
+
+        //- Initialise bin properties
+        virtual void initialise();
+
+        //- Apply the binning to field fieldi
+        template<class Type>
+        bool processField(const label fieldi);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("singleDirectionUniformBin");
+
+
+    // Constructors
+
+        //- Construct from components
+        singleDirectionUniformBin
+        (
+            const dictionary& dict,
+            const fvMesh& mesh,
+            const word& outputPrefix
+        );
+
+        //- No copy construct
+        singleDirectionUniformBin(const singleDirectionUniformBin&) = delete;
+
+        //- No copy assignment
+        void operator=(const singleDirectionUniformBin&) = delete;
+
+
+    //- Destructor
+    virtual ~singleDirectionUniformBin() = default;
+
+
+    // Member Functions
+
+        //- Read the dictionary
+        virtual bool read(const dictionary& dict);
+
+        //- Apply bins
+        virtual void apply();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace binModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "singleDirectionUniformBinTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/binField/binModels/singleDirectionUniformBin/singleDirectionUniformBinTemplates.C b/src/functionObjects/field/binField/binModels/singleDirectionUniformBin/singleDirectionUniformBinTemplates.C
new file mode 100644
index 0000000000000000000000000000000000000000..5fedd25a11483a38da197ee0f3afecd3c8920000
--- /dev/null
+++ b/src/functionObjects/field/binField/binModels/singleDirectionUniformBin/singleDirectionUniformBinTemplates.C
@@ -0,0 +1,195 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021-2022 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+
+template<class Type>
+void Foam::binModels::singleDirectionUniformBin::writeFileHeader
+(
+    OFstream& os
+) const
+{
+    writeHeaderValue(os, "bins", nBin_);
+    writeHeaderValue(os, "start", binMin_);
+    writeHeaderValue(os, "end", binMax_);
+    writeHeaderValue(os, "delta", binDx_);
+    writeHeaderValue(os, "direction", binDir_);
+
+    // Compute and print bin end points in the binning direction
+    vectorField binPoints(nBin_);
+    writeCommented(os, "x co-ords  :");
+    forAll(binPoints, pointi)
+    {
+        binPoints[pointi] = (binMin_ + (pointi + 1)*binDx_)*binDir_;
+        os  << tab << binPoints[pointi].x();
+    }
+    os  << nl;
+
+    writeCommented(os, "y co-ords  :");
+    forAll(binPoints, pointi)
+    {
+        os  << tab << binPoints[pointi].y();
+    }
+    os  << nl;
+
+    writeCommented(os, "z co-ords  :");
+    forAll(binPoints, pointi)
+    {
+        os  << tab << binPoints[pointi].z();
+    }
+    os  << nl;
+
+    writeHeader(os, "");
+    writeCommented(os, "Time");
+
+    for (label i = 0; i < nBin_; ++i)
+    {
+        const word ibin("_" + Foam::name(i));
+        writeTabbed(os, writeComponents<Type>("total" + ibin));
+        writeTabbed(os, writeComponents<Type>("internal" + ibin));
+
+        if (decomposePatchValues_)
+        {
+            writeTabbed(os, writeComponents<Type>("normal" + ibin));
+            writeTabbed(os, writeComponents<Type>("tangenial" + ibin));
+        }
+        else
+        {
+            writeTabbed(os, writeComponents<Type>("patch" + ibin));
+        }
+
+    }
+
+    os  << endl;
+}
+
+
+template<class Type>
+bool Foam::binModels::singleDirectionUniformBin::processField
+(
+    const label fieldi
+)
+{
+    const word& fieldName = fieldNames_[fieldi];
+
+    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
+
+    const VolFieldType* fieldPtr = mesh_.findObject<VolFieldType>(fieldName);
+
+    if (!fieldPtr)
+    {
+        return false;
+    }
+
+    if (Pstream::master() && !writtenHeader_)
+    {
+        writeFileHeader<Type>(filePtrs_[fieldi]);
+    }
+
+    const VolFieldType& fld = *fieldPtr;
+
+    // Total number of fields
+    //
+    // 0: internal
+    // 1: patch total
+    //
+    // OR
+    //
+    // 0: internal
+    // 1: patch normal
+    // 2: patch tangential
+    label nField = 2;
+    if (decomposePatchValues_)
+    {
+        nField += 1;
+    }
+
+    List<List<Type>> data(nField);
+    for (auto& binList : data)
+    {
+        binList.setSize(nBin_, Zero);
+    }
+
+    for (const label zonei : cellZoneIDs_)
+    {
+        const cellZone& cZone = mesh_.cellZones()[zonei];
+
+        for (const label celli : cZone)
+        {
+            const scalar dd = mesh_.C()[celli] & binDir_;
+
+            if (dd < binMin_ || dd > binMax_)
+            {
+                continue;
+            }
+
+            // Find the bin division corresponding to the cell
+            const label bini =
+                min(max(floor((dd - binMin_)/binDx_), 0), nBin_ - 1);
+
+            data[0][bini] += fld[celli];
+        }
+    }
+
+    forAllIters(patchSet_, iter)
+    {
+        const label patchi = iter();
+        const polyPatch& pp = mesh_.boundaryMesh()[patchi];
+        const vectorField np(mesh_.boundary()[patchi].nf());
+
+        const scalarField dd(pp.faceCentres() & binDir_);
+
+        forAll(dd, facei)
+        {
+            // Avoid faces outside of the bin
+            if (dd[facei] < binMin_ || dd[facei] > binMax_)
+            {
+                continue;
+            }
+
+            // Find the bin division corresponding to the face
+            const label bini =
+                min(max(floor((dd[facei] - binMin_)/binDx_), 0), nBin_ - 1);
+
+            const Type& v = fld.boundaryField()[patchi][facei];
+
+            if (!decomposePatchValues(data, bini, v, np[facei]))
+            {
+                data[1][bini] += v;
+            }
+        }
+    }
+
+    if (Pstream::master())
+    {
+        writeBinnedData(data, filePtrs_[fieldi]);
+    }
+
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/binField/binModels/uniformBin/uniformBin.C b/src/functionObjects/field/binField/binModels/uniformBin/uniformBin.C
new file mode 100644
index 0000000000000000000000000000000000000000..77836bf690b9ec932a4a2f21117972f10db38c6c
--- /dev/null
+++ b/src/functionObjects/field/binField/binModels/uniformBin/uniformBin.C
@@ -0,0 +1,315 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021-2022 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 "uniformBin.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace binModels
+{
+    defineTypeNameAndDebug(uniformBin, 0);
+    addToRunTimeSelectionTable(binModel, uniformBin, dictionary);
+}
+}
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void Foam::binModels::uniformBin::initialise()
+{
+    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+
+    // Determine extents of patches in a given coordinate system
+    vector geomMin(GREAT, GREAT, GREAT);
+    vector geomMax(-GREAT, -GREAT, -GREAT);
+
+    for (const label patchi : patchSet_)
+    {
+        const polyPatch& pp = pbm[patchi];
+        const vectorField ppcs(coordSysPtr_->localPosition(pp.faceCentres()));
+
+        for (direction i = 0; i < vector::nComponents; ++i)
+        {
+            geomMin[i] = min(min(ppcs.component(i)), geomMin[i]);
+            geomMax[i] = max(max(ppcs.component(i)), geomMax[i]);
+        }
+    }
+
+    for (const label zonei : cellZoneIDs_)
+    {
+        const cellZone& cZone = mesh_.cellZones()[zonei];
+        const vectorField d
+        (
+            coordSysPtr_->localPosition(vectorField(mesh_.C(), cZone))
+        );
+
+        for (direction i = 0; i < vector::nComponents; ++i)
+        {
+            geomMin[i] = min(min(d.component(i)), geomMin[i]);
+            geomMax[i] = max(max(d.component(i)), geomMax[i]);
+        }
+    }
+
+    reduce(geomMin, minOp<vector>());
+    reduce(geomMax, maxOp<vector>());
+
+    for (direction i = 0; i < vector::nComponents; ++i)
+    {
+        // Slightly boost max so that region of interest is fully within bounds
+        geomMax[i] = 1.0001*(geomMax[i] - geomMin[i]) + geomMin[i];
+
+        // Use geometry limits if not specified by the user
+        if (binMinMax_[i][0] == GREAT) binMinMax_[i][0] = geomMin[i];
+        if (binMinMax_[i][1] == GREAT) binMinMax_[i][1] = geomMax[i];
+
+        if (binMinMax_[i][0] > binMinMax_[i][1])
+        {
+            FatalErrorInFunction
+                << "Max bounds must be greater than min bounds" << nl
+                << "    direction   = " << i << nl
+                << "    min         = " << binMinMax_[i][0] << nl
+                << "    max         = " << binMinMax_[i][1] << nl
+                << exit(FatalError);
+        }
+
+        //- Compute bin widths in binning directions
+        binW_[i] = (binMinMax_[i][1] - binMinMax_[i][0])/scalar(nBins_[i]);
+
+        if (binW_[i] <= 0)
+        {
+            FatalErrorInFunction
+                << "Bin widths must be greater than zero" << nl
+                << "    direction = " << i << nl
+                << "    min bound = " << binMinMax_[i][0] << nl
+                << "    max bound = " << binMinMax_[i][1] << nl
+                << "    bin width = " << binW_[i]
+                << exit(FatalError);
+        }
+    }
+
+    setBinsAddressing();
+}
+
+
+Foam::labelList Foam::binModels::uniformBin::binAddr(const vectorField& d) const
+{
+    labelList binIndices(d.size(), -1);
+
+    forAll(d, i)
+    {
+        // Avoid elements outside of the bin
+        bool faceInside = true;
+        for (direction j = 0; j < vector::nComponents; ++j)
+        {
+            if (d[i][j] < binMinMax_[j][0] || d[i][j] > binMinMax_[j][1])
+            {
+                faceInside = false;
+                break;
+            }
+        }
+
+        if (faceInside)
+        {
+            // Find the bin division corresponding to the element
+            Vector<label> n(Zero);
+            for (direction j = 0; j < vector::nComponents; ++j)
+            {
+                n[j] = floor((d[i][j] - binMinMax_[j][0])/binW_[j]);
+                n[j] = min(max(n[j], 0), nBins_[j] - 1);
+            }
+
+            // Order: (e1, e2, e3), the first varies the fastest
+            binIndices[i] = n[0] + nBins_[0]*n[1] + nBins_[0]*nBins_[1]*n[2];
+        }
+        else
+        {
+            binIndices[i] = -1;
+        }
+    }
+
+    return binIndices;
+}
+
+
+void Foam::binModels::uniformBin::setBinsAddressing()
+{
+    faceToBin_.setSize(mesh_.nBoundaryFaces());
+    faceToBin_ = -1;
+
+    forAllIters(patchSet_, iter)
+    {
+        const polyPatch& pp = mesh_.boundaryMesh()[iter()];
+        const label i0 = pp.start() - mesh_.nInternalFaces();
+
+        SubList<label>(faceToBin_, pp.size(), i0) =
+            binAddr(coordSysPtr_->localPosition(pp.faceCentres()));
+    }
+
+    cellToBin_.setSize(mesh_.nCells());
+    cellToBin_ = -1;
+
+    for (const label zonei : cellZoneIDs_)
+    {
+        const cellZone& cZone = mesh_.cellZones()[zonei];
+        labelList bins
+        (
+            binAddr(coordSysPtr_->localPosition(vectorField(mesh_.C(), cZone)))
+        );
+
+        forAll(cZone, i)
+        {
+            const label celli = cZone[i];
+            cellToBin_[celli] = bins[i];
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::binModels::uniformBin::uniformBin
+(
+    const dictionary& dict,
+    const fvMesh& mesh,
+    const word& outputPrefix
+)
+:
+    binModel(dict, mesh, outputPrefix),
+    nBins_(Zero),
+    binW_(Zero),
+    binMinMax_
+    (
+        vector2D(GREAT, GREAT),
+        vector2D(GREAT, GREAT),
+        vector2D(GREAT, GREAT)
+    )
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::binModels::uniformBin::read(const dictionary& dict)
+{
+    if (!binModel::read(dict))
+    {
+        return false;
+    }
+
+    Info<< "    Activating a set of uniform bins" << endl;
+
+    const dictionary& binDict = dict.subDict("binData");
+
+    nBins_ = binDict.get<Vector<label>>("nBin");
+
+    for (const label n : nBins_)
+    {
+        nBin_ *= n;
+    }
+
+    if (nBin_ <= 0)
+    {
+        FatalIOErrorInFunction(binDict)
+            << "Number of bins must be greater than zero" << nl
+            << "    e1 bins = " << nBins_[0] << nl
+            << "    e2 bins = " << nBins_[1] << nl
+            << "    e3 bins = " << nBins_[2]
+            << exit(FatalIOError);
+    }
+
+    Info<< "    - Employing:" << nl
+        << "        " << nBins_[0] << " e1 bins," << nl
+        << "        " << nBins_[1] << " e2 bins," << nl
+        << "        " << nBins_[2] << " e3 bins"
+        << endl;
+
+    cumulative_ = binDict.getOrDefault<bool>("cumulative", false);
+    Info<< "    - cumulative    : " << cumulative_ << endl;
+    Info<< "    - decomposePatchValues    : " << decomposePatchValues_ << endl;
+
+    if (binDict.found("minMax"))
+    {
+        const dictionary& minMaxDict = binDict.subDict("minMax");
+
+        for (direction i = 0; i < vector::nComponents; ++i)
+        {
+            const word ei("e" + Foam::name(i));
+
+            if (minMaxDict.readIfPresent(ei, binMinMax_[i]))
+            {
+                Info<< "    - " << ei << " min        : "
+                    << binMinMax_[i][0] << nl
+                    << "    - " << ei << " max        : "
+                    << binMinMax_[i][1] << endl;
+            }
+        }
+    }
+    Info<< endl;
+
+    initialise();
+
+    return true;
+}
+
+
+void Foam::binModels::uniformBin::apply()
+{
+    forAll(fieldNames_, i)
+    {
+        const bool ok =
+            processField<scalar>(i)
+         || processField<vector>(i)
+         || processField<sphericalTensor>(i)
+         || processField<symmTensor>(i)
+         || processField<tensor>(i);
+
+        if (!ok)
+        {
+            WarningInFunction
+                << "Unable to find field " << fieldNames_[i]
+                << endl;
+        }
+    }
+
+    writtenHeader_ = true;
+}
+
+
+void Foam::binModels::uniformBin::updateMesh(const mapPolyMesh& mpm)
+{}
+
+
+void Foam::binModels::uniformBin::movePoints(const polyMesh& mesh)
+{
+    setBinsAddressing();
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/binField/binModels/uniformBin/uniformBin.H b/src/functionObjects/field/binField/binModels/uniformBin/uniformBin.H
new file mode 100644
index 0000000000000000000000000000000000000000..bbc043b6daead1b9973ed66f28a7cda82518f749
--- /dev/null
+++ b/src/functionObjects/field/binField/binModels/uniformBin/uniformBin.H
@@ -0,0 +1,222 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021-2022 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::binModels::uniformBin
+
+Description
+    Calculates binned data in multiple segments according to
+    a specified Cartesian or cylindrical coordinate system.
+
+Usage
+    Minimal example by using \c system/controlDict.functions:
+    \verbatim
+    binField1
+    {
+        // Other binField entries
+        ...
+
+        // Mandatory entries
+        binModel          uniformBin;
+
+        binData
+        {
+            // Mandatory entries
+            nBin          <Vector<label>>;
+
+            // Optional entries
+            cumulative    <bool>;
+            minMax
+            {
+                e1        (<scalar> <scalar>); // (min max);
+                e2        (<scalar> <scalar>);
+                e3        (<scalar> <scalar>);
+            }
+        }
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                           | Type   | Reqd | Deflt
+      binModel  | Type name: uniformBin                 | word   | yes  | -
+      binData   | Entries of the chosen bin model       | dict   | yes  | -
+      nBin      | Numbers of bins in specified directions | Vector\<label\> <!--
+                -->                                              | yes  | -
+      cumulative | Flag to bin data accumulated with increasing distance   <!--
+                --> in binning direction                | bool   | no   | false
+      minMax     | Min-max bounds in binning directions with respect to    <!--
+                --> the coordinateSystem's origin       | dict   | no   | -
+    \endtable
+
+Note
+  - The order of bin numbering is (e1, e2, e3), where the first
+    varies the fastest. For example, for a cylindrical bin with
+    \f$ nBin = (radial, azimuth, height) = (2, 4, 2) \f$, the bin indices
+    may look like as follows - note the counterclockwise increments:
+    \verbatim
+                           |-------------------|
+                           | 12 |         | 14 |
+                                | 11 | 13 |
+                                | 9  | 15 |
+                           | 10 |         | 16 |
+                           |-------------------|
+                          /    /         /    /
+                         /    /         /    /
+                        |-------------------|
+                        | 4  |         | 6  |
+                             | 3  | 5  |
+                             | 1  | 7  |
+                        | 2  |         | 8  |
+                        |-------------------|
+    \endverbatim
+
+SourceFiles
+    uniformBin.C
+    uniformBinTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_binModels_uniformBin_H
+#define Foam_binModels_uniformBin_H
+
+#include "binModel.H"
+#include "writeFile.H"
+#include "coordinateSystem.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace binModels
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class uniformBin Declaration
+\*---------------------------------------------------------------------------*/
+
+class uniformBin
+:
+    public binModel
+{
+protected:
+
+    // Protected Data
+
+        //- Numbers of bins in binning directions
+        Vector<label> nBins_;
+
+        //- Equidistant bin widths in binning directions
+        vector binW_;
+
+        //- Min-max bounds of bins in binning directions
+        Vector<vector2D> binMinMax_;
+
+        //- Face index to bin index addressing
+        labelList faceToBin_;
+
+        //- Cell index to bin index addressing
+        labelList cellToBin_;
+
+
+    // Protected Member Functions
+
+        //- Write header for an binned-data file
+        template<class Type>
+        void writeFileHeader(OFstream& os) const;
+
+        //- Initialise bin properties
+        virtual void initialise();
+
+        //- Return list of bin indices corresponding to positions given by d
+        virtual labelList binAddr(const vectorField& d) const;
+
+        //- Set/cache the bin addressing
+        virtual void setBinsAddressing();
+
+        //- Apply the binning to field fieldi
+        template<class Type>
+        bool processField(const label fieldi);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("uniformBin");
+
+
+    // Constructors
+
+        //- Construct from components
+        uniformBin
+        (
+            const dictionary& dict,
+            const fvMesh& mesh,
+            const word& outputPrefix
+        );
+
+        //- No copy construct
+        uniformBin(const uniformBin&) = delete;
+
+        //- No copy assignment
+        void operator=(const uniformBin&) = delete;
+
+
+    //- Destructor
+    virtual ~uniformBin() = default;
+
+
+    // Member Functions
+
+        //- Read the dictionary
+        virtual bool read(const dictionary& dict);
+
+        //- Apply bins
+        virtual void apply();
+
+        //- Update for changes of mesh
+        virtual void updateMesh(const mapPolyMesh& mpm);
+
+        //- Update for changes of mesh
+        virtual void movePoints(const polyMesh& mesh);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "uniformBinTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace binModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/binField/binModels/uniformBin/uniformBinTemplates.C b/src/functionObjects/field/binField/binModels/uniformBin/uniformBinTemplates.C
new file mode 100644
index 0000000000000000000000000000000000000000..8f7c1c103e180ba82ab9ae198343a0b982f27d57
--- /dev/null
+++ b/src/functionObjects/field/binField/binModels/uniformBin/uniformBinTemplates.C
@@ -0,0 +1,178 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021-2022 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+void Foam::binModels::uniformBin::writeFileHeader
+(
+    OFstream& os
+) const
+{
+    writeHeader(os, "bins");
+
+    const tensor& R = coordSysPtr_->R();
+    for (direction i = 0; i < vector::nComponents; ++i)
+    {
+        writeHeaderValue(os, "e" + Foam::name(i) + " bins", nBins_[i]);
+        writeHeaderValue(os, "    start", binMinMax_[i][0]);
+        writeHeaderValue(os, "    end", binMinMax_[i][1]);
+        writeHeaderValue(os, "    delta", binW_[i]);
+        writeHeaderValue(os, "    direction", R.col(i));
+    }
+    writeCommented(os, "bin end co-ordinates:");
+    os  << nl;
+
+    // Compute and print bin end points in binning directions
+    for (direction i = 0; i < vector::nComponents; ++i)
+    {
+        scalar binEnd = binMinMax_[i][0];
+
+        writeCommented(os, "e"+Foam::name(i)+" co-ords   :");
+        for (label j = 0; j < nBins_[i]; ++j)
+        {
+            binEnd += binW_[i];
+            os  << tab << binEnd;
+        }
+        os  << nl;
+    }
+
+    writeHeader(os, "");
+    writeCommented(os, "Time");
+
+    for (label i = 0; i < nBin_; ++i)
+    {
+        const word ibin(Foam::name(i) + ':');
+        writeTabbed(os, writeComponents<Type>("total" + ibin));
+        writeTabbed(os, writeComponents<Type>("internal" + ibin));
+
+        if (decomposePatchValues_)
+        {
+            writeTabbed(os, writeComponents<Type>("normal" + ibin));
+            writeTabbed(os, writeComponents<Type>("tangential" + ibin));
+        }
+        else
+        {
+            writeTabbed(os, writeComponents<Type>("patch" + ibin));
+        }
+    }
+
+    os  << endl;
+}
+
+template<class Type>
+bool Foam::binModels::uniformBin::processField(const label fieldi)
+{
+    const word& fieldName = fieldNames_[fieldi];
+
+    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
+
+    const VolFieldType* fieldPtr = mesh_.findObject<VolFieldType>(fieldName);
+
+    if (!fieldPtr)
+    {
+        return false;
+    }
+
+    if (Pstream::master() && !writtenHeader_)
+    {
+        writeFileHeader<Type>(filePtrs_[fieldi]);
+    }
+
+    const VolFieldType& fld = *fieldPtr;
+
+    // Total number of fields
+    //
+    // 0: internal
+    // 1: patch total
+    //
+    // OR
+    //
+    // 0: internal
+    // 1: patch normal
+    // 2: patch tangential
+    label nField = 2;
+    if (decomposePatchValues_)
+    {
+        nField += 1;
+    }
+
+    List<List<Type>> data(nField);
+    for (auto& binList : data)
+    {
+        binList.setSize(nBin_, Zero);
+    }
+
+    for (const label zonei : cellZoneIDs_)
+    {
+        const cellZone& cZone = mesh_.cellZones()[zonei];
+
+        for (const label celli : cZone)
+        {
+            const label bini = cellToBin_[celli];
+
+            if (bini != -1)
+            {
+                data[0][bini] += fld[celli];
+            }
+        }
+    }
+
+    forAllIters(patchSet_, iter)
+    {
+        const label patchi = iter();
+        const polyPatch& pp = mesh_.boundaryMesh()[patchi];
+        const vectorField np(mesh_.boundary()[patchi].nf());
+
+        forAll(pp, facei)
+        {
+            const label localFacei =
+                pp.start() - mesh_.nInternalFaces() + facei;
+            const label bini = faceToBin_[localFacei];
+
+            if (bini != -1)
+            {
+                const Type& v = fld.boundaryField()[patchi][facei];
+
+                if (!decomposePatchValues(data, bini, v, np[facei]))
+                {
+                    data[1][bini] += v;
+                }
+            }
+        }
+    }
+
+    if (Pstream::master())
+    {
+        writeBinnedData(data, filePtrs_[fieldi]);
+    }
+
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/forces/forceCoeffs/forceCoeffs.C b/src/functionObjects/forces/forceCoeffs/forceCoeffs.C
index 376df8038e35763a0d963edda18706091708ad0b..d9807ec9764dbb06d803c91bbf41d605219cf9db 100644
--- a/src/functionObjects/forces/forceCoeffs/forceCoeffs.C
+++ b/src/functionObjects/forces/forceCoeffs/forceCoeffs.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,47 +50,121 @@ namespace functionObjects
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-void Foam::functionObjects::forceCoeffs::createFiles()
+void Foam::functionObjects::forceCoeffs::initialise()
 {
-    // Note: Only possible to create bin files after bins have been initialised
+    if (initialised_)
+    {
+        return;
+    }
+
+    initialised_ = true;
+}
+
+
+void Foam::functionObjects::forceCoeffs::reset()
+{
+    Cf_.reset();
+    Cm_.reset();
+
+    forceCoeff_ == dimensionedVector(dimless, Zero);
+    momentCoeff_ == dimensionedVector(dimless, Zero);
+}
+
+
+Foam::HashTable<Foam::functionObjects::forceCoeffs::coeffDesc>
+Foam::functionObjects::forceCoeffs::selectCoeffs() const
+{
+    HashTable<coeffDesc> coeffs(16);
+
+    coeffs.insert("Cd", coeffDesc("Drag force", "Cd", 0, 0));
+    coeffs.insert("Cs", coeffDesc("Side force", "Cs", 1, 2));
+    coeffs.insert("Cl", coeffDesc("Lift force", "Cl", 2, 1));
 
-    if (writeToFile() && !coeffFilePtr_)
+    // Add front/rear options
+    const auto frontRearCoeffs(coeffs);
+    forAllConstIters(frontRearCoeffs, iter)
     {
-        coeffFilePtr_ = createFile("coefficient");
-        writeIntegratedHeader("Coefficients", coeffFilePtr_());
+        const auto& fr = iter.val();
+        coeffs.insert(fr.frontName(), fr.front());
+        coeffs.insert(fr.rearName(), fr.rear());
+    }
 
-        if (nBin_ > 1)
-        {
-            CdBinFilePtr_ = createFile("CdBin");
-            writeBinHeader("Drag coefficient bins", CdBinFilePtr_());
+    // Add moments
+    coeffs.insert("CmRoll", coeffDesc("Roll moment", "CmRoll", 0, -1));
+    coeffs.insert("CmPitch", coeffDesc("Pitch moment", "CmPitch", 1, -1));
+    coeffs.insert("CmYaw", coeffDesc("Yaw moment", "CmYaw", 2, -1));
+
+    return coeffs;
+}
 
-            CsBinFilePtr_ = createFile("CsBin");
-            writeBinHeader("Side coefficient bins", CsBinFilePtr_());
 
-            ClBinFilePtr_ = createFile("ClBin");
-            writeBinHeader("Lift coefficient bins", ClBinFilePtr_());
+void Foam::functionObjects::forceCoeffs::calcForceCoeffs()
+{
+    // Calculate scaling factors
+    const scalar pDyn = 0.5*rhoRef_*sqr(magUInf_);
+    const dimensionedScalar forceScaling
+    (
+        dimless/dimForce,
+        scalar(1)/(Aref_*pDyn + SMALL)
+    );
 
-            CmRollBinFilePtr_ = createFile("CmRollBin");
-            writeBinHeader("Roll moment coefficient bins", CmRollBinFilePtr_());
+    const auto& coordSys = coordSysPtr_();
 
-            CmPitchBinFilePtr_ = createFile("CmPitchBin");
-            writeBinHeader("Moment coefficient bins", CmPitchBinFilePtr_());
+    // Calculate force coefficients
+    forceCoeff_ = forceScaling*force_;
 
-            CmYawBinFilePtr_ = createFile("CmYawBin");
-            writeBinHeader("Yaw moment coefficient bins", CmYawBinFilePtr_());
-        }
+    Cf_.reset
+    (
+        forceScaling.value()*coordSys.localVector(sumPatchForcesV_),
+        forceScaling.value()*coordSys.localVector(sumPatchForcesP_),
+        forceScaling.value()*coordSys.localVector(sumInternalForces_)
+    );
+}
+
+
+void Foam::functionObjects::forceCoeffs::calcMomentCoeffs()
+{
+    // Calculate scaling factors
+    const scalar pDyn = 0.5*rhoRef_*sqr(magUInf_);
+    const dimensionedScalar momentScaling
+    (
+        dimless/(dimForce*dimLength),
+        scalar(1)/(Aref_*pDyn*lRef_ + SMALL)
+    );
+
+    const auto& coordSys = coordSysPtr_();
+
+    // Calculate moment coefficients
+    momentCoeff_ = momentScaling*moment_;
+
+    Cm_.reset
+    (
+        momentScaling.value()*coordSys.localVector(sumPatchMomentsP_),
+        momentScaling.value()*coordSys.localVector(sumPatchMomentsV_),
+        momentScaling.value()*coordSys.localVector(sumInternalMoments_)
+    );
+}
+
+
+void Foam::functionObjects::forceCoeffs::createIntegratedDataFile()
+{
+    if (!coeffFilePtr_.valid())
+    {
+        coeffFilePtr_ = createFile("coefficient");
+        writeIntegratedDataFileHeader("Coefficients", coeffFilePtr_());
     }
 }
 
 
-void Foam::functionObjects::forceCoeffs::writeIntegratedHeader
+void Foam::functionObjects::forceCoeffs::writeIntegratedDataFileHeader
 (
     const word& header,
-    Ostream& os
+    OFstream& os
 ) const
 {
     const auto& coordSys = coordSysPtr_();
-    writeHeader(os, "Force coefficients");
+
+    writeHeader(os, "Force and moment coefficients");
     writeHeaderValue(os, "dragDir", coordSys.e1());
     writeHeaderValue(os, "sideDir", coordSys.e2());
     writeHeaderValue(os, "liftDir", coordSys.e3());
@@ -103,128 +177,38 @@ void Foam::functionObjects::forceCoeffs::writeIntegratedHeader
     writeHeaderValue(os, "CofR", coordSys.origin());
     writeHeader(os, "");
     writeCommented(os, "Time");
-    writeTabbed(os, "Cd");
-    writeTabbed(os, "Cs");
-    writeTabbed(os, "Cl");
-    writeTabbed(os, "CmRoll");
-    writeTabbed(os, "CmPitch");
-    writeTabbed(os, "CmYaw");
-    writeTabbed(os, "Cd(f)");
-    writeTabbed(os, "Cd(r)");
-    writeTabbed(os, "Cs(f)");
-    writeTabbed(os, "Cs(r)");
-    writeTabbed(os, "Cl(f)");
-    writeTabbed(os, "Cl(r)");
-    os  << endl;
-}
 
-
-void Foam::functionObjects::forceCoeffs::writeBinHeader
-(
-    const word& header,
-    Ostream& os
-) const
-{
-    writeHeader(os, header);
-    writeHeaderValue(os, "bins", nBin_);
-    writeHeaderValue(os, "start", binMin_);
-    writeHeaderValue(os, "delta", binDx_);
-    writeHeaderValue(os, "direction", binDir_);
-
-    vectorField binPoints(nBin_);
-    writeCommented(os, "x co-ords  :");
-    forAll(binPoints, pointi)
+    for (const auto& iter : coeffs_.sorted())
     {
-        binPoints[pointi] = (binMin_ + (pointi + 1)*binDx_)*binDir_;
-        os  << tab << binPoints[pointi].x();
-    }
-    os  << nl;
+        const auto& coeff = iter.val();
 
-    writeCommented(os, "y co-ords  :");
-    forAll(binPoints, pointi)
-    {
-        os  << tab << binPoints[pointi].y();
-    }
-    os  << nl;
+        if (!coeff.active_) continue;
 
-    writeCommented(os, "z co-ords  :");
-    forAll(binPoints, pointi)
-    {
-        os  << tab << binPoints[pointi].z();
-    }
-    os  << nl;
-
-    writeHeader(os, "");
-    writeCommented(os, "Time");
-
-    for (label j = 0; j < nBin_; ++j)
-    {
-        const word jn(Foam::name(j) + ':');
-        writeTabbed(os, jn + "total");
-        writeTabbed(os, jn + "pressure");
-        writeTabbed(os, jn + "viscous");
-
-        if (porosity_)
-        {
-            writeTabbed(os, jn + "porous");
-        }
+        writeTabbed(os, coeff.name_);
     }
 
     os  << endl;
 }
 
 
-void Foam::functionObjects::forceCoeffs::writeIntegratedData
-(
-    const word& title,
-    const List<Field<scalar>>& coeff
-) const
+void Foam::functionObjects::forceCoeffs::writeIntegratedDataFile()
 {
-    if (!log)
-    {
-        return;
-    }
-
-    const scalar pressure = sum(coeff[0]);
-    const scalar viscous = sum(coeff[1]);
-    const scalar porous = sum(coeff[2]);
-    const scalar total = pressure + viscous + porous;
-
-    Info<< "        " << title << "       : " << total << token::TAB
-        << '('
-        << "pressure: " << pressure << token::TAB
-        << "viscous: " << viscous;
-
-    if (porosity_)
-    {
-        Info<< token::TAB << "porous: " << porous;
-    }
-
-    Info<< ')' << endl;
-}
+    OFstream& os = coeffFilePtr_();
 
-
-void Foam::functionObjects::forceCoeffs::writeBinData
-(
-    const List<Field<scalar>> coeffs,
-    Ostream& os
-) const
-{
     writeCurrentTime(os);
 
-    for (label bini = 0; bini < nBin_; ++bini)
+    for (const auto& iter : coeffs_.sorted())
     {
-        scalar total = coeffs[0][bini] + coeffs[1][bini] + coeffs[2][bini];
+        const auto& coeff = iter.val();
 
-        os  << tab << total << tab << coeffs[0][bini] << tab << coeffs[1][bini];
+        if (!coeff.active_) continue;
 
-        if (porosity_)
-        {
-            os  << tab << coeffs[2][bini];
-        }
+        const vector coeffValue = coeff.value(Cf_, Cm_);
+
+        os  << tab << (coeffValue.x() + coeffValue.y() + coeffValue.z());
     }
 
-    os  << endl;
+    coeffFilePtr_() << endl;
 }
 
 
@@ -239,16 +223,39 @@ Foam::functionObjects::forceCoeffs::forceCoeffs
 )
 :
     forces(name, runTime, dict, false),
+    Cf_(),
+    Cm_(),
+    forceCoeff_
+    (
+        IOobject
+        (
+            "forceCoeff", // scopedName() is not available at ctor level
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedVector(dimless, Zero)
+    ),
+    momentCoeff_
+    (
+        IOobject
+        (
+            "momentCoeff",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedVector(dimless, Zero)
+    ),
+    coeffFilePtr_(),
     magUInf_(Zero),
     lRef_(Zero),
     Aref_(Zero),
-    coeffFilePtr_(),
-    CdBinFilePtr_(),
-    CsBinFilePtr_(),
-    ClBinFilePtr_(),
-    CmRollBinFilePtr_(),
-    CmPitchBinFilePtr_(),
-    CmYawBinFilePtr_()
+    initialised_(false)
 {
     if (readFields)
     {
@@ -263,249 +270,163 @@ Foam::functionObjects::forceCoeffs::forceCoeffs
 
 bool Foam::functionObjects::forceCoeffs::read(const dictionary& dict)
 {
-    forces::read(dict);
-
-    // Free stream velocity magnitude
-    dict.readEntry("magUInf", magUInf_);
-
-    // If case is compressible we must read rhoInf (store in rhoRef_) to
-    // calculate the reference dynamic pressure
-    // Note: for incompressible, rhoRef_ is already initialised
-    if (rhoName_ != "rhoInf")
+    if (!forces::read(dict))
     {
-        dict.readEntry("rhoInf", rhoRef_);
+        return false;
     }
 
-    // Reference length and area scales
+    initialised_ = false;
+
+    // Reference scales
+    dict.readEntry("magUInf", magUInf_);
     dict.readEntry("lRef", lRef_);
     dict.readEntry("Aref", Aref_);
 
-    if (writeFields_)
+    // If case is compressible we must read rhoInf
+    // (stored in rhoRef_) to calculate the reference dynamic pressure
+    // Note: for incompressible, rhoRef_ is already initialised to 1
+    if (rhoName_ != "rhoInf")
     {
-        volVectorField* forceCoeffPtr
-        (
-            new volVectorField
-            (
-                IOobject
-                (
-                    scopedName("forceCoeff"),
-                    mesh_.time().timeName(),
-                    mesh_,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                mesh_,
-                dimensionedVector(dimless, Zero)
-            )
-        );
-
-        mesh_.objectRegistry::store(forceCoeffPtr);
-
-        volVectorField* momentCoeffPtr
-        (
-            new volVectorField
-            (
-                IOobject
-                (
-                    scopedName("momentCoeff"),
-                    mesh_.time().timeName(),
-                    mesh_,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                mesh_,
-                dimensionedVector(dimless, Zero)
-            )
-        );
-
-        mesh_.objectRegistry::store(momentCoeffPtr);
+        dict.readEntry("rhoInf", rhoRef_);
     }
 
-    return true;
-}
+    Info<< "    magUInf: " << magUInf_ << nl
+        << "    lRef   : " << lRef_ << nl
+        << "    Aref   : " << Aref_ << nl
+        << "    rhoInf : " << rhoRef_ << endl;
 
+    if (min(magUInf_, rhoRef_) < SMALL || min(lRef_, Aref_) < SMALL)
+    {
+        FatalIOErrorInFunction(dict)
+            << "Non-zero values are required for reference scales"
+            << exit(FatalIOError);
 
-bool Foam::functionObjects::forceCoeffs::execute()
-{
-    forces::calcForcesMoment();
+        return false;
+    }
 
-    createFiles();
+    if (!dict.found("coefficients"))
+    {
+        Info<< "    Selecting all coefficients" << nl;
 
-    // Storage for pressure, viscous and porous contributions to coeffs
-    List<Field<scalar>> dragCoeffs(3);
-    List<Field<scalar>> sideCoeffs(3);
-    List<Field<scalar>> liftCoeffs(3);
-    List<Field<scalar>> rollMomentCoeffs(3);
-    List<Field<scalar>> pitchMomentCoeffs(3);
-    List<Field<scalar>> yawMomentCoeffs(3);
+        coeffs_ = selectCoeffs();
 
-    forAll(liftCoeffs, i)
-    {
-        dragCoeffs[i].setSize(nBin_);
-        sideCoeffs[i].setSize(nBin_);
-        liftCoeffs[i].setSize(nBin_);
-        rollMomentCoeffs[i].setSize(nBin_);
-        pitchMomentCoeffs[i].setSize(nBin_);
-        yawMomentCoeffs[i].setSize(nBin_);
+        for (auto& iter : coeffs_.sorted())
+        {
+            auto& coeff = iter.val();
+            coeff.active_ = true;
+            Info<< "    - " << coeff << nl;
+        }
     }
+    else
+    {
+        const wordHashSet coeffs(dict.get<wordHashSet>("coefficients"));
 
-    // Calculate coefficients
-    scalar CdTot = 0;
-    scalar CsTot = 0;
-    scalar ClTot = 0;
-    scalar CmRollTot = 0;
-    scalar CmPitchTot = 0;
-    scalar CmYawTot = 0;
-
-    const scalar pDyn = 0.5*rhoRef_*sqr(magUInf_);
+        coeffs_ = selectCoeffs();
 
-    // Avoid divide by zero in 2D cases
-    const scalar momentScaling = 1.0/(Aref_*pDyn*lRef_ + SMALL);
-    const scalar forceScaling = 1.0/(Aref_*pDyn + SMALL);
+        Info<< "    Selecting coefficients:" << nl;
 
-    const auto& coordSys = coordSysPtr_();
+        for (const word& key : coeffs)
+        {
+            auto coeffIter = coeffs_.find(key);
 
-    forAll(liftCoeffs, i)
-    {
-        const Field<vector> localForce(coordSys.localVector(force_[i]));
-        const Field<vector> localMoment(coordSys.localVector(moment_[i]));
-
-        dragCoeffs[i] = forceScaling*(localForce.component(0));
-        sideCoeffs[i] = forceScaling*(localForce.component(1));
-        liftCoeffs[i] = forceScaling*(localForce.component(2));
-        rollMomentCoeffs[i] = momentScaling*(localMoment.component(0));
-        pitchMomentCoeffs[i] = momentScaling*(localMoment.component(1));
-        yawMomentCoeffs[i] = momentScaling*(localMoment.component(2));
-
-        CdTot += sum(dragCoeffs[i]);
-        CsTot += sum(sideCoeffs[i]);
-        ClTot += sum(liftCoeffs[i]);
-        CmRollTot += sum(rollMomentCoeffs[i]);
-        CmPitchTot += sum(pitchMomentCoeffs[i]);
-        CmYawTot += sum(yawMomentCoeffs[i]);
+            if (!coeffIter.good())
+            {
+                FatalIOErrorInFunction(dict)
+                    << "Unknown coefficient type " << key
+                    << " . Valid entries are : " << coeffs_.sortedToc()
+                    << exit(FatalIOError);
+            }
+            auto& coeff = coeffIter.val();
+            coeff.active_ = true;
+            Info<< "    - " << coeff << nl;
+        }
     }
 
-    // Single contributions to the front and rear
-    const scalar CdfTot = 0.5*CdTot + CmRollTot;
-    const scalar CdrTot = 0.5*CdTot - CmRollTot;
-    const scalar CsfTot = 0.5*CsTot + CmYawTot;
-    const scalar CsrTot = 0.5*CsTot - CmYawTot;
-    const scalar ClfTot = 0.5*ClTot + CmPitchTot;
-    const scalar ClrTot = 0.5*ClTot - CmPitchTot;
+    forceCoeff_.rename(scopedName("forceCoeff"));
+    momentCoeff_.rename(scopedName("momentCoeff"));
 
-    Log << type() << " " << name() << " execute:" << nl
-        << "    Coefficients" << nl;
+    Info<< endl;
 
-    writeIntegratedData("Cd", dragCoeffs);
-    writeIntegratedData("Cs", sideCoeffs);
-    writeIntegratedData("Cl", liftCoeffs);
-    writeIntegratedData("CmRoll", rollMomentCoeffs);
-    writeIntegratedData("CmPitch", pitchMomentCoeffs);
-    writeIntegratedData("CmYaw", yawMomentCoeffs);
+    return true;
+}
 
-    Log << "        Cd(f)    : " << CdfTot << nl
-        << "        Cd(r)    : " << CdrTot << nl;
 
-    Log << "        Cs(f)    : " << CsfTot << nl
-        << "        Cs(r)    : " << CsrTot << nl;
 
-    Log << "        Cl(f)    : " << ClfTot << nl
-        << "        Cl(r)    : " << ClrTot << nl;
+bool Foam::functionObjects::forceCoeffs::execute()
+{
+    forces::calcForcesMoments();
 
-    if (writeToFile())
-    {
-        writeCurrentTime(coeffFilePtr_());
-        coeffFilePtr_()
-            << tab << CdTot << tab << CsTot << tab << ClTot
-            << tab << CmRollTot << tab << CmPitchTot << tab << CmYawTot
-            << tab << CdfTot << tab << CdrTot
-            << tab << CsfTot << tab << CsrTot
-            << tab << ClfTot << tab << ClrTot << endl;
-
-        if (nBin_ > 1)
-        {
-            if (binCumulative_)
-            {
-                forAll(liftCoeffs, i)
-                {
-                    for (label bini = 1; bini < nBin_; ++bini)
-                    {
-                        dragCoeffs[i][bini] += dragCoeffs[i][bini-1];
-                        sideCoeffs[i][bini] += sideCoeffs[i][bini-1];
-                        liftCoeffs[i][bini] += liftCoeffs[i][bini-1];
-                        rollMomentCoeffs[i][bini] +=
-                            rollMomentCoeffs[i][bini-1];
-                        pitchMomentCoeffs[i][bini] +=
-                            pitchMomentCoeffs[i][bini-1];
-                        yawMomentCoeffs[i][bini] += yawMomentCoeffs[i][bini-1];
-                    }
-                }
-            }
+    initialise();
 
-            writeBinData(dragCoeffs, CdBinFilePtr_());
-            writeBinData(sideCoeffs, CsBinFilePtr_());
-            writeBinData(liftCoeffs, ClBinFilePtr_());
-            writeBinData(rollMomentCoeffs, CmRollBinFilePtr_());
-            writeBinData(pitchMomentCoeffs, CmPitchBinFilePtr_());
-            writeBinData(yawMomentCoeffs, CmYawBinFilePtr_());
-        }
-    }
+    reset();
 
-    // Write state/results information
-    {
-        setResult("Cd", CdTot);
-        setResult("Cs", CsTot);
-        setResult("Cl", ClTot);
-        setResult("CmRoll", CmRollTot);
-        setResult("CmPitch", CmPitchTot);
-        setResult("CmYaw", CmYawTot);
-        setResult("Cd(f)", CdfTot);
-        setResult("Cd(r)", CdrTot);
-        setResult("Cs(f)", CsfTot);
-        setResult("Cs(r)", CsrTot);
-        setResult("Cl(f)", ClfTot);
-        setResult("Cl(r)", ClrTot);
-    }
+    Log << type() << " " << name() << " write:" << nl
+        << "    " << "Coefficient"
+        << tab << "Total"
+        << tab << "Pressure"
+        << tab << "Viscous"
+        << tab << "Internal"
+        << nl;
 
-    if (writeFields_)
-    {
-        const volVectorField& force =
-            lookupObject<volVectorField>(scopedName("force"));
+    calcForceCoeffs();
 
-        const volVectorField& moment =
-            lookupObject<volVectorField>(scopedName("moment"));
+    calcMomentCoeffs();
 
-        volVectorField& forceCoeff =
-            lookupObjectRef<volVectorField>(scopedName("forceCoeff"));
+    auto logValues = [](const word& name, const vector& coeff, Ostream& os)
+    {
+        os  << "    " << name << ":"
+            << tab << coeff.x() + coeff.y() + coeff.z()
+            << tab << coeff.x()
+            << tab << coeff.y()
+            << tab << coeff.z()
+            << nl;
+    };
+
+    // Always setting all results
+    for (const auto& iter : coeffs_.sorted())
+    {
+        const word& coeffName = iter.key();
+        const auto& coeff = iter.val();
 
-        volVectorField& momentCoeff =
-            lookupObjectRef<volVectorField>(scopedName("momentCoeff"));
+        // Vectors for x:pressure, y:viscous, z:internal
+        const vector coeffValue = coeff.value(Cf_, Cm_);
 
-        dimensionedScalar f0("f0", dimForce, Aref_*pDyn);
-        dimensionedScalar m0("m0", dimForce*dimLength, Aref_*lRef_*pDyn);
+        if (log && coeff.active_)
+        {
+            logValues(coeffName, coeffValue, Info);
+        }
 
-        forceCoeff == force/f0;
-        momentCoeff == moment/m0;
+        setResult(coeffName, coeffValue.x() + coeffValue.y() + coeffValue.z());
+        setResult(coeffName & "pressure", coeffValue.x());
+        setResult(coeffName & "viscous", coeffValue.y());
+        setResult(coeffName & "internal", coeffValue.z());
     }
 
+    Log  << endl;
+
     return true;
 }
 
 
 bool Foam::functionObjects::forceCoeffs::write()
 {
-    if (writeFields_)
+    if (writeToFile())
     {
-        const volVectorField& forceCoeff =
-            lookupObject<volVectorField>(scopedName("forceCoeff"));
+        Log << "    writing force and moment coefficient files." << endl;
+
+        createIntegratedDataFile();
 
-        const volVectorField& momentCoeff =
-            lookupObject<volVectorField>(scopedName("momentCoeff"));
+        writeIntegratedDataFile();
+    }
 
-        forceCoeff.write();
-        momentCoeff.write();
+    if (writeFields_)
+    {
+        forceCoeff_.write();
+        momentCoeff_.write();
     }
 
+    Log << endl;
+
     return true;
 }
 
diff --git a/src/functionObjects/forces/forceCoeffs/forceCoeffs.H b/src/functionObjects/forces/forceCoeffs/forceCoeffs.H
index 0b218f51ae3e9d61b4617074620931f931b22b05..85fbe23d2ed11ce8184d610ef429ef6d46150950 100644
--- a/src/functionObjects/forces/forceCoeffs/forceCoeffs.H
+++ b/src/functionObjects/forces/forceCoeffs/forceCoeffs.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,142 +31,207 @@ Group
     grpForcesFunctionObjects
 
 Description
-    Extends the \c forces functionObject by providing coefficients for:
-    - drag, side and lift forces (Cd, Cs, and Cl)
-    - roll, pitch and yaw moments (CmRoll, CmPitch, and CmYaw)
-    - front and rear axle force contributions (C(f) and C(r)) wherein
+    Computes force and moment coefficients over a given
+    list of patches, and optionally over given porous zones.
+    The following coefficients can be selected and output:
+    \verbatim
+      Cd      | Drag coefficient
+      Cs      | Side-force coefficient
+      Cl      | Lift coefficient
+      CmRoll  | Roll-moment coefficient
+      CmPitch | Pitch-moment coefficient
+      CmYaw   | Yaw-moment coefficient
+    \endverbatim
 
+    The force coefficients can also be optionally output
+    in terms of their front and rear axle constituents:
     \verbatim
-        Cd(f/r) = 0.5*Cd \pm CmRoll
-        Cs(f/r) = 0.5*Cs \pm CmYaw
-        Cl(f/r) = 0.5*Cl \pm CmPitch
+        Cd{f,r} = 0.5*Cd {+,-} CmRoll
+        Cs{f,r} = 0.5*Cs {+,-} CmYaw
+        Cl{f,r} = 0.5*Cl {+,-} CmPitch
     \endverbatim
+    where \c f and \c r represent front and rear axles, respectively.
 
-    The data can optionally be output into bins, defined in a given direction.
+    Force and moment coefficients are output
+    in their total and constituent components:
+    - total force and moment coefficients
+    - pressure contributions
+    - viscous contributions
+    - porous resistance contributions (optional)
 
-    The binned data provides the total and consitituent components per bin:
-    - total coefficient
-    - pressure coefficient contribution
-    - viscous coefficient contribution
-    - porous coefficient contribution
+    Force and moment coefficients can be computed and output in:
+    - the global Cartesian coordinate system (default)
+    - a user-defined Cartesian coordinate system
+
+    Operands:
+    \table
+      Operand       | Type           | Location
+      input         | -              | -
+      output file   | dat            | postProcessing/\<FO\>/\<time\>/\<file\>s
+      output field  | volVectorField | \<time\>/\<outField\>s
+    \endtable
 
-    Data is written into multiple files in the
-    postProcessing/\<functionObjectName\> directory:
-    - coefficient.dat   : integrated coefficients over all geometries
-    - CdBin.dat         : drag coefficient bins
-    - CsBin.dat         : side coefficient bins
-    - ClBin.dat         : lift coefficient bins
-    - CmRollBin.dat     : roll moment coefficient bins
-    - CmPitchBin.dat    : pitch moment coefficient bins
-    - CmYawBin.dat      : yaw moment coefficient bins
+    where \c \<file\>s:
+    \verbatim
+      coefficient.dat  | Integrated coefficients over all patches
+    \endverbatim
+
+    where \c \<outField\>s:
+    \verbatim
+      <namePrefix>:forceCoeff   | Force coefficient field
+      <namePrefix>:momentCoeff  | Moment coefficient field
+    \endverbatim
 
 Usage
-    Example of function object specification:
+    Minimal example by using \c system/controlDict.functions:
     \verbatim
     forceCoeffs1
     {
-        type        forceCoeffs;
-        libs        (forces);
-        ...
-        log         yes;
-        writeFields yes;
-        patches     (walls);
+        // Mandatory entries
+        type                forceCoeffs;
+        libs                (forces);
+        patches             (<wordRes>); // (wall1 "(wall2|wall3)");
+        magUInf             <scalar>;
+        lRef                <scalar>;
+        Aref                <scalar>;
+
+        // Optional entries
+        coefficients        (<wordHashSet>);
+        directForceDensity  <bool>;
+        porosity            <bool>;
+        writeFields         <bool>;
+        useNamePrefix       <bool>;
+
+        // Conditional mandatory entries
+
+            // Cartesian coordinate system specification when evaluating
+            // force and moment coefficients, either of the below
+
+            // Option-1, i.e. the centre of rotation
+            // by inherently using e3=(0 0 1) and e1=(1 0 0)
+            CofR                (0 0 0); // Centre of rotation
+            dragDir             (1 0 0);
+            liftDir             (0 0 1);
+
+            // Option-2, i.e. local coordinate system specification
+            origin              (0 0 0);
+            e1                  (1 0 0);
+            e3                  (0 0 1); // (e1, e2) or (e2, e3) or (e3, e1)
+
+            // Option-3, i.e. general coordinate system specification
+            coordinateSystem
+            {
+                type            cartesian;
+                origin          (0 0 0);
+                rotation
+                {
+                    type        axes;
+                    e3          (0 0 1);
+                    e1          (1 0 0); // (e1, e2) or (e2, e3) or (e3, e1)
+                }
+            }
 
-        // input keywords for directions of force/moment coefficients
-        // refer below for options, and relations
+        // Conditional optional entries
 
-        magUInf     100;
-        lRef        3.5;
-        Aref        2.2;
-        porosity    no;
+            // if directForceDensity == true
+            fD              <word>;
 
-        binData
-        {
-            nBin        20;
-            direction   (1 0 0);
-            cumulative  yes;
-        }
+            // if directForceDensity == false
+            p               <word>;
+            U               <word>;
+            rho             <word>;
+            rhoInf          <scalar>; // redundant for incompressible-flow cases
+            pRef            <scalar>;
+
+        // Inherited entries
+        ...
     }
     \endverbatim
 
-    Where the entries comprise:
-    \table
-        Property     | Description                          | Required | Default
-        type         | Type name: forceCoeffs               | yes |
-        log          | Write force data to standard output  | no  | no
-        writeFields  | Write force,moment coefficient fields | no | no
-        patches      | Patches included in the forces calculation | yes |
-        magUInf      | Free stream velocity magnitude       | yes |
-        rhoInf       | Free stream density | for compressible cases |
-        lRef         | Reference length scale for moment calculations | yes |
-        Aref         | Reference area                       | yes |
-        porosity     | Include porosity contributions       | no  | false
-    \endtable
-
-    Bin data is optional, but if the dictionary is present, the entries must
-    be defined according to following:
+    where the entries mean:
     \table
-        nBin         | Number of data bins                    | yes |
-        direction    | Direction along which bins are defined | yes |
-        cumulative   | Bin data accumulated with incresing distance | yes |
+      Property   | Description                     | Type     | Reqd | Deflt
+      type       | Type name: forceCoeffs          | word     | yes  | -
+      libs       | Library name: forces            | word     | yes  | -
+      patches    | Names of operand patches        | wordRes  | yes  | -
+      coefficients | Names of operand coefficients | wordHashSet | no | -
+      magUInf    | Reference velocity magnitude       | scalar | yes | -
+      lRef       | Reference length scale for moment calculations   <!--
+                 -->                                  | scalar | yes | -
+      Aref       | Reference area                     | scalar | yes | -
+      directForceDensity | Flag to directly supply force density    <!--
+                 -->                                  | bool   | no  | false
+      porosity   | Flag to include porosity contributions | bool | no | false
+      writeFields | Flag to write force and moment fields | bool | no | false
+      useNamePrefix | Flag to include prefix for field names | bool | no | false
+      CofR    | Centre of rotation                    | vector | cndtnl | -
+      origin  | Origin of coordinate system           | vector | cndtnl | -
+      e3      | e3 coordinate axis                    | vector | cndtnl | -
+      e1      | e1 coordinate axis                    | vector | cndtnl | -
+      coordinateSystem | Coordinate system specifier  | dictionary | cndtnl | -
+      fD      | Name of force density field  | word   | cndtnl-no | fD
+      p       | Name of pressure field       | word   | cndtnl-no | p
+      U       | Name of velocity field       | word   | cndtnl-no | U
+      rho     | Name of density field        | word   | cndtnl-no | rho
+      rhoInf  | Value of reference density   | scalar | cndtnl-yes | -
+      pRef    | Value of reference pressure  | scalar | cndtnl-no | 0
     \endtable
 
-    Input of force/moment coefficient directions:
-    - require an origin, and two orthogonal directions; the remaining orthogonal
-    direction is determined accordingly.
-    - can be added by the three options below.
-
+    Options for the \c coefficients entry:
     \verbatim
-        CofR        (0 0 0); // Centre of rotation
-        dragDir     (1 0 0);
-        liftDir     (0 0 1);
-    \endverbatim
-
-    \verbatim
-        origin (0 0 0);
-        e1     (1 0 0);
-        e3     (0 0 1); // combinations: (e1, e2) or (e2, e3) or (e3, e1)
-    \endverbatim
-
-    \verbatim
-        coordinateSystem
-        {
-            origin  (0 0 0);
-            rotation
-            {
-                type axes;
-                e1 (1 0 0);
-                e3 (0 0 1); // combinations: (e1, e2) or (e2, e3) or (e3, e1)
-            }
-        }
+      Cd      | Drag coefficient
+      Cs      | Side-force coefficient
+      Cl      | Lift coefficient
+      CmRoll  | Roll-moment coefficient
+      CmPitch | Pitch-moment coefficient
+      CmYaw   | Yaw-moment coefficient
+      Cdf     | Front-axle drag coefficient
+      Csf     | Front-axle side-force coefficient
+      Clf     | Front-axle lift coefficient
+      Cdr     | Rear-axle drag coefficient
+      Csr     | Rear-axle side-force coefficient
+      Clr     | Rear-axle lift coefficient
     \endverbatim
 
-    The default direction relations are shown below:
+    The inherited entries are elaborated in:
+      - \link functionObject.H \endlink
+      - \link writeFile.H \endlink
+      - \link coordinateSystem.H \endlink
+      - \link forces.H \endlink
+
+Note
+  - \c rhoInf is always redundant for incompressible computations.
+    That is, \c rhoInf is always equal to 1 in incompressible
+    computations no matter which input value is assigned to \c rhoInf.
+    The value of \c rhoInf is only used for compressible computations.
+  - \c writeControl and \c writeInterval entries of function
+    object do control when to output force and moment files and fields.
+  - Input for force/moment coefficient directions
+    require an origin and two orthogonal directions where
+    the remaining orthogonal direction is automatically computed.
+  - The default direction relations are shown below:
 
     \table
-        Property     | Description           | Alias | Direction
-        dragDir      | Drag direction        | e1    | (1 0 0)
-        sideDir      | Side force direction  | e2    | (0 1 0)
-        liftDir      | Lift direction        | e3    | (0 0 1)
-        rollAxis     | Roll axis             | e1    | (1 0 0)
-        pitchAxis    | Pitch axis            | e2    | (0 1 0)
-        yawAxis      | Yaw axis              | e3    | (0 0 1)
+      Property     | Description           | Alias | Direction
+      dragDir      | Drag direction        | e1    | (1 0 0)
+      sideDir      | Side force direction  | e2    | (0 1 0)
+      liftDir      | Lift direction        | e3    | (0 0 1)
+      rollAxis     | Roll axis             | e1    | (1 0 0)
+      pitchAxis    | Pitch axis            | e2    | (0 1 0)
+      yawAxis      | Yaw axis              | e3    | (0 0 1)
     \endtable
 
-See also
-    Foam::functionObject
-    Foam::functionObjects::timeControl
-    Foam::functionObjects::forces
-
 SourceFiles
     forceCoeffs.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef functionObjects_forceCoeffs_H
-#define functionObjects_forceCoeffs_H
+#ifndef Foam_functionObjects_forceCoeffs_H
+#define Foam_functionObjects_forceCoeffs_H
 
 #include "forces.H"
+#include "HashSet.H"
+#include "Enum.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -183,78 +248,261 @@ class forceCoeffs
 :
     public forces
 {
-    // Private data
+public:
 
-        // Free-stream conditions
+    // Container for force and moment components
+    class forceComponents
+    {
+        // Private Data
 
-            //- Free-stream velocity magnitude
-            scalar magUInf_;
+        //- Pressure force/moment component
+        vector pressure_;
 
+        //- Viscous force/moment component
+        vector viscous_;
 
-        // Reference scales
+        //- Internal force/moment component
+        vector internal_;
 
-            //- Reference length [m]
-            scalar lRef_;
 
-            //- Reference area [m^2]
-            scalar Aref_;
+    public:
 
+        // Constructors
 
-        // File streams
+        //- Default construct
+        forceComponents()
+        :
+            pressure_(Zero),
+            viscous_(Zero),
+            internal_(Zero)
+        {}
 
-            //- Integrated coefficients
-            autoPtr<OFstream> coeffFilePtr_;
 
-            //- Drag coefficient
-            autoPtr<OFstream> CdBinFilePtr_;
+        // Member Functions
 
-            //- Side coefficient
-            autoPtr<OFstream> CsBinFilePtr_;
+        //- Return the sum of the components
+        const vector total() const noexcept
+        {
+            return pressure_ + viscous_ + internal_;
+        }
 
-            //- Lift coefficient
-            autoPtr<OFstream> ClBinFilePtr_;
+        //- Reset the components to zero
+        void reset() noexcept
+        {
+            pressure_ = Zero;
+            viscous_ = Zero;
+            internal_ = Zero;
+        }
 
-            //- Roll moment coefficient
-            autoPtr<OFstream> CmRollBinFilePtr_;
+        //- Reset the components to given values
+        void reset
+        (
+            const vector& pressure,
+            const vector& viscous,
+            const vector& internal
+        ) noexcept
+        {
+            pressure_ = pressure;
+            viscous_ = viscous;
+            internal_ = internal;
+        }
 
-            //- Pitch moment coefficient
-            autoPtr<OFstream> CmPitchBinFilePtr_;
 
-            //- Yaw moment coefficient
-            autoPtr<OFstream> CmYawBinFilePtr_;
+        // Member Operators
 
+        //- Return components in a given direction
+        vector operator[](const label cmpt) const
+        {
+            return vector(pressure_[cmpt], viscous_[cmpt], internal_[cmpt]);
+        }
+    };
 
-    // Private Member Functions
 
-        //- No copy construct
-        forceCoeffs(const forceCoeffs&) = delete;
+    // Coefficients description
+    struct coeffDesc
+    {
+        enum splitType
+        {
+            stFront,
+            stRear,
+            stNone
+        };
+
+        string desc_;
+        word name_;
+        label c0_;
+        label c1_;
+        bool active_;
+        splitType splitType_;
+
+        // Constructors
+
+        //- Construct from components
+        coeffDesc
+        (
+            const string& description,
+            const word& name,
+            const label c0,
+            const label c1 = -1
+        )
+        :
+            desc_(description),
+            name_(name),
+            c0_(c0),
+            c1_(c1),
+            active_(false),
+            splitType_(stNone)
+        {}
+
+        //- Return name with front-name appendix
+        const word frontName() const noexcept
+        {
+            return name_ + "(f)";
+        }
 
-        //- No copy assignment
-        void operator=(const forceCoeffs&) = delete;
+        //- Return name with rear-name appendix
+        const word rearName() const noexcept
+        {
+            return name_ + "(r)";
+        }
+
+        //- Return force/moment components based on the specified split type
+        vector value(const forceComponents& f, const forceComponents& m) const
+        {
+            if (c1_ == -1)
+            {
+                return m[c0_];
+            }
+            else
+            {
+                switch (splitType_)
+                {
+                    case stFront:
+                    {
+                        return 0.5*f[c0_] + m[c1_];
+                    }
+                    case stRear:
+                    {
+                        return 0.5*f[c0_] - m[c1_];
+                    }
+                    case stNone:
+                    {
+                        return f[c0_];
+                    }
+                }
+            }
+
+            FatalErrorInFunction
+                << "Cannot determine value"
+                << abort(FatalError);
+
+            return vector::zero;
+        }
+
+        //- Return front-axle coefficient description
+        coeffDesc front() const
+        {
+            coeffDesc coeff(*this);
+            coeff.name_ = coeff.frontName();
+            coeff.desc_ = coeff.desc_ + " front";
+            coeff.splitType_ = stFront;
+            return coeff;
+        }
+
+        //- Return rear-axle coefficient description
+        coeffDesc rear() const
+        {
+            coeffDesc coeff(*this);
+            coeff.name_ = coeff.rearName();
+            coeff.desc_ = coeff.desc_ + " rear";
+            coeff.splitType_ = stRear;
+            return coeff;
+        }
+    };
+
+
+private:
+
+    // Private Data
+
+        //- Force components
+        forceComponents Cf_;
+
+        //- Moment components
+        forceComponents Cm_;
+
+        //- Force coefficient field
+        volVectorField forceCoeff_;
+
+        //- Moment coefficient field
+        volVectorField momentCoeff_;
+
+        //- Operand force and moment coefficients
+        HashTable<coeffDesc> coeffs_;
+
+
+        // File streams
+
+            //- File stream for integrated operand coefficients
+            autoPtr<OFstream> coeffFilePtr_;
+
+
+        // Reference scales
+
+            //- Reference velocity magnitude [m/s]
+            scalar magUInf_;
+
+            //- Reference length scale [m]
+            scalar lRef_;
+
+            //- Reference area [m^2]
+            scalar Aref_;
+
+
+        //- Flag of initialisation (internal)
+        bool initialised_;
 
 
 protected:
 
     // Protected Member Functions
 
-        //- Create the output files
-        void createFiles();
+        //- Initialise containers and fields
+        void initialise();
+
+        //- Reset containers and fields
+        void reset();
+
+
+    // Evaluation
 
-        //- Write header for integrated data
-        void writeIntegratedHeader(const word& header, Ostream& os) const;
+        //- Return the operand coefficients
+        HashTable<coeffDesc> selectCoeffs() const;
 
-        //- Write header for binned data
-        void writeBinHeader(const word& header, Ostream& os) const;
+        //- Calculate force coefficients
+        void calcForceCoeffs();
 
-        //- Write integrated data
-        void writeIntegratedData
+        //- Calculate moment coefficients
+        void calcMomentCoeffs();
+
+        //- Return integrated {total, pressure, viscous, porous} components
+        List<scalar> integrateData(const List<Field<scalar>>& coeff) const;
+
+
+    // I-O
+
+        //- Create the integrated-coefficient file
+        void createIntegratedDataFile();
+
+        //- Write header to the integrated-coefficient file
+        void writeIntegratedDataFileHeader
         (
-            const word& title,
-            const List<Field<scalar>>& coeff
+            const word& header,
+            OFstream& os
         ) const;
 
-        //- Write binned data
-        void writeBinData(const List<Field<scalar>> coeffs, Ostream& os) const;
+        //- Write integrated coefficients to the integrated-coefficient file
+        void writeIntegratedDataFile();
 
 
 public:
@@ -270,10 +518,16 @@ public:
         (
             const word& name,
             const Time& runTime,
-            const dictionary&,
+            const dictionary& dict,
             const bool readFields = true
         );
 
+        //- No copy construct
+        forceCoeffs(const forceCoeffs&) = delete;
+
+        //- No copy assignment
+        void operator=(const forceCoeffs&) = delete;
+
 
     //- Destructor
     virtual ~forceCoeffs() = default;
@@ -281,19 +535,27 @@ public:
 
     // Member Functions
 
-        //- Read the forces data
-        virtual bool read(const dictionary&);
+        //- Read the dictionary
+        virtual bool read(const dictionary& dict);
 
-        //- Execute
+        //- Execute the function object
         virtual bool execute();
 
-        //- Write the forces
+        //- Write to data files/fields and to streams
         virtual bool write();
 };
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+Ostream& operator<<(Ostream& os, const forceCoeffs::coeffDesc& coeff)
+{
+    os  << coeff.desc_.c_str() << ": " << coeff.name_;
+
+    return os;
+}
+
+
 } // End namespace functionObjects
 } // End namespace Foam
 
diff --git a/src/functionObjects/forces/forces/forces.C b/src/functionObjects/forces/forces/forces.C
index 2f3648012f80230b601b3bb6ff3de12e5fa2f9ce..5916aa2a61fcc6e16837ac5b9cf480dbca28b765 100644
--- a/src/functionObjects/forces/forces/forces.C
+++ b/src/functionObjects/forces/forces/forces.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,8 +31,8 @@ License
 #include "porosityModel.H"
 #include "turbulentTransportModel.H"
 #include "turbulentFluidThermoModel.H"
-#include "addToRunTimeSelectionTable.H"
 #include "cartesianCS.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -48,107 +48,6 @@ namespace functionObjects
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-void Foam::functionObjects::forces::createFiles()
-{
-    // Note: Only possible to create bin files after bins have been initialised
-
-    if (writeToFile() && !forceFilePtr_)
-    {
-        forceFilePtr_ = createFile("force");
-        writeIntegratedHeader("Force", forceFilePtr_());
-        momentFilePtr_ = createFile("moment");
-        writeIntegratedHeader("Moment", momentFilePtr_());
-
-        if (nBin_ > 1)
-        {
-            forceBinFilePtr_ = createFile("forceBin");
-            writeBinHeader("Force", forceBinFilePtr_());
-            momentBinFilePtr_ = createFile("momentBin");
-            writeBinHeader("Moment", momentBinFilePtr_());
-        }
-    }
-}
-
-
-void Foam::functionObjects::forces::writeIntegratedHeader
-(
-    const word& header,
-    Ostream& os
-) const
-{
-    writeHeader(os, header);
-    writeHeaderValue(os, "CofR", coordSysPtr_->origin());
-    writeHeader(os, "");
-    writeCommented(os, "Time");
-    writeTabbed(os, "(total_x total_y total_z)");
-    writeTabbed(os, "(pressure_x pressure_y pressure_z)");
-    writeTabbed(os, "(viscous_x viscous_y viscous_z)");
-
-    if (porosity_)
-    {
-        writeTabbed(os, "(porous_x porous_y porous_z)");
-    }
-
-    os  << endl;
-}
-
-
-void Foam::functionObjects::forces::writeBinHeader
-(
-    const word& header,
-    Ostream& os
-) const
-{
-    writeHeader(os, header + " bins");
-    writeHeaderValue(os, "bins", nBin_);
-    writeHeaderValue(os, "start", binMin_);
-    writeHeaderValue(os, "end", binMax_);
-    writeHeaderValue(os, "delta", binDx_);
-    writeHeaderValue(os, "direction", binDir_);
-
-    vectorField binPoints(nBin_);
-    writeCommented(os, "x co-ords  :");
-    forAll(binPoints, pointi)
-    {
-        binPoints[pointi] = (binMin_ + (pointi + 1)*binDx_)*binDir_;
-        os  << tab << binPoints[pointi].x();
-    }
-    os  << nl;
-
-    writeCommented(os, "y co-ords  :");
-    forAll(binPoints, pointi)
-    {
-        os  << tab << binPoints[pointi].y();
-    }
-    os  << nl;
-
-    writeCommented(os, "z co-ords  :");
-    forAll(binPoints, pointi)
-    {
-        os  << tab << binPoints[pointi].z();
-    }
-    os  << nl;
-
-    writeHeader(os, "");
-    writeCommented(os, "Time");
-
-    for (label j = 0; j < nBin_; j++)
-    {
-        const word jn(Foam::name(j) + ':');
-        os  << tab << jn << "(total_x total_y total_z)"
-            << tab << jn << "(pressure_x pressure_y pressure_z)"
-            << tab << jn << "(viscous_x viscous_y viscous_z)";
-
-        if (porosity_)
-        {
-            os  << tab << jn << "(porous_x porous_y porous_z)";
-        }
-    }
-
-    os << endl;
-}
-
-
 void Foam::functionObjects::forces::setCoordinateSystem
 (
     const dictionary& dict,
@@ -189,7 +88,6 @@ void Foam::functionObjects::forces::setCoordinateSystem
             coordSysPtr_.reset(new coordSystem::cartesian(dict));
         }
     }
-
 }
 
 
@@ -215,126 +113,38 @@ void Foam::functionObjects::forces::initialise()
         (
             !foundObject<volVectorField>(UName_)
          || !foundObject<volScalarField>(pName_)
-
         )
         {
             FatalErrorInFunction
-                << "Could not find U: " << UName_ << " or p:" << pName_
-                << " in database"
+                << "Could not find U: " << UName_
+                << " or p:" << pName_ << " in database"
                 << exit(FatalError);
         }
 
         if (rhoName_ != "rhoInf" && !foundObject<volScalarField>(rhoName_))
         {
             FatalErrorInFunction
-                << "Could not find rho:" << rhoName_
+                << "Could not find rho:" << rhoName_ << " in database"
                 << exit(FatalError);
         }
     }
 
-    initialiseBins();
-
     initialised_ = true;
 }
 
 
-void Foam::functionObjects::forces::initialiseBins()
+void Foam::functionObjects::forces::reset()
 {
-    if (nBin_ > 1)
-    {
-        const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
-
-        // Determine extents of patches
-        scalar geomMin = GREAT;
-        scalar geomMax = -GREAT;
-        for (const label patchi : patchSet_)
-        {
-            const polyPatch& pp = pbm[patchi];
-            scalarField d(pp.faceCentres() & binDir_);
-            geomMin = min(min(d), geomMin);
-            geomMax = max(max(d), geomMax);
-        }
-
-        // Include porosity
-        if (porosity_)
-        {
-            const HashTable<const porosityModel*> models =
-                obr_.lookupClass<porosityModel>();
+    sumPatchForcesP_ = Zero;
+    sumPatchForcesV_ = Zero;
+    sumPatchMomentsP_ = Zero;
+    sumPatchMomentsV_ = Zero;
 
-            const scalarField dd(mesh_.C() & binDir_);
+    sumInternalForces_ = Zero;
+    sumInternalMoments_ = Zero;
 
-            forAllConstIters(models, iter)
-            {
-                const porosityModel& pm = *iter();
-                const labelList& cellZoneIDs = pm.cellZoneIDs();
-
-                for (const label zonei : cellZoneIDs)
-                {
-                    const cellZone& cZone = mesh_.cellZones()[zonei];
-                    const scalarField d(dd, cZone);
-                    geomMin = min(min(d), geomMin);
-                    geomMax = max(max(d), geomMax);
-                }
-            }
-        }
-
-        reduce(geomMin, minOp<scalar>());
-        reduce(geomMax, maxOp<scalar>());
-
-        // Slightly boost max so that region of interest is fully within bounds
-        geomMax = 1.0001*(geomMax - geomMin) + geomMin;
-
-        // Use geometry limits if not specified by the user
-        if (binMin_ == GREAT)
-        {
-            binMin_ = geomMin;
-        }
-        if (binMax_ == GREAT)
-        {
-            binMax_ = geomMax;
-        }
-
-        binDx_ = (binMax_ - binMin_)/scalar(nBin_);
-
-        // Create the bin mid-points used for writing
-        binPoints_.setSize(nBin_);
-        forAll(binPoints_, i)
-        {
-            binPoints_[i] = (i + 0.5)*binDir_*binDx_;
-        }
-    }
-
-    // Allocate storage for forces and moments
-    forAll(force_, i)
-    {
-        force_[i].setSize(nBin_, vector::zero);
-        moment_[i].setSize(nBin_, vector::zero);
-    }
-}
-
-
-void Foam::functionObjects::forces::resetFields()
-{
-    force_[0] = Zero;
-    force_[1] = Zero;
-    force_[2] = Zero;
-
-    moment_[0] = Zero;
-    moment_[1] = Zero;
-    moment_[2] = Zero;
-
-    if (writeFields_)
-    {
-        volVectorField& force =
-            lookupObjectRef<volVectorField>(scopedName("force"));
-
-        force == dimensionedVector(force.dimensions(), Zero);
-
-        volVectorField& moment =
-            lookupObjectRef<volVectorField>(scopedName("moment"));
-
-        moment == dimensionedVector(moment.dimensions(), Zero);
-    }
+    force_ == dimensionedVector(force_.dimensions(), Zero);
+    moment_ == dimensionedVector(moment_.dimensions(), Zero);
 }
 
 
@@ -344,42 +154,45 @@ Foam::functionObjects::forces::devRhoReff() const
     typedef compressible::turbulenceModel cmpTurbModel;
     typedef incompressible::turbulenceModel icoTurbModel;
 
-    const auto& U = lookupObject<volVectorField>(UName_);
-
     if (foundObject<cmpTurbModel>(cmpTurbModel::propertiesName))
     {
-        const cmpTurbModel& turb =
+        const auto& turb =
             lookupObject<cmpTurbModel>(cmpTurbModel::propertiesName);
 
-        return turb.devRhoReff(U);
+        return turb.devRhoReff();
     }
     else if (foundObject<icoTurbModel>(icoTurbModel::propertiesName))
     {
-        const incompressible::turbulenceModel& turb =
+        const auto& turb =
             lookupObject<icoTurbModel>(icoTurbModel::propertiesName);
 
-        return rho()*turb.devReff(U);
+        return rho()*turb.devReff();
     }
     else if (foundObject<fluidThermo>(fluidThermo::dictName))
     {
-        const fluidThermo& thermo =
-            lookupObject<fluidThermo>(fluidThermo::dictName);
+        const auto& thermo = lookupObject<fluidThermo>(fluidThermo::dictName);
+
+        const auto& U = lookupObject<volVectorField>(UName_);
 
         return -thermo.mu()*dev(twoSymm(fvc::grad(U)));
     }
     else if (foundObject<transportModel>("transportProperties"))
     {
-        const transportModel& laminarT =
+        const auto& laminarT =
             lookupObject<transportModel>("transportProperties");
 
+        const auto& U = lookupObject<volVectorField>(UName_);
+
         return -rho()*laminarT.nu()*dev(twoSymm(fvc::grad(U)));
     }
     else if (foundObject<dictionary>("transportProperties"))
     {
-        const dictionary& transportProperties =
+        const auto& transportProperties =
             lookupObject<dictionary>("transportProperties");
 
-        dimensionedScalar nu("nu", dimViscosity, transportProperties);
+        const dimensionedScalar nu("nu", dimViscosity, transportProperties);
+
+        const auto& U = lookupObject<volVectorField>(UName_);
 
         return -rho()*nu*dev(twoSymm(fvc::grad(U)));
     }
@@ -398,24 +211,23 @@ Foam::tmp<Foam::volScalarField> Foam::functionObjects::forces::mu() const
 {
     if (foundObject<fluidThermo>(basicThermo::dictName))
     {
-        const fluidThermo& thermo =
-             lookupObject<fluidThermo>(basicThermo::dictName);
+        const auto& thermo = lookupObject<fluidThermo>(basicThermo::dictName);
 
         return thermo.mu();
     }
     else if (foundObject<transportModel>("transportProperties"))
     {
-        const transportModel& laminarT =
+        const auto& laminarT =
             lookupObject<transportModel>("transportProperties");
 
         return rho()*laminarT.nu();
     }
     else if (foundObject<dictionary>("transportProperties"))
     {
-        const dictionary& transportProperties =
+        const auto& transportProperties =
              lookupObject<dictionary>("transportProperties");
 
-        dimensionedScalar nu("nu", dimViscosity, transportProperties);
+        const dimensionedScalar nu("nu", dimViscosity, transportProperties);
 
         return rho()*nu;
     }
@@ -443,11 +255,11 @@ Foam::tmp<Foam::volScalarField> Foam::functionObjects::forces::rho() const
                 mesh_
             ),
             mesh_,
-            dimensionedScalar("rho", dimDensity, rhoRef_)
+            dimensionedScalar(dimDensity, rhoRef_)
         );
     }
 
-    return(lookupObject<volScalarField>(rhoName_));
+    return (lookupObject<volScalarField>(rhoName_));
 }
 
 
@@ -455,7 +267,7 @@ Foam::scalar Foam::functionObjects::forces::rho(const volScalarField& p) const
 {
     if (p.dimensions() == dimPressure)
     {
-        return 1.0;
+        return 1;
     }
 
     if (rhoName_ != "rhoInf")
@@ -469,225 +281,159 @@ Foam::scalar Foam::functionObjects::forces::rho(const volScalarField& p) const
 }
 
 
-void Foam::functionObjects::forces::applyBins
+void Foam::functionObjects::forces::addToPatchFields
 (
+    const label patchi,
     const vectorField& Md,
-    const vectorField& fN,
-    const vectorField& fT,
     const vectorField& fP,
-    const vectorField& d
+    const vectorField& fV
 )
 {
-    if (nBin_ == 1)
-    {
-        force_[0][0] += sum(fN);
-        force_[1][0] += sum(fT);
-        force_[2][0] += sum(fP);
-        moment_[0][0] += sum(Md^fN);
-        moment_[1][0] += sum(Md^fT);
-        moment_[2][0] += sum(Md^fP);
-    }
-    else
-    {
-        scalarField dd((d & binDir_) - binMin_);
+    sumPatchForcesP_ += sum(fP);
+    sumPatchForcesV_ += sum(fV);
+    force_.boundaryFieldRef()[patchi] += fP + fV;
 
-        forAll(dd, i)
-        {
-            label bini = min(max(floor(dd[i]/binDx_), 0), force_[0].size() - 1);
-
-            force_[0][bini] += fN[i];
-            force_[1][bini] += fT[i];
-            force_[2][bini] += fP[i];
-            moment_[0][bini] += Md[i]^fN[i];
-            moment_[1][bini] += Md[i]^fT[i];
-            moment_[2][bini] += Md[i]^fP[i];
-        }
-    }
+    const vectorField mP(Md^fP);
+    const vectorField mV(Md^fV);
+
+    sumPatchMomentsP_ += sum(mP);
+    sumPatchMomentsV_ += sum(mV);
+    moment_.boundaryFieldRef()[patchi] += mP + mV;
 }
 
 
-void Foam::functionObjects::forces::addToFields
+void Foam::functionObjects::forces::addToInternalField
 (
-    const label patchi,
+    const labelList& cellIDs,
     const vectorField& Md,
-    const vectorField& fN,
-    const vectorField& fT,
-    const vectorField& fP
+    const vectorField& f
 )
 {
-    if (!writeFields_)
+    forAll(cellIDs, i)
     {
-        return;
-    }
+        const label celli = cellIDs[i];
 
-    auto& force = lookupObjectRef<volVectorField>(scopedName("force"));
-    vectorField& pf = force.boundaryFieldRef()[patchi];
-    pf += fN + fT + fP;
+        sumInternalForces_ += f[i];
+        force_[celli] += f[i];
 
-    auto& moment = lookupObjectRef<volVectorField>(scopedName("moment"));
-    vectorField& pm = moment.boundaryFieldRef()[patchi];
-    pm = Md^pf;
+        const vector m(Md[i]^f[i]);
+        sumInternalMoments_ += m;
+        moment_[celli] = m;
+    }
 }
 
 
-void Foam::functionObjects::forces::addToFields
-(
-    const labelList& cellIDs,
-    const vectorField& Md,
-    const vectorField& fN,
-    const vectorField& fT,
-    const vectorField& fP
-)
+void Foam::functionObjects::forces::createIntegratedDataFiles()
 {
-    if (!writeFields_)
+    if (!forceFilePtr_.valid())
     {
-        return;
+        forceFilePtr_ = createFile("force");
+        writeIntegratedDataFileHeader("Force", forceFilePtr_());
     }
 
-    auto& force = lookupObjectRef<volVectorField>(scopedName("force"));
-    auto& moment = lookupObjectRef<volVectorField>(scopedName("moment"));
-
-    forAll(cellIDs, i)
+    if (!momentFilePtr_.valid())
     {
-        label celli = cellIDs[i];
-        force[celli] += fN[i] + fT[i] + fP[i];
-        moment[celli] = Md[i]^force[celli];
+        momentFilePtr_ = createFile("moment");
+        writeIntegratedDataFileHeader("Moment", momentFilePtr_());
     }
 }
 
 
-void Foam::functionObjects::forces::writeIntegratedForceMoment
+void Foam::functionObjects::forces::writeIntegratedDataFileHeader
 (
-    const string& descriptor,
-    const vectorField& fm0,
-    const vectorField& fm1,
-    const vectorField& fm2,
-    autoPtr<OFstream>& osPtr
+    const word& header,
+    OFstream& os
 ) const
 {
-    vector pressure = sum(fm0);
-    vector viscous = sum(fm1);
-    vector porous = sum(fm2);
-    vector total = pressure + viscous + porous;
-
-    Log << "    Sum of " << descriptor.c_str() << nl
-        << "        Total    : " << total << nl
-        << "        Pressure : " << pressure << nl
-        << "        Viscous  : " << viscous << nl;
+    const auto& coordSys = coordSysPtr_();
+    const auto vecDesc = [](const word& root)->string
+    {
+        return root + "_x " + root + "_y " + root + "_z";
+    };
+    writeHeader(os, header);
+    writeHeaderValue(os, "CofR", coordSys.origin());
+    writeHeader(os, "");
+    writeCommented(os, "Time");
+    writeTabbed(os, vecDesc("total"));
+    writeTabbed(os, vecDesc("pressure"));
+    writeTabbed(os, vecDesc("viscous"));
 
     if (porosity_)
     {
-        Log << "        Porous   : " << porous << nl;
+        writeTabbed(os, vecDesc("porous"));
     }
 
-    if (writeToFile())
-    {
-        Ostream& os = osPtr();
-
-        writeCurrentTime(os);
-
-        os  << tab << total
-            << tab << pressure
-            << tab << viscous;
-
-        if (porosity_)
-        {
-            os  << tab << porous;
-        }
-
-        os  << endl;
-    }
+    os  << endl;
 }
 
 
-void Foam::functionObjects::forces::writeForces()
+void Foam::functionObjects::forces::writeIntegratedDataFiles()
 {
-    Log << type() << " " << name() << " write:" << nl;
-
     const auto& coordSys = coordSysPtr_();
 
-    writeIntegratedForceMoment
+    writeIntegratedDataFile
     (
-        "forces",
-        coordSys.localVector(force_[0]),
-        coordSys.localVector(force_[1]),
-        coordSys.localVector(force_[2]),
-        forceFilePtr_
+        coordSys.localVector(sumPatchForcesP_),
+        coordSys.localVector(sumPatchForcesV_),
+        coordSys.localVector(sumInternalForces_),
+        forceFilePtr_()
     );
 
-    writeIntegratedForceMoment
+    writeIntegratedDataFile
     (
-        "moments",
-        coordSys.localVector(moment_[0]),
-        coordSys.localVector(moment_[1]),
-        coordSys.localVector(moment_[2]),
-        momentFilePtr_
+        coordSys.localVector(sumPatchMomentsP_),
+        coordSys.localVector(sumPatchMomentsV_),
+        coordSys.localVector(sumInternalMoments_),
+        momentFilePtr_()
     );
-
-    Log << endl;
 }
 
 
-void Foam::functionObjects::forces::writeBinnedForceMoment
+void Foam::functionObjects::forces::writeIntegratedDataFile
 (
-    const List<Field<vector>>& fm,
-    autoPtr<OFstream>& osPtr
+    const vector& pres,
+    const vector& vis,
+    const vector& internal,
+    OFstream& os
 ) const
 {
-    if ((nBin_ == 1) || !writeToFile())
-    {
-        return;
-    }
+    writeCurrentTime(os);
 
-    List<Field<vector>> f(fm);
+    writeValue(os, pres + vis + internal);
+    writeValue(os, pres);
+    writeValue(os, vis);
 
-    if (binCumulative_)
+    if (porosity_)
     {
-        for (label i = 1; i < f[0].size(); i++)
-        {
-            f[0][i] += f[0][i-1];
-            f[1][i] += f[1][i-1];
-            f[2][i] += f[2][i-1];
-        }
+        writeValue(os, internal);
     }
 
-    Ostream& os = osPtr();
+    os  << endl;
+}
 
-    writeCurrentTime(os);
 
-    forAll(f[0], i)
+void Foam::functionObjects::forces::logIntegratedData
+(
+    const string& descriptor,
+    const vector& pres,
+    const vector& vis,
+    const vector& internal
+) const
+{
+    if (!log)
     {
-        vector total = f[0][i] + f[1][i] + f[2][i];
-
-        os  << tab << total
-            << tab << f[0][i]
-            << tab << f[1][i];
-
-        if (porosity_)
-        {
-            os  << tab << f[2][i];
-        }
+        return;
     }
 
-    os  << nl;
-}
-
-
-void Foam::functionObjects::forces::writeBins()
-{
-    const auto& coordSys = coordSysPtr_();
+    Log << "    Sum of " << descriptor.c_str() << nl
+        << "        Total    : " << (pres + vis + internal) << nl
+        << "        Pressure : " << pres << nl
+        << "        Viscous  : " << vis << nl;
 
-    List<Field<vector>> lf(3);
-    List<Field<vector>> lm(3);
-    lf[0] = coordSys.localVector(force_[0]);
-    lf[1] = coordSys.localVector(force_[1]);
-    lf[2] = coordSys.localVector(force_[2]);
-    lm[0] = coordSys.localVector(moment_[0]);
-    lm[1] = coordSys.localVector(moment_[1]);
-    lm[2] = coordSys.localVector(moment_[2]);
-
-    writeBinnedForceMoment(lf, forceBinFilePtr_);
-    writeBinnedForceMoment(lm, momentBinFilePtr_);
+    if (porosity_)
+    {
+        Log << "        Porous   : " << internal << nl;
+    }
 }
 
 
@@ -703,29 +449,50 @@ Foam::functionObjects::forces::forces
 :
     fvMeshFunctionObject(name, runTime, dict),
     writeFile(mesh_, name),
-    force_(3),
-    moment_(3),
+    force_
+    (
+        IOobject
+        (
+            "force", // scopedName() is not available at ctor level
+            time_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedVector(dimForce, Zero)
+    ),
+    moment_
+    (
+        IOobject
+        (
+            "moment",
+            time_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedVector(dimForce*dimLength, Zero)
+    ),
+    sumPatchForcesP_(Zero),
+    sumPatchForcesV_(Zero),
+    sumPatchMomentsP_(Zero),
+    sumPatchMomentsV_(Zero),
+    sumInternalForces_(Zero),
+    sumInternalMoments_(Zero),
     forceFilePtr_(),
     momentFilePtr_(),
-    forceBinFilePtr_(),
-    momentBinFilePtr_(),
+    coordSysPtr_(nullptr),
     patchSet_(),
+    rhoRef_(VGREAT),
+    pRef_(0),
     pName_("p"),
     UName_("U"),
     rhoName_("rho"),
-    directForceDensity_(false),
     fDName_("fD"),
-    rhoRef_(VGREAT),
-    pRef_(0),
-    coordSysPtr_(nullptr),
+    directForceDensity_(false),
     porosity_(false),
-    nBin_(1),
-    binDir_(Zero),
-    binDx_(0),
-    binMin_(GREAT),
-    binMax_(GREAT),
-    binPoints_(),
-    binCumulative_(true),
     writeFields_(false),
     initialised_(false)
 {
@@ -748,29 +515,50 @@ Foam::functionObjects::forces::forces
 :
     fvMeshFunctionObject(name, obr, dict),
     writeFile(mesh_, name),
-    force_(3),
-    moment_(3),
+    force_
+    (
+        IOobject
+        (
+            "force",
+            time_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedVector(dimForce, Zero)
+    ),
+    moment_
+    (
+        IOobject
+        (
+            "moment",
+            time_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedVector(dimForce*dimLength, Zero)
+    ),
+    sumPatchForcesP_(Zero),
+    sumPatchForcesV_(Zero),
+    sumPatchMomentsP_(Zero),
+    sumPatchMomentsV_(Zero),
+    sumInternalForces_(Zero),
+    sumInternalMoments_(Zero),
     forceFilePtr_(),
     momentFilePtr_(),
-    forceBinFilePtr_(),
-    momentBinFilePtr_(),
+    coordSysPtr_(nullptr),
     patchSet_(),
+    rhoRef_(VGREAT),
+    pRef_(0),
     pName_("p"),
     UName_("U"),
     rhoName_("rho"),
-    directForceDensity_(false),
     fDName_("fD"),
-    rhoRef_(VGREAT),
-    pRef_(0),
-    coordSysPtr_(nullptr),
+    directForceDensity_(false),
     porosity_(false),
-    nBin_(1),
-    binDir_(Zero),
-    binDx_(0),
-    binMin_(GREAT),
-    binMax_(GREAT),
-    binPoints_(),
-    binCumulative_(true),
     writeFields_(false),
     initialised_(false)
 {
@@ -780,17 +568,6 @@ Foam::functionObjects::forces::forces
         setCoordinateSystem(dict);
         Log << endl;
     }
-
-/*
-    // Turn off writing to file
-    writeToFile_ = false;
-
-    forAll(force_, i)
-    {
-        force_[i].setSize(nBin_, vector::zero);
-        moment_[i].setSize(nBin_, vector::zero);
-    }
-*/
 }
 
 
@@ -798,14 +575,14 @@ Foam::functionObjects::forces::forces
 
 bool Foam::functionObjects::forces::read(const dictionary& dict)
 {
-    fvMeshFunctionObject::read(dict);
-    writeFile::read(dict);
+    if (!fvMeshFunctionObject::read(dict) || !writeFile::read(dict))
+    {
+        return false;
+    }
 
     initialised_ = false;
 
-    Info<< type() << " " << name() << ":" << nl;
-
-    directForceDensity_ = dict.getOrDefault("directForceDensity", false);
+    Info<< type() << " " << name() << ":" << endl;
 
     patchSet_ =
         mesh_.boundaryMesh().patchSet
@@ -813,9 +590,10 @@ bool Foam::functionObjects::forces::read(const dictionary& dict)
             dict.get<wordRes>("patches")
         );
 
+    dict.readIfPresent("directForceDensity", directForceDensity_);
     if (directForceDensity_)
     {
-        // Optional entry for fDName
+        // Optional name entry for fD
         if (dict.readIfPresent<word>("fD", fDName_))
         {
             Info<< "    fD: " << fDName_ << endl;
@@ -840,7 +618,7 @@ bool Foam::functionObjects::forces::read(const dictionary& dict)
         // Reference density needed for incompressible calculations
         if (rhoName_ == "rhoInf")
         {
-            rhoRef_ = dict.get<scalar>("rhoInf");
+            rhoRef_ = dict.getCheck<scalar>("rhoInf", scalarMinMax::ge(SMALL));
             Info<< "    Freestream density (rhoInf) set to " << rhoRef_ << endl;
         }
 
@@ -861,115 +639,43 @@ bool Foam::functionObjects::forces::read(const dictionary& dict)
         Info<< "    Not including porosity effects" << endl;
     }
 
-    if (dict.found("binData"))
-    {
-        Info<< "    Activated data bins" << endl;
-        const dictionary& binDict(dict.subDict("binData"));
-        nBin_ = binDict.get<label>("nBin");
-
-        if (nBin_ < 0)
-        {
-            FatalIOErrorInFunction(dict)
-                << "Number of bins (nBin) must be zero or greater"
-                << exit(FatalIOError);
-        }
-        else if (nBin_ == 0)
-        {
-            // Case of no bins equates to a single bin to collect all data
-            nBin_ = 1;
-        }
-        else
-        {
-            Info<< "    Employing " << nBin_ << " bins" << endl;
-            if (binDict.readIfPresent("min", binMin_))
-            {
-                Info<< "    - min         : " << binMin_ << endl;
-            }
-            if (binDict.readIfPresent("max", binMax_))
-            {
-                Info<< "    - max         : " << binMax_ << endl;
-            }
-
-            binCumulative_ = binDict.get<bool>("cumulative");
-            Info<< "    - cumuluative : " << binCumulative_ << endl;
-
-            binDir_ = binDict.get<vector>("direction");
-            binDir_.normalise();
-            Info<< "    - direction   : " << binDir_ << endl;
-        }
-    }
-
     writeFields_ = dict.getOrDefault("writeFields", false);
-
     if (writeFields_)
     {
         Info<< "    Fields will be written" << endl;
-
-        volVectorField* forcePtr
-        (
-            new volVectorField
-            (
-                IOobject
-                (
-                    scopedName("force"),
-                    time_.timeName(),
-                    mesh_,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                mesh_,
-                dimensionedVector(dimForce, Zero)
-            )
-        );
-
-        mesh_.objectRegistry::store(forcePtr);
-
-        volVectorField* momentPtr
-        (
-            new volVectorField
-            (
-                IOobject
-                (
-                    scopedName("moment"),
-                    time_.timeName(),
-                    mesh_,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                mesh_,
-                dimensionedVector(dimForce*dimLength, Zero)
-            )
-        );
-
-        mesh_.objectRegistry::store(momentPtr);
     }
 
+    force_.rename(scopedName("force"));
+    moment_.rename(scopedName("moment"));
+
     return true;
 }
 
 
-void Foam::functionObjects::forces::calcForcesMoment()
+void Foam::functionObjects::forces::calcForcesMoments()
 {
     initialise();
 
-    resetFields();
+    reset();
 
     const point& origin = coordSysPtr_->origin();
 
     if (directForceDensity_)
     {
-        const volVectorField& fD = lookupObject<volVectorField>(fDName_);
+        const auto& fD = lookupObject<volVectorField>(fDName_);
 
-        const surfaceVectorField::Boundary& Sfb = mesh_.Sf().boundaryField();
+        const auto& Sfb = mesh_.Sf().boundaryField();
 
         for (const label patchi : patchSet_)
         {
-            vectorField Md(mesh_.C().boundaryField()[patchi] - origin);
+            const vectorField& d = mesh_.C().boundaryField()[patchi];
 
-            scalarField sA(mag(Sfb[patchi]));
+            const vectorField Md(d - origin);
 
-            // Normal force = surfaceUnitNormal*(surfaceNormal & forceDensity)
-            vectorField fN
+            const scalarField sA(mag(Sfb[patchi]));
+
+            // Pressure force = surfaceUnitNormal*(surfaceNormal & forceDensity)
+            const vectorField fP
             (
                 Sfb[patchi]/sA
                *(
@@ -977,72 +683,64 @@ void Foam::functionObjects::forces::calcForcesMoment()
                 )
             );
 
-            // Tangential force (total force minus normal fN)
-            vectorField fT(sA*fD.boundaryField()[patchi] - fN);
-
-            // Porous force
-            vectorField fP(Md.size(), Zero);
+            // Viscous force (total force minus pressure fP)
+            const vectorField fV(sA*fD.boundaryField()[patchi] - fP);
 
-            addToFields(patchi, Md, fN, fT, fP);
-
-            applyBins(Md, fN, fT, fP, mesh_.C().boundaryField()[patchi]);
+            addToPatchFields(patchi, Md, fP, fV);
         }
     }
     else
     {
-        const volScalarField& p = lookupObject<volScalarField>(pName_);
+        const auto& p = lookupObject<volScalarField>(pName_);
 
-        const surfaceVectorField::Boundary& Sfb = mesh_.Sf().boundaryField();
+        const auto& Sfb = mesh_.Sf().boundaryField();
 
         tmp<volSymmTensorField> tdevRhoReff = devRhoReff();
-        const volSymmTensorField::Boundary& devRhoReffb
-            = tdevRhoReff().boundaryField();
+        const auto& devRhoReffb = tdevRhoReff().boundaryField();
 
         // Scale pRef by density for incompressible simulations
-        scalar pRef = pRef_/rho(p);
+        const scalar rhoRef = rho(p);
+        const scalar pRef = pRef_/rhoRef;
 
         for (const label patchi : patchSet_)
         {
-            vectorField Md(mesh_.C().boundaryField()[patchi] - origin);
+            const vectorField& d = mesh_.C().boundaryField()[patchi];
+
+            const vectorField Md(d - origin);
 
-            vectorField fN
+            const vectorField fP
             (
-                rho(p)*Sfb[patchi]*(p.boundaryField()[patchi] - pRef)
+                rhoRef*Sfb[patchi]*(p.boundaryField()[patchi] - pRef)
             );
 
-            vectorField fT(Sfb[patchi] & devRhoReffb[patchi]);
-
-            vectorField fP(Md.size(), Zero);
-
-            addToFields(patchi, Md, fN, fT, fP);
+            const vectorField fV(Sfb[patchi] & devRhoReffb[patchi]);
 
-            applyBins(Md, fN, fT, fP, mesh_.C().boundaryField()[patchi]);
+            addToPatchFields(patchi, Md, fP, fV);
         }
     }
 
     if (porosity_)
     {
-        const volVectorField& U = lookupObject<volVectorField>(UName_);
+        const auto& U = lookupObject<volVectorField>(UName_);
         const volScalarField rho(this->rho());
         const volScalarField mu(this->mu());
 
-        const HashTable<const porosityModel*> models =
-            obr_.lookupClass<porosityModel>();
+        const auto models = obr_.lookupClass<porosityModel>();
 
         if (models.empty())
         {
             WarningInFunction
-                << "Porosity effects requested, but no porosity models found "
-                << "in the database"
+                << "Porosity effects requested, "
+                << "but no porosity models found in the database"
                 << endl;
         }
 
         forAllConstIters(models, iter)
         {
             // Non-const access required if mesh is changing
-            porosityModel& pm = const_cast<porosityModel&>(*iter());
+            auto& pm = const_cast<porosityModel&>(*iter());
 
-            vectorField fPTot(pm.force(U, rho, mu));
+            const vectorField fPTot(pm.force(U, rho, mu));
 
             const labelList& cellZoneIDs = pm.cellZoneIDs();
 
@@ -1054,55 +752,58 @@ void Foam::functionObjects::forces::calcForcesMoment()
                 const vectorField fP(fPTot, cZone);
                 const vectorField Md(d - origin);
 
-                const vectorField fDummy(Md.size(), Zero);
-
-                addToFields(cZone, Md, fDummy, fDummy, fP);
-
-                applyBins(Md, fDummy, fDummy, fP, d);
+                addToInternalField(cZone, Md, fP);
             }
         }
     }
 
-    Pstream::listCombineAllGather(force_, plusEqOp<vectorField>());
-    Pstream::listCombineAllGather(moment_, plusEqOp<vectorField>());
+    reduce(sumPatchForcesP_, sumOp<vector>());
+    reduce(sumPatchForcesV_, sumOp<vector>());
+    reduce(sumPatchMomentsP_, sumOp<vector>());
+    reduce(sumPatchMomentsV_, sumOp<vector>());
+    reduce(sumInternalForces_, sumOp<vector>());
+    reduce(sumInternalMoments_, sumOp<vector>());
 }
 
 
 Foam::vector Foam::functionObjects::forces::forceEff() const
 {
-    return sum(force_[0]) + sum(force_[1]) + sum(force_[2]);
+    return sumPatchForcesP_ + sumPatchForcesV_ + sumInternalForces_;
 }
 
 
 Foam::vector Foam::functionObjects::forces::momentEff() const
 {
-    return sum(moment_[0]) + sum(moment_[1]) + sum(moment_[2]);
+    return sumPatchMomentsP_ + sumPatchMomentsV_ + sumInternalMoments_;
 }
 
 
 bool Foam::functionObjects::forces::execute()
 {
-    calcForcesMoment();
+    calcForcesMoments();
 
-    if (Pstream::master())
-    {
-        createFiles();
+    Log << type() << " " << name() << " write:" << nl;
 
-        writeForces();
+    const auto& coordSys = coordSysPtr_();
 
-        writeBins();
+    const auto localFp(coordSys.localVector(sumPatchForcesP_));
+    const auto localFv(coordSys.localVector(sumPatchForcesV_));
+    const auto localFi(coordSys.localVector(sumInternalForces_));
 
-        Log << endl;
-    }
+    logIntegratedData("forces", localFp, localFv, localFi);
+
+    const auto localMp(coordSys.localVector(sumPatchMomentsP_));
+    const auto localMv(coordSys.localVector(sumPatchMomentsV_));
+    const auto localMi(coordSys.localVector(sumInternalMoments_));
 
-    // Write state/results information
-    setResult("normalForce", sum(force_[0]));
-    setResult("tangentialForce", sum(force_[1]));
-    setResult("porousForce", sum(force_[2]));
+    logIntegratedData("moments", localMp, localMv, localMi);
 
-    setResult("normalMoment", sum(moment_[0]));
-    setResult("tangentialMoment", sum(moment_[1]));
-    setResult("porousMoment", sum(moment_[2]));
+    setResult("pressureForce", localFp);
+    setResult("viscousForce", localFv);
+    setResult("internalForce", localFi);
+    setResult("pressureMoment", localMp);
+    setResult("viscousMoment", localMv);
+    setResult("internalMoment", localMi);
 
     return true;
 }
@@ -1110,12 +811,24 @@ bool Foam::functionObjects::forces::execute()
 
 bool Foam::functionObjects::forces::write()
 {
+    if (writeToFile())
+    {
+        Log << "    writing force and moment files." << endl;
+
+        createIntegratedDataFiles();
+        writeIntegratedDataFiles();
+    }
+
     if (writeFields_)
     {
-        lookupObject<volVectorField>(scopedName("force")).write();
-        lookupObject<volVectorField>(scopedName("moment")).write();
+        Log << "    writing force and moment fields." << endl;
+
+        force_.write();
+        moment_.write();
     }
 
+    Log << endl;
+
     return true;
 }
 
diff --git a/src/functionObjects/forces/forces/forces.H b/src/functionObjects/forces/forces/forces.H
index 23dc50d05fd2980b3d4a75d2db286f909d228122..da01eff4af34119614338c52b5ba40dd610f0d5f 100644
--- a/src/functionObjects/forces/forces/forces.H
+++ b/src/functionObjects/forces/forces/forces.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,127 +31,155 @@ Group
     grpForcesFunctionObjects
 
 Description
-    Calculates the forces and moments by integrating the pressure and
-    skin-friction forces over a given list of patches, and the resistance
-    from porous zones.
-
-    Forces and moments are calculated in a global Cartesian coordinate system
-    by default, or using a user-defined system.  Contributions can be 'binned'
-    according to a user-defined number of uniform-width collection zones (bins)
-    that span the input geometries, oriented by a user-defined direction vector.
-
-    Results are written to multiple files as a function of time in the
-    postProcessing/\<functionObjectName\> directory:
-    - force.dat          : forces
-    - moment.dat         : moments
-    - forceBin.dat       : force bins
-    - momentBin.dat      : moment bins
+    Computes forces and moments over a given list of patches by integrating
+    pressure and viscous forces and moments, and optionally resistance forces
+    and moments from porous zones.
 
-Usage
-    Example of function object specification:
-    \verbatim
-    forces1
-    {
-        type        forces;
-        libs        (forces);
-        ...
-        log         yes;
-        writeFields yes;
-        patches     (walls);
-
-        binData
-        {
-            nBin        20;
-            direction   (1 0 0);
-            cumulative  yes;
-        }
-    }
-    \endverbatim
+    Forces and moments are output in their total and constituent components:
+    - total forces and moments
+    - pressure contributions
+    - viscous contributions
+    - porous resistance contributions (optional)
 
-    Where the entries comprise:
-    \table
-        Property     | Description             | Required    | Default value
-        type         | Type name: forces       | yes         |
-        log          | Write force data to standard output | no | no
-        writeFields  | Write the force and moment fields | no | no
-        patches      | Patches included in the forces calculation | yes |
-        p            | Pressure field name     | no          | p
-        U            | Velocity field name     | no          | U
-        rho          | Density field name (see below) | no   | rho
-        CofR         | Centre of rotation (see below) | no   |
-        porosity     | flag to include porosity contributions | no | no
-        directForceDensity | Force density supplied directly (see below)|no|no
-        fD           | Name of force density field (see below) | no | fD
-    \endtable
+    Forces and moments can be computed and output in:
+    - the global Cartesian coordinate system (default)
+    - a user-defined Cartesian coordinate system
 
-    Bin data is optional, but if the dictionary is present, the entries must
-    be defined according o
+    Operands:
     \table
-        nBin         | number of data bins     | yes         |
-        direction    | direction along which bins are defined | yes |
-        cumulative   | bin data accumulated with incresing distance | yes |
+      Operand       | Type           | Location
+      input         | -              | -
+      output file   | dat            | postProcessing/\<FO\>/\<time\>/\<file\>s
+      output field  | volVectorField | \<time\>/\<outField\>s
     \endtable
 
-Note
-  - For incompressible cases, set \c rho to \c rhoInf.  You will then be
-    required to provide a \c rhoInf value corresponding to the free-stream
-    constant density.
-  - If the force density is supplied directly, set the \c directForceDensity
-    flag to 'yes', and supply the force density field using the \c
-    fDName entry
-  - The centre of rotation (CofR) for moment calculations can either be
-    specified by an \c CofR entry, or be taken from origin of the local
-    coordinate system.  For example,
+    where \c \<file\>s:
     \verbatim
-        CofR        (0 0 0);
+      force.dat        | Forces
+      moment.dat       | Moments
     \endverbatim
-    or
+
+    where \c \<outField\>s:
     \verbatim
-        origin  (0 0 0);
-        e1      (0 1 0);
-        e3      (0 0 1);
+      <namePrefix>:force   | Force field
+      <namePrefix>:moment  | Moment field
     \endverbatim
-    or
+
+Usage
+    Minimal example by using \c system/controlDict.functions:
     \verbatim
-        coordinateSystem
-        {
-            origin  (0 0 0);
-            rotation
+    <namePrefix>
+    {
+        // Mandatory entries
+        type                forces;
+        libs                (forces);
+        patches             (<wordRes>);
+
+        // Optional entries
+        directForceDensity  <bool>;
+        porosity            <bool>;
+        writeFields         <bool>;
+        useNamePrefix       <bool>;
+
+        // Conditional mandatory entries
+
+            // if directForceDensity == true
+            fD              <word>;
+
+
+            // Cartesian coordinate system specification when
+            // evaluating forces and moments, either of the below
+
+            // Option-1, i.e. the centre of rotation
+            // by inherently using e3=(0 0 1) and e1=(1 0 0)
+            CofR                (0 0 0); // Centre of rotation
+
+            // Option-2, i.e. local coordinate system specification
+            origin              (0 0 0);
+            e1                  (1 0 0);
+            e3                  (0 0 1); // (e1, e2) or (e2, e3) or (e3, e1)
+
+            // Option-3, i.e. general coordinate system specification
+            coordinateSystem
             {
-                type    axes;
-                e3      (0 0 1);
-                e1      (1 0 0);
+                type            cartesian;
+                origin          (0 0 0);
+                rotation
+                {
+                    type        axes;
+                    e3          (0 0 1);
+                    e1          (1 0 0); // (e1, e2) or (e2, e3) or (e3, e1)
+                }
             }
-        }
+
+        // Conditional optional entries
+
+            // if directForceDensity == false
+            p               <word>;
+            U               <word>;
+            rho             <word>;
+            rhoInf          <scalar>; // enabled if rho=rhoInf
+            pRef            <scalar>;
+
+        // Inherited entries
+        ...
+    }
     \endverbatim
 
-See also
-    Foam::functionObject
-    Foam::functionObjects::fvMeshFunctionObject
-    Foam::functionObjects::writeFile
-    Foam::functionObjects::timeControl
+    where the entries mean:
+    \table
+      Property   | Description               | Type | Reqd    | Deflt
+      type       | Type name: forces         | word | yes     | -
+      libs       | Library name: forces      | word | yes     | -
+      patches    | Names of operand patches  | wordRes | yes  | -
+      directForceDensity | Flag to directly supply force density <!--
+                 -->                         | bool | no      | false
+      porosity   | Flag to include porosity contributions | bool | no | false
+      writeFields | Flag to write force and moment fields | bool | no | false
+      useNamePrefix | Flag to include prefix for field names | bool | no | false
+      CofR    | Centre of rotation          | vector | cndtnl   | -
+      origin  | Origin of coordinate system | vector | cndtnl   | -
+      e3      | e3 coordinate axis          | vector | cndtnl   | -
+      e1      | e1 coordinate axis          | vector | cndtnl   | -
+      coordinateSystem | Coordinate system specifier | dictionary | cndtnl | -
+      fD      | Name of force density field | word   | cndtnl   | -
+      p       | Name of pressure field      | word   | cndtnl   | p
+      U       | Name of velocity field      | word   | cndtnl   | U
+      rho     | Name of density field       | word   | cndtnl   | rho
+      rhoInf  | Value of reference density  | scalar | cndtnl   | -
+      pRef    | Value of reference pressure | scalar | cndtnl   | 0
+    \endtable
+
+    The inherited entries are elaborated in:
+      - \link functionObject.H \endlink
+      - \link writeFile.H \endlink
+      - \link coordinateSystem.H \endlink
+
+Note
+  - For incompressible cases, set \c rho to \c rhoInf.
+    You will then be required to provide a \c rhoInf
+    value corresponding to the constant freestream density.
+  - \c writeControl and \c writeInterval entries of function
+    object do control when to output force and moment files and fields.
 
 SourceFiles
     forces.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef functionObjects_forces_H
-#define functionObjects_forces_H
+#ifndef Foam_functionObjects_forces_H
+#define Foam_functionObjects_forces_H
 
 #include "fvMeshFunctionObject.H"
 #include "writeFile.H"
 #include "coordinateSystem.H"
-#include "volFieldsFwd.H"
-#include "HashSet.H"
-#include "Tuple2.H"
-#include "OFstream.H"
-#include "Switch.H"
+#include "volFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
+
 namespace functionObjects
 {
 
@@ -164,107 +192,87 @@ class forces
     public fvMeshFunctionObject,
     public writeFile
 {
-
 protected:
 
-    // Protected data
-
-        //- Pressure, viscous and porous force per bin
-        List<Field<vector>> force_;
-
-        //- Pressure, viscous and porous moment per bin
-        List<Field<vector>> moment_;
+    // Protected Data
 
-        // File streams
+        // Fields
 
             //- Forces
-            autoPtr<OFstream> forceFilePtr_;
+            volVectorField force_;
 
             //- Moments
-            autoPtr<OFstream> momentFilePtr_;
+            volVectorField moment_;
 
-            //- Force bins
-            autoPtr<OFstream> forceBinFilePtr_;
+            //- Sum of patch pressure forces
+            vector sumPatchForcesP_;
 
-            //- Moment bins
-            autoPtr<OFstream> momentBinFilePtr_;
+            //- Sum of patch viscous forces
+            vector sumPatchForcesV_;
 
+            //- Sum of patch pressure moments
+            vector sumPatchMomentsP_;
 
-        // Read from dictionary
+            //- Sum of patch viscous moments
+            vector sumPatchMomentsV_;
 
-            //- Patches to integrate forces over
-            labelHashSet patchSet_;
+            //- Sum of internal forces
+            vector sumInternalForces_;
 
-            //- Name of pressure field
-            word pName_;
+            //- Sum of internal moments
+            vector sumInternalMoments_;
 
-            //- Name of velocity field
-            word UName_;
 
-            //- Name of density field (optional)
-            word rhoName_;
+        // File streams
 
-            //- Is the force density being supplied directly?
-            Switch directForceDensity_;
+            //- File stream for forces
+            autoPtr<OFstream> forceFilePtr_;
 
-            //- The name of the force density (fD) field
-            word fDName_;
+            //- File stream for moments
+            autoPtr<OFstream> momentFilePtr_;
 
-            //- Reference density needed for incompressible calculations
-            scalar rhoRef_;
 
-            //- Reference pressure
-            scalar pRef_;
+        // Read from dictionary
 
-            //- Coordinate system used when evaluating forces/moments
+            //- Coordinate system used when evaluating forces and moments
             autoPtr<coordinateSystem> coordSysPtr_;
 
-            //- Flag to include porosity effects
-            bool porosity_;
-
-
-            // Bin information
+            //- Names of operand patches
+            labelHashSet patchSet_;
 
-                //- Number of bins
-                label nBin_;
+            //- Reference density needed for incompressible calculations
+            scalar rhoRef_;
 
-                //- Direction used to determine bin orientation
-                vector binDir_;
+            //- Reference pressure
+            scalar pRef_;
 
-                //- Distance between bin divisions
-                scalar binDx_;
+            //- Name of pressure field
+            word pName_;
 
-                //- Minimum bin bounds
-                scalar binMin_;
+            //- Name of velocity field
+            word UName_;
 
-                //- Maximum bin bounds
-                scalar binMax_;
+            //- Name of density field
+            word rhoName_;
 
-                //- Bin positions along binDir
-                List<point> binPoints_;
+            //- Name of force density field
+            word fDName_;
 
-                //- Should bin data be cumulative?
-                bool binCumulative_;
+            //- Flag to directly supply force density
+            bool directForceDensity_;
 
+            //- Flag to include porosity effects
+            bool porosity_;
 
-            //- Write fields flag
+            //- Flag to write force and moment fields
             bool writeFields_;
 
-            //- Initialised flag
+            //- Flag of initialisation (internal)
             bool initialised_;
 
 
     // Protected Member Functions
 
-        //- Create the output files
-        void createFiles();
-
-        //- Write header for integrated data
-        void writeIntegratedHeader(const word& header, Ostream& os) const;
-
-        //- Write header for binned data
-        void writeBinHeader(const word& header, Ostream& os) const;
-
         //- Set the co-ordinate system from dictionary and axes names
         void setCoordinateSystem
         (
@@ -273,86 +281,79 @@ protected:
             const word& e1Name = word::null
         );
 
-        //- Initialise the fields
+        //- Initialise containers and fields
         void initialise();
 
-        //- Initialise the collection bins
-        void initialiseBins();
+        //- Reset containers and fields
+        void reset();
+
 
-        //- Reset the fields prior to accumulation of force/moments
-        void resetFields();
+    // Evaluation
 
-        //- Return the effective viscous stress (laminar + turbulent).
+        //- Return the effective stress (viscous + turbulent)
         tmp<volSymmTensorField> devRhoReff() const;
 
-        //- Dynamic viscosity field
+        //- Return dynamic viscosity field
         tmp<volScalarField> mu() const;
 
         //- Return rho if specified otherwise rhoRef
         tmp<volScalarField> rho() const;
 
-        //- Return rhoRef if the pressure field is dynamic, i.e. p/rho
-        //  otherwise return 1
+        //- Return rhoRef if the pressure field is
+        //- dynamic (i.e. p/rho), otherwise return 1
         scalar rho(const volScalarField& p) const;
 
-        //- Accumulate bin data
-        void applyBins
-        (
-            const vectorField& Md,
-            const vectorField& fN,
-            const vectorField& fT,
-            const vectorField& fP,
-            const vectorField& d
-        );
-
         //- Add patch contributions to force and moment fields
-        void addToFields
+        void addToPatchFields
         (
             const label patchi,
             const vectorField& Md,
-            const vectorField& fN,
-            const vectorField& fT,
-            const vectorField& fP
+            const vectorField& fP,
+            const vectorField& fV
         );
 
-        //- Add cell contributions to force and moment fields
-        void addToFields
+        //- Add cell contributions to force and
+        //- moment fields, and include porosity effects
+        void addToInternalField
         (
             const labelList& cellIDs,
             const vectorField& Md,
-            const vectorField& fN,
-            const vectorField& fT,
-            const vectorField& fP
+            const vectorField& f
         );
 
-        //- Helper function to write integrated forces and moments
-        void writeIntegratedForceMoment
-        (
-            const string& descriptor,
-            const vectorField& fm0,
-            const vectorField& fm1,
-            const vectorField& fm2,
-            autoPtr<OFstream>& osPtr
-        ) const;
 
-        //- Write force data
-        void writeForces();
+    // I-O
+
+        //- Create the integrated-data files
+        void createIntegratedDataFiles();
 
-        //- Helper function to write binned forces and moments
-        void writeBinnedForceMoment
+        //- Write header for an integrated-data file
+        void writeIntegratedDataFileHeader
         (
-            const List<Field<vector>>& fm,
-            autoPtr<OFstream>& osPtr
+            const word& header,
+            OFstream& os
         ) const;
 
-        //- Write binned data
-        void writeBins();
+        //- Write integrated data to files
+        void writeIntegratedDataFiles();
 
-        //- No copy construct
-        forces(const forces&) = delete;
+        //- Write integrated data to a file
+        void writeIntegratedDataFile
+        (
+            const vector& pres,
+            const vector& vis,
+            const vector& internal,
+            OFstream& os
+        ) const;
 
-        //- No copy assignment
-        void operator=(const forces&) = delete;
+        //- Write integrated data to stream
+        void logIntegratedData
+        (
+            const string& descriptor,
+            const vector& pres,
+            const vector& vis,
+            const vector& internal
+        ) const;
 
 
 public:
@@ -381,6 +382,12 @@ public:
             const bool readFields = true
         );
 
+        //- No copy construct
+        forces(const forces&) = delete;
+
+        //- No copy assignment
+        void operator=(const forces&) = delete;
+
 
     //- Destructor
     virtual ~forces() = default;
@@ -388,11 +395,11 @@ public:
 
     // Member Functions
 
-        //- Read the forces data
-        virtual bool read(const dictionary&);
+        //- Read the dictionary
+        virtual bool read(const dictionary& dict);
 
-        //- Calculate the forces and moments
-        virtual void calcForcesMoment();
+        //- Calculate forces and moments
+        virtual void calcForcesMoments();
 
         //- Return the total force
         virtual vector forceEff() const;
@@ -400,10 +407,10 @@ public:
         //- Return the total moment
         virtual vector momentEff() const;
 
-        //- Execute, currently does nothing
+        //- Execute the function object
         virtual bool execute();
 
-        //- Write the forces
+        //- Write to data files/fields and to streams
         virtual bool write();
 };
 
diff --git a/src/functionObjects/forces/propellerInfo/propellerInfo.C b/src/functionObjects/forces/propellerInfo/propellerInfo.C
index e74a2472dbcb7927e0a0e04297e2f90d6e05ee19..a8f5e88b6c80cace426cb6c827dd2cecfb5ced31 100644
--- a/src/functionObjects/forces/propellerInfo/propellerInfo.C
+++ b/src/functionObjects/forces/propellerInfo/propellerInfo.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -560,8 +560,8 @@ void Foam::functionObjects::propellerInfo::writePropellerPerformance()
     // Update n_
     setRotationalSpeed();
 
-    const vector sumForce(sum(force_[0]) + sum(force_[1]) + sum(force_[2]));
-    const vector sumMoment(sum(moment_[0]) + sum(moment_[1]) + sum(moment_[2]));
+    const vector sumForce = forceEff();
+    const vector sumMoment = momentEff();
 
     const scalar diameter = 2*radius_;
     const scalar URef = URefPtr_->value(time_.timeOutputValue());
@@ -866,7 +866,7 @@ bool Foam::functionObjects::propellerInfo::read(const dictionary& dict)
 
 bool Foam::functionObjects::propellerInfo::execute()
 {
-    calcForcesMoment();
+    calcForcesMoments();
 
     createFiles();
 
diff --git a/src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.C b/src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.C
index 71767b5a993d683aab10830d5dc2591e4c76b24f..359cd99a06d9d2c580905e0962b8219fe4d8e4c1 100644
--- a/src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.C
+++ b/src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.C
@@ -304,7 +304,7 @@ void Foam::rigidBodyMeshMotion::solve()
                 forcesDict.add("CofR", vector::zero);
 
                 functionObjects::forces f("forces", db(), forcesDict);
-                f.calcForcesMoment();
+                f.calcForcesMoments();
 
                 fx[bodyID] = ramp*spatialVector(f.momentEff(), f.forceEff());
             }
diff --git a/src/rigidBodyMeshMotion/rigidBodyMeshMotionSolver/rigidBodyMeshMotionSolver.C b/src/rigidBodyMeshMotion/rigidBodyMeshMotionSolver/rigidBodyMeshMotionSolver.C
index aebf1f763ad37c42a3028aefad7234f0f293b113..3517495f2c8865e8e71baaa4239db46b0e754f7b 100644
--- a/src/rigidBodyMeshMotion/rigidBodyMeshMotionSolver/rigidBodyMeshMotionSolver.C
+++ b/src/rigidBodyMeshMotion/rigidBodyMeshMotionSolver/rigidBodyMeshMotionSolver.C
@@ -228,7 +228,7 @@ void Foam::rigidBodyMeshMotionSolver::solve()
             forcesDict.add("CofR", vector::zero);
 
             functionObjects::forces f("forces", db(), forcesDict);
-            f.calcForcesMoment();
+            f.calcForcesMoments();
 
             fx[bodyID] = spatialVector(f.momentEff(), f.forceEff());
         }
diff --git a/src/sixDoFRigidBodyMotion/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C b/src/sixDoFRigidBodyMotion/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
index f56b85a746966cf30940530d5c6cc2a750e89778..3c6da0ce0cd76c2a4c262e79bd6ccfe523c9d852 100644
--- a/src/sixDoFRigidBodyMotion/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
+++ b/src/sixDoFRigidBodyMotion/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
@@ -221,7 +221,7 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
 
     functionObjects::forces f("forces", db(), forcesDict);
 
-    f.calcForcesMoment();
+    f.calcForcesMoments();
 
     // Get the forces on the patch faces at the current positions
 
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C
index af05c9b282efd3d34e441290efe83e91ce8edde6..eb6d00851c5039a19f4392a0964f889f14ffb506 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C
@@ -247,7 +247,7 @@ void Foam::sixDoFRigidBodyMotionSolver::solve()
 
         functionObjects::forces f("forces", db(), forcesDict);
 
-        f.calcForcesMoment();
+        f.calcForcesMoments();
 
         motion_.update
         (
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/forceCoeffs b/tutorials/incompressible/simpleFoam/motorBike/system/forceCoeffs
index b500bf64822166a5977764e18a22c6679f1c291a..b09c484fab4748014d411071cf3d6df2441fd77e 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/system/forceCoeffs
+++ b/tutorials/incompressible/simpleFoam/motorBike/system/forceCoeffs
@@ -27,15 +27,28 @@ forceCoeffs1
     magUInf         20;
     lRef            1.42;        // Wheelbase length
     Aref            0.75;        // Estimated
-    /*
+}
+
+/*
+binField1
+{
+    type                    binField;
+    libs                    (fieldFunctionObjects);
+    binModel                singleDirectionUniformBin;
+    fields                  (forceCoeff);
+    patches                 (motorBikeGroup);
+    decomposePatchValues    true;
+    CofR                    ${..forceCoeffs1.CofR};
+
     binData
     {
         nBin        20;          // output data into 20 bins
         direction   (1 0 0);     // bin direction
         cumulative  yes;
     }
-    */
+    writeControl            timeStep;
 }
+*/
 
 
 // ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/simpleCar/system/forceCoeffs b/tutorials/incompressible/simpleFoam/simpleCar/system/forceCoeffs
index 14185879cc0d9903d1bd8fd8de744ca954438b7d..56e9a063895400bf9d3987b435ced11aa660cfa0 100644
--- a/tutorials/incompressible/simpleFoam/simpleCar/system/forceCoeffs
+++ b/tutorials/incompressible/simpleFoam/simpleCar/system/forceCoeffs
@@ -27,13 +27,27 @@ forceCoeffs1
     lRef            4;          // Wheelbase length
     Aref            1;          // Estimated
     porosity        on;
+}
+
+
+binField1
+{
+    type                    binField;
+    libs                    (fieldFunctionObjects);
+    binModel                singleDirectionUniformBin;
+    fields                  (forceCoeff);
+    patches                 (body);
+    decomposePatchValues    yes;
+    CofR                    ${..forceCoeffs1.CofR};
+    cellZones               (porousZone);
 
     binData
     {
-        nBin            20;          // output data into 20 bins
-        direction       (1 0 0);     // bin direction
-        cumulative      yes;
+        nBin        20;          // output data into 20 bins
+        direction   (1 0 0);     // bin direction
+        cumulative  yes;
     }
+    writeControl            writeTime;
 }