diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.C index f4265f155ab08f7619c1222e28fd110761953c2f..5b3d1bc483bdc3fd8ce07a02e4fc52377822dee8 100644 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.C +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.C @@ -162,7 +162,7 @@ bool Foam::AMIMethod<SourcePatch, TargetPatch>::initialise } // Reset the octree - resetTree(); + treePtr_.reset(createTree<TargetPatch>(tgtPatch())); // Find initial face match using brute force/octree search if ((srcFacei == -1) || (tgtFacei == -1)) @@ -256,35 +256,31 @@ void Foam::AMIMethod<SourcePatch, TargetPatch>::writeIntersectionOBJ template<class SourcePatch, class TargetPatch> -void Foam::AMIMethod<SourcePatch, TargetPatch>::resetTree() +template<class PatchType> +Foam::autoPtr<Foam::indexedOctree<Foam::treeDataPrimitivePatch<PatchType>>> +Foam::AMIMethod<SourcePatch, TargetPatch>::createTree +( + const PatchType& patch +) const { - const auto& tgt = tgtPatch(); - - // Clear the old octree - treePtr_.clear(); + typedef treeDataPrimitivePatch<PatchType> PatchTreeType; - treeBoundBox bb(tgt.points(), tgt.meshPoints()); + treeBoundBox bb(patch.points(), patch.meshPoints()); bb.inflate(0.01); - if (!treePtr_.valid()) - { - treePtr_.reset + return autoPtr<indexedOctree<PatchTreeType>>::New + ( + PatchTreeType ( - new indexedOctree<treeType> - ( - treeType - ( - false, - tgt, - indexedOctree<treeType>::perturbTol() - ), - bb, // overall search domain - 8, // maxLevel - 10, // leaf size - 3.0 // duplicity - ) - ); - } + false, + patch, + indexedOctree<PatchTreeType>::perturbTol() + ), + bb, // overall search domain + 8, // maxLevel + 10, // leaf size + 3.0 // duplicity + ); } diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.H index e702f8366ab38747eef861057f593da17257c86f..62efb5d48684080300a7c8d7ce7315740201dec0 100644 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.H +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.H @@ -199,7 +199,12 @@ protected: // Common AMI method functions //- Reset the octree for the target patch face search - void resetTree(); + template<class PatchType> + autoPtr<indexedOctree<treeDataPrimitivePatch<PatchType>>> + createTree + ( + const PatchType& patch + ) const; label findTargetFace ( diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.C deleted file mode 100644 index e2d57b8599d2e591aa84f4fffe0660f5cb5a7cc3..0000000000000000000000000000000000000000 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.C +++ /dev/null @@ -1,444 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | www.openfoam.com - \\/ M anipulation | -------------------------------------------------------------------------------- - Copyright (C) 2013-2016 OpenFOAM Foundation - Copyright (C) 2016 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 "mapNearestAMI.H" - -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - -template<class SourcePatch, class TargetPatch> -void Foam::mapNearestAMI<SourcePatch, TargetPatch>::findNearestFace -( - const SourcePatch& srcPatch, - const TargetPatch& tgtPatch, - const label& srcFacei, - label& tgtFacei -) const -{ - const vectorField& srcCf = srcPatch.faceCentres(); - const vectorField& tgtCf = tgtPatch.faceCentres(); - - const vector srcP = srcCf[srcFacei]; - - DynamicList<label> tgtFaces(10); - tgtFaces.append(tgtFacei); - - DynamicList<label> visitedFaces(10); - - scalar d = GREAT; - - do - { - label tgtI = tgtFaces.remove(); - visitedFaces.append(tgtI); - - scalar dTest = magSqr(tgtCf[tgtI] - srcP); - if (dTest < d) - { - tgtFacei = tgtI; - d = dTest; - - this->appendNbrFaces - ( - tgtFacei, - tgtPatch, - visitedFaces, - tgtFaces - ); - } - - } while (tgtFaces.size() > 0); -} - - -template<class SourcePatch, class TargetPatch> -void Foam::mapNearestAMI<SourcePatch, TargetPatch>::setNextNearestFaces -( - boolList& mapFlag, - label& startSeedI, - label& srcFacei, - label& tgtFacei -) const -{ - const labelList& srcNbr = this->srcPatch().faceFaces()[srcFacei]; - - srcFacei = -1; - - for (const label facei : srcNbr) - { - if (mapFlag[facei]) - { - srcFacei = facei; - startSeedI = facei + 1; - - return; - } - } - - forAll(mapFlag, facei) - { - if (mapFlag[facei]) - { - srcFacei = facei; - tgtFacei = this->findTargetFace(facei); - - if (tgtFacei == -1) - { - const vectorField& srcCf = this->srcPatch().faceCentres(); - - FatalErrorInFunction - << "Unable to find target face for source face " - << srcFacei << " with face centre " << srcCf[srcFacei] - << abort(FatalError); - } - - break; - } - } -} - - -template<class SourcePatch, class TargetPatch> -Foam::label Foam::mapNearestAMI<SourcePatch, TargetPatch>::findMappedSrcFace -( - const label tgtFacei, - const List<DynamicList<label>>& tgtToSrc -) const -{ - DynamicList<label> testFaces(10); - DynamicList<label> visitedFaces(10); - - testFaces.append(tgtFacei); - - do - { - // search target tgtFacei neighbours for match with source face - label tgtI = testFaces.remove(); - - if (!visitedFaces.found(tgtI)) - { - visitedFaces.append(tgtI); - - if (tgtToSrc[tgtI].size()) - { - return tgtToSrc[tgtI][0]; - } - else - { - const labelList& nbrFaces = this->tgtPatch().faceFaces()[tgtI]; - - for (const label nbrFacei : nbrFaces) - { - if (!visitedFaces.found(nbrFacei)) - { - testFaces.append(nbrFacei); - } - } - } - } - } while (testFaces.size()); - - // did not find any match - should not be possible to get here! - return -1; -} - - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -template<class SourcePatch, class TargetPatch> -Foam::mapNearestAMI<SourcePatch, TargetPatch>::mapNearestAMI -( - const SourcePatch& srcPatch, - const TargetPatch& tgtPatch, - const faceAreaIntersect::triangulationMode& triMode, - const bool reverseTarget, - const bool requireMatch -) -: - AMIMethod<SourcePatch, TargetPatch> - ( - srcPatch, - tgtPatch, - triMode, - reverseTarget, - requireMatch - ) -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -template<class SourcePatch, class TargetPatch> -bool Foam::mapNearestAMI<SourcePatch, TargetPatch>::calculate -( - labelListList& srcAddress, - scalarListList& srcWeights, - pointListList& srcCentroids, - labelListList& tgtAddress, - scalarListList& tgtWeights, - scalarList& srcMagSf, - scalarList& tgtMagSf, - autoPtr<mapDistribute>& srcMapPtr, - autoPtr<mapDistribute>& tgtMapPtr, - label srcFacei, - label tgtFacei -) -{ - bool ok = - this->initialise - ( - srcAddress, - srcWeights, - tgtAddress, - tgtWeights, - srcFacei, - tgtFacei - ); - - if (!ok) - { - return false; - } - - - // temporary storage for addressing and weights - List<DynamicList<label>> srcAddr(this->srcPatch().size()); - List<DynamicList<label>> tgtAddr(this->tgtPatch().size()); - - - // construct weights and addressing - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - // list to keep track of whether src face can be mapped - boolList mapFlag(srcAddr.size(), true); - - // reset starting seed - label startSeedI = 0; - - DynamicList<label> nonOverlapFaces; - do - { - findNearestFace(this->srcPatch(), this->tgtPatch(), srcFacei, tgtFacei); - - srcAddr[srcFacei].append(tgtFacei); - tgtAddr[tgtFacei].append(srcFacei); - - mapFlag[srcFacei] = false; - - // Do advancing front starting from srcFacei, tgtFacei - setNextNearestFaces - ( - mapFlag, - startSeedI, - srcFacei, - tgtFacei - ); - } while (srcFacei >= 0); - - - // for the case of multiple source faces per target face, select the - // nearest source face only and discard the others - const vectorField& srcCf = this->srcPatch().faceCentres(); - const vectorField& tgtCf = this->tgtPatch().faceCentres(); - - forAll(tgtAddr, targetFacei) - { - if (tgtAddr[targetFacei].size() > 1) - { - const vector& tgtC = tgtCf[tgtFacei]; - - DynamicList<label>& srcFaces = tgtAddr[targetFacei]; - - label srcFacei = srcFaces[0]; - scalar d = magSqr(tgtC - srcCf[srcFacei]); - - for (label i = 1; i < srcFaces.size(); ++i) - { - label srcI = srcFaces[i]; - scalar dNew = magSqr(tgtC - srcCf[srcI]); - if (dNew < d) - { - d = dNew; - srcFacei = srcI; - } - } - - srcFaces.clear(); - srcFaces.append(srcFacei); - } - } - - // If there are more target faces than source faces, some target faces - // might not yet be mapped - forAll(tgtAddr, tgtFacei) - { - if (tgtAddr[tgtFacei].empty()) - { - label srcFacei = findMappedSrcFace(tgtFacei, tgtAddr); - - if (srcFacei >= 0) - { - // note - reversed search from src->tgt to tgt->src - findNearestFace - ( - this->tgtPatch(), - this->srcPatch(), - tgtFacei, - srcFacei - ); - - tgtAddr[tgtFacei].append(srcFacei); - } - } - } - - - // transfer data to persistent storage - const pointField& srcFc = this->srcPatch().faceCentres(); - const pointField& tgtFc = this->tgtPatch().faceCentres(); - - forAll(srcAddr, srcI) - { - srcAddress[srcI].transfer(srcAddr[srcI]); - - const labelList& addr = srcAddress[srcI]; - srcWeights[srcI].setSize(addr.size()); - const point& srcPt = srcFc[srcI]; - forAll(addr, i) - { - srcWeights[srcI][i] = magSqr(srcPt-tgtFc[addr[i]]); - } - } - forAll(tgtAddr, tgtI) - { - tgtAddress[tgtI].transfer(tgtAddr[tgtI]); - - const labelList& addr = tgtAddress[tgtI]; - tgtWeights[tgtI].setSize(addr.size()); - const point& tgtPt = tgtFc[tgtI]; - forAll(addr, i) - { - tgtWeights[tgtI][i] = magSqr(tgtPt-srcFc[addr[i]]); - } - } - - return true; -} - - -template<class SourcePatch, class TargetPatch> -void Foam::mapNearestAMI<SourcePatch, TargetPatch>::setMagSf -( - const TargetPatch& tgtPatch, - const mapDistribute& map, - scalarList& srcMagSf, - scalarList& tgtMagSf -) const -{ - srcMagSf = std::move(this->srcMagSf_); - tgtMagSf = scalarList(tgtPatch.size(), 1.0); -} - - -template<class SourcePatch, class TargetPatch> -void Foam::mapNearestAMI<SourcePatch, TargetPatch>::normaliseWeights -( - const bool verbose, - AMIInterpolation<SourcePatch, TargetPatch>& inter -) -{ - { - labelListList& srcAddress = inter.srcAddress(); - scalarListList& srcWeights = inter.srcWeights(); - - forAll(srcAddress, srcI) - { - labelList& addr = srcAddress[srcI]; - scalarList& wghts = srcWeights[srcI]; - - // Choose one with smallest weight (since calculate above returns - // distance) - if (addr.size()) - { - label minFaceI = addr[0]; - scalar minWeight = wghts[0]; - - for (label i = 0; i < addr.size(); ++i) - { - if (wghts[i] < minWeight) - { - minWeight = wghts[i]; - minFaceI = addr[i]; - } - } - - wghts.setSize(1); - wghts[0] = this->srcMagSf_[srcI]; - addr.setSize(1); - addr[0] = minFaceI; - } - } - } - - { - labelListList& tgtAddress = inter.tgtAddress(); - scalarListList& tgtWeights = inter.tgtWeights(); - - forAll(tgtAddress, tgtI) - { - labelList& addr = tgtAddress[tgtI]; - scalarList& wghts = tgtWeights[tgtI]; - - // Choose one with smallest weight (since calculate above returns - // distance) - if (addr.size()) - { - label minFaceI = addr[0]; - scalar minWeight = wghts[0]; - - forAll(addr, i) - { - if (wghts[i] < minWeight) - { - minWeight = wghts[i]; - minFaceI = addr[i]; - } - } - - wghts.setSize(1); - wghts[0] = inter.tgtMagSf()[tgtI]; - addr.setSize(1); - addr[0] = minFaceI; - } - } - } - - inter.normaliseWeights(this->conformal(), verbose); -} - - -// ************************************************************************* // diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/nearestFaceAMI.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/nearestFaceAMI.C new file mode 100644 index 0000000000000000000000000000000000000000..0aa96a35c4469d14ca93efcc4a1ec474b0e12bc2 --- /dev/null +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/nearestFaceAMI.C @@ -0,0 +1,361 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2020 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 "nearestFaceAMI.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template<class SourcePatch, class TargetPatch> +Foam::autoPtr<Foam::mapDistribute> +Foam::nearestFaceAMI<SourcePatch, TargetPatch>::calcFaceMap +( + const List<nearestAndDist>& localInfo, + const SourcePatch& srcPatch, + const TargetPatch& tgtPatch +) const +{ + // Generate the list of processor bounding boxes + List<boundBox> procBbs(Pstream::nProcs()); + procBbs[Pstream::myProcNo()] = + boundBox(srcPatch.points(), srcPatch.meshPoints(), true); + Pstream::gatherList(procBbs); + Pstream::scatterList(procBbs); + + // Identify which of my local tgt faces intersect with each processor bb + // within the current match's search distance + const pointField& tgtCcs = tgtPatch.faceCentres(); + List<DynamicList<label>> dynSendMap(Pstream::nProcs()); + + forAll(localInfo, tgtFacei) + { + const scalar r2 = localInfo[tgtFacei].second(); + + // Construct local bounding box to test against processor bb + forAll(procBbs, proci) + { + if (proci != Pstream::myProcNo()) + { + if (procBbs[proci].overlaps(tgtCcs[tgtFacei], r2)) + { + dynSendMap[proci].append(tgtFacei); + } + } + } + } + + // Convert dynamicList to labelList + labelListList sendMap(Pstream::nProcs()); + forAll(sendMap, proci) + { + dynSendMap[proci].shrink(); + sendMap[proci].transfer(dynSendMap[proci]); + + if (debug) + { + Pout<< "send map - to proc " << proci << " sending " + << sendMap[proci].size() << " elements" << endl; + } + } + + return autoPtr<mapDistribute>::New(std::move(sendMap)); +} + + +template<class SourcePatch, class TargetPatch> +Foam::autoPtr<Foam::mapDistribute> +Foam::nearestFaceAMI<SourcePatch, TargetPatch>::calcDistributed +( + const SourcePatch& src, + const TargetPatch& tgt, + labelListList& srcToTgtAddr, + scalarListList& srcToTgtWght +) const +{ + const auto tgtTreePtr = this->createTree(tgt); + const auto& tgtTree = tgtTreePtr(); + + // Create global indexing for each patch + globalIndex globalTgtCells(src.size()); + + // First pass: determine local matches + + // Identify local nearest matches + pointField srcCcs(src.faceCentres()); + + List<nearestAndDist> localInfo(src.size()); + forAll(srcCcs, srcCelli) + { + const point& srcCc = srcCcs[srcCelli]; + + pointIndexHit& test = localInfo[srcCelli].first(); + test = tgtTree.findNearest(srcCc, GREAT); + + if (test.hit()) + { + // With a search radius2 of GREAT all cells should receive a hit + localInfo[srcCelli].second() = magSqr(srcCc - test.hitPoint()); + test.setIndex(globalTgtCells.toGlobal(test.index())); + } + } + + // Second pass: determine remote matches + + autoPtr<mapDistribute> mapPtr = calcFaceMap(localInfo, src, tgt); + mapDistribute& map = mapPtr(); + + List<nearestAndDist> remoteInfo(localInfo); + map.distribute(remoteInfo); + + // Note: re-using srcCcs + map.distribute(srcCcs); + + // Test remote target cells against local source cells + nearestAndDist testInfo; + pointIndexHit& test = testInfo.first(); + forAll(remoteInfo, i) + { + test = tgtTree.findNearest(srcCcs[i], remoteInfo[i].second()); + if (test.hit()) + { + testInfo.first().setIndex + ( + globalTgtCells.toGlobal(test.index()) + ); + testInfo.second() = magSqr(test.hitPoint() - srcCcs[i]); + nearestEqOp()(remoteInfo[i], testInfo); + } + } + + // Send back to originating processor. Choose best if sent to multiple + // processors. Note that afterwards all unused entries have the unique + // value nearestZero (distance < 0). This is used later on to see if + // the sample was sent to any processor. + const nearestAndDist nearestZero(pointIndexHit(), -GREAT); + mapDistributeBase::distribute + ( + Pstream::commsTypes::nonBlocking, + List<labelPair>(0), + src.size(), + map.constructMap(), + map.constructHasFlip(), + map.subMap(), + map.subHasFlip(), + remoteInfo, + nearestEqOp(), + noOp(), // no flipping + nearestZero + ); + + // Third pass: combine local and remote info and filter out any + // connections that are further away than threshold distance squared + srcToTgtAddr.setSize(src.size()); + srcToTgtWght.setSize(src.size()); + forAll(srcToTgtAddr, srcCelli) + { + nearestEqOp()(localInfo[srcCelli], remoteInfo[srcCelli]); + if (localInfo[srcCelli].second() < maxDistance2_) + { + const label tgtCelli = localInfo[srcCelli].first().index(); + srcToTgtAddr[srcCelli] = labelList(1, tgtCelli); + srcToTgtWght[srcCelli] = scalarList(1, 1.0); + } + } + + List<Map<label>> cMap; + return autoPtr<mapDistribute>::New(globalTgtCells, srcToTgtAddr, cMap); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template<class SourcePatch, class TargetPatch> +Foam::nearestFaceAMI<SourcePatch, TargetPatch>::nearestFaceAMI +( + const SourcePatch& srcPatch, + const TargetPatch& tgtPatch, + const faceAreaIntersect::triangulationMode& triMode, + const bool reverseTarget, + const bool requireMatch +) +: + AMIMethod<SourcePatch, TargetPatch> + ( + srcPatch, + tgtPatch, + triMode, + reverseTarget, + requireMatch + ), + maxDistance2_(GREAT) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template<class SourcePatch, class TargetPatch> +bool Foam::nearestFaceAMI<SourcePatch, TargetPatch>::calculate +( + labelListList& srcAddress, + scalarListList& srcWeights, + pointListList& srcCentroids, + labelListList& tgtAddress, + scalarListList& tgtWeights, + scalarList& srcMagSf, + scalarList& tgtMagSf, + autoPtr<mapDistribute>& srcMapPtr, + autoPtr<mapDistribute>& tgtMapPtr, + label srcFacei, + label tgtFacei +) +{ +bool symmetric_ = true; + + if (this->distributed()) + { + tgtMapPtr = + calcDistributed + ( + this->srcPatch0_, + this->tgtPatch0_, + srcAddress, + srcWeights + ); + + if (symmetric_) + { + srcMapPtr = + calcDistributed + ( + this->tgtPatch0_, + this->srcPatch0_, + tgtAddress, + tgtWeights + ); + } + } + else + { + srcAddress.setSize(this->srcPatch0_.size()); + srcWeights.setSize(this->srcPatch0_.size()); + + if (symmetric_) + { + tgtAddress.setSize(this->tgtPatch0_.size()); + tgtWeights.setSize(this->tgtPatch0_.size()); + } + + const pointField& srcCcs = this->srcPatch0_.faceCentres(); + const pointField& tgtCcs = this->tgtPatch0_.faceCentres(); + + const auto tgtTreePtr = this->createTree(this->tgtPatch0_); + const auto& tgtTree = tgtTreePtr(); + + forAll(srcCcs, srcFacei) + { + const point& srcCc = srcCcs[srcFacei]; + const pointIndexHit hit = tgtTree.findNearest(srcCc, GREAT); + + if + ( + hit.hit() + && (magSqr(srcCc - tgtCcs[hit.index()]) < maxDistance2_) + ) + { + label tgtFacei = hit.index(); + srcAddress[srcFacei] = labelList(1, tgtFacei); + srcWeights[srcFacei] = scalarList(1, 1.0); + + if (symmetric_) + { + tgtAddress[tgtFacei] = labelList(1, srcFacei); + tgtWeights[tgtFacei] = scalarList(1, 1.0); + } + } + else + { + if (debug) + { + WarningInFunction + << "Unable to find target face for source face " + << srcFacei << endl; + } + } + } + + if (symmetric_) + { + const auto srcTreePtr = this->createTree(this->srcPatch0_); + const auto& srcTree = srcTreePtr(); + + // Check that all source cells have connections and populate any + // missing entries + forAll(tgtWeights, tgtCelli) + { + if (tgtAddress[tgtCelli].empty()) + { + const point& tgtCc = tgtCcs[tgtCelli]; + pointIndexHit hit = srcTree.findNearest(tgtCc, GREAT); + + if + ( + hit.hit() + && (magSqr(tgtCc - srcCcs[hit.index()]) < maxDistance2_) + ) + { + tgtAddress[tgtCelli] = labelList(1, hit.index()); + tgtWeights[tgtCelli] = scalarList(1, 1.0); + } + } + else + { + if (debug) + { + WarningInFunction + << "Unable to find source face for target face " + << tgtCelli << endl; + } + } + } + } + } + + return true; +} + + +template<class SourcePatch, class TargetPatch> +void Foam::nearestFaceAMI<SourcePatch, TargetPatch>::normaliseWeights +( + const bool verbose, + AMIInterpolation<SourcePatch, TargetPatch>& inter +) +{ + // Do nothing - weights already 1-to-1 and normalised +} + + +// ************************************************************************* // diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/nearestFaceAMI.H similarity index 64% rename from src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.H rename to src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/nearestFaceAMI.H index 791c370332d0d097339a64f9340e916f5c75b93b..78bc71d63693c9c25444c7e9e6ee934ac4d9630d 100644 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.H +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/nearestFaceAMI.H @@ -5,8 +5,7 @@ \\ / A nd | www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2013-2016 OpenFOAM Foundation - Copyright (C) 2016 OpenCFD Ltd. + Copyright (C) 2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -25,20 +24,21 @@ License along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. Class - Foam::mapNearestAMI + Foam::nearestFaceAMI Description - Nearest-mapping Arbitrary Mesh Interface (AMI) method + Nearest-face Arbitrary Mesh Interface (AMI) method SourceFiles - mapNearestAMI.C + nearestFaceAMI.C \*---------------------------------------------------------------------------*/ -#ifndef mapNearestAMI_H -#define mapNearestAMI_H +#ifndef nearestFaceAMI_H +#define nearestFaceAMI_H #include "AMIMethod.H" +#include "pointIndexHit.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -46,73 +46,83 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class mapNearestAMI Declaration + Class nearestFaceAMI Declaration \*---------------------------------------------------------------------------*/ template<class SourcePatch, class TargetPatch> -class mapNearestAMI +class nearestFaceAMI : public AMIMethod<SourcePatch, TargetPatch> { +public: + + typedef Tuple2<pointIndexHit, scalar> nearestAndDist; + + //- Helper class for finding nearest + class nearestEqOp + { + + public: + + void operator()(nearestAndDist& x, const nearestAndDist& y) const + { + if (y.first().hit()) + { + if (!x.first().hit()) + { + x = y; + } + else if (y.second() < x.second()) + { + x = y; + } + } + } + }; + private: - // Private Member Functions + // Private Data - //- No copy construct - mapNearestAMI(const mapNearestAMI&) = delete; + //- Maximum squared distance + scalar maxDistance2_; - //- No copy assignment - void operator=(const mapNearestAMI&) = delete; - // Marching front + // Private Member Functions - //- Find nearest target face for source face srcFacei - void findNearestFace - ( - const SourcePatch& srcPatch, - const TargetPatch& tgtPatch, - const label& srcFacei, - label& tgtFacei - ) const; - - //- Determine next source-target face pair - void setNextNearestFaces - ( - boolList& mapFlag, - label& startSeedI, - label& srcFacei, - label& tgtFacei - ) const; - - //- Find mapped source face - label findMappedSrcFace - ( - const label tgtFacei, - const List<DynamicList<label>>& tgtToSrc - ) const; + //- No copy construct + nearestFaceAMI(const nearestFaceAMI&) = delete; + //- No copy assignment + void operator=(const nearestFaceAMI&) = delete; - // Evaluation + autoPtr<mapDistribute> calcFaceMap + ( + const List<nearestAndDist>& localInfo, + const SourcePatch& srcPatch, + const TargetPatch& tgtPatch + ) const; - //- Area of intersection between source and target faces - scalar interArea - ( - const label srcFacei, - const label tgtFacei - ) const; + autoPtr<mapDistribute> calcDistributed + ( + const SourcePatch& src, + const TargetPatch& tgt, + labelListList& srcToTgtAddr, + scalarListList& srcToTgtWght + ) const; public: //- Runtime type information - TypeName("mapNearestAMI"); + TypeName("nearestFaceAMI"); // Constructors //- Construct from components - mapNearestAMI + nearestFaceAMI ( const SourcePatch& srcPatch, const TargetPatch& tgtPatch, @@ -123,7 +133,7 @@ public: //- Destructor - virtual ~mapNearestAMI() = default; + virtual ~nearestFaceAMI() = default; // Member Functions @@ -146,15 +156,6 @@ public: label tgtFacei = -1 ); - //- Set the face areas for parallel runs - virtual void setMagSf - ( - const TargetPatch& tgtPatch, - const mapDistribute& map, - scalarList& srcMagSf, - scalarList& tgtMagSf - ) const; - //- Normalise the weight. Can optionally subset addressing //- (e.g. for mapNearest) virtual void normaliseWeights @@ -172,7 +173,7 @@ public: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #ifdef NoRepository - #include "mapNearestAMI.C" + #include "nearestFaceAMI.C" #endif // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIPatchToPatchInterpolation.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIPatchToPatchInterpolation.C index 3ea1622bed68b5104889157962713e26dfe5c318..3db4e09a0e6a0998043c82c8048b918455d5c942 100644 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIPatchToPatchInterpolation.C +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIPatchToPatchInterpolation.C @@ -28,7 +28,7 @@ License #include "AMIPatchToPatchInterpolation.H" #include "AMIMethod.H" #include "directAMI.H" -#include "mapNearestAMI.H" +#include "nearestFaceAMI.H" #include "faceAreaWeightAMI.H" #include "partialFaceAreaWeightAMI.H" @@ -39,7 +39,7 @@ namespace Foam makeAMIMethod(AMIPatchToPatchInterpolation); makeAMIMethodType(AMIPatchToPatchInterpolation, directAMI); - makeAMIMethodType(AMIPatchToPatchInterpolation, mapNearestAMI); + makeAMIMethodType(AMIPatchToPatchInterpolation, nearestFaceAMI); makeAMIMethodType(AMIPatchToPatchInterpolation, faceAreaWeightAMI); makeAMIMethodType(AMIPatchToPatchInterpolation, partialFaceAreaWeightAMI); }