ENH: forces function object - added data output into bins

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

@@ -43,7 +43,78 @@ License
 defineTypeNameAndDebug(Foam::forces, 0);
 
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+Foam::fileName Foam::forces::baseFileDir() const
+{
+    fileName baseDir;
+
+    if (Pstream::parRun())
+    {
+        // Put in undecomposed case (Note: gives problems for
+        // distributed data running)
+        baseDir = obr_.time().path()/".."/name_;
+    }
+    else
+    {
+        baseDir = obr_.time().path()/name_;
+    }
+
+    return baseDir;
+}
+
+
+void Foam::forces::makeFile()
+{
+    // Create the forces file if not already created
+    if (forcesFilePtr_.empty())
+    {
+        if (debug)
+        {
+            Info<< "Creating forces file" << endl;
+        }
+
+        // File update
+        if (Pstream::master())
+        {
+            word startTimeName =
+                obr_.time().timeName(obr_.time().startTime().value());
+
+            fileName forcesDir = baseFileDir()/startTimeName;
+
+            // Create directory if does not exist.
+            mkDir(forcesDir);
+
+            // Open new file at start up
+            forcesFilePtr_.reset(new OFstream(forcesDir/(type() + ".dat")));
+
+            // Add headers to output data
+            writeFileHeader();
+        }
+    }
+}
+
+
+void Foam::forces::writeFileHeader()
+{
+    if (forcesFilePtr_.valid())
+    {
+        forcesFilePtr_()
+            << "# Time" << tab
+            << "forces(pressure, viscous) moment(pressure, viscous)";
+
+        if (localSystem_)
+        {
+            forcesFilePtr_()
+                << tab
+                << "local forces(pressure, viscous) "
+                << "local moment(pressure, viscous)";
+        }
+
+        forcesFilePtr_()<< endl;
+    }
+}
+
 
 Foam::tmp<Foam::volSymmTensorField> Foam::forces::devRhoReff() const
 {
@@ -167,6 +238,85 @@ Foam::scalar Foam::forces::rho(const volScalarField& p) const
 }
 
 
+void Foam::forces::applyBins
+(
+    const label patchI,
+    const vectorField fN,
+    const vectorField Md,
+    const vectorField fT
+)
+{
+    if (nBin_ == 1)
+    {
+        force_[0][0] = sum(fN);
+        force_[1][0] = sum(fT);
+        moment_[0][0] = sum(Md ^ fN);
+        moment_[1][0] = sum(Md ^ fT);
+    }
+    else
+    {
+        const fvMesh& mesh = refCast<const fvMesh>(obr_);
+        const vectorField& Cf = mesh.C().boundaryField()[patchI];
+        scalarField d((Cf & binDir_) - binMin_);
+
+        forAll(d, i)
+        {
+            label binI = floor(d[i]/binDx_);
+            force_[0][binI] += fN[i];
+            force_[1][binI] += fT[i];
+            moment_[0][binI] += Md[i]^fN[i];
+            moment_[1][binI] += Md[i]^fT[i];
+        }
+    }
+}
+
+
+void Foam::forces::writeBins() const
+{
+    if (nBin_ == 1)
+    {
+        return;
+    }
+
+    autoPtr<writer<vector> > binWriterPtr(writer<vector>::New(binFormat_));
+    coordSet axis("forces", "distance", binPoints_, mag(binPoints_));
+
+    fileName forcesDir = baseFileDir()/"bins"/obr_.time().timeName();
+    mkDir(forcesDir);
+
+    if (log_)
+    {
+        Info<< "    Writing bins to " << forcesDir << endl;
+
+    }
+
+    wordList fieldNames(IStringStream("(pressure viscous)")());
+
+    OFstream osForce(forcesDir/"force");
+    binWriterPtr->write(axis, fieldNames, force_, osForce);
+
+    OFstream osMoment(forcesDir/"moment");
+    binWriterPtr->write(axis, fieldNames, moment_, osMoment);
+
+
+    if (localSystem_)
+    {
+        List<Field<vector> > localForce(2);
+        List<Field<vector> > localMoment(2);
+        localForce[0] = coordSys_.localVector(force_[0]);
+        localForce[1] = coordSys_.localVector(force_[1]);
+        localMoment[0] = coordSys_.localVector(moment_[0]);
+        localMoment[1] = coordSys_.localVector(moment_[1]);
+
+        OFstream osLocalForce(forcesDir/"force_local");
+        binWriterPtr->write(axis, fieldNames, localForce, osLocalForce);
+
+        OFstream osLocalMoment(forcesDir/"moment_local");
+        binWriterPtr->write(axis, fieldNames, localMoment, osLocalMoment);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::forces::forces
@@ -181,6 +331,8 @@ Foam::forces::forces
     obr_(obr),
     active_(true),
     log_(false),
+    force_(2),
+    moment_(2),
     patchSet_(),
     pName_(word::null),
     UName_(word::null),
@@ -191,6 +343,12 @@ Foam::forces::forces
     pRef_(0),
     coordSys_(),
     localSystem_(false),
+    nBin_(1),
+    binDir_(vector::zero),
+    binDx_(0.0),
+    binMin_(GREAT),
+    binPoints_(),
+    binFormat_("undefined"),
     forcesFilePtr_(NULL)
 {
     // Check if the available mesh is an fvMesh otherise deactivate
@@ -231,6 +389,8 @@ Foam::forces::forces
     obr_(obr),
     active_(true),
     log_(false),
+    force_(2),
+    moment_(2),
     patchSet_(patchSet),
     pName_(pName),
     UName_(UName),
@@ -241,6 +401,12 @@ Foam::forces::forces
     pRef_(pRef),
     coordSys_(coordSys),
     localSystem_(false),
+    nBin_(1),
+    binDir_(vector::zero),
+    binDx_(0.0),
+    binMin_(GREAT),
+    binPoints_(),
+    binFormat_("undefined"),
     forcesFilePtr_(NULL)
 {}
 
@@ -261,11 +427,9 @@ void Foam::forces::read(const dictionary& dict)
         directForceDensity_ = dict.lookupOrDefault("directForceDensity", false);
 
         const fvMesh& mesh = refCast<const fvMesh>(obr_);
+        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
 
-        patchSet_ = mesh.boundaryMesh().patchSet
+        patchSet_ = pbm.patchSet(wordReList(dict.lookup("patches")));
-        (
-            wordReList(dict.lookup("patches"))
-        );
 
         if (directForceDensity_)
         {
@@ -334,68 +498,77 @@ void Foam::forces::read(const dictionary& dict)
             coordSys_ = coordinateSystem(dict, obr_);
             localSystem_ = true;
         }
+
+        // read bin information if present
+        if (dict.readIfPresent<label>("nBin", nBin_))
+        {
+            if (nBin_ < 0)
+            {
+                FatalIOErrorIn
+                (
+                    "void Foam::forces::read(const dictionary&)", dict
+                )   << "Number of bins (nBin) must be zero or greater"
+                    << exit(FatalIOError);
             }
-}
+            else if ((nBin_ == 0) || (nBin_ == 1))
-
-
-void Foam::forces::makeFile()
-{
-    // Create the forces file if not already created
-    if (forcesFilePtr_.empty())
             {
-        if (debug)
+                nBin_ = 1;
-        {
+                force_[0].setSize(1);
-            Info<< "Creating forces file." << endl;
+                force_[1].setSize(1);
+                moment_[0].setSize(1);
+                moment_[1].setSize(1);
             }
 
-        // File update
+            if (nBin_ > 1)
-        if (Pstream::master())
             {
-            fileName forcesDir;
+                dict.lookup("binDir") >> binDir_;
-            word startTimeName =
+                binDir_ /= mag(binDir_);
-                obr_.time().timeName(obr_.time().startTime().value());
 
-            if (Pstream::parRun())
+                binMin_ = GREAT;
+                scalar binMax = -GREAT;
+                forAllConstIter(labelHashSet, patchSet_, iter)
                 {
-                // Put in undecomposed case (Note: gives problems for
+                    label patchI = iter.key();
-                // distributed data running)
+                    const polyPatch& pp = pbm[patchI];
-                forcesDir = obr_.time().path()/".."/name_/startTimeName;
+                    scalarField d(pp.faceCentres() & binDir_);
+                    binMin_ = min(min(d), binMin_);
+                    binMax = max(max(d), binMax);
                 }
-            else
+                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)
                 {
-                forcesDir = obr_.time().path()/name_/startTimeName;
+                    binPoints_[i] = (i + 0.5)*binDir_*binDx_;
                 }
 
-            // Create directory if does not exist.
+                dict.lookup("binFormat") >> binFormat_;
-            mkDir(forcesDir);
 
-            // Open new file at start up
+                // allocate storage for forces and moments
-            forcesFilePtr_.reset(new OFstream(forcesDir/(type() + ".dat")));
+                forAll(force_, i)
-
-            // Add headers to output data
-            writeFileHeader();
-        }
-    }
-}
-
-
-void Foam::forces::writeFileHeader()
-{
-    if (forcesFilePtr_.valid())
                 {
-        forcesFilePtr_()
+                    force_[i].setSize(nBin_);
-            << "# Time" << tab
+                    moment_[i].setSize(nBin_);
-            << "forces(pressure, viscous) moment(pressure, viscous)";
+                }
-
+            }
-        if (localSystem_)
-        {
-            forcesFilePtr_()
-                << tab
-                << "local forces(pressure, viscous) "
-                << "local moment(pressure, viscous)";
         }
 
-        forcesFilePtr_()<< endl;
+        if (nBin_ == 1)
+        {
+            // allocate storage for forces and moments
+            force_[0].setSize(1);
+            force_[1].setSize(1);
+            moment_[0].setSize(1);
+            moment_[1].setSize(1);
+        }
     }
 }
 
@@ -414,71 +587,65 @@ void Foam::forces::end()
 
 void Foam::forces::write()
 {
-    if (active_)
+    if (!active_)
     {
+        return;
+    }
+
     // Create the forces file if not already created
     makeFile();
 
-        forcesMoments fm = calcForcesMoment();
+    calcForcesMoment();
 
     if (Pstream::master())
     {
         if (log_)
         {
-                Info<< "forces output:" << nl
+            Info<< type() << " output:" << nl
-                    << "    forces(pressure, viscous)" << fm.first() << nl
+                << "    forces(pressure,viscous)"
-                    << "    moment(pressure, viscous)" << fm.second() << nl;
+                << "(" << sum(force_[0]) << "," << sum(force_[1]) << ")" << nl
+                << "    moment(pressure,viscous)"
+                << "(" << sum(moment_[0]) << "," << sum(moment_[1]) << ")"
+                << nl;
         }
 
-            forcesFilePtr_() << obr_.time().value() << tab << fm;
+        forcesFilePtr_() << obr_.time().value() << tab
+            << "(" << sum(force_[0]) << "," << sum(force_[1]) << ") "
+            << "(" << sum(moment_[0]) << "," << sum(moment_[1]) << ")"
+            << endl;
 
         if (localSystem_)
         {
-                forcesMoments fmLocal;
+            vectorField localForceP(coordSys_.localVector(force_[0]));
+            vectorField localForceV(coordSys_.localVector(force_[1]));
+            vectorField localMomentP(coordSys_.localVector(moment_[0]));
+            vectorField localMomentV(coordSys_.localVector(moment_[1]));
 
-                fmLocal.first().first() =
+            forcesFilePtr_() << obr_.time().value() << tab
-                    coordSys_.localVector(fm.first().first());
+                << "(" << sum(localForceP) << "," << sum(localForceV) << ") "
-
+                << "(" << sum(localMomentP) << "," << sum(localMomentV) << ")"
-                fmLocal.first().second() =
+                << endl;
-                    coordSys_.localVector(fm.first().second());
-
-                fmLocal.second().first() =
-                    coordSys_.localVector(fm.second().first());
-
-                fmLocal.second().second() =
-                    coordSys_.localVector(fm.second().second());
-
-                forcesFilePtr_() << tab << fmLocal;
-
-                if (log_)
-                {
-                    Info<< "  local:" << nl
-                        << "    forces(pressure, viscous)" << fmLocal.first()
-                        << nl
-                        << "    moment(pressure, viscous)" << fmLocal.second()
-                        << nl;
-                }
         }
 
-
+        writeBins();
-            forcesFilePtr_() << endl;
 
         if (log_)
         {
             Info<< endl;
         }
-        }
+
+        forcesFilePtr_() << endl;
     }
 }
 
 
-Foam::forces::forcesMoments Foam::forces::calcForcesMoment() const
+void Foam::forces::calcForcesMoment()
 {
-    forcesMoments fm
+    force_[0] = vector::zero;
-    (
+    force_[1] = vector::zero;
-        pressureViscous(vector::zero, vector::zero),
+
-        pressureViscous(vector::zero, vector::zero)
+    moment_[0] = vector::zero;
-    );
+    moment_[1] = vector::zero;
 
     if (directForceDensity_)
     {
@@ -491,32 +658,28 @@ Foam::forces::forcesMoments Foam::forces::calcForcesMoment() const
 
         forAllConstIter(labelHashSet, patchSet_, iter)
         {
-            label patchi = iter.key();
+            label patchI = iter.key();
 
             vectorField Md
             (
-                mesh.C().boundaryField()[patchi] - coordSys_.origin()
+                mesh.C().boundaryField()[patchI] - coordSys_.origin()
             );
 
-            scalarField sA(mag(Sfb[patchi]));
+            scalarField sA(mag(Sfb[patchI]));
 
-            // Normal force = surfaceUnitNormal * (surfaceNormal & forceDensity)
+            // Normal force = surfaceUnitNormal*(surfaceNormal & forceDensity)
             vectorField fN
             (
-                Sfb[patchi]/sA
+                Sfb[patchI]/sA
                *(
-                    Sfb[patchi] & fD.boundaryField()[patchi]
+                    Sfb[patchI] & fD.boundaryField()[patchI]
                 )
             );
 
-            fm.first().first() += sum(fN);
-            fm.second().first() += sum(Md ^ fN);
-
             // Tangential force (total force minus normal fN)
-            vectorField fT(sA*fD.boundaryField()[patchi] - fN);
+            vectorField fT(sA*fD.boundaryField()[patchI] - fN);
 
-            fm.first().second() += sum(fT);
+            applyBins(patchI, fN, Md, fT);
-            fm.second().second() += sum(Md ^ fT);
         }
     }
     else
@@ -538,28 +701,38 @@ Foam::forces::forcesMoments Foam::forces::calcForcesMoment() const
 
         forAllConstIter(labelHashSet, patchSet_, iter)
         {
-            label patchi = iter.key();
+            label patchI = iter.key();
 
             vectorField Md
             (
-                mesh.C().boundaryField()[patchi] - coordSys_.origin()
+                mesh.C().boundaryField()[patchI] - coordSys_.origin()
             );
 
-            vectorField pf(Sfb[patchi]*(p.boundaryField()[patchi] - pRef));
+            vectorField pf
+            (
+                rho(p)*Sfb[patchI]*(p.boundaryField()[patchI] - pRef)
+            );
 
-            fm.first().first() += rho(p)*sum(pf);
+            vectorField vf(Sfb[patchI] & devRhoReffb[patchI]);
-            fm.second().first() += rho(p)*sum(Md ^ pf);
 
-            vectorField vf(Sfb[patchi] & devRhoReffb[patchi]);
+            applyBins(patchI, pf, Md, vf);
-
-            fm.first().second() += sum(vf);
-            fm.second().second() += sum(Md ^ vf);
         }
     }
 
-    reduce(fm, sumOp());
+    Pstream::listCombineGather(force_, plusEqOp<vectorField>());
+    Pstream::listCombineGather(moment_, plusEqOp<vectorField>());
+}
 
-    return fm;
+
+Foam::vector Foam::forces::forceEff() const
+{
+    return sum(force_[0]) + sum(force_[1]);
+}
+
+
+Foam::vector Foam::forces::momentEff() const
+{
+    return sum(moment_[0]) + sum(moment_[1]);
 }
 
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,6 +31,13 @@ Description
     Member function forces::write() calculates the forces/moments and
     writes the forces/moments into the file \<timeDir\>/forces.dat
 
+    The data can optionally be collected into bins, using e.g.
+    \verbaitim
+        nBin        20;          // output data into bins
+        binDir      (1 0 0);     // bin direction
+        binFormat   gnuplot;
+    \endverbatim
+
 Note
     The centre of rotation for moment calculations can either be specified
     by an \c CofR entry, or be taken from origin of the local coordinateSystem.
@@ -66,6 +73,7 @@ SourceFiles
 #include "OFstream.H"
 #include "Switch.H"
 #include "pointFieldFwd.H"
+#include "writer.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -83,46 +91,9 @@ class mapPolyMesh;
 
 class forces
 {
-public:
-
-    // Tuple for pressure (.first()) and viscous (.second()) forces
-    typedef Tuple2<vector, vector> pressureViscous;
-
-    // Tuple for forces (.first()) and moment (.second())
-    // pressure/viscous forces Tuples.
-    typedef Tuple2<pressureViscous, pressureViscous> forcesMoments;
-
-    //- Sum operation class to accumulate pressure/viscous forces and moments
-    class sumOp
-    {
-    public:
-
-        forcesMoments operator()
-        (
-            const forcesMoments& fm1,
-            const forcesMoments& fm2
-        ) const
-        {
-            return forcesMoments
-            (
-                pressureViscous
-                (
-                    fm1.first().first()  + fm2.first().first(),
-                    fm1.first().second() + fm2.first().second()
-                ),
-                pressureViscous
-                (
-                    fm1.second().first()  + fm2.second().first(),
-                    fm1.second().second() + fm2.second().second()
-                )
-            );
-        }
-    };
-
-
 protected:
 
-    // Private data
+    // Protected data
 
         //- Name of this set of forces,
         //  Also used as the name of the probes directory.
@@ -130,12 +101,19 @@ protected:
 
         const objectRegistry& obr_;
 
-        //- on/off switch
+        //- On/off switch
         bool active_;
 
         //- Switch to send output to Info as well as to file
         Switch log_;
 
+        //- Pressure and viscous force per bin
+        List<Field<vector> > force_;
+
+        //- Pressure and viscous pressure per bin
+        List<Field<vector> > moment_;
+
+
         // Read from dictionary
 
             //- Patches to integrate forces over
@@ -169,11 +147,35 @@ protected:
             bool localSystem_;
 
 
+            // Bin information
+
+                //- Number of bins
+                label nBin_;
+
+                //- Direction used to determine bin orientation
+                vector binDir_;
+
+                //- Distance between bin divisions
+                scalar binDx_;
+
+                //- Minimum bin bounds
+                scalar binMin_;
+                
+                //- Bin positions along binDir
+                List<point> binPoints_;
+
+                //- Write format for bin data
+                word binFormat_;
+
+
         //- Forces/moment file ptr
         autoPtr<OFstream> forcesFilePtr_;
 
 
-    // Private Member Functions
+    // Protected Member Functions
+
+        //- Return the base file directory for output
+        fileName baseFileDir() const;
 
         //- If the forces file has not been created create it
         void makeFile();
@@ -191,6 +193,18 @@ protected:
         //  otherwise return 1
         scalar rho(const volScalarField& p) const;
 
+        //- Accumulate bin data
+        void applyBins
+        (
+            const label patchI,
+            const vectorField fN,
+            const vectorField Md,
+            const vectorField fT
+        );
+
+        //- Helper function to write bin data
+        void writeBins() const;
+
         //- Disallow default bitwise copy construct
         forces(const forces&);
 
@@ -256,8 +270,14 @@ public:
         //- Write the forces
         virtual void write();
 
-        //- Calculate and return forces and moment
+        //- Calculate the forces and moments
-        virtual forcesMoments calcForcesMoment() const;
+        virtual void calcForcesMoment();
+
+        //- Return the total force
+        virtual vector forceEff() const;
+
+        //- Return the total moment
+        virtual vector momentEff() const;
 
         //- Update for changes of mesh
         virtual void updateMesh(const mapPolyMesh&)