From 45881823522b587d59f0e0a944e0a509cf62fa93 Mon Sep 17 00:00:00 2001 From: Mark Olesen Date: Fri, 11 Aug 2017 14:53:24 +0200 Subject: [PATCH 01/13] ENH: add absolute weighting for surfaceFieldValue (issue #567) - can be useful either for flow-rate weighting where backflow is to be ignored in the average, or for flow-rate weighting on surfaces with inconsistent orientation. Reworked to code to make better use of Enum (the NamedEnum replacement). Enum doesn't require contiguous enumeration values, which lets us use bitmasking of similar operations to reduce duplicate code. --- .../surfaceFieldValue/surfaceFieldValue.C | 129 ++++++++-------- .../surfaceFieldValue/surfaceFieldValue.H | 138 ++++++++++++----- .../surfaceFieldValue/surfaceFieldValueI.H | 17 +- .../surfaceFieldValueTemplates.C | 138 +++++++++-------- .../fieldValues/volFieldValue/volFieldValue.C | 71 ++++----- .../fieldValues/volFieldValue/volFieldValue.H | 71 ++++++--- .../volFieldValue/volFieldValueTemplates.C | 145 ++++++++---------- 7 files changed, 405 insertions(+), 304 deletions(-) diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C index f0eb481dd7..db92766bae 100644 --- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C +++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C @@ -69,23 +69,32 @@ const Foam::Enum > Foam::functionObjects::fieldValues::surfaceFieldValue::operationTypeNames_ { + // Normal operations { operationType::opNone, "none" }, + { operationType::opMin, "min" }, + { operationType::opMax, "max" }, { operationType::opSum, "sum" }, - { operationType::opWeightedSum, "weightedSum" }, { operationType::opSumMag, "sumMag" }, { operationType::opSumDirection, "sumDirection" }, { operationType::opSumDirectionBalance, "sumDirectionBalance" }, { operationType::opAverage, "average" }, - { operationType::opWeightedAverage, "weightedAverage" }, { operationType::opAreaAverage, "areaAverage" }, - { operationType::opWeightedAreaAverage, "weightedAreaAverage" }, { operationType::opAreaIntegrate, "areaIntegrate" }, - { operationType::opWeightedAreaIntegrate, "weightedAreaIntegrate" }, - { operationType::opMin, "min" }, - { operationType::opMax, "max" }, { operationType::opCoV, "CoV" }, { operationType::opAreaNormalAverage, "areaNormalAverage" }, { operationType::opAreaNormalIntegrate, "areaNormalIntegrate" }, + + // Using weighting + { operationType::opWeightedSum, "weightedSum" }, + { operationType::opWeightedAverage, "weightedAverage" }, + { operationType::opWeightedAreaAverage, "weightedAreaAverage" }, + { operationType::opWeightedAreaIntegrate, "weightedAreaIntegrate" }, + + // Using absolute weighting + { operationType::opAbsWeightedSum, "absWeightedSum" }, + { operationType::opAbsWeightedAverage, "absWeightedAverage" }, + { operationType::opAbsWeightedAreaAverage, "absWeightedAreaAverage" }, + { operationType::opAbsWeightedAreaIntegrate, "absWeightedAreaIntegrate" }, }; const Foam::Enum @@ -244,7 +253,7 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::combineMeshGeometry { if (facePatchId_[i] != -1) { - label patchi = facePatchId_[i]; + const label patchi = facePatchId_[i]; globalFacesIs[i] += mesh_.boundaryMesh()[patchi].start(); } } @@ -279,9 +288,8 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::combineMeshGeometry // My own data first { const faceList& fcs = allFaces[Pstream::myProcNo()]; - forAll(fcs, i) + for (const face& f : fcs) { - const face& f = fcs[i]; face& newF = faces[nFaces++]; newF.setSize(f.size()); forAll(f, fp) @@ -291,9 +299,9 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::combineMeshGeometry } const pointField& pts = allPoints[Pstream::myProcNo()]; - forAll(pts, i) + for (const point& pt : pts) { - points[nPoints++] = pts[i]; + points[nPoints++] = pt; } } @@ -303,9 +311,8 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::combineMeshGeometry if (proci != Pstream::myProcNo()) { const faceList& fcs = allFaces[proci]; - forAll(fcs, i) + for (const face& f : fcs) { - const face& f = fcs[i]; face& newF = faces[nFaces++]; newF.setSize(f.size()); forAll(f, fp) @@ -315,9 +322,9 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::combineMeshGeometry } const pointField& pts = allPoints[proci]; - forAll(pts, i) + for (const point& pt : pts) { - points[nPoints++] = pts[i]; + points[nPoints++] = pt; } } } @@ -343,9 +350,9 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::combineMeshGeometry } points.transfer(newPoints); - forAll(faces, i) + for (face& f : faces) { - inplaceRenumber(oldToNew, faces[i]); + inplaceRenumber(oldToNew, f); } } } @@ -447,17 +454,19 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::totalArea() const // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // -bool Foam::functionObjects::fieldValues::surfaceFieldValue::needsSf() const +bool Foam::functionObjects::fieldValues::surfaceFieldValue::usesSf() const { - // Many operations use the Sf field + // Only a few operations do not require the Sf field switch (operation_) { case opNone: + case opMin: + case opMax: case opSum: case opSumMag: case opAverage: - case opMin: - case opMax: + case opSumDirection: + case opSumDirectionBalance: { return false; } @@ -469,26 +478,6 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::needsSf() const } -bool Foam::functionObjects::fieldValues::surfaceFieldValue::needsWeight() const -{ - // Only a few operations use weight field - switch (operation_) - { - case opWeightedSum: - case opWeightedAverage: - case opWeightedAreaAverage: - case opWeightedAreaIntegrate: - { - return true; - } - default: - { - return false; - } - } -} - - void Foam::functionObjects::fieldValues::surfaceFieldValue::initialise ( const dictionary& dict @@ -582,19 +571,16 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::initialise weightFieldName_ = "none"; - if (needsWeight()) + if (usesWeight() && dict.readIfPresent("weightField", weightFieldName_)) { - if (dict.readIfPresent("weightField", weightFieldName_)) + if (regionType_ == stSampledSurface) { - if (regionType_ == stSampledSurface) - { - FatalIOErrorInFunction(dict) - << "Cannot use weightField for sampledSurface" - << exit(FatalIOError); - } - - Info<< " weight field = " << weightFieldName_ << nl; + FatalIOErrorInFunction(dict) + << "Cannot use weightField for sampledSurface" + << exit(FatalIOError); } + + Info<< " weight field = " << weightFieldName_ << nl; } // Backwards compatibility for v1612+ and older @@ -660,10 +646,10 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::writeFileHeader os << tab << "Area"; } - forAll(fields_, i) + for (const word& fieldName : fields_) { os << tab << operationTypeNames_[operation_] - << "(" << fields_[i] << ")"; + << "(" << fieldName << ")"; } os << endl; @@ -684,12 +670,12 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::processValues { case opSumDirection: { - vector n(dict_.lookup("direction")); + const vector n(dict_.lookup("direction")); return gSum(pos(values*(Sf & n))*mag(values)); } case opSumDirectionBalance: { - vector n(dict_.lookup("direction")); + const vector n(dict_.lookup("direction")); const scalarField nv(values*(Sf & n)); return gSum(pos(nv)*mag(values) - neg(nv)*mag(values)); @@ -754,10 +740,17 @@ Foam::tmp Foam::functionObjects::fieldValues::surfaceFieldValue::weightingFactor ( const Field& weightField -) +) const { - // pass through - return weightField; + if (usesMag()) + { + return mag(weightField); + } + else + { + // pass through + return weightField; + } } @@ -767,16 +760,21 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::weightingFactor ( const Field& weightField, const vectorField& Sf -) +) const { // scalar * Area if (returnReduce(weightField.empty(), andOp())) { + // No weight field - revert to unweighted form return mag(Sf); } + else if (usesMag()) + { + return mag(weightField * mag(Sf)); + } else { - return weightField * mag(Sf); + return (weightField * mag(Sf)); } } @@ -787,16 +785,21 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::weightingFactor ( const Field& weightField, const vectorField& Sf -) +) const { // vector (dot) Area if (returnReduce(weightField.empty(), andOp())) { + // No weight field - revert to unweighted form return mag(Sf); } + else if (usesMag()) + { + return mag(weightField & Sf); + } else { - return weightField & Sf; + return (weightField & Sf); } } @@ -915,7 +918,7 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::write() // Many operations use the Sf field vectorField Sf; - if (needsSf()) + if (usesSf()) { if (regionType_ == stSurface) { diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H index 04697bdb34..b34a052916 100644 --- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H +++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H @@ -103,22 +103,26 @@ Usage The \c operation is one of: \plaintable none | no operation - sum | sum - weightedSum | weighted sum - sumMag | sum of component magnitudes - sumDirection | sum values which are positive in given direction - sumDirectionBalance | sum of balance of values in given direction - average | ensemble average - weightedAverage | weighted average - areaAverage | area-weighted average - weightedAreaAverage | weighted area average - areaIntegrate | area integral - weightedAreaIntegrate | weighted area integral min | minimum max | maximum + sum | sum + sumMag | sum of component magnitudes + sumDirection | sum values that are positive in given direction + sumDirectionBalance | sum of balance of values in given direction + average | ensemble average + areaAverage | area-weighted average + areaIntegrate | area integral CoV | coefficient of variation: standard deviation/mean areaNormalAverage| area-weighted average in face normal direction areaNormalIntegrate | area-weighted integral in face normal directon + weightedSum | weighted sum + weightedAverage | weighted average + weightedAreaAverage | weighted area average + weightedAreaIntegrate | weighted area integral + absWeightedSum | sum using absolute weighting + absWeightedAverage | average using absolute weighting + absWeightedAreaAverage | area average using absolute weighting + absWeightedAreaIntegrate | area integral using absolute weighting \endplaintable Note @@ -190,36 +194,73 @@ public: //- Region type enumeration enum regionTypes { - stFaceZone, - stPatch, - stSurface, - stSampledSurface + stFaceZone, //!< Calculate on a faceZone + stPatch, //!< Calculate on a patch + stSurface, //!< Calculate with fields on a surfMesh + stSampledSurface //!< Sample onto surface and calculate }; //- Region type names static const Enum regionTypeNames_; + //- Bitmask values for operation variants + enum operationVariant + { + typeBase = 0, //!< Base operation + typeWeighted = 0x100, //!< Operation using weighting + typeAbsolute = 0x200, //!< Operation using mag (eg, for weighting) + }; //- Operation type enumeration enum operationType { - opNone, //!< None - opSum, //!< Sum - opWeightedSum, //!< Weighted sum - opSumMag, //!< Magnitude of sum + // Normal operations + + opNone = 0, //!< No operation + opMin, //!< Minimum value + opMax, //!< Maximum value + opSum, //!< Sum of values + opSumMag, //!< Sum of component magnitudes opSumDirection, //!< Sum in a given direction - opSumDirectionBalance, //!< Sum in a given direction for multiple - opAverage, //!< Average - opWeightedAverage, //!< Weighted average + + //! Sum of balance of values in given direction + opSumDirectionBalance, + + opAverage, //!< Ensemble average opAreaAverage, //!< Area average - opWeightedAreaAverage, //!< Weighted area average opAreaIntegrate, //!< Area integral - opWeightedAreaIntegrate, //!< Weighted area integral - opMin, //!< Minimum - opMax, //!< Maximum opCoV, //!< Coefficient of variation opAreaNormalAverage, //!< Area average in normal direction - opAreaNormalIntegrate //!< Area integral in normal direction + opAreaNormalIntegrate, //!< Area integral in normal direction + + // Weighted variants + + //! Weighted sum + opWeightedSum = (opSum | typeWeighted), + + //! Weighted average + opWeightedAverage = (opAverage | typeWeighted), + + //! Weighted area average + opWeightedAreaAverage = (opAreaAverage | typeWeighted), + + //! Weighted area integral + opWeightedAreaIntegrate = (opAreaIntegrate | typeWeighted), + + // Variants using absolute weighting + + //! Sum using abs weighting + opAbsWeightedSum = (opWeightedSum | typeAbsolute), + + //! Average using abs weighting + opAbsWeightedAverage = (opWeightedAverage | typeAbsolute), + + //! Area average using abs weighting + opAbsWeightedAreaAverage = (opWeightedAreaAverage | typeAbsolute), + + //! Area integral using abs weighting + opAbsWeightedAreaIntegrate = + (opWeightedAreaIntegrate | typeAbsolute), }; //- Operation type names @@ -229,8 +270,8 @@ public: //- Post-operation type enumeration enum postOperationType { - postOpNone, - postOpSqrt + postOpNone, //!< No additional operation after calculation + postOpSqrt //!< Perform sqrt after normal operation }; //- Operation type names @@ -325,11 +366,19 @@ protected: //- Return the local true/false list representing the face flip map inline const boolList& faceFlip() const; - //- True if the specified operation needs a surface Sf - bool needsSf() const; + //- True if the operation needs a surface Sf + bool usesSf() const; - //- True if the specified operation needs a weight-field - bool needsWeight() const; + //- True if the operation variant uses mag + inline bool usesMag() const; + + //- True if the operation variant uses a weight-field + inline bool usesWeight() const; + + //- True if operation variant uses a weight-field that is available. + // Checks for availability on any processor. + template + inline bool canWeight(const Field& weightField) const; //- Initialise, e.g. face addressing void initialise(const dictionary& dict); @@ -365,6 +414,7 @@ protected: const Field& weightField ) const; + //- Filter a surface field according to faceIds template tmp> filterField @@ -379,20 +429,24 @@ protected: const GeometricField& field ) const; - //- Weighting factor + + //- Weighting factor. + // Possibly applies mag() depending on the operation type. template - static tmp weightingFactor + tmp weightingFactor ( const Field& weightField - ); + ) const; - //- Weighting factor, weight field with the area + //- Weighting factor, weight field with the area. + // Possibly applies mag() depending on the operation type. + // Reverts to mag(Sf) if the weight field is not available. template - static tmp weightingFactor + tmp weightingFactor ( const Field& weightField, const vectorField& Sf - ); + ) const; //- Templated helper function to output field values @@ -489,7 +543,7 @@ template<> tmp surfaceFieldValue::weightingFactor ( const Field& weightField -); +) const; //- Specialisation for scalar - scalar * Area @@ -498,7 +552,7 @@ tmp surfaceFieldValue::weightingFactor ( const Field& weightField, const vectorField& Sf -); +) const; //- Specialisation for vector - vector (dot) Area @@ -507,7 +561,7 @@ tmp surfaceFieldValue::weightingFactor ( const Field& weightField, const vectorField& Sf -); +) const; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueI.H b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueI.H index e6dacee4df..8623229b6e 100644 --- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueI.H +++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueI.H @@ -23,7 +23,6 @@ License \*---------------------------------------------------------------------------*/ -#include "surfaceFieldValue.H" #include "Time.H" // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // @@ -51,6 +50,22 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::faceFlip() const // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +inline bool +Foam::functionObjects::fieldValues::surfaceFieldValue::usesMag() const +{ + // Operation specifically tagged to use mag + return (operation_ & typeAbsolute); +} + + +inline bool +Foam::functionObjects::fieldValues::surfaceFieldValue::usesWeight() const +{ + // Operation specifically tagged to require a weight field + return (operation_ & typeWeighted); +} + + inline const Foam::functionObjects::fieldValues::surfaceFieldValue::regionTypes& Foam::functionObjects::fieldValues::surfaceFieldValue::regionType() const { diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C index 806a068454..5852047488 100644 --- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C +++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C @@ -33,6 +33,20 @@ License // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // +template +inline bool Foam::functionObjects::fieldValues::surfaceFieldValue::canWeight +( + const Field& weightField +) const +{ + return + ( + usesWeight() + && returnReduce(!weightField.empty(), orOp()) // On some processor + ); +} + + template bool Foam::functionObjects::fieldValues::surfaceFieldValue::validField ( @@ -142,23 +156,14 @@ processSameTypeValues { break; } - case opSum: + case opMin: { - result = gSum(values); + result = gMin(values); break; } - case opWeightedSum: + case opMax: { - if (returnReduce(weightField.empty(), andOp())) - { - result = gSum(values); - } - else - { - tmp weight(weightingFactor(weightField)); - - result = gSum(weight*values); - } + result = gMax(values); break; } case opSumMag: @@ -166,6 +171,23 @@ processSameTypeValues result = gSum(cmptMag(values)); break; } + case opSum: + case opWeightedSum: + case opAbsWeightedSum: + { + if (canWeight(weightField)) + { + tmp weight(weightingFactor(weightField)); + + result = gSum(weight*values); + } + else + { + // Unweighted form + result = gSum(values); + } + break; + } case opSumDirection: case opSumDirectionBalance: { @@ -178,62 +200,56 @@ processSameTypeValues break; } case opAverage: - { - const label n = returnReduce(values.size(), sumOp