nearestFaceAMI.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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"
#include "addToRunTimeSelectionTable.H"
 
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
namespace Foam
{
    defineTypeNameAndDebug(nearestFaceAMI, 0);
    addToRunTimeSelectionTable(AMIInterpolation, nearestFaceAMI, dict);
    addToRunTimeSelectionTable(AMIInterpolation, nearestFaceAMI, component);
}
 
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
Foam::autoPtr<Foam::mapDistribute> Foam::nearestFaceAMI::calcFaceMap
(
    const List<nearestAndDist>& localInfo,
    const primitivePatch& srcPatch,
    const primitivePatch& tgtPatch
) const
{
    // Generate the list of processor bounding boxes for tgtPatch
    List<boundBox> procBbs(Pstream::nProcs());
    procBbs[Pstream::myProcNo()] =
        boundBox(tgtPatch.points(), tgtPatch.meshPoints(), true);
    Pstream::gatherList(procBbs);
    Pstream::scatterList(procBbs);
 
    // Identify which of my local src faces intersect with each processor
    // tgtPatch bb within the current match's search distance
    const pointField& srcCcs = srcPatch.faceCentres();
    List<DynamicList<label>> dynSendMap(Pstream::nProcs());
 
    forAll(localInfo, srcFacei)
    {
        // Test local srcPatch face centres against remote processor tgtPatch bb
        // using distance from initial pass
 
        const scalar r2 = localInfo[srcFacei].second();
 
        forAll(procBbs, proci)
        {
            if (proci != Pstream::myProcNo())
            {
                if (procBbs[proci].overlaps(srcCcs[srcFacei], r2))
                {
                    dynSendMap[proci].append(srcFacei);
                }
            }
        }
    }
 
    // 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));
}
 
 
Foam::autoPtr<Foam::mapDistribute> Foam::nearestFaceAMI::calcDistributed
(
    const primitivePatch& src,
    const primitivePatch& tgt,
    labelListList& srcToTgtAddr,
    scalarListList& srcToTgtWght
) const
{
    autoPtr<indexedOctree<treeType>> tgtTreePtr;
    if (tgt.size())
    {
        tgtTreePtr = this->createTree(tgt);
    }
 
    // Create global indexing for tgtPatch
    globalIndex globalTgtCells(tgt.size());
 
 
    // First pass
    // ==========
    // For each srcPatch face, determine local match on tgtPatch
 
    // Identify local nearest matches
    pointField srcCcs(src.faceCentres());
 
    List<nearestAndDist> localInfo(src.size());
    if (tgt.size())
    {
        const auto& tgtTree = tgtTreePtr();
 
        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()));
            }
        }
    }
    else
    {
        // No local tgtPatch faces - initialise nearest distance for all
        // srcPatch faces to GREAT so that they [should be] set by remote
        // tgtPatch faces
        for (auto& info : localInfo)
        {
            info.second() = GREAT;
        }
    }
 
 
    // Second pass
    // ===========
    // Determine remote matches
 
    // Map returns labels of src patch faces to send to each proc
    autoPtr<mapDistribute> mapPtr = calcFaceMap(localInfo, src, tgt);
    mapDistribute& map = mapPtr();
 
    List<nearestAndDist> remoteInfo(localInfo);
    map.distribute(remoteInfo);
 
    // Note: re-using srcCcs
    map.distribute(srcCcs);
 
    if (tgt.size())
    {
        const auto& tgtTree = tgtTreePtr();
 
        // Test remote srcPatch faces against local tgtPatch faces
        nearestAndDist testInfo;
        pointIndexHit& test = testInfo.first();
        forAll(remoteInfo, i)
        {
            test = tgtTree.findNearest(srcCcs[i], remoteInfo[i].second());
            if (test.hit())
            {
                test.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,
        nearestZero,
        nearestEqOp(),
        noOp()              // no flipping
    );
 
 
    // 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, srcFacei)
    {
        nearestEqOp()(localInfo[srcFacei], remoteInfo[srcFacei]);
        if (localInfo[srcFacei].second() < maxDistance2_)
        {
            const label tgtFacei = localInfo[srcFacei].first().index();
            srcToTgtAddr[srcFacei] = labelList(1, tgtFacei);
            srcToTgtWght[srcFacei] = scalarList(1, 1.0);
        }
    }
 
    List<Map<label>> cMap;
    return autoPtr<mapDistribute>::New(globalTgtCells, srcToTgtAddr, cMap);
}
 
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
Foam::nearestFaceAMI::nearestFaceAMI
(
    const dictionary& dict,
    const bool reverseTarget
)
:
    AMIInterpolation(dict, reverseTarget),
    maxDistance2_(dict.getOrDefault<scalar>("maxDistance2", GREAT))
{}
 
 
Foam::nearestFaceAMI::nearestFaceAMI
(
    const bool requireMatch,
    const bool reverseTarget,
    const scalar lowWeightCorrection
)
:
    AMIInterpolation(requireMatch, reverseTarget, lowWeightCorrection),
    maxDistance2_(GREAT)
{}
 
 
Foam::nearestFaceAMI::nearestFaceAMI(const nearestFaceAMI& ami)
:
    AMIInterpolation(ami),
    maxDistance2_(ami.maxDistance2_)
{}
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
bool Foam::nearestFaceAMI::calculate
(
    const primitivePatch& srcPatch,
    const primitivePatch& tgtPatch,
    const autoPtr<searchableSurface>& surfPtr
)
{
    if (upToDate_)
    {
        return false;
    }
 
    AMIInterpolation::calculate(srcPatch, tgtPatch, surfPtr);
 
    const auto& src = this->srcPatch0();
    const auto& tgt = this->tgtPatch0();
 
    // Set face area magnitudes
    srcMagSf_ = mag(src.faceAreas());
    tgtMagSf_ = mag(tgt.faceAreas());
 
    // TODO: implement symmetric calculation controls; assume yes for now
    bool symmetric_ = true;
 
    if (this->distributed())
    {
        tgtMapPtr_ =
            calcDistributed
            (
                src,
                tgt,
                srcAddress_,
                srcWeights_
            );
 
        if (symmetric_)
        {
            srcMapPtr_ =
                calcDistributed
                (
                    tgt,
                    src,
                    tgtAddress_,
                    tgtWeights_
                );
            }
    }
    else
    {
        srcAddress_.setSize(src.size());
        srcWeights_.setSize(src.size());
 
        if (symmetric_)
        {
            tgtAddress_.setSize(tgt.size());
            tgtWeights_.setSize(tgt.size());
        }
 
        const pointField& srcCcs = src.faceCentres();
        const pointField& tgtCcs = tgt.faceCentres();
 
        const auto tgtTreePtr = this->createTree(tgtPatch);
        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(srcPatch);
            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;
                    }
                }
            }
        }
    }
 
    srcWeightsSum_.setSize(srcWeights_.size(), 1);
    tgtWeightsSum_.setSize(tgtWeights_.size(), 1);
 
    upToDate_ = true;
 
    return upToDate_;
}
 
 
// ************************************************************************* //