Added to the forces functionObject the ability to read a force density (fD) volVectorField directly rather than a p and U field. To use, the following entries go in the controlDict directForceDensity true; fDName fDMean;

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

@@ -151,6 +151,8 @@ Foam::forces::forces
     patchSet_(),
     pName_(""),
     UName_(""),
+    directForceDensity_(false),
+    fDName_(""),
     rhoRef_(0),
     CofR_(vector::zero),
     forcesFilePtr_(NULL)
@@ -189,6 +191,28 @@ void Foam::forces::read(const dictionary& dict)
         patchSet_ =
             mesh.boundaryMesh().patchSet(wordList(dict.lookup("patches")));
 
+        dict.readIfPresent("directForceDensity", directForceDensity_);
+
+        if (directForceDensity_)
+        {
+            // Optional entry for fDName
+            fDName_ = dict.lookupOrDefault<word>("fDName", "fD");
+
+            // Check whether fDName exists, if not deactivate forces
+            if
+            (
+                !obr_.foundObject<volVectorField>(fDName_)
+            )
+            {
+                active_ = false;
+                WarningIn("void forces::read(const dictionary& dict)")
+                << "Could not find " << fDName_ << " in database." << nl
+                    << "    De-activating forces."
+                    << endl;
+            }
+        }
+        else
+        {
             // Optional entries U and p
             pName_ = dict.lookupOrDefault<word>("pName", "p");
             UName_ = dict.lookupOrDefault<word>("UName", "U");
@@ -210,6 +234,7 @@ void Foam::forces::read(const dictionary& dict)
 
             // Reference density needed for incompressible calculations
             rhoRef_ = readScalar(dict.lookup("rhoInf"));
+        }
 
         // Centre of rotation for moment calculations
         CofR_ = dict.lookup("CofR");
@@ -307,17 +332,53 @@ void Foam::forces::write()
 
 Foam::forces::forcesMoments Foam::forces::calcForcesMoment() const
 {
-    const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
-    const volScalarField& p = obr_.lookupObject<volScalarField>(pName_);
-
-    const fvMesh& mesh = U.mesh();
-
     forcesMoments fm
     (
         pressureViscous(vector::zero, vector::zero),
         pressureViscous(vector::zero, vector::zero)
     );
 
+    if (directForceDensity_)
+    {
+        const volVectorField& fD = obr_.lookupObject<volVectorField>(fDName_);
+
+        const fvMesh& mesh = fD.mesh();
+
+        const surfaceVectorField::GeometricBoundaryField& Sfb =
+            mesh.Sf().boundaryField();
+
+        forAllConstIter(labelHashSet, patchSet_, iter)
+        {
+            label patchi = iter.key();
+
+            vectorField Md = mesh.C().boundaryField()[patchi] - CofR_;
+
+            scalarField sA = mag(Sfb[patchi]);
+
+            // Normal force = surfaceUnitNormal * (surfaceNormal & forceDensity)
+            vectorField fN =
+                Sfb[patchi]/sA
+               *(
+                    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;
+
+            fm.first().second() += sum(fT);
+            fm.second().second() += sum(Md ^ fT);
+        }
+    }
+    else
+    {
+        const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
+        const volScalarField& p = obr_.lookupObject<volScalarField>(pName_);
+
+        const fvMesh& mesh = U.mesh();
+
         const surfaceVectorField::GeometricBoundaryField& Sfb =
             mesh.Sf().boundaryField();
 
@@ -331,8 +392,7 @@ Foam::forces::forcesMoments Foam::forces::calcForcesMoment() const
 
             vectorField Md = mesh.C().boundaryField()[patchi] - CofR_;
 
-        vectorField pf =
+            vectorField pf = Sfb[patchi]*p.boundaryField()[patchi];
-            mesh.Sf().boundaryField()[patchi]*p.boundaryField()[patchi];
 
             fm.first().first() += rho(p)*sum(pf);
             fm.second().first() += rho(p)*sum(Md ^ pf);
@@ -342,6 +402,7 @@ Foam::forces::forcesMoments Foam::forces::calcForcesMoment() const
             fm.first().second() += sum(vf);
             fm.second().second() += sum(Md ^ vf);
         }
+    }
 
     reduce(fm, sumOp());
 

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

@@ -121,7 +121,7 @@ protected:
         //- Switch to send output to Info as well as to file
         Switch log_;
 
-        // Read from dictonary
+        // Read from dictionary
 
             //- Patches to integrate forces over
             labelHashSet patchSet_;
@@ -132,6 +132,12 @@ protected:
             //- Name of velocity field
             word UName_;
 
+            //- Is the force density being supplied directly?
+            Switch directForceDensity_;
+
+            //- The name of the force density (fD) field
+            word fDName_;
+
             //- Reference density needed for incompressible calculations
             scalar rhoRef_;