diff --git a/src/postProcessing/functionObjects/forces/forces/forces.C b/src/postProcessing/functionObjects/forces/forces/forces.C index cf835d66ef..318e732b31 100644 --- a/src/postProcessing/functionObjects/forces/forces/forces.C +++ b/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::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(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 > binWriterPtr(writer::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 > localForce(2); + List > 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(obr_); + const polyBoundaryMesh& pbm = mesh.boundaryMesh(); - patchSet_ = mesh.boundaryMesh().patchSet - ( - wordReList(dict.lookup("patches")) - ); + patchSet_ = pbm.patchSet(wordReList(dict.lookup("patches"))); if (directForceDensity_) { @@ -334,68 +498,77 @@ void Foam::forces::read(const dictionary& dict) coordSys_ = coordinateSystem(dict, obr_); localSystem_ = true; } - } -} - -void Foam::forces::makeFile() -{ - // Create the forces file if not already created - if (forcesFilePtr_.empty()) - { - if (debug) + // read bin information if present + if (dict.readIfPresent