Skip to content
Snippets Groups Projects
mappedPatchBase.C 55.2 KiB
Newer Older
  • Learn to ignore specific revisions
  • /*---------------------------------------------------------------------------*\
      =========                 |
      \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
       \\    /   O peration     |
    
    OpenFOAM bot's avatar
    OpenFOAM bot committed
        \\  /    A nd           | www.openfoam.com
    
         \\/     M anipulation  |
    -------------------------------------------------------------------------------
    
    OpenFOAM bot's avatar
    OpenFOAM bot committed
        Copyright (C) 2011-2016 OpenFOAM Foundation
    
        Copyright (C) 2015-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 "mappedPatchBase.H"
    #include "addToRunTimeSelectionTable.H"
    #include "ListListOps.H"
    
    #include "meshSearchMeshObject.H"
    
    #include "meshTools.H"
    #include "OFstream.H"
    #include "Random.H"
    #include "treeDataFace.H"
    
    #include "treeDataPoint.H"
    
    #include "indexedOctree.H"
    #include "polyMesh.H"
    #include "polyPatch.H"
    #include "Time.H"
    #include "mapDistribute.H"
    #include "SubField.H"
    
    #include "triPointRef.H"
    
    #include "syncTools.H"
    
    #include "treeDataCell.H"
    
    #include "DynamicField.H"
    
    #include "faceAreaWeightAMI.H"
    
    #include "OTstream.H"
    
    
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    
    namespace Foam
    {
        defineTypeNameAndDebug(mappedPatchBase, 0);
    }
    
    
    
    const Foam::Enum
    <
        Foam::mappedPatchBase::sampleMode
    >
    Foam::mappedPatchBase::sampleModeNames_
    
        { sampleMode::NEARESTCELL, "nearestCell" },
        { sampleMode::NEARESTPATCHFACE, "nearestPatchFace" },
        { sampleMode::NEARESTPATCHFACEAMI, "nearestPatchFaceAMI" },
        { sampleMode::NEARESTPATCHPOINT, "nearestPatchPoint" },
        { sampleMode::NEARESTFACE, "nearestFace" },
        { sampleMode::NEARESTONLYCELL, "nearestOnlyCell" },
    
    Mark OLESEN's avatar
    Mark OLESEN committed
    });
    
    
    
    const Foam::Enum
    <
        Foam::mappedPatchBase::offsetMode
    >
    Foam::mappedPatchBase::offsetModeNames_
    
        { offsetMode::UNIFORM, "uniform" },
        { offsetMode::NONUNIFORM, "nonuniform" },
        { offsetMode::NORMAL, "normal" },
    
    Mark OLESEN's avatar
    Mark OLESEN committed
    });
    
    
    
    // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
    
    
    Foam::autoPtr<Foam::fileName> Foam::mappedPatchBase::readDatabase
    (
        const dictionary& dict
    )
    {
        autoPtr<fileName> dbNamePtr_;
    
        if (dict.found("sampleDatabase"))
        {
            const bool useDb = dict.get<bool>("sampleDatabase");
            if (useDb)
            {
                dbNamePtr_.set
                (
                    new fileName
                    (
                        dict.lookupOrDefault<fileName>
                        (
                            "sampleDatabasePath",
                            fileName::null
                        )
                    )
                );
            }
        }
        else if (dict.found("sampleDatabasePath"))
        {
            dbNamePtr_.set(new fileName(dict.get<fileName>("sampleDatabasePath")));
        }
    
        return dbNamePtr_;
    }
    
    
    Foam::label Foam::mappedPatchBase::communicator
    (
        const word& sampleWorld
    )
    {
        // Start off with local world
        label comm = UPstream::worldComm;
    
        if (!sampleWorld.empty() && Pstream::parRun())
        {
            if (!UPstream::allWorlds().found(sampleWorld))
            {
                FatalErrorInFunction << "Cannot find sampleWorld " << sampleWorld
                    << " in set of worlds " << UPstream::allWorlds()
                    << exit(FatalError);
            }
    
            const labelList& worldIDs = UPstream::worldIDs();
    
            DynamicList<label> subRanks(worldIDs.size());
            forAll(worldIDs, proci)
            {
                const label worldi = worldIDs[proci];
                if
                (
                    worldi == UPstream::myWorldID()
                 || UPstream::allWorlds()[worldi] == sampleWorld
                )
                {
                    subRanks.append(proci);
                }
            }
    
            // Allocate new communicator with parent 0 (= world)
            comm = UPstream::allocateCommunicator(0, subRanks, true);
    
            if (debug)
            {
                Pout<< "mappedPatchBase::communicator :"
                    << " myWorld:" << UPstream::myWorld()
                    << " sampleWorld:" << sampleWorld
                    << " using subRanks:" << subRanks
                    << " new comm:" << comm << endl;
            }
        }
    
        return comm;
    }
    
    
    
    Foam::tmp<Foam::pointField> Foam::mappedPatchBase::facePoints
    (
        const polyPatch& pp
    ) const
    {
        const polyMesh& mesh = pp.boundaryMesh().mesh();
    
    
        // Force construction of min-tet decomp
        (void)mesh.tetBasePtIs();
    
    
        // Initialise to face-centre
    
        tmp<pointField> tfacePoints(new pointField(patch_.size()));
    
        pointField& facePoints = tfacePoints.ref();
    
            facePoints[facei] = facePoint
    
    void Foam::mappedPatchBase::collectSamples
    (
    
        const label mySampleWorld,      // Wanted world
    
        const pointField& facePoints,
    
        pointField& samples,            // All samples
        labelList& patchFaceWorlds,     // Per sample: wanted world
    
    Mark OLESEN's avatar
    Mark OLESEN committed
        labelList& patchFaceProcs,      // Per sample: originating processor
    
        labelList& patchFaces,          // Per sample: originating patchFace index
        pointField& patchFc             // Per sample: originating centre
    
        const label oldComm(Pstream::warnComm);
        Pstream::warnComm = comm_;
    
        const label myRank = Pstream::myProcNo(comm_);
        const label nProcs = Pstream::nProcs(comm_);
    
    
        // Collect all sample points and the faces they come from.
    
            List<pointField> globalFc(nProcs);
            globalFc[myRank] = facePoints;
            Pstream::gatherList(globalFc, Pstream::msgType(), comm_);
            Pstream::scatterList(globalFc, Pstream::msgType(), comm_);
    
            // Rework into straight list
            patchFc = ListListOps::combine<pointField>
            (
                globalFc,
                accessOp<pointField>()
            );
        }
    
            List<pointField> globalSamples(nProcs);
            globalSamples[myRank] = samplePoints(facePoints);
            Pstream::gatherList(globalSamples, Pstream::msgType(), comm_);
            Pstream::scatterList(globalSamples, Pstream::msgType(), comm_);
    
            // Rework into straight list
            samples = ListListOps::combine<pointField>
            (
                globalSamples,
                accessOp<pointField>()
            );
        }
    
        {
    
            labelListList globalFaces(nProcs);
            globalFaces[myRank] = identity(patch_.size());
    
            // Distribute to all processors
    
            Pstream::gatherList(globalFaces, Pstream::msgType(), comm_);
            Pstream::scatterList(globalFaces, Pstream::msgType(), comm_);
    
    
            patchFaces = ListListOps::combine<labelList>
    
            (
                globalFaces,
                accessOp<labelList>()
    
            labelList procToWorldIndex(nProcs);
            procToWorldIndex[myRank] = mySampleWorld;
            Pstream::gatherList(procToWorldIndex, Pstream::msgType(), comm_);
            Pstream::scatterList(procToWorldIndex, Pstream::msgType(), comm_);
    
            labelList nPerProc(nProcs);
            nPerProc[myRank] = patch_.size();
            Pstream::gatherList(nPerProc, Pstream::msgType(), comm_);
            Pstream::scatterList(nPerProc, Pstream::msgType(), comm_);
    
            patchFaceWorlds.setSize(patchFaces.size());
    
            patchFaceProcs.setSize(patchFaces.size());
    
            label sampleI = 0;
    
                for (label i = 0; i < nPerProc[proci]; i++)
    
                    patchFaceWorlds[sampleI] = procToWorldIndex[proci];
                    patchFaceProcs[sampleI] = proci;
                    sampleI++;
    
        Pstream::warnComm = oldComm;
    
    void Foam::mappedPatchBase::findLocalSamples
    
        const sampleMode mode,
    
    
        const label mySampleWorld,  // local world to sample == my own world
        const word& sampleRegion,   // local region to sample
        const word& samplePatch,    // local patch to sample
    
    
        const pointField& samples,
    
        List<nearInfoWorld>& nearest
    
        // Find the local cell containing the samples
    
        const label myRank = Pstream::myProcNo(comm_);
    
    
        // Lookup the correct region
    
        const polyMesh& mesh = lookupMesh(sampleRegion);
    
    
        // All the info for nearest. Construct to miss
    
        nearest.setSize(samples.size());
        nearInfoWorld miss;
        {
            miss.first().second() = Tuple2<scalar, label>(Foam::sqr(GREAT), -1);
            miss.second() = -1; // set world to be ignored
        }
        nearest = miss;
    
        switch (mode)
    
                if (samplePatch.size() && samplePatch != "none")
    
                    FatalErrorInFunction
                        << "No need to supply a patch name when in "
    
                        << sampleModeNames_[mode] << " mode." << exit(FatalError);
    
                //- Note: face-diagonal decomposition
    
                const indexedOctree<Foam::treeDataCell>& tree = mesh.cellTree();
    
    
                forAll(samples, sampleI)
                {
                    const point& sample = samples[sampleI];
    
                    nearInfoWorld& near = nearest[sampleI];
    
                    label celli = tree.findInside(sample);
    
                        near.first().second().first() = Foam::sqr(GREAT);
                        near.first().second().second() = myRank;
                        near.second() = mySampleWorld;
    
                        const point& cc = mesh.cellCentres()[celli];
    
                        near.first().first() = pointIndexHit
    
                        near.first().second().first() = magSqr(cc-sample);
                        near.first().second().second() = myRank;
                        near.second() = mySampleWorld;
    
            case NEARESTONLYCELL:
            {
    
                if (samplePatch.size() && samplePatch != "none")
    
                    FatalErrorInFunction
                        << "No need to supply a patch name when in "
    
                        << sampleModeNames_[mode] << " mode." << exit(FatalError);
                }
    
                //- Note: face-diagonal decomposition
                const indexedOctree<Foam::treeDataCell>& tree = mesh.cellTree();
    
                forAll(samples, sampleI)
                {
                    const point& sample = samples[sampleI];
    
                    nearInfoWorld& near = nearest[sampleI];
    
                    near.first().first() = tree.findNearest(sample, sqr(GREAT));
                    near.first().second().first() = magSqr
    
                        near.first().first().hitPoint()
    
                    near.first().second().second() = myRank;
                    near.second() = mySampleWorld;
    
            case NEARESTPATCHFACE:
            {
                Random rndGen(123456);
    
    
                const polyPatch& pp = lookupPatch(sampleRegion, samplePatch);
    
    
                if (pp.empty())
                {
                    forAll(samples, sampleI)
                    {
    
                        nearInfoWorld& near = nearest[sampleI];
                        near.first().second().first() = Foam::sqr(GREAT);
                        near.first().second().second() = myRank;
                        near.second() = mySampleWorld;
    
                    }
                }
                else
                {
                    treeBoundBox patchBb
                    (
                        treeBoundBox(pp.points(), pp.meshPoints()).extend
                        (
                            rndGen,
    
                    patchBb.min() -= point::uniform(ROOTVSMALL);
                    patchBb.max() += point::uniform(ROOTVSMALL);
    
    
                    indexedOctree<treeDataFace> boundaryTree
                    (
                        treeDataFace    // all information needed to search faces
                        (
                            false,      // do not cache bb
                            mesh,
    
                            identity(pp.range())  // boundary faces only
    
                        ),
                        patchBb,        // overall search domain
                        8,              // maxLevel
                        10,             // leafsize
                        3.0             // duplicity
                    );
    
                    forAll(samples, sampleI)
                    {
                        const point& sample = samples[sampleI];
    
    
                        nearInfoWorld& near = nearest[sampleI];
                        pointIndexHit& nearInfo = near.first().first();
    
                        nearInfo = boundaryTree.findNearest
                        (
                            sample,
                            magSqr(patchBb.span())
                        );
    
                        if (!nearInfo.hit())
                        {
    
                            near.first().second().first() = Foam::sqr(GREAT);
                            near.first().second().second() = myRank;
                            near.second() = mySampleWorld;
    
                        }
                        else
                        {
                            point fc(pp[nearInfo.index()].centre(pp.points()));
                            nearInfo.setPoint(fc);
    
                            near.first().second().first() = magSqr(fc-sample);
                            near.first().second().second() = myRank;
                            near.second() = mySampleWorld;
    
            case NEARESTPATCHPOINT:
            {
                Random rndGen(123456);
    
    
                const polyPatch& pp = lookupPatch(sampleRegion, samplePatch);
    
    
                if (pp.empty())
                {
                    forAll(samples, sampleI)
                    {
    
                        nearInfoWorld& near = nearest[sampleI];
                        near.first().second().first() = Foam::sqr(GREAT);
                        near.first().second().second() = myRank;
                        near.second() = mySampleWorld;
    
                    }
                }
                else
                {
                    // patch (local) points
                    treeBoundBox patchBb
                    (
                        treeBoundBox(pp.points(), pp.meshPoints()).extend
                        (
                            rndGen,
                            1e-4
                        )
                    );
    
                    patchBb.min() -= point::uniform(ROOTVSMALL);
                    patchBb.max() += point::uniform(ROOTVSMALL);
    
    
                    indexedOctree<treeDataPoint> boundaryTree
                    (
                        treeDataPoint   // all information needed to search faces
                        (
                            mesh.points(),
                            pp.meshPoints() // selection of points to search on
                        ),
                        patchBb,        // overall search domain
                        8,              // maxLevel
                        10,             // leafsize
                        3.0             // duplicity
                    );
    
                    forAll(samples, sampleI)
                    {
                        const point& sample = samples[sampleI];
    
    
                        nearInfoWorld& near = nearest[sampleI];
                        pointIndexHit& nearInfo = near.first().first();
    
                        nearInfo = boundaryTree.findNearest
                        (
                            sample,
                            magSqr(patchBb.span())
                        );
    
                        if (!nearInfo.hit())
                        {
    
                            near.first().second().first() = Foam::sqr(GREAT);
                            near.first().second().second() = myRank;
                            near.second() = mySampleWorld;
    
                        }
                        else
                        {
                            const point& pt = nearInfo.hitPoint();
    
    
                            near.first().second().first() = magSqr(pt-sample);
                            near.first().second().second() = myRank;
                            near.second() = mySampleWorld;
    
                if (samplePatch.size() && samplePatch != "none")
    
                    FatalErrorInFunction
                        << "No need to supply a patch name when in "
    
                        << sampleModeNames_[mode] << " mode." << exit(FatalError);
    
                //- Note: face-diagonal decomposition
                const meshSearchMeshObject& meshSearchEngine =
                    meshSearchMeshObject::New(mesh);
    
    
                forAll(samples, sampleI)
                {
                    const point& sample = samples[sampleI];
    
                    nearInfoWorld& near = nearest[sampleI];
    
                    label facei = meshSearchEngine.findNearestFace(sample);
    
                        near.first().second().first() = Foam::sqr(GREAT);
                        near.first().second().second() = myRank;
                        near.second() = mySampleWorld;
    
                        const point& fc = mesh.faceCentres()[facei];
    
                        near.first().first() = pointIndexHit(true, fc, facei);
                        near.first().second().first() = magSqr(fc-sample);
                        near.first().second().second() = myRank;
                        near.second() = mySampleWorld;
    
                    }
                }
                break;
            }
    
            case NEARESTPATCHFACEAMI:
            {
                // nothing to do here
                return;
            }
    
            default:
            {
    
                FatalErrorInFunction
    
                    << "problem." << abort(FatalError);
            }
        }
    
    void Foam::mappedPatchBase::findSamples
    (
        const sampleMode mode,
        const label myWorld,
        const pointField& samples,
        const labelList& wantedWorlds,
        const labelList& origProcs,
    
        labelList& sampleProcs,
        labelList& sampleIndices,
        pointField& sampleLocations
    ) const
    {
        // Find the processor/cell containing the samples. Does not account
        // for samples being found in two processors.
    
        const label myRank = Pstream::myProcNo(comm_);
        const label nProcs = Pstream::nProcs(comm_);
    
        wordList samplePatches(nProcs);
        {
            const label oldComm(Pstream::warnComm);
            Pstream::warnComm = comm_;
            samplePatches[myRank] = samplePatch_;
            Pstream::gatherList(samplePatches, Pstream::msgType(), comm_);
            Pstream::scatterList(samplePatches, Pstream::msgType(), comm_);
            Pstream::warnComm = oldComm;
        }
        wordList sampleRegions(nProcs);
    
            const label oldComm(Pstream::warnComm);
            Pstream::warnComm = comm_;
            sampleRegions[myRank] = sampleRegion_;
            Pstream::gatherList(sampleRegions, Pstream::msgType(), comm_);
            Pstream::scatterList(sampleRegions, Pstream::msgType(), comm_);
            Pstream::warnComm = oldComm;
        }
    
        // Find all the info for nearest
        List<nearInfoWorld> nearest(samples.size());
        forAll(nearest, samplei)
        {
            nearest[samplei].first() = nearInfo
            (
                pointIndexHit(),
                Tuple2<scalar, label>(Foam::sqr(GREAT), -1)
            );
            nearest[samplei].second() = wantedWorlds[samplei];
        }
    
    
        // Extract samples to search for locally
        {
            DynamicList<label> localMap(samples.size());
            forAll(wantedWorlds, samplei)
    
                if (wantedWorlds[samplei] == myWorld)
                {
                    localMap.append(samplei);
                }
            }
    
            if (localMap.size())
            {
                pointField localSamples(samples, localMap);
                labelList localOrigProcs(origProcs, localMap);
    
                //Assume single patch to sample for now
                const word localOrigPatch(samplePatches[localOrigProcs[0]]);
                const word localOrigRegion(sampleRegions[localOrigProcs[0]]);
                List<nearInfoWorld> localNearest(localSamples.size());
    
                if (debug)
                {
                    Pout<< "Searching locally for " << localSamples.size()
                        << " samples on region:" << localOrigRegion
                        << " on patch:" << localOrigPatch << endl;
                }
                findLocalSamples
                (
                    mode,
                    myWorld,
                    localOrigRegion,
                    localOrigPatch,
                    localSamples,
                    localNearest
                );
                UIndirectList<nearInfoWorld>(nearest, localMap) = localNearest;
    
        const label oldComm(Pstream::warnComm);
        Pstream::warnComm = comm_;
    
        // Find nearest. Combine on master.
        Pstream::listCombineGather
        (
            nearest,
            nearestWorldEqOp(),
            Pstream::msgType(),
            comm_
        );
        Pstream::listCombineScatter(nearest, Pstream::msgType(), comm_);
    
        //if (debug)
        //{
        //    Pout<< "** AFter combining:" << endl;
        //    forAll(nearest, samplei)
        //    {
        //        Pout<< "  sample:" << samples[samplei]
        //            << " origating from proc:" << origProcs[samplei]
        //            << " to be found on world:" << wantedWorlds[samplei] << nl
        //            << "    found on world:" << nearest[samplei].second() << nl
        //            << "    found on proc:"
        //            << nearest[samplei].first().second().second() << nl
        //            << "    found on patchfacei:"
        //            << nearest[samplei].first().first().index() << nl
        //            << "    found at location:"
        //            << nearest[samplei].first().first().rawPoint() << nl;
        //    }
        //    Pout<< endl;
        //}
    
    
        // Convert back into proc+local index
        sampleProcs.setSize(samples.size());
        sampleIndices.setSize(samples.size());
        sampleLocations.setSize(samples.size());
    
        forAll(nearest, sampleI)
        {
    
            const nearInfo& ni = nearest[sampleI].first();
    
            if (!ni.first().hit())
    
            {
                sampleProcs[sampleI] = -1;
                sampleIndices[sampleI] = -1;
                sampleLocations[sampleI] = vector::max;
            }
            else
            {
    
                sampleProcs[sampleI] = ni.second().second();
                sampleIndices[sampleI] = ni.first().index();
                sampleLocations[sampleI] = ni.first().hitPoint();
    
    
        Pstream::warnComm = oldComm;
    
    }
    
    
    void Foam::mappedPatchBase::calcMapping() const
    {
    
        static bool hasWarned = false;
    
            FatalErrorInFunction
    
                << "Mapping already calculated" << exit(FatalError);
        }
    
    
        //// Make sure if running in database that there is a syncObjects FO
        //if (sampleDatabase() && !sameWorld())
        //{
        //    const word syncName("syncObjects");
        //    const polyMesh& mesh = patch_.boundaryMesh().mesh();
        //    const Time& runTime = mesh.time();
        //    const functionObjectList& fos = runTime.functionObjects();
        //
        //    forAll(fos, i)
        //    {
        //        Pout<< "** FO:" << fos[i].name() << " tpye:" << fos[i].type()
        //            << endl;
        //    }
        //
        //    if (fos.findObjectID(syncName) == -1)
        //    {
        //        //FatalErrorInFunction
        //        WarningInFunction
        //            << "Did not detect functionObject " << syncName
        //            << ". This is used to synchronise data inbetween worlds"
        //            << endl;
        //    }
        //}
    
    
    
        // Get points on face (since cannot use face-centres - might be off
        // face-diagonal decomposed tets.
        tmp<pointField> patchPoints(facePoints(patch_));
    
        // Get offsetted points
    
        const pointField offsettedPoints(samplePoints(patchPoints()));
    
        // Do a sanity check - am I sampling my own patch?
        // This only makes sense for a non-zero offset.
    
        bool sampleMyself =
        (
            mode_ == NEARESTPATCHFACE
    
         && sampleWorld() == UPstream::myWorld()
    
         && sampleRegion() == patch_.boundaryMesh().mesh().name()
         && samplePatch() == patch_.name()
    
        if (sampleMyself)
    
            // Check offset
            vectorField d(offsettedPoints-patchPoints());
            bool coincident = (gAverage(mag(d)) <= ROOTVSMALL);
    
            if (sampleMyself && coincident)
            {
                WarningInFunction
                    << "Invalid offset " << d << endl
                    << "Offset is the vector added to the patch face centres to"
                    << " find the patch face supplying the data." << endl
                    << "Setting it to " << d
                    << " on the same patch, on the same region"
                    << " will find the faces themselves which does not make sense"
                    << " for anything but testing." << endl
                    << "patch_:" << patch_.name() << endl
                    << "sampleRegion_:" << sampleRegion() << endl
                    << "mode_:" << sampleModeNames_[mode_] << endl
                    << "samplePatch_:" << samplePatch() << endl
                    << "offsetMode_:" << offsetModeNames_[offsetMode_] << endl;
            }
    
        // Get local world and world-to-sample in index form
        const label myWorld = Pstream::myWorldID();
        const label mySampleWorld = Pstream::allWorlds().find(sampleWorld_);
    
    
    
        // Get global list of all samples and the processor and face they come from.
    
        pointField samples;         // coordinates
        labelList patchFaceWorlds;  // world to sample
        labelList patchFaceProcs;   // originating processor
        labelList patchFaces;       // originating face on processor patch
        pointField patchFc;         // originating face centre
    
            mySampleWorld,          // world I want to sample
    
            patchFaceWorlds,
    
            patchFaceProcs,
            patchFaces,
            patchFc
        );
    
    
        //if (debug)
        //{
        //    forAll(samples, samplei)
        //    {
        //        Pout<< "    sample:" << samples[samplei]
        //            << " origating from proc:" << patchFaceProcs[samplei]
        //            << "  face:" << patchFaces[samplei]
        //            << " to be found on world:" << patchFaceWorlds[samplei] << nl;
        //    }
        //}
    
    
        // Find processor and cell/face samples are in and actual location.
        labelList sampleProcs;
        labelList sampleIndices;
        pointField sampleLocations;
    
        findSamples
        (
            mode_,
            myWorld,        // my world (in index form)
            samples,
            patchFaceWorlds,
            patchFaceProcs,
    
            sampleProcs,
            sampleIndices,
            sampleLocations
        );
    
        // Check for samples that were not found. This will only happen for
        // NEARESTCELL since finds cell containing a location
        if (mode_ == NEARESTCELL)
        {
            label nNotFound = 0;
            forAll(sampleProcs, sampleI)
            {
                if (sampleProcs[sampleI] == -1)
                {
                    nNotFound++;
                }
            }
    
            reduce(nNotFound, sumOp<label>(), Pstream::msgType(), comm_);
    
    
            if (nNotFound > 0)
            {
                if (!hasWarned)
                {
    
                    WarningInFunction
                        << "Did not find " << nNotFound
    
                        << " out of " << sampleProcs.size() << " total samples."
                        << " Sampling these on owner cell centre instead." << endl
    
                        << "On patch " << patch_.name()
    
                        << " on region " << sampleRegion()
    
                        << " in mode " << sampleModeNames_[mode_] << endl
    
                        << "with offset mode " << offsetModeNames_[offsetMode_]
                        << ". Suppressing further warnings from " << type() << endl;
    
                // Collect the samples that cannot be found
                DynamicList<label> subMap;
    
                forAll(sampleProcs, sampleI)
                {
                    if (sampleProcs[sampleI] == -1)
                    {
    
                        subMap.append(sampleI);
    
                // And re-search for pure nearest (should not fail)
                labelList subSampleProcs;
                labelList subSampleIndices;
                pointField subSampleLocations;
                findSamples
                (
                    NEARESTONLYCELL,
    
                    myWorld,        // my world (in index form)
    
                    pointField(samples, subMap),
                    UIndirectList<label>(patchFaceWorlds, subMap)(),
                    UIndirectList<label>(patchFaceProcs, subMap)(),
    
    
                    subSampleProcs,
                    subSampleIndices,
                    subSampleLocations
                );
    
                // Insert
    
                labelUIndList(sampleProcs, subMap) = subSampleProcs;
                labelUIndList(sampleIndices, subMap) = subSampleIndices;
    
                UIndirectList<point>(sampleLocations, subMap) = subSampleLocations;
    
        // Now we have all the data we need:
        // - where sample originates from (so destination when mapping):
        //   patchFaces, patchFaceProcs.
        // - cell/face sample is in (so source when mapping)
        //   sampleIndices, sampleProcs.
    
    
        if (Pstream::master(comm_))
        {
            forAll(samples, i)
            {
                if (sampleProcs[i] == -1)
                {
                    FatalErrorInFunction << "did not find sample "
                        << patchFc[i] << " on patch " << patch_.name()
                        << " on region "
                        << patch_.boundaryMesh().mesh().name()
                        << " on processor " << patchFaceProcs[i]
                        << exit(FatalError);
                }
            }
        }
    
        if (debug && Pstream::master(comm_))
    
            //forAll(samples, i)
            //{
            //    Pout<< i << " need data in region "
            //        << patch_.boundaryMesh().mesh().name()
            //        << " for proc:" << patchFaceProcs[i]
            //        << " face:" << patchFaces[i]
            //        << " at:" << patchFc[i] << endl
            //        << "Found data in region " << sampleRegion()
            //        << " at proc:" << sampleProcs[i]
            //        << " face:" << sampleIndices[i]
            //        << " at:" << sampleLocations[i]
            //        << nl << endl;
            //}
    
    
            OFstream str
            (
                patch_.boundaryMesh().mesh().time().path()
              / patch_.name()
              + "_mapped.obj"
            );
            Pout<< "Dumping mapping as lines from patch faceCentres to"
    
                << " sampled cell/faceCentres/points to file " << str.name()
                << endl;
    
    
            label vertI = 0;
    
            forAll(patchFc, i)
            {
                meshTools::writeOBJ(str, patchFc[i]);
                vertI++;
                meshTools::writeOBJ(str, sampleLocations[i]);
                vertI++;
                str << "l " << vertI-1 << ' ' << vertI << nl;
            }
        }
    
        // Determine schedule.
    
        mapPtr_.reset(new mapDistribute(sampleProcs, patchFaceProcs, comm_));
    
    
        // Rework the schedule from indices into samples to cell data to send,
        // face data to receive.