From 3cd4bc9c09116dbc1ccdf46f393dba06ceb0209d Mon Sep 17 00:00:00 2001
From: Andrew Heather <>
Date: Wed, 8 Apr 2020 07:59:57 +0100
Subject: [PATCH] ENH: mapNearestAMI - re-written to use tree search instead of
advancing front
---
.../AMIMethod/AMIMethod/AMIMethod.C | 46 +-
.../AMIMethod/AMIMethod/AMIMethod.H | 7 +-
.../AMIMethod/mapNearestAMI/mapNearestAMI.C | 444 ------------------
.../AMIMethod/mapNearestAMI/nearestFaceAMI.C | 361 ++++++++++++++
.../{mapNearestAMI.H => nearestFaceAMI.H} | 119 ++---
.../AMIPatchToPatchInterpolation.C | 4 +-
6 files changed, 450 insertions(+), 531 deletions(-)
delete mode 100644 src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.C
create mode 100644 src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/nearestFaceAMI.C
rename src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/{mapNearestAMI.H => nearestFaceAMI.H} (64%)
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.C
index f4265f155ab..5b3d1bc483b 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 e702f8366ab..62efb5d4868 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 e2d57b8599d..00000000000
--- 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 00000000000..0aa96a35c44
--- /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 791c370332d..78bc71d6369 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 3ea1622bed6..3db4e09a0e6 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);
}
--
GitLab