Adding post-processing routines to compute surface forces (#659)

* Adding new post-processing routines to compute surface forces

* Replacing previous version of forcePerS

src/utilities/postProcessing/calcForcePerSCompressible/calcForcePerS.C

@@ -11,7 +11,11 @@
 
 #include "fvCFD.H"
 #include "argList.H"
+#include "autoPtr.H"
 #include "Time.H"
+#include "timeSelector.H"
+#include "TimePaths.H"
+#include "ListOps.H"
 #include "fvMesh.H"
 #include "OFstream.H"
 #include "simpleControl.H"
@@ -24,40 +28,31 @@ using namespace Foam;
 
 int main(int argc, char* argv[])
 {
-    Info << "Computing forcePerS...." << endl;
+    Info << "Computing forces...." << endl;
 
     argList::addOption(
         "patchNames",
-        "'(inlet)'",
+        "'(wing)'",
         "List of patch names to compute");
 
-    argList::addOption(
-        "forceDir",
-        "'(0 0 1)'",
-        "Force direction");
-
     argList::addOption(
         "time",
         "1000",
-        "Tme instance to compute");
+        "Time instance to compute, if not provided runs all times");
 
 #include "setRootCase.H"
 #include "createTime.H"
 
-    scalar time;
-    if (args.optionFound("time"))
+    if (!args.optionFound("time"))
     {
-        time = readScalar(args.optionLookup("time")());
+        Info << "Time not set, running all times." << endl;
     }
-    else
-    {
-        Info << "time not set! Exit." << endl;
-        return 1;
-    }
-    runTime.setTime(time, 0);
 
-#include "createMesh.H"
-#include "createFields.H"
+    // Create the processor databases
+    autoPtr<TimePaths> timePaths;
+    timePaths = autoPtr<TimePaths>::New(args.rootPath(), args.caseName());
+
+    const instantList timeDirs(timeSelector::select(timePaths->times(), args));
 
     List<wordRe> patchNames;
     if (args.optionFound("patchNames"))
@@ -70,83 +65,82 @@ int main(int argc, char* argv[])
         return 1;
     }
 
-    List<scalar> forceDir1;
-    if (args.optionFound("forceDir"))
-    {
-        forceDir1 = scalarList(args.optionLookup("forceDir")());
-    }
-    else
+    forAll(timeDirs, iTime)
     {
-        Info << "forceDir not set! Exit." << endl;
-        return 1;
-    }
-    vector forceDir(vector::zero);
-    forceDir.x() = forceDir1[0];
-    forceDir.y() = forceDir1[1];
-    forceDir.z() = forceDir1[2];
-
-    volScalarField forcePerS(
-        IOobject(
-            "forcePerS",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE),
-        mesh,
-        dimensionedScalar("forcePerS", dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0),
-        fixedValueFvPatchScalarField::typeName);
-
-    // this code is pulled from:
-    // src/functionObjects/forcces/forces.C
-    // modified slightly
-    vector forces(vector::zero);
-
-    const surfaceVectorField::Boundary& Sfb = mesh.Sf().boundaryField();
-    const surfaceScalarField::Boundary& magSfb = mesh.magSf().boundaryField();
-
-    volSymmTensorField devRhoReff(
-        IOobject(
-            IOobject::groupName("devRhoReff", U.group()),
-            mesh.time().timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE),
-        (-rho * nuEff) * dev(twoSymm(fvc::grad(U))));
+        runTime.setTime(timeDirs[iTime].value(), 0);
 
-    const volSymmTensorField::Boundary& devRhoReffb = devRhoReff.boundaryField();
+#include "createMesh.H"
+#include "createFields.H"
 
-    forAll(patchNames, cI)
-    {
-        // get the patch id label
-        label patchI = mesh.boundaryMesh().findPatchID(patchNames[cI]);
-        if (patchI < 0)
-        {
-            Info << "ERROR: Patch name " << patchNames[cI] << " not found in constant/polyMesh/boundary! Exit!" << endl;
-            return 1;
-        }
-        // create a shorter handle for the boundary patch
-        const fvPatch& patch = mesh.boundary()[patchI];
-        // normal force
-        vectorField fN(
-            Sfb[patchI] * p.boundaryField()[patchI]);
-        // tangential force
-        vectorField fT(Sfb[patchI] & devRhoReffb[patchI]);
-        // sum them up
-        forAll(patch, faceI)
+        // Initialize surface field for face-centered forces
+        volVectorField forcePerS(
+            IOobject(
+                "forcePerS",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE),
+            mesh,
+            dimensionedVector("forcePerS", dimensionSet(1, -1, -2, 0, 0, 0, 0), vector::zero),
+            fixedValueFvPatchScalarField::typeName);
+
+        // this code is pulled from:
+        // src/functionObjects/forcces/forces.C
+        // modified slightly
+        vector forces(vector::zero);
+
+        const surfaceVectorField::Boundary& Sfb = mesh.Sf().boundaryField();
+        const surfaceScalarField::Boundary& magSfb = mesh.magSf().boundaryField();
+
+        volSymmTensorField devRhoReff(
+            IOobject(
+                IOobject::groupName("devRhoReff", U.group()),
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE),
+            (-rho * nuEff) * dev(twoSymm(fvc::grad(U))));
+
+        const volSymmTensorField::Boundary& devRhoReffb = devRhoReff.boundaryField();
+
+        vector totalForce(vector::zero);
+        forAll(patchNames, cI)
         {
-            forces.x() = fN[faceI].x() + fT[faceI].x();
-            forces.y() = fN[faceI].y() + fT[faceI].y();
-            forces.z() = fN[faceI].z() + fT[faceI].z();
-            scalar force = forces & forceDir;
-            forcePerS.boundaryFieldRef()[patchI][faceI] = force / magSfb[patchI][faceI];
+            // get the patch id label
+            label patchI = mesh.boundaryMesh().findPatchID(patchNames[cI]);
+            if (patchI < 0)
+            {
+                Info << "ERROR: Patch name " << patchNames[cI] << " not found in constant/polyMesh/boundary! Exit!" << endl;
+                return 1;
+            }
+            // create a shorter handle for the boundary patch
+            const fvPatch& patch = mesh.boundary()[patchI];
+            // normal force
+            vectorField fN(Sfb[patchI] * p.boundaryField()[patchI]);
+            // tangential force
+            vectorField fT(Sfb[patchI] & devRhoReffb[patchI]);
+            // sum them up
+            forAll(patch, faceI)
+            {
+                // compute forces
+                forces.x() = fN[faceI].x() + fT[faceI].x();
+                forces.y() = fN[faceI].y() + fT[faceI].y();
+                forces.z() = fN[faceI].z() + fT[faceI].z();
+
+                // project force direction
+                forcePerS.boundaryFieldRef()[patchI][faceI] = forces / magSfb[patchI][faceI];
+
+                totalForce.x() += forces.x();
+                totalForce.y() += forces.y();
+                totalForce.z() += forces.z();
+            }
         }
-    }
-    forcePerS.write();
-
-    Info << "Force: " << forces << endl;
+        forcePerS.write();
 
-    Info << "Computing forcePerS.... Completed!" << endl;
+        Info << "Total force: " << totalForce << endl;
 
+        Info << "Computing forces.... Completed!" << endl;
+    }
     return 0;
 }