From 55af2fc2c65fefcf5c9dfb8488f498044a4ca76e Mon Sep 17 00:00:00 2001 From: Mark Olesen Date: Mon, 22 Nov 2021 15:59:46 +0100 Subject: [PATCH] BUG: missing unit-normal weighting for surfaceFieldValue (fixes #2273) - when using a vector field for weighting, it either used mag() or mag * area, but did not have a unit-normal projection version --- .../field/fieldValues/fieldValue/fieldValue.C | 2 +- .../field/fieldValues/fieldValue/fieldValue.H | 10 +- .../fieldValues/fieldValue/fieldValueI.H | 13 +- .../surfaceFieldValue/surfaceFieldValue.C | 149 ++++++++++++++---- .../surfaceFieldValue/surfaceFieldValue.H | 115 +++++++++----- .../surfaceFieldValue/surfaceFieldValueI.H | 25 +-- .../surfaceFieldValueTemplates.C | 77 +++++---- .../fieldValues/volFieldValue/volFieldValue.C | 34 +--- .../fieldValues/volFieldValue/volFieldValue.H | 17 +- .../volFieldValue/volFieldValueI.H | 56 +++++++ .../volFieldValue/volFieldValueTemplates.C | 8 +- 11 files changed, 343 insertions(+), 163 deletions(-) create mode 100644 src/functionObjects/field/fieldValues/volFieldValue/volFieldValueI.H diff --git a/src/functionObjects/field/fieldValues/fieldValue/fieldValue.C b/src/functionObjects/field/fieldValues/fieldValue/fieldValue.C index 21be1ff6cd..dc9746466e 100644 --- a/src/functionObjects/field/fieldValues/fieldValue/fieldValue.C +++ b/src/functionObjects/field/fieldValues/fieldValue/fieldValue.C @@ -114,7 +114,7 @@ bool Foam::functionObjects::fieldValue::execute() bool Foam::functionObjects::fieldValue::write() { - Log << type() << " " << name() << " write:" << nl; + Log << type() << ' ' << name() << " write:" << nl; return true; } diff --git a/src/functionObjects/field/fieldValues/fieldValue/fieldValue.H b/src/functionObjects/field/fieldValues/fieldValue/fieldValue.H index a182e30022..9964736369 100644 --- a/src/functionObjects/field/fieldValues/fieldValue/fieldValue.H +++ b/src/functionObjects/field/fieldValues/fieldValue/fieldValue.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2016-2020 OpenCFD Ltd. + Copyright (C) 2016-2021 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -186,16 +186,16 @@ public: // Member Functions //- Return the reference to the construction dictionary - inline const dictionary& dict() const; + inline const dictionary& dict() const noexcept; //- Return the region name - inline const word& regionName() const; + inline const word& regionName() const noexcept; //- Return the list of field names - inline const wordList& fields() const; + inline const wordList& fields() const noexcept; //- Return the output field values flag - inline bool writeFields() const; + inline bool writeFields() const noexcept; //- Read from dictionary virtual bool read(const dictionary& dict); diff --git a/src/functionObjects/field/fieldValues/fieldValue/fieldValueI.H b/src/functionObjects/field/fieldValues/fieldValue/fieldValueI.H index f667fe087f..2dcacd60a5 100644 --- a/src/functionObjects/field/fieldValues/fieldValue/fieldValueI.H +++ b/src/functionObjects/field/fieldValues/fieldValue/fieldValueI.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2020 OpenCFD Ltd. + Copyright (C) 2020-2021 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -28,25 +28,28 @@ License // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -inline const Foam::dictionary& Foam::functionObjects::fieldValue::dict() const +inline const Foam::dictionary& Foam::functionObjects::fieldValue::dict() +const noexcept { return dict_; } -inline const Foam::word& Foam::functionObjects::fieldValue::regionName() const +inline const Foam::word& Foam::functionObjects::fieldValue::regionName() +const noexcept { return regionName_; } -inline const Foam::wordList& Foam::functionObjects::fieldValue::fields() const +inline const Foam::wordList& Foam::functionObjects::fieldValue::fields() +const noexcept { return fields_; } -inline bool Foam::functionObjects::fieldValue::writeFields() const +inline bool Foam::functionObjects::fieldValue::writeFields() const noexcept { return writeFields_; } diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C index 3400e1b66b..92a5839027 100644 --- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C +++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C @@ -543,9 +543,10 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::totalArea() const // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // -bool Foam::functionObjects::fieldValues::surfaceFieldValue::usesSf() const +bool Foam::functionObjects::fieldValues::surfaceFieldValue::usesSf() +const noexcept { - // Only a few operations do not require the Sf field + // Few operations do not require the Sf field switch (operation_) { case opNone: @@ -554,13 +555,10 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::usesSf() const case opSum: case opSumMag: case opAverage: - { return false; - } + default: - { return true; - } } } @@ -703,11 +701,14 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::processValues scalar mean, numer; - if (canWeight(weightField)) + if (is_weightedOp() && canWeight(weightField)) { // Weighted quantity = (Weight * phi * dA) - tmp weight(weightingFactor(weightField)); + tmp weight + ( + weightingFactor(weightField, is_magOp()) + ); // Mean weighted value (area-averaged) mean = gSum(weight()*areaVal()) / areaTotal; @@ -786,11 +787,14 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::processValues scalar mean, numer; - if (canWeight(weightField)) + if (is_weightedOp() && canWeight(weightField)) { // Weighted quantity = (Weight * phi . dA) - tmp weight(weightingFactor(weightField)); + tmp weight + ( + weightingFactor(weightField, is_magOp()) + ); // Mean weighted value (area-averaged) mean = gSum(weight()*areaVal()) / areaTotal; @@ -824,14 +828,17 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::processValues } +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + template<> Foam::tmp Foam::functionObjects::fieldValues::surfaceFieldValue::weightingFactor ( - const Field& weightField -) const + const Field& weightField, + const bool useMag +) { - if (usesMag()) + if (useMag) { return mag(weightField); } @@ -846,17 +853,48 @@ Foam::tmp Foam::functionObjects::fieldValues::surfaceFieldValue::weightingFactor ( const Field& weightField, - const vectorField& Sf -) const + const vectorField& Sf, + const bool useMag +) +{ + // scalar * unit-normal + + // Can skip this check - already used canWeight() + /// if (returnReduce(weightField.empty(), andOp())) + /// { + /// // No weight field - revert to unweighted form? + /// return tmp::New(Sf.size(), scalar(1)); + /// } + + if (useMag) + { + return mag(weightField); + } + + // pass through + return weightField; +} + + +template<> +Foam::tmp +Foam::functionObjects::fieldValues::surfaceFieldValue::areaWeightingFactor +( + const Field& weightField, + const vectorField& Sf, + const bool useMag +) { // scalar * Area - if (returnReduce(weightField.empty(), andOp())) - { - // No weight field - revert to unweighted form - return mag(Sf); - } - else if (usesMag()) + // Can skip this check - already used canWeight() + /// if (returnReduce(weightField.empty(), andOp())) + /// { + /// // No weight field - revert to unweighted form + /// return mag(Sf); + /// } + + if (useMag) { return mag(weightField * mag(Sf)); } @@ -870,17 +908,61 @@ Foam::tmp Foam::functionObjects::fieldValues::surfaceFieldValue::weightingFactor ( const Field& weightField, - const vectorField& Sf -) const + const vectorField& Sf, + const bool useMag +) +{ + // vector (dot) unit-normal + + // Can skip this check - already used canWeight() + /// if (returnReduce(weightField.empty(), andOp())) + /// { + /// // No weight field - revert to unweighted form + /// return tmp::New(Sf.size(), scalar(1)); + /// } + + const label len = weightField.size(); + + auto tresult = tmp::New(weightField.size()); + auto& result = tresult.ref(); + + for (label facei=0; facei < len; ++facei) + { + const vector unitNormal(normalised(Sf[facei])); + result[facei] = (weightField[facei] & unitNormal); + } + + if (useMag) + { + for (scalar& val : result) + { + val = mag(val); + } + } + + return tresult; +} + + +template<> +Foam::tmp +Foam::functionObjects::fieldValues::surfaceFieldValue::areaWeightingFactor +( + const Field& weightField, + const vectorField& Sf, + const bool useMag +) { // vector (dot) Area - if (returnReduce(weightField.empty(), andOp())) - { - // No weight field - revert to unweighted form - return mag(Sf); - } - else if (usesMag()) + // Can skip this check - already used canWeight() + /// if (returnReduce(weightField.empty(), andOp())) + /// { + /// // No weight field - revert to unweighted form + /// return mag(Sf); + /// } + + if (useMag) { return mag(weightField & Sf); } @@ -959,6 +1041,13 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::surfaceFieldValue } +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +// Needs completed sampledSurface, surfaceWriter +Foam::functionObjects::fieldValues::surfaceFieldValue::~surfaceFieldValue() +{} + + // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // bool Foam::functionObjects::fieldValues::surfaceFieldValue::read @@ -1048,7 +1137,7 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::read Info<< operationTypeNames_[operation_] << nl; } - if (usesWeight()) + if (is_weightedOp()) { // Can have "weightFields" or "weightField" diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H index 71a4543db7..3782e878a7 100644 --- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H +++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H @@ -435,6 +435,40 @@ protected: autoPtr surfaceWriterPtr_; + // Static Member Functions + + //- Weighting factor. + // Possibly applies mag() depending on the operation type. + template + static tmp weightingFactor + ( + const Field& weightField, + const bool useMag + ); + + //- Weighting factor, weight field projected onto unit-normal. + // Possibly applies mag() depending on the operation type. + // Reverts to 'one' if the weight field is unavailable. + template + static tmp weightingFactor + ( + const Field& weightField, + const vectorField& Sf, + const bool useMag + ); + + //- Weighting factor, weight field with area factor. + // Possibly applies mag() depending on the operation type. + // Reverts to mag(Sf) if the weight field is unavailable. + template + static tmp areaWeightingFactor + ( + const Field& weightField, + const vectorField& Sf, + const bool useMag + ); + + // Protected Member Functions //- The volume mesh or surface registry being used @@ -444,30 +478,29 @@ protected: inline bool withSurfaceFields() const; //- Can use mesh topological merge? - inline bool withTopologicalMerge() const; + inline bool withTopologicalMerge() const noexcept; //- Return the local list of face IDs - inline const labelList& faceId() const; + inline const labelList& faceId() const noexcept; //- Return the local list of patch ID per face - inline const labelList& facePatch() const; + inline const labelList& facePatch() const noexcept; //- Return the local true/false list representing the face flip map - inline const boolList& faceFlip() const; + inline const boolList& faceFlip() const noexcept; //- True if the operation needs a surface Sf - bool usesSf() const; + bool usesSf() const noexcept; //- True if the operation variant uses mag - inline bool usesMag() const; + inline bool is_magOp() const noexcept; //- True if the operation variant uses a weight-field - inline bool usesWeight() const; + inline bool is_weightedOp() const noexcept; - //- True if operation variant uses a weight-field that is available. - // Checks for availability on any processor. + //- True if field is non-empty on any processor. template - inline bool canWeight(const Field& weightField) const; + inline bool canWeight(const Field& fld) const; //- Update the surface and surface information as required. // Do nothing (and return false) if no update was required @@ -519,26 +552,6 @@ protected: const GeometricField& field ) const; - - //- Weighting factor. - // Possibly applies mag() depending on the operation type. - template - tmp weightingFactor - ( - const Field& weightField - ) const; - - //- 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 - tmp weightingFactor - ( - const Field& weightField, - const vectorField& Sf - ) const; - - //- Templated helper function to output field values template label writeAll @@ -597,13 +610,13 @@ public: //- Destructor - virtual ~surfaceFieldValue() = default; + virtual ~surfaceFieldValue(); // Member Functions //- Return the region type - inline regionTypes regionType() const; + inline regionTypes regionType() const noexcept; //- Return the output directory inline fileName outputDir() const; @@ -646,26 +659,46 @@ vector surfaceFieldValue::processValues template<> tmp surfaceFieldValue::weightingFactor ( - const Field& weightField -) const; + const Field& weightField, + const bool useMag +); - -//- Specialisation for scalar - scalar * Area +//- Specialisation for scalar - pass through template<> tmp surfaceFieldValue::weightingFactor ( const Field& weightField, - const vectorField& Sf -) const; + const vectorField& Sf /* unused */, + const bool useMag +); + +//- Specialisation for scalar - scalar * Area +template<> +tmp surfaceFieldValue::areaWeightingFactor +( + const Field& weightField, + const vectorField& Sf, + const bool useMag +); -//- Specialisation for vector - vector (dot) Area +//- Specialisation for vector - vector (dot) unit-normal template<> tmp surfaceFieldValue::weightingFactor ( const Field& weightField, - const vectorField& Sf -) const; + const vectorField& Sf, + const bool useMag +); + +//- Specialisation for vector - vector (dot) Area +template<> +tmp surfaceFieldValue::areaWeightingFactor +( + const Field& weightField, + const vectorField& Sf, + const bool useMag +); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueI.H b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueI.H index f1ce22a48a..4e6a22d49e 100644 --- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueI.H +++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueI.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2017-2020 OpenCFD Ltd. + Copyright (C) 2017-2021 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -44,28 +44,31 @@ withSurfaceFields() const inline bool Foam::functionObjects::fieldValues::surfaceFieldValue:: -withTopologicalMerge() const +withTopologicalMerge() const noexcept { return (stFaceZone == regionType_ || stPatch == regionType_); } inline const Foam::labelList& -Foam::functionObjects::fieldValues::surfaceFieldValue::faceId() const +Foam::functionObjects::fieldValues::surfaceFieldValue::faceId() +const noexcept { return faceId_; } inline const Foam::labelList& -Foam::functionObjects::fieldValues::surfaceFieldValue::facePatch() const +Foam::functionObjects::fieldValues::surfaceFieldValue::facePatch() +const noexcept { return facePatchId_; } inline const Foam::boolList& -Foam::functionObjects::fieldValues::surfaceFieldValue::faceFlip() const +Foam::functionObjects::fieldValues::surfaceFieldValue::faceFlip() +const noexcept { return faceFlip_; } @@ -73,16 +76,17 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::faceFlip() const // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -inline bool -Foam::functionObjects::fieldValues::surfaceFieldValue::usesMag() const +inline bool Foam::functionObjects::fieldValues::surfaceFieldValue:: +is_magOp() +const noexcept { // Operation specifically tagged to use mag return (operation_ & typeAbsolute); } -inline bool -Foam::functionObjects::fieldValues::surfaceFieldValue::usesWeight() const +inline bool Foam::functionObjects::fieldValues::surfaceFieldValue:: +is_weightedOp() const noexcept { // Operation specifically tagged to require a weight field return (operation_ & typeWeighted); @@ -90,7 +94,8 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::usesWeight() const inline Foam::functionObjects::fieldValues::surfaceFieldValue::regionTypes -Foam::functionObjects::fieldValues::surfaceFieldValue::regionType() const +Foam::functionObjects::fieldValues::surfaceFieldValue::regionType() +const noexcept { return regionType_; } diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C index 9db9386e64..b4c18d8e64 100644 --- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C +++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C @@ -35,19 +35,32 @@ License #include "interpolationCell.H" #include "interpolationCellPoint.H" +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + +template +Foam::tmp +Foam::functionObjects::fieldValues::surfaceFieldValue::weightingFactor +( + const Field& weightField, + const bool useMag /* ignore */ +) +{ + // The scalar form is specialized. + // Other types: use mag() to generate a scalar field. + return mag(weightField); +} + + // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // template inline bool Foam::functionObjects::fieldValues::surfaceFieldValue::canWeight ( - const Field& weightField + const Field& fld ) const { - return - ( - usesWeight() - && returnReduce(!weightField.empty(), orOp()) // On some processor - ); + // Non-empty on some processor + return returnReduce(!fld.empty(), orOp()); } @@ -156,9 +169,13 @@ processSameTypeValues case opWeightedSum: case opAbsWeightedSum: { - if (canWeight(weightField)) + if (is_weightedOp() && canWeight(weightField)) { - tmp weight(weightingFactor(weightField)); + tmp weight + ( + weightingFactor(weightField, Sf, is_magOp()) + ); + result = gSum(weight*values); } else @@ -183,9 +200,12 @@ processSameTypeValues case opWeightedAverage: case opAbsWeightedAverage: { - if (canWeight(weightField)) + if (is_weightedOp() && canWeight(weightField)) { - const scalarField factor(weightingFactor(weightField)); + const scalarField factor + ( + weightingFactor(weightField, Sf, is_magOp()) + ); result = gSum(factor*values)/(gSum(factor) + ROOTVSMALL); } @@ -193,6 +213,7 @@ processSameTypeValues { // Unweighted form const label n = returnReduce(values.size(), sumOp