diff --git a/src/finiteVolume/functionObjects/volRegion/volRegion.C b/src/finiteVolume/functionObjects/volRegion/volRegion.C
index 1f4779afebd605ed4f2bac46f1d72722fb36814c..c39c9cf8c5f9ad45c2124ba2e63c85da98f441be 100644
--- a/src/finiteVolume/functionObjects/volRegion/volRegion.C
+++ b/src/finiteVolume/functionObjects/volRegion/volRegion.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,6 +25,7 @@ License
 
 #include "volRegion.H"
 #include "volMesh.H"
+#include "cellSet.H"
 #include "globalMeshData.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -44,8 +45,9 @@ const Foam::Enum
 >
 Foam::functionObjects::volRegion::regionTypeNames_
 ({
-    { regionTypes::vrtCellZone, "cellZone" },
     { regionTypes::vrtAll, "all" },
+    { regionTypes::vrtCellSet, "cellSet" },
+    { regionTypes::vrtCellZone, "cellZone" },
 });
 
 
@@ -73,7 +75,7 @@ Foam::functionObjects::volRegion::volRegion
     const dictionary& dict
 )
 :
-    mesh_(mesh),
+    volMesh_(mesh),
     regionType_
     (
         regionTypeNames_.lookupOrDefault
@@ -83,7 +85,7 @@ Foam::functionObjects::volRegion::volRegion
             regionTypes::vrtAll
         )
     ),
-    regionName_(polyMesh::defaultRegion),
+    regionName_(volMesh_.name()),
     regionID_(-1)
 {
     read(dict);
@@ -94,12 +96,6 @@ Foam::functionObjects::volRegion::volRegion
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::functionObjects::volRegion::~volRegion()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool Foam::functionObjects::volRegion::read
@@ -107,19 +103,41 @@ bool Foam::functionObjects::volRegion::read
     const dictionary& dict
 )
 {
+    regionID_ = -1;
+    cellIds_.clear();
+
     switch (regionType_)
     {
+        case vrtCellSet:
+        {
+            dict.readEntry("name", regionName_);
+
+            cellIds_ = cellSet(volMesh_, regionName_).sortedToc();
+
+            if (nCells() == 0)
+            {
+                FatalIOErrorInFunction(dict)
+                    << regionTypeNames_[regionType_]
+                    << "(" << regionName_ << "):" << nl
+                    << "    Region has no cells"
+                    << exit(FatalIOError);
+            }
+
+            break;
+        }
+
         case vrtCellZone:
         {
             dict.readEntry("name", regionName_);
 
-            regionID_ = mesh_.cellZones().findZoneID(regionName_);
+            regionID_ = volMesh_.cellZones().findZoneID(regionName_);
 
             if (regionID_ < 0)
             {
                 FatalIOErrorInFunction(dict)
                     << "Unknown cell zone name: " << regionName_
-                    << ". Valid cell zones are: " << mesh_.cellZones().names()
+                    << ". Valid cell zones    : "
+                    << flatOutput(volMesh_.cellZones().names())
                     << exit(FatalIOError);
             }
 
@@ -137,6 +155,7 @@ bool Foam::functionObjects::volRegion::read
 
         case vrtAll:
         {
+            regionName_= volMesh_.name();
             break;
         }
 
@@ -144,7 +163,7 @@ bool Foam::functionObjects::volRegion::read
         {
             FatalIOErrorInFunction(dict)
                 << "Unknown region type. Valid region types are:"
-                << regionTypeNames_.sortedToc()
+                << flatOutput(regionTypeNames_.names()) << nl
                 << exit(FatalIOError);
         }
     }
@@ -155,14 +174,21 @@ bool Foam::functionObjects::volRegion::read
 
 const Foam::labelList& Foam::functionObjects::volRegion::cellIDs() const
 {
-    if (regionType_ == vrtAll)
-    {
-        return labelList::null();
-    }
-    else
+    switch (regionType_)
     {
-        return mesh_.cellZones()[regionID_];
+        case vrtCellSet:
+            return cellIds_;
+            break;
+
+        case vrtCellZone:
+            return volMesh_.cellZones()[regionID_];
+            break;
+
+        default:
+            break;
     }
+
+    return labelList::null();
 }
 
 
@@ -170,7 +196,7 @@ Foam::label Foam::functionObjects::volRegion::nCells() const
 {
     if (regionType_ == vrtAll)
     {
-        return mesh_.globalData().nTotalCells();
+        return volMesh_.globalData().nTotalCells();
     }
     else
     {
@@ -183,11 +209,17 @@ Foam::scalar Foam::functionObjects::volRegion::V() const
 {
     if (regionType_ == vrtAll)
     {
-        return gSum(mesh_.V());
+        return gSum(volMesh_.V());
     }
     else
     {
-        return gSum(scalarField(mesh_.V(), cellIDs()));
+        scalar vol = 0;
+        for (const label celli : cellIDs())
+        {
+            vol += volMesh_.V()[celli];
+        }
+
+        return returnReduce(vol, sumOp<scalar>());
     }
 }
 
diff --git a/src/finiteVolume/functionObjects/volRegion/volRegion.H b/src/finiteVolume/functionObjects/volRegion/volRegion.H
index f95504684d9d5e922c35b7a7169f02d40b76eab2..ef33682b94c8d13114363c9d932b85b62cd0db63 100644
--- a/src/finiteVolume/functionObjects/volRegion/volRegion.H
+++ b/src/finiteVolume/functionObjects/volRegion/volRegion.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -54,9 +54,9 @@ Description
 
 Usage
     \table
-        Property     | Description                | Required     | Default value
-        regionType   | cellZone or all              | no | all
-        name         | Name of cellZone if required | no |
+        Property     | Description                          | Required | Default
+        regionType   | Selection type: all/cellSet/cellZone | no | all
+        name         | Name of cellSet/cellZone if required | no |
     \endtable
 
 See also
@@ -78,7 +78,7 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of classes
+// Forward declarations
 class fvMesh;
 
 namespace functionObjects
@@ -90,12 +90,16 @@ namespace functionObjects
 
 class volRegion
 {
-    // Private member data
+    // Private Member Data
 
-        const fvMesh& mesh_;
+        const fvMesh& volMesh_;
+
+        //- The cell ids, from cellSet
+        labelList cellIds_;
 
         // Cache integral properties of the region for writeFileHeader
         label nCells_;
+
         scalar V_;
 
 
@@ -106,8 +110,9 @@ public:
         //- Region type enumeration
         enum regionTypes
         {
-            vrtCellZone,    //!< cell zone
-            vrtAll          //!< all cells
+            vrtAll,             //!< All cells
+            vrtCellSet,         //!< A cellSet
+            vrtCellZone         //!< A cellZone
         };
 
         //- Region type names
@@ -116,15 +121,15 @@ public:
 
 protected:
 
-    // Protected data
+    // Protected Data
 
         //- Region type
         regionTypes regionType_;
 
-        //- Region name (patch, zone, etc.)
+        //- Region name (cellSet, cellZone, ...)
         word regionName_;
 
-        //- Region ID (patch ID, zone ID, etc.)
+        //- Region ID (zone ID, ...)
         label regionID_;
 
 
@@ -147,13 +152,13 @@ public:
 
 
     //- Destructor
-    virtual ~volRegion();
+    virtual ~volRegion() = default;
 
 
     // Public Member Functions
 
         //- Read from dictionary
-        bool read(const dictionary&);
+        bool read(const dictionary& dict);
 
         //- Return the region type
         inline const regionTypes& regionType() const;
@@ -161,7 +166,7 @@ public:
         //- Return the local list of cell IDs
         const labelList& cellIDs() const;
 
-        //- Return the number of cells in the region
+        //- Return the number of cells selected in the region
         label nCells() const;
 
         //- Return total volume of the region
diff --git a/src/functionObjects/field/Make/files b/src/functionObjects/field/Make/files
index 1afdfb2ea7d2daafd0526756deb8c98eddab3761..d5aa38ee3dddab5e461818f15f825f0874bd423d 100644
--- a/src/functionObjects/field/Make/files
+++ b/src/functionObjects/field/Make/files
@@ -59,6 +59,7 @@ flowType/flowType.C
 CourantNo/CourantNo.C
 PecletNo/PecletNo.C
 blendingFactor/blendingFactor.C
+momentum/momentum.C
 pressure/pressure.C
 MachNo/MachNo.C
 Curle/Curle.C
diff --git a/src/functionObjects/field/momentum/momentum.C b/src/functionObjects/field/momentum/momentum.C
new file mode 100644
index 0000000000000000000000000000000000000000..53438c91af637b765758f06379647021a0a49824
--- /dev/null
+++ b/src/functionObjects/field/momentum/momentum.C
@@ -0,0 +1,573 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 "momentum.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "cellSet.H"
+#include "cylindricalRotation.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(momentum, 0);
+    addToRunTimeSelectionTable(functionObject, momentum, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
+
+template<class GeoField>
+Foam::autoPtr<GeoField>
+Foam::functionObjects::momentum::newField
+(
+    const word& baseName,
+    const dimensionSet& dims,
+    bool registerObject
+) const
+{
+    return
+        autoPtr<GeoField>::New
+        (
+            IOobject
+            (
+                scopedName(baseName),
+                time_.timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                registerObject
+            ),
+            mesh_,
+            dimensioned<typename GeoField::value_type>(dims, Zero)
+        );
+}
+
+
+void Foam::functionObjects::momentum::calc()
+{
+    initialise();
+
+    // When field writing is not enabled we need our local storage
+    // for the momentum and angular velocity fields
+    autoPtr<volVectorField> tmomentum, tAngularMom, tAngularVel;
+
+
+    // The base fields required
+    const auto& U = lookupObject<volVectorField>(UName_);
+    const auto* rhoPtr = findObject<volScalarField>(rhoName_);
+
+    const dimensionedScalar rhoRef("rho", dimDensity, rhoRef_);
+
+    // For quantities such as the mass-averaged angular velocity,
+    // we would calculate the mass per-cell here.
+
+    // tmp<volScalarField::Internal> tmass =
+    // (
+    //     rhoPtr
+    //   ? (mesh_.V() * (*rhoPtr))
+    //   : (mesh_.V() * rhoRef)
+    // );
+
+
+    // Linear momentum
+    // ~~~~~~~~~~~~~~~
+
+    auto* pmomentum = getObjectPtr<volVectorField>(scopedName("momentum"));
+
+    if (!pmomentum)
+    {
+        tmomentum = newField<volVectorField>("momentum", dimVelocity*dimMass);
+        pmomentum = tmomentum.get();  // get(), not release()
+    }
+    auto& momentum = *pmomentum;
+
+    if (rhoPtr)
+    {
+        momentum.ref() = (U * mesh_.V() * (*rhoPtr));
+    }
+    else
+    {
+        momentum.ref() = (U * mesh_.V() * rhoRef);
+    }
+    momentum.correctBoundaryConditions();
+
+
+    // Angular momentum
+    // ~~~~~~~~~~~~~~~~
+
+    auto* pAngularMom =
+        getObjectPtr<volVectorField>(scopedName("angularMomentum"));
+
+    if (hasCsys_ && !pAngularMom)
+    {
+        tAngularMom =
+            newField<volVectorField>("angularMomentum", dimVelocity*dimMass);
+        pAngularMom = tAngularMom.get();  // get(), not release()
+    }
+    else if (!pAngularMom)
+    {
+        // Angular momentum not requested, but alias to normal momentum
+        // to simplify logic when calculating the summations
+        pAngularMom = pmomentum;
+    }
+    auto& angularMom = *pAngularMom;
+
+
+    // Angular velocity
+    // ~~~~~~~~~~~~~~~~
+
+    auto* pAngularVel =
+        getObjectPtr<volVectorField>(scopedName("angularVelocity"));
+
+    if (hasCsys_)
+    {
+        if (!pAngularVel)
+        {
+            tAngularVel =
+                newField<volVectorField>("angularVelocity", dimVelocity);
+            pAngularVel = tAngularVel.get();  // get(), not release()
+        }
+        auto& angularVel = *pAngularVel;
+
+
+        // Global to local
+
+        angularVel.primitiveFieldRef() =
+            csys_.invTransform(mesh_.cellCentres(), U.internalField());
+
+        angularVel.correctBoundaryConditions();
+
+        if (rhoPtr)
+        {
+            angularMom.ref() = (angularVel * mesh_.V() * (*rhoPtr));
+        }
+        else
+        {
+            angularMom.ref() = (angularVel * mesh_.V() * rhoRef);
+        }
+
+        angularMom.correctBoundaryConditions();
+    }
+
+
+    // Integrate the selection
+
+    sumMomentum_ = Zero;
+    sumAngularMom_ = Zero;
+
+    switch (regionType_)
+    {
+        case vrtCellSet:
+        case vrtCellZone:
+        {
+            for (const label celli : cellIDs())
+            {
+                sumMomentum_ += momentum[celli];
+                sumAngularMom_ += angularMom[celli];
+            }
+            break;
+        }
+        case vrtAll:
+        {
+            for (label celli=0; celli < mesh_.nCells(); ++celli)
+            {
+                sumMomentum_ += momentum[celli];
+                sumAngularMom_ += angularMom[celli];
+            }
+            break;
+        }
+    }
+
+    reduce(sumMomentum_, sumOp<vector>());
+    reduce(sumAngularMom_, sumOp<vector>());
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void Foam::functionObjects::momentum::writeFileHeader(Ostream& os)
+{
+    if (!writeToFile() || writtenHeader_)
+    {
+        return;
+    }
+
+    if (hasCsys_)
+    {
+        writeHeader(os, "Momentum, Angular Momentum");
+        writeHeaderValue(os, "origin", csys_.origin());
+        writeHeaderValue(os, "axis", csys_.e3());
+    }
+    else
+    {
+        writeHeader(os, "Momentum");
+    }
+
+    if (regionType_ != vrtAll)
+    {
+        writeHeader
+        (
+            os,
+            "Selection " + regionTypeNames_[regionType_]
+          + " = " + regionName_
+        );
+    }
+
+    writeHeader(os, "");
+    writeCommented(os, "Time");
+    writeTabbed(os, "(momentum_x momentum_y momentum_z)");
+
+    if (hasCsys_)
+    {
+        writeTabbed(os, "(momentum_r momentum_rtheta momentum_axis)");
+    }
+    os  << endl;
+
+
+    writtenHeader_ = true;
+}
+
+
+void Foam::functionObjects::momentum::initialise()
+{
+    if (initialised_)
+    {
+        return;
+    }
+
+    if (!foundObject<volVectorField>(UName_))
+    {
+        FatalErrorInFunction
+            << "Could not find U: " << UName_ << " in database"
+            << exit(FatalError);
+    }
+
+
+    const auto* pPtr = cfindObject<volScalarField>(pName_);
+
+    if (pPtr && pPtr->dimensions() == dimPressure)
+    {
+        // Compressible - rho is mandatory
+
+        if (!foundObject<volScalarField>(rhoName_))
+        {
+            FatalErrorInFunction
+                << "Could not find rho:" << rhoName_
+                << exit(FatalError);
+        }
+    }
+
+    initialised_ = true;
+}
+
+
+void Foam::functionObjects::momentum::writeValues(Ostream& os)
+{
+    Log << type() << " " << name() << " write:" << nl;
+
+    Log << "    Sum of Momentum";
+
+    if (regionType_ != vrtAll)
+    {
+        Log << ' ' << regionTypeNames_[regionType_]
+            << ' ' << regionName_;
+    }
+
+    Log << nl
+        << "        linear  : " << sumMomentum_ << nl;
+
+    if (hasCsys_)
+    {
+        Log << "        angular : " << sumAngularMom_ << nl;
+    }
+
+    if (writeToFile())
+    {
+        writeTime(os);
+
+        os << tab << sumMomentum_;
+
+        if (hasCsys_)
+        {
+            os << tab << sumAngularMom_;
+        }
+        os << endl;
+    }
+
+    Log << endl;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::momentum::momentum
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict,
+    bool readFields
+)
+:
+    fvMeshFunctionObject(name, runTime, dict),
+    volRegion(fvMeshFunctionObject::mesh_, dict),
+    writeFile(mesh_, name, typeName, dict),
+    sumMomentum_(Zero),
+    sumAngularMom_(Zero),
+    UName_(),
+    pName_(),
+    rhoName_(),
+    rhoRef_(1.0),
+    csys_(),
+    hasCsys_(false),
+    writeMomentum_(false),
+    writeVelocity_(false),
+    writePosition_(false),
+    initialised_(false)
+{
+    if (readFields)
+    {
+        read(dict);
+        Log << endl;
+    }
+}
+
+
+Foam::functionObjects::momentum::momentum
+(
+    const word& name,
+    const objectRegistry& obr,
+    const dictionary& dict,
+    bool readFields
+)
+:
+    fvMeshFunctionObject(name, obr, dict),
+    volRegion(fvMeshFunctionObject::mesh_, dict),
+    writeFile(mesh_, name, typeName, dict),
+    sumMomentum_(Zero),
+    sumAngularMom_(Zero),
+    UName_(),
+    pName_(),
+    rhoName_(),
+    rhoRef_(1.0),
+    csys_(),
+    hasCsys_(false),
+    writeMomentum_(false),
+    writeVelocity_(false),
+    writePosition_(false),
+    initialised_(false)
+{
+    if (readFields)
+    {
+        read(dict);
+        Log << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::momentum::read(const dictionary& dict)
+{
+    fvMeshFunctionObject::read(dict);
+    volRegion::read(dict);
+    writeFile::read(dict);
+
+    initialised_ = false;
+
+    Info<< type() << " " << name() << ":" << nl;
+
+    // Optional entries U and p
+    UName_ = dict.lookupOrDefault<word>("U", "U");
+    pName_ = dict.lookupOrDefault<word>("p", "p");
+    rhoName_ = dict.lookupOrDefault<word>("rho", "rho");
+    rhoRef_ = dict.lookupOrDefault<scalar>("rhoRef", 1.0);
+
+    rhoRef_ = dict.lookupOrDefault<scalar>("rhoRef", 1.0);
+
+    hasCsys_ = dict.lookupOrDefault("cylindrical", false);
+
+    if (hasCsys_)
+    {
+        csys_ = coordSystem::cylindrical(dict);
+    }
+
+    writeMomentum_ = dict.lookupOrDefault("writeMomentum", false);
+    writeVelocity_ = dict.lookupOrDefault("writeVelocity", false);
+    writePosition_ = dict.lookupOrDefault("writePosition", false);
+
+    Info<<"Integrating for selection: "
+        << regionTypeNames_[regionType_]
+        << " (" << regionName_ << ")" << nl;
+
+    if (writeMomentum_)
+    {
+        Info<< "    Momentum fields will be written" << endl;
+
+        mesh_.objectRegistry::store
+        (
+            newField<volVectorField>("momentum", dimVelocity*dimMass)
+        );
+
+        if (hasCsys_)
+        {
+            mesh_.objectRegistry::store
+            (
+                newField<volVectorField>("angularMomentum", dimVelocity*dimMass)
+            );
+        }
+    }
+
+    if (hasCsys_)
+    {
+        if (writeVelocity_)
+        {
+            Info<< "    Angular velocity will be written" << endl;
+
+            mesh_.objectRegistry::store
+            (
+                newField<volVectorField>("angularVelocity", dimVelocity)
+            );
+        }
+
+        if (writePosition_)
+        {
+            Info<< "    Angular position will be written" << endl;
+        }
+    }
+
+    return true;
+}
+
+
+bool Foam::functionObjects::momentum::execute()
+{
+    calc();
+
+    if (Pstream::master())
+    {
+        writeFileHeader(file());
+
+        writeValues(file());
+
+        Log << endl;
+    }
+
+    // Write state/results information
+    setResult("momentum_x", sumMomentum_[0]);
+    setResult("momentum_y", sumMomentum_[1]);
+    setResult("momentum_z", sumMomentum_[2]);
+
+    setResult("momentum_r", sumAngularMom_[0]);
+    setResult("momentum_rtheta", sumAngularMom_[1]);
+    setResult("momentum_axis", sumAngularMom_[2]);
+
+    return true;
+}
+
+
+bool Foam::functionObjects::momentum::write()
+{
+    if (writeMomentum_ || (hasCsys_ && (writeVelocity_ || writePosition_)))
+    {
+        Log <<"Writing fields" << nl;
+
+        const volVectorField* fieldPtr;
+
+        fieldPtr = findObject<volVectorField>(scopedName("momentum"));
+        if (fieldPtr) fieldPtr->write();
+
+        fieldPtr = findObject<volVectorField>(scopedName("angularMomentum"));
+        if (fieldPtr) fieldPtr->write();
+
+        fieldPtr = findObject<volVectorField>(scopedName("angularVelocity"));
+        if (fieldPtr) fieldPtr->write();
+
+        if (hasCsys_ && writePosition_)
+        {
+            // Clunky, but currently no simple means of handling
+            // component-wise conversion and output
+
+            auto cyl_r = newField<volScalarField>("cyl_r", dimLength);
+            auto cyl_t = newField<volScalarField>("cyl_theta", dimless);
+            auto cyl_z = newField<volScalarField>("cyl_z", dimLength);
+
+            // Internal
+            {
+                const auto& pts = mesh_.cellCentres();
+                const label len = pts.size();
+
+                UList<scalar>& r = cyl_r->primitiveFieldRef(false);
+                UList<scalar>& t = cyl_t->primitiveFieldRef(false);
+                UList<scalar>& z = cyl_z->primitiveFieldRef(false);
+
+                for (label i=0; i < len; ++i)
+                {
+                    point p(csys_.localPosition(pts[i]));
+
+                    r[i] = p.x();
+                    t[i] = p.y();
+                    z[i] = p.z();
+                }
+            }
+
+            // Boundary
+            const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+
+            forAll(pbm, patchi)
+            {
+                const auto& pts = pbm[patchi].faceCentres();
+                const label len = pts.size();
+
+                UList<scalar>& r = cyl_r->boundaryFieldRef(false)[patchi];
+                UList<scalar>& t = cyl_t->boundaryFieldRef(false)[patchi];
+                UList<scalar>& z = cyl_z->boundaryFieldRef(false)[patchi];
+
+                for (label i=0; i < len; ++i)
+                {
+                    point p(csys_.localPosition(pts[i]));
+
+                    r[i] = p.x();
+                    t[i] = p.y();
+                    z[i] = p.z();
+                }
+            }
+
+            cyl_r->write();
+            cyl_t->write();
+            cyl_z->write();
+        }
+    }
+
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/momentum/momentum.H b/src/functionObjects/field/momentum/momentum.H
new file mode 100644
index 0000000000000000000000000000000000000000..2da62e1ea0b7e187a3f1957c5134fcd498c49965
--- /dev/null
+++ b/src/functionObjects/field/momentum/momentum.H
@@ -0,0 +1,266 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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::momentum
+
+Group
+    grpFieldFunctionObjects
+
+Description
+    Calculates linear/angular momentum, reporting integral values
+    and optionally writing the fields.
+
+    Data is written into momentum.dat in the
+    postProcessing/\<functionObjectName\> directory.
+
+Usage
+    Example of function object specification:
+    \verbatim
+    momentum1
+    {
+        type        momentum;
+        libs        ("libfieldFunctionObjects.so");
+        ...
+        log         yes;
+
+        regionType      all;
+        writeMomentum   yes;
+        writePosition   yes;
+        writeVelocity   yes;
+
+        cylindrical     true;
+
+        origin  (0 0 0);
+        e1      (1 0 0);
+        e3      (0 0 1);
+    }
+    \endverbatim
+
+    Where the entries comprise:
+    \table
+        Property     | Description                          | Required | Default
+        type         | Type name: momentum                  | yes |
+        log          | Log information to standard output   | no  | no
+        writeMomentum | Write (linear, angular) momentum  fields  | no | no
+        writePosition | Write angular position component fields   | no | no
+        writeVelocity | Write angular velocity fields       | no  | no
+        p            | Pressure field name                  | no  | p
+        U            | Velocity field name                  | no  | U
+        rho          | Density field name                   | no  | rho
+        rhoRef       | Reference density (incompressible)   | no  | 1.0
+        cylindrical  | Use cylindrical coordinates          | no  | false
+        origin       | Origin for cylindrical coordinates   | no  |
+        regionType   | Selection type: all/cellSet/cellZone | no  | all
+        name         | Name of cellSet/cellZone if required | no  |
+    \endtable
+
+Note
+  - For incompressible cases, the value of \c rhoRef is used.
+  - When specifying the cylindrical coordinate system, the rotation
+    can be specified directly with e1/e2/e3 axes, or via a \c rotation
+    sub-dictionary
+    For example,
+    \verbatim
+        origin      (0 0 0);
+        rotation
+        {
+            type    cylindrical;
+            axis    (0 0 1);
+        }
+    \endverbatim
+
+See also
+    Foam::functionObject
+    Foam::functionObjects::fvMeshFunctionObject
+    Foam::functionObjects::volRegion
+    Foam::functionObjects::writeFile
+    Foam::functionObjects::timeControl
+
+SourceFiles
+    momentum.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjects_momentum_H
+#define functionObjects_momentum_H
+
+#include "fvMeshFunctionObject.H"
+#include "writeFile.H"
+#include "cylindricalCS.H"
+#include "volFieldsFwd.H"
+#include "volRegion.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declarations
+class dimensionSet;
+
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class momentum Declaration
+\*---------------------------------------------------------------------------*/
+
+class momentum
+:
+    public fvMeshFunctionObject,
+    public volRegion,
+    public writeFile
+{
+    // Private Member Functions
+
+        //- Calculate the fields and integral values
+        void calc();
+
+        //- Allocate a new zero geometric field
+        template<class GeoField>
+        autoPtr<GeoField> newField
+        (
+            const word& baseName,
+            const dimensionSet& dims,
+            bool registerObject=true
+        ) const;
+
+
+protected:
+
+    // Protected data
+
+        //- Integral (linear) momentum
+        vector sumMomentum_;
+
+        //- Integral angular momentum
+        vector sumAngularMom_;
+
+
+    // Read from dictionary
+
+        //- The velocity field name (optional)
+        word UName_;
+
+        //- The pressure field name (optional)
+        //  Only used to determine incompressible/compressible
+        word pName_;
+
+        //- The density field name (optional)
+        word rhoName_;
+
+        //- Reference density (for incompressible)
+        scalar rhoRef_;
+
+        //- Coordinate system for evaluating angular momentum
+        coordSystem::cylindrical csys_;
+
+        //- Are we using the cylindrical coordinate system?
+        bool hasCsys_;
+
+        //- Write fields flag
+        bool writeMomentum_;
+
+        //- Write fields flag
+        bool writeVelocity_;
+
+        //- Write fields flag
+        bool writePosition_;
+
+        //- Initialised flag
+        bool initialised_;
+
+
+    // Protected Member Functions
+
+        //- Initialise the fields
+        void initialise();
+
+        //- Output file header information
+        virtual void writeFileHeader(Ostream& os);
+
+        //- Write momentum data
+        void writeValues(Ostream& os);
+
+        //- No copy construct
+        momentum(const momentum&) = delete;
+
+        //- No copy assignment
+        void operator=(const momentum&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("momentum");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        momentum
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict,
+            const bool readFields = true
+        );
+
+        //- Construct from objectRegistry and dictionary
+        momentum
+        (
+            const word& name,
+            const objectRegistry& obr,
+            const dictionary& dict,
+            const bool readFields = true
+        );
+
+
+    //- Destructor
+    virtual ~momentum() = default;
+
+
+    // Member Functions
+
+        //- Read the momentum data
+        virtual bool read(const dictionary&);
+
+        //- Calculate and report the integral momentum
+        virtual bool execute();
+
+        //- Write the momentum, possibly angular momentum and velocity
+        virtual bool write();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/pipeCyclic/system/controlDict b/tutorials/incompressible/simpleFoam/pipeCyclic/system/controlDict
index fae83278cb6cfaef6d6e55823dc0d0a9e238fc04..0beed5f5bd3ff006a9c90363eabc5d6a1b17abbd 100644
--- a/tutorials/incompressible/simpleFoam/pipeCyclic/system/controlDict
+++ b/tutorials/incompressible/simpleFoam/pipeCyclic/system/controlDict
@@ -48,6 +48,7 @@ runTimeModifiable true;
 functions
 {
     #include "coordinateTransform"
+    #include "momentum"
 }
 
 
diff --git a/tutorials/incompressible/simpleFoam/pipeCyclic/system/momentum b/tutorials/incompressible/simpleFoam/pipeCyclic/system/momentum
new file mode 100644
index 0000000000000000000000000000000000000000..6f05784e69fe58f92b9a0c0abcad98bbd88c7e88
--- /dev/null
+++ b/tutorials/incompressible/simpleFoam/pipeCyclic/system/momentum
@@ -0,0 +1,33 @@
+// -*- C++ -*-
+// Calculate momentum fields
+momentum
+{
+    type    momentum;
+    libs    ("libfieldFunctionObjects.so");
+    log     true;
+
+    executeInterval 10;
+    writeControl    writeTime;
+
+    // writeToFile     true;
+
+    writeMomentum   true;
+    writePosition   true;
+    writeVelocity   true;
+
+    // Cells to select (all/cellSet/cellZone)
+    regionType  all;
+    // name     c0;
+
+    cylindrical true;
+
+    origin  (0 0 0);
+    rotation
+    {
+        type cylindrical;
+        axis (1 0 0);   //< local Z
+    }
+}
+
+
+// ************************************************************************* //