diff --git a/src/utilities/postProcessing/calcForcePerSCompressible/calcForcePerS.C b/src/utilities/postProcessing/calcForcePerSCompressible/calcForcePerS.C index d2021d0d..750283ff 100755 --- a/src/utilities/postProcessing/calcForcePerSCompressible/calcForcePerS.C +++ b/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 = autoPtr::New(args.rootPath(), args.caseName()); + + const instantList timeDirs(timeSelector::select(timePaths->times(), args)); List patchNames; if (args.optionFound("patchNames")) @@ -70,83 +65,82 @@ int main(int argc, char* argv[]) return 1; } - List 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; } diff --git a/src/utilities/postProcessing/calcForcePerSIncompressible/calcForcePerS.C b/src/utilities/postProcessing/calcForcePerSIncompressible/calcForcePerS.C index f004e293..781891b2 100755 --- a/src/utilities/postProcessing/calcForcePerSIncompressible/calcForcePerS.C +++ b/src/utilities/postProcessing/calcForcePerSIncompressible/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,36 @@ 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"); + + argList::addOption( + "rho", + "1.0", + "Freestream density"); #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 = autoPtr::New(args.rootPath(), args.caseName()); + + const instantList timeDirs(timeSelector::select(timePaths->times(), args)); List patchNames; if (args.optionFound("patchNames")) @@ -70,83 +70,92 @@ int main(int argc, char* argv[]) return 1; } - List forceDir1; - if (args.optionFound("forceDir")) + scalar rho = 1.0; + if (args.optionFound("rho")) { - forceDir1 = scalarList(args.optionLookup("forceDir")()); + rho = readScalar(args.optionLookup("rho")()); } else { - Info << "forceDir not set! Exit." << endl; - return 1; + Info << "Density not set! Defaulting to 1.0." << endl; } - 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), - -nuEff * dev(twoSymm(fvc::grad(U)))); - - const volSymmTensorField::Boundary& devRhoReffb = devRhoReff.boundaryField(); - forAll(patchNames, cI) + forAll(timeDirs, iTime) { - // 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) + runTime.setTime(timeDirs[iTime].value(), 0); + +#include "createMesh.H" +#include "createFields.H" + + // 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), + -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] * rho); + // tangential force + vectorField fT(Sfb[patchI] & devRhoReffb[patchI] * rho); + // 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; }