diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C index 36887ad927e793a2d522c044ffea9f9da863bbb1..4818000b3a0c973c8590a5bc93ed80e420dc9ce7 100644 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C @@ -798,6 +798,20 @@ bool Foam::AMIInterpolation::calculate } +bool Foam::AMIInterpolation::calculate +( + const polyMesh& mesh, + const label srcPatchi, + const primitivePatch& srcPatch, + const label tgtPatchi, + const primitivePatch& tgtPatch, + const autoPtr<searchableSurface>& surfPtr +) +{ + return calculate(srcPatch, tgtPatch, surfPtr); +} + + void Foam::AMIInterpolation::reset ( autoPtr<mapDistribute>&& srcToTgtMap, diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H index 5111a50ff2c521bcd4418cf6bc1e590ae84dec7c..549fe30869edf9205a4f8329c1b6d35feaece33e 100644 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H @@ -73,6 +73,8 @@ SourceFiles namespace Foam { +class polyMesh; + /*---------------------------------------------------------------------------*\ Class AMIInterpolation Declaration \*---------------------------------------------------------------------------*/ @@ -198,15 +200,6 @@ protected: ) const; - // Access - - //- Return the orginal src patch with optionally updated points - inline const primitivePatch& srcPatch0() const; - - //- Return the orginal tgt patch with optionally updated points - inline const primitivePatch& tgtPatch0() const; - - // Evaluation //- Normalise the (area) weights - suppresses numerical error in @@ -348,6 +341,15 @@ public: // Member Functions + // Access + + //- Return the orginal src patch with optionally updated points + inline const primitivePatch& srcPatch0() const; + + //- Return the orginal tgt patch with optionally updated points + inline const primitivePatch& tgtPatch0() const; + + // Access //- Access to the up-to-date flag @@ -446,6 +448,12 @@ public: //- patch weights (i.e. the sum before normalisation) inline const scalarField& tgtWeightsSum() const; + //- Return const access to target patch face centroids + inline const pointListList& tgtCentroids() const; + + //- Return access to target patch face centroids + inline pointListList& tgtCentroids(); + //- Return access to normalisation factor of target //- patch weights (i.e. the sum before normalisation) inline scalarField& tgtWeightsSum(); @@ -466,6 +474,17 @@ public: const autoPtr<searchableSurface>& surfPtr = nullptr ); + //- Update addressing, weights and (optional) centroids + virtual bool calculate + ( + const polyMesh& mesh, + const label srcPatchi, + const primitivePatch& srcPatch, + const label tgtPatchi, + const primitivePatch& tgtPatch, + const autoPtr<searchableSurface>& surfPtr = nullptr + ); + //- Set the maps, addresses and weights from an external source void reset ( diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H index cb9401d88593f44883e439aec1089e8f2d2e5cfc..d99323eb53de755851c27da1e611ab646787c036 100644 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H @@ -222,6 +222,18 @@ inline const Foam::scalarField& Foam::AMIInterpolation::tgtWeightsSum() const } +inline const Foam::pointListList& Foam::AMIInterpolation::tgtCentroids() const +{ + return tgtCentroids_; +} + + +inline Foam::pointListList& Foam::AMIInterpolation::tgtCentroids() +{ + return tgtCentroids_; +} + + inline Foam::scalarField& Foam::AMIInterpolation::tgtWeightsSum() { return tgtWeightsSum_; diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/lowWeightCorrection/edgeTopoDistancesData.C b/src/meshTools/AMIInterpolation/AMIInterpolation/lowWeightCorrection/edgeTopoDistancesData.C new file mode 100644 index 0000000000000000000000000000000000000000..8ad0dbb740e2e852b79c3ffcf8ca8d92ca650bb8 --- /dev/null +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/lowWeightCorrection/edgeTopoDistancesData.C @@ -0,0 +1,52 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2022 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // + +template<class Type, class PrimitivePatchType> +Foam::Ostream& Foam::operator<< +( + Foam::Ostream& os, + const Foam::edgeTopoDistancesData<Type, PrimitivePatchType>& wDist +) +{ + return os << wDist.distance_ << token::SPACE << wDist.data_; +} + + +template<class Type, class PrimitivePatchType> +Foam::Istream& Foam::operator>> +( + Foam::Istream& is, + Foam::edgeTopoDistancesData<Type, PrimitivePatchType>& wDist +) +{ + return is >> wDist.distance_ >> wDist.data_; +} + + +// ************************************************************************* // diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/lowWeightCorrection/edgeTopoDistancesData.H b/src/meshTools/AMIInterpolation/AMIInterpolation/lowWeightCorrection/edgeTopoDistancesData.H new file mode 100644 index 0000000000000000000000000000000000000000..119d671f69d72a2cb14b6c2f909b5a680f88117a --- /dev/null +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/lowWeightCorrection/edgeTopoDistancesData.H @@ -0,0 +1,261 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2022 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +Class + Foam::edgeTopoDistancesData + +Description + For use with PatchEdgeFaceWave. Determines topological distance to + starting edges. Templated on passive transported data. + +SourceFiles + edgeTopoDistancesDataI.H + edgeTopoDistancesData.C + +\*---------------------------------------------------------------------------*/ + +#ifndef edgeTopoDistancesData_H +#define edgeTopoDistancesData_H + +#include "point.H" +#include "tensor.H" +#include "indirectPrimitivePatch.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward Declarations +class polyPatch; +class polyMesh; +template<class Type, class PrimitivePatchType> +class edgeTopoDistancesData; + +template<class Type, class PrimitivePatchType> +Istream& operator>> +( + Istream&, + edgeTopoDistancesData<Type, PrimitivePatchType>& +); +template<class Type, class PrimitivePatchType> +Ostream& operator<< +( + Ostream&, + const edgeTopoDistancesData<Type, PrimitivePatchType>& +); + + +/*---------------------------------------------------------------------------*\ + Class edgeTopoDistancesData Declaration +\*---------------------------------------------------------------------------*/ + +template<class Type, class PrimitivePatchType = indirectPrimitivePatch> +class edgeTopoDistancesData +{ +public: + + //- Class used to pass additional data in + struct trackData + { + //- Max number of entries stored per element + label n_; + }; + + +protected: + + // Protected data + + //- Distance (up to size n_) + labelList distance_; + + //- Starting data (up to size n_) + List<Type> data_; + + +public: + + typedef Type dataType; + + // Constructors + + //- Construct null with null constructor for distance and data + inline edgeTopoDistancesData(); + + //- Construct from distance, data + inline edgeTopoDistancesData + ( + const labelList& distance, + const List<Type>& data + ); + + + // Member Functions + + // Access + + inline const labelList& distance() const + { + return distance_; + } + + inline const List<Type>& data() const + { + return data_; + } + + + // Needed by PatchEdgeFaceWave + + //- Check whether origin has been changed at all or + // still contains original (invalid) value. + template<class TrackingData> + inline bool valid(TrackingData& td) const; + + //- Apply rotation matrix + template<class TrackingData> + inline void transform + ( + const polyMesh& mesh, + const PrimitivePatchType& patch, + const tensor& rotTensor, + const scalar tol, + TrackingData& td + ); + + //- Influence of face on edge + template<class TrackingData> + inline bool updateEdge + ( + const polyMesh& mesh, + const PrimitivePatchType& patch, + const label edgeI, + const label facei, + const edgeTopoDistancesData<Type, PrimitivePatchType>& faceInfo, + const scalar tol, + TrackingData& td + ); + + //- New information for edge (from e.g. coupled edge) + template<class TrackingData> + inline bool updateEdge + ( + const polyMesh& mesh, + const PrimitivePatchType& patch, + const edgeTopoDistancesData<Type, PrimitivePatchType>& edgeInfo, + const bool sameOrientation, + const scalar tol, + TrackingData& td + ); + + //- Influence of edge on face. + template<class TrackingData> + inline bool updateFace + ( + const polyMesh& mesh, + const PrimitivePatchType& patch, + const label facei, + const label edgeI, + const edgeTopoDistancesData<Type, PrimitivePatchType>& edgeInfo, + const scalar tol, + TrackingData& td + ); + + //- Same (like operator==) + template<class TrackingData> + inline bool equal + ( + const edgeTopoDistancesData<Type, PrimitivePatchType>&, + TrackingData& + ) const; + + + // Member Operators + + // Needed for List IO + inline bool operator== + ( + const edgeTopoDistancesData<Type, PrimitivePatchType>& + ) const; + inline bool operator!= + ( + const edgeTopoDistancesData<Type, PrimitivePatchType>& + ) const; + + + // IOstream Operators + + friend Ostream& operator<< <Type, PrimitivePatchType> + ( + Ostream&, + const edgeTopoDistancesData<Type, PrimitivePatchType>& + ); + friend Istream& operator>> <Type, PrimitivePatchType> + ( + Istream&, + edgeTopoDistancesData<Type, PrimitivePatchType>& + ); +}; + + +// * * * * * * * * * * * * * * * * * Traits * * * * * * * * * * * * * * * * // + + +//- Data are contiguous if data type is contiguous +template<class Type, class PrimitivePatchType> +struct is_contiguous<edgeTopoDistancesData<Type, PrimitivePatchType>> : + std::false_type {}; + +//- Data are contiguous label if data type is label +template<class Type, class PrimitivePatchType> +struct is_contiguous_label<edgeTopoDistancesData<Type, PrimitivePatchType>> : + std::false_type {}; + +//- Data are contiguous scalar if data type is scalar +template<class Type, class PrimitivePatchType> +struct is_contiguous_scalar<edgeTopoDistancesData<Type, PrimitivePatchType>> : + std::false_type {}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository + #include "edgeTopoDistancesData.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "edgeTopoDistancesDataI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/lowWeightCorrection/edgeTopoDistancesDataI.H b/src/meshTools/AMIInterpolation/AMIInterpolation/lowWeightCorrection/edgeTopoDistancesDataI.H new file mode 100644 index 0000000000000000000000000000000000000000..e6fc2e8b30d5dbb38c37886069870d734b687825 --- /dev/null +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/lowWeightCorrection/edgeTopoDistancesDataI.H @@ -0,0 +1,224 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2022 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "polyMesh.H" +#include "transform.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template<class Type, class PrimitivePatchType> +inline +Foam::edgeTopoDistancesData<Type, PrimitivePatchType>::edgeTopoDistancesData() +: + distance_(), + data_() +{} + + +template<class Type, class PrimitivePatchType> +inline +Foam::edgeTopoDistancesData<Type, PrimitivePatchType>::edgeTopoDistancesData +( + const labelList& distance, + const List<Type>& data +) +: + distance_(distance), + data_(data) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template<class Type, class PrimitivePatchType> +template<class TrackingData> +inline bool Foam::edgeTopoDistancesData<Type, PrimitivePatchType>::valid +( + TrackingData& td +) const +{ + return distance_.size(); +} + + +template<class Type, class PrimitivePatchType> +template<class TrackingData> +inline void Foam::edgeTopoDistancesData<Type, PrimitivePatchType>::transform +( + const polyMesh& mesh, + const PrimitivePatchType& patch, + const tensor& rotTensor, + const scalar tol, + TrackingData& td +) +{} + + +template<class Type, class PrimitivePatchType> +template<class TrackingData> +inline bool Foam::edgeTopoDistancesData<Type, PrimitivePatchType>::updateEdge +( + const polyMesh& mesh, + const PrimitivePatchType& patch, + const label edgeI, + const label facei, + const edgeTopoDistancesData<Type, PrimitivePatchType>& faceInfo, + const scalar tol, + TrackingData& td +) +{ + // From face to edge + + if (distance_.size() >= td.n_) + { + return false; + } + + const label oldSize = data_.size(); + + forAll(faceInfo.data_, i) + { + const auto& d = faceInfo.data_[i]; + + if (!data_.found(d)) + { + data_.append(d); + distance_.append(faceInfo.distance_[i] + 1); + } + } + return data_.size() > oldSize; +} + + +template<class Type, class PrimitivePatchType> +template<class TrackingData> +inline bool Foam::edgeTopoDistancesData<Type, PrimitivePatchType>::updateEdge +( + const polyMesh& mesh, + const PrimitivePatchType& patch, + const edgeTopoDistancesData<Type, PrimitivePatchType>& edgeInfo, + const bool sameOrientation, + const scalar tol, + TrackingData& td +) +{ + // From edge to edge (e.g. coupled edges) + + if (distance_.size() >= td.n_) + { + return false; + } + + const label oldSize = data_.size(); + + forAll(edgeInfo.data_, i) + { + const auto& d = edgeInfo.data_[i]; + + if (!data_.found(d)) + { + data_.append(d); + distance_.append(edgeInfo.distance_[i]); + } + } + + return data_.size() > oldSize; +} + + +template<class Type, class PrimitivePatchType> +template<class TrackingData> +inline bool Foam::edgeTopoDistancesData<Type, PrimitivePatchType>::updateFace +( + const polyMesh& mesh, + const PrimitivePatchType& patch, + const label facei, + const label edgeI, + const edgeTopoDistancesData<Type, PrimitivePatchType>& edgeInfo, + const scalar tol, + TrackingData& td +) +{ + // From edge to edge (e.g. coupled edges) + + if (distance_.size() >= td.n_) + { + return false; + } + + const label oldSize = data_.size(); + + forAll(edgeInfo.data_, i) + { + const auto& d = edgeInfo.data_[i]; + + if (!data_.found(d)) + { + data_.append(d); + distance_.append(edgeInfo.distance_[i]); + } + } + + return data_.size() > oldSize; +} + + +template<class Type, class PrimitivePatchType> +template<class TrackingData> +inline bool Foam::edgeTopoDistancesData<Type, PrimitivePatchType>::equal +( + const edgeTopoDistancesData<Type, PrimitivePatchType>& rhs, + TrackingData& td +) const +{ + return operator==(rhs); +} + + +// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // + +template<class Type, class PrimitivePatchType> +inline bool Foam::edgeTopoDistancesData<Type, PrimitivePatchType>::operator== +( + const Foam::edgeTopoDistancesData<Type, PrimitivePatchType>& rhs +) const +{ + return distance() == rhs.distance() && data() == rhs.data(); +} + + +template<class Type, class PrimitivePatchType> +inline bool Foam::edgeTopoDistancesData<Type, PrimitivePatchType>::operator!= +( + const Foam::edgeTopoDistancesData<Type, PrimitivePatchType>& rhs +) const +{ + return !(*this == rhs); +} + + +// ************************************************************************* // diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/lowWeightCorrection/lowWeightCorrection.C b/src/meshTools/AMIInterpolation/AMIInterpolation/lowWeightCorrection/lowWeightCorrection.C new file mode 100644 index 0000000000000000000000000000000000000000..697e45e6eeca3fb88c48a8cad11277c3b5e8a134 --- /dev/null +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/lowWeightCorrection/lowWeightCorrection.C @@ -0,0 +1,578 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2022 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "lowWeightCorrection.H" +#include "addToRunTimeSelectionTable.H" +#include "edgeTopoDistancesData.H" +#include "PatchEdgeFaceWave.H" +#include "OBJstream.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(lowWeightCorrection, 0); + addToRunTimeSelectionTable(AMIInterpolation, lowWeightCorrection, dict); + addToRunTimeSelectionTable + ( + AMIInterpolation, + lowWeightCorrection, + component + ); +} + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::lowWeightCorrection::findNearest +( + const polyMesh& mesh, + const primitivePatch& pp, + const globalIndex& globalFaces, + const labelList& uncoveredFaces, + labelListList& nearestCoveredFaces +) const +{ + // Data on all edges and faces + typedef edgeTopoDistancesData<label, primitivePatch> Type; + + List<Type> allEdgeInfo(pp.nEdges()); + List<Type> allFaceInfo(pp.size()); + + + // Seed all covered faces + bitSet isUncovered(pp.size(), uncoveredFaces); + + DynamicList<label> changedEdges(pp.nEdges()); + DynamicList<Type> changedInfo(pp.nEdges()); + forAll(pp, facei) + { + if (!isUncovered(facei)) + { + const labelList& fEdges = pp.faceEdges()[facei]; + + const Type seed + ( + labelList(1, 0), + labelList(1, globalFaces.toGlobal(facei)) + ); + + allFaceInfo[facei] = seed; + for (const label edgei : fEdges) + { + changedEdges.append(edgei); + changedInfo.append(seed); + } + } + } + + + // Walk + Type::trackData td; + td.n_ = nDonors_; + + PatchEdgeFaceWave + < + primitivePatch, + Type, + Type::trackData + > calc + ( + mesh, + pp, + changedEdges, + changedInfo, + allEdgeInfo, + allFaceInfo, + 0, //returnReduce(pp.nEdges(), sumOp<label>()) + td + ); + calc.iterate(nIters_); // should be enough iterations? + + + nearestCoveredFaces.setSize(uncoveredFaces.size()); + forAll(uncoveredFaces, i) + { + const label facei = uncoveredFaces[i]; + if (allFaceInfo[facei].valid(calc.data())) + { + // Collect donor faces + nearestCoveredFaces[i] = allFaceInfo[facei].data(); + } + } +} + + +void Foam::lowWeightCorrection::calculateStencil +( + const label facei, + const point& sample, + const labelList& covered, + const labelListList& stencilAddressing, + + labelListList& addressing, + scalarListList& weights, + scalarField& sumWeights, + const pointField& faceCentres +) const +{ + addressing[facei].clear(); + weights[facei].clear(); + sumWeights[facei] = Zero; + + // The (uncovered) facei has local donors + + // Add the donors + for (const label coveredSlot : covered) + { + // Get remote slots for the local covered face + const labelList& slots = stencilAddressing[coveredSlot]; + + for (const label sloti : slots) + { + const point& donorPt = faceCentres[sloti]; + const label index = addressing[facei].find(sloti); + if (index == -1) + { + addressing[facei].append(sloti); + weights[facei].append(1.0/mag(sample-donorPt)); + } + else + { + weights[facei][index] += 1.0/mag(sample-donorPt); + } + } + } + + scalarList& wghts = weights[facei]; + const scalar w = sum(wghts); + forAll(wghts, i) + { + wghts[i] /= w; + } + sumWeights[facei] = 1.0; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::lowWeightCorrection::lowWeightCorrection +( + const dictionary& dict, + const bool reverseTarget +) +: + AMIInterpolation(dict, reverseTarget), + lowWeightCorrection_(dict.get<scalar>("lowWeightCorrection")), + nDonors_(dict.get<label>("nDonors")), + nIters_(dict.get<label>("nIters")), + AMIPtr_ + ( + AMIInterpolation::New + ( + dict.subDict(type() + "Coeffs").get<word>("AMIMethod"), + dict.subDict(type() + "Coeffs"), + reverseTarget + ) + ) +{ + if (nDonors_ <= 0 || nIters_ <= 0) + { + WarningInFunction << "Disabled low-weight correction. Using " + << AMIPtr_->type() << " AMI interpolation" << endl; + } +} + + +Foam::lowWeightCorrection::lowWeightCorrection +( + const bool requireMatch, + const bool reverseTarget, + const scalar lowWeightCorrection +) +: + AMIInterpolation(requireMatch, reverseTarget, lowWeightCorrection), + lowWeightCorrection_(-1), + nDonors_(0), + nIters_(0), + AMIPtr_() +{} + + +Foam::lowWeightCorrection::lowWeightCorrection(const lowWeightCorrection& ami) +: + AMIInterpolation(ami), + lowWeightCorrection_(ami.lowWeightCorrection_), + nDonors_(ami.nDonors_), + nIters_(ami.nIters_), + AMIPtr_(ami.AMIPtr_.clone()) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +bool Foam::lowWeightCorrection::calculate +( + const polyMesh& mesh, + const label srcPatchi, + const primitivePatch& srcPatch, + const label tgtPatchi, + const primitivePatch& tgtPatch, + const autoPtr<searchableSurface>& surfPtr +) +{ + if (upToDate_) + { + return false; + } + + + AMIInterpolation::calculate(srcPatch, tgtPatch, surfPtr); + AMIPtr_->calculate + ( + mesh, + srcPatchi, + srcPatch, + tgtPatchi, + tgtPatch, + surfPtr + ); + + // Take over AMI data + // ~~~~~~~~~~~~~~~~~~ + + upToDate_ = AMIPtr_->upToDate(); + // distributed / singlePatchProc + singlePatchProc_ = AMIPtr_->singlePatchProc();; + + // Source patch + srcMagSf_ = AMIPtr_->srcMagSf(); + srcAddress_ = AMIPtr_->srcAddress(); + srcWeights_ = AMIPtr_->srcWeights(); + srcWeightsSum_ = AMIPtr_->srcWeightsSum(); + srcCentroids_ = AMIPtr_->srcCentroids(); + //TBD: srcPatchPts_ = AMIPtr_->srcPatchPts(); + tsrcPatch0_ = AMIPtr_->srcPatch0(); + if (AMIPtr_->distributed()) + { + srcMapPtr_.reset(new mapDistribute(AMIPtr_->srcMap())); + } + + // Target + tgtMagSf_ = AMIPtr_->tgtMagSf(); + tgtAddress_ = AMIPtr_->tgtAddress(); + tgtWeights_ = AMIPtr_->tgtWeights(); + tgtWeightsSum_ = AMIPtr_->tgtWeightsSum(); + tgtCentroids_ = AMIPtr_->tgtCentroids(); + //TBD: tgtPatchPts_ = AMIPtr_->tgtPatchPts(); + ttgtPatch0_ = AMIPtr_->tgtPatch0(); + if (AMIPtr_->distributed()) + { + tgtMapPtr_.reset(new mapDistribute(AMIPtr_->tgtMap())); + } + + + // Walk out valid donors + // ~~~~~~~~~~~~~~~~~~~~~ + + if (nDonors_ > 0 && nIters_ > 0) + { + // Extend source side + { + // Low weight faces + DynamicList<label> uncoveredFaces; + forAll(srcWeightsSum(), facei) + { + if (srcWeightsSum()[facei] < lowWeightCorrection_) + { + uncoveredFaces.append(facei); + } + } + uncoveredFaces.shrink(); + + + // Global indexing for src faces + const globalIndex globalFaces(srcPatch.size()); + + + // Find nearest face with high weight + labelListList nearestCoveredFaces; + findNearest + ( + mesh, + mesh.boundaryMesh()[srcPatchi], + globalFaces, + uncoveredFaces, + nearestCoveredFaces + ); + + + // Create map to get remote data into nearestCoveredFaces + // order + List<Map<label>> compactMap; + const mapDistribute uncoveredToPatchMap + ( + globalFaces, + nearestCoveredFaces, + compactMap + ); + // Get srcPatch faceCentres in stencil ordering + pointField stencilFcs(srcPatch.faceCentres()); + uncoveredToPatchMap.distribute(stencilFcs); + // Get addressing (to donors) over in stencil ordering + labelListList stencilAddressing(srcAddress()); + uncoveredToPatchMap.distribute(stencilAddressing); + //scalarListList stencilWeights(srcWeights()); + //uncoveredToPatchMap.distribute(stencilWeights); + + + // Target side patch centres + tmp<pointField> totherFcs; + if (AMIPtr_->distributed()) + { + totherFcs = new pointField(tgtPatch.faceCentres()); + AMIPtr_->tgtMap().distribute(totherFcs.ref()); + } + else + { + totherFcs = tgtPatch.faceCentres(); + } + const pointField& otherFcs = totherFcs(); + + + //forAll(uncoveredFaces, i) + //{ + // const label facei = uncoveredFaces[i]; + // const labelList& covered = nearestCoveredFaces[i]; + // Pout<< "Uncovered face:" << facei + // << " low weihgt:" << srcWeightsSum()[facei] + // << " at:" << srcPatch.faceCentres()[facei] + // << " has local donors:" << endl; + // for (const label coveredSlot : covered) + // { + // Pout<< " fc:" << stencilFcs[coveredSlot] << nl + // << " which has remotes:" + // << stencilAddressing[coveredSlot] << nl + // << " at:" + // << UIndirectList<point> + // ( + // otherFcs, + // stencilAddressing[coveredSlot] + // ) << nl + // //<< " with weights:" + // //<< stencilWeights[coveredSlot] << nl + // << endl; + // } + //} + + + // Re-do interpolation on uncovered faces + forAll(uncoveredFaces, i) + { + const label facei = uncoveredFaces[i]; + + calculateStencil + ( + facei, + srcPatch.faceCentres()[facei], + nearestCoveredFaces[i], + stencilAddressing, + + srcAddress_, + srcWeights_, + srcWeightsSum_, + otherFcs + ); + } + + //OBJstream os + //( + // mesh.time().path() + // /mesh.boundaryMesh()[srcPatchi].name()+"_src.obj" + //); + //forAll(uncoveredFaces, i) + //{ + // const label facei = uncoveredFaces[i]; + // const point& fc = srcPatch.faceCentres()[facei]; + // const labelList& slots = srcAddress_[facei]; + // for (const label sloti : slots) + // { + // os.write(linePointRef(fc, otherFcs[sloti])); + // } + //} + } + + + // Extend target side + { + // Low weight faces + DynamicList<label> uncoveredFaces; + forAll(tgtWeightsSum(), facei) + { + if (tgtWeightsSum()[facei] < lowWeightCorrection_) + { + uncoveredFaces.append(facei); + } + } + uncoveredFaces.shrink(); + + + // Global indexing for tgt faces + const globalIndex globalFaces(tgtPatch.size()); + + + // Find nearest face with high weight + labelListList nearestCoveredFaces; + findNearest + ( + mesh, + mesh.boundaryMesh()[tgtPatchi], + globalFaces, + uncoveredFaces, + nearestCoveredFaces + ); + + + // Create map to get remote data into nearestCoveredFaces + // order + List<Map<label>> compactMap; + const mapDistribute uncoveredToPatchMap + ( + globalFaces, + nearestCoveredFaces, + compactMap + ); + // Get tgtPatch faceCentres in stencil ordering + pointField stencilFcs(tgtPatch.faceCentres()); + uncoveredToPatchMap.distribute(stencilFcs); + // Get addressing (to donors) over in stencil ordering + labelListList stencilAddressing(tgtAddress()); + uncoveredToPatchMap.distribute(stencilAddressing); + scalarListList stencilWeights(tgtWeights()); + uncoveredToPatchMap.distribute(stencilWeights); + + + // Target side patch centres + tmp<pointField> totherFcs; + if (AMIPtr_->distributed()) + { + totherFcs = new pointField(srcPatch.faceCentres()); + AMIPtr_->srcMap().distribute(totherFcs.ref()); + } + else + { + totherFcs = srcPatch.faceCentres(); + } + const pointField& otherFcs = totherFcs(); + + + //forAll(uncoveredFaces, i) + //{ + // const label facei = uncoveredFaces[i]; + // const labelList& covered = nearestCoveredFaces[i]; + // Pout<< "Uncovered face:" << facei + // << " low weihgt:" << tgtWeights()[facei] + // << " at:" << tgtPatch.faceCentres()[facei] + // << " has local donors:" << endl; + // for (const label coveredSlot : covered) + // { + // Pout<< " fc:" << stencilFcs[coveredSlot] << nl + // << " which has remotes:" + // << stencilAddressing[coveredSlot] << nl + // << " at:" + // << UIndirectList<point> + // ( + // otherFcs, + // stencilAddressing[coveredSlot] + // ) << nl + // << " with weights:" + // << stencilWeights[coveredSlot] << nl + // << endl; + // } + //} + + // Re-do interpolation on uncovered faces + forAll(uncoveredFaces, i) + { + const label facei = uncoveredFaces[i]; + + calculateStencil + ( + facei, + tgtPatch.faceCentres()[facei], + nearestCoveredFaces[i], + stencilAddressing, + + tgtAddress_, + tgtWeights_, + tgtWeightsSum_, + otherFcs + ); + } + + //OBJstream os + //( + // mesh.time().path() + // /mesh.boundaryMesh()[srcPatchi].name()+"_tgt.obj" + //); + //forAll(uncoveredFaces, i) + //{ + // const label facei = uncoveredFaces[i]; + // const point& fc = tgtPatch.faceCentres()[facei]; + // const labelList& slots = tgtAddress_[facei]; + // for (const label sloti : slots) + // { + // os.write(linePointRef(fc, otherFcs[sloti])); + // } + //} + } + + + ////- Write face connectivity as OBJ file + //writeFaceConnectivity + //( + // srcPatch, + // tgtPatch, + // srcAddress_ + //); + } + + return upToDate_; +} + + +void Foam::lowWeightCorrection::write(Ostream& os) const +{ + AMIInterpolation::write(os); + os.writeEntry("nDonors", nDonors_); + os.writeEntry("nIters", nIters_); + os.beginBlock(word(this->type() + "Coeffs")); + AMIPtr_->write(os); + os.endBlock(); +} + + +// ************************************************************************* // diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/lowWeightCorrection/lowWeightCorrection.H b/src/meshTools/AMIInterpolation/AMIInterpolation/lowWeightCorrection/lowWeightCorrection.H new file mode 100644 index 0000000000000000000000000000000000000000..8572603de6e54cffb9b7ab5606929867c23c6d44 --- /dev/null +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/lowWeightCorrection/lowWeightCorrection.H @@ -0,0 +1,216 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2022 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +Class + Foam::lowWeightCorrection + +Description + Arbitrary Mesh Interface (AMI) method wrapper that extends low-weight + faces to use donors of local neighbouring (non-low weight) faces. + + It collects valid nearby faces and uses their stencil instead. + + Example of the boundary specification: + \verbatim + <patchName> + { + type cyclicAMI; + + AMIMethod lowWeightCorrection; + + // Low weight value; which faces to override + lowWeightCorrection 0.1; + // How far to walk out from valid faces. 0 means no change. + nIters 2; + // Limit amount of valid faces to store. + nDonors 2; + // Underlying AMI method to use. + lowWeightCorrectionCoeffs + { + AMIMethod faceAreaWeightAMI; + } + } + \endverbatim + + This version walks out the local faces and then collects all their + remote donors. Hence even for e.g. 2 local faces (nDonors) the stencil + might still include lots of remote faces. + + The number of iterations to walk out might be set to infinite if one + wants a value on all faces. + + +SourceFiles + lowWeightCorrection.C + +\*---------------------------------------------------------------------------*/ + +#ifndef lowWeightCorrection_H +#define lowWeightCorrection_H + +#include "AMIInterpolation.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class lowWeightCorrection Declaration +\*---------------------------------------------------------------------------*/ + +class lowWeightCorrection +: + public AMIInterpolation +{ +private: + + // Private Data + + //- Threshold weight below which interpolation switches to distance + // weighted + const scalar lowWeightCorrection_; + + //- Number of donors to weigh + const label nDonors_; + + //- Number of iterations to walk out to search for donors + const label nIters_; + + //- Uncorrected AMI + autoPtr<AMIInterpolation> AMIPtr_; + + + // Private Member Functions + + //- Find (local) donor faces for all uncovered faces + void findNearest + ( + const polyMesh& mesh, + const primitivePatch& pp, + const globalIndex& globalFaces, + const labelList& uncoveredFaces, + labelListList& nearestCoveredFaces + ) const; + + //- Determine stencil for single face. Gets list of (local) faces + // whose addressing can be used for calculating the stencil + void calculateStencil + ( + const label facei, + const point& sample, + const labelList& nearestCovered, + const labelListList& stencilAddressing, + + labelListList& addressing, + scalarListList& weights, + scalarField& sumWeights, + const pointField& faceCentres + ) const; + + + //- No copy assignment + void operator=(const lowWeightCorrection&) = delete; + + +public: + + //- Runtime type information + TypeName("lowWeightCorrection"); + + + // Constructors + + //- Construct from dictionary + lowWeightCorrection + ( + const dictionary& dict, + const bool reverseTarget = false + ); + + //- Construct from components + lowWeightCorrection + ( + const bool requireMatch = true, + const bool reverseTarget = false, + const scalar lowWeightCorrection = -1 + ); + + //- Construct as copy + lowWeightCorrection(const lowWeightCorrection& ami); + + //- Construct and return a clone + virtual autoPtr<AMIInterpolation> clone() const + { + return autoPtr<AMIInterpolation>(new lowWeightCorrection(*this)); + } + + + //- Destructor + virtual ~lowWeightCorrection() = default; + + + // Member Functions + + //- Update addressing and weights + virtual bool calculate + ( + const primitivePatch& srcPatch, + const primitivePatch& tgtPatch, + const autoPtr<searchableSurface>& surfPtr = nullptr + ) + { + NotImplemented; + return false; + } + + //- Update addressing and weights + virtual bool calculate + ( + const polyMesh& mesh, + const label srcPatchi, + const primitivePatch& srcPatch, + const label tgtPatchi, + const primitivePatch& tgtPatch, + const autoPtr<searchableSurface>& surfPtr = nullptr + ); + + + // I-O + + //- Write + virtual void write(Ostream& os) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C b/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C index b81b905323b20d32ac80cfd3c24365d0dbe08735..23fdb1f946528ed9b56fcba15d3f3dbe4e99bf72 100644 --- a/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C +++ b/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C @@ -435,7 +435,15 @@ void Foam::cyclicAMIPolyPatch::resetAMI(const UList<point>& points) const // Construct/apply AMI interpolation to determine addressing and weights AMIPtr_->upToDate() = false; - AMIPtr_->calculate(patch0, nbrPatch0, surfPtr()); + AMIPtr_->calculate + ( + boundaryMesh().mesh(), + index(), + patch0, + nbr.index(), + nbrPatch0, + surfPtr() + ); if (debug) { diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files index 4309290e20b48c54b0d295b64e24528e7419fd0b..da1fd19384dafc134a4e5c794524a3b2f6538dc9 100644 --- a/src/meshTools/Make/files +++ b/src/meshTools/Make/files @@ -293,6 +293,7 @@ $(AMI)/GAMG/interfaceFields/cyclicACMIGAMGInterfaceField/cyclicACMIGAMGInterface $(AMI)/triangle2D/triangle2D.C $(AMI)/AMIInterpolation/faceAreaWeightAMI2D/faceAreaWeightAMI2D.C +$(AMI)/AMIInterpolation/lowWeightCorrection/lowWeightCorrection.C AMICycPatches=$(AMI)/patches/cyclicAMI