Commit 3cd4bc9c authored by Andrew Heather's avatar Andrew Heather Committed by Andrew Heather
Browse files

ENH: mapNearestAMI - re-written to use tree search instead of advancing front

parent c357d7a7
......@@ -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
);
}
......
......@@ -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
(
......
/*---------------------------------------------------------------------------*\
========= |
\\ / 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);
}
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / 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();