ENH: forces/forceCoeffs FOs updated follwing changes to functionObjectFile

src/postProcessing/functionObjects/forces/forceCoeffs/forceCoeffs.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,6 +28,9 @@ License
 #include "Time.H"
 #include "Pstream.H"
 #include "IOmanip.H"
+#include "fvMesh.H"
+#include "dimensionedTypes.H"
+#include "volFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -39,81 +42,157 @@ namespace Foam
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-void Foam::forceCoeffs::writeFileHeader(const label i)
+void Foam::forceCoeffs::createFiles()
 {
-    if (i == 0)
+    // Note: Only possible to create bin files after bins have been initialised
-    {
-        // force coeff data
 
-        writeHeader(file(i), "Force coefficients");
+    if (writeToFile() && !coeffFilePtr_.valid())
-        writeHeaderValue(file(i), "liftDir", liftDir_);
+    {
-        writeHeaderValue(file(i), "dragDir", dragDir_);
+        coeffFilePtr_ = createFile("coefficient");
-        writeHeaderValue(file(i), "pitchAxis", pitchAxis_);
+        writeIntegratedHeader("Coefficients", coeffFilePtr_());
-        writeHeaderValue(file(i), "magUInf", magUInf_);
+
-        writeHeaderValue(file(i), "lRef", lRef_);
+        if (nBin_ > 1)
-        writeHeaderValue(file(i), "Aref", Aref_);
+        {
-        writeHeaderValue(file(i), "CofR", coordSys_.origin());
+            CmBinFilePtr_ = createFile("CmBin");
-        writeCommented(file(i), "Time");
+            writeBinHeader("Moment coefficient bins", CmBinFilePtr_());
-        writeTabbed(file(i), "Cm");
+            CdBinFilePtr_ = createFile("CdBin");
-        writeTabbed(file(i), "Cd");
+            writeBinHeader("Drag coefficient bins", CdBinFilePtr_());
-        writeTabbed(file(i), "Cl");
+            ClBinFilePtr_ = createFile("ClBin");
-        writeTabbed(file(i), "Cl(f)");
+            writeBinHeader("Lift coefficient bins", ClBinFilePtr_());
-        writeTabbed(file(i), "Cl(r)");
-        file(i)
-            << tab << "Cm" << tab << "Cd" << tab << "Cl" << tab << "Cl(f)"
-            << tab << "Cl(r)";
         }
-    else if (i == 1)
+    }
-    {
+}
-        // bin coeff data
 
-        writeHeader(file(i), "Force coefficient bins");
+
-        writeHeaderValue(file(i), "bins", nBin_);
+void Foam::forceCoeffs::writeIntegratedHeader
-        writeHeaderValue(file(i), "start", binMin_);
+(
-        writeHeaderValue(file(i), "delta", binDx_);
+    const word& header,
-        writeHeaderValue(file(i), "direction", binDir_);
+    Ostream& os
+) const
+{
+    writeHeader(os, "Force coefficients");
+    writeHeaderValue(os, "liftDir", liftDir_);
+    writeHeaderValue(os, "dragDir", dragDir_);
+    writeHeaderValue(os, "pitchAxis", pitchAxis_);
+    writeHeaderValue(os, "magUInf", magUInf_);
+    writeHeaderValue(os, "lRef", lRef_);
+    writeHeaderValue(os, "Aref", Aref_);
+    writeHeaderValue(os, "CofR", coordSys_.origin());
+    writeHeader(os, "");
+    writeCommented(os, "Time");
+    writeTabbed(os, "Cm");
+    writeTabbed(os, "Cd");
+    writeTabbed(os, "Cl");
+    writeTabbed(os, "Cl(f)");
+    writeTabbed(os, "Cl(r)");
+    os  << endl;
+}
+
+
+void Foam::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(file(i), "x co-ords  :");
+    writeCommented(os, "x co-ords  :");
     forAll(binPoints, pointI)
     {
         binPoints[pointI] = (binMin_ + (pointI + 1)*binDx_)*binDir_;
-            file(i) << tab << binPoints[pointI].x();
+        os << tab << binPoints[pointI].x();
     }
-        file(i) << nl;
+    os << nl;
 
-        writeCommented(file(i), "y co-ords  :");
+    writeCommented(os, "y co-ords  :");
     forAll(binPoints, pointI)
     {
-            file(i) << tab << binPoints[pointI].y();
+        os << tab << binPoints[pointI].y();
     }
-        file(i) << nl;
+    os << nl;
 
-        writeCommented(file(i), "z co-ords  :");
+    writeCommented(os, "z co-ords  :");
     forAll(binPoints, pointI)
     {
-            file(i) << tab << binPoints[pointI].z();
+        os << tab << binPoints[pointI].z();
     }
-        file(i) << nl;
+    os << nl;
 
-        writeCommented(file(i), "Time");
+    writeHeader(os, "");
+    writeCommented(os, "Time");
 
     for (label j = 0; j < nBin_; j++)
     {
-            const word jn('(' + Foam::name(j) + ')');
+        word jn(Foam::name(j) + ':');
-            writeTabbed(file(i), "Cm" + jn);
+        writeTabbed(os, jn + "total");
-            writeTabbed(file(i), "Cd" + jn);
+        writeTabbed(os, jn + "pressure");
-            writeTabbed(file(i), "Cl" + jn);
+        writeTabbed(os, jn + "viscous");
-        }
+
-    }
+        if (porosity_)
-    else
         {
-        FatalErrorIn("void Foam::forces::writeFileHeader(const label)")
+            writeTabbed(os, jn + "porous");
-            << "Unhandled file index: " << i
+        }
-            << abort(FatalError);
     }
 
-    file(i)<< endl;
+    os  << endl;
+}
+
+
+void Foam::forceCoeffs::writeIntegratedData
+(
+    const word& title,
+    const List<Field<scalar> >& coeff
+) const
+{
+    scalar pressure = sum(coeff[0]);
+    scalar viscous = sum(coeff[1]);
+    scalar porous = sum(coeff[2]);
+    scalar total = pressure + viscous + porous;
+
+    if (log_)
+    {
+        Info<< "        " << title << "       : " << total << token::TAB
+            << "("
+            << "pressure: " << pressure << token::TAB
+            << "viscous: " << viscous;
+
+        if (porosity_)
+        {
+            Info<< token::TAB << "porous: " << porous;
+        }
+
+        Info<< ")" << endl;
+    }
+}
+
+
+void Foam::forceCoeffs::writeBinData
+(
+    const List<Field<scalar> > coeffs,
+    Ostream& os
+) const
+{
+    os  << obr_.time().value();
+
+    for (label binI = 0; binI < nBin_; binI++)
+    {
+        scalar total = coeffs[0][binI] + coeffs[1][binI] + coeffs[2][binI];
+
+        os  << tab << total << tab << coeffs[0][binI] << tab << coeffs[1][binI];
+
+        if (porosity_)
+        {
+            os  << tab << coeffs[2][binI];
+        }
+    }
+
+    os  << endl;
 }
 
 
@@ -124,7 +203,8 @@ Foam::forceCoeffs::forceCoeffs
     const word& name,
     const objectRegistry& obr,
     const dictionary& dict,
-    const bool loadFromFiles
+    const bool loadFromFiles,
+    const bool readFields
 )
 :
     forces(name, obr, dict, loadFromFiles, false),
@@ -133,11 +213,17 @@ Foam::forceCoeffs::forceCoeffs
     pitchAxis_(vector::zero),
     magUInf_(0.0),
     lRef_(0.0),
-    Aref_(0.0)
+    Aref_(0.0),
+    coeffFilePtr_(),
+    CmBinFilePtr_(),
+    CdBinFilePtr_(),
+    ClBinFilePtr_()
 {
+    if (readFields)
+    {
         read(dict);
-
+        if (log_) Info << endl;
-    Info<< endl;
+    }
 }
 
 
@@ -151,8 +237,11 @@ Foam::forceCoeffs::~forceCoeffs()
 
 void Foam::forceCoeffs::read(const dictionary& dict)
 {
-    if (active_)
+    if (!active_)
     {
+        return;
+    }
+
     forces::read(dict);
 
     // Directions for lift and drag forces, and pitch moment
@@ -166,100 +255,198 @@ void Foam::forceCoeffs::read(const dictionary& dict)
     // Reference length and area scales
     dict.lookup("lRef") >> lRef_;
     dict.lookup("Aref") >> Aref_;
+
+    if (writeFields_)
+    {
+        const fvMesh& mesh = refCast<const fvMesh>(obr_);
+
+        tmp<volVectorField> tforceCoeff
+        (
+            new volVectorField
+            (
+                IOobject
+                (
+                    fieldName("forceCoeff"),
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                dimensionedVector("0", dimless, vector::zero)
+            )
+        );
+
+        obr_.store(tforceCoeff.ptr());
+
+        tmp<volVectorField> tmomentCoeff
+        (
+            new volVectorField
+            (
+                IOobject
+                (
+                    fieldName("momentCoeff"),
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                dimensionedVector("0", dimless, vector::zero)
+            )
+        );
+
+        obr_.store(tmomentCoeff.ptr());
     }
 }
 
 
 void Foam::forceCoeffs::execute()
 {
-    // Do nothing - only valid on write
-}
-
-
-void Foam::forceCoeffs::end()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::forceCoeffs::timeSet()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::forceCoeffs::write()
-{
-    forces::calcForcesMoment();
-
     if (!active_)
     {
         return;
     }
 
-    if (Pstream::master())
+    forces::calcForcesMoment();
-    {
+
-        functionObjectFile::write();
+    createFiles();
 
     scalar pDyn = 0.5*rhoRef_*magUInf_*magUInf_;
 
-        Field<vector> totForce(force_[0] + force_[1] + force_[2]);
+    // Storage for pressure, viscous and porous contributions to coeffs
-        Field<vector> totMoment(moment_[0] + moment_[1] + moment_[2]);
+    List<Field<scalar> > momentCoeffs(3);
+    List<Field<scalar> > dragCoeffs(3);
+    List<Field<scalar> > liftCoeffs(3);
+    forAll(liftCoeffs, i)
+    {
+        momentCoeffs[i].setSize(nBin_);
+        dragCoeffs[i].setSize(nBin_);
+        liftCoeffs[i].setSize(nBin_);
+    }
 
-        List<Field<scalar> > coeffs(3);
+    // Calculate coefficients
-        coeffs[0].setSize(nBin_);
+    scalar CmTot = 0;
-        coeffs[1].setSize(nBin_);
+    scalar CdTot = 0;
-        coeffs[2].setSize(nBin_);
+    scalar ClTot = 0;
+    forAll(liftCoeffs, i)
+    {
+        momentCoeffs[i] = (moment_[i] & pitchAxis_)/(Aref_*pDyn*lRef_);
+        dragCoeffs[i] = (force_[i] & dragDir_)/(Aref_*pDyn);
+        liftCoeffs[i] = (force_[i] & liftDir_)/(Aref_*pDyn);
 
-        // lift, drag and moment
+        CmTot += sum(momentCoeffs[i]);
-        coeffs[0] = (totForce & liftDir_)/(Aref_*pDyn);
+        CdTot += sum(dragCoeffs[i]);
-        coeffs[1] = (totForce & dragDir_)/(Aref_*pDyn);
+        ClTot += sum(liftCoeffs[i]);
-        coeffs[2] = (totMoment & pitchAxis_)/(Aref_*lRef_*pDyn);
+    }
 
-        scalar Cl = sum(coeffs[0]);
+    scalar ClfTot = ClTot/2.0 + CmTot;
-        scalar Cd = sum(coeffs[1]);
+    scalar ClrTot = ClTot/2.0 - CmTot;
-        scalar Cm = sum(coeffs[2]);
 
-        scalar Clf = Cl/2.0 + Cm;
+    if (log_)
-        scalar Clr = Cl/2.0 - Cm;
+    {
+        Info<< type() << " " << name_ << " output:" << nl
+            << "    Coefficients" << nl;
+    }
+    writeIntegratedData("Cm", momentCoeffs);
+    writeIntegratedData("Cd", dragCoeffs);
+    writeIntegratedData("Cl", liftCoeffs);
+    if (log_)
+    {
+        Info<< "        Cl(f)    : " << ClfTot << nl
+            << "        Cl(r)    : " << ClrTot << nl
+            << endl;
+    }
 
-        file(0)
+    if (writeToFile())
-            << obr_.time().value() << tab << Cm << tab  << Cd
+    {
-            << tab << Cl << tab << Clf << tab << Clr << endl;
+        coeffFilePtr_()
+            << obr_.time().value() << tab << CmTot << tab  << CdTot
+            << tab << ClTot << tab << ClfTot << tab << ClrTot << endl;
 
-        if (log_) Info<< type() << " " << name_ << " output:" << nl
-            << "    Cm    = " << Cm << nl
-            << "    Cd    = " << Cd << nl
-            << "    Cl    = " << Cl << nl
-            << "    Cl(f) = " << Clf << nl
-            << "    Cl(r) = " << Clr << endl;
 
         if (nBin_ > 1)
         {
             if (binCumulative_)
             {
-                for (label i = 1; i < coeffs[0].size(); i++)
+                forAll(liftCoeffs, i)
                 {
-                    coeffs[0][i] += coeffs[0][i-1];
+                    for (label binI = 1; binI < nBin_; binI++)
-                    coeffs[1][i] += coeffs[1][i-1];
-                    coeffs[2][i] += coeffs[2][i-1];
-                }
-            }
-
-            file(1)<< obr_.time().value();
-
-            forAll(coeffs[0], i)
                     {
-                file(1)
+                        liftCoeffs[i][binI] += liftCoeffs[i][binI-1];
-                    << tab << coeffs[2][i]
+                        dragCoeffs[i][binI] += dragCoeffs[i][binI-1];
-                    << tab << coeffs[1][i]
+                        momentCoeffs[i][binI] += momentCoeffs[i][binI-1];
-                    << tab << coeffs[0][i];
+                    }
+                }
             }
 
-            file(1) << endl;
+            writeBinData(dragCoeffs, CdBinFilePtr_());
+            writeBinData(liftCoeffs, ClBinFilePtr_());
+            writeBinData(momentCoeffs, CmBinFilePtr_());
+        }
     }
 
-        if (log_) Info<< endl;
+    if (writeFields_)
+    {
+        const volVectorField& force =
+            obr_.lookupObject<volVectorField>(fieldName("force"));
+
+        const volVectorField& moment =
+            obr_.lookupObject<volVectorField>(fieldName("moment"));
+
+        volVectorField& forceCoeff =
+            const_cast<volVectorField&>
+            (
+                obr_.lookupObject<volVectorField>(fieldName("forceCoeff"))
+            );
+
+        volVectorField& momentCoeff =
+            const_cast<volVectorField&>
+            (
+                obr_.lookupObject<volVectorField>(fieldName("momentCoeff"))
+            );
+
+        dimensionedScalar f0("f0", dimForce, Aref_*pDyn);
+        dimensionedScalar m0("m0", dimForce*dimLength, Aref_*lRef_*pDyn);
+
+        forceCoeff == force/f0;
+        momentCoeff == moment/m0;
+    }
+}
+
+
+void Foam::forceCoeffs::end()
+{
+    if (active_)
+    {
+        execute();
+    }
+}
+
+
+void Foam::forceCoeffs::timeSet()
+{
+    // Do nothing
+}
+
+
+void Foam::forceCoeffs::write()
+{
+    if (!active_)
+    {
+        return;
+    }
+
+    if (writeFields_)
+    {
+        const volVectorField& forceCoeff =
+            obr_.lookupObject<volVectorField>(fieldName("forceCoeff"));
+
+        const volVectorField& momentCoeff =
+            obr_.lookupObject<volVectorField>(fieldName("momentCoeff"));
+
+        forceCoeff.write();
+        momentCoeff.write();
     }
 }
 

src/postProcessing/functionObjects/forces/forceCoeffs/forceCoeffs.H

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,6 +32,19 @@ Description
     lift, drag and moment coefficients.  The data can optionally be output into
     bins, defined in a given direction.
 
+    The binned data provides the total and consitituentcomponents per bin:
+    - total coefficient
+    - pressure coefficient contribution
+    - viscous coefficient contribution
+    - porous coefficient contribution
+
+    Data is written into multiple files in the
+    postProcessing/\<functionObjectName\> directory:
+    - coefficient.dat   : integrated coefficients over all geometries
+    - CmBin.dat         : moment coefficient bins
+    - CdBin.dat         : drag coefficient bins
+    - ClBin.dat         : lift coefficient bins
+
     Example of function object specification:
     \verbatim
     forceCoeffs1
@@ -40,6 +53,7 @@ Description
         functionObjectLibs ("libforces.so");
         ...
         log         yes;
+        writeFields yes;
         patches     (walls);
         liftDir     (0 1 0);
         dragDir     (-1 0 0);
@@ -62,6 +76,7 @@ Description
         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 coefficient fields | no | no
         patches      | patches included in the forces calculation | yes |
         liftDir      | lift direction          | yes         |
         dragDir      | drag direction          | yes         |
@@ -137,6 +152,21 @@ class forceCoeffs
             scalar Aref_;
 
 
+        // File streams
+
+            //- Integrated coefficients
+            autoPtr<OFstream> coeffFilePtr_;
+
+            //- Moment coefficient
+            autoPtr<OFstream> CmBinFilePtr_;
+
+            //- Drag coefficient
+            autoPtr<OFstream> CdBinFilePtr_;
+
+            //- Lift coefficient
+            autoPtr<OFstream> ClBinFilePtr_;
+
+
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
@@ -148,8 +178,26 @@ class forceCoeffs
 
 protected:
 
-    //- Output file header information
+    // Protected Member Functions
-    virtual void writeFileHeader(const label i);
+
+        //- 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;
+
+        //- Write integrated data
+        void writeIntegratedData
+        (
+            const word& title,
+            const List<Field<scalar> >& coeff
+        ) const;
+
+        //- Write binned data
+        void writeBinData(const List<Field<scalar> > coeffs, Ostream& os) const;
 
 
 public:
@@ -167,7 +215,8 @@ public:
             const word& name,
             const objectRegistry&,
             const dictionary&,
-            const bool loadFromFiles = false
+            const bool loadFromFiles = false,
+            const bool readFields = true
         );
 
 

src/postProcessing/functionObjects/forces/forces/forces.C

@@ -43,113 +43,121 @@ namespace Foam
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-Foam::wordList Foam::forces::createFileNames(const dictionary& dict) const
+Foam::word Foam::forces::fieldName(const word& name) const
 {
-    DynamicList<word> names(1);
+    return name_ + ":" + name;
-
-    const word forceType(dict.lookup("type"));
-
-    if (dict.found("binData"))
-    {
-        const dictionary& binDict(dict.subDict("binData"));
-        label nb = readLabel(binDict.lookup("nBin"));
-        if (nb > 0)
-        {
-            names.append(forceType + "_bins");
-        }
-    }
-
-    names.append(forceType);
-
-    return names;
 }
 
 
-void Foam::forces::writeFileHeader(const label i)
+void Foam::forces::createFiles()
 {
-    if (i == 0)
+    // Note: Only possible to create bin files after bins have been initialised
+
+    if (writeToFile() && !forceFilePtr_.valid())
     {
-        // force data
+        forceFilePtr_ = createFile("force");
+        writeIntegratedHeader("Force", forceFilePtr_());
+        momentFilePtr_ = createFile("moment");
+        writeIntegratedHeader("Moment", momentFilePtr_());
 
-        writeHeader(file(i), "Forces");
+        if (nBin_ > 1)
-        writeHeaderValue(file(i), "CofR", coordSys_.origin());
+        {
-        writeCommented(file(i), "Time");
+            forceBinFilePtr_ = createFile("forceBin");
-
+            writeBinHeader("Force", forceBinFilePtr_());
-        file(i)
+            momentBinFilePtr_ = createFile("momentBin");
-            << "forces(pressure viscous porous) "
+            writeBinHeader("Moment", momentBinFilePtr_());
-            << "moment(pressure viscous porous)";
+        }
 
         if (localSystem_)
         {
-            file(i)
+            localForceFilePtr_ = createFile("localForce");
-                << tab
+            writeIntegratedHeader("Force", localForceFilePtr_());
-                << "localForces(pressure,viscous,porous) "
+            localMomentFilePtr_ = createFile("localMoment");
-                << "localMoments(pressure,viscous,porous)";
+            writeIntegratedHeader("Moment", localMomentFilePtr_());
-        }
-    }
-    else if (i == 1)
-    {
-        // bin data
 
-        writeHeader(file(i), "Force bins");
+            if (nBin_ > 1)
-        writeHeaderValue(file(i), "bins", nBin_);
+            {
-        writeHeaderValue(file(i), "start", binMin_);
+                localForceBinFilePtr_ = createFile("localForceBin");
-        writeHeaderValue(file(i), "delta", binDx_);
+                writeBinHeader("Force", localForceBinFilePtr_());
-        writeHeaderValue(file(i), "direction", binDir_);
+                localMomentBinFilePtr_ = createFile("localMomentBin");
+                writeBinHeader("Moment", localMomentBinFilePtr_());
+            }
+        }
+    }
+}
+
+
+void Foam::forces::writeIntegratedHeader
+(
+    const word& header,
+    Ostream& os
+) const
+{
+    writeHeader(os, header);
+    writeHeaderValue(os, "CofR", coordSys_.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::forces::writeBinHeader(const word& header, Ostream& os) const
+{
+    writeHeader(os, header + " bins");
+    writeHeaderValue(os, "bins", nBin_);
+    writeHeaderValue(os, "start", binMin_);
+    writeHeaderValue(os, "delta", binDx_);
+    writeHeaderValue(os, "direction", binDir_);
 
     vectorField binPoints(nBin_);
-        writeCommented(file(i), "x co-ords  :");
+    writeCommented(os, "x co-ords  :");
     forAll(binPoints, pointI)
     {
         binPoints[pointI] = (binMin_ + (pointI + 1)*binDx_)*binDir_;
-            file(i) << tab << binPoints[pointI].x();
+        os  << tab << binPoints[pointI].x();
     }
-        file(i) << nl;
+    os  << nl;
 
-        writeCommented(file(i), "y co-ords  :");
+    writeCommented(os, "y co-ords  :");
     forAll(binPoints, pointI)
     {
-            file(i) << tab << binPoints[pointI].y();
+        os  << tab << binPoints[pointI].y();
     }
-        file(i) << nl;
+    os  << nl;
 
-        writeCommented(file(i), "z co-ords  :");
+    writeCommented(os, "z co-ords  :");
     forAll(binPoints, pointI)
     {
-            file(i) << tab << binPoints[pointI].z();
+        os  << tab << binPoints[pointI].z();
     }
-        file(i) << nl;
+    os  << nl;
 
-        writeCommented(file(i), "Time");
+    writeHeader(os, "");
+    writeCommented(os, "Time");
 
     for (label j = 0; j < nBin_; j++)
     {
-            const word jn('(' + Foam::name(j) + ')');
+        const word jn(Foam::name(j) + ':');
-            const word f("forces" + jn + "[pressure,viscous,porous]");
+        os  << tab << jn << "(total_x total_y total_z)"
-            const word m("moments" + jn + "[pressure,viscous,porous]");
+            << tab << jn << "(pressure_x pressure_y pressure_z)"
+            << tab << jn << "(viscous_x viscous_y viscous_z)";
 
-            file(i)<< tab << f << tab << m;
+        if (porosity_)
-        }
-        if (localSystem_)
         {
-            for (label j = 0; j < nBin_; j++)
+            os  << tab << jn << "(porous_x porous_y porous_z)";
-            {
-                const word jn('(' + Foam::name(j) + ')');
-                const word f("localForces" + jn + "[pressure,viscous,porous]");
-                const word m("localMoments" + jn + "[pressure,viscous,porous]");
-
-                file(i)<< tab << f << tab << m;
         }
     }
-    }
-    else
-    {
-        FatalErrorIn("void Foam::forces::writeFileHeader(const label)")
-            << "Unhandled file index: " << i
-            << abort(FatalError);
-    }
 
-    file(i)<< endl;
+    os << endl;
 }
 
 
@@ -198,10 +206,117 @@ void Foam::forces::initialise()
         }
     }
 
+    initialiseBins();
+
     initialised_ = true;
 }
 
 
+void Foam::forces::initialiseBins()
+{
+    if (!active_)
+    {
+        return;
+    }
+
+    if (nBin_ > 1)
+    {
+        const fvMesh& mesh = refCast<const fvMesh>(obr_);
+        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+        // Determine extents of patches
+        binMin_ = GREAT;
+        scalar binMax = -GREAT;
+        forAllConstIter(labelHashSet, patchSet_, iter)
+        {
+            label patchI = iter.key();
+            const polyPatch& pp = pbm[patchI];
+            scalarField d(pp.faceCentres() & binDir_);
+            binMin_ = min(min(d), binMin_);
+            binMax = max(max(d), binMax);
+        }
+
+        // Include porosity
+        if (porosity_)
+        {
+            const HashTable<const porosityModel*> models =
+                obr_.lookupClass<porosityModel>();
+
+            const scalarField dd(mesh.C() & binDir_);
+
+            forAllConstIter(HashTable<const porosityModel*>, models, iter)
+            {
+                const porosityModel& pm = *iter();
+                const labelList& cellZoneIDs = pm.cellZoneIDs();
+
+                forAll(cellZoneIDs, i)
+                {
+                    label zoneI = cellZoneIDs[i];
+                    const cellZone& cZone = mesh.cellZones()[zoneI];
+                    const scalarField d(dd, cZone);
+                    binMin_ = min(min(d), binMin_);
+                    binMax = max(max(d), binMax);
+                }
+            }
+        }
+
+        reduce(binMin_, minOp<scalar>());
+        reduce(binMax, maxOp<scalar>());
+
+        // Slightly boost binMax so that region of interest is fully
+        // within bounds
+        binMax = 1.0001*(binMax - binMin_) + binMin_;
+
+        binDx_ = (binMax - binMin_)/scalar(nBin_);
+
+        // Create the bin 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_);
+            moment_[i].setSize(nBin_);
+        }
+    }
+}
+
+
+void Foam::forces::resetFields()
+{
+    force_[0] = vector::zero;
+    force_[1] = vector::zero;
+    force_[2] = vector::zero;
+
+    moment_[0] = vector::zero;
+    moment_[1] = vector::zero;
+    moment_[2] = vector::zero;
+
+    if (writeFields_)
+    {
+        volVectorField& force =
+            const_cast<volVectorField&>
+            (
+                obr_.lookupObject<volVectorField>(fieldName("force"))
+            );
+
+        force == dimensionedVector("0", force.dimensions(), vector::zero);
+
+        volVectorField& moment =
+            const_cast<volVectorField&>
+            (
+                obr_.lookupObject<volVectorField>(fieldName("moment"))
+            );
+
+        moment == dimensionedVector("0", moment.dimensions(), vector::zero);
+    }
+}
+
+
 Foam::tmp<Foam::volSymmTensorField> Foam::forces::devRhoReff() const
 {
     typedef compressible::turbulenceModel cmpTurbModel;
@@ -388,59 +503,180 @@ void Foam::forces::applyBins
 }
 
 
-void Foam::forces::writeForces()
+void Foam::forces::addToFields
+(
+    const label patchI,
+    const vectorField& Md,
+    const vectorField& fN,
+    const vectorField& fT,
+    const vectorField& fP
+)
 {
-    if (log_) Info
+    if (!writeFields_)
-        << type() << " " << name_ << " output:" << nl
-        << "    sum of forces:" << nl
-        << "        pressure : " << sum(force_[0]) << nl
-        << "        viscous  : " << sum(force_[1]) << nl
-        << "        porous   : " << sum(force_[2]) << nl
-        << "    sum of moments:" << nl
-        << "        pressure : " << sum(moment_[0]) << nl
-        << "        viscous  : " << sum(moment_[1]) << nl
-        << "        porous   : " << sum(moment_[2])
-        << endl;
-
-    file(0) << obr_.time().value() << tab << setw(1) << '('
-        << sum(force_[0]) << setw(1) << ' '
-        << sum(force_[1]) << setw(1) << ' '
-        << sum(force_[2]) << setw(3) << ") ("
-        << sum(moment_[0]) << setw(1) << ' '
-        << sum(moment_[1]) << setw(1) << ' '
-        << sum(moment_[2]) << setw(1) << ')'
-        << endl;
-
-    if (localSystem_)
-    {
-        vectorField localForceN(coordSys_.localVector(force_[0]));
-        vectorField localForceT(coordSys_.localVector(force_[1]));
-        vectorField localForceP(coordSys_.localVector(force_[2]));
-        vectorField localMomentN(coordSys_.localVector(moment_[0]));
-        vectorField localMomentT(coordSys_.localVector(moment_[1]));
-        vectorField localMomentP(coordSys_.localVector(moment_[2]));
-
-        file(0) << obr_.time().value() << tab << setw(1) << '('
-            << sum(localForceN) << setw(1) << ' '
-            << sum(localForceT) << setw(1) << ' '
-            << sum(localForceP) << setw(3) << ") ("
-            << sum(localMomentN) << setw(1) << ' '
-            << sum(localMomentT) << setw(1) << ' '
-            << sum(localMomentP) << setw(1) << ')'
-            << endl;
-    }
-}
-
-
-void Foam::forces::writeBins()
-{
-    if (nBin_ == 1)
     {
         return;
     }
 
-    List<Field<vector> > f(force_);
+    volVectorField& force =
-    List<Field<vector> > m(moment_);
+        const_cast<volVectorField&>
+        (
+            obr_.lookupObject<volVectorField>(fieldName("force"))
+        );
+
+    vectorField& pf = force.boundaryField()[patchI];
+    pf += fN + fT + fP;
+
+    volVectorField& moment =
+        const_cast<volVectorField&>
+        (
+            obr_.lookupObject<volVectorField>(fieldName("moment"))
+        );
+
+    vectorField& pm = moment.boundaryField()[patchI];
+    pm += Md;
+}
+
+
+void Foam::forces::addToFields
+(
+    const labelList& cellIDs,
+    const vectorField& Md,
+    const vectorField& fN,
+    const vectorField& fT,
+    const vectorField& fP
+)
+{
+    if (!writeFields_)
+    {
+        return;
+    }
+
+    volVectorField& force =
+        const_cast<volVectorField&>
+        (
+            obr_.lookupObject<volVectorField>(fieldName("force"))
+        );
+
+    volVectorField& moment =
+        const_cast<volVectorField&>
+        (
+            obr_.lookupObject<volVectorField>(fieldName("moment"))
+        );
+
+    forAll(cellIDs, i)
+    {
+        label cellI = cellIDs[i];
+        force[cellI] += fN[i] + fT[i] + fP[i];
+        moment[cellI] += Md[i];
+    }
+}
+
+
+void Foam::forces::writeIntegratedForceMoment
+(
+    const string& descriptor,
+    const vectorField& fm0,
+    const vectorField& fm1,
+    const vectorField& fm2,
+    autoPtr<OFstream>& osPtr
+) const
+{
+    vector pressure = sum(fm0);
+    vector viscous = sum(fm1);
+    vector porous = sum(fm2);
+    vector total = pressure + viscous + porous;
+
+    if (log_)
+    {
+        Info<< "    Sum of " << descriptor.c_str() << nl
+            << "        Total    : " << total << nl
+            << "        Pressure : " << pressure << nl
+            << "        Viscous  : " << viscous << nl;
+
+        if (porosity_)
+        {
+            Info<< "        Porous   : " << porous << nl;
+        }
+    }
+
+    if (writeToFile())
+    {
+        Ostream& os = osPtr();
+
+        os  << obr_.time().value()
+            << tab << total
+            << tab << pressure
+            << tab << viscous;
+
+        if (porosity_)
+        {
+            os  << tab << porous;
+        }
+
+        os  << endl;
+    }
+}
+
+
+void Foam::forces::writeForces()
+{
+    if (log_) Info << type() << " " << name_ << " output:" << nl;
+
+    writeIntegratedForceMoment
+    (
+        "forces",
+        force_[0],
+        force_[1],
+        force_[2],
+        forceFilePtr_
+    );
+
+    writeIntegratedForceMoment
+    (
+        "moments",
+        moment_[0],
+        moment_[1],
+        moment_[2],
+        momentFilePtr_
+    );
+
+    if (localSystem_)
+    {
+        writeIntegratedForceMoment
+        (
+            "local forces",
+            coordSys_.localVector(force_[0]),
+            coordSys_.localVector(force_[1]),
+            coordSys_.localVector(force_[2]),
+            localForceFilePtr_
+        );
+
+        writeIntegratedForceMoment
+        (
+            "local moments",
+            coordSys_.localVector(moment_[0]),
+            coordSys_.localVector(moment_[1]),
+            coordSys_.localVector(moment_[2]),
+            localMomentFilePtr_
+        );
+    }
+
+    if (log_) Info << endl;
+}
+
+
+void Foam::forces::writeBinnedForceMoment
+(
+    const List<Field<vector> >& fm,
+    autoPtr<OFstream>& osPtr
+) const
+{
+    if ((nBin_ == 1) || !writeToFile())
+    {
+        return;
+    }
+
+    List<Field<vector> > f(fm);
 
     if (binCumulative_)
     {
@@ -449,26 +685,35 @@ void Foam::forces::writeBins()
             f[0][i] += f[0][i-1];
             f[1][i] += f[1][i-1];
             f[2][i] += f[2][i-1];
-
-            m[0][i] += m[0][i-1];
-            m[1][i] += m[1][i-1];
-            m[2][i] += m[2][i-1];
         }
     }
 
-    file(1) << obr_.time().value();
+    Ostream& os = osPtr();
+
+    os  << obr_.time().value();
 
     forAll(f[0], i)
     {
-        file(1)
+        vector total = f[0][i] + f[1][i] + f[2][i];
-            << tab << setw(1) << '('
+
-            << f[0][i] << setw(1) << ' '
+        os  << tab << total
-            << f[1][i] << setw(1) << ' '
+            << tab << f[0][i]
-            << f[2][i] << setw(3) << ") ("
+            << tab << f[1][i];
-            << m[0][i] << setw(1) << ' '
+
-            << m[1][i] << setw(1) << ' '
+        if (porosity_)
-            << m[2][i] << setw(1) << ')';
+        {
+            os  << tab << f[2][i];
         }
+    }
+
+    os  << nl;
+}
+
+
+void Foam::forces::writeBins()
+{
+    writeBinnedForceMoment(force_, forceBinFilePtr_);
+    writeBinnedForceMoment(moment_, momentBinFilePtr_);
 
     if (localSystem_)
     {
@@ -481,33 +726,9 @@ void Foam::forces::writeBins()
         lm[1] = coordSys_.localVector(moment_[1]);
         lm[2] = coordSys_.localVector(moment_[2]);
 
-        if (binCumulative_)
+        writeBinnedForceMoment(lf, localForceBinFilePtr_);
-        {
+        writeBinnedForceMoment(lm, localMomentBinFilePtr_);
-            for (label i = 1; i < lf[0].size(); i++)
-            {
-                lf[0][i] += lf[0][i-1];
-                lf[1][i] += lf[1][i-1];
-                lf[2][i] += lf[2][i-1];
-                lm[0][i] += lm[0][i-1];
-                lm[1][i] += lm[1][i-1];
-                lm[2][i] += lm[2][i-1];
     }
-        }
-
-        forAll(lf[0], i)
-        {
-            file(1)
-                << tab << setw(1) << '('
-                << lf[0][i] << setw(1) << ' '
-                << lf[1][i] << setw(1) << ' '
-                << lf[2][i] << setw(3) << ") ("
-                << lm[0][i] << setw(1) << ' '
-                << lm[1][i] << setw(1) << ' '
-                << lm[2][i] << setw(1) << ')';
-        }
-    }
-
-    file(1) << endl;
 }
 
 
@@ -522,13 +743,21 @@ Foam::forces::forces
     const bool readFields
 )
 :
-    functionObjectFile(obr, name, createFileNames(dict)),
+    functionObjectFile(obr, name),
     name_(name),
     obr_(obr),
     active_(true),
     log_(true),
     force_(3),
     moment_(3),
+    forceFilePtr_(),
+    momentFilePtr_(),
+    forceBinFilePtr_(),
+    momentBinFilePtr_(),
+    localForceFilePtr_(),
+    localMomentFilePtr_(),
+    localForceBinFilePtr_(),
+    localMomentBinFilePtr_(),
     patchSet_(),
     pName_(word::null),
     UName_(word::null),
@@ -546,6 +775,7 @@ Foam::forces::forces
     binMin_(GREAT),
     binPoints_(),
     binCumulative_(true),
+    writeFields_(false),
     initialised_(false)
 {
     // Check if the available mesh is an fvMesh otherise deactivate
@@ -572,7 +802,6 @@ Foam::forces::forces
         )   << "No fvMesh available, deactivating " << name_
             << endl;
     }
-
 }
 
 
@@ -589,13 +818,20 @@ Foam::forces::forces
     const coordinateSystem& coordSys
 )
 :
-    functionObjectFile(obr, name, typeName),
+    functionObjectFile(obr, name),
-    name_(name),
     obr_(obr),
     active_(true),
     log_(true),
     force_(3),
     moment_(3),
+    forceFilePtr_(),
+    momentFilePtr_(),
+    forceBinFilePtr_(),
+    momentBinFilePtr_(),
+    localForceFilePtr_(),
+    localMomentFilePtr_(),
+    localForceBinFilePtr_(),
+    localMomentBinFilePtr_(),
     patchSet_(patchSet),
     pName_(pName),
     UName_(UName),
@@ -613,8 +849,12 @@ Foam::forces::forces
     binMin_(GREAT),
     binPoints_(),
     binCumulative_(true),
+    writeFields_(false),
     initialised_(false)
 {
+    // Turn off writing to file
+    writeToFile_ = false;
+
     forAll(force_, i)
     {
         force_[i].setSize(nBin_);
@@ -633,13 +873,18 @@ Foam::forces::~forces()
 
 void Foam::forces::read(const dictionary& dict)
 {
-    if (active_)
+    if (!active_)
     {
+        return;
+    }
+
+    functionObjectFile::read(dict);
+
     initialised_ = false;
 
     log_ = dict.lookupOrDefault<Switch>("log", false);
 
-        if (log_) Info<< type() << " " << name_ << ":" << nl;
+    if (log_) Info << type() << " " << name_ << ":" << nl;
 
     directForceDensity_ = dict.lookupOrDefault("directForceDensity", false);
 
@@ -680,11 +925,11 @@ void Foam::forces::read(const dictionary& dict)
     dict.readIfPresent("porosity", porosity_);
     if (porosity_)
     {
-            if (log_) Info<< "    Including porosity effects" << endl;
+        if (log_) Info << "    Including porosity effects" << endl;
     }
     else
     {
-            if (log_) Info<< "    Not including porosity effects" << endl;
+        if (log_) Info << "    Not including porosity effects" << endl;
     }
 
     if (dict.found("binData"))
@@ -700,9 +945,21 @@ void Foam::forces::read(const dictionary& dict)
             )   << "Number of bins (nBin) must be zero or greater"
                 << exit(FatalIOError);
         }
-            else if ((nBin_ == 0) || (nBin_ == 1))
+        else if (nBin_ == 0)
         {
             nBin_ = 1;
+        }
+        else
+        {
+            binDict.lookup("cumulative") >> binCumulative_;
+            binDict.lookup("direction") >> binDir_;
+            binDir_ /= mag(binDir_);
+        }
+    }
+
+    if (nBin_ == 1)
+    {
+        // Allocate storage for forces and moments
         forAll(force_, i)
         {
             force_[i].setSize(1);
@@ -710,82 +967,61 @@ void Foam::forces::read(const dictionary& dict)
         }
     }
 
-            if (nBin_ > 1)
+    writeFields_ = dict.lookupOrDefault("writeFields", false);
+
+    if (writeFields_)
     {
-                binDict.lookup("direction") >> binDir_;
+        if (log_) Info << "    Fields will be written" << endl;
-                binDir_ /= mag(binDir_);
 
-                binMin_ = GREAT;
+        tmp<volVectorField> tforce
-                scalar binMax = -GREAT;
+        (
-                forAllConstIter(labelHashSet, patchSet_, iter)
+            new volVectorField
-                {
+            (
-                    label patchI = iter.key();
+                IOobject
-                    const polyPatch& pp = pbm[patchI];
+                (
-                    scalarField d(pp.faceCentres() & binDir_);
+                    fieldName("force"),
-                    binMin_ = min(min(d), binMin_);
+                    mesh.time().timeName(),
-                    binMax = max(max(d), binMax);
+                    mesh,
-                }
+                    IOobject::NO_READ,
-                reduce(binMin_, minOp<scalar>());
+                    IOobject::NO_WRITE
-                reduce(binMax, maxOp<scalar>());
+                ),
+                mesh,
+                dimensionedVector("0", dimForce, vector::zero)
+            )
+        );
 
-                // slightly boost binMax so that region of interest is fully
+        obr_.store(tforce.ptr());
-                // within bounds
-                binMax = 1.0001*(binMax - binMin_) + binMin_;
 
-                binDx_ = (binMax - binMin_)/scalar(nBin_);
+        tmp<volVectorField> tmoment
+        (
+            new volVectorField
+            (
+                IOobject
+                (
+                    fieldName("moment"),
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                dimensionedVector("0", dimForce*dimLength, vector::zero)
+            )
+        );
 
-                // create the bin points used for writing
+        obr_.store(tmoment.ptr());
-                binPoints_.setSize(nBin_);
-                forAll(binPoints_, i)
-                {
-                    binPoints_[i] = (i + 0.5)*binDir_*binDx_;
-                }
-
-                binDict.lookup("cumulative") >> binCumulative_;
-
-                // allocate storage for forces and moments
-                forAll(force_, i)
-                {
-                    force_[i].setSize(nBin_);
-                    moment_[i].setSize(nBin_);
-                }
-            }
-        }
-
-        if (nBin_ == 1)
-        {
-            // allocate storage for forces and moments
-            force_[0].setSize(1);
-            force_[1].setSize(1);
-            force_[2].setSize(1);
-            moment_[0].setSize(1);
-            moment_[1].setSize(1);
-            moment_[2].setSize(1);
-        }
     }
 }
 
 
 void Foam::forces::execute()
 {
-    // Do nothing - only valid on write
+    if (!active_)
-}
+    {
+        return;
+    }
 
-
+    // calcForcesMoment may have reset the active flag - need to re-check
-void Foam::forces::end()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::forces::timeSet()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::forces::write()
-{
     calcForcesMoment();
 
     if (!active_)
@@ -795,33 +1031,64 @@ void Foam::forces::write()
 
     if (Pstream::master())
     {
-        functionObjectFile::write();
+        createFiles();
 
         writeForces();
 
         writeBins();
 
-        if (log_) Info<< endl;
+        if (log_) Info << endl;
+    }
+
+}
+
+
+void Foam::forces::end()
+{
+    if (active_)
+    {
+        execute();
+    }
+}
+
+
+void Foam::forces::timeSet()
+{
+    // Do nothing
+}
+
+
+void Foam::forces::write()
+{
+    if (!active_)
+    {
+        return;
+    }
+
+    if (writeFields_)
+    {
+        obr_.lookupObject<volVectorField>(fieldName("force")).write();
+        obr_.lookupObject<volVectorField>(fieldName("moment")).write();
     }
 }
 
 
 void Foam::forces::calcForcesMoment()
 {
-    initialise();
-
     if (!active_)
     {
         return;
     }
 
-    force_[0] = vector::zero;
+    initialise();
-    force_[1] = vector::zero;
-    force_[2] = vector::zero;
 
-    moment_[0] = vector::zero;
+    // Initialise may have reset the active flag - need to re-check
-    moment_[1] = vector::zero;
+    if (!active_)
-    moment_[2] = vector::zero;
+    {
+        return;
+    }
+
+    resetFields();
 
     if (directForceDensity_)
     {
@@ -858,6 +1125,8 @@ void Foam::forces::calcForcesMoment()
             //- Porous force
             vectorField fP(Md.size(), vector::zero);
 
+            addToFields(patchI, Md, fN, fT, fP);
+
             applyBins(Md, fN, fT, fP, mesh.C().boundaryField()[patchI]);
         }
     }
@@ -896,6 +1165,8 @@ void Foam::forces::calcForcesMoment()
 
             vectorField fP(Md.size(), vector::zero);
 
+            addToFields(patchI, Md, fN, fT, fP);
+
             applyBins(Md, fN, fT, fP, mesh.C().boundaryField()[patchI]);
         }
     }
@@ -939,6 +1210,8 @@ void Foam::forces::calcForcesMoment()
 
                 const vectorField fDummy(Md.size(), vector::zero);
 
+                addToFields(cZone, Md, fDummy, fDummy, fP);
+
                 applyBins(Md, fDummy, fDummy, fP, d);
             }
         }

src/postProcessing/functionObjects/forces/forces/forces.H

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,11 +29,24 @@ Group
 
 Description
     This function object calculates the forces and moments by integrating the
-    pressure and skin-friction forces over a given list of patches.
+    pressure and skin-friction forces over a given list of patches, and
+    the resistance from porous zones.
 
-    Member function forces::write() calculates the forces/moments and
+    Forces and moments are calculated, with optional co-ordinate system and
-    writes the forces/moments into the file \<timeDir\>/forces.dat and bin
+    writing of binned data, where force and moment contributions are collected
-    data (if selected) to the file \<timeDir\>/forces_bin.dat
+    into a user-defined number of bins that span the input geometries for a
+    user-defined direction vector.
+
+    Data is written into multiple files in the
+    postProcessing/\<functionObjectName\> directory:
+    - force.dat          : forces in global Cartesian co-ordinate system
+    - moment.dat         : moments in global Cartesian co-ordinate system
+    - forceBin.dat       : force bins in global Cartesian co-ordinate system
+    - momentBin.dat      : moment bins in global Cartesian co-ordinate system
+    - localForce.dat     : forces in local co-ordinate system
+    - localMoment.dat    : moments in local co-ordinate system
+    - localForceBin.dat  : force bins in local co-ordinate system
+    - localMomentBin.dat : moment bins in local co-ordinate system
 
     Example of function object specification:
     \verbatim
@@ -43,6 +56,7 @@ Description
         functionObjectLibs ("libforces.so");
         ...
         log         yes;
+        writeFields yes;
         patches     (walls);
 
         binData
@@ -59,6 +73,7 @@ Description
         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 |
         pName        | pressure field name     | no          | p
         UName        | velocity field name     | no          | U
@@ -101,6 +116,7 @@ Note
 
 SeeAlso
     Foam::functionObject
+    Foam::functionObjectFile
     Foam::OutputFilterFunctionObject
     Foam::forceCoeffs
 
@@ -122,7 +138,6 @@ SourceFiles
 #include "Tuple2.H"
 #include "OFstream.H"
 #include "Switch.H"
-#include "writer.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -165,6 +180,32 @@ protected:
         //- Pressure, viscous and porous moment per bin
         List<Field<vector> > moment_;
 
+        // File streams
+
+            //- Forces
+            autoPtr<OFstream> forceFilePtr_;
+
+            //- Moments
+            autoPtr<OFstream> momentFilePtr_;
+
+            //- Force bins
+            autoPtr<OFstream> forceBinFilePtr_;
+
+            //- Moment bins
+            autoPtr<OFstream> momentBinFilePtr_;
+
+            //- Local force
+            autoPtr<OFstream> localForceFilePtr_;
+
+            //- Local moment
+            autoPtr<OFstream> localMomentFilePtr_;
+
+            //- Local force bins
+            autoPtr<OFstream> localForceBinFilePtr_;
+
+            //- Local moment bins
+            autoPtr<OFstream> localMomentBinFilePtr_;
+
 
         // Read from dictionary
 
@@ -223,21 +264,36 @@ protected:
                 bool binCumulative_;
 
 
+            //- Write fields flag
+            bool writeFields_;
+
             //- Initialised flag
             bool initialised_;
 
 
     // Protected Member Functions
 
-        //- Create file names for forces and bins
+        //- Create a field name
-        wordList createFileNames(const dictionary& dict) const;
+        word fieldName(const word& name) const;
 
-        //- Output file header information
+        //- Create the output files
-        virtual void writeFileHeader(const label i);
+        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;
 
         //- Initialise the fields
         void initialise();
 
+        //- Initialise the collection bins
+        void initialiseBins();
+
+        //- Reset the fields prior to accumulation of force/moments
+        void resetFields();
+
         //- Return the effective viscous stress (laminar + turbulent).
         tmp<volSymmTensorField> devRhoReff() const;
 
@@ -261,10 +317,47 @@ protected:
             const vectorField& d
         );
 
-        //- Helper function to write force data
+        //- Add patch contributions to force and moment fields
+        void addToFields
+        (
+            const label patchI,
+            const vectorField& Md,
+            const vectorField& fN,
+            const vectorField& fT,
+            const vectorField& fP
+        );
+
+        //- Add cell contributions to force and moment fields
+        void addToFields
+        (
+            const labelList& cellIDs,
+            const vectorField& Md,
+            const vectorField& fN,
+            const vectorField& fT,
+            const vectorField& fP
+        );
+
+        //- 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();
 
-        //- Helper function to write bin data
+        //- Helper function to write binned forces and moments
+        void writeBinnedForceMoment
+        (
+            const List<Field<vector> >& fm,
+            autoPtr<OFstream>& osPtr
+        ) const;
+
+        //- Write binned data
         void writeBins();
 
         //- Disallow default bitwise copy construct