/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  dev                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 1;

vertices
(
    (0 0 0)
    (4 0 0)
    (0 1 0)
    (4 1 0)
    (0 2 0)
    (4 2 0)
    (0 0 2)
    (4 0 2)
    (0 1 2)
    (4 1 2)
    (0 2 2)
    (4 2 2)
);

blocks
(
    hex (0 1 3 2 6 7 9 8) (400 250 30) simpleGrading (1 10.7028 1)
    hex (2 3 5 4 8 9 11 10) (400 250 30) simpleGrading (1 0.0934 1)
);

edges
(
);

boundary
(
    bottomWall
    {
        type            wall;
        faces           ((0 1 7 6));
    }
    topWall
    {
        type            wall;
        faces           ((4 10 11 5));
    }

    sides1_half0
    {
        type            cyclic;
        neighbourPatch  sides1_half1;
        faces           ((0 2 3 1));
    }
    sides1_half1
    {
        type            cyclic;
        neighbourPatch  sides1_half0;
        faces           ((6 7 9 8));
    }

    sides2_half0
    {
        type            cyclic;
        neighbourPatch  sides2_half1;
        faces           ((2 4 5 3));
    }
    sides2_half1
    {
        type            cyclic;
        neighbourPatch  sides2_half0;
        faces           ((8 9 11 10));
    }

    inout1_half0
    {
        type            cyclic;
        neighbourPatch  inout1_half1;
        faces           ((1 3 9 7));
    }
    inout1_half1
    {
        type            cyclic;
        neighbourPatch  inout1_half0;
        faces           ((0 6 8 2));
    }

    inout2_half0
    {
        type            cyclic;
        neighbourPatch  inout2_half1;
        faces           ((3 5 11 9));
    }
    inout2_half1
    {
        type            cyclic;
        neighbourPatch  inout2_half0;
        faces           ((2 8 10 4));
    }
);

mergePatchPairs
(
);

// ************************************************************************* //

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
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 "FaceCellWave.H"
#include "polyMesh.H"
#include "processorPolyPatch.H"
#include "cyclicPolyPatch.H"
#include "cyclicAMIPolyPatch.H"
#include "OPstream.H"
#include "IPstream.H"
#include "PstreamReduceOps.H"
#include "debug.H"
#include "typeInfo.H"
#include "SubField.H"
#include "globalMeshData.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

template<class Type, class TrackingData>
const Foam::scalar Foam::FaceCellWave<Type, TrackingData>::geomTol_ = 1e-6;

template<class Type, class TrackingData>
Foam::scalar Foam::FaceCellWave<Type, TrackingData>::propagationTol_ = 0.01;

template<class Type, class TrackingData>
int Foam::FaceCellWave<Type, TrackingData>::dummyTrackData_ = 12345;

namespace Foam
{
    template<class Type, class TrackingData>
    class combine
    {
        //- Combine operator for AMIInterpolation

        FaceCellWave<Type, TrackingData>& solver_;

        const cyclicAMIPolyPatch& patch_;

        public:

            combine
            (
                FaceCellWave<Type, TrackingData>& solver,
                const cyclicAMIPolyPatch& patch
            )
            :
                solver_(solver),
                patch_(patch)
            {}


            void operator()
            (
                Type& x,
                const label facei,
                const Type& y,
                const scalar weight
            ) const
            {
                if (y.valid(solver_.data()))
                {
                    label meshFacei = -1;
                    if (patch_.owner())
                    {
                        meshFacei = patch_.start() + facei;
                    }
                    else
                    {
                        meshFacei = patch_.neighbPatch().start() + facei;
                    }
                    x.updateFace
                    (
                        solver_.mesh(),
                        meshFacei,
                        y,
                        solver_.propagationTol(),
                        solver_.data()
                    );
                }
            }
    };
}


// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

template<class Type, class TrackingData>
bool Foam::FaceCellWave<Type, TrackingData>::updateCell
(
    const label celli,
    const label neighbourFacei,
    const Type& neighbourInfo,
    const scalar tol,
    Type& cellInfo
)
{
    // Update info for celli, at position pt, with information from
    // neighbouring face/cell.
    // Updates:
    //      - changedCell_, changedCells_
    //      - statistics: nEvals_, nUnvisitedCells_

    nEvals_++;

    bool wasValid = cellInfo.valid(td_);

    bool propagate =
        cellInfo.updateCell
        (
            mesh_,
            celli,
            neighbourFacei,
            neighbourInfo,
            tol,
            td_
        );

    if (propagate)
    {
        if (!changedCell_[celli])
        {
            changedCell_[celli] = true;
            changedCells_.append(celli);
        }
    }

    if (!wasValid && cellInfo.valid(td_))
    {
        --nUnvisitedCells_;
    }

    return propagate;
}


template<class Type, class TrackingData>
bool Foam::FaceCellWave<Type, TrackingData>::updateFace
(
    const label facei,
    const label neighbourCelli,
    const Type& neighbourInfo,
    const scalar tol,
    Type& faceInfo
)
{
    // Update info for facei, at position pt, with information from
    // neighbouring face/cell.
    // Updates:
    //      - changedFace_, changedFaces_,
    //      - statistics: nEvals_, nUnvisitedFaces_

    nEvals_++;

    bool wasValid = faceInfo.valid(td_);

    bool propagate =
        faceInfo.updateFace
        (
            mesh_,
            facei,
            neighbourCelli,
            neighbourInfo,
            tol,
            td_
        );

    if (propagate)
    {
        if (!changedFace_[facei])
        {
            changedFace_[facei] = true;
            changedFaces_.append(facei);
        }
    }

    if (!wasValid && faceInfo.valid(td_))
    {
        --nUnvisitedFaces_;
    }

    return propagate;
}


template<class Type, class TrackingData>
bool Foam::FaceCellWave<Type, TrackingData>::updateFace
(
    const label facei,
    const Type& neighbourInfo,
    const scalar tol,
    Type& faceInfo
)
{
    // Update info for facei, at position pt, with information from
    // same face.
    // Updates:
    //      - changedFace_, changedFaces_,
    //      - statistics: nEvals_, nUnvisitedFaces_

    nEvals_++;

    bool wasValid = faceInfo.valid(td_);

    bool propagate =
        faceInfo.updateFace
        (
            mesh_,
            facei,
            neighbourInfo,
            tol,
            td_
        );

    if (propagate)
    {
        if (!changedFace_[facei])
        {
            changedFace_[facei] = true;
            changedFaces_.append(facei);
        }
    }

    if (!wasValid && faceInfo.valid(td_))
    {
        --nUnvisitedFaces_;
    }

    return propagate;
}


template<class Type, class TrackingData>
void Foam::FaceCellWave<Type, TrackingData>::checkCyclic
(
    const polyPatch& patch
) const
{
    // For debugging: check status on both sides of cyclic

    const cyclicPolyPatch& nbrPatch =
        refCast<const cyclicPolyPatch>(patch).neighbPatch();

    forAll(patch, patchFacei)
    {
        label i1 = patch.start() + patchFacei;
        label i2 = nbrPatch.start() + patchFacei;

        if
        (
           !allFaceInfo_[i1].sameGeometry
            (
                mesh_,
                allFaceInfo_[i2],
                geomTol_,
                td_
            )
        )
        {
            FatalErrorInFunction
                << "   faceInfo:" << allFaceInfo_[i1]
                << "   otherfaceInfo:" << allFaceInfo_[i2]
                << abort(FatalError);
        }

        if (changedFace_[i1] != changedFace_[i2])
        {
            FatalErrorInFunction
                << "   faceInfo:" << allFaceInfo_[i1]
                << "   otherfaceInfo:" << allFaceInfo_[i2]
                << "   changedFace:" << changedFace_[i1]
                << "   otherchangedFace:" << changedFace_[i2]
                << abort(FatalError);
        }
    }
}


template<class Type, class TrackingData>
template<class PatchType>
bool Foam::FaceCellWave<Type, TrackingData>::hasPatch() const
{
    forAll(mesh_.boundaryMesh(), patchi)
    {
        if (isA<PatchType>(mesh_.boundaryMesh()[patchi]))
        {
            return true;
        }
    }
    return false;
}


template<class Type, class TrackingData>
void Foam::FaceCellWave<Type, TrackingData>::setFaceInfo
(
    const labelList& changedFaces,
    const List<Type>& changedFacesInfo
)
{
    forAll(changedFaces, changedFacei)
    {
        label facei = changedFaces[changedFacei];

        bool wasValid = allFaceInfo_[facei].valid(td_);

        // Copy info for facei
        allFaceInfo_[facei] = changedFacesInfo[changedFacei];

        // Maintain count of unset faces
        if (!wasValid && allFaceInfo_[facei].valid(td_))
        {
            --nUnvisitedFaces_;
        }

        // Mark facei as changed, both on list and on face itself.

        changedFace_[facei] = true;
        changedFaces_.append(facei);
    }
}


template<class Type, class TrackingData>
void Foam::FaceCellWave<Type, TrackingData>::mergeFaceInfo
(
    const polyPatch& patch,
    const label nFaces,
    const labelList& changedFaces,
    const List<Type>& changedFacesInfo
)
{
    // Merge face information into member data

    for (label changedFacei = 0; changedFacei < nFaces; changedFacei++)
    {
        const Type& neighbourWallInfo = changedFacesInfo[changedFacei];
        label patchFacei = changedFaces[changedFacei];

        label meshFacei = patch.start() + patchFacei;

        Type& currentWallInfo = allFaceInfo_[meshFacei];

        if (!currentWallInfo.equal(neighbourWallInfo, td_))
        {
            updateFace
            (
                meshFacei,
                neighbourWallInfo,
                propagationTol_,
                currentWallInfo
            );
        }
    }
}


template<class Type, class TrackingData>
Foam::label Foam::FaceCellWave<Type, TrackingData>::getChangedPatchFaces
(
    const polyPatch& patch,
    const label startFacei,
    const label nFaces,
    labelList& changedPatchFaces,
    List<Type>& changedPatchFacesInfo
) const
{
    // Construct compact patchFace change arrays for a (slice of a) single
    // patch. changedPatchFaces in local patch numbering.
    // Return length of arrays.
    label nChangedPatchFaces = 0;

    for (label i = 0; i < nFaces; i++)
    {
        label patchFacei = i + startFacei;

        label meshFacei = patch.start() + patchFacei;

        if (changedFace_[meshFacei])
        {
            changedPatchFaces[nChangedPatchFaces] = patchFacei;
            changedPatchFacesInfo[nChangedPatchFaces] = allFaceInfo_[meshFacei];
            nChangedPatchFaces++;
        }
    }
    return nChangedPatchFaces;
}


template<class Type, class TrackingData>
void Foam::FaceCellWave<Type, TrackingData>::leaveDomain
(
    const polyPatch& patch,
    const label nFaces,
    const labelList& faceLabels,
    List<Type>& faceInfo
) const
{
    // Handle leaving domain. Implementation referred to Type

    const vectorField& fc = mesh_.faceCentres();

    for (label i = 0; i < nFaces; i++)
    {
        label patchFacei = faceLabels[i];

        label meshFacei = patch.start() + patchFacei;
        faceInfo[i].leaveDomain(mesh_, patch, patchFacei, fc[meshFacei], td_);
    }
}


template<class Type, class TrackingData>
void Foam::FaceCellWave<Type, TrackingData>::enterDomain
(
    const polyPatch& patch,
    const label nFaces,
    const labelList& faceLabels,
    List<Type>& faceInfo
) const
{
    // Handle entering domain. Implementation referred to Type

    const vectorField& fc = mesh_.faceCentres();

    for (label i = 0; i < nFaces; i++)
    {
        label patchFacei = faceLabels[i];

        label meshFacei = patch.start() + patchFacei;
        faceInfo[i].enterDomain(mesh_, patch, patchFacei, fc[meshFacei], td_);
    }
}


template<class Type, class TrackingData>
void Foam::FaceCellWave<Type, TrackingData>::transform
(
    const tensorField& rotTensor,
    const label nFaces,
    List<Type>& faceInfo
)
{
    // Transform. Implementation referred to Type

    if (rotTensor.size() == 1)
    {
        const tensor& T = rotTensor[0];

        for (label facei = 0; facei < nFaces; facei++)
        {
            faceInfo[facei].transform(mesh_, T, td_);
        }
    }
    else
    {
        for (label facei = 0; facei < nFaces; facei++)
        {
            faceInfo[facei].transform(mesh_, rotTensor[facei], td_);
        }
    }
}


template<class Type, class TrackingData>
void Foam::FaceCellWave<Type, TrackingData>::offset
(
    const polyPatch&,
    const label cycOffset,
    const label nFaces,
    labelList& faces
)
{
    // Offset mesh face. Used for transferring from one cyclic half to the
    // other.

    for (label facei = 0; facei < nFaces; facei++)
    {
        faces[facei] += cycOffset;
    }
}


template<class Type, class TrackingData>
void Foam::FaceCellWave<Type, TrackingData>::handleProcPatches()
{
    // Tranfer all the information to/from neighbouring processors

    const globalMeshData& pData = mesh_.globalData();

    // Which patches are processor patches
    const labelList& procPatches = pData.processorPatches();

    // Send all

    PstreamBuffers pBufs(Pstream::nonBlocking);

    forAll(procPatches, i)
    {
        label patchi = procPatches[i];

        const processorPolyPatch& procPatch =
            refCast<const processorPolyPatch>(mesh_.boundaryMesh()[patchi]);

        // Allocate buffers
        label nSendFaces;
        labelList sendFaces(procPatch.size());
        List<Type> sendFacesInfo(procPatch.size());

        // Determine which faces changed on current patch
        nSendFaces = getChangedPatchFaces
        (
            procPatch,
            0,
            procPatch.size(),
            sendFaces,
            sendFacesInfo
        );

        // Adapt wallInfo for leaving domain
        leaveDomain
        (
            procPatch,
            nSendFaces,
            sendFaces,
            sendFacesInfo
        );

        if (debug & 2)
        {
            Pout<< " Processor patch " << patchi << ' ' << procPatch.name()
                << " communicating with " << procPatch.neighbProcNo()
                << "  Sending:" << nSendFaces
                << endl;
        }

        UOPstream toNeighbour(procPatch.neighbProcNo(), pBufs);
        //writeFaces(nSendFaces, sendFaces, sendFacesInfo, toNeighbour);
        toNeighbour
            << SubList<label>(sendFaces, nSendFaces)
            << SubList<Type>(sendFacesInfo, nSendFaces);
    }

    pBufs.finishedSends();

    // Receive all

    forAll(procPatches, i)
    {
        label patchi = procPatches[i];

        const processorPolyPatch& procPatch =
            refCast<const processorPolyPatch>(mesh_.boundaryMesh()[patchi]);

        // Allocate buffers
        labelList receiveFaces;
        List<Type> receiveFacesInfo;

        {
            UIPstream fromNeighbour(procPatch.neighbProcNo(), pBufs);
            fromNeighbour >> receiveFaces >> receiveFacesInfo;
        }

        if (debug & 2)
        {
            Pout<< " Processor patch " << patchi << ' ' << procPatch.name()
                << " communicating with " << procPatch.neighbProcNo()
                << "  Receiving:" << receiveFaces.size()
                << endl;
        }

        // Apply transform to received data for non-parallel planes
        if (!procPatch.parallel())
        {
            transform
            (
                procPatch.forwardT(),
                receiveFaces.size(),
                receiveFacesInfo
            );
        }

        // Adapt wallInfo for entering domain
        enterDomain
        (
            procPatch,
            receiveFaces.size(),
            receiveFaces,
            receiveFacesInfo
        );

        // Merge received info
        mergeFaceInfo
        (
            procPatch,
            receiveFaces.size(),
            receiveFaces,
            receiveFacesInfo
        );
    }
}


template<class Type, class TrackingData>
void Foam::FaceCellWave<Type, TrackingData>::handleCyclicPatches()
{
    // Transfer information across cyclic halves.

    forAll(mesh_.boundaryMesh(), patchi)
    {
        const polyPatch& patch = mesh_.boundaryMesh()[patchi];

        if (isA<cyclicPolyPatch>(patch))
        {
            const cyclicPolyPatch& nbrPatch =
                refCast<const cyclicPolyPatch>(patch).neighbPatch();

            // Allocate buffers
            label nReceiveFaces;
            labelList receiveFaces(patch.size());
            List<Type> receiveFacesInfo(patch.size());

            // Determine which faces changed
            nReceiveFaces = getChangedPatchFaces
            (
                nbrPatch,
                0,
                nbrPatch.size(),
                receiveFaces,
                receiveFacesInfo
            );

            // Adapt wallInfo for leaving domain
            leaveDomain
            (
                nbrPatch,
                nReceiveFaces,
                receiveFaces,
                receiveFacesInfo
            );

            const cyclicPolyPatch& cycPatch =
                refCast<const cyclicPolyPatch>(patch);

            if (!cycPatch.parallel())
            {
                // received data from other half
                transform
                (
                    cycPatch.forwardT(),
                    nReceiveFaces,
                    receiveFacesInfo
                );
            }

            if (debug & 2)
            {
                Pout<< " Cyclic patch " << patchi << ' ' << cycPatch.name()
                    << "  Changed : " << nReceiveFaces
                    << endl;
            }

            // Half2: Adapt wallInfo for entering domain
            enterDomain
            (
                cycPatch,
                nReceiveFaces,
                receiveFaces,
                receiveFacesInfo
            );

            // Merge into global storage
            mergeFaceInfo
            (
                cycPatch,
                nReceiveFaces,
                receiveFaces,
                receiveFacesInfo
            );

            if (debug)
            {
                checkCyclic(cycPatch);
            }
        }
    }
}


template<class Type, class TrackingData>
void Foam::FaceCellWave<Type, TrackingData>::handleAMICyclicPatches()
{
    // Transfer information across cyclicAMI halves.

    forAll(mesh_.boundaryMesh(), patchi)
    {
        const polyPatch& patch = mesh_.boundaryMesh()[patchi];

        if (isA<cyclicAMIPolyPatch>(patch))
        {
            const cyclicAMIPolyPatch& cycPatch =
                refCast<const cyclicAMIPolyPatch>(patch);

            List<Type> receiveInfo;

            {
                const cyclicAMIPolyPatch& nbrPatch =
                    refCast<const cyclicAMIPolyPatch>(patch).neighbPatch();

                // Get nbrPatch data (so not just changed faces)
                typename List<Type>::subList sendInfo
                (
                    nbrPatch.patchSlice
                    (
                        allFaceInfo_
                    )
                );

                if (!nbrPatch.parallel() || nbrPatch.separated())
                {
                    // Adapt sendInfo for leaving domain
                    const vectorField::subField fc = nbrPatch.faceCentres();
                    forAll(sendInfo, i)
                    {
                        sendInfo[i].leaveDomain(mesh_, nbrPatch, i, fc[i], td_);
                    }
                }

                // Transfer sendInfo to cycPatch
                combine<Type, TrackingData> cmb(*this, cycPatch);

                if (cycPatch.applyLowWeightCorrection())
                {
                    List<Type> defVals
                    (
                        cycPatch.patchInternalList(allCellInfo_)
                    );

                    cycPatch.interpolate(sendInfo, cmb, receiveInfo, defVals);
                }
                else
                {
                    cycPatch.interpolate(sendInfo, cmb, receiveInfo);
                }
            }

            // Apply transform to received data for non-parallel planes
            if (!cycPatch.parallel())
            {
                transform
                (
                    cycPatch.forwardT(),
                    receiveInfo.size(),
                    receiveInfo
                );
            }

            if (!cycPatch.parallel() || cycPatch.separated())
            {
                // Adapt receiveInfo for entering domain
                const vectorField::subField fc = cycPatch.faceCentres();
                forAll(receiveInfo, i)
                {
                    receiveInfo[i].enterDomain(mesh_, cycPatch, i, fc[i], td_);
                }
            }

            // Merge into global storage
            forAll(receiveInfo, i)
            {
                label meshFacei = cycPatch.start()+i;

                Type& currentWallInfo = allFaceInfo_[meshFacei];

                if
                (
                    receiveInfo[i].valid(td_)
                && !currentWallInfo.equal(receiveInfo[i], td_)
                )
                {
                    updateFace
                    (
                        meshFacei,
                        receiveInfo[i],
                        propagationTol_,
                        currentWallInfo
                    );
                }
            }
        }
    }
}


template<class Type, class TrackingData>
void Foam::FaceCellWave<Type, TrackingData>::handleExplicitConnections()
{
    // Collect changed information

    DynamicList<label> f0Baffle(explicitConnections_.size());
    DynamicList<Type> f0Info(explicitConnections_.size());

    DynamicList<label> f1Baffle(explicitConnections_.size());
    DynamicList<Type> f1Info(explicitConnections_.size());

    forAll(explicitConnections_, connI)
    {
        const labelPair& baffle = explicitConnections_[connI];

        label f0 = baffle[0];
        if (changedFace_[f0])
        {
            f0Baffle.append(connI);
            f0Info.append(allFaceInfo_[f0]);
        }

        label f1 = baffle[1];
        if (changedFace_[f1])
        {
            f1Baffle.append(connI);
            f1Info.append(allFaceInfo_[f1]);
        }
    }


    // Update other side with changed information

    forAll(f1Info, index)
    {
        const labelPair& baffle = explicitConnections_[f1Baffle[index]];

        label f0 = baffle[0];
        Type& currentWallInfo = allFaceInfo_[f0];

        if (!currentWallInfo.equal(f1Info[index], td_))
        {
            updateFace
            (
                f0,
                f1Info[index],
                propagationTol_,
                currentWallInfo
            );
        }
    }

    forAll(f0Info, index)
    {
        const labelPair& baffle = explicitConnections_[f0Baffle[index]];

        label f1 = baffle[1];
        Type& currentWallInfo = allFaceInfo_[f1];

        if (!currentWallInfo.equal(f0Info[index], td_))
        {
            updateFace
            (
                f1,
                f0Info[index],
                propagationTol_,
                currentWallInfo
            );
        }
    }
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class Type, class TrackingData>
Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
(
    const polyMesh& mesh,
    UList<Type>& allFaceInfo,
    UList<Type>& allCellInfo,
    TrackingData& td
)
:
    mesh_(mesh),
    explicitConnections_(0),
    allFaceInfo_(allFaceInfo),
    allCellInfo_(allCellInfo),
    td_(td),
    changedFace_(mesh_.nFaces(), false),
    changedFaces_(mesh_.nFaces()),
    changedCell_(mesh_.nCells(), false),
    changedCells_(mesh_.nCells()),
    hasCyclicPatches_(hasPatch<cyclicPolyPatch>()),
    hasCyclicAMIPatches_
    (
        returnReduce(hasPatch<cyclicAMIPolyPatch>(), orOp<bool>())
    ),
    nEvals_(0),
    nUnvisitedCells_(mesh_.nCells()),
    nUnvisitedFaces_(mesh_.nFaces())
{
    if
    (
        allFaceInfo.size() != mesh_.nFaces()
     || allCellInfo.size() != mesh_.nCells()
    )
    {
        FatalErrorInFunction
            << "face and cell storage not the size of mesh faces, cells:"
            << endl
            << "    allFaceInfo   :" << allFaceInfo.size() << endl
            << "    mesh_.nFaces():" << mesh_.nFaces() << endl
            << "    allCellInfo   :" << allCellInfo.size() << endl
            << "    mesh_.nCells():" << mesh_.nCells()
            << exit(FatalError);
    }
}


template<class Type, class TrackingData>
Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
(
    const polyMesh& mesh,
    const labelList& changedFaces,
    const List<Type>& changedFacesInfo,
    UList<Type>& allFaceInfo,
    UList<Type>& allCellInfo,
    const label maxIter,
    TrackingData& td
)
:
    mesh_(mesh),
    explicitConnections_(0),
    allFaceInfo_(allFaceInfo),
    allCellInfo_(allCellInfo),
    td_(td),
    changedFace_(mesh_.nFaces(), false),
    changedFaces_(mesh_.nFaces()),
    changedCell_(mesh_.nCells(), false),
    changedCells_(mesh_.nCells()),
    hasCyclicPatches_(hasPatch<cyclicPolyPatch>()),
    hasCyclicAMIPatches_
    (
        returnReduce(hasPatch<cyclicAMIPolyPatch>(), orOp<bool>())
    ),
    nEvals_(0),
    nUnvisitedCells_(mesh_.nCells()),
    nUnvisitedFaces_(mesh_.nFaces())
{
    if
    (
        allFaceInfo.size() != mesh_.nFaces()
     || allCellInfo.size() != mesh_.nCells()
    )
    {
        FatalErrorInFunction
            << "face and cell storage not the size of mesh faces, cells:"
            << endl
            << "    allFaceInfo   :" << allFaceInfo.size() << endl
            << "    mesh_.nFaces():" << mesh_.nFaces() << endl
            << "    allCellInfo   :" << allCellInfo.size() << endl
            << "    mesh_.nCells():" << mesh_.nCells()
            << exit(FatalError);
    }

    // Copy initial changed faces data
    setFaceInfo(changedFaces, changedFacesInfo);

    // Iterate until nothing changes
    label iter = iterate(maxIter);

    if ((maxIter > 0) && (iter >= maxIter))
    {
        FatalErrorInFunction
            << "Maximum number of iterations reached. Increase maxIter." << endl
            << "    maxIter:" << maxIter << endl
            << "    nChangedCells:" << changedCells_.size() << endl
            << "    nChangedFaces:" << changedFaces_.size() << endl
            << exit(FatalError);
    }
}


template<class Type, class TrackingData>
Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
(
    const polyMesh& mesh,
    const labelPairList& explicitConnections,
    const bool handleCyclicAMI,
    const labelList& changedFaces,
    const List<Type>& changedFacesInfo,
    UList<Type>& allFaceInfo,
    UList<Type>& allCellInfo,
    const label maxIter,
    TrackingData& td
)
:
    mesh_(mesh),
    explicitConnections_(explicitConnections),
    allFaceInfo_(allFaceInfo),
    allCellInfo_(allCellInfo),
    td_(td),
    changedFace_(mesh_.nFaces(), false),
    changedFaces_(mesh_.nFaces()),
    changedCell_(mesh_.nCells(), false),
    changedCells_(mesh_.nCells()),
    hasCyclicPatches_(hasPatch<cyclicPolyPatch>()),
    hasCyclicAMIPatches_
    (
        handleCyclicAMI
     && returnReduce(hasPatch<cyclicAMIPolyPatch>(), orOp<bool>())
    ),
    nEvals_(0),
    nUnvisitedCells_(mesh_.nCells()),
    nUnvisitedFaces_(mesh_.nFaces())
{
    if
    (
        allFaceInfo.size() != mesh_.nFaces()
     || allCellInfo.size() != mesh_.nCells()
    )
    {
        FatalErrorInFunction
            << "face and cell storage not the size of mesh faces, cells:"
            << endl
            << "    allFaceInfo   :" << allFaceInfo.size() << endl
            << "    mesh_.nFaces():" << mesh_.nFaces() << endl
            << "    allCellInfo   :" << allCellInfo.size() << endl
            << "    mesh_.nCells():" << mesh_.nCells()
            << exit(FatalError);
    }

    // Copy initial changed faces data
    setFaceInfo(changedFaces, changedFacesInfo);

    // Iterate until nothing changes
    label iter = iterate(maxIter);

    if ((maxIter > 0) && (iter >= maxIter))
    {
        FatalErrorInFunction
            << "Maximum number of iterations reached. Increase maxIter." << endl
            << "    maxIter:" << maxIter << endl
            << "    nChangedCells:" << changedCells_.size() << endl
            << "    nChangedFaces:" << changedFaces_.size() << endl
            << exit(FatalError);
    }
}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class Type, class TrackingData>
Foam::label Foam::FaceCellWave<Type, TrackingData>::getUnsetCells() const
{
    return nUnvisitedCells_;
}


template<class Type, class TrackingData>
Foam::label Foam::FaceCellWave<Type, TrackingData>::getUnsetFaces() const
{
    return nUnvisitedFaces_;
}


template<class Type, class TrackingData>
Foam::label Foam::FaceCellWave<Type, TrackingData>::faceToCell()
{
    // Propagate face to cell

    const labelList& owner = mesh_.faceOwner();
    const labelList& neighbour = mesh_.faceNeighbour();
    label nInternalFaces = mesh_.nInternalFaces();

    forAll(changedFaces_, changedFacei)
    {
        label facei = changedFaces_[changedFacei];
        if (!changedFace_[facei])
        {
            FatalErrorInFunction
                << "Face " << facei
                << " not marked as having been changed"
                << abort(FatalError);
        }


        const Type& neighbourWallInfo = allFaceInfo_[facei];

        // Evaluate all connected cells

        // Owner
        label celli = owner[facei];
        Type& currentWallInfo = allCellInfo_[celli];

        if (!currentWallInfo.equal(neighbourWallInfo, td_))
        {
            updateCell
            (
                celli,
                facei,
                neighbourWallInfo,
                propagationTol_,
                currentWallInfo
            );
        }

        // Neighbour.
        if (facei < nInternalFaces)
        {
            celli = neighbour[facei];
            Type& currentWallInfo2 = allCellInfo_[celli];

            if (!currentWallInfo2.equal(neighbourWallInfo, td_))
            {
                updateCell
                (
                    celli,
                    facei,
                    neighbourWallInfo,
                    propagationTol_,
                    currentWallInfo2
                );
            }
        }

        // Reset status of face
        changedFace_[facei] = false;
    }

    // Handled all changed faces by now
    changedFaces_.clear();

    if (debug & 2)
    {
        Pout<< " Changed cells            : " << changedCells_.size() << endl;
    }

    // Sum changedCells over all procs
    label totNChanged = changedCells_.size();

    reduce(totNChanged, sumOp<label>());

    return totNChanged;
}


template<class Type, class TrackingData>
Foam::label Foam::FaceCellWave<Type, TrackingData>::cellToFace()
{
    // Propagate cell to face

    const cellList& cells = mesh_.cells();

    forAll(changedCells_, changedCelli)
    {
        label celli = changedCells_[changedCelli];
        if (!changedCell_[celli])
        {
            FatalErrorInFunction
                << "Cell " << celli << " not marked as having been changed"
                << abort(FatalError);
        }

        const Type& neighbourWallInfo = allCellInfo_[celli];

        // Evaluate all connected faces

        const labelList& faceLabels = cells[celli];
        forAll(faceLabels, faceLabelI)
        {
            label facei = faceLabels[faceLabelI];
            Type& currentWallInfo = allFaceInfo_[facei];

            if (!currentWallInfo.equal(neighbourWallInfo, td_))
            {
                updateFace
                (
                    facei,
                    celli,
                    neighbourWallInfo,
                    propagationTol_,
                    currentWallInfo
                );
            }
        }

        // Reset status of cell
        changedCell_[celli] = false;
    }

    // Handled all changed cells by now
    changedCells_.clear();


    // Transfer across any explicitly provided internal connections
    handleExplicitConnections();

    if (hasCyclicPatches_)
    {
        // Transfer changed faces across cyclic halves
        handleCyclicPatches();
    }

    if (hasCyclicAMIPatches_)
    {
        handleAMICyclicPatches();
    }

    if (Pstream::parRun())
    {
        // Transfer changed faces from neighbouring processors.
        handleProcPatches();
    }

    if (debug & 2)
    {
        Pout<< " Changed faces            : " << changedFaces_.size() << endl;
    }

    // Sum nChangedFaces over all procs
    label totNChanged = changedFaces_.size();

    reduce(totNChanged, sumOp<label>());

    return totNChanged;
}


// Iterate
template<class Type, class TrackingData>
Foam::label Foam::FaceCellWave<Type, TrackingData>::iterate(const label maxIter)
{
    if (hasCyclicPatches_)
    {
        // Transfer changed faces across cyclic halves
        handleCyclicPatches();
    }

    if (hasCyclicAMIPatches_)
    {
        handleAMICyclicPatches();
    }

    if (Pstream::parRun())
    {
        // Transfer changed faces from neighbouring processors.
        handleProcPatches();
    }

    label iter = 0;

    while (iter < maxIter)
    {
        if (debug)
        {
            Info<< " Iteration " << iter << endl;
        }

        nEvals_ = 0;

        label nCells = faceToCell();

        if (debug)
        {
            Info<< " Total changed cells      : " << nCells << endl;
        }

        if (nCells == 0)
        {
            break;
        }

        label nFaces = cellToFace();

        if (debug)
        {
            Info<< " Total changed faces      : " << nFaces << nl
                << " Total evaluations        : " << nEvals_ << nl
                << " Remaining unvisited cells: " << nUnvisitedCells_ << nl
                << " Remaining unvisited faces: " << nUnvisitedFaces_ << endl;
        }

        if (nFaces == 0)
        {
            break;
        }

        ++iter;
    }

    return iter;
}


// ************************************************************************* //

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
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/>.

Class
    Foam::FaceCellWave

Description
    Wave propagation of information through grid. Every iteration
    information goes through one layer of cells. Templated on information
    that is transferred.

    Handles parallel and cyclics and non-parallel cyclics.

    Note: whether to propagate depends on the return value of Type::update
    which returns true (i.e. propagate) if the value changes by more than a
    certain tolerance.
    This tolerance can be very strict for normal face-cell and parallel
    cyclics (we use a value of 0.01 just to limit propagation of small changes)
    but for non-parallel cyclics this tolerance can be critical and if chosen
    too small can lead to non-convergence.

SourceFiles
    FaceCellWave.C

\*---------------------------------------------------------------------------*/

#ifndef FaceCellWave_H
#define FaceCellWave_H

#include "PackedBoolList.H"
#include "DynamicList.H"
#include "primitiveFieldsFwd.H"
#include "labelPair.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

// Forward declaration of classes
class polyMesh;
class polyPatch;

/*---------------------------------------------------------------------------*\
                        Class FaceCellWaveName Declaration
\*---------------------------------------------------------------------------*/

TemplateName(FaceCellWave);


/*---------------------------------------------------------------------------*\
                           Class FaceCellWave Declaration
\*---------------------------------------------------------------------------*/

template<class Type, class TrackingData = int>
class FaceCellWave
:
    public FaceCellWaveName
{
private:

    // Private Member Functions

        //- Disallow default bitwise copy construct
        FaceCellWave(const FaceCellWave&);

        //- Disallow default bitwise assignment
        void operator=(const FaceCellWave&);


protected:

    // Private data

        //- Reference to mesh
        const polyMesh& mesh_;

        //- Optional boundary faces that information should travel through
        const labelPairList explicitConnections_;

        //- Information for all faces
        UList<Type>& allFaceInfo_;

        //- Information for all cells
        UList<Type>& allCellInfo_;

        //- Additional data to be passed into container
        TrackingData& td_;

        //- Has face changed
        PackedBoolList changedFace_;

        //- List of changed faces
        DynamicList<label> changedFaces_;

        //- Has cell changed
        PackedBoolList changedCell_;

        // Cells that have changed
        DynamicList<label> changedCells_;


        //- Contains cyclics
        const bool hasCyclicPatches_;

        //- Contains cyclicAMI
        const bool hasCyclicAMIPatches_;

        //- Number of evaluations
        label nEvals_;

        //- Number of unvisited cells/faces
        label nUnvisitedCells_;
        label nUnvisitedFaces_;


        //- Updates cellInfo with information from neighbour. Updates all
        //  statistics.
        bool updateCell
        (
            const label celli,
            const label neighbourFacei,
            const Type& neighbourInfo,
            const scalar tol,
            Type& cellInfo
        );

        //- Updates faceInfo with information from neighbour. Updates all
        //  statistics.
        bool updateFace
        (
            const label facei,
            const label neighbourCelli,
            const Type& neighbourInfo,
            const scalar tol,
            Type& faceInfo
        );

        //- Updates faceInfo with information from same face. Updates all
        //  statistics.
        bool updateFace
        (
            const label facei,
            const Type& neighbourInfo,
            const scalar tol,
            Type& faceInfo
        );


        // Parallel, cyclic

            //- Debugging: check info on both sides of cyclic
            void checkCyclic(const polyPatch& pPatch) const;

            //- Has cyclic patch?
            template<class PatchType>
            bool hasPatch() const;

            //- Merge received patch data into global data
            void mergeFaceInfo
            (
                const polyPatch& patch,
                const label nFaces,
                const labelList&,
                const List<Type>&
            );

            //- Extract info for single patch only
            label getChangedPatchFaces
            (
                const polyPatch& patch,
                const label startFacei,
                const label nFaces,
                labelList& changedPatchFaces,
                List<Type>& changedPatchFacesInfo
            ) const;

            //- Handle leaving domain. Implementation referred to Type
            void leaveDomain
            (
                const polyPatch& patch,
                const label nFaces,
                const labelList& faceLabels,
                List<Type>& faceInfo
            ) const;

            //- Handle leaving domain. Implementation referred to Type
            void enterDomain
            (
                const polyPatch& patch,
                const label nFaces,
                const labelList& faceLabels,
                List<Type>& faceInfo
            ) const;

            //- Offset face labels by constant value
            static void offset
            (
                const polyPatch& patch,
                const label off,
                const label nFaces,
                labelList& faces
            );

            //- Apply transformation to Type
            void transform
            (
                const tensorField& rotTensor,
                const label nFaces,
                List<Type>& faceInfo
            );

            //- Merge data from across processor boundaries
            void handleProcPatches();

            //- Merge data from across cyclics
            void handleCyclicPatches();

            //- Merge data from across AMI cyclics
            void handleAMICyclicPatches();

            //- Merge data across explicitly provided local connections (usually
            //  baffles)
            void handleExplicitConnections();


      // Private static data

            static const scalar geomTol_;
            static scalar propagationTol_;

            //- Used as default trackdata value to satisfy default template
            //  argument.
            static int dummyTrackData_;


public:

    // Static Functions

        //- Access to tolerance
        static scalar propagationTol()
        {
            return propagationTol_;
        }

        //- Change tolerance
        static void setPropagationTol(const scalar tol)
        {
            propagationTol_ = tol;
        }


    // Constructors

        // Construct from mesh. Use setFaceInfo and iterate() to do actual
        // calculation.
        FaceCellWave
        (
            const polyMesh&,
            UList<Type>& allFaceInfo,
            UList<Type>& allCellInfo,
            TrackingData& td = dummyTrackData_
        );

        //- Construct from mesh and list of changed faces with the Type
        //  for these faces. Iterates until nothing changes or maxIter reached.
        //  (maxIter can be 0)
        FaceCellWave
        (
            const polyMesh&,
            const labelList& initialChangedFaces,
            const List<Type>& changedFacesInfo,
            UList<Type>& allFaceInfo,
            UList<Type>& allCellInfo,
            const label maxIter,
            TrackingData& td = dummyTrackData_
        );

        //- Construct from mesh and explicitly connected boundary faces
        //  and list of changed faces with the Type
        //  for these faces. Iterates until nothing changes or maxIter reached.
        //  (maxIter can be 0)
        FaceCellWave
        (
            const polyMesh&,
            const labelPairList& explicitConnections,
            const bool handleCyclicAMI,
            const labelList& initialChangedFaces,
            const List<Type>& changedFacesInfo,
            UList<Type>& allFaceInfo,
            UList<Type>& allCellInfo,
            const label maxIter,
            TrackingData& td = dummyTrackData_
        );


    //- Destructor
    virtual ~FaceCellWave()
    {}


    // Member Functions

        // Access

            //- Access allFaceInfo
            UList<Type>& allFaceInfo()
            {
                return allFaceInfo_;
            }

            //- Access allCellInfo
            UList<Type>& allCellInfo()
            {
                return allCellInfo_;
            }

            //- Additional data to be passed into container
            const TrackingData& data() const
            {
                return td_;
            }

            //- Access mesh
            const polyMesh& mesh() const
            {
                return mesh_;
            }

            //- Get number of unvisited cells, i.e. cells that were not (yet)
            //  reached from walking across mesh. This can happen from
            //  - not enough iterations done
            //  - a disconnected mesh
            //  - a mesh without walls in it
            label getUnsetCells() const;

            //- Get number of unvisited faces
            label getUnsetFaces() const;


        // Edit

            //- Set initial changed faces
            void setFaceInfo
            (
                const labelList& changedFaces,
                const List<Type>& changedFacesInfo
            );

            //- Propagate from face to cell. Returns total number of cells
            //  (over all processors) changed.
           virtual label faceToCell();

            //- Propagate from cell to face. Returns total number of faces
            //  (over all processors) changed. (Faces on processorpatches are
            //  counted double)
            virtual label cellToFace();

            //- Iterate until no changes or maxIter reached.  Returns actual
            //  number of iterations.
            virtual label iterate(const label maxIter);
};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#ifdef NoRepository
    #include "FaceCellWave.C"
#endif

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
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 "FaceCellWave.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
defineTypeNameAndDebug(FaceCellWaveName, 0);
}


// ************************************************************************* //

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
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 "regionSplit.H"
#include "FaceCellWave.H"
#include "cyclicPolyPatch.H"
#include "processorPolyPatch.H"
#include "globalIndex.H"
#include "syncTools.H"
#include "minData.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{

defineTypeNameAndDebug(regionSplit, 0);

}

// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

void Foam::regionSplit::calcNonCompactRegionSplit
(
    const globalIndex& globalFaces,
    const boolList& blockedFace,
    const List<labelPair>& explicitConnections,

    labelList& cellRegion
) const
{
    // Field on cells and faces.
    List<minData> cellData(mesh().nCells());
    List<minData> faceData(mesh().nFaces());

    // Take over blockedFaces by seeding a negative number
    // (so is always less than the decomposition)
    label nUnblocked = 0;
    forAll(faceData, facei)
    {
        if (blockedFace.size() && blockedFace[facei])
        {
            faceData[facei] = minData(-2);
        }
        else
        {
            nUnblocked++;
        }
    }

    // Seed unblocked faces
    labelList seedFaces(nUnblocked);
    List<minData> seedData(nUnblocked);
    nUnblocked = 0;


    forAll(faceData, facei)
    {
        if (blockedFace.empty() || !blockedFace[facei])
        {
            seedFaces[nUnblocked] = facei;
            // Seed face with globally unique number
            seedData[nUnblocked] = minData(globalFaces.toGlobal(facei));
            nUnblocked++;
        }
    }


    // Propagate information inwards
    FaceCellWave<minData> deltaCalc
    (
        mesh(),
        explicitConnections,
        false,  // disable walking through cyclicAMI for backwards compatibility
        seedFaces,
        seedData,
        faceData,
        cellData,
        mesh().globalData().nTotalCells()+1
    );


    // And extract
    cellRegion.setSize(mesh().nCells());
    forAll(cellRegion, celli)
    {
        if (cellData[celli].valid(deltaCalc.data()))
        {
            cellRegion[celli] = cellData[celli].data();
        }
        else
        {
            // Unvisited cell -> only possible if surrounded by blocked faces.
            // If so make up region from any of the faces
            const cell& cFaces = mesh().cells()[celli];
            label facei = cFaces[0];

            if (blockedFace.size() && !blockedFace[facei])
            {
                FatalErrorIn("regionSplit::calcNonCompactRegionSplit(..)")
                    << "Problem: unblocked face " << facei
                    << " at " << mesh().faceCentres()[facei]
                    << " on unassigned cell " << celli
                    << mesh().cellCentres()[celli]
                    << exit(FatalError);
            }
            cellRegion[celli] = globalFaces.toGlobal(facei);
        }
    }
}


Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
(
    const bool doGlobalRegions,
    const boolList& blockedFace,
    const List<labelPair>& explicitConnections,

    labelList& cellRegion
) const
{
    // See header in regionSplit.H


    if (!doGlobalRegions)
    {
        // Block all parallel faces to avoid comms across
        boolList coupledOrBlockedFace(blockedFace);
        const polyBoundaryMesh& pbm = mesh().boundaryMesh();

        if (coupledOrBlockedFace.size())
        {
            forAll(pbm, patchi)
            {
                const polyPatch& pp = pbm[patchi];
                if (isA<processorPolyPatch>(pp))
                {
                    label facei = pp.start();
                    forAll(pp, i)
                    {
                        coupledOrBlockedFace[facei++] = true;
                    }
                }
            }
        }

        // Create dummy (local only) globalIndex
        labelList offsets(Pstream::nProcs()+1, 0);
        for (label i = Pstream::myProcNo()+1; i < offsets.size(); i++)
        {
            offsets[i] = mesh().nFaces();
        }
        const globalIndex globalRegions(offsets.xfer());

        // Minimise regions across connected cells
        // Note: still uses global decisions so all processors are running
        //       in lock-step, i.e. slowest determines overall time.
        //       To avoid this we could switch off Pstream::parRun.
        calcNonCompactRegionSplit
        (
            globalRegions,
            coupledOrBlockedFace,
            explicitConnections,
            cellRegion
        );

        // Compact
        Map<label> globalToCompact(mesh().nCells()/8);
        forAll(cellRegion, celli)
        {
            label region = cellRegion[celli];

            label globalRegion;

            Map<label>::const_iterator fnd = globalToCompact.find(region);
            if (fnd == globalToCompact.end())
            {
                globalRegion = globalRegions.toGlobal(globalToCompact.size());
                globalToCompact.insert(region, globalRegion);
            }
            else
            {
                globalRegion = fnd();
            }
            cellRegion[celli] = globalRegion;
        }


        // Return globalIndex with size = localSize and all regions local
        labelList compactOffsets(Pstream::nProcs()+1, 0);
        for (label i = Pstream::myProcNo()+1; i < compactOffsets.size(); i++)
        {
            compactOffsets[i] = globalToCompact.size();
        }

        return autoPtr<globalIndex>(new globalIndex(compactOffsets.xfer()));
    }



    // Initial global region numbers
    const globalIndex globalRegions(mesh().nFaces());

    // Minimise regions across connected cells (including parallel)
    calcNonCompactRegionSplit
    (
        globalRegions,
        blockedFace,
        explicitConnections,
        cellRegion
    );


    // Now our cellRegion will have
    // - non-local regions (i.e. originating from other processors)
    // - non-compact locally originating regions
    // so we'll need to compact

    // 4a: count per originating processor the number of regions
    labelList nOriginating(Pstream::nProcs(), 0);
    {
        labelHashSet haveRegion(mesh().nCells()/8);

        forAll(cellRegion, celli)
        {
            label region = cellRegion[celli];

            // Count originating processor. Use isLocal as efficiency since
            // most cells are locally originating.
            if (globalRegions.isLocal(region))
            {
                if (haveRegion.insert(region))
                {
                    nOriginating[Pstream::myProcNo()]++;
                }
            }
            else
            {
                label proci = globalRegions.whichProcID(region);
                if (haveRegion.insert(region))
                {
                    nOriginating[proci]++;
                }
            }
        }
    }

    if (debug)
    {
        Pout<< "Counted " << nOriginating[Pstream::myProcNo()]
            << " local regions." << endl;
    }


    // Global numbering for compacted local regions
    autoPtr<globalIndex> globalCompactPtr
    (
        new globalIndex(nOriginating[Pstream::myProcNo()])
    );
    const globalIndex& globalCompact = globalCompactPtr();


    // 4b: renumber
    // Renumber into compact indices. Note that since we've already made
    // all regions global we now need a Map to store the compacting information
    // instead of a labelList - otherwise we could have used a straight
    // labelList.

    // Local compaction map
    Map<label> globalToCompact(2*nOriginating[Pstream::myProcNo()]);
    // Remote regions we want the compact number for
    List<labelHashSet> nonLocal(Pstream::nProcs());
    forAll(nonLocal, proci)
    {
        if (proci != Pstream::myProcNo())
        {
            nonLocal[proci].resize(2*nOriginating[proci]);
        }
    }

    forAll(cellRegion, celli)
    {
        label region = cellRegion[celli];
        if (globalRegions.isLocal(region))
        {
            // Insert new compact region (if not yet present)
            globalToCompact.insert
            (
                region,
                globalCompact.toGlobal(globalToCompact.size())
            );
        }
        else
        {
            nonLocal[globalRegions.whichProcID(region)].insert(region);
        }
    }


    // Now we have all the local regions compacted. Now we need to get the
    // non-local ones from the processors to whom they are local.
    // Convert the nonLocal (labelHashSets) to labelLists.

    labelListList sendNonLocal(Pstream::nProcs());
    forAll(sendNonLocal, proci)
    {
        sendNonLocal[proci] = nonLocal[proci].toc();
    }

    if (debug)
    {
        forAll(sendNonLocal, proci)
        {
            Pout<< "    from processor " << proci
                << " want " << sendNonLocal[proci].size()
                << " region numbers."
                << endl;
        }
        Pout<< endl;
    }


    // Get the wanted region labels into recvNonLocal
    labelListList recvNonLocal;
    Pstream::exchange<labelList, label>(sendNonLocal, recvNonLocal);

    // Now we have the wanted compact region labels that proci wants in
    // recvNonLocal[proci]. Construct corresponding list of compact
    // region labels to send back.

    labelListList sendWantedLocal(Pstream::nProcs());
    forAll(recvNonLocal, proci)
    {
        const labelList& nonLocal = recvNonLocal[proci];
        sendWantedLocal[proci].setSize(nonLocal.size());

        forAll(nonLocal, i)
        {
            sendWantedLocal[proci][i] = globalToCompact[nonLocal[i]];
        }
    }


    // Send back (into recvNonLocal)
    recvNonLocal.clear();
    Pstream::exchange<labelList, label>(sendWantedLocal, recvNonLocal);
    sendWantedLocal.clear();

    // Now recvNonLocal contains for every element in setNonLocal the
    // corresponding compact number. Insert these into the local compaction
    // map.

    forAll(recvNonLocal, proci)
    {
        const labelList& wantedRegions = sendNonLocal[proci];
        const labelList& compactRegions = recvNonLocal[proci];

        forAll(wantedRegions, i)
        {
            globalToCompact.insert(wantedRegions[i], compactRegions[i]);
        }
    }

    // Finally renumber the regions
    forAll(cellRegion, celli)
    {
        cellRegion[celli] = globalToCompact[cellRegion[celli]];
    }

    return globalCompactPtr;
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::regionSplit::regionSplit(const polyMesh& mesh, const bool doGlobalRegions)
:
    MeshObject<polyMesh, Foam::TopologicalMeshObject, regionSplit>(mesh),
    labelList(mesh.nCells(), -1)
{
    globalNumberingPtr_ = calcRegionSplit
    (
        doGlobalRegions,    //do global regions
        boolList(0, false), //blockedFaces
        List<labelPair>(0), //explicitConnections,
        *this
    );
}


Foam::regionSplit::regionSplit
(
    const polyMesh& mesh,
    const boolList& blockedFace,
    const bool doGlobalRegions
)
:
    MeshObject<polyMesh, Foam::TopologicalMeshObject, regionSplit>(mesh),
    labelList(mesh.nCells(), -1)
{
    globalNumberingPtr_ = calcRegionSplit
    (
        doGlobalRegions,
        blockedFace,        //blockedFaces
        List<labelPair>(0), //explicitConnections,
        *this
    );
}


Foam::regionSplit::regionSplit
(
    const polyMesh& mesh,
    const boolList& blockedFace,
    const List<labelPair>& explicitConnections,
    const bool doGlobalRegions
)
:
    MeshObject<polyMesh, Foam::TopologicalMeshObject, regionSplit>(mesh),
    labelList(mesh.nCells(), -1)
{
    globalNumberingPtr_ = calcRegionSplit
    (
        doGlobalRegions,
        blockedFace,            //blockedFaces
        explicitConnections,    //explicitConnections,
        *this
    );
}


// ************************************************************************* //

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
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/>.

Class
    Foam::regionSplit

Description
    This class separates the mesh into distinct unconnected regions,
    each of which is then given a label according to globalNumbering().


    Say 6 cells, 3 processors, with single baffle on proc1.

              baffle
                |
    +---+---+---+---+---+---+
    |   |   |   |   |   |   |
    +---+---+---+---+---+---+
      proc0 | proc1 | proc2



    1: determine local regions (uncoupled)

    +---+---+---+---+---+---+
    | 0 | 0 | 0 | 1 | 0 | 0 |
    +---+---+---+---+---+---+
      proc0 | proc1 | proc2



    2: make global

    +---+---+---+---+---+---+
    | 0 | 0 | 1 | 2 | 3 | 3 |
    +---+---+---+---+---+---+
      proc0 | proc1 | proc2



    3: merge connected across procs

    +---+---+---+---+---+---+
    | 0 | 0 | 0 | 2 | 2 | 2 |
    +---+---+---+---+---+---+
      proc0 | proc1 | proc2



    4. determine locally owner regions. determine compact numbering for the
    local regions and send these to all processors that need them:

    proc0 uses regions:
        - 0 which is local to it.
    proc1 uses regions
        - 0 which originates from proc0
        - 2 which is local to it
    proc2 uses regions
        - 2 which originates from proc1

    So proc1 needs to get the compact number for region 0 from proc0 and proc2
    needs to get the compact number for region 2 from proc1:

    +---+---+---+---+---+---+
    | 0 | 0 | 0 | 1 | 1 | 1 |
    +---+---+---+---+---+---+
      proc0 | proc1 | proc2


    Can optionally keep all regions local to the processor.

    Note: does not walk across cyclicAMI/cyclicACMI - since not 'coupled()'
          at the patch level.


SourceFiles
    regionSplit.C

\*---------------------------------------------------------------------------*/

#ifndef regionSplit_H
#define regionSplit_H

#include "globalIndex.H"
#include "labelPair.H"
#include "boolList.H"
#include "MeshObject.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

class polyMesh;

/*---------------------------------------------------------------------------*\
                         Class regionSplit Declaration
\*---------------------------------------------------------------------------*/

class regionSplit
:
    public MeshObject<polyMesh, TopologicalMeshObject, regionSplit>,
    public labelList
{
    // Private data

        mutable autoPtr<globalIndex> globalNumberingPtr_;


    // Private Member Functions

        //- Calculate region split in non-compact (global) numbering.
        void calcNonCompactRegionSplit
        (
            const globalIndex& globalFaces,
            const boolList& blockedFace,
            const List<labelPair>& explicitConnections,
            labelList& cellRegion
        ) const;

        //- Calculate global region split. Return globalIndex.
        autoPtr<globalIndex> calcRegionSplit
        (
            const bool doGlobalRegions,
            const boolList& blockedFace,
            const List<labelPair>& explicitConnections,
            labelList& cellRegion
        ) const;


public:

    //- Runtime type information
    ClassName("regionSplit");


    // Constructors

        //- Construct from mesh
        regionSplit
        (
            const polyMesh&,
            const bool doGlobalRegions = Pstream::parRun()
        );

        //- Construct from mesh and whether face is blocked
        //  NOTE: blockedFace has to be consistent across coupled faces!
        regionSplit
        (
            const polyMesh&,
            const boolList& blockedFace,
            const bool doGlobalRegions = Pstream::parRun()
        );

        //- Construct from mesh and whether face is blocked. Additional explicit
        //  connections between normal boundary faces.
        //  NOTE: blockedFace has to be consistent across coupled faces!
        regionSplit
        (
            const polyMesh&,
            const boolList& blockedFace,
            const List<labelPair>&,
            const bool doGlobalRegions = Pstream::parRun()
        );


    // Member Functions

        //- Return global region numbering
        const globalIndex& globalNumbering() const
        {
            return globalNumberingPtr_();
        }

        //- Return local number of regions
        label nLocalRegions() const
        {
            return globalNumbering().localSize(Pstream::myProcNo());
        }

        //- Return total number of regions
        label nRegions() const
        {
            return globalNumbering().size();
        }
};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //